COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/bipartite_matching.h @ 2550:f26368148b9c

Last change on this file since 2550:f26368148b9c was 2550:f26368148b9c, checked in by Balazs Dezso, 11 years ago

Changing degree of tournament tree
Bug fix in union find
Small efficiency improvment in bipartite matchings

File size: 53.4 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 MaxWeightedBipartiteMatchingDefaultTraits {
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::ANodeMap<int>.
595    typedef typename BpUGraph::template ANodeMap<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 =
648            MaxWeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
649#endif
650  class MaxWeightedBipartiteMatching {
651  public:
652
653    typedef _Traits Traits;
654    typedef typename Traits::BpUGraph BpUGraph;
655    typedef typename Traits::WeightMap WeightMap;
656    typedef typename Traits::Value Value;
657
658  protected:
659
660    typedef typename Traits::HeapCrossRef HeapCrossRef;
661    typedef typename Traits::Heap Heap;
662
663   
664    typedef typename BpUGraph::Node Node;
665    typedef typename BpUGraph::ANodeIt ANodeIt;
666    typedef typename BpUGraph::BNodeIt BNodeIt;
667    typedef typename BpUGraph::UEdge UEdge;
668    typedef typename BpUGraph::UEdgeIt UEdgeIt;
669    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
670
671    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
672    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
673
674    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
675    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
676
677
678  public:
679
680    /// \brief \ref Exception for uninitialized parameters.
681    ///
682    /// This error represents problems in the initialization
683    /// of the parameters of the algorithms.
684    class UninitializedParameter : public lemon::UninitializedParameter {
685    public:
686      virtual const char* what() const throw() {
687        return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
688      }
689    };
690
691    ///\name Named template parameters
692
693    ///@{
694
695    template <class H, class CR>
696    struct DefHeapTraits : public Traits {
697      typedef CR HeapCrossRef;
698      typedef H Heap;
699      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
700        throw UninitializedParameter();
701      }
702      static Heap *createHeap(HeapCrossRef &) {
703        throw UninitializedParameter();
704      }
705    };
706
707    /// \brief \ref named-templ-param "Named parameter" for setting heap
708    /// and cross reference type
709    ///
710    /// \ref named-templ-param "Named parameter" for setting heap and cross
711    /// reference type
712    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
713    struct DefHeap
714      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
715                                            DefHeapTraits<H, CR> > {
716      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
717                                           DefHeapTraits<H, CR> > Create;
718    };
719
720    template <class H, class CR>
721    struct DefStandardHeapTraits : public Traits {
722      typedef CR HeapCrossRef;
723      typedef H Heap;
724      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
725        return new HeapCrossRef(graph);
726      }
727      static Heap *createHeap(HeapCrossRef &crossref) {
728        return new Heap(crossref);
729      }
730    };
731
732    /// \brief \ref named-templ-param "Named parameter" for setting heap and
733    /// cross reference type with automatic allocation
734    ///
735    /// \ref named-templ-param "Named parameter" for setting heap and cross
736    /// reference type. It can allocate the heap and the cross reference
737    /// object if the cross reference's constructor waits for the graph as
738    /// parameter and the heap's constructor waits for the cross reference.
739    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
740    struct DefStandardHeap
741      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
742                                            DefStandardHeapTraits<H, CR> > {
743      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
744                                           DefStandardHeapTraits<H, CR> >
745      Create;
746    };
747
748    ///@}
749
750
751    /// \brief Constructor.
752    ///
753    /// Constructor of the algorithm.
754    MaxWeightedBipartiteMatching(const BpUGraph& _graph,
755                                 const WeightMap& _weight)
756      : graph(&_graph), weight(&_weight),
757        anode_matching(_graph), bnode_matching(_graph),
758        anode_potential(_graph), bnode_potential(_graph),
759        _heap_cross_ref(0), local_heap_cross_ref(false),
760        _heap(0), local_heap(0) {}
761
762    /// \brief Destructor.
763    ///
764    /// Destructor of the algorithm.
765    ~MaxWeightedBipartiteMatching() {
766      destroyStructures();
767    }
768
769    /// \brief Sets the heap and the cross reference used by algorithm.
770    ///
771    /// Sets the heap and the cross reference used by algorithm.
772    /// If you don't use this function before calling \ref run(),
773    /// it will allocate one. The destuctor deallocates this
774    /// automatically allocated map, of course.
775    /// \return \c (*this)
776    MaxWeightedBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
777      if(local_heap_cross_ref) {
778        delete _heap_cross_ref;
779        local_heap_cross_ref = false;
780      }
781      _heap_cross_ref = &cr;
782      if(local_heap) {
783        delete _heap;
784        local_heap = false;
785      }
786      _heap = &hp;
787      return *this;
788    }
789
790    /// \name Execution control
791    /// The simplest way to execute the algorithm is to use
792    /// one of the member functions called \c run().
793    /// \n
794    /// If you need more control on the execution,
795    /// first you must call \ref init() or one alternative for it.
796    /// Finally \ref start() will perform the matching computation or
797    /// with step-by-step execution you can augment the solution.
798
799    /// @{
800
801    /// \brief Initalize the data structures.
802    ///
803    /// It initalizes the data structures and creates an empty matching.
804    void init() {
805      initStructures();
806      for (ANodeIt it(*graph); it != INVALID; ++it) {
807        anode_matching[it] = INVALID;
808        anode_potential[it] = 0;
809      }
810      for (BNodeIt it(*graph); it != INVALID; ++it) {
811        bnode_matching[it] = INVALID;
812        bnode_potential[it] = 0;
813        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
814          if ((*weight)[jt] > bnode_potential[it]) {
815            bnode_potential[it] = (*weight)[jt];
816          }
817        }
818      }
819      matching_value = 0;
820      matching_size = 0;
821    }
822
823
824    /// \brief An augmenting phase of the weighted matching algorithm
825    ///
826    /// It runs an augmenting phase of the weighted matching
827    /// algorithm. This phase finds the best augmenting path and
828    /// augments only on this paths.
829    ///
830    /// The algorithm consists at most
831    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
832    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
833    /// with binary heap.
834    /// \param decrease If the given parameter true the matching value
835    /// can be decreased in the augmenting phase. If we would like
836    /// to calculate the maximum cardinality maximum weighted matching
837    /// then we should let the algorithm to decrease the matching
838    /// value in order to increase the number of the matching edges.
839    bool augment(bool decrease = false) {
840
841      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
842      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
843
844      Node bestNode = INVALID;
845      Value bestValue = 0;
846
847      _heap->clear();
848      for (ANodeIt it(*graph); it != INVALID; ++it) {
849        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
850      }
851
852      for (ANodeIt it(*graph); it != INVALID; ++it) {
853        if (anode_matching[it] == INVALID) {
854          _heap->push(it, 0);
855        }
856      }
857
858      Value bdistMax = 0;
859      while (!_heap->empty()) {
860        Node anode = _heap->top();
861        Value avalue = _heap->prio();
862        _heap->pop();
863        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
864          if (jt == anode_matching[anode]) continue;
865          Node bnode = graph->bNode(jt);
866          Value bvalue = avalue  - (*weight)[jt] +
867            anode_potential[anode] + bnode_potential[bnode];
868          if (bvalue > bdistMax) {
869            bdistMax = bvalue;
870          }
871          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
872            bdist[bnode] = bvalue;
873            bpred[bnode] = jt;
874          } else continue;
875          if (bnode_matching[bnode] != INVALID) {
876            Node newanode = graph->aNode(bnode_matching[bnode]);
877            switch (_heap->state(newanode)) {
878            case Heap::PRE_HEAP:
879              _heap->push(newanode, bvalue);
880              break;
881            case Heap::IN_HEAP:
882              if (bvalue < (*_heap)[newanode]) {
883                _heap->decrease(newanode, bvalue);
884              }
885              break;
886            case Heap::POST_HEAP:
887              break;
888            }
889          } else {
890            if (bestNode == INVALID ||
891                bnode_potential[bnode] - bvalue > bestValue) {
892              bestValue = bnode_potential[bnode] - bvalue;
893              bestNode = bnode;
894            }
895          }
896        }
897      }
898
899      if (bestNode == INVALID || (!decrease && bestValue < 0)) {
900        return false;
901      }
902
903      matching_value += bestValue;
904      ++matching_size;
905
906      for (BNodeIt it(*graph); it != INVALID; ++it) {
907        if (bpred[it] != INVALID) {
908          bnode_potential[it] -= bdist[it];
909        } else {
910          bnode_potential[it] -= bdistMax;
911        }
912      }
913      for (ANodeIt it(*graph); it != INVALID; ++it) {
914        if (anode_matching[it] != INVALID) {
915          Node bnode = graph->bNode(anode_matching[it]);
916          if (bpred[bnode] != INVALID) {
917            anode_potential[it] += bdist[bnode];
918          } else {
919            anode_potential[it] += bdistMax;
920          }
921        }
922      }
923
924      while (bestNode != INVALID) {
925        UEdge uedge = bpred[bestNode];
926        Node anode = graph->aNode(uedge);
927       
928        bnode_matching[bestNode] = uedge;
929        if (anode_matching[anode] != INVALID) {
930          bestNode = graph->bNode(anode_matching[anode]);
931        } else {
932          bestNode = INVALID;
933        }
934        anode_matching[anode] = uedge;
935      }
936
937
938      return true;
939    }
940
941    /// \brief Starts the algorithm.
942    ///
943    /// Starts the algorithm. It runs augmenting phases until the
944    /// optimal solution reached.
945    ///
946    /// \param maxCardinality If the given value is true it will
947    /// calculate the maximum cardinality maximum matching instead of
948    /// the maximum matching.
949    void start(bool maxCardinality = false) {
950      while (augment(maxCardinality)) {}
951    }
952
953    /// \brief Runs the algorithm.
954    ///
955    /// It just initalize the algorithm and then start it.
956    ///
957    /// \param maxCardinality If the given value is true it will
958    /// calculate the maximum cardinality maximum matching instead of
959    /// the maximum matching.
960    void run(bool maxCardinality = false) {
961      init();
962      start(maxCardinality);
963    }
964
965    /// @}
966
967    /// \name Query Functions
968    /// The result of the %Matching algorithm can be obtained using these
969    /// functions.\n
970    /// Before the use of these functions,
971    /// either run() or start() must be called.
972   
973    ///@{
974
975    /// \brief Gives back the potential in the NodeMap
976    ///
977    /// Gives back the potential in the NodeMap. The matching is optimal
978    /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
979    /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
980    /// for each edges.
981    template <typename PotentialMap>
982    void potential(PotentialMap& pt) const {
983      for (ANodeIt it(*graph); it != INVALID; ++it) {
984        pt.set(it, anode_potential[it]);
985      }
986      for (BNodeIt it(*graph); it != INVALID; ++it) {
987        pt.set(it, bnode_potential[it]);
988      }
989    }
990
991    /// \brief Set true all matching uedge in the map.
992    ///
993    /// Set true all matching uedge in the map. It does not change the
994    /// value mapped to the other uedges.
995    /// \return The number of the matching edges.
996    template <typename MatchingMap>
997    int quickMatching(MatchingMap& mm) const {
998      for (ANodeIt it(*graph); it != INVALID; ++it) {
999        if (anode_matching[it] != INVALID) {
1000          mm.set(anode_matching[it], true);
1001        }
1002      }
1003      return matching_size;
1004    }
1005
1006    /// \brief Set true all matching uedge in the map and the others to false.
1007    ///
1008    /// Set true all matching uedge in the map and the others to false.
1009    /// \return The number of the matching edges.
1010    template <typename MatchingMap>
1011    int matching(MatchingMap& mm) const {
1012      for (UEdgeIt it(*graph); it != INVALID; ++it) {
1013        mm.set(it, it == anode_matching[graph->aNode(it)]);
1014      }
1015      return matching_size;
1016    }
1017
1018    ///Gives back the matching in an ANodeMap.
1019
1020    ///Gives back the matching in an ANodeMap. The parameter should
1021    ///be a write ANodeMap of UEdge values.
1022    ///\return The number of the matching edges.
1023    template<class MatchingMap>
1024    int aMatching(MatchingMap& mm) const {
1025      for (ANodeIt it(*graph); it != INVALID; ++it) {
1026        mm.set(it, anode_matching[it]);
1027      }
1028      return matching_size;
1029    }
1030
1031    ///Gives back the matching in a BNodeMap.
1032
1033    ///Gives back the matching in a BNodeMap. The parameter should
1034    ///be a write BNodeMap of UEdge values.
1035    ///\return The number of the matching edges.
1036    template<class MatchingMap>
1037    int bMatching(MatchingMap& mm) const {
1038      for (BNodeIt it(*graph); it != INVALID; ++it) {
1039        mm.set(it, bnode_matching[it]);
1040      }
1041      return matching_size;
1042    }
1043
1044
1045    /// \brief Return true if the given uedge is in the matching.
1046    ///
1047    /// It returns true if the given uedge is in the matching.
1048    bool matchingEdge(const UEdge& edge) const {
1049      return anode_matching[graph->aNode(edge)] == edge;
1050    }
1051
1052    /// \brief Returns the matching edge from the node.
1053    ///
1054    /// Returns the matching edge from the node. If there is not such
1055    /// edge it gives back \c INVALID.
1056    UEdge matchingEdge(const Node& node) const {
1057      if (graph->aNode(node)) {
1058        return anode_matching[node];
1059      } else {
1060        return bnode_matching[node];
1061      }
1062    }
1063
1064    /// \brief Gives back the sum of weights of the matching edges.
1065    ///
1066    /// Gives back the sum of weights of the matching edges.
1067    Value matchingValue() const {
1068      return matching_value;
1069    }
1070
1071    /// \brief Gives back the number of the matching edges.
1072    ///
1073    /// Gives back the number of the matching edges.
1074    int matchingSize() const {
1075      return matching_size;
1076    }
1077
1078    /// @}
1079
1080  private:
1081
1082    void initStructures() {
1083      if (!_heap_cross_ref) {
1084        local_heap_cross_ref = true;
1085        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1086      }
1087      if (!_heap) {
1088        local_heap = true;
1089        _heap = Traits::createHeap(*_heap_cross_ref);
1090      }
1091    }
1092
1093    void destroyStructures() {
1094      if (local_heap_cross_ref) delete _heap_cross_ref;
1095      if (local_heap) delete _heap;
1096    }
1097
1098
1099  private:
1100   
1101    const BpUGraph *graph;
1102    const WeightMap* weight;
1103
1104    ANodeMatchingMap anode_matching;
1105    BNodeMatchingMap bnode_matching;
1106
1107    ANodePotentialMap anode_potential;
1108    BNodePotentialMap bnode_potential;
1109
1110    Value matching_value;
1111    int matching_size;
1112
1113    HeapCrossRef *_heap_cross_ref;
1114    bool local_heap_cross_ref;
1115
1116    Heap *_heap;
1117    bool local_heap;
1118 
1119  };
1120
1121  /// \ingroup matching
1122  ///
1123  /// \brief Maximum weighted bipartite matching
1124  ///
1125  /// This function calculates the maximum weighted matching
1126  /// in a bipartite graph. It gives back the matching in an undirected
1127  /// edge map.
1128  ///
1129  /// \param graph The bipartite graph.
1130  /// \param weight The undirected edge map which contains the weights.
1131  /// \retval matching The undirected edge map which will be set to
1132  /// the matching.
1133  /// \return The value of the matching.
1134  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1135  typename WeightMap::Value
1136  maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
1137                               MatchingMap& matching) {
1138    MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1139      bpmatching(graph, weight);
1140    bpmatching.run();
1141    bpmatching.matching(matching);
1142    return bpmatching.matchingValue();
1143  }
1144
1145  /// \ingroup matching
1146  ///
1147  /// \brief Maximum weighted maximum cardinality bipartite matching
1148  ///
1149  /// This function calculates the maximum weighted of the maximum cardinality
1150  /// matchings of a bipartite graph. It gives back the matching in an
1151  /// undirected edge map.
1152  ///
1153  /// \param graph The bipartite graph.
1154  /// \param weight The undirected edge map which contains the weights.
1155  /// \retval matching The undirected edge map which will be set to
1156  /// the matching.
1157  /// \return The value of the matching.
1158  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1159  typename WeightMap::Value
1160  maxWeightedMaxBipartiteMatching(const BpUGraph& graph,
1161                                  const WeightMap& weight,
1162                                  MatchingMap& matching) {
1163    MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1164      bpmatching(graph, weight);
1165    bpmatching.run(true);
1166    bpmatching.matching(matching);
1167    return bpmatching.matchingValue();
1168  }
1169
1170  /// \brief Default traits class for minimum cost bipartite matching
1171  /// algoritms.
1172  ///
1173  /// Default traits class for minimum cost bipartite matching
1174  /// algoritms. 
1175  ///
1176  /// \param _BpUGraph The bipartite undirected graph
1177  /// type. 
1178  ///
1179  /// \param _CostMap Type of cost map.
1180  template <typename _BpUGraph, typename _CostMap>
1181  struct MinCostMaxBipartiteMatchingDefaultTraits {
1182    /// \brief The type of the cost of the undirected edges.
1183    typedef typename _CostMap::Value Value;
1184
1185    /// The undirected bipartite graph type the algorithm runs on.
1186    typedef _BpUGraph BpUGraph;
1187
1188    /// The map of the edges costs
1189    typedef _CostMap CostMap;
1190
1191    /// \brief The cross reference type used by heap.
1192    ///
1193    /// The cross reference type used by heap.
1194    /// Usually it is \c Graph::NodeMap<int>.
1195    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1196
1197    /// \brief Instantiates a HeapCrossRef.
1198    ///
1199    /// This function instantiates a \ref HeapCrossRef.
1200    /// \param graph is the graph, to which we would like to define the
1201    /// HeapCrossRef.
1202    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1203      return new HeapCrossRef(graph);
1204    }
1205   
1206    /// \brief The heap type used by costed matching algorithms.
1207    ///
1208    /// The heap type used by costed matching algorithms. It should
1209    /// minimize the priorities and the heap's key type is the graph's
1210    /// anode graph's node.
1211    ///
1212    /// \sa BinHeap
1213    typedef BinHeap<Value, HeapCrossRef> Heap;
1214   
1215    /// \brief Instantiates a Heap.
1216    ///
1217    /// This function instantiates a \ref Heap.
1218    /// \param crossref The cross reference of the heap.
1219    static Heap *createHeap(HeapCrossRef& crossref) {
1220      return new Heap(crossref);
1221    }
1222
1223  };
1224
1225
1226  /// \ingroup matching
1227  ///
1228  /// \brief Bipartite Min Cost Matching algorithm
1229  ///
1230  /// This class implements the bipartite Min Cost Matching algorithm.
1231  /// It uses the successive shortest path algorithm to calculate the
1232  /// minimum cost maximum matching in the bipartite graph. The time
1233  /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
1234  /// default binary heap implementation but this can be improved to
1235  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1236  ///
1237  /// The algorithm also provides a potential function on the nodes
1238  /// which a dual solution of the matching algorithm and it can be
1239  /// used to proof the optimality of the given pimal solution.
1240#ifdef DOXYGEN
1241  template <typename _BpUGraph, typename _CostMap, typename _Traits>
1242#else
1243  template <typename _BpUGraph,
1244            typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
1245            typename _Traits =
1246            MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
1247#endif
1248  class MinCostMaxBipartiteMatching {
1249  public:
1250
1251    typedef _Traits Traits;
1252    typedef typename Traits::BpUGraph BpUGraph;
1253    typedef typename Traits::CostMap CostMap;
1254    typedef typename Traits::Value Value;
1255
1256  protected:
1257
1258    typedef typename Traits::HeapCrossRef HeapCrossRef;
1259    typedef typename Traits::Heap Heap;
1260
1261   
1262    typedef typename BpUGraph::Node Node;
1263    typedef typename BpUGraph::ANodeIt ANodeIt;
1264    typedef typename BpUGraph::BNodeIt BNodeIt;
1265    typedef typename BpUGraph::UEdge UEdge;
1266    typedef typename BpUGraph::UEdgeIt UEdgeIt;
1267    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1268
1269    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1270    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1271
1272    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1273    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1274
1275
1276  public:
1277
1278    /// \brief \ref Exception for uninitialized parameters.
1279    ///
1280    /// This error represents problems in the initialization
1281    /// of the parameters of the algorithms.
1282    class UninitializedParameter : public lemon::UninitializedParameter {
1283    public:
1284      virtual const char* what() const throw() {
1285        return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
1286      }
1287    };
1288
1289    ///\name Named template parameters
1290
1291    ///@{
1292
1293    template <class H, class CR>
1294    struct DefHeapTraits : public Traits {
1295      typedef CR HeapCrossRef;
1296      typedef H Heap;
1297      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1298        throw UninitializedParameter();
1299      }
1300      static Heap *createHeap(HeapCrossRef &) {
1301        throw UninitializedParameter();
1302      }
1303    };
1304
1305    /// \brief \ref named-templ-param "Named parameter" for setting heap
1306    /// and cross reference type
1307    ///
1308    /// \ref named-templ-param "Named parameter" for setting heap and cross
1309    /// reference type
1310    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1311    struct DefHeap
1312      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1313                                            DefHeapTraits<H, CR> > {
1314      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1315                                           DefHeapTraits<H, CR> > Create;
1316    };
1317
1318    template <class H, class CR>
1319    struct DefStandardHeapTraits : public Traits {
1320      typedef CR HeapCrossRef;
1321      typedef H Heap;
1322      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1323        return new HeapCrossRef(graph);
1324      }
1325      static Heap *createHeap(HeapCrossRef &crossref) {
1326        return new Heap(crossref);
1327      }
1328    };
1329
1330    /// \brief \ref named-templ-param "Named parameter" for setting heap and
1331    /// cross reference type with automatic allocation
1332    ///
1333    /// \ref named-templ-param "Named parameter" for setting heap and cross
1334    /// reference type. It can allocate the heap and the cross reference
1335    /// object if the cross reference's constructor waits for the graph as
1336    /// parameter and the heap's constructor waits for the cross reference.
1337    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1338    struct DefStandardHeap
1339      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1340                                            DefStandardHeapTraits<H, CR> > {
1341      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1342                                           DefStandardHeapTraits<H, CR> >
1343      Create;
1344    };
1345
1346    ///@}
1347
1348
1349    /// \brief Constructor.
1350    ///
1351    /// Constructor of the algorithm.
1352    MinCostMaxBipartiteMatching(const BpUGraph& _graph,
1353                                 const CostMap& _cost)
1354      : graph(&_graph), cost(&_cost),
1355        anode_matching(_graph), bnode_matching(_graph),
1356        anode_potential(_graph), bnode_potential(_graph),
1357        _heap_cross_ref(0), local_heap_cross_ref(false),
1358        _heap(0), local_heap(0) {}
1359
1360    /// \brief Destructor.
1361    ///
1362    /// Destructor of the algorithm.
1363    ~MinCostMaxBipartiteMatching() {
1364      destroyStructures();
1365    }
1366
1367    /// \brief Sets the heap and the cross reference used by algorithm.
1368    ///
1369    /// Sets the heap and the cross reference used by algorithm.
1370    /// If you don't use this function before calling \ref run(),
1371    /// it will allocate one. The destuctor deallocates this
1372    /// automatically allocated map, of course.
1373    /// \return \c (*this)
1374    MinCostMaxBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
1375      if(local_heap_cross_ref) {
1376        delete _heap_cross_ref;
1377        local_heap_cross_ref = false;
1378      }
1379      _heap_cross_ref = &cr;
1380      if(local_heap) {
1381        delete _heap;
1382        local_heap = false;
1383      }
1384      _heap = &hp;
1385      return *this;
1386    }
1387
1388    /// \name Execution control
1389    /// The simplest way to execute the algorithm is to use
1390    /// one of the member functions called \c run().
1391    /// \n
1392    /// If you need more control on the execution,
1393    /// first you must call \ref init() or one alternative for it.
1394    /// Finally \ref start() will perform the matching computation or
1395    /// with step-by-step execution you can augment the solution.
1396
1397    /// @{
1398
1399    /// \brief Initalize the data structures.
1400    ///
1401    /// It initalizes the data structures and creates an empty matching.
1402    void init() {
1403      initStructures();
1404      for (ANodeIt it(*graph); it != INVALID; ++it) {
1405        anode_matching[it] = INVALID;
1406        anode_potential[it] = 0;
1407      }
1408      for (BNodeIt it(*graph); it != INVALID; ++it) {
1409        bnode_matching[it] = INVALID;
1410        bnode_potential[it] = 0;
1411      }
1412      matching_cost = 0;
1413      matching_size = 0;
1414    }
1415
1416
1417    /// \brief An augmenting phase of the costed matching algorithm
1418    ///
1419    /// It runs an augmenting phase of the matching algorithm. The
1420    /// phase finds the best augmenting path and augments only on this
1421    /// paths.
1422    ///
1423    /// The algorithm consists at most
1424    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1425    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1426    /// with binary heap.
1427    bool augment() {
1428
1429      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1430      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1431
1432      Node bestNode = INVALID;
1433      Value bestValue = 0;
1434
1435      _heap->clear();
1436      for (ANodeIt it(*graph); it != INVALID; ++it) {
1437        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1438      }
1439
1440      for (ANodeIt it(*graph); it != INVALID; ++it) {
1441        if (anode_matching[it] == INVALID) {
1442          _heap->push(it, 0);
1443        }
1444      }
1445      Value bdistMax = 0;
1446
1447      while (!_heap->empty()) {
1448        Node anode = _heap->top();
1449        Value avalue = _heap->prio();
1450        _heap->pop();
1451        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1452          if (jt == anode_matching[anode]) continue;
1453          Node bnode = graph->bNode(jt);
1454          Value bvalue = avalue + (*cost)[jt] +
1455            anode_potential[anode] - bnode_potential[bnode];
1456          if (bvalue > bdistMax) {
1457            bdistMax = bvalue;
1458          }
1459          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1460            bdist[bnode] = bvalue;
1461            bpred[bnode] = jt;
1462          } else continue;
1463          if (bnode_matching[bnode] != INVALID) {
1464            Node newanode = graph->aNode(bnode_matching[bnode]);
1465            switch (_heap->state(newanode)) {
1466            case Heap::PRE_HEAP:
1467              _heap->push(newanode, bvalue);
1468              break;
1469            case Heap::IN_HEAP:
1470              if (bvalue < (*_heap)[newanode]) {
1471                _heap->decrease(newanode, bvalue);
1472              }
1473              break;
1474            case Heap::POST_HEAP:
1475              break;
1476            }
1477          } else {
1478            if (bestNode == INVALID ||
1479                bvalue + bnode_potential[bnode] < bestValue) {
1480              bestValue = bvalue + bnode_potential[bnode];
1481              bestNode = bnode;
1482            }
1483          }
1484        }
1485      }
1486
1487      if (bestNode == INVALID) {
1488        return false;
1489      }
1490
1491      matching_cost += bestValue;
1492      ++matching_size;
1493
1494      for (BNodeIt it(*graph); it != INVALID; ++it) {
1495        if (bpred[it] != INVALID) {
1496          bnode_potential[it] += bdist[it];
1497        } else {
1498          bnode_potential[it] += bdistMax;
1499        }
1500      }
1501      for (ANodeIt it(*graph); it != INVALID; ++it) {
1502        if (anode_matching[it] != INVALID) {
1503          Node bnode = graph->bNode(anode_matching[it]);
1504          if (bpred[bnode] != INVALID) {
1505            anode_potential[it] += bdist[bnode];
1506          } else {
1507            anode_potential[it] += bdistMax;
1508          }
1509        }
1510      }
1511
1512      while (bestNode != INVALID) {
1513        UEdge uedge = bpred[bestNode];
1514        Node anode = graph->aNode(uedge);
1515       
1516        bnode_matching[bestNode] = uedge;
1517        if (anode_matching[anode] != INVALID) {
1518          bestNode = graph->bNode(anode_matching[anode]);
1519        } else {
1520          bestNode = INVALID;
1521        }
1522        anode_matching[anode] = uedge;
1523      }
1524
1525
1526      return true;
1527    }
1528
1529    /// \brief Starts the algorithm.
1530    ///
1531    /// Starts the algorithm. It runs augmenting phases until the
1532    /// optimal solution reached.
1533    void start() {
1534      while (augment()) {}
1535    }
1536
1537    /// \brief Runs the algorithm.
1538    ///
1539    /// It just initalize the algorithm and then start it.
1540    void run() {
1541      init();
1542      start();
1543    }
1544
1545    /// @}
1546
1547    /// \name Query Functions
1548    /// The result of the %Matching algorithm can be obtained using these
1549    /// functions.\n
1550    /// Before the use of these functions,
1551    /// either run() or start() must be called.
1552   
1553    ///@{
1554
1555    /// \brief Gives back the potential in the NodeMap
1556    ///
1557    /// Gives back the potential in the NodeMap. The matching is optimal
1558    /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
1559    /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
1560    /// for each edges.
1561    template <typename PotentialMap>
1562    void potential(PotentialMap& pt) const {
1563      for (ANodeIt it(*graph); it != INVALID; ++it) {
1564        pt.set(it, anode_potential[it]);
1565      }
1566      for (BNodeIt it(*graph); it != INVALID; ++it) {
1567        pt.set(it, bnode_potential[it]);
1568      }
1569    }
1570
1571    /// \brief Set true all matching uedge in the map.
1572    ///
1573    /// Set true all matching uedge in the map. It does not change the
1574    /// value mapped to the other uedges.
1575    /// \return The number of the matching edges.
1576    template <typename MatchingMap>
1577    int quickMatching(MatchingMap& mm) const {
1578      for (ANodeIt it(*graph); it != INVALID; ++it) {
1579        if (anode_matching[it] != INVALID) {
1580          mm.set(anode_matching[it], true);
1581        }
1582      }
1583      return matching_size;
1584    }
1585
1586    /// \brief Set true all matching uedge in the map and the others to false.
1587    ///
1588    /// Set true all matching uedge in the map and the others to false.
1589    /// \return The number of the matching edges.
1590    template <typename MatchingMap>
1591    int matching(MatchingMap& mm) const {
1592      for (UEdgeIt it(*graph); it != INVALID; ++it) {
1593        mm.set(it, it == anode_matching[graph->aNode(it)]);
1594      }
1595      return matching_size;
1596    }
1597
1598    /// \brief Gives back the matching in an ANodeMap.
1599    ///
1600    /// Gives back the matching in an ANodeMap. The parameter should
1601    /// be a write ANodeMap of UEdge values.
1602    /// \return The number of the matching edges.
1603    template<class MatchingMap>
1604    int aMatching(MatchingMap& mm) const {
1605      for (ANodeIt it(*graph); it != INVALID; ++it) {
1606        mm.set(it, anode_matching[it]);
1607      }
1608      return matching_size;
1609    }
1610
1611    /// \brief Gives back the matching in a BNodeMap.
1612    ///
1613    /// Gives back the matching in a BNodeMap. The parameter should
1614    /// be a write BNodeMap of UEdge values.
1615    /// \return The number of the matching edges.
1616    template<class MatchingMap>
1617    int bMatching(MatchingMap& mm) const {
1618      for (BNodeIt it(*graph); it != INVALID; ++it) {
1619        mm.set(it, bnode_matching[it]);
1620      }
1621      return matching_size;
1622    }
1623
1624    /// \brief Return true if the given uedge is in the matching.
1625    ///
1626    /// It returns true if the given uedge is in the matching.
1627    bool matchingEdge(const UEdge& edge) const {
1628      return anode_matching[graph->aNode(edge)] == edge;
1629    }
1630
1631    /// \brief Returns the matching edge from the node.
1632    ///
1633    /// Returns the matching edge from the node. If there is not such
1634    /// edge it gives back \c INVALID.
1635    UEdge matchingEdge(const Node& node) const {
1636      if (graph->aNode(node)) {
1637        return anode_matching[node];
1638      } else {
1639        return bnode_matching[node];
1640      }
1641    }
1642
1643    /// \brief Gives back the sum of costs of the matching edges.
1644    ///
1645    /// Gives back the sum of costs of the matching edges.
1646    Value matchingCost() const {
1647      return matching_cost;
1648    }
1649
1650    /// \brief Gives back the number of the matching edges.
1651    ///
1652    /// Gives back the number of the matching edges.
1653    int matchingSize() const {
1654      return matching_size;
1655    }
1656
1657    /// @}
1658
1659  private:
1660
1661    void initStructures() {
1662      if (!_heap_cross_ref) {
1663        local_heap_cross_ref = true;
1664        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1665      }
1666      if (!_heap) {
1667        local_heap = true;
1668        _heap = Traits::createHeap(*_heap_cross_ref);
1669      }
1670    }
1671
1672    void destroyStructures() {
1673      if (local_heap_cross_ref) delete _heap_cross_ref;
1674      if (local_heap) delete _heap;
1675    }
1676
1677
1678  private:
1679   
1680    const BpUGraph *graph;
1681    const CostMap* cost;
1682
1683    ANodeMatchingMap anode_matching;
1684    BNodeMatchingMap bnode_matching;
1685
1686    ANodePotentialMap anode_potential;
1687    BNodePotentialMap bnode_potential;
1688
1689    Value matching_cost;
1690    int matching_size;
1691
1692    HeapCrossRef *_heap_cross_ref;
1693    bool local_heap_cross_ref;
1694
1695    Heap *_heap;
1696    bool local_heap;
1697 
1698  };
1699
1700  /// \ingroup matching
1701  ///
1702  /// \brief Minimum cost maximum cardinality bipartite matching
1703  ///
1704  /// This function calculates the maximum cardinality matching with
1705  /// minimum cost of a bipartite graph. It gives back the matching in
1706  /// an undirected edge map.
1707  ///
1708  /// \param graph The bipartite graph.
1709  /// \param cost The undirected edge map which contains the costs.
1710  /// \retval matching The undirected edge map which will be set to
1711  /// the matching.
1712  /// \return The cost of the matching.
1713  template <typename BpUGraph, typename CostMap, typename MatchingMap>
1714  typename CostMap::Value
1715  minCostMaxBipartiteMatching(const BpUGraph& graph,
1716                              const CostMap& cost,
1717                              MatchingMap& matching) {
1718    MinCostMaxBipartiteMatching<BpUGraph, CostMap>
1719      bpmatching(graph, cost);
1720    bpmatching.run();
1721    bpmatching.matching(matching);
1722    return bpmatching.matchingCost();
1723  }
1724
1725}
1726
1727#endif
Note: See TracBrowser for help on using the repository browser.