1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2009
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_MATCHING_H
20 #define LEMON_MATCHING_H
27 #include <lemon/core.h>
28 #include <lemon/unionfind.h>
29 #include <lemon/bin_heap.h>
30 #include <lemon/maps.h>
34 ///\brief Maximum matching algorithms in general graphs.
40 /// \brief Maximum cardinality matching in general graphs
42 /// This class implements Edmonds' alternating forest matching algorithm
43 /// for finding a maximum cardinality matching in a general undirected graph.
44 /// It can be started from an arbitrary initial matching
45 /// (the default is the empty one).
47 /// The dual solution of the problem is a map of the nodes to
48 /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
49 /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
50 /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
51 /// with factor-critical components, the nodes in \c ODD/A form the
52 /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
53 /// a perfect matching. The number of the factor-critical components
54 /// minus the number of barrier nodes is a lower bound on the
55 /// unmatched nodes, and the matching is optimal if and only if this bound is
56 /// tight. This decomposition can be obtained using \ref status() or
57 /// \ref statusMap() after running the algorithm.
59 /// \tparam GR The undirected graph type the algorithm runs on.
60 template <typename GR>
64 /// The graph type of the algorithm
66 /// The type of the matching map
67 typedef typename Graph::template NodeMap<typename Graph::Arc>
70 ///\brief Status constants for Gallai-Edmonds decomposition.
72 ///These constants are used for indicating the Gallai-Edmonds
73 ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
74 ///induce a subgraph with factor-critical components, the nodes with
75 ///status \c ODD (or \c A) form the canonical barrier, and the nodes
76 ///with status \c MATCHED (or \c C) induce a subgraph having a
79 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.)
81 MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.)
83 ODD = -1, ///< = -1. (\c A is an alias for \c ODD.)
85 UNMATCHED = -2 ///< = -2.
88 /// The type of the status map
89 typedef typename Graph::template NodeMap<Status> StatusMap;
93 TEMPLATE_GRAPH_TYPEDEFS(Graph);
95 typedef UnionFindEnum<IntNodeMap> BlossomSet;
96 typedef ExtendFindEnum<IntNodeMap> TreeSet;
97 typedef RangeMap<Node> NodeIntMap;
98 typedef MatchingMap EarMap;
99 typedef std::vector<Node> NodeQueue;
102 MatchingMap* _matching;
107 IntNodeMap* _blossom_set_index;
108 BlossomSet* _blossom_set;
109 NodeIntMap* _blossom_rep;
111 IntNodeMap* _tree_set_index;
114 NodeQueue _node_queue;
115 int _process, _postpone, _last;
121 void createStructures() {
122 _node_num = countNodes(_graph);
124 _matching = new MatchingMap(_graph);
127 _status = new StatusMap(_graph);
130 _ear = new EarMap(_graph);
133 _blossom_set_index = new IntNodeMap(_graph);
134 _blossom_set = new BlossomSet(*_blossom_set_index);
137 _blossom_rep = new NodeIntMap(_node_num);
140 _tree_set_index = new IntNodeMap(_graph);
141 _tree_set = new TreeSet(*_tree_set_index);
143 _node_queue.resize(_node_num);
146 void destroyStructures() {
158 delete _blossom_set_index;
164 delete _tree_set_index;
169 void processDense(const Node& n) {
170 _process = _postpone = _last = 0;
171 _node_queue[_last++] = n;
173 while (_process != _last) {
174 Node u = _node_queue[_process++];
175 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
176 Node v = _graph.target(a);
177 if ((*_status)[v] == MATCHED) {
179 } else if ((*_status)[v] == UNMATCHED) {
186 while (_postpone != _last) {
187 Node u = _node_queue[_postpone++];
189 for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
190 Node v = _graph.target(a);
192 if ((*_status)[v] == EVEN) {
193 if (_blossom_set->find(u) != _blossom_set->find(v)) {
198 while (_process != _last) {
199 Node w = _node_queue[_process++];
200 for (OutArcIt b(_graph, w); b != INVALID; ++b) {
201 Node x = _graph.target(b);
202 if ((*_status)[x] == MATCHED) {
204 } else if ((*_status)[x] == UNMATCHED) {
214 void processSparse(const Node& n) {
215 _process = _last = 0;
216 _node_queue[_last++] = n;
217 while (_process != _last) {
218 Node u = _node_queue[_process++];
219 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
220 Node v = _graph.target(a);
222 if ((*_status)[v] == EVEN) {
223 if (_blossom_set->find(u) != _blossom_set->find(v)) {
226 } else if ((*_status)[v] == MATCHED) {
228 } else if ((*_status)[v] == UNMATCHED) {
236 void shrinkOnEdge(const Edge& e) {
240 std::set<Node> left_set, right_set;
242 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
243 left_set.insert(left);
245 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
246 right_set.insert(right);
249 if ((*_matching)[left] == INVALID) break;
250 left = _graph.target((*_matching)[left]);
251 left = (*_blossom_rep)[_blossom_set->
252 find(_graph.target((*_ear)[left]))];
253 if (right_set.find(left) != right_set.end()) {
257 left_set.insert(left);
259 if ((*_matching)[right] == INVALID) break;
260 right = _graph.target((*_matching)[right]);
261 right = (*_blossom_rep)[_blossom_set->
262 find(_graph.target((*_ear)[right]))];
263 if (left_set.find(right) != left_set.end()) {
267 right_set.insert(right);
270 if (nca == INVALID) {
271 if ((*_matching)[left] == INVALID) {
273 while (left_set.find(nca) == left_set.end()) {
274 nca = _graph.target((*_matching)[nca]);
275 nca =(*_blossom_rep)[_blossom_set->
276 find(_graph.target((*_ear)[nca]))];
280 while (right_set.find(nca) == right_set.end()) {
281 nca = _graph.target((*_matching)[nca]);
282 nca = (*_blossom_rep)[_blossom_set->
283 find(_graph.target((*_ear)[nca]))];
291 Node node = _graph.u(e);
292 Arc arc = _graph.direct(e, true);
293 Node base = (*_blossom_rep)[_blossom_set->find(node)];
295 while (base != nca) {
300 n = _graph.target((*_matching)[n]);
302 n = _graph.target(a);
303 (*_ear)[n] = _graph.oppositeArc(a);
305 node = _graph.target((*_matching)[base]);
306 _tree_set->erase(base);
307 _tree_set->erase(node);
308 _blossom_set->insert(node, _blossom_set->find(base));
309 (*_status)[node] = EVEN;
310 _node_queue[_last++] = node;
311 arc = _graph.oppositeArc((*_ear)[node]);
312 node = _graph.target((*_ear)[node]);
313 base = (*_blossom_rep)[_blossom_set->find(node)];
314 _blossom_set->join(_graph.target(arc), base);
318 (*_blossom_rep)[_blossom_set->find(nca)] = nca;
322 Node node = _graph.v(e);
323 Arc arc = _graph.direct(e, false);
324 Node base = (*_blossom_rep)[_blossom_set->find(node)];
326 while (base != nca) {
331 n = _graph.target((*_matching)[n]);
333 n = _graph.target(a);
334 (*_ear)[n] = _graph.oppositeArc(a);
336 node = _graph.target((*_matching)[base]);
337 _tree_set->erase(base);
338 _tree_set->erase(node);
339 _blossom_set->insert(node, _blossom_set->find(base));
340 (*_status)[node] = EVEN;
341 _node_queue[_last++] = node;
342 arc = _graph.oppositeArc((*_ear)[node]);
343 node = _graph.target((*_ear)[node]);
344 base = (*_blossom_rep)[_blossom_set->find(node)];
345 _blossom_set->join(_graph.target(arc), base);
349 (*_blossom_rep)[_blossom_set->find(nca)] = nca;
352 void extendOnArc(const Arc& a) {
353 Node base = _graph.source(a);
354 Node odd = _graph.target(a);
356 (*_ear)[odd] = _graph.oppositeArc(a);
357 Node even = _graph.target((*_matching)[odd]);
358 (*_blossom_rep)[_blossom_set->insert(even)] = even;
359 (*_status)[odd] = ODD;
360 (*_status)[even] = EVEN;
361 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
362 _tree_set->insert(odd, tree);
363 _tree_set->insert(even, tree);
364 _node_queue[_last++] = even;
368 void augmentOnArc(const Arc& a) {
369 Node even = _graph.source(a);
370 Node odd = _graph.target(a);
372 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
374 (*_matching)[odd] = _graph.oppositeArc(a);
375 (*_status)[odd] = MATCHED;
377 Arc arc = (*_matching)[even];
378 (*_matching)[even] = a;
380 while (arc != INVALID) {
381 odd = _graph.target(arc);
383 even = _graph.target(arc);
384 (*_matching)[odd] = arc;
385 arc = (*_matching)[even];
386 (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
389 for (typename TreeSet::ItemIt it(*_tree_set, tree);
390 it != INVALID; ++it) {
391 if ((*_status)[it] == ODD) {
392 (*_status)[it] = MATCHED;
394 int blossom = _blossom_set->find(it);
395 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
396 jt != INVALID; ++jt) {
397 (*_status)[jt] = MATCHED;
399 _blossom_set->eraseClass(blossom);
402 _tree_set->eraseClass(tree);
408 /// \brief Constructor
411 MaxMatching(const Graph& graph)
412 : _graph(graph), _matching(0), _status(0), _ear(0),
413 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
414 _tree_set_index(0), _tree_set(0) {}
420 /// \name Execution Control
421 /// The simplest way to execute the algorithm is to use the
422 /// \c run() member function.\n
423 /// If you need better control on the execution, you have to call
424 /// one of the functions \ref init(), \ref greedyInit() or
425 /// \ref matchingInit() first, then you can start the algorithm with
426 /// \ref startSparse() or \ref startDense().
430 /// \brief Set the initial matching to the empty matching.
432 /// This function sets the initial matching to the empty matching.
435 for(NodeIt n(_graph); n != INVALID; ++n) {
436 (*_matching)[n] = INVALID;
437 (*_status)[n] = UNMATCHED;
441 /// \brief Find an initial matching in a greedy way.
443 /// This function finds an initial matching in a greedy way.
446 for (NodeIt n(_graph); n != INVALID; ++n) {
447 (*_matching)[n] = INVALID;
448 (*_status)[n] = UNMATCHED;
450 for (NodeIt n(_graph); n != INVALID; ++n) {
451 if ((*_matching)[n] == INVALID) {
452 for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
453 Node v = _graph.target(a);
454 if ((*_matching)[v] == INVALID && v != n) {
456 (*_status)[n] = MATCHED;
457 (*_matching)[v] = _graph.oppositeArc(a);
458 (*_status)[v] = MATCHED;
467 /// \brief Initialize the matching from a map.
469 /// This function initializes the matching from a \c bool valued edge
470 /// map. This map should have the property that there are no two incident
471 /// edges with \c true value, i.e. it really contains a matching.
472 /// \return \c true if the map contains a matching.
473 template <typename MatchingMap>
474 bool matchingInit(const MatchingMap& matching) {
477 for (NodeIt n(_graph); n != INVALID; ++n) {
478 (*_matching)[n] = INVALID;
479 (*_status)[n] = UNMATCHED;
481 for(EdgeIt e(_graph); e!=INVALID; ++e) {
484 Node u = _graph.u(e);
485 if ((*_matching)[u] != INVALID) return false;
486 (*_matching)[u] = _graph.direct(e, true);
487 (*_status)[u] = MATCHED;
489 Node v = _graph.v(e);
490 if ((*_matching)[v] != INVALID) return false;
491 (*_matching)[v] = _graph.direct(e, false);
492 (*_status)[v] = MATCHED;
498 /// \brief Start Edmonds' algorithm
500 /// This function runs the original Edmonds' algorithm.
502 /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
503 /// called before using this function.
505 for(NodeIt n(_graph); n != INVALID; ++n) {
506 if ((*_status)[n] == UNMATCHED) {
507 (*_blossom_rep)[_blossom_set->insert(n)] = n;
508 _tree_set->insert(n);
509 (*_status)[n] = EVEN;
515 /// \brief Start Edmonds' algorithm with a heuristic improvement
518 /// This function runs Edmonds' algorithm with a heuristic of postponing
519 /// shrinks, therefore resulting in a faster algorithm for dense graphs.
521 /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
522 /// called before using this function.
524 for(NodeIt n(_graph); n != INVALID; ++n) {
525 if ((*_status)[n] == UNMATCHED) {
526 (*_blossom_rep)[_blossom_set->insert(n)] = n;
527 _tree_set->insert(n);
528 (*_status)[n] = EVEN;
535 /// \brief Run Edmonds' algorithm
537 /// This function runs Edmonds' algorithm. An additional heuristic of
538 /// postponing shrinks is used for relatively dense graphs
539 /// (for which <tt>m>=2*n</tt> holds).
541 if (countEdges(_graph) < 2 * countNodes(_graph)) {
552 /// \name Primal Solution
553 /// Functions to get the primal solution, i.e. the maximum matching.
557 /// \brief Return the size (cardinality) of the matching.
559 /// This function returns the size (cardinality) of the current matching.
560 /// After run() it returns the size of the maximum matching in the graph.
561 int matchingSize() const {
563 for (NodeIt n(_graph); n != INVALID; ++n) {
564 if ((*_matching)[n] != INVALID) {
571 /// \brief Return \c true if the given edge is in the matching.
573 /// This function returns \c true if the given edge is in the current
575 bool matching(const Edge& edge) const {
576 return edge == (*_matching)[_graph.u(edge)];
579 /// \brief Return the matching arc (or edge) incident to the given node.
581 /// This function returns the matching arc (or edge) incident to the
582 /// given node in the current matching or \c INVALID if the node is
583 /// not covered by the matching.
584 Arc matching(const Node& n) const {
585 return (*_matching)[n];
588 /// \brief Return a const reference to the matching map.
590 /// This function returns a const reference to a node map that stores
591 /// the matching arc (or edge) incident to each node.
592 const MatchingMap& matchingMap() const {
596 /// \brief Return the mate of the given node.
598 /// This function returns the mate of the given node in the current
599 /// matching or \c INVALID if the node is not covered by the matching.
600 Node mate(const Node& n) const {
601 return (*_matching)[n] != INVALID ?
602 _graph.target((*_matching)[n]) : INVALID;
607 /// \name Dual Solution
608 /// Functions to get the dual solution, i.e. the Gallai-Edmonds
613 /// \brief Return the status of the given node in the Edmonds-Gallai
616 /// This function returns the \ref Status "status" of the given node
617 /// in the Edmonds-Gallai decomposition.
618 Status status(const Node& n) const {
619 return (*_status)[n];
622 /// \brief Return a const reference to the status map, which stores
623 /// the Edmonds-Gallai decomposition.
625 /// This function returns a const reference to a node map that stores the
626 /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
627 const StatusMap& statusMap() const {
631 /// \brief Return \c true if the given node is in the barrier.
633 /// This function returns \c true if the given node is in the barrier.
634 bool barrier(const Node& n) const {
635 return (*_status)[n] == ODD;
642 /// \ingroup matching
644 /// \brief Weighted matching in general graphs
646 /// This class provides an efficient implementation of Edmond's
647 /// maximum weighted matching algorithm. The implementation is based
648 /// on extensive use of priority queues and provides
649 /// \f$O(nm\log n)\f$ time complexity.
651 /// The maximum weighted matching problem is to find a subset of the
652 /// edges in an undirected graph with maximum overall weight for which
653 /// each node has at most one incident edge.
654 /// It can be formulated with the following linear program.
655 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
656 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
657 \quad \forall B\in\mathcal{O}\f] */
658 /// \f[x_e \ge 0\quad \forall e\in E\f]
659 /// \f[\max \sum_{e\in E}x_ew_e\f]
660 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
661 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
662 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
663 /// subsets of the nodes.
665 /// The algorithm calculates an optimal matching and a proof of the
666 /// optimality. The solution of the dual problem can be used to check
667 /// the result of the algorithm. The dual linear problem is the
669 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
670 z_B \ge w_{uv} \quad \forall uv\in E\f] */
671 /// \f[y_u \ge 0 \quad \forall u \in V\f]
672 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
673 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
674 \frac{\vert B \vert - 1}{2}z_B\f] */
676 /// The algorithm can be executed with the run() function.
677 /// After it the matching (the primal solution) and the dual solution
678 /// can be obtained using the query functions and the
679 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
680 /// which is able to iterate on the nodes of a blossom.
681 /// If the value type is integer, then the dual solution is multiplied
682 /// by \ref MaxWeightedMatching::dualScale "4".
684 /// \tparam GR The undirected graph type the algorithm runs on.
685 /// \tparam WM The type edge weight map. The default type is
686 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
688 template <typename GR, typename WM>
690 template <typename GR,
691 typename WM = typename GR::template EdgeMap<int> >
693 class MaxWeightedMatching {
696 /// The graph type of the algorithm
698 /// The type of the edge weight map
699 typedef WM WeightMap;
700 /// The value type of the edge weights
701 typedef typename WeightMap::Value Value;
703 /// The type of the matching map
704 typedef typename Graph::template NodeMap<typename Graph::Arc>
707 /// \brief Scaling factor for dual solution
709 /// Scaling factor for dual solution. It is equal to 4 or 1
710 /// according to the value type.
711 static const int dualScale =
712 std::numeric_limits<Value>::is_integer ? 4 : 1;
716 TEMPLATE_GRAPH_TYPEDEFS(Graph);
718 typedef typename Graph::template NodeMap<Value> NodePotential;
719 typedef std::vector<Node> BlossomNodeList;
721 struct BlossomVariable {
725 BlossomVariable(int _begin, int _end, Value _value)
726 : begin(_begin), end(_end), value(_value) {}
730 typedef std::vector<BlossomVariable> BlossomPotential;
733 const WeightMap& _weight;
735 MatchingMap* _matching;
737 NodePotential* _node_potential;
739 BlossomPotential _blossom_potential;
740 BlossomNodeList _blossom_node_list;
745 typedef RangeMap<int> IntIntMap;
748 EVEN = -1, MATCHED = 0, ODD = 1
751 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
760 IntNodeMap *_blossom_index;
761 BlossomSet *_blossom_set;
762 RangeMap<BlossomData>* _blossom_data;
764 IntNodeMap *_node_index;
765 IntArcMap *_node_heap_index;
769 NodeData(IntArcMap& node_heap_index)
770 : heap(node_heap_index) {}
774 BinHeap<Value, IntArcMap> heap;
775 std::map<int, Arc> heap_index;
780 RangeMap<NodeData>* _node_data;
782 typedef ExtendFindEnum<IntIntMap> TreeSet;
784 IntIntMap *_tree_set_index;
787 IntNodeMap *_delta1_index;
788 BinHeap<Value, IntNodeMap> *_delta1;
790 IntIntMap *_delta2_index;
791 BinHeap<Value, IntIntMap> *_delta2;
793 IntEdgeMap *_delta3_index;
794 BinHeap<Value, IntEdgeMap> *_delta3;
796 IntIntMap *_delta4_index;
797 BinHeap<Value, IntIntMap> *_delta4;
801 void createStructures() {
802 _node_num = countNodes(_graph);
803 _blossom_num = _node_num * 3 / 2;
806 _matching = new MatchingMap(_graph);
808 if (!_node_potential) {
809 _node_potential = new NodePotential(_graph);
812 _blossom_index = new IntNodeMap(_graph);
813 _blossom_set = new BlossomSet(*_blossom_index);
814 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
818 _node_index = new IntNodeMap(_graph);
819 _node_heap_index = new IntArcMap(_graph);
820 _node_data = new RangeMap<NodeData>(_node_num,
821 NodeData(*_node_heap_index));
825 _tree_set_index = new IntIntMap(_blossom_num);
826 _tree_set = new TreeSet(*_tree_set_index);
829 _delta1_index = new IntNodeMap(_graph);
830 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
833 _delta2_index = new IntIntMap(_blossom_num);
834 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
837 _delta3_index = new IntEdgeMap(_graph);
838 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
841 _delta4_index = new IntIntMap(_blossom_num);
842 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
846 void destroyStructures() {
850 if (_node_potential) {
851 delete _node_potential;
854 delete _blossom_index;
856 delete _blossom_data;
861 delete _node_heap_index;
866 delete _tree_set_index;
870 delete _delta1_index;
874 delete _delta2_index;
878 delete _delta3_index;
882 delete _delta4_index;
887 void matchedToEven(int blossom, int tree) {
888 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
889 _delta2->erase(blossom);
892 if (!_blossom_set->trivial(blossom)) {
893 (*_blossom_data)[blossom].pot -=
894 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
897 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
900 _blossom_set->increase(n, std::numeric_limits<Value>::max());
901 int ni = (*_node_index)[n];
903 (*_node_data)[ni].heap.clear();
904 (*_node_data)[ni].heap_index.clear();
906 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
908 _delta1->push(n, (*_node_data)[ni].pot);
910 for (InArcIt e(_graph, n); e != INVALID; ++e) {
911 Node v = _graph.source(e);
912 int vb = _blossom_set->find(v);
913 int vi = (*_node_index)[v];
915 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
916 dualScale * _weight[e];
918 if ((*_blossom_data)[vb].status == EVEN) {
919 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
920 _delta3->push(e, rw / 2);
923 typename std::map<int, Arc>::iterator it =
924 (*_node_data)[vi].heap_index.find(tree);
926 if (it != (*_node_data)[vi].heap_index.end()) {
927 if ((*_node_data)[vi].heap[it->second] > rw) {
928 (*_node_data)[vi].heap.replace(it->second, e);
929 (*_node_data)[vi].heap.decrease(e, rw);
933 (*_node_data)[vi].heap.push(e, rw);
934 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
937 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
938 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
940 if ((*_blossom_data)[vb].status == MATCHED) {
941 if (_delta2->state(vb) != _delta2->IN_HEAP) {
942 _delta2->push(vb, _blossom_set->classPrio(vb) -
943 (*_blossom_data)[vb].offset);
944 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
945 (*_blossom_data)[vb].offset) {
946 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
947 (*_blossom_data)[vb].offset);
954 (*_blossom_data)[blossom].offset = 0;
957 void matchedToOdd(int blossom) {
958 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
959 _delta2->erase(blossom);
961 (*_blossom_data)[blossom].offset += _delta_sum;
962 if (!_blossom_set->trivial(blossom)) {
963 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
964 (*_blossom_data)[blossom].offset);
968 void evenToMatched(int blossom, int tree) {
969 if (!_blossom_set->trivial(blossom)) {
970 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
973 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
975 int ni = (*_node_index)[n];
976 (*_node_data)[ni].pot -= _delta_sum;
980 for (InArcIt e(_graph, n); e != INVALID; ++e) {
981 Node v = _graph.source(e);
982 int vb = _blossom_set->find(v);
983 int vi = (*_node_index)[v];
985 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
986 dualScale * _weight[e];
989 if (_delta3->state(e) == _delta3->IN_HEAP) {
992 } else if ((*_blossom_data)[vb].status == EVEN) {
994 if (_delta3->state(e) == _delta3->IN_HEAP) {
998 int vt = _tree_set->find(vb);
1002 Arc r = _graph.oppositeArc(e);
1004 typename std::map<int, Arc>::iterator it =
1005 (*_node_data)[ni].heap_index.find(vt);
1007 if (it != (*_node_data)[ni].heap_index.end()) {
1008 if ((*_node_data)[ni].heap[it->second] > rw) {
1009 (*_node_data)[ni].heap.replace(it->second, r);
1010 (*_node_data)[ni].heap.decrease(r, rw);
1014 (*_node_data)[ni].heap.push(r, rw);
1015 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1018 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1019 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1021 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1022 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1023 (*_blossom_data)[blossom].offset);
1024 } else if ((*_delta2)[blossom] >
1025 _blossom_set->classPrio(blossom) -
1026 (*_blossom_data)[blossom].offset){
1027 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1028 (*_blossom_data)[blossom].offset);
1034 typename std::map<int, Arc>::iterator it =
1035 (*_node_data)[vi].heap_index.find(tree);
1037 if (it != (*_node_data)[vi].heap_index.end()) {
1038 (*_node_data)[vi].heap.erase(it->second);
1039 (*_node_data)[vi].heap_index.erase(it);
1040 if ((*_node_data)[vi].heap.empty()) {
1041 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1042 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1043 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1046 if ((*_blossom_data)[vb].status == MATCHED) {
1047 if (_blossom_set->classPrio(vb) ==
1048 std::numeric_limits<Value>::max()) {
1050 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1051 (*_blossom_data)[vb].offset) {
1052 _delta2->increase(vb, _blossom_set->classPrio(vb) -
1053 (*_blossom_data)[vb].offset);
1062 void oddToMatched(int blossom) {
1063 (*_blossom_data)[blossom].offset -= _delta_sum;
1065 if (_blossom_set->classPrio(blossom) !=
1066 std::numeric_limits<Value>::max()) {
1067 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1068 (*_blossom_data)[blossom].offset);
1071 if (!_blossom_set->trivial(blossom)) {
1072 _delta4->erase(blossom);
1076 void oddToEven(int blossom, int tree) {
1077 if (!_blossom_set->trivial(blossom)) {
1078 _delta4->erase(blossom);
1079 (*_blossom_data)[blossom].pot -=
1080 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1083 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1084 n != INVALID; ++n) {
1085 int ni = (*_node_index)[n];
1087 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1089 (*_node_data)[ni].heap.clear();
1090 (*_node_data)[ni].heap_index.clear();
1091 (*_node_data)[ni].pot +=
1092 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1094 _delta1->push(n, (*_node_data)[ni].pot);
1096 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1097 Node v = _graph.source(e);
1098 int vb = _blossom_set->find(v);
1099 int vi = (*_node_index)[v];
1101 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1102 dualScale * _weight[e];
1104 if ((*_blossom_data)[vb].status == EVEN) {
1105 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1106 _delta3->push(e, rw / 2);
1110 typename std::map<int, Arc>::iterator it =
1111 (*_node_data)[vi].heap_index.find(tree);
1113 if (it != (*_node_data)[vi].heap_index.end()) {
1114 if ((*_node_data)[vi].heap[it->second] > rw) {
1115 (*_node_data)[vi].heap.replace(it->second, e);
1116 (*_node_data)[vi].heap.decrease(e, rw);
1120 (*_node_data)[vi].heap.push(e, rw);
1121 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1124 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1125 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1127 if ((*_blossom_data)[vb].status == MATCHED) {
1128 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1129 _delta2->push(vb, _blossom_set->classPrio(vb) -
1130 (*_blossom_data)[vb].offset);
1131 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1132 (*_blossom_data)[vb].offset) {
1133 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1134 (*_blossom_data)[vb].offset);
1141 (*_blossom_data)[blossom].offset = 0;
1144 void alternatePath(int even, int tree) {
1147 evenToMatched(even, tree);
1148 (*_blossom_data)[even].status = MATCHED;
1150 while ((*_blossom_data)[even].pred != INVALID) {
1151 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1152 (*_blossom_data)[odd].status = MATCHED;
1154 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1156 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1157 (*_blossom_data)[even].status = MATCHED;
1158 evenToMatched(even, tree);
1159 (*_blossom_data)[even].next =
1160 _graph.oppositeArc((*_blossom_data)[odd].pred);
1165 void destroyTree(int tree) {
1166 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1167 if ((*_blossom_data)[b].status == EVEN) {
1168 (*_blossom_data)[b].status = MATCHED;
1169 evenToMatched(b, tree);
1170 } else if ((*_blossom_data)[b].status == ODD) {
1171 (*_blossom_data)[b].status = MATCHED;
1175 _tree_set->eraseClass(tree);
1179 void unmatchNode(const Node& node) {
1180 int blossom = _blossom_set->find(node);
1181 int tree = _tree_set->find(blossom);
1183 alternatePath(blossom, tree);
1186 (*_blossom_data)[blossom].base = node;
1187 (*_blossom_data)[blossom].next = INVALID;
1190 void augmentOnEdge(const Edge& edge) {
1192 int left = _blossom_set->find(_graph.u(edge));
1193 int right = _blossom_set->find(_graph.v(edge));
1195 int left_tree = _tree_set->find(left);
1196 alternatePath(left, left_tree);
1197 destroyTree(left_tree);
1199 int right_tree = _tree_set->find(right);
1200 alternatePath(right, right_tree);
1201 destroyTree(right_tree);
1203 (*_blossom_data)[left].next = _graph.direct(edge, true);
1204 (*_blossom_data)[right].next = _graph.direct(edge, false);
1207 void augmentOnArc(const Arc& arc) {
1209 int left = _blossom_set->find(_graph.source(arc));
1210 int right = _blossom_set->find(_graph.target(arc));
1212 (*_blossom_data)[left].status = MATCHED;
1214 int right_tree = _tree_set->find(right);
1215 alternatePath(right, right_tree);
1216 destroyTree(right_tree);
1218 (*_blossom_data)[left].next = arc;
1219 (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1222 void extendOnArc(const Arc& arc) {
1223 int base = _blossom_set->find(_graph.target(arc));
1224 int tree = _tree_set->find(base);
1226 int odd = _blossom_set->find(_graph.source(arc));
1227 _tree_set->insert(odd, tree);
1228 (*_blossom_data)[odd].status = ODD;
1230 (*_blossom_data)[odd].pred = arc;
1232 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1233 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1234 _tree_set->insert(even, tree);
1235 (*_blossom_data)[even].status = EVEN;
1236 matchedToEven(even, tree);
1239 void shrinkOnEdge(const Edge& edge, int tree) {
1241 std::vector<int> left_path, right_path;
1244 std::set<int> left_set, right_set;
1245 int left = _blossom_set->find(_graph.u(edge));
1246 left_path.push_back(left);
1247 left_set.insert(left);
1249 int right = _blossom_set->find(_graph.v(edge));
1250 right_path.push_back(right);
1251 right_set.insert(right);
1255 if ((*_blossom_data)[left].pred == INVALID) break;
1258 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1259 left_path.push_back(left);
1261 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1262 left_path.push_back(left);
1264 left_set.insert(left);
1266 if (right_set.find(left) != right_set.end()) {
1271 if ((*_blossom_data)[right].pred == INVALID) break;
1274 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1275 right_path.push_back(right);
1277 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1278 right_path.push_back(right);
1280 right_set.insert(right);
1282 if (left_set.find(right) != left_set.end()) {
1290 if ((*_blossom_data)[left].pred == INVALID) {
1292 while (left_set.find(nca) == left_set.end()) {
1294 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1295 right_path.push_back(nca);
1297 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1298 right_path.push_back(nca);
1302 while (right_set.find(nca) == right_set.end()) {
1304 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1305 left_path.push_back(nca);
1307 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1308 left_path.push_back(nca);
1314 std::vector<int> subblossoms;
1317 prev = _graph.direct(edge, true);
1318 for (int i = 0; left_path[i] != nca; i += 2) {
1319 subblossoms.push_back(left_path[i]);
1320 (*_blossom_data)[left_path[i]].next = prev;
1321 _tree_set->erase(left_path[i]);
1323 subblossoms.push_back(left_path[i + 1]);
1324 (*_blossom_data)[left_path[i + 1]].status = EVEN;
1325 oddToEven(left_path[i + 1], tree);
1326 _tree_set->erase(left_path[i + 1]);
1327 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1331 while (right_path[k] != nca) ++k;
1333 subblossoms.push_back(nca);
1334 (*_blossom_data)[nca].next = prev;
1336 for (int i = k - 2; i >= 0; i -= 2) {
1337 subblossoms.push_back(right_path[i + 1]);
1338 (*_blossom_data)[right_path[i + 1]].status = EVEN;
1339 oddToEven(right_path[i + 1], tree);
1340 _tree_set->erase(right_path[i + 1]);
1342 (*_blossom_data)[right_path[i + 1]].next =
1343 (*_blossom_data)[right_path[i + 1]].pred;
1345 subblossoms.push_back(right_path[i]);
1346 _tree_set->erase(right_path[i]);
1350 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1352 for (int i = 0; i < int(subblossoms.size()); ++i) {
1353 if (!_blossom_set->trivial(subblossoms[i])) {
1354 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1356 (*_blossom_data)[subblossoms[i]].status = MATCHED;
1359 (*_blossom_data)[surface].pot = -2 * _delta_sum;
1360 (*_blossom_data)[surface].offset = 0;
1361 (*_blossom_data)[surface].status = EVEN;
1362 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1363 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1365 _tree_set->insert(surface, tree);
1366 _tree_set->erase(nca);
1369 void splitBlossom(int blossom) {
1370 Arc next = (*_blossom_data)[blossom].next;
1371 Arc pred = (*_blossom_data)[blossom].pred;
1373 int tree = _tree_set->find(blossom);
1375 (*_blossom_data)[blossom].status = MATCHED;
1376 oddToMatched(blossom);
1377 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1378 _delta2->erase(blossom);
1381 std::vector<int> subblossoms;
1382 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1384 Value offset = (*_blossom_data)[blossom].offset;
1385 int b = _blossom_set->find(_graph.source(pred));
1386 int d = _blossom_set->find(_graph.source(next));
1388 int ib = -1, id = -1;
1389 for (int i = 0; i < int(subblossoms.size()); ++i) {
1390 if (subblossoms[i] == b) ib = i;
1391 if (subblossoms[i] == d) id = i;
1393 (*_blossom_data)[subblossoms[i]].offset = offset;
1394 if (!_blossom_set->trivial(subblossoms[i])) {
1395 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1397 if (_blossom_set->classPrio(subblossoms[i]) !=
1398 std::numeric_limits<Value>::max()) {
1399 _delta2->push(subblossoms[i],
1400 _blossom_set->classPrio(subblossoms[i]) -
1401 (*_blossom_data)[subblossoms[i]].offset);
1405 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1406 for (int i = (id + 1) % subblossoms.size();
1407 i != ib; i = (i + 2) % subblossoms.size()) {
1408 int sb = subblossoms[i];
1409 int tb = subblossoms[(i + 1) % subblossoms.size()];
1410 (*_blossom_data)[sb].next =
1411 _graph.oppositeArc((*_blossom_data)[tb].next);
1414 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1415 int sb = subblossoms[i];
1416 int tb = subblossoms[(i + 1) % subblossoms.size()];
1417 int ub = subblossoms[(i + 2) % subblossoms.size()];
1419 (*_blossom_data)[sb].status = ODD;
1421 _tree_set->insert(sb, tree);
1422 (*_blossom_data)[sb].pred = pred;
1423 (*_blossom_data)[sb].next =
1424 _graph.oppositeArc((*_blossom_data)[tb].next);
1426 pred = (*_blossom_data)[ub].next;
1428 (*_blossom_data)[tb].status = EVEN;
1429 matchedToEven(tb, tree);
1430 _tree_set->insert(tb, tree);
1431 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1434 (*_blossom_data)[subblossoms[id]].status = ODD;
1435 matchedToOdd(subblossoms[id]);
1436 _tree_set->insert(subblossoms[id], tree);
1437 (*_blossom_data)[subblossoms[id]].next = next;
1438 (*_blossom_data)[subblossoms[id]].pred = pred;
1442 for (int i = (ib + 1) % subblossoms.size();
1443 i != id; i = (i + 2) % subblossoms.size()) {
1444 int sb = subblossoms[i];
1445 int tb = subblossoms[(i + 1) % subblossoms.size()];
1446 (*_blossom_data)[sb].next =
1447 _graph.oppositeArc((*_blossom_data)[tb].next);
1450 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1451 int sb = subblossoms[i];
1452 int tb = subblossoms[(i + 1) % subblossoms.size()];
1453 int ub = subblossoms[(i + 2) % subblossoms.size()];
1455 (*_blossom_data)[sb].status = ODD;
1457 _tree_set->insert(sb, tree);
1458 (*_blossom_data)[sb].next = next;
1459 (*_blossom_data)[sb].pred =
1460 _graph.oppositeArc((*_blossom_data)[tb].next);
1462 (*_blossom_data)[tb].status = EVEN;
1463 matchedToEven(tb, tree);
1464 _tree_set->insert(tb, tree);
1465 (*_blossom_data)[tb].pred =
1466 (*_blossom_data)[tb].next =
1467 _graph.oppositeArc((*_blossom_data)[ub].next);
1468 next = (*_blossom_data)[ub].next;
1471 (*_blossom_data)[subblossoms[ib]].status = ODD;
1472 matchedToOdd(subblossoms[ib]);
1473 _tree_set->insert(subblossoms[ib], tree);
1474 (*_blossom_data)[subblossoms[ib]].next = next;
1475 (*_blossom_data)[subblossoms[ib]].pred = pred;
1477 _tree_set->erase(blossom);
1480 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1481 if (_blossom_set->trivial(blossom)) {
1482 int bi = (*_node_index)[base];
1483 Value pot = (*_node_data)[bi].pot;
1485 (*_matching)[base] = matching;
1486 _blossom_node_list.push_back(base);
1487 (*_node_potential)[base] = pot;
1490 Value pot = (*_blossom_data)[blossom].pot;
1491 int bn = _blossom_node_list.size();
1493 std::vector<int> subblossoms;
1494 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1495 int b = _blossom_set->find(base);
1497 for (int i = 0; i < int(subblossoms.size()); ++i) {
1498 if (subblossoms[i] == b) { ib = i; break; }
1501 for (int i = 1; i < int(subblossoms.size()); i += 2) {
1502 int sb = subblossoms[(ib + i) % subblossoms.size()];
1503 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1505 Arc m = (*_blossom_data)[tb].next;
1506 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1507 extractBlossom(tb, _graph.source(m), m);
1509 extractBlossom(subblossoms[ib], base, matching);
1511 int en = _blossom_node_list.size();
1513 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1517 void extractMatching() {
1518 std::vector<int> blossoms;
1519 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1520 blossoms.push_back(c);
1523 for (int i = 0; i < int(blossoms.size()); ++i) {
1524 if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1526 Value offset = (*_blossom_data)[blossoms[i]].offset;
1527 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1528 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1529 n != INVALID; ++n) {
1530 (*_node_data)[(*_node_index)[n]].pot -= offset;
1533 Arc matching = (*_blossom_data)[blossoms[i]].next;
1534 Node base = _graph.source(matching);
1535 extractBlossom(blossoms[i], base, matching);
1537 Node base = (*_blossom_data)[blossoms[i]].base;
1538 extractBlossom(blossoms[i], base, INVALID);
1545 /// \brief Constructor
1548 MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1549 : _graph(graph), _weight(weight), _matching(0),
1550 _node_potential(0), _blossom_potential(), _blossom_node_list(),
1551 _node_num(0), _blossom_num(0),
1553 _blossom_index(0), _blossom_set(0), _blossom_data(0),
1554 _node_index(0), _node_heap_index(0), _node_data(0),
1555 _tree_set_index(0), _tree_set(0),
1557 _delta1_index(0), _delta1(0),
1558 _delta2_index(0), _delta2(0),
1559 _delta3_index(0), _delta3(0),
1560 _delta4_index(0), _delta4(0),
1564 ~MaxWeightedMatching() {
1565 destroyStructures();
1568 /// \name Execution Control
1569 /// The simplest way to execute the algorithm is to use the
1570 /// \ref run() member function.
1574 /// \brief Initialize the algorithm
1576 /// This function initializes the algorithm.
1580 for (ArcIt e(_graph); e != INVALID; ++e) {
1581 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1583 for (NodeIt n(_graph); n != INVALID; ++n) {
1584 (*_delta1_index)[n] = _delta1->PRE_HEAP;
1586 for (EdgeIt e(_graph); e != INVALID; ++e) {
1587 (*_delta3_index)[e] = _delta3->PRE_HEAP;
1589 for (int i = 0; i < _blossom_num; ++i) {
1590 (*_delta2_index)[i] = _delta2->PRE_HEAP;
1591 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1595 for (NodeIt n(_graph); n != INVALID; ++n) {
1597 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1598 if (_graph.target(e) == n) continue;
1599 if ((dualScale * _weight[e]) / 2 > max) {
1600 max = (dualScale * _weight[e]) / 2;
1603 (*_node_index)[n] = index;
1604 (*_node_data)[index].pot = max;
1605 _delta1->push(n, max);
1607 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1609 _tree_set->insert(blossom);
1611 (*_blossom_data)[blossom].status = EVEN;
1612 (*_blossom_data)[blossom].pred = INVALID;
1613 (*_blossom_data)[blossom].next = INVALID;
1614 (*_blossom_data)[blossom].pot = 0;
1615 (*_blossom_data)[blossom].offset = 0;
1618 for (EdgeIt e(_graph); e != INVALID; ++e) {
1619 int si = (*_node_index)[_graph.u(e)];
1620 int ti = (*_node_index)[_graph.v(e)];
1621 if (_graph.u(e) != _graph.v(e)) {
1622 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1623 dualScale * _weight[e]) / 2);
1628 /// \brief Start the algorithm
1630 /// This function starts the algorithm.
1632 /// \pre \ref init() must be called before using this function.
1638 int unmatched = _node_num;
1639 while (unmatched > 0) {
1640 Value d1 = !_delta1->empty() ?
1641 _delta1->prio() : std::numeric_limits<Value>::max();
1643 Value d2 = !_delta2->empty() ?
1644 _delta2->prio() : std::numeric_limits<Value>::max();
1646 Value d3 = !_delta3->empty() ?
1647 _delta3->prio() : std::numeric_limits<Value>::max();
1649 Value d4 = !_delta4->empty() ?
1650 _delta4->prio() : std::numeric_limits<Value>::max();
1652 _delta_sum = d3; OpType ot = D3;
1653 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1654 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1655 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1660 Node n = _delta1->top();
1667 int blossom = _delta2->top();
1668 Node n = _blossom_set->classTop(blossom);
1669 Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1670 if ((*_blossom_data)[blossom].next == INVALID) {
1680 Edge e = _delta3->top();
1682 int left_blossom = _blossom_set->find(_graph.u(e));
1683 int right_blossom = _blossom_set->find(_graph.v(e));
1685 if (left_blossom == right_blossom) {
1688 int left_tree = _tree_set->find(left_blossom);
1689 int right_tree = _tree_set->find(right_blossom);
1691 if (left_tree == right_tree) {
1692 shrinkOnEdge(e, left_tree);
1700 splitBlossom(_delta4->top());
1707 /// \brief Run the algorithm.
1709 /// This method runs the \c %MaxWeightedMatching algorithm.
1711 /// \note mwm.run() is just a shortcut of the following code.
1723 /// \name Primal Solution
1724 /// Functions to get the primal solution, i.e. the maximum weighted
1726 /// Either \ref run() or \ref start() function should be called before
1731 /// \brief Return the weight of the matching.
1733 /// This function returns the weight of the found matching.
1735 /// \pre Either run() or start() must be called before using this function.
1736 Value matchingWeight() const {
1738 for (NodeIt n(_graph); n != INVALID; ++n) {
1739 if ((*_matching)[n] != INVALID) {
1740 sum += _weight[(*_matching)[n]];
1746 /// \brief Return the size (cardinality) of the matching.
1748 /// This function returns the size (cardinality) of the found matching.
1750 /// \pre Either run() or start() must be called before using this function.
1751 int matchingSize() const {
1753 for (NodeIt n(_graph); n != INVALID; ++n) {
1754 if ((*_matching)[n] != INVALID) {
1761 /// \brief Return \c true if the given edge is in the matching.
1763 /// This function returns \c true if the given edge is in the found
1766 /// \pre Either run() or start() must be called before using this function.
1767 bool matching(const Edge& edge) const {
1768 return edge == (*_matching)[_graph.u(edge)];
1771 /// \brief Return the matching arc (or edge) incident to the given node.
1773 /// This function returns the matching arc (or edge) incident to the
1774 /// given node in the found matching or \c INVALID if the node is
1775 /// not covered by the matching.
1777 /// \pre Either run() or start() must be called before using this function.
1778 Arc matching(const Node& node) const {
1779 return (*_matching)[node];
1782 /// \brief Return a const reference to the matching map.
1784 /// This function returns a const reference to a node map that stores
1785 /// the matching arc (or edge) incident to each node.
1786 const MatchingMap& matchingMap() const {
1790 /// \brief Return the mate of the given node.
1792 /// This function returns the mate of the given node in the found
1793 /// matching or \c INVALID if the node is not covered by the matching.
1795 /// \pre Either run() or start() must be called before using this function.
1796 Node mate(const Node& node) const {
1797 return (*_matching)[node] != INVALID ?
1798 _graph.target((*_matching)[node]) : INVALID;
1803 /// \name Dual Solution
1804 /// Functions to get the dual solution.\n
1805 /// Either \ref run() or \ref start() function should be called before
1810 /// \brief Return the value of the dual solution.
1812 /// This function returns the value of the dual solution.
1813 /// It should be equal to the primal value scaled by \ref dualScale
1816 /// \pre Either run() or start() must be called before using this function.
1817 Value dualValue() const {
1819 for (NodeIt n(_graph); n != INVALID; ++n) {
1820 sum += nodeValue(n);
1822 for (int i = 0; i < blossomNum(); ++i) {
1823 sum += blossomValue(i) * (blossomSize(i) / 2);
1828 /// \brief Return the dual value (potential) of the given node.
1830 /// This function returns the dual value (potential) of the given node.
1832 /// \pre Either run() or start() must be called before using this function.
1833 Value nodeValue(const Node& n) const {
1834 return (*_node_potential)[n];
1837 /// \brief Return the number of the blossoms in the basis.
1839 /// This function returns the number of the blossoms in the basis.
1841 /// \pre Either run() or start() must be called before using this function.
1843 int blossomNum() const {
1844 return _blossom_potential.size();
1847 /// \brief Return the number of the nodes in the given blossom.
1849 /// This function returns the number of the nodes in the given blossom.
1851 /// \pre Either run() or start() must be called before using this function.
1853 int blossomSize(int k) const {
1854 return _blossom_potential[k].end - _blossom_potential[k].begin;
1857 /// \brief Return the dual value (ptential) of the given blossom.
1859 /// This function returns the dual value (ptential) of the given blossom.
1861 /// \pre Either run() or start() must be called before using this function.
1862 Value blossomValue(int k) const {
1863 return _blossom_potential[k].value;
1866 /// \brief Iterator for obtaining the nodes of a blossom.
1868 /// This class provides an iterator for obtaining the nodes of the
1869 /// given blossom. It lists a subset of the nodes.
1870 /// Before using this iterator, you must allocate a
1871 /// MaxWeightedMatching class and execute it.
1875 /// \brief Constructor.
1877 /// Constructor to get the nodes of the given variable.
1879 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1880 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1881 /// called before initializing this iterator.
1882 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1883 : _algorithm(&algorithm)
1885 _index = _algorithm->_blossom_potential[variable].begin;
1886 _last = _algorithm->_blossom_potential[variable].end;
1889 /// \brief Conversion to \c Node.
1891 /// Conversion to \c Node.
1892 operator Node() const {
1893 return _algorithm->_blossom_node_list[_index];
1896 /// \brief Increment operator.
1898 /// Increment operator.
1899 BlossomIt& operator++() {
1904 /// \brief Validity checking
1906 /// Checks whether the iterator is invalid.
1907 bool operator==(Invalid) const { return _index == _last; }
1909 /// \brief Validity checking
1911 /// Checks whether the iterator is valid.
1912 bool operator!=(Invalid) const { return _index != _last; }
1915 const MaxWeightedMatching* _algorithm;
1924 /// \ingroup matching
1926 /// \brief Weighted perfect matching in general graphs
1928 /// This class provides an efficient implementation of Edmond's
1929 /// maximum weighted perfect matching algorithm. The implementation
1930 /// is based on extensive use of priority queues and provides
1931 /// \f$O(nm\log n)\f$ time complexity.
1933 /// The maximum weighted perfect matching problem is to find a subset of
1934 /// the edges in an undirected graph with maximum overall weight for which
1935 /// each node has exactly one incident edge.
1936 /// It can be formulated with the following linear program.
1937 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1938 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1939 \quad \forall B\in\mathcal{O}\f] */
1940 /// \f[x_e \ge 0\quad \forall e\in E\f]
1941 /// \f[\max \sum_{e\in E}x_ew_e\f]
1942 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1943 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1944 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1945 /// subsets of the nodes.
1947 /// The algorithm calculates an optimal matching and a proof of the
1948 /// optimality. The solution of the dual problem can be used to check
1949 /// the result of the algorithm. The dual linear problem is the
1951 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1952 w_{uv} \quad \forall uv\in E\f] */
1953 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1954 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1955 \frac{\vert B \vert - 1}{2}z_B\f] */
1957 /// The algorithm can be executed with the run() function.
1958 /// After it the matching (the primal solution) and the dual solution
1959 /// can be obtained using the query functions and the
1960 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
1961 /// which is able to iterate on the nodes of a blossom.
1962 /// If the value type is integer, then the dual solution is multiplied
1963 /// by \ref MaxWeightedMatching::dualScale "4".
1965 /// \tparam GR The undirected graph type the algorithm runs on.
1966 /// \tparam WM The type edge weight map. The default type is
1967 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1969 template <typename GR, typename WM>
1971 template <typename GR,
1972 typename WM = typename GR::template EdgeMap<int> >
1974 class MaxWeightedPerfectMatching {
1977 /// The graph type of the algorithm
1979 /// The type of the edge weight map
1980 typedef WM WeightMap;
1981 /// The value type of the edge weights
1982 typedef typename WeightMap::Value Value;
1984 /// \brief Scaling factor for dual solution
1986 /// Scaling factor for dual solution, it is equal to 4 or 1
1987 /// according to the value type.
1988 static const int dualScale =
1989 std::numeric_limits<Value>::is_integer ? 4 : 1;
1991 /// The type of the matching map
1992 typedef typename Graph::template NodeMap<typename Graph::Arc>
1997 TEMPLATE_GRAPH_TYPEDEFS(Graph);
1999 typedef typename Graph::template NodeMap<Value> NodePotential;
2000 typedef std::vector<Node> BlossomNodeList;
2002 struct BlossomVariable {
2006 BlossomVariable(int _begin, int _end, Value _value)
2007 : begin(_begin), end(_end), value(_value) {}
2011 typedef std::vector<BlossomVariable> BlossomPotential;
2013 const Graph& _graph;
2014 const WeightMap& _weight;
2016 MatchingMap* _matching;
2018 NodePotential* _node_potential;
2020 BlossomPotential _blossom_potential;
2021 BlossomNodeList _blossom_node_list;
2026 typedef RangeMap<int> IntIntMap;
2029 EVEN = -1, MATCHED = 0, ODD = 1
2032 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2033 struct BlossomData {
2040 IntNodeMap *_blossom_index;
2041 BlossomSet *_blossom_set;
2042 RangeMap<BlossomData>* _blossom_data;
2044 IntNodeMap *_node_index;
2045 IntArcMap *_node_heap_index;
2049 NodeData(IntArcMap& node_heap_index)
2050 : heap(node_heap_index) {}
2054 BinHeap<Value, IntArcMap> heap;
2055 std::map<int, Arc> heap_index;
2060 RangeMap<NodeData>* _node_data;
2062 typedef ExtendFindEnum<IntIntMap> TreeSet;
2064 IntIntMap *_tree_set_index;
2067 IntIntMap *_delta2_index;
2068 BinHeap<Value, IntIntMap> *_delta2;
2070 IntEdgeMap *_delta3_index;
2071 BinHeap<Value, IntEdgeMap> *_delta3;
2073 IntIntMap *_delta4_index;
2074 BinHeap<Value, IntIntMap> *_delta4;
2078 void createStructures() {
2079 _node_num = countNodes(_graph);
2080 _blossom_num = _node_num * 3 / 2;
2083 _matching = new MatchingMap(_graph);
2085 if (!_node_potential) {
2086 _node_potential = new NodePotential(_graph);
2088 if (!_blossom_set) {
2089 _blossom_index = new IntNodeMap(_graph);
2090 _blossom_set = new BlossomSet(*_blossom_index);
2091 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2095 _node_index = new IntNodeMap(_graph);
2096 _node_heap_index = new IntArcMap(_graph);
2097 _node_data = new RangeMap<NodeData>(_node_num,
2098 NodeData(*_node_heap_index));
2102 _tree_set_index = new IntIntMap(_blossom_num);
2103 _tree_set = new TreeSet(*_tree_set_index);
2106 _delta2_index = new IntIntMap(_blossom_num);
2107 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2110 _delta3_index = new IntEdgeMap(_graph);
2111 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2114 _delta4_index = new IntIntMap(_blossom_num);
2115 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2119 void destroyStructures() {
2123 if (_node_potential) {
2124 delete _node_potential;
2127 delete _blossom_index;
2128 delete _blossom_set;
2129 delete _blossom_data;
2134 delete _node_heap_index;
2139 delete _tree_set_index;
2143 delete _delta2_index;
2147 delete _delta3_index;
2151 delete _delta4_index;
2156 void matchedToEven(int blossom, int tree) {
2157 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2158 _delta2->erase(blossom);
2161 if (!_blossom_set->trivial(blossom)) {
2162 (*_blossom_data)[blossom].pot -=
2163 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2166 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2167 n != INVALID; ++n) {
2169 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2170 int ni = (*_node_index)[n];
2172 (*_node_data)[ni].heap.clear();
2173 (*_node_data)[ni].heap_index.clear();
2175 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2177 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2178 Node v = _graph.source(e);
2179 int vb = _blossom_set->find(v);
2180 int vi = (*_node_index)[v];
2182 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2183 dualScale * _weight[e];
2185 if ((*_blossom_data)[vb].status == EVEN) {
2186 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2187 _delta3->push(e, rw / 2);
2190 typename std::map<int, Arc>::iterator it =
2191 (*_node_data)[vi].heap_index.find(tree);
2193 if (it != (*_node_data)[vi].heap_index.end()) {
2194 if ((*_node_data)[vi].heap[it->second] > rw) {
2195 (*_node_data)[vi].heap.replace(it->second, e);
2196 (*_node_data)[vi].heap.decrease(e, rw);
2200 (*_node_data)[vi].heap.push(e, rw);
2201 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2204 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2205 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2207 if ((*_blossom_data)[vb].status == MATCHED) {
2208 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2209 _delta2->push(vb, _blossom_set->classPrio(vb) -
2210 (*_blossom_data)[vb].offset);
2211 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2212 (*_blossom_data)[vb].offset){
2213 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2214 (*_blossom_data)[vb].offset);
2221 (*_blossom_data)[blossom].offset = 0;
2224 void matchedToOdd(int blossom) {
2225 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2226 _delta2->erase(blossom);
2228 (*_blossom_data)[blossom].offset += _delta_sum;
2229 if (!_blossom_set->trivial(blossom)) {
2230 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2231 (*_blossom_data)[blossom].offset);
2235 void evenToMatched(int blossom, int tree) {
2236 if (!_blossom_set->trivial(blossom)) {
2237 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2240 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2241 n != INVALID; ++n) {
2242 int ni = (*_node_index)[n];
2243 (*_node_data)[ni].pot -= _delta_sum;
2245 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2246 Node v = _graph.source(e);
2247 int vb = _blossom_set->find(v);
2248 int vi = (*_node_index)[v];
2250 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2251 dualScale * _weight[e];
2253 if (vb == blossom) {
2254 if (_delta3->state(e) == _delta3->IN_HEAP) {
2257 } else if ((*_blossom_data)[vb].status == EVEN) {
2259 if (_delta3->state(e) == _delta3->IN_HEAP) {
2263 int vt = _tree_set->find(vb);
2267 Arc r = _graph.oppositeArc(e);
2269 typename std::map<int, Arc>::iterator it =
2270 (*_node_data)[ni].heap_index.find(vt);
2272 if (it != (*_node_data)[ni].heap_index.end()) {
2273 if ((*_node_data)[ni].heap[it->second] > rw) {
2274 (*_node_data)[ni].heap.replace(it->second, r);
2275 (*_node_data)[ni].heap.decrease(r, rw);
2279 (*_node_data)[ni].heap.push(r, rw);
2280 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2283 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2284 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2286 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2287 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2288 (*_blossom_data)[blossom].offset);
2289 } else if ((*_delta2)[blossom] >
2290 _blossom_set->classPrio(blossom) -
2291 (*_blossom_data)[blossom].offset){
2292 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2293 (*_blossom_data)[blossom].offset);
2299 typename std::map<int, Arc>::iterator it =
2300 (*_node_data)[vi].heap_index.find(tree);
2302 if (it != (*_node_data)[vi].heap_index.end()) {
2303 (*_node_data)[vi].heap.erase(it->second);
2304 (*_node_data)[vi].heap_index.erase(it);
2305 if ((*_node_data)[vi].heap.empty()) {
2306 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2307 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2308 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2311 if ((*_blossom_data)[vb].status == MATCHED) {
2312 if (_blossom_set->classPrio(vb) ==
2313 std::numeric_limits<Value>::max()) {
2315 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2316 (*_blossom_data)[vb].offset) {
2317 _delta2->increase(vb, _blossom_set->classPrio(vb) -
2318 (*_blossom_data)[vb].offset);
2327 void oddToMatched(int blossom) {
2328 (*_blossom_data)[blossom].offset -= _delta_sum;
2330 if (_blossom_set->classPrio(blossom) !=
2331 std::numeric_limits<Value>::max()) {
2332 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2333 (*_blossom_data)[blossom].offset);
2336 if (!_blossom_set->trivial(blossom)) {
2337 _delta4->erase(blossom);
2341 void oddToEven(int blossom, int tree) {
2342 if (!_blossom_set->trivial(blossom)) {
2343 _delta4->erase(blossom);
2344 (*_blossom_data)[blossom].pot -=
2345 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2348 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2349 n != INVALID; ++n) {
2350 int ni = (*_node_index)[n];
2352 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2354 (*_node_data)[ni].heap.clear();
2355 (*_node_data)[ni].heap_index.clear();
2356 (*_node_data)[ni].pot +=
2357 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2359 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2360 Node v = _graph.source(e);
2361 int vb = _blossom_set->find(v);
2362 int vi = (*_node_index)[v];
2364 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2365 dualScale * _weight[e];
2367 if ((*_blossom_data)[vb].status == EVEN) {
2368 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2369 _delta3->push(e, rw / 2);
2373 typename std::map<int, Arc>::iterator it =
2374 (*_node_data)[vi].heap_index.find(tree);
2376 if (it != (*_node_data)[vi].heap_index.end()) {
2377 if ((*_node_data)[vi].heap[it->second] > rw) {
2378 (*_node_data)[vi].heap.replace(it->second, e);
2379 (*_node_data)[vi].heap.decrease(e, rw);
2383 (*_node_data)[vi].heap.push(e, rw);
2384 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2387 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2388 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2390 if ((*_blossom_data)[vb].status == MATCHED) {
2391 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2392 _delta2->push(vb, _blossom_set->classPrio(vb) -
2393 (*_blossom_data)[vb].offset);
2394 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2395 (*_blossom_data)[vb].offset) {
2396 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2397 (*_blossom_data)[vb].offset);
2404 (*_blossom_data)[blossom].offset = 0;
2407 void alternatePath(int even, int tree) {
2410 evenToMatched(even, tree);
2411 (*_blossom_data)[even].status = MATCHED;
2413 while ((*_blossom_data)[even].pred != INVALID) {
2414 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2415 (*_blossom_data)[odd].status = MATCHED;
2417 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2419 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2420 (*_blossom_data)[even].status = MATCHED;
2421 evenToMatched(even, tree);
2422 (*_blossom_data)[even].next =
2423 _graph.oppositeArc((*_blossom_data)[odd].pred);
2428 void destroyTree(int tree) {
2429 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2430 if ((*_blossom_data)[b].status == EVEN) {
2431 (*_blossom_data)[b].status = MATCHED;
2432 evenToMatched(b, tree);
2433 } else if ((*_blossom_data)[b].status == ODD) {
2434 (*_blossom_data)[b].status = MATCHED;
2438 _tree_set->eraseClass(tree);
2441 void augmentOnEdge(const Edge& edge) {
2443 int left = _blossom_set->find(_graph.u(edge));
2444 int right = _blossom_set->find(_graph.v(edge));
2446 int left_tree = _tree_set->find(left);
2447 alternatePath(left, left_tree);
2448 destroyTree(left_tree);
2450 int right_tree = _tree_set->find(right);
2451 alternatePath(right, right_tree);
2452 destroyTree(right_tree);
2454 (*_blossom_data)[left].next = _graph.direct(edge, true);
2455 (*_blossom_data)[right].next = _graph.direct(edge, false);
2458 void extendOnArc(const Arc& arc) {
2459 int base = _blossom_set->find(_graph.target(arc));
2460 int tree = _tree_set->find(base);
2462 int odd = _blossom_set->find(_graph.source(arc));
2463 _tree_set->insert(odd, tree);
2464 (*_blossom_data)[odd].status = ODD;
2466 (*_blossom_data)[odd].pred = arc;
2468 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2469 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2470 _tree_set->insert(even, tree);
2471 (*_blossom_data)[even].status = EVEN;
2472 matchedToEven(even, tree);
2475 void shrinkOnEdge(const Edge& edge, int tree) {
2477 std::vector<int> left_path, right_path;
2480 std::set<int> left_set, right_set;
2481 int left = _blossom_set->find(_graph.u(edge));
2482 left_path.push_back(left);
2483 left_set.insert(left);
2485 int right = _blossom_set->find(_graph.v(edge));
2486 right_path.push_back(right);
2487 right_set.insert(right);
2491 if ((*_blossom_data)[left].pred == INVALID) break;
2494 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2495 left_path.push_back(left);
2497 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2498 left_path.push_back(left);
2500 left_set.insert(left);
2502 if (right_set.find(left) != right_set.end()) {
2507 if ((*_blossom_data)[right].pred == INVALID) break;
2510 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2511 right_path.push_back(right);
2513 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2514 right_path.push_back(right);
2516 right_set.insert(right);
2518 if (left_set.find(right) != left_set.end()) {
2526 if ((*_blossom_data)[left].pred == INVALID) {
2528 while (left_set.find(nca) == left_set.end()) {
2530 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2531 right_path.push_back(nca);
2533 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2534 right_path.push_back(nca);
2538 while (right_set.find(nca) == right_set.end()) {
2540 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2541 left_path.push_back(nca);
2543 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2544 left_path.push_back(nca);
2550 std::vector<int> subblossoms;
2553 prev = _graph.direct(edge, true);
2554 for (int i = 0; left_path[i] != nca; i += 2) {
2555 subblossoms.push_back(left_path[i]);
2556 (*_blossom_data)[left_path[i]].next = prev;
2557 _tree_set->erase(left_path[i]);
2559 subblossoms.push_back(left_path[i + 1]);
2560 (*_blossom_data)[left_path[i + 1]].status = EVEN;
2561 oddToEven(left_path[i + 1], tree);
2562 _tree_set->erase(left_path[i + 1]);
2563 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2567 while (right_path[k] != nca) ++k;
2569 subblossoms.push_back(nca);
2570 (*_blossom_data)[nca].next = prev;
2572 for (int i = k - 2; i >= 0; i -= 2) {
2573 subblossoms.push_back(right_path[i + 1]);
2574 (*_blossom_data)[right_path[i + 1]].status = EVEN;
2575 oddToEven(right_path[i + 1], tree);
2576 _tree_set->erase(right_path[i + 1]);
2578 (*_blossom_data)[right_path[i + 1]].next =
2579 (*_blossom_data)[right_path[i + 1]].pred;
2581 subblossoms.push_back(right_path[i]);
2582 _tree_set->erase(right_path[i]);
2586 _blossom_set->join(subblossoms.begin(), subblossoms.end());
2588 for (int i = 0; i < int(subblossoms.size()); ++i) {
2589 if (!_blossom_set->trivial(subblossoms[i])) {
2590 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2592 (*_blossom_data)[subblossoms[i]].status = MATCHED;
2595 (*_blossom_data)[surface].pot = -2 * _delta_sum;
2596 (*_blossom_data)[surface].offset = 0;
2597 (*_blossom_data)[surface].status = EVEN;
2598 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2599 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2601 _tree_set->insert(surface, tree);
2602 _tree_set->erase(nca);
2605 void splitBlossom(int blossom) {
2606 Arc next = (*_blossom_data)[blossom].next;
2607 Arc pred = (*_blossom_data)[blossom].pred;
2609 int tree = _tree_set->find(blossom);
2611 (*_blossom_data)[blossom].status = MATCHED;
2612 oddToMatched(blossom);
2613 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2614 _delta2->erase(blossom);
2617 std::vector<int> subblossoms;
2618 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2620 Value offset = (*_blossom_data)[blossom].offset;
2621 int b = _blossom_set->find(_graph.source(pred));
2622 int d = _blossom_set->find(_graph.source(next));
2624 int ib = -1, id = -1;
2625 for (int i = 0; i < int(subblossoms.size()); ++i) {
2626 if (subblossoms[i] == b) ib = i;
2627 if (subblossoms[i] == d) id = i;
2629 (*_blossom_data)[subblossoms[i]].offset = offset;
2630 if (!_blossom_set->trivial(subblossoms[i])) {
2631 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2633 if (_blossom_set->classPrio(subblossoms[i]) !=
2634 std::numeric_limits<Value>::max()) {
2635 _delta2->push(subblossoms[i],
2636 _blossom_set->classPrio(subblossoms[i]) -
2637 (*_blossom_data)[subblossoms[i]].offset);
2641 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2642 for (int i = (id + 1) % subblossoms.size();
2643 i != ib; i = (i + 2) % subblossoms.size()) {
2644 int sb = subblossoms[i];
2645 int tb = subblossoms[(i + 1) % subblossoms.size()];
2646 (*_blossom_data)[sb].next =
2647 _graph.oppositeArc((*_blossom_data)[tb].next);
2650 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2651 int sb = subblossoms[i];
2652 int tb = subblossoms[(i + 1) % subblossoms.size()];
2653 int ub = subblossoms[(i + 2) % subblossoms.size()];
2655 (*_blossom_data)[sb].status = ODD;
2657 _tree_set->insert(sb, tree);
2658 (*_blossom_data)[sb].pred = pred;
2659 (*_blossom_data)[sb].next =
2660 _graph.oppositeArc((*_blossom_data)[tb].next);
2662 pred = (*_blossom_data)[ub].next;
2664 (*_blossom_data)[tb].status = EVEN;
2665 matchedToEven(tb, tree);
2666 _tree_set->insert(tb, tree);
2667 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2670 (*_blossom_data)[subblossoms[id]].status = ODD;
2671 matchedToOdd(subblossoms[id]);
2672 _tree_set->insert(subblossoms[id], tree);
2673 (*_blossom_data)[subblossoms[id]].next = next;
2674 (*_blossom_data)[subblossoms[id]].pred = pred;
2678 for (int i = (ib + 1) % subblossoms.size();
2679 i != id; i = (i + 2) % subblossoms.size()) {
2680 int sb = subblossoms[i];
2681 int tb = subblossoms[(i + 1) % subblossoms.size()];
2682 (*_blossom_data)[sb].next =
2683 _graph.oppositeArc((*_blossom_data)[tb].next);
2686 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2687 int sb = subblossoms[i];
2688 int tb = subblossoms[(i + 1) % subblossoms.size()];
2689 int ub = subblossoms[(i + 2) % subblossoms.size()];
2691 (*_blossom_data)[sb].status = ODD;
2693 _tree_set->insert(sb, tree);
2694 (*_blossom_data)[sb].next = next;
2695 (*_blossom_data)[sb].pred =
2696 _graph.oppositeArc((*_blossom_data)[tb].next);
2698 (*_blossom_data)[tb].status = EVEN;
2699 matchedToEven(tb, tree);
2700 _tree_set->insert(tb, tree);
2701 (*_blossom_data)[tb].pred =
2702 (*_blossom_data)[tb].next =
2703 _graph.oppositeArc((*_blossom_data)[ub].next);
2704 next = (*_blossom_data)[ub].next;
2707 (*_blossom_data)[subblossoms[ib]].status = ODD;
2708 matchedToOdd(subblossoms[ib]);
2709 _tree_set->insert(subblossoms[ib], tree);
2710 (*_blossom_data)[subblossoms[ib]].next = next;
2711 (*_blossom_data)[subblossoms[ib]].pred = pred;
2713 _tree_set->erase(blossom);
2716 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2717 if (_blossom_set->trivial(blossom)) {
2718 int bi = (*_node_index)[base];
2719 Value pot = (*_node_data)[bi].pot;
2721 (*_matching)[base] = matching;
2722 _blossom_node_list.push_back(base);
2723 (*_node_potential)[base] = pot;
2726 Value pot = (*_blossom_data)[blossom].pot;
2727 int bn = _blossom_node_list.size();
2729 std::vector<int> subblossoms;
2730 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2731 int b = _blossom_set->find(base);
2733 for (int i = 0; i < int(subblossoms.size()); ++i) {
2734 if (subblossoms[i] == b) { ib = i; break; }
2737 for (int i = 1; i < int(subblossoms.size()); i += 2) {
2738 int sb = subblossoms[(ib + i) % subblossoms.size()];
2739 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2741 Arc m = (*_blossom_data)[tb].next;
2742 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2743 extractBlossom(tb, _graph.source(m), m);
2745 extractBlossom(subblossoms[ib], base, matching);
2747 int en = _blossom_node_list.size();
2749 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2753 void extractMatching() {
2754 std::vector<int> blossoms;
2755 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2756 blossoms.push_back(c);
2759 for (int i = 0; i < int(blossoms.size()); ++i) {
2761 Value offset = (*_blossom_data)[blossoms[i]].offset;
2762 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2763 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2764 n != INVALID; ++n) {
2765 (*_node_data)[(*_node_index)[n]].pot -= offset;
2768 Arc matching = (*_blossom_data)[blossoms[i]].next;
2769 Node base = _graph.source(matching);
2770 extractBlossom(blossoms[i], base, matching);
2776 /// \brief Constructor
2779 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2780 : _graph(graph), _weight(weight), _matching(0),
2781 _node_potential(0), _blossom_potential(), _blossom_node_list(),
2782 _node_num(0), _blossom_num(0),
2784 _blossom_index(0), _blossom_set(0), _blossom_data(0),
2785 _node_index(0), _node_heap_index(0), _node_data(0),
2786 _tree_set_index(0), _tree_set(0),
2788 _delta2_index(0), _delta2(0),
2789 _delta3_index(0), _delta3(0),
2790 _delta4_index(0), _delta4(0),
2794 ~MaxWeightedPerfectMatching() {
2795 destroyStructures();
2798 /// \name Execution Control
2799 /// The simplest way to execute the algorithm is to use the
2800 /// \ref run() member function.
2804 /// \brief Initialize the algorithm
2806 /// This function initializes the algorithm.
2810 for (ArcIt e(_graph); e != INVALID; ++e) {
2811 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2813 for (EdgeIt e(_graph); e != INVALID; ++e) {
2814 (*_delta3_index)[e] = _delta3->PRE_HEAP;
2816 for (int i = 0; i < _blossom_num; ++i) {
2817 (*_delta2_index)[i] = _delta2->PRE_HEAP;
2818 (*_delta4_index)[i] = _delta4->PRE_HEAP;
2822 for (NodeIt n(_graph); n != INVALID; ++n) {
2823 Value max = - std::numeric_limits<Value>::max();
2824 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2825 if (_graph.target(e) == n) continue;
2826 if ((dualScale * _weight[e]) / 2 > max) {
2827 max = (dualScale * _weight[e]) / 2;
2830 (*_node_index)[n] = index;
2831 (*_node_data)[index].pot = max;
2833 _blossom_set->insert(n, std::numeric_limits<Value>::max());
2835 _tree_set->insert(blossom);
2837 (*_blossom_data)[blossom].status = EVEN;
2838 (*_blossom_data)[blossom].pred = INVALID;
2839 (*_blossom_data)[blossom].next = INVALID;
2840 (*_blossom_data)[blossom].pot = 0;
2841 (*_blossom_data)[blossom].offset = 0;
2844 for (EdgeIt e(_graph); e != INVALID; ++e) {
2845 int si = (*_node_index)[_graph.u(e)];
2846 int ti = (*_node_index)[_graph.v(e)];
2847 if (_graph.u(e) != _graph.v(e)) {
2848 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2849 dualScale * _weight[e]) / 2);
2854 /// \brief Start the algorithm
2856 /// This function starts the algorithm.
2858 /// \pre \ref init() must be called before using this function.
2864 int unmatched = _node_num;
2865 while (unmatched > 0) {
2866 Value d2 = !_delta2->empty() ?
2867 _delta2->prio() : std::numeric_limits<Value>::max();
2869 Value d3 = !_delta3->empty() ?
2870 _delta3->prio() : std::numeric_limits<Value>::max();
2872 Value d4 = !_delta4->empty() ?
2873 _delta4->prio() : std::numeric_limits<Value>::max();
2875 _delta_sum = d3; OpType ot = D3;
2876 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
2877 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2879 if (_delta_sum == std::numeric_limits<Value>::max()) {
2886 int blossom = _delta2->top();
2887 Node n = _blossom_set->classTop(blossom);
2888 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2894 Edge e = _delta3->top();
2896 int left_blossom = _blossom_set->find(_graph.u(e));
2897 int right_blossom = _blossom_set->find(_graph.v(e));
2899 if (left_blossom == right_blossom) {
2902 int left_tree = _tree_set->find(left_blossom);
2903 int right_tree = _tree_set->find(right_blossom);
2905 if (left_tree == right_tree) {
2906 shrinkOnEdge(e, left_tree);
2914 splitBlossom(_delta4->top());
2922 /// \brief Run the algorithm.
2924 /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
2926 /// \note mwpm.run() is just a shortcut of the following code.
2938 /// \name Primal Solution
2939 /// Functions to get the primal solution, i.e. the maximum weighted
2940 /// perfect matching.\n
2941 /// Either \ref run() or \ref start() function should be called before
2946 /// \brief Return the weight of the matching.
2948 /// This function returns the weight of the found matching.
2950 /// \pre Either run() or start() must be called before using this function.
2951 Value matchingWeight() const {
2953 for (NodeIt n(_graph); n != INVALID; ++n) {
2954 if ((*_matching)[n] != INVALID) {
2955 sum += _weight[(*_matching)[n]];
2961 /// \brief Return \c true if the given edge is in the matching.
2963 /// This function returns \c true if the given edge is in the found
2966 /// \pre Either run() or start() must be called before using this function.
2967 bool matching(const Edge& edge) const {
2968 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
2971 /// \brief Return the matching arc (or edge) incident to the given node.
2973 /// This function returns the matching arc (or edge) incident to the
2974 /// given node in the found matching or \c INVALID if the node is
2975 /// not covered by the matching.
2977 /// \pre Either run() or start() must be called before using this function.
2978 Arc matching(const Node& node) const {
2979 return (*_matching)[node];
2982 /// \brief Return a const reference to the matching map.
2984 /// This function returns a const reference to a node map that stores
2985 /// the matching arc (or edge) incident to each node.
2986 const MatchingMap& matchingMap() const {
2990 /// \brief Return the mate of the given node.
2992 /// This function returns the mate of the given node in the found
2993 /// matching or \c INVALID if the node is not covered by the matching.
2995 /// \pre Either run() or start() must be called before using this function.
2996 Node mate(const Node& node) const {
2997 return _graph.target((*_matching)[node]);
3002 /// \name Dual Solution
3003 /// Functions to get the dual solution.\n
3004 /// Either \ref run() or \ref start() function should be called before
3009 /// \brief Return the value of the dual solution.
3011 /// This function returns the value of the dual solution.
3012 /// It should be equal to the primal value scaled by \ref dualScale
3015 /// \pre Either run() or start() must be called before using this function.
3016 Value dualValue() const {
3018 for (NodeIt n(_graph); n != INVALID; ++n) {
3019 sum += nodeValue(n);
3021 for (int i = 0; i < blossomNum(); ++i) {
3022 sum += blossomValue(i) * (blossomSize(i) / 2);
3027 /// \brief Return the dual value (potential) of the given node.
3029 /// This function returns the dual value (potential) of the given node.
3031 /// \pre Either run() or start() must be called before using this function.
3032 Value nodeValue(const Node& n) const {
3033 return (*_node_potential)[n];
3036 /// \brief Return the number of the blossoms in the basis.
3038 /// This function returns the number of the blossoms in the basis.
3040 /// \pre Either run() or start() must be called before using this function.
3042 int blossomNum() const {
3043 return _blossom_potential.size();
3046 /// \brief Return the number of the nodes in the given blossom.
3048 /// This function returns the number of the nodes in the given blossom.
3050 /// \pre Either run() or start() must be called before using this function.
3052 int blossomSize(int k) const {
3053 return _blossom_potential[k].end - _blossom_potential[k].begin;
3056 /// \brief Return the dual value (ptential) of the given blossom.
3058 /// This function returns the dual value (ptential) of the given blossom.
3060 /// \pre Either run() or start() must be called before using this function.
3061 Value blossomValue(int k) const {
3062 return _blossom_potential[k].value;
3065 /// \brief Iterator for obtaining the nodes of a blossom.
3067 /// This class provides an iterator for obtaining the nodes of the
3068 /// given blossom. It lists a subset of the nodes.
3069 /// Before using this iterator, you must allocate a
3070 /// MaxWeightedPerfectMatching class and execute it.
3074 /// \brief Constructor.
3076 /// Constructor to get the nodes of the given variable.
3078 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3079 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3080 /// must be called before initializing this iterator.
3081 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3082 : _algorithm(&algorithm)
3084 _index = _algorithm->_blossom_potential[variable].begin;
3085 _last = _algorithm->_blossom_potential[variable].end;
3088 /// \brief Conversion to \c Node.
3090 /// Conversion to \c Node.
3091 operator Node() const {
3092 return _algorithm->_blossom_node_list[_index];
3095 /// \brief Increment operator.
3097 /// Increment operator.
3098 BlossomIt& operator++() {
3103 /// \brief Validity checking
3105 /// This function checks whether the iterator is invalid.
3106 bool operator==(Invalid) const { return _index == _last; }
3108 /// \brief Validity checking
3110 /// This function checks whether the iterator is valid.
3111 bool operator!=(Invalid) const { return _index != _last; }
3114 const MaxWeightedPerfectMatching* _algorithm;
3123 } //END OF NAMESPACE LEMON
3125 #endif //LEMON_MATCHING_H