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_MAX_MATCHING_H
20 #define LEMON_MAX_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 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 by calling \c
57 /// decomposition() after running the algorithm.
59 /// \tparam GR The graph type the algorithm runs on.
60 template <typename GR>
64 /// The graph type of the algorithm
66 typedef typename Graph::template NodeMap<typename Graph::Arc>
69 ///\brief Status constants for Gallai-Edmonds decomposition.
71 ///These constants are used for indicating the Gallai-Edmonds
72 ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
73 ///induce a subgraph with factor-critical components, the nodes with
74 ///status \c ODD (or \c A) form the canonical barrier, and the nodes
75 ///with status \c MATCHED (or \c C) induce a subgraph having a
78 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.)
80 MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.)
82 ODD = -1, ///< = -1. (\c A is an alias for \c ODD.)
84 UNMATCHED = -2 ///< = -2.
87 typedef typename Graph::template NodeMap<Status> StatusMap;
91 TEMPLATE_GRAPH_TYPEDEFS(Graph);
93 typedef UnionFindEnum<IntNodeMap> BlossomSet;
94 typedef ExtendFindEnum<IntNodeMap> TreeSet;
95 typedef RangeMap<Node> NodeIntMap;
96 typedef MatchingMap EarMap;
97 typedef std::vector<Node> NodeQueue;
100 MatchingMap* _matching;
105 IntNodeMap* _blossom_set_index;
106 BlossomSet* _blossom_set;
107 NodeIntMap* _blossom_rep;
109 IntNodeMap* _tree_set_index;
112 NodeQueue _node_queue;
113 int _process, _postpone, _last;
119 void createStructures() {
120 _node_num = countNodes(_graph);
122 _matching = new MatchingMap(_graph);
125 _status = new StatusMap(_graph);
128 _ear = new EarMap(_graph);
131 _blossom_set_index = new IntNodeMap(_graph);
132 _blossom_set = new BlossomSet(*_blossom_set_index);
135 _blossom_rep = new NodeIntMap(_node_num);
138 _tree_set_index = new IntNodeMap(_graph);
139 _tree_set = new TreeSet(*_tree_set_index);
141 _node_queue.resize(_node_num);
144 void destroyStructures() {
156 delete _blossom_set_index;
162 delete _tree_set_index;
167 void processDense(const Node& n) {
168 _process = _postpone = _last = 0;
169 _node_queue[_last++] = n;
171 while (_process != _last) {
172 Node u = _node_queue[_process++];
173 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
174 Node v = _graph.target(a);
175 if ((*_status)[v] == MATCHED) {
177 } else if ((*_status)[v] == UNMATCHED) {
184 while (_postpone != _last) {
185 Node u = _node_queue[_postpone++];
187 for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
188 Node v = _graph.target(a);
190 if ((*_status)[v] == EVEN) {
191 if (_blossom_set->find(u) != _blossom_set->find(v)) {
196 while (_process != _last) {
197 Node w = _node_queue[_process++];
198 for (OutArcIt b(_graph, w); b != INVALID; ++b) {
199 Node x = _graph.target(b);
200 if ((*_status)[x] == MATCHED) {
202 } else if ((*_status)[x] == UNMATCHED) {
212 void processSparse(const Node& n) {
213 _process = _last = 0;
214 _node_queue[_last++] = n;
215 while (_process != _last) {
216 Node u = _node_queue[_process++];
217 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
218 Node v = _graph.target(a);
220 if ((*_status)[v] == EVEN) {
221 if (_blossom_set->find(u) != _blossom_set->find(v)) {
224 } else if ((*_status)[v] == MATCHED) {
226 } else if ((*_status)[v] == UNMATCHED) {
234 void shrinkOnEdge(const Edge& e) {
238 std::set<Node> left_set, right_set;
240 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
241 left_set.insert(left);
243 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
244 right_set.insert(right);
247 if ((*_matching)[left] == INVALID) break;
248 left = _graph.target((*_matching)[left]);
249 left = (*_blossom_rep)[_blossom_set->
250 find(_graph.target((*_ear)[left]))];
251 if (right_set.find(left) != right_set.end()) {
255 left_set.insert(left);
257 if ((*_matching)[right] == INVALID) break;
258 right = _graph.target((*_matching)[right]);
259 right = (*_blossom_rep)[_blossom_set->
260 find(_graph.target((*_ear)[right]))];
261 if (left_set.find(right) != left_set.end()) {
265 right_set.insert(right);
268 if (nca == INVALID) {
269 if ((*_matching)[left] == INVALID) {
271 while (left_set.find(nca) == left_set.end()) {
272 nca = _graph.target((*_matching)[nca]);
273 nca =(*_blossom_rep)[_blossom_set->
274 find(_graph.target((*_ear)[nca]))];
278 while (right_set.find(nca) == right_set.end()) {
279 nca = _graph.target((*_matching)[nca]);
280 nca = (*_blossom_rep)[_blossom_set->
281 find(_graph.target((*_ear)[nca]))];
289 Node node = _graph.u(e);
290 Arc arc = _graph.direct(e, true);
291 Node base = (*_blossom_rep)[_blossom_set->find(node)];
293 while (base != nca) {
298 n = _graph.target((*_matching)[n]);
300 n = _graph.target(a);
301 (*_ear)[n] = _graph.oppositeArc(a);
303 node = _graph.target((*_matching)[base]);
304 _tree_set->erase(base);
305 _tree_set->erase(node);
306 _blossom_set->insert(node, _blossom_set->find(base));
307 (*_status)[node] = EVEN;
308 _node_queue[_last++] = node;
309 arc = _graph.oppositeArc((*_ear)[node]);
310 node = _graph.target((*_ear)[node]);
311 base = (*_blossom_rep)[_blossom_set->find(node)];
312 _blossom_set->join(_graph.target(arc), base);
316 (*_blossom_rep)[_blossom_set->find(nca)] = nca;
320 Node node = _graph.v(e);
321 Arc arc = _graph.direct(e, false);
322 Node base = (*_blossom_rep)[_blossom_set->find(node)];
324 while (base != nca) {
329 n = _graph.target((*_matching)[n]);
331 n = _graph.target(a);
332 (*_ear)[n] = _graph.oppositeArc(a);
334 node = _graph.target((*_matching)[base]);
335 _tree_set->erase(base);
336 _tree_set->erase(node);
337 _blossom_set->insert(node, _blossom_set->find(base));
338 (*_status)[node] = EVEN;
339 _node_queue[_last++] = node;
340 arc = _graph.oppositeArc((*_ear)[node]);
341 node = _graph.target((*_ear)[node]);
342 base = (*_blossom_rep)[_blossom_set->find(node)];
343 _blossom_set->join(_graph.target(arc), base);
347 (*_blossom_rep)[_blossom_set->find(nca)] = nca;
350 void extendOnArc(const Arc& a) {
351 Node base = _graph.source(a);
352 Node odd = _graph.target(a);
354 (*_ear)[odd] = _graph.oppositeArc(a);
355 Node even = _graph.target((*_matching)[odd]);
356 (*_blossom_rep)[_blossom_set->insert(even)] = even;
357 (*_status)[odd] = ODD;
358 (*_status)[even] = EVEN;
359 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
360 _tree_set->insert(odd, tree);
361 _tree_set->insert(even, tree);
362 _node_queue[_last++] = even;
366 void augmentOnArc(const Arc& a) {
367 Node even = _graph.source(a);
368 Node odd = _graph.target(a);
370 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
372 (*_matching)[odd] = _graph.oppositeArc(a);
373 (*_status)[odd] = MATCHED;
375 Arc arc = (*_matching)[even];
376 (*_matching)[even] = a;
378 while (arc != INVALID) {
379 odd = _graph.target(arc);
381 even = _graph.target(arc);
382 (*_matching)[odd] = arc;
383 arc = (*_matching)[even];
384 (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
387 for (typename TreeSet::ItemIt it(*_tree_set, tree);
388 it != INVALID; ++it) {
389 if ((*_status)[it] == ODD) {
390 (*_status)[it] = MATCHED;
392 int blossom = _blossom_set->find(it);
393 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
394 jt != INVALID; ++jt) {
395 (*_status)[jt] = MATCHED;
397 _blossom_set->eraseClass(blossom);
400 _tree_set->eraseClass(tree);
406 /// \brief Constructor
409 MaxMatching(const Graph& graph)
410 : _graph(graph), _matching(0), _status(0), _ear(0),
411 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
412 _tree_set_index(0), _tree_set(0) {}
418 /// \name Execution Control
419 /// The simplest way to execute the algorithm is to use the
420 /// \c run() member function.\n
421 /// If you need better control on the execution, you have to call
422 /// one of the functions \ref init(), \ref greedyInit() or
423 /// \ref matchingInit() first, then you can start the algorithm with
424 /// \ref startSparse() or \ref startDense().
428 /// \brief Set the initial matching to the empty matching.
430 /// This function sets the initial matching to the empty matching.
433 for(NodeIt n(_graph); n != INVALID; ++n) {
434 (*_matching)[n] = INVALID;
435 (*_status)[n] = UNMATCHED;
439 /// \brief Find an initial matching in a greedy way.
441 /// This function finds an initial matching in a greedy way.
444 for (NodeIt n(_graph); n != INVALID; ++n) {
445 (*_matching)[n] = INVALID;
446 (*_status)[n] = UNMATCHED;
448 for (NodeIt n(_graph); n != INVALID; ++n) {
449 if ((*_matching)[n] == INVALID) {
450 for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
451 Node v = _graph.target(a);
452 if ((*_matching)[v] == INVALID && v != n) {
454 (*_status)[n] = MATCHED;
455 (*_matching)[v] = _graph.oppositeArc(a);
456 (*_status)[v] = MATCHED;
465 /// \brief Initialize the matching from a map.
467 /// This function initializes the matching from a \c bool valued edge
468 /// map. This map should have the property that there are no two incident
469 /// edges with \c true value, i.e. it really contains a matching.
470 /// \return \c true if the map contains a matching.
471 template <typename MatchingMap>
472 bool matchingInit(const MatchingMap& matching) {
475 for (NodeIt n(_graph); n != INVALID; ++n) {
476 (*_matching)[n] = INVALID;
477 (*_status)[n] = UNMATCHED;
479 for(EdgeIt e(_graph); e!=INVALID; ++e) {
482 Node u = _graph.u(e);
483 if ((*_matching)[u] != INVALID) return false;
484 (*_matching)[u] = _graph.direct(e, true);
485 (*_status)[u] = MATCHED;
487 Node v = _graph.v(e);
488 if ((*_matching)[v] != INVALID) return false;
489 (*_matching)[v] = _graph.direct(e, false);
490 (*_status)[v] = MATCHED;
496 /// \brief Start Edmonds' algorithm
498 /// This function runs the original Edmonds' algorithm.
500 /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be
501 /// called before using this function.
503 for(NodeIt n(_graph); n != INVALID; ++n) {
504 if ((*_status)[n] == UNMATCHED) {
505 (*_blossom_rep)[_blossom_set->insert(n)] = n;
506 _tree_set->insert(n);
507 (*_status)[n] = EVEN;
513 /// \brief Start Edmonds' algorithm with a heuristic improvement
516 /// This function runs Edmonds' algorithm with a heuristic of postponing
517 /// shrinks, therefore resulting in a faster algorithm for dense graphs.
519 /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be
520 /// called before using this function.
522 for(NodeIt n(_graph); n != INVALID; ++n) {
523 if ((*_status)[n] == UNMATCHED) {
524 (*_blossom_rep)[_blossom_set->insert(n)] = n;
525 _tree_set->insert(n);
526 (*_status)[n] = EVEN;
533 /// \brief Run Edmonds' algorithm
535 /// This function runs Edmonds' algorithm. An additional heuristic of
536 /// postponing shrinks is used for relatively dense graphs
537 /// (for which <tt>m>=2*n</tt> holds).
539 if (countEdges(_graph) < 2 * countNodes(_graph)) {
550 /// \name Primal Solution
551 /// Functions to get the primal solution, i.e. the maximum matching.
555 /// \brief Return the size (cardinality) of the matching.
557 /// This function returns the size (cardinality) of the current matching.
558 /// After run() it returns the size of the maximum matching in the graph.
559 int matchingSize() const {
561 for (NodeIt n(_graph); n != INVALID; ++n) {
562 if ((*_matching)[n] != INVALID) {
569 /// \brief Return \c true if the given edge is in the matching.
571 /// This function returns \c true if the given edge is in the current
573 bool matching(const Edge& edge) const {
574 return edge == (*_matching)[_graph.u(edge)];
577 /// \brief Return the matching arc (or edge) incident to the given node.
579 /// This function returns the matching arc (or edge) incident to the
580 /// given node in the current matching or \c INVALID if the node is
581 /// not covered by the matching.
582 Arc matching(const Node& n) const {
583 return (*_matching)[n];
586 /// \brief Return the mate of the given node.
588 /// This function returns the mate of the given node in the current
589 /// matching or \c INVALID if the node is not covered by the matching.
590 Node mate(const Node& n) const {
591 return (*_matching)[n] != INVALID ?
592 _graph.target((*_matching)[n]) : INVALID;
597 /// \name Dual Solution
598 /// Functions to get the dual solution, i.e. the Gallai-Edmonds
603 /// \brief Return the status of the given node in the Edmonds-Gallai
606 /// This function returns the \ref Status "status" of the given node
607 /// in the Edmonds-Gallai decomposition.
608 Status decomposition(const Node& n) const {
609 return (*_status)[n];
612 /// \brief Return \c true if the given node is in the barrier.
614 /// This function returns \c true if the given node is in the barrier.
615 bool barrier(const Node& n) const {
616 return (*_status)[n] == ODD;
623 /// \ingroup matching
625 /// \brief Weighted matching in general graphs
627 /// This class provides an efficient implementation of Edmond's
628 /// maximum weighted matching algorithm. The implementation is based
629 /// on extensive use of priority queues and provides
630 /// \f$O(nm\log n)\f$ time complexity.
632 /// The maximum weighted matching problem is to find a subset of the
633 /// edges in an undirected graph with maximum overall weight for which
634 /// each node has at most one incident edge.
635 /// It can be formulated with the following linear program.
636 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
637 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
638 \quad \forall B\in\mathcal{O}\f] */
639 /// \f[x_e \ge 0\quad \forall e\in E\f]
640 /// \f[\max \sum_{e\in E}x_ew_e\f]
641 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
642 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
643 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
644 /// subsets of the nodes.
646 /// The algorithm calculates an optimal matching and a proof of the
647 /// optimality. The solution of the dual problem can be used to check
648 /// the result of the algorithm. The dual linear problem is the
650 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
651 z_B \ge w_{uv} \quad \forall uv\in E\f] */
652 /// \f[y_u \ge 0 \quad \forall u \in V\f]
653 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
654 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
655 \frac{\vert B \vert - 1}{2}z_B\f] */
657 /// The algorithm can be executed with the run() function.
658 /// After it the matching (the primal solution) and the dual solution
659 /// can be obtained using the query functions and the
660 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
661 /// which is able to iterate on the nodes of a blossom.
662 /// If the value type is integer, then the dual solution is multiplied
663 /// by \ref MaxWeightedMatching::dualScale "4".
665 /// \tparam GR The graph type the algorithm runs on.
666 /// \tparam WM The type edge weight map. The default type is
667 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
669 template <typename GR, typename WM>
671 template <typename GR,
672 typename WM = typename GR::template EdgeMap<int> >
674 class MaxWeightedMatching {
677 /// The graph type of the algorithm
679 /// The type of the edge weight map
680 typedef WM WeightMap;
681 /// The value type of the edge weights
682 typedef typename WeightMap::Value Value;
684 typedef typename Graph::template NodeMap<typename Graph::Arc>
687 /// \brief Scaling factor for dual solution
689 /// Scaling factor for dual solution. It is equal to 4 or 1
690 /// according to the value type.
691 static const int dualScale =
692 std::numeric_limits<Value>::is_integer ? 4 : 1;
696 TEMPLATE_GRAPH_TYPEDEFS(Graph);
698 typedef typename Graph::template NodeMap<Value> NodePotential;
699 typedef std::vector<Node> BlossomNodeList;
701 struct BlossomVariable {
705 BlossomVariable(int _begin, int _end, Value _value)
706 : begin(_begin), end(_end), value(_value) {}
710 typedef std::vector<BlossomVariable> BlossomPotential;
713 const WeightMap& _weight;
715 MatchingMap* _matching;
717 NodePotential* _node_potential;
719 BlossomPotential _blossom_potential;
720 BlossomNodeList _blossom_node_list;
725 typedef RangeMap<int> IntIntMap;
728 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
731 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
740 IntNodeMap *_blossom_index;
741 BlossomSet *_blossom_set;
742 RangeMap<BlossomData>* _blossom_data;
744 IntNodeMap *_node_index;
745 IntArcMap *_node_heap_index;
749 NodeData(IntArcMap& node_heap_index)
750 : heap(node_heap_index) {}
754 BinHeap<Value, IntArcMap> heap;
755 std::map<int, Arc> heap_index;
760 RangeMap<NodeData>* _node_data;
762 typedef ExtendFindEnum<IntIntMap> TreeSet;
764 IntIntMap *_tree_set_index;
767 IntNodeMap *_delta1_index;
768 BinHeap<Value, IntNodeMap> *_delta1;
770 IntIntMap *_delta2_index;
771 BinHeap<Value, IntIntMap> *_delta2;
773 IntEdgeMap *_delta3_index;
774 BinHeap<Value, IntEdgeMap> *_delta3;
776 IntIntMap *_delta4_index;
777 BinHeap<Value, IntIntMap> *_delta4;
781 void createStructures() {
782 _node_num = countNodes(_graph);
783 _blossom_num = _node_num * 3 / 2;
786 _matching = new MatchingMap(_graph);
788 if (!_node_potential) {
789 _node_potential = new NodePotential(_graph);
792 _blossom_index = new IntNodeMap(_graph);
793 _blossom_set = new BlossomSet(*_blossom_index);
794 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
798 _node_index = new IntNodeMap(_graph);
799 _node_heap_index = new IntArcMap(_graph);
800 _node_data = new RangeMap<NodeData>(_node_num,
801 NodeData(*_node_heap_index));
805 _tree_set_index = new IntIntMap(_blossom_num);
806 _tree_set = new TreeSet(*_tree_set_index);
809 _delta1_index = new IntNodeMap(_graph);
810 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
813 _delta2_index = new IntIntMap(_blossom_num);
814 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
817 _delta3_index = new IntEdgeMap(_graph);
818 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
821 _delta4_index = new IntIntMap(_blossom_num);
822 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
826 void destroyStructures() {
827 _node_num = countNodes(_graph);
828 _blossom_num = _node_num * 3 / 2;
833 if (_node_potential) {
834 delete _node_potential;
837 delete _blossom_index;
839 delete _blossom_data;
844 delete _node_heap_index;
849 delete _tree_set_index;
853 delete _delta1_index;
857 delete _delta2_index;
861 delete _delta3_index;
865 delete _delta4_index;
870 void matchedToEven(int blossom, int tree) {
871 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
872 _delta2->erase(blossom);
875 if (!_blossom_set->trivial(blossom)) {
876 (*_blossom_data)[blossom].pot -=
877 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
880 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
883 _blossom_set->increase(n, std::numeric_limits<Value>::max());
884 int ni = (*_node_index)[n];
886 (*_node_data)[ni].heap.clear();
887 (*_node_data)[ni].heap_index.clear();
889 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
891 _delta1->push(n, (*_node_data)[ni].pot);
893 for (InArcIt e(_graph, n); e != INVALID; ++e) {
894 Node v = _graph.source(e);
895 int vb = _blossom_set->find(v);
896 int vi = (*_node_index)[v];
898 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
899 dualScale * _weight[e];
901 if ((*_blossom_data)[vb].status == EVEN) {
902 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
903 _delta3->push(e, rw / 2);
905 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
906 if (_delta3->state(e) != _delta3->IN_HEAP) {
907 _delta3->push(e, rw);
910 typename std::map<int, Arc>::iterator it =
911 (*_node_data)[vi].heap_index.find(tree);
913 if (it != (*_node_data)[vi].heap_index.end()) {
914 if ((*_node_data)[vi].heap[it->second] > rw) {
915 (*_node_data)[vi].heap.replace(it->second, e);
916 (*_node_data)[vi].heap.decrease(e, rw);
920 (*_node_data)[vi].heap.push(e, rw);
921 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
924 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
925 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
927 if ((*_blossom_data)[vb].status == MATCHED) {
928 if (_delta2->state(vb) != _delta2->IN_HEAP) {
929 _delta2->push(vb, _blossom_set->classPrio(vb) -
930 (*_blossom_data)[vb].offset);
931 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
932 (*_blossom_data)[vb].offset){
933 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
934 (*_blossom_data)[vb].offset);
941 (*_blossom_data)[blossom].offset = 0;
944 void matchedToOdd(int blossom) {
945 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
946 _delta2->erase(blossom);
948 (*_blossom_data)[blossom].offset += _delta_sum;
949 if (!_blossom_set->trivial(blossom)) {
950 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
951 (*_blossom_data)[blossom].offset);
955 void evenToMatched(int blossom, int tree) {
956 if (!_blossom_set->trivial(blossom)) {
957 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
960 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
962 int ni = (*_node_index)[n];
963 (*_node_data)[ni].pot -= _delta_sum;
967 for (InArcIt e(_graph, n); e != INVALID; ++e) {
968 Node v = _graph.source(e);
969 int vb = _blossom_set->find(v);
970 int vi = (*_node_index)[v];
972 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
973 dualScale * _weight[e];
976 if (_delta3->state(e) == _delta3->IN_HEAP) {
979 } else if ((*_blossom_data)[vb].status == EVEN) {
981 if (_delta3->state(e) == _delta3->IN_HEAP) {
985 int vt = _tree_set->find(vb);
989 Arc r = _graph.oppositeArc(e);
991 typename std::map<int, Arc>::iterator it =
992 (*_node_data)[ni].heap_index.find(vt);
994 if (it != (*_node_data)[ni].heap_index.end()) {
995 if ((*_node_data)[ni].heap[it->second] > rw) {
996 (*_node_data)[ni].heap.replace(it->second, r);
997 (*_node_data)[ni].heap.decrease(r, rw);
1001 (*_node_data)[ni].heap.push(r, rw);
1002 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1005 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1006 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1008 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1009 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1010 (*_blossom_data)[blossom].offset);
1011 } else if ((*_delta2)[blossom] >
1012 _blossom_set->classPrio(blossom) -
1013 (*_blossom_data)[blossom].offset){
1014 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1015 (*_blossom_data)[blossom].offset);
1020 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1021 if (_delta3->state(e) == _delta3->IN_HEAP) {
1026 typename std::map<int, Arc>::iterator it =
1027 (*_node_data)[vi].heap_index.find(tree);
1029 if (it != (*_node_data)[vi].heap_index.end()) {
1030 (*_node_data)[vi].heap.erase(it->second);
1031 (*_node_data)[vi].heap_index.erase(it);
1032 if ((*_node_data)[vi].heap.empty()) {
1033 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1034 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1035 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1038 if ((*_blossom_data)[vb].status == MATCHED) {
1039 if (_blossom_set->classPrio(vb) ==
1040 std::numeric_limits<Value>::max()) {
1042 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1043 (*_blossom_data)[vb].offset) {
1044 _delta2->increase(vb, _blossom_set->classPrio(vb) -
1045 (*_blossom_data)[vb].offset);
1054 void oddToMatched(int blossom) {
1055 (*_blossom_data)[blossom].offset -= _delta_sum;
1057 if (_blossom_set->classPrio(blossom) !=
1058 std::numeric_limits<Value>::max()) {
1059 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1060 (*_blossom_data)[blossom].offset);
1063 if (!_blossom_set->trivial(blossom)) {
1064 _delta4->erase(blossom);
1068 void oddToEven(int blossom, int tree) {
1069 if (!_blossom_set->trivial(blossom)) {
1070 _delta4->erase(blossom);
1071 (*_blossom_data)[blossom].pot -=
1072 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1075 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1076 n != INVALID; ++n) {
1077 int ni = (*_node_index)[n];
1079 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1081 (*_node_data)[ni].heap.clear();
1082 (*_node_data)[ni].heap_index.clear();
1083 (*_node_data)[ni].pot +=
1084 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1086 _delta1->push(n, (*_node_data)[ni].pot);
1088 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1089 Node v = _graph.source(e);
1090 int vb = _blossom_set->find(v);
1091 int vi = (*_node_index)[v];
1093 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1094 dualScale * _weight[e];
1096 if ((*_blossom_data)[vb].status == EVEN) {
1097 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1098 _delta3->push(e, rw / 2);
1100 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1101 if (_delta3->state(e) != _delta3->IN_HEAP) {
1102 _delta3->push(e, rw);
1106 typename std::map<int, Arc>::iterator it =
1107 (*_node_data)[vi].heap_index.find(tree);
1109 if (it != (*_node_data)[vi].heap_index.end()) {
1110 if ((*_node_data)[vi].heap[it->second] > rw) {
1111 (*_node_data)[vi].heap.replace(it->second, e);
1112 (*_node_data)[vi].heap.decrease(e, rw);
1116 (*_node_data)[vi].heap.push(e, rw);
1117 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1120 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1121 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1123 if ((*_blossom_data)[vb].status == MATCHED) {
1124 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1125 _delta2->push(vb, _blossom_set->classPrio(vb) -
1126 (*_blossom_data)[vb].offset);
1127 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1128 (*_blossom_data)[vb].offset) {
1129 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1130 (*_blossom_data)[vb].offset);
1137 (*_blossom_data)[blossom].offset = 0;
1141 void matchedToUnmatched(int blossom) {
1142 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1143 _delta2->erase(blossom);
1146 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1147 n != INVALID; ++n) {
1148 int ni = (*_node_index)[n];
1150 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1152 (*_node_data)[ni].heap.clear();
1153 (*_node_data)[ni].heap_index.clear();
1155 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1156 Node v = _graph.target(e);
1157 int vb = _blossom_set->find(v);
1158 int vi = (*_node_index)[v];
1160 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1161 dualScale * _weight[e];
1163 if ((*_blossom_data)[vb].status == EVEN) {
1164 if (_delta3->state(e) != _delta3->IN_HEAP) {
1165 _delta3->push(e, rw);
1172 void unmatchedToMatched(int blossom) {
1173 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1174 n != INVALID; ++n) {
1175 int ni = (*_node_index)[n];
1177 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1178 Node v = _graph.source(e);
1179 int vb = _blossom_set->find(v);
1180 int vi = (*_node_index)[v];
1182 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1183 dualScale * _weight[e];
1185 if (vb == blossom) {
1186 if (_delta3->state(e) == _delta3->IN_HEAP) {
1189 } else if ((*_blossom_data)[vb].status == EVEN) {
1191 if (_delta3->state(e) == _delta3->IN_HEAP) {
1195 int vt = _tree_set->find(vb);
1197 Arc r = _graph.oppositeArc(e);
1199 typename std::map<int, Arc>::iterator it =
1200 (*_node_data)[ni].heap_index.find(vt);
1202 if (it != (*_node_data)[ni].heap_index.end()) {
1203 if ((*_node_data)[ni].heap[it->second] > rw) {
1204 (*_node_data)[ni].heap.replace(it->second, r);
1205 (*_node_data)[ni].heap.decrease(r, rw);
1209 (*_node_data)[ni].heap.push(r, rw);
1210 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1213 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1214 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1216 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1217 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1218 (*_blossom_data)[blossom].offset);
1219 } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1220 (*_blossom_data)[blossom].offset){
1221 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1222 (*_blossom_data)[blossom].offset);
1226 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1227 if (_delta3->state(e) == _delta3->IN_HEAP) {
1235 void alternatePath(int even, int tree) {
1238 evenToMatched(even, tree);
1239 (*_blossom_data)[even].status = MATCHED;
1241 while ((*_blossom_data)[even].pred != INVALID) {
1242 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1243 (*_blossom_data)[odd].status = MATCHED;
1245 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1247 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1248 (*_blossom_data)[even].status = MATCHED;
1249 evenToMatched(even, tree);
1250 (*_blossom_data)[even].next =
1251 _graph.oppositeArc((*_blossom_data)[odd].pred);
1256 void destroyTree(int tree) {
1257 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1258 if ((*_blossom_data)[b].status == EVEN) {
1259 (*_blossom_data)[b].status = MATCHED;
1260 evenToMatched(b, tree);
1261 } else if ((*_blossom_data)[b].status == ODD) {
1262 (*_blossom_data)[b].status = MATCHED;
1266 _tree_set->eraseClass(tree);
1270 void unmatchNode(const Node& node) {
1271 int blossom = _blossom_set->find(node);
1272 int tree = _tree_set->find(blossom);
1274 alternatePath(blossom, tree);
1277 (*_blossom_data)[blossom].status = UNMATCHED;
1278 (*_blossom_data)[blossom].base = node;
1279 matchedToUnmatched(blossom);
1283 void augmentOnEdge(const Edge& edge) {
1285 int left = _blossom_set->find(_graph.u(edge));
1286 int right = _blossom_set->find(_graph.v(edge));
1288 if ((*_blossom_data)[left].status == EVEN) {
1289 int left_tree = _tree_set->find(left);
1290 alternatePath(left, left_tree);
1291 destroyTree(left_tree);
1293 (*_blossom_data)[left].status = MATCHED;
1294 unmatchedToMatched(left);
1297 if ((*_blossom_data)[right].status == EVEN) {
1298 int right_tree = _tree_set->find(right);
1299 alternatePath(right, right_tree);
1300 destroyTree(right_tree);
1302 (*_blossom_data)[right].status = MATCHED;
1303 unmatchedToMatched(right);
1306 (*_blossom_data)[left].next = _graph.direct(edge, true);
1307 (*_blossom_data)[right].next = _graph.direct(edge, false);
1310 void extendOnArc(const Arc& arc) {
1311 int base = _blossom_set->find(_graph.target(arc));
1312 int tree = _tree_set->find(base);
1314 int odd = _blossom_set->find(_graph.source(arc));
1315 _tree_set->insert(odd, tree);
1316 (*_blossom_data)[odd].status = ODD;
1318 (*_blossom_data)[odd].pred = arc;
1320 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1321 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1322 _tree_set->insert(even, tree);
1323 (*_blossom_data)[even].status = EVEN;
1324 matchedToEven(even, tree);
1327 void shrinkOnEdge(const Edge& edge, int tree) {
1329 std::vector<int> left_path, right_path;
1332 std::set<int> left_set, right_set;
1333 int left = _blossom_set->find(_graph.u(edge));
1334 left_path.push_back(left);
1335 left_set.insert(left);
1337 int right = _blossom_set->find(_graph.v(edge));
1338 right_path.push_back(right);
1339 right_set.insert(right);
1343 if ((*_blossom_data)[left].pred == INVALID) break;
1346 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1347 left_path.push_back(left);
1349 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1350 left_path.push_back(left);
1352 left_set.insert(left);
1354 if (right_set.find(left) != right_set.end()) {
1359 if ((*_blossom_data)[right].pred == INVALID) break;
1362 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1363 right_path.push_back(right);
1365 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1366 right_path.push_back(right);
1368 right_set.insert(right);
1370 if (left_set.find(right) != left_set.end()) {
1378 if ((*_blossom_data)[left].pred == INVALID) {
1380 while (left_set.find(nca) == left_set.end()) {
1382 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1383 right_path.push_back(nca);
1385 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1386 right_path.push_back(nca);
1390 while (right_set.find(nca) == right_set.end()) {
1392 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1393 left_path.push_back(nca);
1395 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1396 left_path.push_back(nca);
1402 std::vector<int> subblossoms;
1405 prev = _graph.direct(edge, true);
1406 for (int i = 0; left_path[i] != nca; i += 2) {
1407 subblossoms.push_back(left_path[i]);
1408 (*_blossom_data)[left_path[i]].next = prev;
1409 _tree_set->erase(left_path[i]);
1411 subblossoms.push_back(left_path[i + 1]);
1412 (*_blossom_data)[left_path[i + 1]].status = EVEN;
1413 oddToEven(left_path[i + 1], tree);
1414 _tree_set->erase(left_path[i + 1]);
1415 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1419 while (right_path[k] != nca) ++k;
1421 subblossoms.push_back(nca);
1422 (*_blossom_data)[nca].next = prev;
1424 for (int i = k - 2; i >= 0; i -= 2) {
1425 subblossoms.push_back(right_path[i + 1]);
1426 (*_blossom_data)[right_path[i + 1]].status = EVEN;
1427 oddToEven(right_path[i + 1], tree);
1428 _tree_set->erase(right_path[i + 1]);
1430 (*_blossom_data)[right_path[i + 1]].next =
1431 (*_blossom_data)[right_path[i + 1]].pred;
1433 subblossoms.push_back(right_path[i]);
1434 _tree_set->erase(right_path[i]);
1438 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1440 for (int i = 0; i < int(subblossoms.size()); ++i) {
1441 if (!_blossom_set->trivial(subblossoms[i])) {
1442 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1444 (*_blossom_data)[subblossoms[i]].status = MATCHED;
1447 (*_blossom_data)[surface].pot = -2 * _delta_sum;
1448 (*_blossom_data)[surface].offset = 0;
1449 (*_blossom_data)[surface].status = EVEN;
1450 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1451 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1453 _tree_set->insert(surface, tree);
1454 _tree_set->erase(nca);
1457 void splitBlossom(int blossom) {
1458 Arc next = (*_blossom_data)[blossom].next;
1459 Arc pred = (*_blossom_data)[blossom].pred;
1461 int tree = _tree_set->find(blossom);
1463 (*_blossom_data)[blossom].status = MATCHED;
1464 oddToMatched(blossom);
1465 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1466 _delta2->erase(blossom);
1469 std::vector<int> subblossoms;
1470 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1472 Value offset = (*_blossom_data)[blossom].offset;
1473 int b = _blossom_set->find(_graph.source(pred));
1474 int d = _blossom_set->find(_graph.source(next));
1476 int ib = -1, id = -1;
1477 for (int i = 0; i < int(subblossoms.size()); ++i) {
1478 if (subblossoms[i] == b) ib = i;
1479 if (subblossoms[i] == d) id = i;
1481 (*_blossom_data)[subblossoms[i]].offset = offset;
1482 if (!_blossom_set->trivial(subblossoms[i])) {
1483 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1485 if (_blossom_set->classPrio(subblossoms[i]) !=
1486 std::numeric_limits<Value>::max()) {
1487 _delta2->push(subblossoms[i],
1488 _blossom_set->classPrio(subblossoms[i]) -
1489 (*_blossom_data)[subblossoms[i]].offset);
1493 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1494 for (int i = (id + 1) % subblossoms.size();
1495 i != ib; i = (i + 2) % subblossoms.size()) {
1496 int sb = subblossoms[i];
1497 int tb = subblossoms[(i + 1) % subblossoms.size()];
1498 (*_blossom_data)[sb].next =
1499 _graph.oppositeArc((*_blossom_data)[tb].next);
1502 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1503 int sb = subblossoms[i];
1504 int tb = subblossoms[(i + 1) % subblossoms.size()];
1505 int ub = subblossoms[(i + 2) % subblossoms.size()];
1507 (*_blossom_data)[sb].status = ODD;
1509 _tree_set->insert(sb, tree);
1510 (*_blossom_data)[sb].pred = pred;
1511 (*_blossom_data)[sb].next =
1512 _graph.oppositeArc((*_blossom_data)[tb].next);
1514 pred = (*_blossom_data)[ub].next;
1516 (*_blossom_data)[tb].status = EVEN;
1517 matchedToEven(tb, tree);
1518 _tree_set->insert(tb, tree);
1519 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1522 (*_blossom_data)[subblossoms[id]].status = ODD;
1523 matchedToOdd(subblossoms[id]);
1524 _tree_set->insert(subblossoms[id], tree);
1525 (*_blossom_data)[subblossoms[id]].next = next;
1526 (*_blossom_data)[subblossoms[id]].pred = pred;
1530 for (int i = (ib + 1) % subblossoms.size();
1531 i != id; i = (i + 2) % subblossoms.size()) {
1532 int sb = subblossoms[i];
1533 int tb = subblossoms[(i + 1) % subblossoms.size()];
1534 (*_blossom_data)[sb].next =
1535 _graph.oppositeArc((*_blossom_data)[tb].next);
1538 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1539 int sb = subblossoms[i];
1540 int tb = subblossoms[(i + 1) % subblossoms.size()];
1541 int ub = subblossoms[(i + 2) % subblossoms.size()];
1543 (*_blossom_data)[sb].status = ODD;
1545 _tree_set->insert(sb, tree);
1546 (*_blossom_data)[sb].next = next;
1547 (*_blossom_data)[sb].pred =
1548 _graph.oppositeArc((*_blossom_data)[tb].next);
1550 (*_blossom_data)[tb].status = EVEN;
1551 matchedToEven(tb, tree);
1552 _tree_set->insert(tb, tree);
1553 (*_blossom_data)[tb].pred =
1554 (*_blossom_data)[tb].next =
1555 _graph.oppositeArc((*_blossom_data)[ub].next);
1556 next = (*_blossom_data)[ub].next;
1559 (*_blossom_data)[subblossoms[ib]].status = ODD;
1560 matchedToOdd(subblossoms[ib]);
1561 _tree_set->insert(subblossoms[ib], tree);
1562 (*_blossom_data)[subblossoms[ib]].next = next;
1563 (*_blossom_data)[subblossoms[ib]].pred = pred;
1565 _tree_set->erase(blossom);
1568 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1569 if (_blossom_set->trivial(blossom)) {
1570 int bi = (*_node_index)[base];
1571 Value pot = (*_node_data)[bi].pot;
1573 (*_matching)[base] = matching;
1574 _blossom_node_list.push_back(base);
1575 (*_node_potential)[base] = pot;
1578 Value pot = (*_blossom_data)[blossom].pot;
1579 int bn = _blossom_node_list.size();
1581 std::vector<int> subblossoms;
1582 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1583 int b = _blossom_set->find(base);
1585 for (int i = 0; i < int(subblossoms.size()); ++i) {
1586 if (subblossoms[i] == b) { ib = i; break; }
1589 for (int i = 1; i < int(subblossoms.size()); i += 2) {
1590 int sb = subblossoms[(ib + i) % subblossoms.size()];
1591 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1593 Arc m = (*_blossom_data)[tb].next;
1594 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1595 extractBlossom(tb, _graph.source(m), m);
1597 extractBlossom(subblossoms[ib], base, matching);
1599 int en = _blossom_node_list.size();
1601 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1605 void extractMatching() {
1606 std::vector<int> blossoms;
1607 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1608 blossoms.push_back(c);
1611 for (int i = 0; i < int(blossoms.size()); ++i) {
1612 if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1614 Value offset = (*_blossom_data)[blossoms[i]].offset;
1615 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1616 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1617 n != INVALID; ++n) {
1618 (*_node_data)[(*_node_index)[n]].pot -= offset;
1621 Arc matching = (*_blossom_data)[blossoms[i]].next;
1622 Node base = _graph.source(matching);
1623 extractBlossom(blossoms[i], base, matching);
1625 Node base = (*_blossom_data)[blossoms[i]].base;
1626 extractBlossom(blossoms[i], base, INVALID);
1633 /// \brief Constructor
1636 MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1637 : _graph(graph), _weight(weight), _matching(0),
1638 _node_potential(0), _blossom_potential(), _blossom_node_list(),
1639 _node_num(0), _blossom_num(0),
1641 _blossom_index(0), _blossom_set(0), _blossom_data(0),
1642 _node_index(0), _node_heap_index(0), _node_data(0),
1643 _tree_set_index(0), _tree_set(0),
1645 _delta1_index(0), _delta1(0),
1646 _delta2_index(0), _delta2(0),
1647 _delta3_index(0), _delta3(0),
1648 _delta4_index(0), _delta4(0),
1652 ~MaxWeightedMatching() {
1653 destroyStructures();
1656 /// \name Execution Control
1657 /// The simplest way to execute the algorithm is to use the
1658 /// \ref run() member function.
1662 /// \brief Initialize the algorithm
1664 /// This function initializes the algorithm.
1668 for (ArcIt e(_graph); e != INVALID; ++e) {
1669 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1671 for (NodeIt n(_graph); n != INVALID; ++n) {
1672 (*_delta1_index)[n] = _delta1->PRE_HEAP;
1674 for (EdgeIt e(_graph); e != INVALID; ++e) {
1675 (*_delta3_index)[e] = _delta3->PRE_HEAP;
1677 for (int i = 0; i < _blossom_num; ++i) {
1678 (*_delta2_index)[i] = _delta2->PRE_HEAP;
1679 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1683 for (NodeIt n(_graph); n != INVALID; ++n) {
1685 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1686 if (_graph.target(e) == n) continue;
1687 if ((dualScale * _weight[e]) / 2 > max) {
1688 max = (dualScale * _weight[e]) / 2;
1691 (*_node_index)[n] = index;
1692 (*_node_data)[index].pot = max;
1693 _delta1->push(n, max);
1695 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1697 _tree_set->insert(blossom);
1699 (*_blossom_data)[blossom].status = EVEN;
1700 (*_blossom_data)[blossom].pred = INVALID;
1701 (*_blossom_data)[blossom].next = INVALID;
1702 (*_blossom_data)[blossom].pot = 0;
1703 (*_blossom_data)[blossom].offset = 0;
1706 for (EdgeIt e(_graph); e != INVALID; ++e) {
1707 int si = (*_node_index)[_graph.u(e)];
1708 int ti = (*_node_index)[_graph.v(e)];
1709 if (_graph.u(e) != _graph.v(e)) {
1710 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1711 dualScale * _weight[e]) / 2);
1716 /// \brief Start the algorithm
1718 /// This function starts the algorithm.
1720 /// \pre \ref init() must be called before using this function.
1726 int unmatched = _node_num;
1727 while (unmatched > 0) {
1728 Value d1 = !_delta1->empty() ?
1729 _delta1->prio() : std::numeric_limits<Value>::max();
1731 Value d2 = !_delta2->empty() ?
1732 _delta2->prio() : std::numeric_limits<Value>::max();
1734 Value d3 = !_delta3->empty() ?
1735 _delta3->prio() : std::numeric_limits<Value>::max();
1737 Value d4 = !_delta4->empty() ?
1738 _delta4->prio() : std::numeric_limits<Value>::max();
1740 _delta_sum = d1; OpType ot = D1;
1741 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1742 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1743 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1749 Node n = _delta1->top();
1756 int blossom = _delta2->top();
1757 Node n = _blossom_set->classTop(blossom);
1758 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1764 Edge e = _delta3->top();
1766 int left_blossom = _blossom_set->find(_graph.u(e));
1767 int right_blossom = _blossom_set->find(_graph.v(e));
1769 if (left_blossom == right_blossom) {
1773 if ((*_blossom_data)[left_blossom].status == EVEN) {
1774 left_tree = _tree_set->find(left_blossom);
1780 if ((*_blossom_data)[right_blossom].status == EVEN) {
1781 right_tree = _tree_set->find(right_blossom);
1787 if (left_tree == right_tree) {
1788 shrinkOnEdge(e, left_tree);
1796 splitBlossom(_delta4->top());
1803 /// \brief Run the algorithm.
1805 /// This method runs the \c %MaxWeightedMatching algorithm.
1807 /// \note mwm.run() is just a shortcut of the following code.
1819 /// \name Primal Solution
1820 /// Functions to get the primal solution, i.e. the maximum weighted
1822 /// Either \ref run() or \ref start() function should be called before
1827 /// \brief Return the weight of the matching.
1829 /// This function returns the weight of the found matching.
1831 /// \pre Either run() or start() must be called before using this function.
1832 Value matchingValue() const {
1834 for (NodeIt n(_graph); n != INVALID; ++n) {
1835 if ((*_matching)[n] != INVALID) {
1836 sum += _weight[(*_matching)[n]];
1842 /// \brief Return the size (cardinality) of the matching.
1844 /// This function returns the size (cardinality) of the found matching.
1846 /// \pre Either run() or start() must be called before using this function.
1847 int matchingSize() const {
1849 for (NodeIt n(_graph); n != INVALID; ++n) {
1850 if ((*_matching)[n] != INVALID) {
1857 /// \brief Return \c true if the given edge is in the matching.
1859 /// This function returns \c true if the given edge is in the found
1862 /// \pre Either run() or start() must be called before using this function.
1863 bool matching(const Edge& edge) const {
1864 return edge == (*_matching)[_graph.u(edge)];
1867 /// \brief Return the matching arc (or edge) incident to the given node.
1869 /// This function returns the matching arc (or edge) incident to the
1870 /// given node in the found matching or \c INVALID if the node is
1871 /// not covered by the matching.
1873 /// \pre Either run() or start() must be called before using this function.
1874 Arc matching(const Node& node) const {
1875 return (*_matching)[node];
1878 /// \brief Return the mate of the given node.
1880 /// This function returns the mate of the given node in the found
1881 /// matching or \c INVALID if the node is not covered by the matching.
1883 /// \pre Either run() or start() must be called before using this function.
1884 Node mate(const Node& node) const {
1885 return (*_matching)[node] != INVALID ?
1886 _graph.target((*_matching)[node]) : INVALID;
1891 /// \name Dual Solution
1892 /// Functions to get the dual solution.\n
1893 /// Either \ref run() or \ref start() function should be called before
1898 /// \brief Return the value of the dual solution.
1900 /// This function returns the value of the dual solution.
1901 /// It should be equal to the primal value scaled by \ref dualScale
1904 /// \pre Either run() or start() must be called before using this function.
1905 Value dualValue() const {
1907 for (NodeIt n(_graph); n != INVALID; ++n) {
1908 sum += nodeValue(n);
1910 for (int i = 0; i < blossomNum(); ++i) {
1911 sum += blossomValue(i) * (blossomSize(i) / 2);
1916 /// \brief Return the dual value (potential) of the given node.
1918 /// This function returns the dual value (potential) of the given node.
1920 /// \pre Either run() or start() must be called before using this function.
1921 Value nodeValue(const Node& n) const {
1922 return (*_node_potential)[n];
1925 /// \brief Return the number of the blossoms in the basis.
1927 /// This function returns the number of the blossoms in the basis.
1929 /// \pre Either run() or start() must be called before using this function.
1931 int blossomNum() const {
1932 return _blossom_potential.size();
1935 /// \brief Return the number of the nodes in the given blossom.
1937 /// This function returns the number of the nodes in the given blossom.
1939 /// \pre Either run() or start() must be called before using this function.
1941 int blossomSize(int k) const {
1942 return _blossom_potential[k].end - _blossom_potential[k].begin;
1945 /// \brief Return the dual value (ptential) of the given blossom.
1947 /// This function returns the dual value (ptential) of the given blossom.
1949 /// \pre Either run() or start() must be called before using this function.
1950 Value blossomValue(int k) const {
1951 return _blossom_potential[k].value;
1954 /// \brief Iterator for obtaining the nodes of a blossom.
1956 /// This class provides an iterator for obtaining the nodes of the
1957 /// given blossom. It lists a subset of the nodes.
1958 /// Before using this iterator, you must allocate a
1959 /// MaxWeightedMatching class and execute it.
1963 /// \brief Constructor.
1965 /// Constructor to get the nodes of the given variable.
1967 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1968 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1969 /// called before initializing this iterator.
1970 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1971 : _algorithm(&algorithm)
1973 _index = _algorithm->_blossom_potential[variable].begin;
1974 _last = _algorithm->_blossom_potential[variable].end;
1977 /// \brief Conversion to \c Node.
1979 /// Conversion to \c Node.
1980 operator Node() const {
1981 return _algorithm->_blossom_node_list[_index];
1984 /// \brief Increment operator.
1986 /// Increment operator.
1987 BlossomIt& operator++() {
1992 /// \brief Validity checking
1994 /// Checks whether the iterator is invalid.
1995 bool operator==(Invalid) const { return _index == _last; }
1997 /// \brief Validity checking
1999 /// Checks whether the iterator is valid.
2000 bool operator!=(Invalid) const { return _index != _last; }
2003 const MaxWeightedMatching* _algorithm;
2012 /// \ingroup matching
2014 /// \brief Weighted perfect matching in general graphs
2016 /// This class provides an efficient implementation of Edmond's
2017 /// maximum weighted perfect matching algorithm. The implementation
2018 /// is based on extensive use of priority queues and provides
2019 /// \f$O(nm\log n)\f$ time complexity.
2021 /// The maximum weighted perfect matching problem is to find a subset of
2022 /// the edges in an undirected graph with maximum overall weight for which
2023 /// each node has exactly one incident edge.
2024 /// It can be formulated with the following linear program.
2025 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2026 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2027 \quad \forall B\in\mathcal{O}\f] */
2028 /// \f[x_e \ge 0\quad \forall e\in E\f]
2029 /// \f[\max \sum_{e\in E}x_ew_e\f]
2030 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2031 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2032 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2033 /// subsets of the nodes.
2035 /// The algorithm calculates an optimal matching and a proof of the
2036 /// optimality. The solution of the dual problem can be used to check
2037 /// the result of the algorithm. The dual linear problem is the
2039 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2040 w_{uv} \quad \forall uv\in E\f] */
2041 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2042 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2043 \frac{\vert B \vert - 1}{2}z_B\f] */
2045 /// The algorithm can be executed with the run() function.
2046 /// After it the matching (the primal solution) and the dual solution
2047 /// can be obtained using the query functions and the
2048 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2049 /// which is able to iterate on the nodes of a blossom.
2050 /// If the value type is integer, then the dual solution is multiplied
2051 /// by \ref MaxWeightedMatching::dualScale "4".
2053 /// \tparam GR The graph type the algorithm runs on.
2054 /// \tparam WM The type edge weight map. The default type is
2055 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2057 template <typename GR, typename WM>
2059 template <typename GR,
2060 typename WM = typename GR::template EdgeMap<int> >
2062 class MaxWeightedPerfectMatching {
2065 /// The graph type of the algorithm
2067 /// The type of the edge weight map
2068 typedef WM WeightMap;
2069 /// The value type of the edge weights
2070 typedef typename WeightMap::Value Value;
2072 /// \brief Scaling factor for dual solution
2074 /// Scaling factor for dual solution, it is equal to 4 or 1
2075 /// according to the value type.
2076 static const int dualScale =
2077 std::numeric_limits<Value>::is_integer ? 4 : 1;
2079 typedef typename Graph::template NodeMap<typename Graph::Arc>
2084 TEMPLATE_GRAPH_TYPEDEFS(Graph);
2086 typedef typename Graph::template NodeMap<Value> NodePotential;
2087 typedef std::vector<Node> BlossomNodeList;
2089 struct BlossomVariable {
2093 BlossomVariable(int _begin, int _end, Value _value)
2094 : begin(_begin), end(_end), value(_value) {}
2098 typedef std::vector<BlossomVariable> BlossomPotential;
2100 const Graph& _graph;
2101 const WeightMap& _weight;
2103 MatchingMap* _matching;
2105 NodePotential* _node_potential;
2107 BlossomPotential _blossom_potential;
2108 BlossomNodeList _blossom_node_list;
2113 typedef RangeMap<int> IntIntMap;
2116 EVEN = -1, MATCHED = 0, ODD = 1
2119 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2120 struct BlossomData {
2127 IntNodeMap *_blossom_index;
2128 BlossomSet *_blossom_set;
2129 RangeMap<BlossomData>* _blossom_data;
2131 IntNodeMap *_node_index;
2132 IntArcMap *_node_heap_index;
2136 NodeData(IntArcMap& node_heap_index)
2137 : heap(node_heap_index) {}
2141 BinHeap<Value, IntArcMap> heap;
2142 std::map<int, Arc> heap_index;
2147 RangeMap<NodeData>* _node_data;
2149 typedef ExtendFindEnum<IntIntMap> TreeSet;
2151 IntIntMap *_tree_set_index;
2154 IntIntMap *_delta2_index;
2155 BinHeap<Value, IntIntMap> *_delta2;
2157 IntEdgeMap *_delta3_index;
2158 BinHeap<Value, IntEdgeMap> *_delta3;
2160 IntIntMap *_delta4_index;
2161 BinHeap<Value, IntIntMap> *_delta4;
2165 void createStructures() {
2166 _node_num = countNodes(_graph);
2167 _blossom_num = _node_num * 3 / 2;
2170 _matching = new MatchingMap(_graph);
2172 if (!_node_potential) {
2173 _node_potential = new NodePotential(_graph);
2175 if (!_blossom_set) {
2176 _blossom_index = new IntNodeMap(_graph);
2177 _blossom_set = new BlossomSet(*_blossom_index);
2178 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2182 _node_index = new IntNodeMap(_graph);
2183 _node_heap_index = new IntArcMap(_graph);
2184 _node_data = new RangeMap<NodeData>(_node_num,
2185 NodeData(*_node_heap_index));
2189 _tree_set_index = new IntIntMap(_blossom_num);
2190 _tree_set = new TreeSet(*_tree_set_index);
2193 _delta2_index = new IntIntMap(_blossom_num);
2194 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2197 _delta3_index = new IntEdgeMap(_graph);
2198 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2201 _delta4_index = new IntIntMap(_blossom_num);
2202 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2206 void destroyStructures() {
2207 _node_num = countNodes(_graph);
2208 _blossom_num = _node_num * 3 / 2;
2213 if (_node_potential) {
2214 delete _node_potential;
2217 delete _blossom_index;
2218 delete _blossom_set;
2219 delete _blossom_data;
2224 delete _node_heap_index;
2229 delete _tree_set_index;
2233 delete _delta2_index;
2237 delete _delta3_index;
2241 delete _delta4_index;
2246 void matchedToEven(int blossom, int tree) {
2247 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2248 _delta2->erase(blossom);
2251 if (!_blossom_set->trivial(blossom)) {
2252 (*_blossom_data)[blossom].pot -=
2253 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2256 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2257 n != INVALID; ++n) {
2259 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2260 int ni = (*_node_index)[n];
2262 (*_node_data)[ni].heap.clear();
2263 (*_node_data)[ni].heap_index.clear();
2265 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2267 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2268 Node v = _graph.source(e);
2269 int vb = _blossom_set->find(v);
2270 int vi = (*_node_index)[v];
2272 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2273 dualScale * _weight[e];
2275 if ((*_blossom_data)[vb].status == EVEN) {
2276 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2277 _delta3->push(e, rw / 2);
2280 typename std::map<int, Arc>::iterator it =
2281 (*_node_data)[vi].heap_index.find(tree);
2283 if (it != (*_node_data)[vi].heap_index.end()) {
2284 if ((*_node_data)[vi].heap[it->second] > rw) {
2285 (*_node_data)[vi].heap.replace(it->second, e);
2286 (*_node_data)[vi].heap.decrease(e, rw);
2290 (*_node_data)[vi].heap.push(e, rw);
2291 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2294 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2295 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2297 if ((*_blossom_data)[vb].status == MATCHED) {
2298 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2299 _delta2->push(vb, _blossom_set->classPrio(vb) -
2300 (*_blossom_data)[vb].offset);
2301 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2302 (*_blossom_data)[vb].offset){
2303 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2304 (*_blossom_data)[vb].offset);
2311 (*_blossom_data)[blossom].offset = 0;
2314 void matchedToOdd(int blossom) {
2315 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2316 _delta2->erase(blossom);
2318 (*_blossom_data)[blossom].offset += _delta_sum;
2319 if (!_blossom_set->trivial(blossom)) {
2320 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2321 (*_blossom_data)[blossom].offset);
2325 void evenToMatched(int blossom, int tree) {
2326 if (!_blossom_set->trivial(blossom)) {
2327 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2330 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2331 n != INVALID; ++n) {
2332 int ni = (*_node_index)[n];
2333 (*_node_data)[ni].pot -= _delta_sum;
2335 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2336 Node v = _graph.source(e);
2337 int vb = _blossom_set->find(v);
2338 int vi = (*_node_index)[v];
2340 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2341 dualScale * _weight[e];
2343 if (vb == blossom) {
2344 if (_delta3->state(e) == _delta3->IN_HEAP) {
2347 } else if ((*_blossom_data)[vb].status == EVEN) {
2349 if (_delta3->state(e) == _delta3->IN_HEAP) {
2353 int vt = _tree_set->find(vb);
2357 Arc r = _graph.oppositeArc(e);
2359 typename std::map<int, Arc>::iterator it =
2360 (*_node_data)[ni].heap_index.find(vt);
2362 if (it != (*_node_data)[ni].heap_index.end()) {
2363 if ((*_node_data)[ni].heap[it->second] > rw) {
2364 (*_node_data)[ni].heap.replace(it->second, r);
2365 (*_node_data)[ni].heap.decrease(r, rw);
2369 (*_node_data)[ni].heap.push(r, rw);
2370 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2373 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2374 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2376 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2377 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2378 (*_blossom_data)[blossom].offset);
2379 } else if ((*_delta2)[blossom] >
2380 _blossom_set->classPrio(blossom) -
2381 (*_blossom_data)[blossom].offset){
2382 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2383 (*_blossom_data)[blossom].offset);
2389 typename std::map<int, Arc>::iterator it =
2390 (*_node_data)[vi].heap_index.find(tree);
2392 if (it != (*_node_data)[vi].heap_index.end()) {
2393 (*_node_data)[vi].heap.erase(it->second);
2394 (*_node_data)[vi].heap_index.erase(it);
2395 if ((*_node_data)[vi].heap.empty()) {
2396 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2397 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2398 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2401 if ((*_blossom_data)[vb].status == MATCHED) {
2402 if (_blossom_set->classPrio(vb) ==
2403 std::numeric_limits<Value>::max()) {
2405 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2406 (*_blossom_data)[vb].offset) {
2407 _delta2->increase(vb, _blossom_set->classPrio(vb) -
2408 (*_blossom_data)[vb].offset);
2417 void oddToMatched(int blossom) {
2418 (*_blossom_data)[blossom].offset -= _delta_sum;
2420 if (_blossom_set->classPrio(blossom) !=
2421 std::numeric_limits<Value>::max()) {
2422 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2423 (*_blossom_data)[blossom].offset);
2426 if (!_blossom_set->trivial(blossom)) {
2427 _delta4->erase(blossom);
2431 void oddToEven(int blossom, int tree) {
2432 if (!_blossom_set->trivial(blossom)) {
2433 _delta4->erase(blossom);
2434 (*_blossom_data)[blossom].pot -=
2435 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2438 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2439 n != INVALID; ++n) {
2440 int ni = (*_node_index)[n];
2442 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2444 (*_node_data)[ni].heap.clear();
2445 (*_node_data)[ni].heap_index.clear();
2446 (*_node_data)[ni].pot +=
2447 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2449 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2450 Node v = _graph.source(e);
2451 int vb = _blossom_set->find(v);
2452 int vi = (*_node_index)[v];
2454 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2455 dualScale * _weight[e];
2457 if ((*_blossom_data)[vb].status == EVEN) {
2458 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2459 _delta3->push(e, rw / 2);
2463 typename std::map<int, Arc>::iterator it =
2464 (*_node_data)[vi].heap_index.find(tree);
2466 if (it != (*_node_data)[vi].heap_index.end()) {
2467 if ((*_node_data)[vi].heap[it->second] > rw) {
2468 (*_node_data)[vi].heap.replace(it->second, e);
2469 (*_node_data)[vi].heap.decrease(e, rw);
2473 (*_node_data)[vi].heap.push(e, rw);
2474 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2477 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2478 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2480 if ((*_blossom_data)[vb].status == MATCHED) {
2481 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2482 _delta2->push(vb, _blossom_set->classPrio(vb) -
2483 (*_blossom_data)[vb].offset);
2484 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2485 (*_blossom_data)[vb].offset) {
2486 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2487 (*_blossom_data)[vb].offset);
2494 (*_blossom_data)[blossom].offset = 0;
2497 void alternatePath(int even, int tree) {
2500 evenToMatched(even, tree);
2501 (*_blossom_data)[even].status = MATCHED;
2503 while ((*_blossom_data)[even].pred != INVALID) {
2504 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2505 (*_blossom_data)[odd].status = MATCHED;
2507 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2509 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2510 (*_blossom_data)[even].status = MATCHED;
2511 evenToMatched(even, tree);
2512 (*_blossom_data)[even].next =
2513 _graph.oppositeArc((*_blossom_data)[odd].pred);
2518 void destroyTree(int tree) {
2519 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2520 if ((*_blossom_data)[b].status == EVEN) {
2521 (*_blossom_data)[b].status = MATCHED;
2522 evenToMatched(b, tree);
2523 } else if ((*_blossom_data)[b].status == ODD) {
2524 (*_blossom_data)[b].status = MATCHED;
2528 _tree_set->eraseClass(tree);
2531 void augmentOnEdge(const Edge& edge) {
2533 int left = _blossom_set->find(_graph.u(edge));
2534 int right = _blossom_set->find(_graph.v(edge));
2536 int left_tree = _tree_set->find(left);
2537 alternatePath(left, left_tree);
2538 destroyTree(left_tree);
2540 int right_tree = _tree_set->find(right);
2541 alternatePath(right, right_tree);
2542 destroyTree(right_tree);
2544 (*_blossom_data)[left].next = _graph.direct(edge, true);
2545 (*_blossom_data)[right].next = _graph.direct(edge, false);
2548 void extendOnArc(const Arc& arc) {
2549 int base = _blossom_set->find(_graph.target(arc));
2550 int tree = _tree_set->find(base);
2552 int odd = _blossom_set->find(_graph.source(arc));
2553 _tree_set->insert(odd, tree);
2554 (*_blossom_data)[odd].status = ODD;
2556 (*_blossom_data)[odd].pred = arc;
2558 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2559 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2560 _tree_set->insert(even, tree);
2561 (*_blossom_data)[even].status = EVEN;
2562 matchedToEven(even, tree);
2565 void shrinkOnEdge(const Edge& edge, int tree) {
2567 std::vector<int> left_path, right_path;
2570 std::set<int> left_set, right_set;
2571 int left = _blossom_set->find(_graph.u(edge));
2572 left_path.push_back(left);
2573 left_set.insert(left);
2575 int right = _blossom_set->find(_graph.v(edge));
2576 right_path.push_back(right);
2577 right_set.insert(right);
2581 if ((*_blossom_data)[left].pred == INVALID) break;
2584 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2585 left_path.push_back(left);
2587 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2588 left_path.push_back(left);
2590 left_set.insert(left);
2592 if (right_set.find(left) != right_set.end()) {
2597 if ((*_blossom_data)[right].pred == INVALID) break;
2600 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2601 right_path.push_back(right);
2603 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2604 right_path.push_back(right);
2606 right_set.insert(right);
2608 if (left_set.find(right) != left_set.end()) {
2616 if ((*_blossom_data)[left].pred == INVALID) {
2618 while (left_set.find(nca) == left_set.end()) {
2620 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2621 right_path.push_back(nca);
2623 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2624 right_path.push_back(nca);
2628 while (right_set.find(nca) == right_set.end()) {
2630 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2631 left_path.push_back(nca);
2633 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2634 left_path.push_back(nca);
2640 std::vector<int> subblossoms;
2643 prev = _graph.direct(edge, true);
2644 for (int i = 0; left_path[i] != nca; i += 2) {
2645 subblossoms.push_back(left_path[i]);
2646 (*_blossom_data)[left_path[i]].next = prev;
2647 _tree_set->erase(left_path[i]);
2649 subblossoms.push_back(left_path[i + 1]);
2650 (*_blossom_data)[left_path[i + 1]].status = EVEN;
2651 oddToEven(left_path[i + 1], tree);
2652 _tree_set->erase(left_path[i + 1]);
2653 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2657 while (right_path[k] != nca) ++k;
2659 subblossoms.push_back(nca);
2660 (*_blossom_data)[nca].next = prev;
2662 for (int i = k - 2; i >= 0; i -= 2) {
2663 subblossoms.push_back(right_path[i + 1]);
2664 (*_blossom_data)[right_path[i + 1]].status = EVEN;
2665 oddToEven(right_path[i + 1], tree);
2666 _tree_set->erase(right_path[i + 1]);
2668 (*_blossom_data)[right_path[i + 1]].next =
2669 (*_blossom_data)[right_path[i + 1]].pred;
2671 subblossoms.push_back(right_path[i]);
2672 _tree_set->erase(right_path[i]);
2676 _blossom_set->join(subblossoms.begin(), subblossoms.end());
2678 for (int i = 0; i < int(subblossoms.size()); ++i) {
2679 if (!_blossom_set->trivial(subblossoms[i])) {
2680 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2682 (*_blossom_data)[subblossoms[i]].status = MATCHED;
2685 (*_blossom_data)[surface].pot = -2 * _delta_sum;
2686 (*_blossom_data)[surface].offset = 0;
2687 (*_blossom_data)[surface].status = EVEN;
2688 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2689 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2691 _tree_set->insert(surface, tree);
2692 _tree_set->erase(nca);
2695 void splitBlossom(int blossom) {
2696 Arc next = (*_blossom_data)[blossom].next;
2697 Arc pred = (*_blossom_data)[blossom].pred;
2699 int tree = _tree_set->find(blossom);
2701 (*_blossom_data)[blossom].status = MATCHED;
2702 oddToMatched(blossom);
2703 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2704 _delta2->erase(blossom);
2707 std::vector<int> subblossoms;
2708 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2710 Value offset = (*_blossom_data)[blossom].offset;
2711 int b = _blossom_set->find(_graph.source(pred));
2712 int d = _blossom_set->find(_graph.source(next));
2714 int ib = -1, id = -1;
2715 for (int i = 0; i < int(subblossoms.size()); ++i) {
2716 if (subblossoms[i] == b) ib = i;
2717 if (subblossoms[i] == d) id = i;
2719 (*_blossom_data)[subblossoms[i]].offset = offset;
2720 if (!_blossom_set->trivial(subblossoms[i])) {
2721 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2723 if (_blossom_set->classPrio(subblossoms[i]) !=
2724 std::numeric_limits<Value>::max()) {
2725 _delta2->push(subblossoms[i],
2726 _blossom_set->classPrio(subblossoms[i]) -
2727 (*_blossom_data)[subblossoms[i]].offset);
2731 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2732 for (int i = (id + 1) % subblossoms.size();
2733 i != ib; i = (i + 2) % subblossoms.size()) {
2734 int sb = subblossoms[i];
2735 int tb = subblossoms[(i + 1) % subblossoms.size()];
2736 (*_blossom_data)[sb].next =
2737 _graph.oppositeArc((*_blossom_data)[tb].next);
2740 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2741 int sb = subblossoms[i];
2742 int tb = subblossoms[(i + 1) % subblossoms.size()];
2743 int ub = subblossoms[(i + 2) % subblossoms.size()];
2745 (*_blossom_data)[sb].status = ODD;
2747 _tree_set->insert(sb, tree);
2748 (*_blossom_data)[sb].pred = pred;
2749 (*_blossom_data)[sb].next =
2750 _graph.oppositeArc((*_blossom_data)[tb].next);
2752 pred = (*_blossom_data)[ub].next;
2754 (*_blossom_data)[tb].status = EVEN;
2755 matchedToEven(tb, tree);
2756 _tree_set->insert(tb, tree);
2757 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2760 (*_blossom_data)[subblossoms[id]].status = ODD;
2761 matchedToOdd(subblossoms[id]);
2762 _tree_set->insert(subblossoms[id], tree);
2763 (*_blossom_data)[subblossoms[id]].next = next;
2764 (*_blossom_data)[subblossoms[id]].pred = pred;
2768 for (int i = (ib + 1) % subblossoms.size();
2769 i != id; i = (i + 2) % subblossoms.size()) {
2770 int sb = subblossoms[i];
2771 int tb = subblossoms[(i + 1) % subblossoms.size()];
2772 (*_blossom_data)[sb].next =
2773 _graph.oppositeArc((*_blossom_data)[tb].next);
2776 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2777 int sb = subblossoms[i];
2778 int tb = subblossoms[(i + 1) % subblossoms.size()];
2779 int ub = subblossoms[(i + 2) % subblossoms.size()];
2781 (*_blossom_data)[sb].status = ODD;
2783 _tree_set->insert(sb, tree);
2784 (*_blossom_data)[sb].next = next;
2785 (*_blossom_data)[sb].pred =
2786 _graph.oppositeArc((*_blossom_data)[tb].next);
2788 (*_blossom_data)[tb].status = EVEN;
2789 matchedToEven(tb, tree);
2790 _tree_set->insert(tb, tree);
2791 (*_blossom_data)[tb].pred =
2792 (*_blossom_data)[tb].next =
2793 _graph.oppositeArc((*_blossom_data)[ub].next);
2794 next = (*_blossom_data)[ub].next;
2797 (*_blossom_data)[subblossoms[ib]].status = ODD;
2798 matchedToOdd(subblossoms[ib]);
2799 _tree_set->insert(subblossoms[ib], tree);
2800 (*_blossom_data)[subblossoms[ib]].next = next;
2801 (*_blossom_data)[subblossoms[ib]].pred = pred;
2803 _tree_set->erase(blossom);
2806 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2807 if (_blossom_set->trivial(blossom)) {
2808 int bi = (*_node_index)[base];
2809 Value pot = (*_node_data)[bi].pot;
2811 (*_matching)[base] = matching;
2812 _blossom_node_list.push_back(base);
2813 (*_node_potential)[base] = pot;
2816 Value pot = (*_blossom_data)[blossom].pot;
2817 int bn = _blossom_node_list.size();
2819 std::vector<int> subblossoms;
2820 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2821 int b = _blossom_set->find(base);
2823 for (int i = 0; i < int(subblossoms.size()); ++i) {
2824 if (subblossoms[i] == b) { ib = i; break; }
2827 for (int i = 1; i < int(subblossoms.size()); i += 2) {
2828 int sb = subblossoms[(ib + i) % subblossoms.size()];
2829 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2831 Arc m = (*_blossom_data)[tb].next;
2832 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2833 extractBlossom(tb, _graph.source(m), m);
2835 extractBlossom(subblossoms[ib], base, matching);
2837 int en = _blossom_node_list.size();
2839 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2843 void extractMatching() {
2844 std::vector<int> blossoms;
2845 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2846 blossoms.push_back(c);
2849 for (int i = 0; i < int(blossoms.size()); ++i) {
2851 Value offset = (*_blossom_data)[blossoms[i]].offset;
2852 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2853 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2854 n != INVALID; ++n) {
2855 (*_node_data)[(*_node_index)[n]].pot -= offset;
2858 Arc matching = (*_blossom_data)[blossoms[i]].next;
2859 Node base = _graph.source(matching);
2860 extractBlossom(blossoms[i], base, matching);
2866 /// \brief Constructor
2869 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2870 : _graph(graph), _weight(weight), _matching(0),
2871 _node_potential(0), _blossom_potential(), _blossom_node_list(),
2872 _node_num(0), _blossom_num(0),
2874 _blossom_index(0), _blossom_set(0), _blossom_data(0),
2875 _node_index(0), _node_heap_index(0), _node_data(0),
2876 _tree_set_index(0), _tree_set(0),
2878 _delta2_index(0), _delta2(0),
2879 _delta3_index(0), _delta3(0),
2880 _delta4_index(0), _delta4(0),
2884 ~MaxWeightedPerfectMatching() {
2885 destroyStructures();
2888 /// \name Execution Control
2889 /// The simplest way to execute the algorithm is to use the
2890 /// \ref run() member function.
2894 /// \brief Initialize the algorithm
2896 /// This function initializes the algorithm.
2900 for (ArcIt e(_graph); e != INVALID; ++e) {
2901 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2903 for (EdgeIt e(_graph); e != INVALID; ++e) {
2904 (*_delta3_index)[e] = _delta3->PRE_HEAP;
2906 for (int i = 0; i < _blossom_num; ++i) {
2907 (*_delta2_index)[i] = _delta2->PRE_HEAP;
2908 (*_delta4_index)[i] = _delta4->PRE_HEAP;
2912 for (NodeIt n(_graph); n != INVALID; ++n) {
2913 Value max = - std::numeric_limits<Value>::max();
2914 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2915 if (_graph.target(e) == n) continue;
2916 if ((dualScale * _weight[e]) / 2 > max) {
2917 max = (dualScale * _weight[e]) / 2;
2920 (*_node_index)[n] = index;
2921 (*_node_data)[index].pot = max;
2923 _blossom_set->insert(n, std::numeric_limits<Value>::max());
2925 _tree_set->insert(blossom);
2927 (*_blossom_data)[blossom].status = EVEN;
2928 (*_blossom_data)[blossom].pred = INVALID;
2929 (*_blossom_data)[blossom].next = INVALID;
2930 (*_blossom_data)[blossom].pot = 0;
2931 (*_blossom_data)[blossom].offset = 0;
2934 for (EdgeIt e(_graph); e != INVALID; ++e) {
2935 int si = (*_node_index)[_graph.u(e)];
2936 int ti = (*_node_index)[_graph.v(e)];
2937 if (_graph.u(e) != _graph.v(e)) {
2938 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2939 dualScale * _weight[e]) / 2);
2944 /// \brief Start the algorithm
2946 /// This function starts the algorithm.
2948 /// \pre \ref init() must be called before using this function.
2954 int unmatched = _node_num;
2955 while (unmatched > 0) {
2956 Value d2 = !_delta2->empty() ?
2957 _delta2->prio() : std::numeric_limits<Value>::max();
2959 Value d3 = !_delta3->empty() ?
2960 _delta3->prio() : std::numeric_limits<Value>::max();
2962 Value d4 = !_delta4->empty() ?
2963 _delta4->prio() : std::numeric_limits<Value>::max();
2965 _delta_sum = d2; OpType ot = D2;
2966 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2967 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2969 if (_delta_sum == std::numeric_limits<Value>::max()) {
2976 int blossom = _delta2->top();
2977 Node n = _blossom_set->classTop(blossom);
2978 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2984 Edge e = _delta3->top();
2986 int left_blossom = _blossom_set->find(_graph.u(e));
2987 int right_blossom = _blossom_set->find(_graph.v(e));
2989 if (left_blossom == right_blossom) {
2992 int left_tree = _tree_set->find(left_blossom);
2993 int right_tree = _tree_set->find(right_blossom);
2995 if (left_tree == right_tree) {
2996 shrinkOnEdge(e, left_tree);
3004 splitBlossom(_delta4->top());
3012 /// \brief Run the algorithm.
3014 /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3016 /// \note mwpm.run() is just a shortcut of the following code.
3028 /// \name Primal Solution
3029 /// Functions to get the primal solution, i.e. the maximum weighted
3030 /// perfect matching.\n
3031 /// Either \ref run() or \ref start() function should be called before
3036 /// \brief Return the weight of the matching.
3038 /// This function returns the weight of the found matching.
3040 /// \pre Either run() or start() must be called before using this function.
3041 Value matchingValue() const {
3043 for (NodeIt n(_graph); n != INVALID; ++n) {
3044 if ((*_matching)[n] != INVALID) {
3045 sum += _weight[(*_matching)[n]];
3051 /// \brief Return \c true if the given edge is in the matching.
3053 /// This function returns \c true if the given edge is in the found
3056 /// \pre Either run() or start() must be called before using this function.
3057 bool matching(const Edge& edge) const {
3058 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3061 /// \brief Return the matching arc (or edge) incident to the given node.
3063 /// This function returns the matching arc (or edge) incident to the
3064 /// given node in the found matching or \c INVALID if the node is
3065 /// not covered by the matching.
3067 /// \pre Either run() or start() must be called before using this function.
3068 Arc matching(const Node& node) const {
3069 return (*_matching)[node];
3072 /// \brief Return the mate of the given node.
3074 /// This function returns the mate of the given node in the found
3075 /// matching or \c INVALID if the node is not covered by the matching.
3077 /// \pre Either run() or start() must be called before using this function.
3078 Node mate(const Node& node) const {
3079 return _graph.target((*_matching)[node]);
3084 /// \name Dual Solution
3085 /// Functions to get the dual solution.\n
3086 /// Either \ref run() or \ref start() function should be called before
3091 /// \brief Return the value of the dual solution.
3093 /// This function returns the value of the dual solution.
3094 /// It should be equal to the primal value scaled by \ref dualScale
3097 /// \pre Either run() or start() must be called before using this function.
3098 Value dualValue() const {
3100 for (NodeIt n(_graph); n != INVALID; ++n) {
3101 sum += nodeValue(n);
3103 for (int i = 0; i < blossomNum(); ++i) {
3104 sum += blossomValue(i) * (blossomSize(i) / 2);
3109 /// \brief Return the dual value (potential) of the given node.
3111 /// This function returns the dual value (potential) of the given node.
3113 /// \pre Either run() or start() must be called before using this function.
3114 Value nodeValue(const Node& n) const {
3115 return (*_node_potential)[n];
3118 /// \brief Return the number of the blossoms in the basis.
3120 /// This function returns the number of the blossoms in the basis.
3122 /// \pre Either run() or start() must be called before using this function.
3124 int blossomNum() const {
3125 return _blossom_potential.size();
3128 /// \brief Return the number of the nodes in the given blossom.
3130 /// This function returns the number of the nodes in the given blossom.
3132 /// \pre Either run() or start() must be called before using this function.
3134 int blossomSize(int k) const {
3135 return _blossom_potential[k].end - _blossom_potential[k].begin;
3138 /// \brief Return the dual value (ptential) of the given blossom.
3140 /// This function returns the dual value (ptential) of the given blossom.
3142 /// \pre Either run() or start() must be called before using this function.
3143 Value blossomValue(int k) const {
3144 return _blossom_potential[k].value;
3147 /// \brief Iterator for obtaining the nodes of a blossom.
3149 /// This class provides an iterator for obtaining the nodes of the
3150 /// given blossom. It lists a subset of the nodes.
3151 /// Before using this iterator, you must allocate a
3152 /// MaxWeightedPerfectMatching class and execute it.
3156 /// \brief Constructor.
3158 /// Constructor to get the nodes of the given variable.
3160 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3161 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3162 /// must be called before initializing this iterator.
3163 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3164 : _algorithm(&algorithm)
3166 _index = _algorithm->_blossom_potential[variable].begin;
3167 _last = _algorithm->_blossom_potential[variable].end;
3170 /// \brief Conversion to \c Node.
3172 /// Conversion to \c Node.
3173 operator Node() const {
3174 return _algorithm->_blossom_node_list[_index];
3177 /// \brief Increment operator.
3179 /// Increment operator.
3180 BlossomIt& operator++() {
3185 /// \brief Validity checking
3187 /// This function checks whether the iterator is invalid.
3188 bool operator==(Invalid) const { return _index == _last; }
3190 /// \brief Validity checking
3192 /// This function checks whether the iterator is valid.
3193 bool operator!=(Invalid) const { return _index != _last; }
3196 const MaxWeightedPerfectMatching* _algorithm;
3205 } //END OF NAMESPACE LEMON
3207 #endif //LEMON_MAX_MATCHING_H