COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/bipartite_matching.h @ 2136:4f64d6b3e9ec

Last change on this file since 2136:4f64d6b3e9ec was 2136:4f64d6b3e9ec, checked in by Balazs Dezso, 18 years ago

Bug fix in MinCostMaxBipartiteMatching?
The augmenting phase have not changed the
unreached nodes' potential which caused invalid
dual solution in some cases

File size: 49.7 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_BIPARTITE_MATCHING
20#define LEMON_BIPARTITE_MATCHING
21
22#include <functional>
23
24#include <lemon/bin_heap.h>
25#include <lemon/maps.h>
26
27#include <iostream>
28
29///\ingroup matching
30///\file
31///\brief Maximum matching algorithms in bipartite graphs.
32
33namespace lemon {
34
35  /// \ingroup matching
36  ///
37  /// \brief Bipartite Max Cardinality Matching algorithm
38  ///
39  /// Bipartite Max Cardinality Matching algorithm. This class implements
40  /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
41  /// complexity.
42  template <typename BpUGraph>
43  class MaxBipartiteMatching {
44  protected:
45
46    typedef BpUGraph Graph;
47
48    typedef typename Graph::Node Node;
49    typedef typename Graph::ANodeIt ANodeIt;
50    typedef typename Graph::BNodeIt BNodeIt;
51    typedef typename Graph::UEdge UEdge;
52    typedef typename Graph::UEdgeIt UEdgeIt;
53    typedef typename Graph::IncEdgeIt IncEdgeIt;
54
55    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
56    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
57
58
59  public:
60
61    /// \brief Constructor.
62    ///
63    /// Constructor of the algorithm.
64    MaxBipartiteMatching(const BpUGraph& _graph)
65      : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
66
67    /// \name Execution control
68    /// The simplest way to execute the algorithm is to use
69    /// one of the member functions called \c run().
70    /// \n
71    /// If you need more control on the execution,
72    /// first you must call \ref init() or one alternative for it.
73    /// Finally \ref start() will perform the matching computation or
74    /// with step-by-step execution you can augment the solution.
75
76    /// @{
77
78    /// \brief Initalize the data structures.
79    ///
80    /// It initalizes the data structures and creates an empty matching.
81    void init() {
82      for (ANodeIt it(*graph); it != INVALID; ++it) {
83        anode_matching[it] = INVALID;
84      }
85      for (BNodeIt it(*graph); it != INVALID; ++it) {
86        bnode_matching[it] = INVALID;
87      }
88      matching_size = 0;
89    }
90
91    /// \brief Initalize the data structures.
92    ///
93    /// It initalizes the data structures and creates a greedy
94    /// matching.  From this matching sometimes it is faster to get
95    /// the matching than from the initial empty matching.
96    void greedyInit() {
97      matching_size = 0;
98      for (BNodeIt it(*graph); it != INVALID; ++it) {
99        bnode_matching[it] = INVALID;
100      }
101      for (ANodeIt it(*graph); it != INVALID; ++it) {
102        anode_matching[it] = INVALID;
103        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
104          if (bnode_matching[graph->bNode(jt)] == INVALID) {
105            anode_matching[it] = jt;
106            bnode_matching[graph->bNode(jt)] = jt;
107            ++matching_size;
108            break;
109          }
110        }
111      }
112    }
113
114    /// \brief Initalize the data structures with an initial matching.
115    ///
116    /// It initalizes the data structures with an initial matching.
117    template <typename MatchingMap>
118    void matchingInit(const MatchingMap& matching) {
119      for (ANodeIt it(*graph); it != INVALID; ++it) {
120        anode_matching[it] = INVALID;
121      }
122      for (BNodeIt it(*graph); it != INVALID; ++it) {
123        bnode_matching[it] = INVALID;
124      }
125      matching_size = 0;
126      for (UEdgeIt it(*graph); it != INVALID; ++it) {
127        if (matching[it]) {
128          ++matching_size;
129          anode_matching[graph->aNode(it)] = it;
130          bnode_matching[graph->bNode(it)] = it;
131        }
132      }
133    }
134
135    /// \brief Initalize the data structures with an initial matching.
136    ///
137    /// It initalizes the data structures with an initial matching.
138    /// \return %True when the given map contains really a matching.
139    template <typename MatchingMap>
140    void checkedMatchingInit(const MatchingMap& matching) {
141      for (ANodeIt it(*graph); it != INVALID; ++it) {
142        anode_matching[it] = INVALID;
143      }
144      for (BNodeIt it(*graph); it != INVALID; ++it) {
145        bnode_matching[it] = INVALID;
146      }
147      matching_size = 0;
148      for (UEdgeIt it(*graph); it != INVALID; ++it) {
149        if (matching[it]) {
150          ++matching_size;
151          if (anode_matching[graph->aNode(it)] != INVALID) {
152            return false;
153          }
154          anode_matching[graph->aNode(it)] = it;
155          if (bnode_matching[graph->aNode(it)] != INVALID) {
156            return false;
157          }
158          bnode_matching[graph->bNode(it)] = it;
159        }
160      }
161      return false;
162    }
163
164    /// \brief An augmenting phase of the Hopcroft-Karp algorithm
165    ///
166    /// It runs an augmenting phase of the Hopcroft-Karp
167    /// algorithm. The phase finds maximum count of edge disjoint
168    /// augmenting paths and augments on these paths. The algorithm
169    /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
170    /// \f$ O(e) \f$ long.
171    bool augment() {
172
173      typename Graph::template ANodeMap<bool> areached(*graph, false);
174      typename Graph::template BNodeMap<bool> breached(*graph, false);
175
176      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
177
178      std::vector<Node> queue, bqueue;
179      for (ANodeIt it(*graph); it != INVALID; ++it) {
180        if (anode_matching[it] == INVALID) {
181          queue.push_back(it);
182          areached[it] = true;
183        }
184      }
185
186      bool success = false;
187
188      while (!success && !queue.empty()) {
189        std::vector<Node> newqueue;
190        for (int i = 0; i < (int)queue.size(); ++i) {
191          Node anode = queue[i];
192          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
193            Node bnode = graph->bNode(jt);
194            if (breached[bnode]) continue;
195            breached[bnode] = true;
196            bpred[bnode] = jt;
197            if (bnode_matching[bnode] == INVALID) {
198              bqueue.push_back(bnode);
199              success = true;
200            } else {           
201              Node newanode = graph->aNode(bnode_matching[bnode]);
202              if (!areached[newanode]) {
203                areached[newanode] = true;
204                newqueue.push_back(newanode);
205              }
206            }
207          }
208        }
209        queue.swap(newqueue);
210      }
211
212      if (success) {
213
214        typename Graph::template ANodeMap<bool> aused(*graph, false);
215       
216        for (int i = 0; i < (int)bqueue.size(); ++i) {
217          Node bnode = bqueue[i];
218
219          bool used = false;
220
221          while (bnode != INVALID) {
222            UEdge uedge = bpred[bnode];
223            Node anode = graph->aNode(uedge);
224           
225            if (aused[anode]) {
226              used = true;
227              break;
228            }
229           
230            bnode = anode_matching[anode] != INVALID ?
231              graph->bNode(anode_matching[anode]) : INVALID;
232           
233          }
234         
235          if (used) continue;
236
237          bnode = bqueue[i];
238          while (bnode != INVALID) {
239            UEdge uedge = bpred[bnode];
240            Node anode = graph->aNode(uedge);
241           
242            bnode_matching[bnode] = uedge;
243
244            bnode = anode_matching[anode] != INVALID ?
245              graph->bNode(anode_matching[anode]) : INVALID;
246           
247            anode_matching[anode] = uedge;
248
249            aused[anode] = true;
250          }
251          ++matching_size;
252
253        }
254      }
255      return success;
256    }
257
258    /// \brief An augmenting phase of the Ford-Fulkerson algorithm
259    ///
260    /// It runs an augmenting phase of the Ford-Fulkerson
261    /// algorithm. The phase finds only one augmenting path and
262    /// augments only on this paths. The algorithm consists at most
263    /// of \f$ O(n) \f$ simple phase and one phase is at most
264    /// \f$ O(e) \f$ long.
265    bool simpleAugment() {
266
267      typename Graph::template ANodeMap<bool> areached(*graph, false);
268      typename Graph::template BNodeMap<bool> breached(*graph, false);
269
270      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
271
272      std::vector<Node> queue;
273      for (ANodeIt it(*graph); it != INVALID; ++it) {
274        if (anode_matching[it] == INVALID) {
275          queue.push_back(it);
276          areached[it] = true;
277        }
278      }
279
280      while (!queue.empty()) {
281        std::vector<Node> newqueue;
282        for (int i = 0; i < (int)queue.size(); ++i) {
283          Node anode = queue[i];
284          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
285            Node bnode = graph->bNode(jt);
286            if (breached[bnode]) continue;
287            breached[bnode] = true;
288            bpred[bnode] = jt;
289            if (bnode_matching[bnode] == INVALID) {
290              while (bnode != INVALID) {
291                UEdge uedge = bpred[bnode];
292                anode = graph->aNode(uedge);
293               
294                bnode_matching[bnode] = uedge;
295               
296                bnode = anode_matching[anode] != INVALID ?
297                  graph->bNode(anode_matching[anode]) : INVALID;
298               
299                anode_matching[anode] = uedge;
300               
301              }
302              ++matching_size;
303              return true;
304            } else {           
305              Node newanode = graph->aNode(bnode_matching[bnode]);
306              if (!areached[newanode]) {
307                areached[newanode] = true;
308                newqueue.push_back(newanode);
309              }
310            }
311          }
312        }
313        queue.swap(newqueue);
314      }
315     
316      return false;
317    }
318
319    /// \brief Starts the algorithm.
320    ///
321    /// Starts the algorithm. It runs augmenting phases until the optimal
322    /// solution reached.
323    void start() {
324      while (augment()) {}
325    }
326
327    /// \brief Runs the algorithm.
328    ///
329    /// It just initalize the algorithm and then start it.
330    void run() {
331      greedyInit();
332      start();
333    }
334
335    /// @}
336
337    /// \name Query Functions
338    /// The result of the %Matching algorithm can be obtained using these
339    /// functions.\n
340    /// Before the use of these functions,
341    /// either run() or start() must be called.
342   
343    ///@{
344
345    /// \brief Returns an minimum covering of the nodes.
346    ///
347    /// The minimum covering set problem is the dual solution of the
348    /// maximum bipartite matching. It provides an solution for this
349    /// problem what is proof of the optimality of the matching.
350    /// \return The size of the cover set.
351    template <typename CoverMap>
352    int coverSet(CoverMap& covering) const {
353
354      typename Graph::template ANodeMap<bool> areached(*graph, false);
355      typename Graph::template BNodeMap<bool> breached(*graph, false);
356     
357      std::vector<Node> queue;
358      for (ANodeIt it(*graph); it != INVALID; ++it) {
359        if (anode_matching[it] == INVALID) {
360          queue.push_back(it);
361        }
362      }
363
364      while (!queue.empty()) {
365        std::vector<Node> newqueue;
366        for (int i = 0; i < (int)queue.size(); ++i) {
367          Node anode = queue[i];
368          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
369            Node bnode = graph->bNode(jt);
370            if (breached[bnode]) continue;
371            breached[bnode] = true;
372            if (bnode_matching[bnode] != INVALID) {
373              Node newanode = graph->aNode(bnode_matching[bnode]);
374              if (!areached[newanode]) {
375                areached[newanode] = true;
376                newqueue.push_back(newanode);
377              }
378            }
379          }
380        }
381        queue.swap(newqueue);
382      }
383
384      int size = 0;
385      for (ANodeIt it(*graph); it != INVALID; ++it) {
386        covering[it] = !areached[it] && anode_matching[it] != INVALID;
387        if (!areached[it] && anode_matching[it] != INVALID) {
388          ++size;
389        }
390      }
391      for (BNodeIt it(*graph); it != INVALID; ++it) {
392        covering[it] = breached[it];
393        if (breached[it]) {
394          ++size;
395        }
396      }
397      return size;
398    }
399
400    /// \brief Set true all matching uedge in the map.
401    ///
402    /// Set true all matching uedge in the map. It does not change the
403    /// value mapped to the other uedges.
404    /// \return The number of the matching edges.
405    template <typename MatchingMap>
406    int quickMatching(MatchingMap& matching) const {
407      for (ANodeIt it(*graph); it != INVALID; ++it) {
408        if (anode_matching[it] != INVALID) {
409          matching[anode_matching[it]] = true;
410        }
411      }
412      return matching_size;
413    }
414
415    /// \brief Set true all matching uedge in the map and the others to false.
416    ///
417    /// Set true all matching uedge in the map and the others to false.
418    /// \return The number of the matching edges.
419    template <typename MatchingMap>
420    int matching(MatchingMap& matching) const {
421      for (UEdgeIt it(*graph); it != INVALID; ++it) {
422        matching[it] = it == anode_matching[graph->aNode(it)];
423      }
424      return matching_size;
425    }
426
427
428    /// \brief Return true if the given uedge is in the matching.
429    ///
430    /// It returns true if the given uedge is in the matching.
431    bool matchingEdge(const UEdge& edge) const {
432      return anode_matching[graph->aNode(edge)] == edge;
433    }
434
435    /// \brief Returns the matching edge from the node.
436    ///
437    /// Returns the matching edge from the node. If there is not such
438    /// edge it gives back \c INVALID.
439    UEdge matchingEdge(const Node& node) const {
440      if (graph->aNode(node)) {
441        return anode_matching[node];
442      } else {
443        return bnode_matching[node];
444      }
445    }
446
447    /// \brief Gives back the number of the matching edges.
448    ///
449    /// Gives back the number of the matching edges.
450    int matchingSize() const {
451      return matching_size;
452    }
453
454    /// @}
455
456  private:
457
458    ANodeMatchingMap anode_matching;
459    BNodeMatchingMap bnode_matching;
460    const Graph *graph;
461
462    int matching_size;
463 
464  };
465
466  /// \ingroup matching
467  ///
468  /// \brief Maximum cardinality bipartite matching
469  ///
470  /// This function calculates the maximum cardinality matching
471  /// in a bipartite graph. It gives back the matching in an undirected
472  /// edge map.
473  ///
474  /// \param graph The bipartite graph.
475  /// \retval matching The undirected edge map which will be set to
476  /// the matching.
477  /// \return The size of the matching.
478  template <typename BpUGraph, typename MatchingMap>
479  int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
480    MaxBipartiteMatching<BpUGraph> bpmatching(graph);
481    bpmatching.run();
482    bpmatching.matching(matching);
483    return bpmatching.matchingSize();
484  }
485
486  /// \brief Default traits class for weighted bipartite matching algoritms.
487  ///
488  /// Default traits class for weighted bipartite matching algoritms.
489  /// \param _BpUGraph The bipartite undirected graph type.
490  /// \param _WeightMap Type of weight map.
491  template <typename _BpUGraph, typename _WeightMap>
492  struct WeightedBipartiteMatchingDefaultTraits {
493    /// \brief The type of the weight of the undirected edges.
494    typedef typename _WeightMap::Value Value;
495
496    /// The undirected bipartite graph type the algorithm runs on.
497    typedef _BpUGraph BpUGraph;
498
499    /// The map of the edges weights
500    typedef _WeightMap WeightMap;
501
502    /// \brief The cross reference type used by heap.
503    ///
504    /// The cross reference type used by heap.
505    /// Usually it is \c Graph::NodeMap<int>.
506    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
507
508    /// \brief Instantiates a HeapCrossRef.
509    ///
510    /// This function instantiates a \ref HeapCrossRef.
511    /// \param graph is the graph, to which we would like to define the
512    /// HeapCrossRef.
513    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
514      return new HeapCrossRef(graph);
515    }
516   
517    /// \brief The heap type used by weighted matching algorithms.
518    ///
519    /// The heap type used by weighted matching algorithms. It should
520    /// minimize the priorities and the heap's key type is the graph's
521    /// anode graph's node.
522    ///
523    /// \sa BinHeap
524    typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
525   
526    /// \brief Instantiates a Heap.
527    ///
528    /// This function instantiates a \ref Heap.
529    /// \param crossref The cross reference of the heap.
530    static Heap *createHeap(HeapCrossRef& crossref) {
531      return new Heap(crossref);
532    }
533
534  };
535
536
537  /// \ingroup matching
538  ///
539  /// \brief Bipartite Max Weighted Matching algorithm
540  ///
541  /// This class implements the bipartite Max Weighted Matching
542  /// algorithm.  It uses the successive shortest path algorithm to
543  /// calculate the maximum weighted matching in the bipartite
544  /// graph. The algorithm can be used also to calculate the maximum
545  /// cardinality maximum weighted matching. The time complexity
546  /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
547  /// heap implementation but this can be improved to
548  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
549  ///
550  /// The algorithm also provides a potential function on the nodes
551  /// which a dual solution of the matching algorithm and it can be
552  /// used to proof the optimality of the given pimal solution.
553#ifdef DOXYGEN
554  template <typename _BpUGraph, typename _WeightMap, typename _Traits>
555#else
556  template <typename _BpUGraph,
557            typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
558            typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
559#endif
560  class MaxWeightedBipartiteMatching {
561  public:
562
563    typedef _Traits Traits;
564    typedef typename Traits::BpUGraph BpUGraph;
565    typedef typename Traits::WeightMap WeightMap;
566    typedef typename Traits::Value Value;
567
568  protected:
569
570    typedef typename Traits::HeapCrossRef HeapCrossRef;
571    typedef typename Traits::Heap Heap;
572
573   
574    typedef typename BpUGraph::Node Node;
575    typedef typename BpUGraph::ANodeIt ANodeIt;
576    typedef typename BpUGraph::BNodeIt BNodeIt;
577    typedef typename BpUGraph::UEdge UEdge;
578    typedef typename BpUGraph::UEdgeIt UEdgeIt;
579    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
580
581    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
582    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
583
584    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
585    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
586
587
588  public:
589
590    /// \brief \ref Exception for uninitialized parameters.
591    ///
592    /// This error represents problems in the initialization
593    /// of the parameters of the algorithms.
594    class UninitializedParameter : public lemon::UninitializedParameter {
595    public:
596      virtual const char* exceptionName() const {
597        return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
598      }
599    };
600
601    ///\name Named template parameters
602
603    ///@{
604
605    template <class H, class CR>
606    struct DefHeapTraits : public Traits {
607      typedef CR HeapCrossRef;
608      typedef H Heap;
609      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
610        throw UninitializedParameter();
611      }
612      static Heap *createHeap(HeapCrossRef &) {
613        throw UninitializedParameter();
614      }
615    };
616
617    /// \brief \ref named-templ-param "Named parameter" for setting heap
618    /// and cross reference type
619    ///
620    /// \ref named-templ-param "Named parameter" for setting heap and cross
621    /// reference type
622    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
623    struct DefHeap
624      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
625                                            DefHeapTraits<H, CR> > {
626      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
627                                           DefHeapTraits<H, CR> > Create;
628    };
629
630    template <class H, class CR>
631    struct DefStandardHeapTraits : public Traits {
632      typedef CR HeapCrossRef;
633      typedef H Heap;
634      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
635        return new HeapCrossRef(graph);
636      }
637      static Heap *createHeap(HeapCrossRef &crossref) {
638        return new Heap(crossref);
639      }
640    };
641
642    /// \brief \ref named-templ-param "Named parameter" for setting heap and
643    /// cross reference type with automatic allocation
644    ///
645    /// \ref named-templ-param "Named parameter" for setting heap and cross
646    /// reference type. It can allocate the heap and the cross reference
647    /// object if the cross reference's constructor waits for the graph as
648    /// parameter and the heap's constructor waits for the cross reference.
649    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
650    struct DefStandardHeap
651      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
652                                            DefStandardHeapTraits<H, CR> > {
653      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
654                                           DefStandardHeapTraits<H, CR> >
655      Create;
656    };
657
658    ///@}
659
660
661    /// \brief Constructor.
662    ///
663    /// Constructor of the algorithm.
664    MaxWeightedBipartiteMatching(const BpUGraph& _graph,
665                                 const WeightMap& _weight)
666      : graph(&_graph), weight(&_weight),
667        anode_matching(_graph), bnode_matching(_graph),
668        anode_potential(_graph), bnode_potential(_graph),
669        _heap_cross_ref(0), local_heap_cross_ref(false),
670        _heap(0), local_heap(0) {}
671
672    /// \brief Destructor.
673    ///
674    /// Destructor of the algorithm.
675    ~MaxWeightedBipartiteMatching() {
676      destroyStructures();
677    }
678
679    /// \brief Sets the heap and the cross reference used by algorithm.
680    ///
681    /// Sets the heap and the cross reference used by algorithm.
682    /// If you don't use this function before calling \ref run(),
683    /// it will allocate one. The destuctor deallocates this
684    /// automatically allocated map, of course.
685    /// \return \c (*this)
686    MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
687      if(local_heap_cross_ref) {
688        delete _heap_cross_ref;
689        local_heap_cross_ref = false;
690      }
691      _heap_cross_ref = &crossRef;
692      if(local_heap) {
693        delete _heap;
694        local_heap = false;
695      }
696      _heap = &heap;
697      return *this;
698    }
699
700    /// \name Execution control
701    /// The simplest way to execute the algorithm is to use
702    /// one of the member functions called \c run().
703    /// \n
704    /// If you need more control on the execution,
705    /// first you must call \ref init() or one alternative for it.
706    /// Finally \ref start() will perform the matching computation or
707    /// with step-by-step execution you can augment the solution.
708
709    /// @{
710
711    /// \brief Initalize the data structures.
712    ///
713    /// It initalizes the data structures and creates an empty matching.
714    void init() {
715      initStructures();
716      for (ANodeIt it(*graph); it != INVALID; ++it) {
717        anode_matching[it] = INVALID;
718        anode_potential[it] = 0;
719      }
720      for (BNodeIt it(*graph); it != INVALID; ++it) {
721        bnode_matching[it] = INVALID;
722        bnode_potential[it] = 0;
723        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
724          if ((*weight)[jt] > bnode_potential[it]) {
725            bnode_potential[it] = (*weight)[jt];
726          }
727        }
728      }
729      matching_value = 0;
730      matching_size = 0;
731    }
732
733
734    /// \brief An augmenting phase of the weighted matching algorithm
735    ///
736    /// It runs an augmenting phase of the weighted matching
737    /// algorithm. The phase finds the best augmenting path and
738    /// augments only on this paths.
739    ///
740    /// The algorithm consists at most
741    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
742    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
743    /// with binary heap.
744    /// \param decrease If the given parameter true the matching value
745    /// can be decreased in the augmenting phase. If we would like
746    /// to calculate the maximum cardinality maximum weighted matching
747    /// then we should let the algorithm to decrease the matching
748    /// value in order to increase the number of the matching edges.
749    bool augment(bool decrease = false) {
750
751      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
752      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
753
754      Node bestNode = INVALID;
755      Value bestValue = 0;
756
757      _heap->clear();
758      for (ANodeIt it(*graph); it != INVALID; ++it) {
759        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
760      }
761
762      for (ANodeIt it(*graph); it != INVALID; ++it) {
763        if (anode_matching[it] == INVALID) {
764          _heap->push(it, 0);
765        }
766      }
767
768      Value bdistMax = 0;
769      while (!_heap->empty()) {
770        Node anode = _heap->top();
771        Value avalue = _heap->prio();
772        _heap->pop();
773        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
774          if (jt == anode_matching[anode]) continue;
775          Node bnode = graph->bNode(jt);
776          Value bvalue = avalue  - (*weight)[jt] +
777            anode_potential[anode] + bnode_potential[bnode];
778          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
779            bdist[bnode] = bvalue;
780            bpred[bnode] = jt;
781          }
782          if (bvalue > bdistMax) {
783            bdistMax = bvalue;
784          }
785          if (bnode_matching[bnode] != INVALID) {
786            Node newanode = graph->aNode(bnode_matching[bnode]);
787            switch (_heap->state(newanode)) {
788            case Heap::PRE_HEAP:
789              _heap->push(newanode, bvalue);
790              break;
791            case Heap::IN_HEAP:
792              if (bvalue < (*_heap)[newanode]) {
793                _heap->decrease(newanode, bvalue);
794              }
795              break;
796            case Heap::POST_HEAP:
797              break;
798            }
799          } else {
800            if (bestNode == INVALID ||
801                bnode_potential[bnode] - bvalue > bestValue) {
802              bestValue = bnode_potential[bnode] - bvalue;
803              bestNode = bnode;
804            }
805          }
806        }
807      }
808
809      if (bestNode == INVALID || (!decrease && bestValue < 0)) {
810        return false;
811      }
812
813      matching_value += bestValue;
814      ++matching_size;
815
816      for (BNodeIt it(*graph); it != INVALID; ++it) {
817        if (bpred[it] != INVALID) {
818          bnode_potential[it] -= bdist[it];
819        } else {
820          bnode_potential[it] -= bdistMax;
821        }
822      }
823      for (ANodeIt it(*graph); it != INVALID; ++it) {
824        if (anode_matching[it] != INVALID) {
825          Node bnode = graph->bNode(anode_matching[it]);
826          if (bpred[bnode] != INVALID) {
827            anode_potential[it] += bdist[bnode];
828          } else {
829            anode_potential[it] += bdistMax;
830          }
831        }
832      }
833
834      while (bestNode != INVALID) {
835        UEdge uedge = bpred[bestNode];
836        Node anode = graph->aNode(uedge);
837       
838        bnode_matching[bestNode] = uedge;
839        if (anode_matching[anode] != INVALID) {
840          bestNode = graph->bNode(anode_matching[anode]);
841        } else {
842          bestNode = INVALID;
843        }
844        anode_matching[anode] = uedge;
845      }
846
847
848      return true;
849    }
850
851    /// \brief Starts the algorithm.
852    ///
853    /// Starts the algorithm. It runs augmenting phases until the
854    /// optimal solution reached.
855    ///
856    /// \param maxCardinality If the given value is true it will
857    /// calculate the maximum cardinality maximum matching instead of
858    /// the maximum matching.
859    void start(bool maxCardinality = false) {
860      while (augment(maxCardinality)) {}
861    }
862
863    /// \brief Runs the algorithm.
864    ///
865    /// It just initalize the algorithm and then start it.
866    ///
867    /// \param maxCardinality If the given value is true it will
868    /// calculate the maximum cardinality maximum matching instead of
869    /// the maximum matching.
870    void run(bool maxCardinality = false) {
871      init();
872      start(maxCardinality);
873    }
874
875    /// @}
876
877    /// \name Query Functions
878    /// The result of the %Matching algorithm can be obtained using these
879    /// functions.\n
880    /// Before the use of these functions,
881    /// either run() or start() must be called.
882   
883    ///@{
884
885    /// \brief Gives back the potential in the NodeMap
886    ///
887    /// Gives back the potential in the NodeMap. The matching is optimal
888    /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
889    /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
890    /// for each edges.
891    template <typename PotentialMap>
892    void potential(PotentialMap& potential) const {
893      for (ANodeIt it(*graph); it != INVALID; ++it) {
894        potential[it] = anode_potential[it];
895      }
896      for (BNodeIt it(*graph); it != INVALID; ++it) {
897        potential[it] = bnode_potential[it];
898      }
899    }
900
901    /// \brief Set true all matching uedge in the map.
902    ///
903    /// Set true all matching uedge in the map. It does not change the
904    /// value mapped to the other uedges.
905    /// \return The number of the matching edges.
906    template <typename MatchingMap>
907    int quickMatching(MatchingMap& matching) const {
908      for (ANodeIt it(*graph); it != INVALID; ++it) {
909        if (anode_matching[it] != INVALID) {
910          matching[anode_matching[it]] = true;
911        }
912      }
913      return matching_size;
914    }
915
916    /// \brief Set true all matching uedge in the map and the others to false.
917    ///
918    /// Set true all matching uedge in the map and the others to false.
919    /// \return The number of the matching edges.
920    template <typename MatchingMap>
921    int matching(MatchingMap& matching) const {
922      for (UEdgeIt it(*graph); it != INVALID; ++it) {
923        matching[it] = it == anode_matching[graph->aNode(it)];
924      }
925      return matching_size;
926    }
927
928
929    /// \brief Return true if the given uedge is in the matching.
930    ///
931    /// It returns true if the given uedge is in the matching.
932    bool matchingEdge(const UEdge& edge) const {
933      return anode_matching[graph->aNode(edge)] == edge;
934    }
935
936    /// \brief Returns the matching edge from the node.
937    ///
938    /// Returns the matching edge from the node. If there is not such
939    /// edge it gives back \c INVALID.
940    UEdge matchingEdge(const Node& node) const {
941      if (graph->aNode(node)) {
942        return anode_matching[node];
943      } else {
944        return bnode_matching[node];
945      }
946    }
947
948    /// \brief Gives back the sum of weights of the matching edges.
949    ///
950    /// Gives back the sum of weights of the matching edges.
951    Value matchingValue() const {
952      return matching_value;
953    }
954
955    /// \brief Gives back the number of the matching edges.
956    ///
957    /// Gives back the number of the matching edges.
958    int matchingSize() const {
959      return matching_size;
960    }
961
962    /// @}
963
964  private:
965
966    void initStructures() {
967      if (!_heap_cross_ref) {
968        local_heap_cross_ref = true;
969        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
970      }
971      if (!_heap) {
972        local_heap = true;
973        _heap = Traits::createHeap(*_heap_cross_ref);
974      }
975    }
976
977    void destroyStructures() {
978      if (local_heap_cross_ref) delete _heap_cross_ref;
979      if (local_heap) delete _heap;
980    }
981
982
983  private:
984   
985    const BpUGraph *graph;
986    const WeightMap* weight;
987
988    ANodeMatchingMap anode_matching;
989    BNodeMatchingMap bnode_matching;
990
991    ANodePotentialMap anode_potential;
992    BNodePotentialMap bnode_potential;
993
994    Value matching_value;
995    int matching_size;
996
997    HeapCrossRef *_heap_cross_ref;
998    bool local_heap_cross_ref;
999
1000    Heap *_heap;
1001    bool local_heap;
1002 
1003  };
1004
1005  /// \ingroup matching
1006  ///
1007  /// \brief Maximum weighted bipartite matching
1008  ///
1009  /// This function calculates the maximum weighted matching
1010  /// in a bipartite graph. It gives back the matching in an undirected
1011  /// edge map.
1012  ///
1013  /// \param graph The bipartite graph.
1014  /// \param weight The undirected edge map which contains the weights.
1015  /// \retval matching The undirected edge map which will be set to
1016  /// the matching.
1017  /// \return The value of the matching.
1018  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1019  typename WeightMap::Value
1020  maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
1021                               MatchingMap& matching) {
1022    MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1023      bpmatching(graph, weight);
1024    bpmatching.run();
1025    bpmatching.matching(matching);
1026    return bpmatching.matchingValue();
1027  }
1028
1029  /// \ingroup matching
1030  ///
1031  /// \brief Maximum weighted maximum cardinality bipartite matching
1032  ///
1033  /// This function calculates the maximum weighted of the maximum cardinality
1034  /// matchings of a bipartite graph. It gives back the matching in an
1035  /// undirected edge map.
1036  ///
1037  /// \param graph The bipartite graph.
1038  /// \param weight The undirected edge map which contains the weights.
1039  /// \retval matching The undirected edge map which will be set to
1040  /// the matching.
1041  /// \return The value of the matching.
1042  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1043  typename WeightMap::Value
1044  maxWeightedMaxBipartiteMatching(const BpUGraph& graph,
1045                                  const WeightMap& weight,
1046                                  MatchingMap& matching) {
1047    MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1048      bpmatching(graph, weight);
1049    bpmatching.run(true);
1050    bpmatching.matching(matching);
1051    return bpmatching.matchingValue();
1052  }
1053
1054  /// \brief Default traits class for minimum cost bipartite matching
1055  /// algoritms.
1056  ///
1057  /// Default traits class for minimum cost bipartite matching
1058  /// algoritms. 
1059  ///
1060  /// \param _BpUGraph The bipartite undirected graph
1061  /// type. 
1062  ///
1063  /// \param _CostMap Type of cost map.
1064  template <typename _BpUGraph, typename _CostMap>
1065  struct MinCostMaxBipartiteMatchingDefaultTraits {
1066    /// \brief The type of the cost of the undirected edges.
1067    typedef typename _CostMap::Value Value;
1068
1069    /// The undirected bipartite graph type the algorithm runs on.
1070    typedef _BpUGraph BpUGraph;
1071
1072    /// The map of the edges costs
1073    typedef _CostMap CostMap;
1074
1075    /// \brief The cross reference type used by heap.
1076    ///
1077    /// The cross reference type used by heap.
1078    /// Usually it is \c Graph::NodeMap<int>.
1079    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1080
1081    /// \brief Instantiates a HeapCrossRef.
1082    ///
1083    /// This function instantiates a \ref HeapCrossRef.
1084    /// \param graph is the graph, to which we would like to define the
1085    /// HeapCrossRef.
1086    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1087      return new HeapCrossRef(graph);
1088    }
1089   
1090    /// \brief The heap type used by costed matching algorithms.
1091    ///
1092    /// The heap type used by costed matching algorithms. It should
1093    /// minimize the priorities and the heap's key type is the graph's
1094    /// anode graph's node.
1095    ///
1096    /// \sa BinHeap
1097    typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
1098   
1099    /// \brief Instantiates a Heap.
1100    ///
1101    /// This function instantiates a \ref Heap.
1102    /// \param crossref The cross reference of the heap.
1103    static Heap *createHeap(HeapCrossRef& crossref) {
1104      return new Heap(crossref);
1105    }
1106
1107  };
1108
1109
1110  /// \ingroup matching
1111  ///
1112  /// \brief Bipartite Min Cost Matching algorithm
1113  ///
1114  /// This class implements the bipartite Min Cost Matching algorithm.
1115  /// It uses the successive shortest path algorithm to calculate the
1116  /// minimum cost maximum matching in the bipartite graph. The time
1117  /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
1118  /// default binary heap implementation but this can be improved to
1119  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1120  ///
1121  /// The algorithm also provides a potential function on the nodes
1122  /// which a dual solution of the matching algorithm and it can be
1123  /// used to proof the optimality of the given pimal solution.
1124#ifdef DOXYGEN
1125  template <typename _BpUGraph, typename _CostMap, typename _Traits>
1126#else
1127  template <typename _BpUGraph,
1128            typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
1129            typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
1130#endif
1131  class MinCostMaxBipartiteMatching {
1132  public:
1133
1134    typedef _Traits Traits;
1135    typedef typename Traits::BpUGraph BpUGraph;
1136    typedef typename Traits::CostMap CostMap;
1137    typedef typename Traits::Value Value;
1138
1139  protected:
1140
1141    typedef typename Traits::HeapCrossRef HeapCrossRef;
1142    typedef typename Traits::Heap Heap;
1143
1144   
1145    typedef typename BpUGraph::Node Node;
1146    typedef typename BpUGraph::ANodeIt ANodeIt;
1147    typedef typename BpUGraph::BNodeIt BNodeIt;
1148    typedef typename BpUGraph::UEdge UEdge;
1149    typedef typename BpUGraph::UEdgeIt UEdgeIt;
1150    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1151
1152    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1153    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1154
1155    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1156    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1157
1158
1159  public:
1160
1161    /// \brief \ref Exception for uninitialized parameters.
1162    ///
1163    /// This error represents problems in the initialization
1164    /// of the parameters of the algorithms.
1165    class UninitializedParameter : public lemon::UninitializedParameter {
1166    public:
1167      virtual const char* exceptionName() const {
1168        return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
1169      }
1170    };
1171
1172    ///\name Named template parameters
1173
1174    ///@{
1175
1176    template <class H, class CR>
1177    struct DefHeapTraits : public Traits {
1178      typedef CR HeapCrossRef;
1179      typedef H Heap;
1180      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1181        throw UninitializedParameter();
1182      }
1183      static Heap *createHeap(HeapCrossRef &) {
1184        throw UninitializedParameter();
1185      }
1186    };
1187
1188    /// \brief \ref named-templ-param "Named parameter" for setting heap
1189    /// and cross reference type
1190    ///
1191    /// \ref named-templ-param "Named parameter" for setting heap and cross
1192    /// reference type
1193    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1194    struct DefHeap
1195      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1196                                            DefHeapTraits<H, CR> > {
1197      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1198                                           DefHeapTraits<H, CR> > Create;
1199    };
1200
1201    template <class H, class CR>
1202    struct DefStandardHeapTraits : public Traits {
1203      typedef CR HeapCrossRef;
1204      typedef H Heap;
1205      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1206        return new HeapCrossRef(graph);
1207      }
1208      static Heap *createHeap(HeapCrossRef &crossref) {
1209        return new Heap(crossref);
1210      }
1211    };
1212
1213    /// \brief \ref named-templ-param "Named parameter" for setting heap and
1214    /// cross reference type with automatic allocation
1215    ///
1216    /// \ref named-templ-param "Named parameter" for setting heap and cross
1217    /// reference type. It can allocate the heap and the cross reference
1218    /// object if the cross reference's constructor waits for the graph as
1219    /// parameter and the heap's constructor waits for the cross reference.
1220    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1221    struct DefStandardHeap
1222      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1223                                            DefStandardHeapTraits<H, CR> > {
1224      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1225                                           DefStandardHeapTraits<H, CR> >
1226      Create;
1227    };
1228
1229    ///@}
1230
1231
1232    /// \brief Constructor.
1233    ///
1234    /// Constructor of the algorithm.
1235    MinCostMaxBipartiteMatching(const BpUGraph& _graph,
1236                                 const CostMap& _cost)
1237      : graph(&_graph), cost(&_cost),
1238        anode_matching(_graph), bnode_matching(_graph),
1239        anode_potential(_graph), bnode_potential(_graph),
1240        _heap_cross_ref(0), local_heap_cross_ref(false),
1241        _heap(0), local_heap(0) {}
1242
1243    /// \brief Destructor.
1244    ///
1245    /// Destructor of the algorithm.
1246    ~MinCostMaxBipartiteMatching() {
1247      destroyStructures();
1248    }
1249
1250    /// \brief Sets the heap and the cross reference used by algorithm.
1251    ///
1252    /// Sets the heap and the cross reference used by algorithm.
1253    /// If you don't use this function before calling \ref run(),
1254    /// it will allocate one. The destuctor deallocates this
1255    /// automatically allocated map, of course.
1256    /// \return \c (*this)
1257    MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
1258      if(local_heap_cross_ref) {
1259        delete _heap_cross_ref;
1260        local_heap_cross_ref = false;
1261      }
1262      _heap_cross_ref = &crossRef;
1263      if(local_heap) {
1264        delete _heap;
1265        local_heap = false;
1266      }
1267      _heap = &heap;
1268      return *this;
1269    }
1270
1271    /// \name Execution control
1272    /// The simplest way to execute the algorithm is to use
1273    /// one of the member functions called \c run().
1274    /// \n
1275    /// If you need more control on the execution,
1276    /// first you must call \ref init() or one alternative for it.
1277    /// Finally \ref start() will perform the matching computation or
1278    /// with step-by-step execution you can augment the solution.
1279
1280    /// @{
1281
1282    /// \brief Initalize the data structures.
1283    ///
1284    /// It initalizes the data structures and creates an empty matching.
1285    void init() {
1286      initStructures();
1287      for (ANodeIt it(*graph); it != INVALID; ++it) {
1288        anode_matching[it] = INVALID;
1289        anode_potential[it] = 0;
1290      }
1291      for (BNodeIt it(*graph); it != INVALID; ++it) {
1292        bnode_matching[it] = INVALID;
1293        bnode_potential[it] = 0;
1294      }
1295      matching_cost = 0;
1296      matching_size = 0;
1297    }
1298
1299
1300    /// \brief An augmenting phase of the costed matching algorithm
1301    ///
1302    /// It runs an augmenting phase of the matching algorithm. The
1303    /// phase finds the best augmenting path and augments only on this
1304    /// paths.
1305    ///
1306    /// The algorithm consists at most
1307    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1308    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1309    /// with binary heap.
1310    bool augment() {
1311
1312      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1313      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1314
1315      Node bestNode = INVALID;
1316      Value bestValue = 0;
1317
1318      _heap->clear();
1319      for (ANodeIt it(*graph); it != INVALID; ++it) {
1320        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1321      }
1322
1323      for (ANodeIt it(*graph); it != INVALID; ++it) {
1324        if (anode_matching[it] == INVALID) {
1325          _heap->push(it, 0);
1326        }
1327      }
1328      Value bdistMax = 0;
1329
1330      while (!_heap->empty()) {
1331        Node anode = _heap->top();
1332        Value avalue = _heap->prio();
1333        _heap->pop();
1334        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1335          if (jt == anode_matching[anode]) continue;
1336          Node bnode = graph->bNode(jt);
1337          Value bvalue = avalue + (*cost)[jt] +
1338            anode_potential[anode] - bnode_potential[bnode];
1339          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1340            bdist[bnode] = bvalue;
1341            bpred[bnode] = jt;
1342          }
1343          if (bvalue > bdistMax) {
1344            bdistMax = bvalue;
1345          }
1346          if (bnode_matching[bnode] != INVALID) {
1347            Node newanode = graph->aNode(bnode_matching[bnode]);
1348            switch (_heap->state(newanode)) {
1349            case Heap::PRE_HEAP:
1350              _heap->push(newanode, bvalue);
1351              break;
1352            case Heap::IN_HEAP:
1353              if (bvalue < (*_heap)[newanode]) {
1354                _heap->decrease(newanode, bvalue);
1355              }
1356              break;
1357            case Heap::POST_HEAP:
1358              break;
1359            }
1360          } else {
1361            if (bestNode == INVALID ||
1362                bvalue + bnode_potential[bnode] < bestValue) {
1363              bestValue = bvalue + bnode_potential[bnode];
1364              bestNode = bnode;
1365            }
1366          }
1367        }
1368      }
1369
1370      if (bestNode == INVALID) {
1371        return false;
1372      }
1373
1374      matching_cost += bestValue;
1375      ++matching_size;
1376
1377      for (BNodeIt it(*graph); it != INVALID; ++it) {
1378        if (bpred[it] != INVALID) {
1379          bnode_potential[it] += bdist[it];
1380        } else {
1381          bnode_potential[it] += bdistMax;
1382        }
1383      }
1384      for (ANodeIt it(*graph); it != INVALID; ++it) {
1385        if (anode_matching[it] != INVALID) {
1386          Node bnode = graph->bNode(anode_matching[it]);
1387          if (bpred[bnode] != INVALID) {
1388            anode_potential[it] += bdist[bnode];
1389          } else {
1390            anode_potential[it] += bdistMax;
1391          }
1392        }
1393      }
1394
1395      while (bestNode != INVALID) {
1396        UEdge uedge = bpred[bestNode];
1397        Node anode = graph->aNode(uedge);
1398       
1399        bnode_matching[bestNode] = uedge;
1400        if (anode_matching[anode] != INVALID) {
1401          bestNode = graph->bNode(anode_matching[anode]);
1402        } else {
1403          bestNode = INVALID;
1404        }
1405        anode_matching[anode] = uedge;
1406      }
1407
1408
1409      return true;
1410    }
1411
1412    /// \brief Starts the algorithm.
1413    ///
1414    /// Starts the algorithm. It runs augmenting phases until the
1415    /// optimal solution reached.
1416    void start() {
1417      while (augment()) {}
1418    }
1419
1420    /// \brief Runs the algorithm.
1421    ///
1422    /// It just initalize the algorithm and then start it.
1423    void run() {
1424      init();
1425      start();
1426    }
1427
1428    /// @}
1429
1430    /// \name Query Functions
1431    /// The result of the %Matching algorithm can be obtained using these
1432    /// functions.\n
1433    /// Before the use of these functions,
1434    /// either run() or start() must be called.
1435   
1436    ///@{
1437
1438    /// \brief Gives back the potential in the NodeMap
1439    ///
1440    /// Gives back the potential in the NodeMap. The potential is optimal with
1441    /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
1442    /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
1443    /// for each edges.
1444    template <typename PotentialMap>
1445    void potential(PotentialMap& potential) const {
1446      for (ANodeIt it(*graph); it != INVALID; ++it) {
1447        potential[it] = anode_potential[it];
1448      }
1449      for (BNodeIt it(*graph); it != INVALID; ++it) {
1450        potential[it] = bnode_potential[it];
1451      }
1452    }
1453
1454    /// \brief Set true all matching uedge in the map.
1455    ///
1456    /// Set true all matching uedge in the map. It does not change the
1457    /// value mapped to the other uedges.
1458    /// \return The number of the matching edges.
1459    template <typename MatchingMap>
1460    int quickMatching(MatchingMap& matching) const {
1461      for (ANodeIt it(*graph); it != INVALID; ++it) {
1462        if (anode_matching[it] != INVALID) {
1463          matching[anode_matching[it]] = true;
1464        }
1465      }
1466      return matching_size;
1467    }
1468
1469    /// \brief Set true all matching uedge in the map and the others to false.
1470    ///
1471    /// Set true all matching uedge in the map and the others to false.
1472    /// \return The number of the matching edges.
1473    template <typename MatchingMap>
1474    int matching(MatchingMap& matching) const {
1475      for (UEdgeIt it(*graph); it != INVALID; ++it) {
1476        matching[it] = it == anode_matching[graph->aNode(it)];
1477      }
1478      return matching_size;
1479    }
1480
1481
1482    /// \brief Return true if the given uedge is in the matching.
1483    ///
1484    /// It returns true if the given uedge is in the matching.
1485    bool matchingEdge(const UEdge& edge) const {
1486      return anode_matching[graph->aNode(edge)] == edge;
1487    }
1488
1489    /// \brief Returns the matching edge from the node.
1490    ///
1491    /// Returns the matching edge from the node. If there is not such
1492    /// edge it gives back \c INVALID.
1493    UEdge matchingEdge(const Node& node) const {
1494      if (graph->aNode(node)) {
1495        return anode_matching[node];
1496      } else {
1497        return bnode_matching[node];
1498      }
1499    }
1500
1501    /// \brief Gives back the sum of costs of the matching edges.
1502    ///
1503    /// Gives back the sum of costs of the matching edges.
1504    Value matchingCost() const {
1505      return matching_cost;
1506    }
1507
1508    /// \brief Gives back the number of the matching edges.
1509    ///
1510    /// Gives back the number of the matching edges.
1511    int matchingSize() const {
1512      return matching_size;
1513    }
1514
1515    /// @}
1516
1517  private:
1518
1519    void initStructures() {
1520      if (!_heap_cross_ref) {
1521        local_heap_cross_ref = true;
1522        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1523      }
1524      if (!_heap) {
1525        local_heap = true;
1526        _heap = Traits::createHeap(*_heap_cross_ref);
1527      }
1528    }
1529
1530    void destroyStructures() {
1531      if (local_heap_cross_ref) delete _heap_cross_ref;
1532      if (local_heap) delete _heap;
1533    }
1534
1535
1536  private:
1537   
1538    const BpUGraph *graph;
1539    const CostMap* cost;
1540
1541    ANodeMatchingMap anode_matching;
1542    BNodeMatchingMap bnode_matching;
1543
1544    ANodePotentialMap anode_potential;
1545    BNodePotentialMap bnode_potential;
1546
1547    Value matching_cost;
1548    int matching_size;
1549
1550    HeapCrossRef *_heap_cross_ref;
1551    bool local_heap_cross_ref;
1552
1553    Heap *_heap;
1554    bool local_heap;
1555 
1556  };
1557
1558  /// \ingroup matching
1559  ///
1560  /// \brief Minimum cost maximum cardinality bipartite matching
1561  ///
1562  /// This function calculates the minimum cost matching of the maximum
1563  /// cardinality matchings of a bipartite graph. It gives back the matching
1564  /// in an undirected edge map.
1565  ///
1566  /// \param graph The bipartite graph.
1567  /// \param cost The undirected edge map which contains the costs.
1568  /// \retval matching The undirected edge map which will be set to
1569  /// the matching.
1570  /// \return The cost of the matching.
1571  template <typename BpUGraph, typename CostMap, typename MatchingMap>
1572  typename CostMap::Value
1573  minCostMaxBipartiteMatching(const BpUGraph& graph,
1574                              const CostMap& cost,
1575                              MatchingMap& matching) {
1576    MinCostMaxBipartiteMatching<BpUGraph, CostMap>
1577      bpmatching(graph, cost);
1578    bpmatching.run();
1579    bpmatching.matching(matching);
1580    return bpmatching.matchingCost();
1581  }
1582
1583}
1584
1585#endif
Note: See TracBrowser for help on using the repository browser.