Interface to the cplex MIP solver: it is little, a bit sour but it is ours.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2006
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.
37 /// \brief Bipartite Max Cardinality Matching algorithm
39 /// Bipartite Max Cardinality Matching algorithm. This class implements
40 /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
42 template <typename BpUGraph>
43 class MaxBipartiteMatching {
46 typedef BpUGraph Graph;
48 typedef typename Graph::Node Node;
49 typedef typename Graph::ANodeIt ANodeIt;
50 typedef typename Graph::BNodeIt BNodeIt;
51 typedef typename Graph::UEdge UEdge;
52 typedef typename Graph::UEdgeIt UEdgeIt;
53 typedef typename Graph::IncEdgeIt IncEdgeIt;
55 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
56 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
61 /// \brief Constructor.
63 /// Constructor of the algorithm.
64 MaxBipartiteMatching(const BpUGraph& _graph)
65 : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
67 /// \name Execution control
68 /// The simplest way to execute the algorithm is to use
69 /// one of the member functions called \c run().
71 /// If you need more control on the execution,
72 /// first you must call \ref init() or one alternative for it.
73 /// Finally \ref start() will perform the matching computation or
74 /// with step-by-step execution you can augment the solution.
78 /// \brief Initalize the data structures.
80 /// It initalizes the data structures and creates an empty matching.
82 for (ANodeIt it(*graph); it != INVALID; ++it) {
83 anode_matching[it] = INVALID;
85 for (BNodeIt it(*graph); it != INVALID; ++it) {
86 bnode_matching[it] = INVALID;
91 /// \brief Initalize the data structures.
93 /// It initalizes the data structures and creates a greedy
94 /// matching. From this matching sometimes it is faster to get
95 /// the matching than from the initial empty matching.
98 for (BNodeIt it(*graph); it != INVALID; ++it) {
99 bnode_matching[it] = INVALID;
101 for (ANodeIt it(*graph); it != INVALID; ++it) {
102 anode_matching[it] = INVALID;
103 for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
104 if (bnode_matching[graph->bNode(jt)] == INVALID) {
105 anode_matching[it] = jt;
106 bnode_matching[graph->bNode(jt)] = jt;
114 /// \brief Initalize the data structures with an initial matching.
116 /// It initalizes the data structures with an initial matching.
117 template <typename MatchingMap>
118 void matchingInit(const MatchingMap& matching) {
119 for (ANodeIt it(*graph); it != INVALID; ++it) {
120 anode_matching[it] = INVALID;
122 for (BNodeIt it(*graph); it != INVALID; ++it) {
123 bnode_matching[it] = INVALID;
126 for (UEdgeIt it(*graph); it != INVALID; ++it) {
129 anode_matching[graph->aNode(it)] = it;
130 bnode_matching[graph->bNode(it)] = it;
135 /// \brief Initalize the data structures with an initial matching.
137 /// It initalizes the data structures with an initial matching.
138 /// \return %True when the given map contains really a matching.
139 template <typename MatchingMap>
140 void checkedMatchingInit(const MatchingMap& matching) {
141 for (ANodeIt it(*graph); it != INVALID; ++it) {
142 anode_matching[it] = INVALID;
144 for (BNodeIt it(*graph); it != INVALID; ++it) {
145 bnode_matching[it] = INVALID;
148 for (UEdgeIt it(*graph); it != INVALID; ++it) {
151 if (anode_matching[graph->aNode(it)] != INVALID) {
154 anode_matching[graph->aNode(it)] = it;
155 if (bnode_matching[graph->aNode(it)] != INVALID) {
158 bnode_matching[graph->bNode(it)] = it;
164 /// \brief An augmenting phase of the Hopcroft-Karp algorithm
166 /// It runs an augmenting phase of the Hopcroft-Karp
167 /// algorithm. The phase finds maximum count of edge disjoint
168 /// augmenting paths and augments on these paths. The algorithm
169 /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
170 /// \f$ O(e) \f$ long.
173 typename Graph::template ANodeMap<bool> areached(*graph, false);
174 typename Graph::template BNodeMap<bool> breached(*graph, false);
176 typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
178 std::vector<Node> queue, bqueue;
179 for (ANodeIt it(*graph); it != INVALID; ++it) {
180 if (anode_matching[it] == INVALID) {
186 bool success = false;
188 while (!success && !queue.empty()) {
189 std::vector<Node> newqueue;
190 for (int i = 0; i < (int)queue.size(); ++i) {
191 Node anode = queue[i];
192 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
193 Node bnode = graph->bNode(jt);
194 if (breached[bnode]) continue;
195 breached[bnode] = true;
197 if (bnode_matching[bnode] == INVALID) {
198 bqueue.push_back(bnode);
201 Node newanode = graph->aNode(bnode_matching[bnode]);
202 if (!areached[newanode]) {
203 areached[newanode] = true;
204 newqueue.push_back(newanode);
209 queue.swap(newqueue);
214 typename Graph::template ANodeMap<bool> aused(*graph, false);
216 for (int i = 0; i < (int)bqueue.size(); ++i) {
217 Node bnode = bqueue[i];
221 while (bnode != INVALID) {
222 UEdge uedge = bpred[bnode];
223 Node anode = graph->aNode(uedge);
230 bnode = anode_matching[anode] != INVALID ?
231 graph->bNode(anode_matching[anode]) : INVALID;
238 while (bnode != INVALID) {
239 UEdge uedge = bpred[bnode];
240 Node anode = graph->aNode(uedge);
242 bnode_matching[bnode] = uedge;
244 bnode = anode_matching[anode] != INVALID ?
245 graph->bNode(anode_matching[anode]) : INVALID;
247 anode_matching[anode] = uedge;
258 /// \brief An augmenting phase of the Ford-Fulkerson algorithm
260 /// It runs an augmenting phase of the Ford-Fulkerson
261 /// algorithm. The phase finds only one augmenting path and
262 /// augments only on this paths. The algorithm consists at most
263 /// of \f$ O(n) \f$ simple phase and one phase is at most
264 /// \f$ O(e) \f$ long.
265 bool simpleAugment() {
267 typename Graph::template ANodeMap<bool> areached(*graph, false);
268 typename Graph::template BNodeMap<bool> breached(*graph, false);
270 typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
272 std::vector<Node> queue;
273 for (ANodeIt it(*graph); it != INVALID; ++it) {
274 if (anode_matching[it] == INVALID) {
280 while (!queue.empty()) {
281 std::vector<Node> newqueue;
282 for (int i = 0; i < (int)queue.size(); ++i) {
283 Node anode = queue[i];
284 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
285 Node bnode = graph->bNode(jt);
286 if (breached[bnode]) continue;
287 breached[bnode] = true;
289 if (bnode_matching[bnode] == INVALID) {
290 while (bnode != INVALID) {
291 UEdge uedge = bpred[bnode];
292 anode = graph->aNode(uedge);
294 bnode_matching[bnode] = uedge;
296 bnode = anode_matching[anode] != INVALID ?
297 graph->bNode(anode_matching[anode]) : INVALID;
299 anode_matching[anode] = uedge;
305 Node newanode = graph->aNode(bnode_matching[bnode]);
306 if (!areached[newanode]) {
307 areached[newanode] = true;
308 newqueue.push_back(newanode);
313 queue.swap(newqueue);
319 /// \brief Starts the algorithm.
321 /// Starts the algorithm. It runs augmenting phases until the optimal
322 /// solution reached.
327 /// \brief Runs the algorithm.
329 /// It just initalize the algorithm and then start it.
337 /// \name Query Functions
338 /// The result of the %Matching algorithm can be obtained using these
340 /// Before the use of these functions,
341 /// either run() or start() must be called.
345 /// \brief Returns an minimum covering of the nodes.
347 /// The minimum covering set problem is the dual solution of the
348 /// maximum bipartite matching. It provides an solution for this
349 /// problem what is proof of the optimality of the matching.
350 /// \return The size of the cover set.
351 template <typename CoverMap>
352 int coverSet(CoverMap& covering) const {
354 typename Graph::template ANodeMap<bool> areached(*graph, false);
355 typename Graph::template BNodeMap<bool> breached(*graph, false);
357 std::vector<Node> queue;
358 for (ANodeIt it(*graph); it != INVALID; ++it) {
359 if (anode_matching[it] == INVALID) {
364 while (!queue.empty()) {
365 std::vector<Node> newqueue;
366 for (int i = 0; i < (int)queue.size(); ++i) {
367 Node anode = queue[i];
368 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
369 Node bnode = graph->bNode(jt);
370 if (breached[bnode]) continue;
371 breached[bnode] = true;
372 if (bnode_matching[bnode] != INVALID) {
373 Node newanode = graph->aNode(bnode_matching[bnode]);
374 if (!areached[newanode]) {
375 areached[newanode] = true;
376 newqueue.push_back(newanode);
381 queue.swap(newqueue);
385 for (ANodeIt it(*graph); it != INVALID; ++it) {
386 covering[it] = !areached[it] && anode_matching[it] != INVALID;
387 if (!areached[it] && anode_matching[it] != INVALID) {
391 for (BNodeIt it(*graph); it != INVALID; ++it) {
392 covering[it] = breached[it];
400 /// \brief Set true all matching uedge in the map.
402 /// Set true all matching uedge in the map. It does not change the
403 /// value mapped to the other uedges.
404 /// \return The number of the matching edges.
405 template <typename MatchingMap>
406 int quickMatching(MatchingMap& matching) const {
407 for (ANodeIt it(*graph); it != INVALID; ++it) {
408 if (anode_matching[it] != INVALID) {
409 matching[anode_matching[it]] = true;
412 return matching_size;
415 /// \brief Set true all matching uedge in the map and the others to false.
417 /// Set true all matching uedge in the map and the others to false.
418 /// \return The number of the matching edges.
419 template <typename MatchingMap>
420 int matching(MatchingMap& matching) const {
421 for (UEdgeIt it(*graph); it != INVALID; ++it) {
422 matching[it] = it == anode_matching[graph->aNode(it)];
424 return matching_size;
428 /// \brief Return true if the given uedge is in the matching.
430 /// It returns true if the given uedge is in the matching.
431 bool matchingEdge(const UEdge& edge) const {
432 return anode_matching[graph->aNode(edge)] == edge;
435 /// \brief Returns the matching edge from the node.
437 /// Returns the matching edge from the node. If there is not such
438 /// edge it gives back \c INVALID.
439 UEdge matchingEdge(const Node& node) const {
440 if (graph->aNode(node)) {
441 return anode_matching[node];
443 return bnode_matching[node];
447 /// \brief Gives back the number of the matching edges.
449 /// Gives back the number of the matching edges.
450 int matchingSize() const {
451 return matching_size;
458 ANodeMatchingMap anode_matching;
459 BNodeMatchingMap bnode_matching;
466 /// \ingroup matching
468 /// \brief Maximum cardinality bipartite matching
470 /// This function calculates the maximum cardinality matching
471 /// in a bipartite graph. It gives back the matching in an undirected
474 /// \param graph The bipartite graph.
475 /// \retval matching The undirected edge map which will be set to
477 /// \return The size of the matching.
478 template <typename BpUGraph, typename MatchingMap>
479 int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
480 MaxBipartiteMatching<BpUGraph> bpmatching(graph);
482 bpmatching.matching(matching);
483 return bpmatching.matchingSize();
486 /// \brief Default traits class for weighted bipartite matching algoritms.
488 /// Default traits class for weighted bipartite matching algoritms.
489 /// \param _BpUGraph The bipartite undirected graph type.
490 /// \param _WeightMap Type of weight map.
491 template <typename _BpUGraph, typename _WeightMap>
492 struct WeightedBipartiteMatchingDefaultTraits {
493 /// \brief The type of the weight of the undirected edges.
494 typedef typename _WeightMap::Value Value;
496 /// The undirected bipartite graph type the algorithm runs on.
497 typedef _BpUGraph BpUGraph;
499 /// The map of the edges weights
500 typedef _WeightMap WeightMap;
502 /// \brief The cross reference type used by heap.
504 /// The cross reference type used by heap.
505 /// Usually it is \c Graph::NodeMap<int>.
506 typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
508 /// \brief Instantiates a HeapCrossRef.
510 /// This function instantiates a \ref HeapCrossRef.
511 /// \param graph is the graph, to which we would like to define the
513 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
514 return new HeapCrossRef(graph);
517 /// \brief The heap type used by weighted matching algorithms.
519 /// The heap type used by weighted matching algorithms. It should
520 /// minimize the priorities and the heap's key type is the graph's
521 /// anode graph's node.
524 typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
526 /// \brief Instantiates a Heap.
528 /// This function instantiates a \ref Heap.
529 /// \param crossref The cross reference of the heap.
530 static Heap *createHeap(HeapCrossRef& crossref) {
531 return new Heap(crossref);
537 /// \ingroup matching
539 /// \brief Bipartite Max Weighted Matching algorithm
541 /// This class implements the bipartite Max Weighted Matching
542 /// algorithm. It uses the successive shortest path algorithm to
543 /// calculate the maximum weighted matching in the bipartite
544 /// graph. The algorithm can be used also to calculate the maximum
545 /// cardinality maximum weighted matching. The time complexity
546 /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
547 /// heap implementation but this can be improved to
548 /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
550 /// The algorithm also provides a potential function on the nodes
551 /// which a dual solution of the matching algorithm and it can be
552 /// used to proof the optimality of the given pimal solution.
554 template <typename _BpUGraph, typename _WeightMap, typename _Traits>
556 template <typename _BpUGraph,
557 typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
558 typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
560 class MaxWeightedBipartiteMatching {
563 typedef _Traits Traits;
564 typedef typename Traits::BpUGraph BpUGraph;
565 typedef typename Traits::WeightMap WeightMap;
566 typedef typename Traits::Value Value;
570 typedef typename Traits::HeapCrossRef HeapCrossRef;
571 typedef typename Traits::Heap Heap;
574 typedef typename BpUGraph::Node Node;
575 typedef typename BpUGraph::ANodeIt ANodeIt;
576 typedef typename BpUGraph::BNodeIt BNodeIt;
577 typedef typename BpUGraph::UEdge UEdge;
578 typedef typename BpUGraph::UEdgeIt UEdgeIt;
579 typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
581 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
582 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
584 typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
585 typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
590 /// \brief \ref Exception for uninitialized parameters.
592 /// This error represents problems in the initialization
593 /// of the parameters of the algorithms.
594 class UninitializedParameter : public lemon::UninitializedParameter {
596 virtual const char* what() const throw() {
597 return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
601 ///\name Named template parameters
605 template <class H, class CR>
606 struct DefHeapTraits : public Traits {
607 typedef CR HeapCrossRef;
609 static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
610 throw UninitializedParameter();
612 static Heap *createHeap(HeapCrossRef &) {
613 throw UninitializedParameter();
617 /// \brief \ref named-templ-param "Named parameter" for setting heap
618 /// and cross reference type
620 /// \ref named-templ-param "Named parameter" for setting heap and cross
622 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
624 : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
625 DefHeapTraits<H, CR> > {
626 typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
627 DefHeapTraits<H, CR> > Create;
630 template <class H, class CR>
631 struct DefStandardHeapTraits : public Traits {
632 typedef CR HeapCrossRef;
634 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
635 return new HeapCrossRef(graph);
637 static Heap *createHeap(HeapCrossRef &crossref) {
638 return new Heap(crossref);
642 /// \brief \ref named-templ-param "Named parameter" for setting heap and
643 /// cross reference type with automatic allocation
645 /// \ref named-templ-param "Named parameter" for setting heap and cross
646 /// reference type. It can allocate the heap and the cross reference
647 /// object if the cross reference's constructor waits for the graph as
648 /// parameter and the heap's constructor waits for the cross reference.
649 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
650 struct DefStandardHeap
651 : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
652 DefStandardHeapTraits<H, CR> > {
653 typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
654 DefStandardHeapTraits<H, CR> >
661 /// \brief Constructor.
663 /// Constructor of the algorithm.
664 MaxWeightedBipartiteMatching(const BpUGraph& _graph,
665 const WeightMap& _weight)
666 : graph(&_graph), weight(&_weight),
667 anode_matching(_graph), bnode_matching(_graph),
668 anode_potential(_graph), bnode_potential(_graph),
669 _heap_cross_ref(0), local_heap_cross_ref(false),
670 _heap(0), local_heap(0) {}
672 /// \brief Destructor.
674 /// Destructor of the algorithm.
675 ~MaxWeightedBipartiteMatching() {
679 /// \brief Sets the heap and the cross reference used by algorithm.
681 /// Sets the heap and the cross reference used by algorithm.
682 /// If you don't use this function before calling \ref run(),
683 /// it will allocate one. The destuctor deallocates this
684 /// automatically allocated map, of course.
685 /// \return \c (*this)
686 MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
687 if(local_heap_cross_ref) {
688 delete _heap_cross_ref;
689 local_heap_cross_ref = false;
691 _heap_cross_ref = &crossRef;
700 /// \name Execution control
701 /// The simplest way to execute the algorithm is to use
702 /// one of the member functions called \c run().
704 /// If you need more control on the execution,
705 /// first you must call \ref init() or one alternative for it.
706 /// Finally \ref start() will perform the matching computation or
707 /// with step-by-step execution you can augment the solution.
711 /// \brief Initalize the data structures.
713 /// It initalizes the data structures and creates an empty matching.
716 for (ANodeIt it(*graph); it != INVALID; ++it) {
717 anode_matching[it] = INVALID;
718 anode_potential[it] = 0;
720 for (BNodeIt it(*graph); it != INVALID; ++it) {
721 bnode_matching[it] = INVALID;
722 bnode_potential[it] = 0;
723 for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
724 if ((*weight)[jt] > bnode_potential[it]) {
725 bnode_potential[it] = (*weight)[jt];
734 /// \brief An augmenting phase of the weighted matching algorithm
736 /// It runs an augmenting phase of the weighted matching
737 /// algorithm. The phase finds the best augmenting path and
738 /// augments only on this paths.
740 /// The algorithm consists at most
741 /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
742 /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
743 /// with binary heap.
744 /// \param decrease If the given parameter true the matching value
745 /// can be decreased in the augmenting phase. If we would like
746 /// to calculate the maximum cardinality maximum weighted matching
747 /// then we should let the algorithm to decrease the matching
748 /// value in order to increase the number of the matching edges.
749 bool augment(bool decrease = false) {
751 typename BpUGraph::template BNodeMap<Value> bdist(*graph);
752 typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
754 Node bestNode = INVALID;
758 for (ANodeIt it(*graph); it != INVALID; ++it) {
759 (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
762 for (ANodeIt it(*graph); it != INVALID; ++it) {
763 if (anode_matching[it] == INVALID) {
769 while (!_heap->empty()) {
770 Node anode = _heap->top();
771 Value avalue = _heap->prio();
773 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
774 if (jt == anode_matching[anode]) continue;
775 Node bnode = graph->bNode(jt);
776 Value bvalue = avalue - (*weight)[jt] +
777 anode_potential[anode] + bnode_potential[bnode];
778 if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
779 bdist[bnode] = bvalue;
782 if (bvalue > bdistMax) {
785 if (bnode_matching[bnode] != INVALID) {
786 Node newanode = graph->aNode(bnode_matching[bnode]);
787 switch (_heap->state(newanode)) {
789 _heap->push(newanode, bvalue);
792 if (bvalue < (*_heap)[newanode]) {
793 _heap->decrease(newanode, bvalue);
796 case Heap::POST_HEAP:
800 if (bestNode == INVALID ||
801 bnode_potential[bnode] - bvalue > bestValue) {
802 bestValue = bnode_potential[bnode] - bvalue;
809 if (bestNode == INVALID || (!decrease && bestValue < 0)) {
813 matching_value += bestValue;
816 for (BNodeIt it(*graph); it != INVALID; ++it) {
817 if (bpred[it] != INVALID) {
818 bnode_potential[it] -= bdist[it];
820 bnode_potential[it] -= bdistMax;
823 for (ANodeIt it(*graph); it != INVALID; ++it) {
824 if (anode_matching[it] != INVALID) {
825 Node bnode = graph->bNode(anode_matching[it]);
826 if (bpred[bnode] != INVALID) {
827 anode_potential[it] += bdist[bnode];
829 anode_potential[it] += bdistMax;
834 while (bestNode != INVALID) {
835 UEdge uedge = bpred[bestNode];
836 Node anode = graph->aNode(uedge);
838 bnode_matching[bestNode] = uedge;
839 if (anode_matching[anode] != INVALID) {
840 bestNode = graph->bNode(anode_matching[anode]);
844 anode_matching[anode] = uedge;
851 /// \brief Starts the algorithm.
853 /// Starts the algorithm. It runs augmenting phases until the
854 /// optimal solution reached.
856 /// \param maxCardinality If the given value is true it will
857 /// calculate the maximum cardinality maximum matching instead of
858 /// the maximum matching.
859 void start(bool maxCardinality = false) {
860 while (augment(maxCardinality)) {}
863 /// \brief Runs the algorithm.
865 /// It just initalize the algorithm and then start it.
867 /// \param maxCardinality If the given value is true it will
868 /// calculate the maximum cardinality maximum matching instead of
869 /// the maximum matching.
870 void run(bool maxCardinality = false) {
872 start(maxCardinality);
877 /// \name Query Functions
878 /// The result of the %Matching algorithm can be obtained using these
880 /// Before the use of these functions,
881 /// either run() or start() must be called.
885 /// \brief Gives back the potential in the NodeMap
887 /// Gives back the potential in the NodeMap. The matching is optimal
888 /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
889 /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
891 template <typename PotentialMap>
892 void potential(PotentialMap& potential) const {
893 for (ANodeIt it(*graph); it != INVALID; ++it) {
894 potential[it] = anode_potential[it];
896 for (BNodeIt it(*graph); it != INVALID; ++it) {
897 potential[it] = bnode_potential[it];
901 /// \brief Set true all matching uedge in the map.
903 /// Set true all matching uedge in the map. It does not change the
904 /// value mapped to the other uedges.
905 /// \return The number of the matching edges.
906 template <typename MatchingMap>
907 int quickMatching(MatchingMap& matching) const {
908 for (ANodeIt it(*graph); it != INVALID; ++it) {
909 if (anode_matching[it] != INVALID) {
910 matching[anode_matching[it]] = true;
913 return matching_size;
916 /// \brief Set true all matching uedge in the map and the others to false.
918 /// Set true all matching uedge in the map and the others to false.
919 /// \return The number of the matching edges.
920 template <typename MatchingMap>
921 int matching(MatchingMap& matching) const {
922 for (UEdgeIt it(*graph); it != INVALID; ++it) {
923 matching[it] = it == anode_matching[graph->aNode(it)];
925 return matching_size;
929 /// \brief Return true if the given uedge is in the matching.
931 /// It returns true if the given uedge is in the matching.
932 bool matchingEdge(const UEdge& edge) const {
933 return anode_matching[graph->aNode(edge)] == edge;
936 /// \brief Returns the matching edge from the node.
938 /// Returns the matching edge from the node. If there is not such
939 /// edge it gives back \c INVALID.
940 UEdge matchingEdge(const Node& node) const {
941 if (graph->aNode(node)) {
942 return anode_matching[node];
944 return bnode_matching[node];
948 /// \brief Gives back the sum of weights of the matching edges.
950 /// Gives back the sum of weights of the matching edges.
951 Value matchingValue() const {
952 return matching_value;
955 /// \brief Gives back the number of the matching edges.
957 /// Gives back the number of the matching edges.
958 int matchingSize() const {
959 return matching_size;
966 void initStructures() {
967 if (!_heap_cross_ref) {
968 local_heap_cross_ref = true;
969 _heap_cross_ref = Traits::createHeapCrossRef(*graph);
973 _heap = Traits::createHeap(*_heap_cross_ref);
977 void destroyStructures() {
978 if (local_heap_cross_ref) delete _heap_cross_ref;
979 if (local_heap) delete _heap;
985 const BpUGraph *graph;
986 const WeightMap* weight;
988 ANodeMatchingMap anode_matching;
989 BNodeMatchingMap bnode_matching;
991 ANodePotentialMap anode_potential;
992 BNodePotentialMap bnode_potential;
994 Value matching_value;
997 HeapCrossRef *_heap_cross_ref;
998 bool local_heap_cross_ref;
1005 /// \ingroup matching
1007 /// \brief Maximum weighted bipartite matching
1009 /// This function calculates the maximum weighted matching
1010 /// in a bipartite graph. It gives back the matching in an undirected
1013 /// \param graph The bipartite graph.
1014 /// \param weight The undirected edge map which contains the weights.
1015 /// \retval matching The undirected edge map which will be set to
1017 /// \return The value of the matching.
1018 template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1019 typename WeightMap::Value
1020 maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
1021 MatchingMap& matching) {
1022 MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1023 bpmatching(graph, weight);
1025 bpmatching.matching(matching);
1026 return bpmatching.matchingValue();
1029 /// \ingroup matching
1031 /// \brief Maximum weighted maximum cardinality bipartite matching
1033 /// This function calculates the maximum weighted of the maximum cardinality
1034 /// matchings of a bipartite graph. It gives back the matching in an
1035 /// undirected edge map.
1037 /// \param graph The bipartite graph.
1038 /// \param weight The undirected edge map which contains the weights.
1039 /// \retval matching The undirected edge map which will be set to
1041 /// \return The value of the matching.
1042 template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1043 typename WeightMap::Value
1044 maxWeightedMaxBipartiteMatching(const BpUGraph& graph,
1045 const WeightMap& weight,
1046 MatchingMap& matching) {
1047 MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1048 bpmatching(graph, weight);
1049 bpmatching.run(true);
1050 bpmatching.matching(matching);
1051 return bpmatching.matchingValue();
1054 /// \brief Default traits class for minimum cost bipartite matching
1057 /// Default traits class for minimum cost bipartite matching
1060 /// \param _BpUGraph The bipartite undirected graph
1063 /// \param _CostMap Type of cost map.
1064 template <typename _BpUGraph, typename _CostMap>
1065 struct MinCostMaxBipartiteMatchingDefaultTraits {
1066 /// \brief The type of the cost of the undirected edges.
1067 typedef typename _CostMap::Value Value;
1069 /// The undirected bipartite graph type the algorithm runs on.
1070 typedef _BpUGraph BpUGraph;
1072 /// The map of the edges costs
1073 typedef _CostMap CostMap;
1075 /// \brief The cross reference type used by heap.
1077 /// The cross reference type used by heap.
1078 /// Usually it is \c Graph::NodeMap<int>.
1079 typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1081 /// \brief Instantiates a HeapCrossRef.
1083 /// This function instantiates a \ref HeapCrossRef.
1084 /// \param graph is the graph, to which we would like to define the
1086 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1087 return new HeapCrossRef(graph);
1090 /// \brief The heap type used by costed matching algorithms.
1092 /// The heap type used by costed matching algorithms. It should
1093 /// minimize the priorities and the heap's key type is the graph's
1094 /// anode graph's node.
1097 typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
1099 /// \brief Instantiates a Heap.
1101 /// This function instantiates a \ref Heap.
1102 /// \param crossref The cross reference of the heap.
1103 static Heap *createHeap(HeapCrossRef& crossref) {
1104 return new Heap(crossref);
1110 /// \ingroup matching
1112 /// \brief Bipartite Min Cost Matching algorithm
1114 /// This class implements the bipartite Min Cost Matching algorithm.
1115 /// It uses the successive shortest path algorithm to calculate the
1116 /// minimum cost maximum matching in the bipartite graph. The time
1117 /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
1118 /// default binary heap implementation but this can be improved to
1119 /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1121 /// The algorithm also provides a potential function on the nodes
1122 /// which a dual solution of the matching algorithm and it can be
1123 /// used to proof the optimality of the given pimal solution.
1125 template <typename _BpUGraph, typename _CostMap, typename _Traits>
1127 template <typename _BpUGraph,
1128 typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
1129 typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
1131 class MinCostMaxBipartiteMatching {
1134 typedef _Traits Traits;
1135 typedef typename Traits::BpUGraph BpUGraph;
1136 typedef typename Traits::CostMap CostMap;
1137 typedef typename Traits::Value Value;
1141 typedef typename Traits::HeapCrossRef HeapCrossRef;
1142 typedef typename Traits::Heap Heap;
1145 typedef typename BpUGraph::Node Node;
1146 typedef typename BpUGraph::ANodeIt ANodeIt;
1147 typedef typename BpUGraph::BNodeIt BNodeIt;
1148 typedef typename BpUGraph::UEdge UEdge;
1149 typedef typename BpUGraph::UEdgeIt UEdgeIt;
1150 typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1152 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1153 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1155 typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1156 typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1161 /// \brief \ref Exception for uninitialized parameters.
1163 /// This error represents problems in the initialization
1164 /// of the parameters of the algorithms.
1165 class UninitializedParameter : public lemon::UninitializedParameter {
1167 virtual const char* what() const throw() {
1168 return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
1172 ///\name Named template parameters
1176 template <class H, class CR>
1177 struct DefHeapTraits : public Traits {
1178 typedef CR HeapCrossRef;
1180 static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1181 throw UninitializedParameter();
1183 static Heap *createHeap(HeapCrossRef &) {
1184 throw UninitializedParameter();
1188 /// \brief \ref named-templ-param "Named parameter" for setting heap
1189 /// and cross reference type
1191 /// \ref named-templ-param "Named parameter" for setting heap and cross
1193 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1195 : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1196 DefHeapTraits<H, CR> > {
1197 typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1198 DefHeapTraits<H, CR> > Create;
1201 template <class H, class CR>
1202 struct DefStandardHeapTraits : public Traits {
1203 typedef CR HeapCrossRef;
1205 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1206 return new HeapCrossRef(graph);
1208 static Heap *createHeap(HeapCrossRef &crossref) {
1209 return new Heap(crossref);
1213 /// \brief \ref named-templ-param "Named parameter" for setting heap and
1214 /// cross reference type with automatic allocation
1216 /// \ref named-templ-param "Named parameter" for setting heap and cross
1217 /// reference type. It can allocate the heap and the cross reference
1218 /// object if the cross reference's constructor waits for the graph as
1219 /// parameter and the heap's constructor waits for the cross reference.
1220 template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1221 struct DefStandardHeap
1222 : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1223 DefStandardHeapTraits<H, CR> > {
1224 typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1225 DefStandardHeapTraits<H, CR> >
1232 /// \brief Constructor.
1234 /// Constructor of the algorithm.
1235 MinCostMaxBipartiteMatching(const BpUGraph& _graph,
1236 const CostMap& _cost)
1237 : graph(&_graph), cost(&_cost),
1238 anode_matching(_graph), bnode_matching(_graph),
1239 anode_potential(_graph), bnode_potential(_graph),
1240 _heap_cross_ref(0), local_heap_cross_ref(false),
1241 _heap(0), local_heap(0) {}
1243 /// \brief Destructor.
1245 /// Destructor of the algorithm.
1246 ~MinCostMaxBipartiteMatching() {
1247 destroyStructures();
1250 /// \brief Sets the heap and the cross reference used by algorithm.
1252 /// Sets the heap and the cross reference used by algorithm.
1253 /// If you don't use this function before calling \ref run(),
1254 /// it will allocate one. The destuctor deallocates this
1255 /// automatically allocated map, of course.
1256 /// \return \c (*this)
1257 MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
1258 if(local_heap_cross_ref) {
1259 delete _heap_cross_ref;
1260 local_heap_cross_ref = false;
1262 _heap_cross_ref = &crossRef;
1271 /// \name Execution control
1272 /// The simplest way to execute the algorithm is to use
1273 /// one of the member functions called \c run().
1275 /// If you need more control on the execution,
1276 /// first you must call \ref init() or one alternative for it.
1277 /// Finally \ref start() will perform the matching computation or
1278 /// with step-by-step execution you can augment the solution.
1282 /// \brief Initalize the data structures.
1284 /// It initalizes the data structures and creates an empty matching.
1287 for (ANodeIt it(*graph); it != INVALID; ++it) {
1288 anode_matching[it] = INVALID;
1289 anode_potential[it] = 0;
1291 for (BNodeIt it(*graph); it != INVALID; ++it) {
1292 bnode_matching[it] = INVALID;
1293 bnode_potential[it] = 0;
1300 /// \brief An augmenting phase of the costed matching algorithm
1302 /// It runs an augmenting phase of the matching algorithm. The
1303 /// phase finds the best augmenting path and augments only on this
1306 /// The algorithm consists at most
1307 /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1308 /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1309 /// with binary heap.
1312 typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1313 typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1315 Node bestNode = INVALID;
1316 Value bestValue = 0;
1319 for (ANodeIt it(*graph); it != INVALID; ++it) {
1320 (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1323 for (ANodeIt it(*graph); it != INVALID; ++it) {
1324 if (anode_matching[it] == INVALID) {
1330 while (!_heap->empty()) {
1331 Node anode = _heap->top();
1332 Value avalue = _heap->prio();
1334 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1335 if (jt == anode_matching[anode]) continue;
1336 Node bnode = graph->bNode(jt);
1337 Value bvalue = avalue + (*cost)[jt] +
1338 anode_potential[anode] - bnode_potential[bnode];
1339 if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1340 bdist[bnode] = bvalue;
1343 if (bvalue > bdistMax) {
1346 if (bnode_matching[bnode] != INVALID) {
1347 Node newanode = graph->aNode(bnode_matching[bnode]);
1348 switch (_heap->state(newanode)) {
1349 case Heap::PRE_HEAP:
1350 _heap->push(newanode, bvalue);
1353 if (bvalue < (*_heap)[newanode]) {
1354 _heap->decrease(newanode, bvalue);
1357 case Heap::POST_HEAP:
1361 if (bestNode == INVALID ||
1362 bvalue + bnode_potential[bnode] < bestValue) {
1363 bestValue = bvalue + bnode_potential[bnode];
1370 if (bestNode == INVALID) {
1374 matching_cost += bestValue;
1377 for (BNodeIt it(*graph); it != INVALID; ++it) {
1378 if (bpred[it] != INVALID) {
1379 bnode_potential[it] += bdist[it];
1381 bnode_potential[it] += bdistMax;
1384 for (ANodeIt it(*graph); it != INVALID; ++it) {
1385 if (anode_matching[it] != INVALID) {
1386 Node bnode = graph->bNode(anode_matching[it]);
1387 if (bpred[bnode] != INVALID) {
1388 anode_potential[it] += bdist[bnode];
1390 anode_potential[it] += bdistMax;
1395 while (bestNode != INVALID) {
1396 UEdge uedge = bpred[bestNode];
1397 Node anode = graph->aNode(uedge);
1399 bnode_matching[bestNode] = uedge;
1400 if (anode_matching[anode] != INVALID) {
1401 bestNode = graph->bNode(anode_matching[anode]);
1405 anode_matching[anode] = uedge;
1412 /// \brief Starts the algorithm.
1414 /// Starts the algorithm. It runs augmenting phases until the
1415 /// optimal solution reached.
1417 while (augment()) {}
1420 /// \brief Runs the algorithm.
1422 /// It just initalize the algorithm and then start it.
1430 /// \name Query Functions
1431 /// The result of the %Matching algorithm can be obtained using these
1433 /// Before the use of these functions,
1434 /// either run() or start() must be called.
1438 /// \brief Gives back the potential in the NodeMap
1440 /// Gives back the potential in the NodeMap. The potential is optimal with
1441 /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
1442 /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
1444 template <typename PotentialMap>
1445 void potential(PotentialMap& potential) const {
1446 for (ANodeIt it(*graph); it != INVALID; ++it) {
1447 potential[it] = anode_potential[it];
1449 for (BNodeIt it(*graph); it != INVALID; ++it) {
1450 potential[it] = bnode_potential[it];
1454 /// \brief Set true all matching uedge in the map.
1456 /// Set true all matching uedge in the map. It does not change the
1457 /// value mapped to the other uedges.
1458 /// \return The number of the matching edges.
1459 template <typename MatchingMap>
1460 int quickMatching(MatchingMap& matching) const {
1461 for (ANodeIt it(*graph); it != INVALID; ++it) {
1462 if (anode_matching[it] != INVALID) {
1463 matching[anode_matching[it]] = true;
1466 return matching_size;
1469 /// \brief Set true all matching uedge in the map and the others to false.
1471 /// Set true all matching uedge in the map and the others to false.
1472 /// \return The number of the matching edges.
1473 template <typename MatchingMap>
1474 int matching(MatchingMap& matching) const {
1475 for (UEdgeIt it(*graph); it != INVALID; ++it) {
1476 matching[it] = it == anode_matching[graph->aNode(it)];
1478 return matching_size;
1482 /// \brief Return true if the given uedge is in the matching.
1484 /// It returns true if the given uedge is in the matching.
1485 bool matchingEdge(const UEdge& edge) const {
1486 return anode_matching[graph->aNode(edge)] == edge;
1489 /// \brief Returns the matching edge from the node.
1491 /// Returns the matching edge from the node. If there is not such
1492 /// edge it gives back \c INVALID.
1493 UEdge matchingEdge(const Node& node) const {
1494 if (graph->aNode(node)) {
1495 return anode_matching[node];
1497 return bnode_matching[node];
1501 /// \brief Gives back the sum of costs of the matching edges.
1503 /// Gives back the sum of costs of the matching edges.
1504 Value matchingCost() const {
1505 return matching_cost;
1508 /// \brief Gives back the number of the matching edges.
1510 /// Gives back the number of the matching edges.
1511 int matchingSize() const {
1512 return matching_size;
1519 void initStructures() {
1520 if (!_heap_cross_ref) {
1521 local_heap_cross_ref = true;
1522 _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1526 _heap = Traits::createHeap(*_heap_cross_ref);
1530 void destroyStructures() {
1531 if (local_heap_cross_ref) delete _heap_cross_ref;
1532 if (local_heap) delete _heap;
1538 const BpUGraph *graph;
1539 const CostMap* cost;
1541 ANodeMatchingMap anode_matching;
1542 BNodeMatchingMap bnode_matching;
1544 ANodePotentialMap anode_potential;
1545 BNodePotentialMap bnode_potential;
1547 Value matching_cost;
1550 HeapCrossRef *_heap_cross_ref;
1551 bool local_heap_cross_ref;
1558 /// \ingroup matching
1560 /// \brief Minimum cost maximum cardinality bipartite matching
1562 /// This function calculates the minimum cost matching of the maximum
1563 /// cardinality matchings of a bipartite graph. It gives back the matching
1564 /// in an undirected edge map.
1566 /// \param graph The bipartite graph.
1567 /// \param cost The undirected edge map which contains the costs.
1568 /// \retval matching The undirected edge map which will be set to
1570 /// \return The cost of the matching.
1571 template <typename BpUGraph, typename CostMap, typename MatchingMap>
1572 typename CostMap::Value
1573 minCostMaxBipartiteMatching(const BpUGraph& graph,
1574 const CostMap& cost,
1575 MatchingMap& matching) {
1576 MinCostMaxBipartiteMatching<BpUGraph, CostMap>
1577 bpmatching(graph, cost);
1579 bpmatching.matching(matching);
1580 return bpmatching.matchingCost();