3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
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.
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
19 #ifndef LEMON_BIPARTITE_MATCHING
20 #define LEMON_BIPARTITE_MATCHING
24 #include <lemon/bin_heap.h>
25 #include <lemon/maps.h>
31 ///\brief Maximum matching algorithms in bipartite graphs.
33 ///\note The pr_bipartite_matching.h file also contains algorithms to
34 ///solve maximum cardinality bipartite matching problems.
40 /// \brief Bipartite Max Cardinality Matching algorithm
42 /// Bipartite Max Cardinality Matching algorithm. This class implements
43 /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
46 /// \note In several cases the push-relabel based algorithms have
47 /// better runtime performance than the augmenting path based ones.
49 /// \see PrBipartiteMatching
50 template <typename BpUGraph>
51 class MaxBipartiteMatching {
54 typedef BpUGraph Graph;
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;
63 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
64 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
69 /// \brief Constructor.
71 /// Constructor of the algorithm.
72 MaxBipartiteMatching(const BpUGraph& _graph)
73 : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
75 /// \name Execution control
76 /// The simplest way to execute the algorithm is to use
77 /// one of the member functions called \c run().
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.
86 /// \brief Initalize the data structures.
88 /// It initalizes the data structures and creates an empty matching.
90 for (ANodeIt it(*graph); it != INVALID; ++it) {
91 anode_matching[it] = INVALID;
93 for (BNodeIt it(*graph); it != INVALID; ++it) {
94 bnode_matching[it] = INVALID;
99 /// \brief Initalize the data structures.
101 /// It initalizes the data structures and creates a greedy
102 /// matching. From this matching sometimes it is faster to get
103 /// the matching than from the initial empty matching.
106 for (BNodeIt it(*graph); it != INVALID; ++it) {
107 bnode_matching[it] = INVALID;
109 for (ANodeIt it(*graph); it != INVALID; ++it) {
110 anode_matching[it] = INVALID;
111 for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
112 if (bnode_matching[graph->bNode(jt)] == INVALID) {
113 anode_matching[it] = jt;
114 bnode_matching[graph->bNode(jt)] = jt;
122 /// \brief Initalize the data structures with an initial matching.
124 /// It initalizes the data structures with an initial matching.
125 template <typename MatchingMap>
126 void matchingInit(const MatchingMap& mm) {
127 for (ANodeIt it(*graph); it != INVALID; ++it) {
128 anode_matching[it] = INVALID;
130 for (BNodeIt it(*graph); it != INVALID; ++it) {
131 bnode_matching[it] = INVALID;
134 for (UEdgeIt it(*graph); it != INVALID; ++it) {
137 anode_matching[graph->aNode(it)] = it;
138 bnode_matching[graph->bNode(it)] = it;
143 /// \brief Initalize the data structures with an initial matching.
145 /// It initalizes the data structures with an initial matching.
146 /// \return %True when the given map contains really a matching.
147 template <typename MatchingMap>
148 void checkedMatchingInit(const MatchingMap& mm) {
149 for (ANodeIt it(*graph); it != INVALID; ++it) {
150 anode_matching[it] = INVALID;
152 for (BNodeIt it(*graph); it != INVALID; ++it) {
153 bnode_matching[it] = INVALID;
156 for (UEdgeIt it(*graph); it != INVALID; ++it) {
159 if (anode_matching[graph->aNode(it)] != INVALID) {
162 anode_matching[graph->aNode(it)] = it;
163 if (bnode_matching[graph->aNode(it)] != INVALID) {
166 bnode_matching[graph->bNode(it)] = it;
172 /// \brief An augmenting phase of the Hopcroft-Karp algorithm
174 /// It runs an augmenting phase of the Hopcroft-Karp
175 /// algorithm. This phase finds maximum count of edge disjoint
176 /// augmenting paths and augments on these paths. The algorithm
177 /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
178 /// \f$ O(e) \f$ long.
181 typename Graph::template ANodeMap<bool> areached(*graph, false);
182 typename Graph::template BNodeMap<bool> breached(*graph, false);
184 typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
186 std::vector<Node> queue, bqueue;
187 for (ANodeIt it(*graph); it != INVALID; ++it) {
188 if (anode_matching[it] == INVALID) {
194 bool success = false;
196 while (!success && !queue.empty()) {
197 std::vector<Node> newqueue;
198 for (int i = 0; i < int(queue.size()); ++i) {
199 Node anode = queue[i];
200 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
201 Node bnode = graph->bNode(jt);
202 if (breached[bnode]) continue;
203 breached[bnode] = true;
205 if (bnode_matching[bnode] == INVALID) {
206 bqueue.push_back(bnode);
209 Node newanode = graph->aNode(bnode_matching[bnode]);
210 if (!areached[newanode]) {
211 areached[newanode] = true;
212 newqueue.push_back(newanode);
217 queue.swap(newqueue);
222 typename Graph::template ANodeMap<bool> aused(*graph, false);
224 for (int i = 0; i < int(bqueue.size()); ++i) {
225 Node bnode = bqueue[i];
229 while (bnode != INVALID) {
230 UEdge uedge = bpred[bnode];
231 Node anode = graph->aNode(uedge);
238 bnode = anode_matching[anode] != INVALID ?
239 graph->bNode(anode_matching[anode]) : INVALID;
246 while (bnode != INVALID) {
247 UEdge uedge = bpred[bnode];
248 Node anode = graph->aNode(uedge);
250 bnode_matching[bnode] = uedge;
252 bnode = anode_matching[anode] != INVALID ?
253 graph->bNode(anode_matching[anode]) : INVALID;
255 anode_matching[anode] = uedge;
266 /// \brief An augmenting phase of the Ford-Fulkerson algorithm
268 /// It runs an augmenting phase of the Ford-Fulkerson
269 /// algorithm. This phase finds only one augmenting path and
270 /// augments only on this paths. The algorithm consists at most
271 /// of \f$ O(n) \f$ simple phase and one phase is at most
272 /// \f$ O(e) \f$ long.
273 bool simpleAugment() {
275 typename Graph::template ANodeMap<bool> areached(*graph, false);
276 typename Graph::template BNodeMap<bool> breached(*graph, false);
278 typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
280 std::vector<Node> queue;
281 for (ANodeIt it(*graph); it != INVALID; ++it) {
282 if (anode_matching[it] == INVALID) {
288 while (!queue.empty()) {
289 std::vector<Node> newqueue;
290 for (int i = 0; i < int(queue.size()); ++i) {
291 Node anode = queue[i];
292 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
293 Node bnode = graph->bNode(jt);
294 if (breached[bnode]) continue;
295 breached[bnode] = true;
297 if (bnode_matching[bnode] == INVALID) {
298 while (bnode != INVALID) {
299 UEdge uedge = bpred[bnode];
300 anode = graph->aNode(uedge);
302 bnode_matching[bnode] = uedge;
304 bnode = anode_matching[anode] != INVALID ?
305 graph->bNode(anode_matching[anode]) : INVALID;
307 anode_matching[anode] = uedge;
313 Node newanode = graph->aNode(bnode_matching[bnode]);
314 if (!areached[newanode]) {
315 areached[newanode] = true;
316 newqueue.push_back(newanode);
321 queue.swap(newqueue);
327 /// \brief Starts the algorithm.
329 /// Starts the algorithm. It runs augmenting phases until the optimal
330 /// solution reached.
335 /// \brief Runs the algorithm.
337 /// It just initalize the algorithm and then start it.
345 /// \name Query Functions
346 /// The result of the %Matching algorithm can be obtained using these
348 /// Before the use of these functions,
349 /// either run() or start() must be called.
353 /// \brief Set true all matching uedge in the map.
355 /// Set true all matching uedge in the map. It does not change the
356 /// value mapped to the other uedges.
357 /// \return The number of the matching edges.
358 template <typename MatchingMap>
359 int quickMatching(MatchingMap& mm) const {
360 for (ANodeIt it(*graph); it != INVALID; ++it) {
361 if (anode_matching[it] != INVALID) {
362 mm[anode_matching[it]] = true;
365 return matching_size;
368 /// \brief Set true all matching uedge in the map and the others to false.
370 /// Set true all matching uedge in the map and the others to false.
371 /// \return The number of the matching edges.
372 template <typename MatchingMap>
373 int matching(MatchingMap& mm) const {
374 for (UEdgeIt it(*graph); it != INVALID; ++it) {
375 mm[it] = it == anode_matching[graph->aNode(it)];
377 return matching_size;
381 /// \brief Return true if the given uedge is in the matching.
383 /// It returns true if the given uedge is in the matching.
384 bool matchingEdge(const UEdge& edge) const {
385 return anode_matching[graph->aNode(edge)] == edge;
388 /// \brief Returns the matching edge from the node.
390 /// Returns the matching edge from the node. If there is not such
391 /// edge it gives back \c INVALID.
392 UEdge matchingEdge(const Node& node) const {
393 if (graph->aNode(node)) {
394 return anode_matching[node];
396 return bnode_matching[node];
400 /// \brief Gives back the number of the matching edges.
402 /// Gives back the number of the matching edges.
403 int matchingSize() const {
404 return matching_size;
407 /// \brief Returns a minimum covering of the nodes.
409 /// The minimum covering set problem is the dual solution of the
410 /// maximum bipartite matching. It provides a solution for this
411 /// problem what is proof of the optimality of the matching.
412 /// \return The size of the cover set.
413 template <typename CoverMap>
414 int coverSet(CoverMap& covering) const {
416 typename Graph::template ANodeMap<bool> areached(*graph, false);
417 typename Graph::template BNodeMap<bool> breached(*graph, false);
419 std::vector<Node> queue;
420 for (ANodeIt it(*graph); it != INVALID; ++it) {
421 if (anode_matching[it] == INVALID) {
426 while (!queue.empty()) {
427 std::vector<Node> newqueue;
428 for (int i = 0; i < int(queue.size()); ++i) {
429 Node anode = queue[i];
430 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
431 Node bnode = graph->bNode(jt);
432 if (breached[bnode]) continue;
433 breached[bnode] = true;
434 if (bnode_matching[bnode] != INVALID) {
435 Node newanode = graph->aNode(bnode_matching[bnode]);
436 if (!areached[newanode]) {
437 areached[newanode] = true;
438 newqueue.push_back(newanode);
443 queue.swap(newqueue);
447 for (ANodeIt it(*graph); it != INVALID; ++it) {
448 covering[it] = !areached[it] && anode_matching[it] != INVALID;
449 if (!areached[it] && anode_matching[it] != INVALID) {
453 for (BNodeIt it(*graph); it != INVALID; ++it) {
454 covering[it] = breached[it];
462 /// \brief Gives back a barrier on the A-nodes
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 {
471 typename Graph::template ANodeMap<bool> areached(*graph, false);
472 typename Graph::template BNodeMap<bool> breached(*graph, false);
474 std::vector<Node> queue;
475 for (ANodeIt it(*graph); it != INVALID; ++it) {
476 if (anode_matching[it] == INVALID) {
481 while (!queue.empty()) {
482 std::vector<Node> newqueue;
483 for (int i = 0; i < int(queue.size()); ++i) {
484 Node anode = queue[i];
485 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
486 Node bnode = graph->bNode(jt);
487 if (breached[bnode]) continue;
488 breached[bnode] = true;
489 if (bnode_matching[bnode] != INVALID) {
490 Node newanode = graph->aNode(bnode_matching[bnode]);
491 if (!areached[newanode]) {
492 areached[newanode] = true;
493 newqueue.push_back(newanode);
498 queue.swap(newqueue);
501 for (ANodeIt it(*graph); it != INVALID; ++it) {
502 barrier[it] = areached[it] || anode_matching[it] == INVALID;
506 /// \brief Gives back a barrier on the B-nodes
508 /// The barrier is s subset of the nodes on the same side of the
509 /// graph, which size minus its neighbours is exactly the
510 /// unmatched nodes on the B-side.
511 /// \retval barrier A WriteMap on the BNodes with bool value.
512 template <typename BarrierMap>
513 void bBarrier(BarrierMap& barrier) const {
515 typename Graph::template ANodeMap<bool> areached(*graph, false);
516 typename Graph::template BNodeMap<bool> breached(*graph, false);
518 std::vector<Node> queue;
519 for (ANodeIt it(*graph); it != INVALID; ++it) {
520 if (anode_matching[it] == INVALID) {
525 while (!queue.empty()) {
526 std::vector<Node> newqueue;
527 for (int i = 0; i < int(queue.size()); ++i) {
528 Node anode = queue[i];
529 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
530 Node bnode = graph->bNode(jt);
531 if (breached[bnode]) continue;
532 breached[bnode] = true;
533 if (bnode_matching[bnode] != INVALID) {
534 Node newanode = graph->aNode(bnode_matching[bnode]);
535 if (!areached[newanode]) {
536 areached[newanode] = true;
537 newqueue.push_back(newanode);
542 queue.swap(newqueue);
545 for (BNodeIt it(*graph); it != INVALID; ++it) {
546 barrier[it] = !breached[it];
554 ANodeMatchingMap anode_matching;
555 BNodeMatchingMap bnode_matching;
562 /// \ingroup matching
564 /// \brief Maximum cardinality bipartite matching
566 /// This function calculates the maximum cardinality matching
567 /// in a bipartite graph. It gives back the matching in an undirected
570 /// \param graph The bipartite graph.
571 /// \retval matching The undirected edge map which will be set to
573 /// \return The size of the matching.
574 template <typename BpUGraph, typename MatchingMap>
575 int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
576 MaxBipartiteMatching<BpUGraph> bpmatching(graph);
578 bpmatching.matching(matching);
579 return bpmatching.matchingSize();
582 /// \brief Default traits class for weighted bipartite matching algoritms.
584 /// Default traits class for weighted bipartite matching algoritms.
585 /// \param _BpUGraph The bipartite undirected graph type.
586 /// \param _WeightMap Type of weight map.
587 template <typename _BpUGraph, typename _WeightMap>
588 struct WeightedBipartiteMatchingDefaultTraits {
589 /// \brief The type of the weight of the undirected edges.
590 typedef typename _WeightMap::Value Value;
592 /// The undirected bipartite graph type the algorithm runs on.
593 typedef _BpUGraph BpUGraph;
595 /// The map of the edges weights
596 typedef _WeightMap WeightMap;
598 /// \brief The cross reference type used by heap.
600 /// The cross reference type used by heap.
601 /// Usually it is \c Graph::NodeMap<int>.
602 typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
604 /// \brief Instantiates a HeapCrossRef.
606 /// This function instantiates a \ref HeapCrossRef.
607 /// \param graph is the graph, to which we would like to define the
609 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
610 return new HeapCrossRef(graph);
613 /// \brief The heap type used by weighted matching algorithms.
615 /// The heap type used by weighted matching algorithms. It should
616 /// minimize the priorities and the heap's key type is the graph's
617 /// anode graph's node.
620 typedef BinHeap<Value, HeapCrossRef> Heap;
622 /// \brief Instantiates a Heap.
624 /// This function instantiates a \ref Heap.
625 /// \param crossref The cross reference of the heap.
626 static Heap *createHeap(HeapCrossRef& crossref) {
627 return new Heap(crossref);
633 /// \ingroup matching
635 /// \brief Bipartite Max Weighted Matching algorithm
637 /// This class implements the bipartite Max Weighted Matching
638 /// algorithm. It uses the successive shortest path algorithm to
639 /// calculate the maximum weighted matching in the bipartite
640 /// graph. The algorithm can be used also to calculate the maximum
641 /// cardinality maximum weighted matching. The time complexity
642 /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
643 /// heap implementation but this can be improved to
644 /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
646 /// The algorithm also provides a potential function on the nodes
647 /// which a dual solution of the matching algorithm and it can be
648 /// used to proof the optimality of the given pimal solution.
650 template <typename _BpUGraph, typename _WeightMap, typename _Traits>
652 template <typename _BpUGraph,
653 typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
654 typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
656 class MaxWeightedBipartiteMatching {
659 typedef _Traits Traits;
660 typedef typename Traits::BpUGraph BpUGraph;
661 typedef typename Traits::WeightMap WeightMap;
662 typedef typename Traits::Value Value;
666 typedef typename Traits::HeapCrossRef HeapCrossRef;
667 typedef typename Traits::Heap Heap;
670 typedef typename BpUGraph::Node Node;
671 typedef typename BpUGraph::ANodeIt ANodeIt;
672 typedef typename BpUGraph::BNodeIt BNodeIt;
673 typedef typename BpUGraph::UEdge UEdge;
674 typedef typename BpUGraph::UEdgeIt UEdgeIt;
675 typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
677 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
678 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
680 typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
681 typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
686 /// \brief \ref Exception for uninitialized parameters.
688 /// This error represents problems in the initialization
689 /// of the parameters of the algorithms.
690 class UninitializedParameter : public lemon::UninitializedParameter {
692 virtual const char* what() const throw() {
693 return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
697 ///\name Named template parameters
701 template <class H, class CR>
702 struct DefHeapTraits : public Traits {
703 typedef CR HeapCrossRef;
705 static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
706 throw UninitializedParameter();
708 static Heap *createHeap(HeapCrossRef &) {
709 throw UninitializedParameter();
713 /// \brief \ref named-templ-param "Named parameter" for setting heap
714 /// and cross reference type
716 /// \ref named-templ-param "Named parameter" for setting heap and cross
718 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
720 : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
721 DefHeapTraits<H, CR> > {
722 typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
723 DefHeapTraits<H, CR> > Create;
726 template <class H, class CR>
727 struct DefStandardHeapTraits : public Traits {
728 typedef CR HeapCrossRef;
730 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
731 return new HeapCrossRef(graph);
733 static Heap *createHeap(HeapCrossRef &crossref) {
734 return new Heap(crossref);
738 /// \brief \ref named-templ-param "Named parameter" for setting heap and
739 /// cross reference type with automatic allocation
741 /// \ref named-templ-param "Named parameter" for setting heap and cross
742 /// reference type. It can allocate the heap and the cross reference
743 /// object if the cross reference's constructor waits for the graph as
744 /// parameter and the heap's constructor waits for the cross reference.
745 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
746 struct DefStandardHeap
747 : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
748 DefStandardHeapTraits<H, CR> > {
749 typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
750 DefStandardHeapTraits<H, CR> >
757 /// \brief Constructor.
759 /// Constructor of the algorithm.
760 MaxWeightedBipartiteMatching(const BpUGraph& _graph,
761 const WeightMap& _weight)
762 : graph(&_graph), weight(&_weight),
763 anode_matching(_graph), bnode_matching(_graph),
764 anode_potential(_graph), bnode_potential(_graph),
765 _heap_cross_ref(0), local_heap_cross_ref(false),
766 _heap(0), local_heap(0) {}
768 /// \brief Destructor.
770 /// Destructor of the algorithm.
771 ~MaxWeightedBipartiteMatching() {
775 /// \brief Sets the heap and the cross reference used by algorithm.
777 /// Sets the heap and the cross reference used by algorithm.
778 /// If you don't use this function before calling \ref run(),
779 /// it will allocate one. The destuctor deallocates this
780 /// automatically allocated map, of course.
781 /// \return \c (*this)
782 MaxWeightedBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
783 if(local_heap_cross_ref) {
784 delete _heap_cross_ref;
785 local_heap_cross_ref = false;
787 _heap_cross_ref = &cr;
796 /// \name Execution control
797 /// The simplest way to execute the algorithm is to use
798 /// one of the member functions called \c run().
800 /// If you need more control on the execution,
801 /// first you must call \ref init() or one alternative for it.
802 /// Finally \ref start() will perform the matching computation or
803 /// with step-by-step execution you can augment the solution.
807 /// \brief Initalize the data structures.
809 /// It initalizes the data structures and creates an empty matching.
812 for (ANodeIt it(*graph); it != INVALID; ++it) {
813 anode_matching[it] = INVALID;
814 anode_potential[it] = 0;
816 for (BNodeIt it(*graph); it != INVALID; ++it) {
817 bnode_matching[it] = INVALID;
818 bnode_potential[it] = 0;
819 for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
820 if ((*weight)[jt] > bnode_potential[it]) {
821 bnode_potential[it] = (*weight)[jt];
830 /// \brief An augmenting phase of the weighted matching algorithm
832 /// It runs an augmenting phase of the weighted matching
833 /// algorithm. This phase finds the best augmenting path and
834 /// augments only on this paths.
836 /// The algorithm consists at most
837 /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
838 /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
839 /// with binary heap.
840 /// \param decrease If the given parameter true the matching value
841 /// can be decreased in the augmenting phase. If we would like
842 /// to calculate the maximum cardinality maximum weighted matching
843 /// then we should let the algorithm to decrease the matching
844 /// value in order to increase the number of the matching edges.
845 bool augment(bool decrease = false) {
847 typename BpUGraph::template BNodeMap<Value> bdist(*graph);
848 typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
850 Node bestNode = INVALID;
854 for (ANodeIt it(*graph); it != INVALID; ++it) {
855 (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
858 for (ANodeIt it(*graph); it != INVALID; ++it) {
859 if (anode_matching[it] == INVALID) {
865 while (!_heap->empty()) {
866 Node anode = _heap->top();
867 Value avalue = _heap->prio();
869 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
870 if (jt == anode_matching[anode]) continue;
871 Node bnode = graph->bNode(jt);
872 Value bvalue = avalue - (*weight)[jt] +
873 anode_potential[anode] + bnode_potential[bnode];
874 if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
875 bdist[bnode] = bvalue;
878 if (bvalue > bdistMax) {
881 if (bnode_matching[bnode] != INVALID) {
882 Node newanode = graph->aNode(bnode_matching[bnode]);
883 switch (_heap->state(newanode)) {
885 _heap->push(newanode, bvalue);
888 if (bvalue < (*_heap)[newanode]) {
889 _heap->decrease(newanode, bvalue);
892 case Heap::POST_HEAP:
896 if (bestNode == INVALID ||
897 bnode_potential[bnode] - bvalue > bestValue) {
898 bestValue = bnode_potential[bnode] - bvalue;
905 if (bestNode == INVALID || (!decrease && bestValue < 0)) {
909 matching_value += bestValue;
912 for (BNodeIt it(*graph); it != INVALID; ++it) {
913 if (bpred[it] != INVALID) {
914 bnode_potential[it] -= bdist[it];
916 bnode_potential[it] -= bdistMax;
919 for (ANodeIt it(*graph); it != INVALID; ++it) {
920 if (anode_matching[it] != INVALID) {
921 Node bnode = graph->bNode(anode_matching[it]);
922 if (bpred[bnode] != INVALID) {
923 anode_potential[it] += bdist[bnode];
925 anode_potential[it] += bdistMax;
930 while (bestNode != INVALID) {
931 UEdge uedge = bpred[bestNode];
932 Node anode = graph->aNode(uedge);
934 bnode_matching[bestNode] = uedge;
935 if (anode_matching[anode] != INVALID) {
936 bestNode = graph->bNode(anode_matching[anode]);
940 anode_matching[anode] = uedge;
947 /// \brief Starts the algorithm.
949 /// Starts the algorithm. It runs augmenting phases until the
950 /// optimal solution reached.
952 /// \param maxCardinality If the given value is true it will
953 /// calculate the maximum cardinality maximum matching instead of
954 /// the maximum matching.
955 void start(bool maxCardinality = false) {
956 while (augment(maxCardinality)) {}
959 /// \brief Runs the algorithm.
961 /// It just initalize the algorithm and then start it.
963 /// \param maxCardinality If the given value is true it will
964 /// calculate the maximum cardinality maximum matching instead of
965 /// the maximum matching.
966 void run(bool maxCardinality = false) {
968 start(maxCardinality);
973 /// \name Query Functions
974 /// The result of the %Matching algorithm can be obtained using these
976 /// Before the use of these functions,
977 /// either run() or start() must be called.
981 /// \brief Gives back the potential in the NodeMap
983 /// Gives back the potential in the NodeMap. The matching is optimal
984 /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
985 /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
987 template <typename PotentialMap>
988 void potential(PotentialMap& pt) const {
989 for (ANodeIt it(*graph); it != INVALID; ++it) {
990 pt[it] = anode_potential[it];
992 for (BNodeIt it(*graph); it != INVALID; ++it) {
993 pt[it] = bnode_potential[it];
997 /// \brief Set true all matching uedge in the map.
999 /// Set true all matching uedge in the map. It does not change the
1000 /// value mapped to the other uedges.
1001 /// \return The number of the matching edges.
1002 template <typename MatchingMap>
1003 int quickMatching(MatchingMap& mm) const {
1004 for (ANodeIt it(*graph); it != INVALID; ++it) {
1005 if (anode_matching[it] != INVALID) {
1006 mm[anode_matching[it]] = true;
1009 return matching_size;
1012 /// \brief Set true all matching uedge in the map and the others to false.
1014 /// Set true all matching uedge in the map and the others to false.
1015 /// \return The number of the matching edges.
1016 template <typename MatchingMap>
1017 int matching(MatchingMap& mm) const {
1018 for (UEdgeIt it(*graph); it != INVALID; ++it) {
1019 mm[it] = it == anode_matching[graph->aNode(it)];
1021 return matching_size;
1025 /// \brief Return true if the given uedge is in the matching.
1027 /// It returns true if the given uedge is in the matching.
1028 bool matchingEdge(const UEdge& edge) const {
1029 return anode_matching[graph->aNode(edge)] == edge;
1032 /// \brief Returns the matching edge from the node.
1034 /// Returns the matching edge from the node. If there is not such
1035 /// edge it gives back \c INVALID.
1036 UEdge matchingEdge(const Node& node) const {
1037 if (graph->aNode(node)) {
1038 return anode_matching[node];
1040 return bnode_matching[node];
1044 /// \brief Gives back the sum of weights of the matching edges.
1046 /// Gives back the sum of weights of the matching edges.
1047 Value matchingValue() const {
1048 return matching_value;
1051 /// \brief Gives back the number of the matching edges.
1053 /// Gives back the number of the matching edges.
1054 int matchingSize() const {
1055 return matching_size;
1062 void initStructures() {
1063 if (!_heap_cross_ref) {
1064 local_heap_cross_ref = true;
1065 _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1069 _heap = Traits::createHeap(*_heap_cross_ref);
1073 void destroyStructures() {
1074 if (local_heap_cross_ref) delete _heap_cross_ref;
1075 if (local_heap) delete _heap;
1081 const BpUGraph *graph;
1082 const WeightMap* weight;
1084 ANodeMatchingMap anode_matching;
1085 BNodeMatchingMap bnode_matching;
1087 ANodePotentialMap anode_potential;
1088 BNodePotentialMap bnode_potential;
1090 Value matching_value;
1093 HeapCrossRef *_heap_cross_ref;
1094 bool local_heap_cross_ref;
1101 /// \ingroup matching
1103 /// \brief Maximum weighted bipartite matching
1105 /// This function calculates the maximum weighted matching
1106 /// in a bipartite graph. It gives back the matching in an undirected
1109 /// \param graph The bipartite graph.
1110 /// \param weight The undirected edge map which contains the weights.
1111 /// \retval matching The undirected edge map which will be set to
1113 /// \return The value of the matching.
1114 template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1115 typename WeightMap::Value
1116 maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
1117 MatchingMap& matching) {
1118 MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1119 bpmatching(graph, weight);
1121 bpmatching.matching(matching);
1122 return bpmatching.matchingValue();
1125 /// \ingroup matching
1127 /// \brief Maximum weighted maximum cardinality bipartite matching
1129 /// This function calculates the maximum weighted of the maximum cardinality
1130 /// matchings of a bipartite graph. It gives back the matching in an
1131 /// undirected edge map.
1133 /// \param graph The bipartite graph.
1134 /// \param weight The undirected edge map which contains the weights.
1135 /// \retval matching The undirected edge map which will be set to
1137 /// \return The value of the matching.
1138 template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1139 typename WeightMap::Value
1140 maxWeightedMaxBipartiteMatching(const BpUGraph& graph,
1141 const WeightMap& weight,
1142 MatchingMap& matching) {
1143 MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1144 bpmatching(graph, weight);
1145 bpmatching.run(true);
1146 bpmatching.matching(matching);
1147 return bpmatching.matchingValue();
1150 /// \brief Default traits class for minimum cost bipartite matching
1153 /// Default traits class for minimum cost bipartite matching
1156 /// \param _BpUGraph The bipartite undirected graph
1159 /// \param _CostMap Type of cost map.
1160 template <typename _BpUGraph, typename _CostMap>
1161 struct MinCostMaxBipartiteMatchingDefaultTraits {
1162 /// \brief The type of the cost of the undirected edges.
1163 typedef typename _CostMap::Value Value;
1165 /// The undirected bipartite graph type the algorithm runs on.
1166 typedef _BpUGraph BpUGraph;
1168 /// The map of the edges costs
1169 typedef _CostMap CostMap;
1171 /// \brief The cross reference type used by heap.
1173 /// The cross reference type used by heap.
1174 /// Usually it is \c Graph::NodeMap<int>.
1175 typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1177 /// \brief Instantiates a HeapCrossRef.
1179 /// This function instantiates a \ref HeapCrossRef.
1180 /// \param graph is the graph, to which we would like to define the
1182 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1183 return new HeapCrossRef(graph);
1186 /// \brief The heap type used by costed matching algorithms.
1188 /// The heap type used by costed matching algorithms. It should
1189 /// minimize the priorities and the heap's key type is the graph's
1190 /// anode graph's node.
1193 typedef BinHeap<Value, HeapCrossRef> Heap;
1195 /// \brief Instantiates a Heap.
1197 /// This function instantiates a \ref Heap.
1198 /// \param crossref The cross reference of the heap.
1199 static Heap *createHeap(HeapCrossRef& crossref) {
1200 return new Heap(crossref);
1206 /// \ingroup matching
1208 /// \brief Bipartite Min Cost Matching algorithm
1210 /// This class implements the bipartite Min Cost Matching algorithm.
1211 /// It uses the successive shortest path algorithm to calculate the
1212 /// minimum cost maximum matching in the bipartite graph. The time
1213 /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
1214 /// default binary heap implementation but this can be improved to
1215 /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1217 /// The algorithm also provides a potential function on the nodes
1218 /// which a dual solution of the matching algorithm and it can be
1219 /// used to proof the optimality of the given pimal solution.
1221 template <typename _BpUGraph, typename _CostMap, typename _Traits>
1223 template <typename _BpUGraph,
1224 typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
1225 typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
1227 class MinCostMaxBipartiteMatching {
1230 typedef _Traits Traits;
1231 typedef typename Traits::BpUGraph BpUGraph;
1232 typedef typename Traits::CostMap CostMap;
1233 typedef typename Traits::Value Value;
1237 typedef typename Traits::HeapCrossRef HeapCrossRef;
1238 typedef typename Traits::Heap Heap;
1241 typedef typename BpUGraph::Node Node;
1242 typedef typename BpUGraph::ANodeIt ANodeIt;
1243 typedef typename BpUGraph::BNodeIt BNodeIt;
1244 typedef typename BpUGraph::UEdge UEdge;
1245 typedef typename BpUGraph::UEdgeIt UEdgeIt;
1246 typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1248 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1249 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1251 typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1252 typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1257 /// \brief \ref Exception for uninitialized parameters.
1259 /// This error represents problems in the initialization
1260 /// of the parameters of the algorithms.
1261 class UninitializedParameter : public lemon::UninitializedParameter {
1263 virtual const char* what() const throw() {
1264 return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
1268 ///\name Named template parameters
1272 template <class H, class CR>
1273 struct DefHeapTraits : public Traits {
1274 typedef CR HeapCrossRef;
1276 static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1277 throw UninitializedParameter();
1279 static Heap *createHeap(HeapCrossRef &) {
1280 throw UninitializedParameter();
1284 /// \brief \ref named-templ-param "Named parameter" for setting heap
1285 /// and cross reference type
1287 /// \ref named-templ-param "Named parameter" for setting heap and cross
1289 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1291 : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1292 DefHeapTraits<H, CR> > {
1293 typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1294 DefHeapTraits<H, CR> > Create;
1297 template <class H, class CR>
1298 struct DefStandardHeapTraits : public Traits {
1299 typedef CR HeapCrossRef;
1301 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1302 return new HeapCrossRef(graph);
1304 static Heap *createHeap(HeapCrossRef &crossref) {
1305 return new Heap(crossref);
1309 /// \brief \ref named-templ-param "Named parameter" for setting heap and
1310 /// cross reference type with automatic allocation
1312 /// \ref named-templ-param "Named parameter" for setting heap and cross
1313 /// reference type. It can allocate the heap and the cross reference
1314 /// object if the cross reference's constructor waits for the graph as
1315 /// parameter and the heap's constructor waits for the cross reference.
1316 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1317 struct DefStandardHeap
1318 : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1319 DefStandardHeapTraits<H, CR> > {
1320 typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1321 DefStandardHeapTraits<H, CR> >
1328 /// \brief Constructor.
1330 /// Constructor of the algorithm.
1331 MinCostMaxBipartiteMatching(const BpUGraph& _graph,
1332 const CostMap& _cost)
1333 : graph(&_graph), cost(&_cost),
1334 anode_matching(_graph), bnode_matching(_graph),
1335 anode_potential(_graph), bnode_potential(_graph),
1336 _heap_cross_ref(0), local_heap_cross_ref(false),
1337 _heap(0), local_heap(0) {}
1339 /// \brief Destructor.
1341 /// Destructor of the algorithm.
1342 ~MinCostMaxBipartiteMatching() {
1343 destroyStructures();
1346 /// \brief Sets the heap and the cross reference used by algorithm.
1348 /// Sets the heap and the cross reference used by algorithm.
1349 /// If you don't use this function before calling \ref run(),
1350 /// it will allocate one. The destuctor deallocates this
1351 /// automatically allocated map, of course.
1352 /// \return \c (*this)
1353 MinCostMaxBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
1354 if(local_heap_cross_ref) {
1355 delete _heap_cross_ref;
1356 local_heap_cross_ref = false;
1358 _heap_cross_ref = &cr;
1367 /// \name Execution control
1368 /// The simplest way to execute the algorithm is to use
1369 /// one of the member functions called \c run().
1371 /// If you need more control on the execution,
1372 /// first you must call \ref init() or one alternative for it.
1373 /// Finally \ref start() will perform the matching computation or
1374 /// with step-by-step execution you can augment the solution.
1378 /// \brief Initalize the data structures.
1380 /// It initalizes the data structures and creates an empty matching.
1383 for (ANodeIt it(*graph); it != INVALID; ++it) {
1384 anode_matching[it] = INVALID;
1385 anode_potential[it] = 0;
1387 for (BNodeIt it(*graph); it != INVALID; ++it) {
1388 bnode_matching[it] = INVALID;
1389 bnode_potential[it] = 0;
1396 /// \brief An augmenting phase of the costed matching algorithm
1398 /// It runs an augmenting phase of the matching algorithm. The
1399 /// phase finds the best augmenting path and augments only on this
1402 /// The algorithm consists at most
1403 /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1404 /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1405 /// with binary heap.
1408 typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1409 typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1411 Node bestNode = INVALID;
1412 Value bestValue = 0;
1415 for (ANodeIt it(*graph); it != INVALID; ++it) {
1416 (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1419 for (ANodeIt it(*graph); it != INVALID; ++it) {
1420 if (anode_matching[it] == INVALID) {
1426 while (!_heap->empty()) {
1427 Node anode = _heap->top();
1428 Value avalue = _heap->prio();
1430 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1431 if (jt == anode_matching[anode]) continue;
1432 Node bnode = graph->bNode(jt);
1433 Value bvalue = avalue + (*cost)[jt] +
1434 anode_potential[anode] - bnode_potential[bnode];
1435 if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1436 bdist[bnode] = bvalue;
1439 if (bvalue > bdistMax) {
1442 if (bnode_matching[bnode] != INVALID) {
1443 Node newanode = graph->aNode(bnode_matching[bnode]);
1444 switch (_heap->state(newanode)) {
1445 case Heap::PRE_HEAP:
1446 _heap->push(newanode, bvalue);
1449 if (bvalue < (*_heap)[newanode]) {
1450 _heap->decrease(newanode, bvalue);
1453 case Heap::POST_HEAP:
1457 if (bestNode == INVALID ||
1458 bvalue + bnode_potential[bnode] < bestValue) {
1459 bestValue = bvalue + bnode_potential[bnode];
1466 if (bestNode == INVALID) {
1470 matching_cost += bestValue;
1473 for (BNodeIt it(*graph); it != INVALID; ++it) {
1474 if (bpred[it] != INVALID) {
1475 bnode_potential[it] += bdist[it];
1477 bnode_potential[it] += bdistMax;
1480 for (ANodeIt it(*graph); it != INVALID; ++it) {
1481 if (anode_matching[it] != INVALID) {
1482 Node bnode = graph->bNode(anode_matching[it]);
1483 if (bpred[bnode] != INVALID) {
1484 anode_potential[it] += bdist[bnode];
1486 anode_potential[it] += bdistMax;
1491 while (bestNode != INVALID) {
1492 UEdge uedge = bpred[bestNode];
1493 Node anode = graph->aNode(uedge);
1495 bnode_matching[bestNode] = uedge;
1496 if (anode_matching[anode] != INVALID) {
1497 bestNode = graph->bNode(anode_matching[anode]);
1501 anode_matching[anode] = uedge;
1508 /// \brief Starts the algorithm.
1510 /// Starts the algorithm. It runs augmenting phases until the
1511 /// optimal solution reached.
1513 while (augment()) {}
1516 /// \brief Runs the algorithm.
1518 /// It just initalize the algorithm and then start it.
1526 /// \name Query Functions
1527 /// The result of the %Matching algorithm can be obtained using these
1529 /// Before the use of these functions,
1530 /// either run() or start() must be called.
1534 /// \brief Gives back the potential in the NodeMap
1536 /// Gives back the potential in the NodeMap. The potential is optimal with
1537 /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
1538 /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
1540 template <typename PotentialMap>
1541 void potential(PotentialMap& pt) const {
1542 for (ANodeIt it(*graph); it != INVALID; ++it) {
1543 pt[it] = anode_potential[it];
1545 for (BNodeIt it(*graph); it != INVALID; ++it) {
1546 pt[it] = bnode_potential[it];
1550 /// \brief Set true all matching uedge in the map.
1552 /// Set true all matching uedge in the map. It does not change the
1553 /// value mapped to the other uedges.
1554 /// \return The number of the matching edges.
1555 template <typename MatchingMap>
1556 int quickMatching(MatchingMap& mm) const {
1557 for (ANodeIt it(*graph); it != INVALID; ++it) {
1558 if (anode_matching[it] != INVALID) {
1559 mm[anode_matching[it]] = true;
1562 return matching_size;
1565 /// \brief Set true all matching uedge in the map and the others to false.
1567 /// Set true all matching uedge in the map and the others to false.
1568 /// \return The number of the matching edges.
1569 template <typename MatchingMap>
1570 int matching(MatchingMap& mm) const {
1571 for (UEdgeIt it(*graph); it != INVALID; ++it) {
1572 mm[it] = it == anode_matching[graph->aNode(it)];
1574 return matching_size;
1578 /// \brief Return true if the given uedge is in the matching.
1580 /// It returns true if the given uedge is in the matching.
1581 bool matchingEdge(const UEdge& edge) const {
1582 return anode_matching[graph->aNode(edge)] == edge;
1585 /// \brief Returns the matching edge from the node.
1587 /// Returns the matching edge from the node. If there is not such
1588 /// edge it gives back \c INVALID.
1589 UEdge matchingEdge(const Node& node) const {
1590 if (graph->aNode(node)) {
1591 return anode_matching[node];
1593 return bnode_matching[node];
1597 /// \brief Gives back the sum of costs of the matching edges.
1599 /// Gives back the sum of costs of the matching edges.
1600 Value matchingCost() const {
1601 return matching_cost;
1604 /// \brief Gives back the number of the matching edges.
1606 /// Gives back the number of the matching edges.
1607 int matchingSize() const {
1608 return matching_size;
1615 void initStructures() {
1616 if (!_heap_cross_ref) {
1617 local_heap_cross_ref = true;
1618 _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1622 _heap = Traits::createHeap(*_heap_cross_ref);
1626 void destroyStructures() {
1627 if (local_heap_cross_ref) delete _heap_cross_ref;
1628 if (local_heap) delete _heap;
1634 const BpUGraph *graph;
1635 const CostMap* cost;
1637 ANodeMatchingMap anode_matching;
1638 BNodeMatchingMap bnode_matching;
1640 ANodePotentialMap anode_potential;
1641 BNodePotentialMap bnode_potential;
1643 Value matching_cost;
1646 HeapCrossRef *_heap_cross_ref;
1647 bool local_heap_cross_ref;
1654 /// \ingroup matching
1656 /// \brief Minimum cost maximum cardinality bipartite matching
1658 /// This function calculates the minimum cost matching of the maximum
1659 /// cardinality matchings of a bipartite graph. It gives back the matching
1660 /// in an undirected edge map.
1662 /// \param graph The bipartite graph.
1663 /// \param cost The undirected edge map which contains the costs.
1664 /// \retval matching The undirected edge map which will be set to
1666 /// \return The cost of the matching.
1667 template <typename BpUGraph, typename CostMap, typename MatchingMap>
1668 typename CostMap::Value
1669 minCostMaxBipartiteMatching(const BpUGraph& graph,
1670 const CostMap& cost,
1671 MatchingMap& matching) {
1672 MinCostMaxBipartiteMatching<BpUGraph, CostMap>
1673 bpmatching(graph, cost);
1675 bpmatching.matching(matching);
1676 return bpmatching.matchingCost();