COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/bipartite_matching.h @ 2542:faaa54ec4520

Last change on this file since 2542:faaa54ec4520 was 2488:da94e3b332f3, checked in by Balazs Dezso, 16 years ago

Bug fix in MaxMatching?

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