COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/bipartite_matching.h @ 2117:96efb4fa0736

Last change on this file since 2117:96efb4fa0736 was 2058:0b1fc1566fdb, checked in by Balazs Dezso, 14 years ago

Refinements in bipartite matching algorithms

File size: 49.5 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
1329      while (!_heap->empty()) {
1330        Node anode = _heap->top();
1331        Value avalue = _heap->prio();
1332        _heap->pop();
1333        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1334          if (jt == anode_matching[anode]) continue;
1335          Node bnode = graph->bNode(jt);
1336          Value bvalue = avalue + (*cost)[jt] +
1337            anode_potential[anode] - bnode_potential[bnode];
1338          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1339            bdist[bnode] = bvalue;
1340            bpred[bnode] = jt;
1341          }
1342          if (bnode_matching[bnode] != INVALID) {
1343            Node newanode = graph->aNode(bnode_matching[bnode]);
1344            switch (_heap->state(newanode)) {
1345            case Heap::PRE_HEAP:
1346              _heap->push(newanode, bvalue);
1347              break;
1348            case Heap::IN_HEAP:
1349              if (bvalue < (*_heap)[newanode]) {
1350                _heap->decrease(newanode, bvalue);
1351              }
1352              break;
1353            case Heap::POST_HEAP:
1354              break;
1355            }
1356          } else {
1357            if (bestNode == INVALID ||
1358                bvalue + bnode_potential[bnode] < bestValue) {
1359              bestValue = bvalue + bnode_potential[bnode];
1360              bestNode = bnode;
1361            }
1362          }
1363        }
1364      }
1365
1366      if (bestNode == INVALID) {
1367        return false;
1368      }
1369
1370      matching_cost += bestValue;
1371      ++matching_size;
1372
1373      for (BNodeIt it(*graph); it != INVALID; ++it) {
1374        if (bpred[it] != INVALID) {
1375          bnode_potential[it] += bdist[it];
1376        }
1377      }
1378      for (ANodeIt it(*graph); it != INVALID; ++it) {
1379        if (anode_matching[it] != INVALID) {
1380          Node bnode = graph->bNode(anode_matching[it]);
1381          if (bpred[bnode] != INVALID) {
1382            anode_potential[it] += bdist[bnode];
1383          }
1384        }
1385      }
1386
1387      while (bestNode != INVALID) {
1388        UEdge uedge = bpred[bestNode];
1389        Node anode = graph->aNode(uedge);
1390       
1391        bnode_matching[bestNode] = uedge;
1392        if (anode_matching[anode] != INVALID) {
1393          bestNode = graph->bNode(anode_matching[anode]);
1394        } else {
1395          bestNode = INVALID;
1396        }
1397        anode_matching[anode] = uedge;
1398      }
1399
1400
1401      return true;
1402    }
1403
1404    /// \brief Starts the algorithm.
1405    ///
1406    /// Starts the algorithm. It runs augmenting phases until the
1407    /// optimal solution reached.
1408    void start() {
1409      while (augment()) {}
1410    }
1411
1412    /// \brief Runs the algorithm.
1413    ///
1414    /// It just initalize the algorithm and then start it.
1415    void run() {
1416      init();
1417      start();
1418    }
1419
1420    /// @}
1421
1422    /// \name Query Functions
1423    /// The result of the %Matching algorithm can be obtained using these
1424    /// functions.\n
1425    /// Before the use of these functions,
1426    /// either run() or start() must be called.
1427   
1428    ///@{
1429
1430    /// \brief Gives back the potential in the NodeMap
1431    ///
1432    /// Gives back the potential in the NodeMap. The potential is optimal with
1433    /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
1434    /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
1435    /// for each edges.
1436    template <typename PotentialMap>
1437    void potential(PotentialMap& potential) const {
1438      for (ANodeIt it(*graph); it != INVALID; ++it) {
1439        potential[it] = anode_potential[it];
1440      }
1441      for (BNodeIt it(*graph); it != INVALID; ++it) {
1442        potential[it] = bnode_potential[it];
1443      }
1444    }
1445
1446    /// \brief Set true all matching uedge in the map.
1447    ///
1448    /// Set true all matching uedge in the map. It does not change the
1449    /// value mapped to the other uedges.
1450    /// \return The number of the matching edges.
1451    template <typename MatchingMap>
1452    int quickMatching(MatchingMap& matching) const {
1453      for (ANodeIt it(*graph); it != INVALID; ++it) {
1454        if (anode_matching[it] != INVALID) {
1455          matching[anode_matching[it]] = true;
1456        }
1457      }
1458      return matching_size;
1459    }
1460
1461    /// \brief Set true all matching uedge in the map and the others to false.
1462    ///
1463    /// Set true all matching uedge in the map and the others to false.
1464    /// \return The number of the matching edges.
1465    template <typename MatchingMap>
1466    int matching(MatchingMap& matching) const {
1467      for (UEdgeIt it(*graph); it != INVALID; ++it) {
1468        matching[it] = it == anode_matching[graph->aNode(it)];
1469      }
1470      return matching_size;
1471    }
1472
1473
1474    /// \brief Return true if the given uedge is in the matching.
1475    ///
1476    /// It returns true if the given uedge is in the matching.
1477    bool matchingEdge(const UEdge& edge) const {
1478      return anode_matching[graph->aNode(edge)] == edge;
1479    }
1480
1481    /// \brief Returns the matching edge from the node.
1482    ///
1483    /// Returns the matching edge from the node. If there is not such
1484    /// edge it gives back \c INVALID.
1485    UEdge matchingEdge(const Node& node) const {
1486      if (graph->aNode(node)) {
1487        return anode_matching[node];
1488      } else {
1489        return bnode_matching[node];
1490      }
1491    }
1492
1493    /// \brief Gives back the sum of costs of the matching edges.
1494    ///
1495    /// Gives back the sum of costs of the matching edges.
1496    Value matchingCost() const {
1497      return matching_cost;
1498    }
1499
1500    /// \brief Gives back the number of the matching edges.
1501    ///
1502    /// Gives back the number of the matching edges.
1503    int matchingSize() const {
1504      return matching_size;
1505    }
1506
1507    /// @}
1508
1509  private:
1510
1511    void initStructures() {
1512      if (!_heap_cross_ref) {
1513        local_heap_cross_ref = true;
1514        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1515      }
1516      if (!_heap) {
1517        local_heap = true;
1518        _heap = Traits::createHeap(*_heap_cross_ref);
1519      }
1520    }
1521
1522    void destroyStructures() {
1523      if (local_heap_cross_ref) delete _heap_cross_ref;
1524      if (local_heap) delete _heap;
1525    }
1526
1527
1528  private:
1529   
1530    const BpUGraph *graph;
1531    const CostMap* cost;
1532
1533    ANodeMatchingMap anode_matching;
1534    BNodeMatchingMap bnode_matching;
1535
1536    ANodePotentialMap anode_potential;
1537    BNodePotentialMap bnode_potential;
1538
1539    Value matching_cost;
1540    int matching_size;
1541
1542    HeapCrossRef *_heap_cross_ref;
1543    bool local_heap_cross_ref;
1544
1545    Heap *_heap;
1546    bool local_heap;
1547 
1548  };
1549
1550  /// \ingroup matching
1551  ///
1552  /// \brief Minimum cost maximum cardinality bipartite matching
1553  ///
1554  /// This function calculates the minimum cost matching of the maximum
1555  /// cardinality matchings of a bipartite graph. It gives back the matching
1556  /// in an undirected edge map.
1557  ///
1558  /// \param graph The bipartite graph.
1559  /// \param cost The undirected edge map which contains the costs.
1560  /// \retval matching The undirected edge map which will be set to
1561  /// the matching.
1562  /// \return The cost of the matching.
1563  template <typename BpUGraph, typename CostMap, typename MatchingMap>
1564  typename CostMap::Value
1565  minCostMaxBipartiteMatching(const BpUGraph& graph,
1566                              const CostMap& cost,
1567                              MatchingMap& matching) {
1568    MinCostMaxBipartiteMatching<BpUGraph, CostMap>
1569      bpmatching(graph, cost);
1570    bpmatching.run();
1571    bpmatching.matching(matching);
1572    return bpmatching.matchingCost();
1573  }
1574
1575}
1576
1577#endif
Note: See TracBrowser for help on using the repository browser.