1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2008
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 Edmonds' alternating forest maximum matching algorithm.
42 /// This class provides Edmonds' alternating forest matching
43 /// algorithm. The starting matching (if any) can be passed to the
44 /// algorithm using some of init functions.
46 /// The dual side of a matching is a map of the nodes to
47 /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
48 /// MATCHED/C showing the Gallai-Edmonds decomposition of the
49 /// graph. The nodes in \c EVEN/D induce a graph with
50 /// factor-critical components, the nodes in \c ODD/A form the
51 /// barrier, and the nodes in \c MATCHED/C induce a graph having a
52 /// perfect matching. The number of the fractor critical components
53 /// minus the number of barrier nodes is a lower bound on the
54 /// unmatched nodes, and if the matching is optimal this bound is
55 /// tight. This decomposition can be attained by calling \c
56 /// decomposition() after running the algorithm.
58 /// \param _Graph The graph type the algorithm runs on.
59 template <typename _Graph>
64 typedef typename Graph::template NodeMap<typename Graph::Arc>
67 ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
69 ///Indicates the Gallai-Edmonds decomposition of the graph, which
70 ///shows an upper bound on the size of a maximum matching. The
71 ///nodes with Status \c EVEN/D induce a graph with factor-critical
72 ///components, the nodes in \c ODD/A form the canonical barrier,
73 ///and the nodes in \c MATCHED/C induce a graph having a perfect
76 EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
79 typedef typename Graph::template NodeMap<Status> StatusMap;
83 TEMPLATE_GRAPH_TYPEDEFS(Graph);
85 typedef UnionFindEnum<IntNodeMap> BlossomSet;
86 typedef ExtendFindEnum<IntNodeMap> TreeSet;
87 typedef RangeMap<Node> NodeIntMap;
88 typedef MatchingMap EarMap;
89 typedef std::vector<Node> NodeQueue;
92 MatchingMap* _matching;
97 IntNodeMap* _blossom_set_index;
98 BlossomSet* _blossom_set;
99 NodeIntMap* _blossom_rep;
101 IntNodeMap* _tree_set_index;
104 NodeQueue _node_queue;
105 int _process, _postpone, _last;
111 void createStructures() {
112 _node_num = countNodes(_graph);
114 _matching = new MatchingMap(_graph);
117 _status = new StatusMap(_graph);
120 _ear = new EarMap(_graph);
123 _blossom_set_index = new IntNodeMap(_graph);
124 _blossom_set = new BlossomSet(*_blossom_set_index);
127 _blossom_rep = new NodeIntMap(_node_num);
130 _tree_set_index = new IntNodeMap(_graph);
131 _tree_set = new TreeSet(*_tree_set_index);
133 _node_queue.resize(_node_num);
136 void destroyStructures() {
148 delete _blossom_set_index;
154 delete _tree_set_index;
159 void processDense(const Node& n) {
160 _process = _postpone = _last = 0;
161 _node_queue[_last++] = n;
163 while (_process != _last) {
164 Node u = _node_queue[_process++];
165 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
166 Node v = _graph.target(a);
167 if ((*_status)[v] == MATCHED) {
169 } else if ((*_status)[v] == UNMATCHED) {
176 while (_postpone != _last) {
177 Node u = _node_queue[_postpone++];
179 for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
180 Node v = _graph.target(a);
182 if ((*_status)[v] == EVEN) {
183 if (_blossom_set->find(u) != _blossom_set->find(v)) {
188 while (_process != _last) {
189 Node w = _node_queue[_process++];
190 for (OutArcIt b(_graph, w); b != INVALID; ++b) {
191 Node x = _graph.target(b);
192 if ((*_status)[x] == MATCHED) {
194 } else if ((*_status)[x] == UNMATCHED) {
204 void processSparse(const Node& n) {
205 _process = _last = 0;
206 _node_queue[_last++] = n;
207 while (_process != _last) {
208 Node u = _node_queue[_process++];
209 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
210 Node v = _graph.target(a);
212 if ((*_status)[v] == EVEN) {
213 if (_blossom_set->find(u) != _blossom_set->find(v)) {
216 } else if ((*_status)[v] == MATCHED) {
218 } else if ((*_status)[v] == UNMATCHED) {
226 void shrinkOnEdge(const Edge& e) {
230 std::set<Node> left_set, right_set;
232 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
233 left_set.insert(left);
235 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
236 right_set.insert(right);
239 if ((*_matching)[left] == INVALID) break;
240 left = _graph.target((*_matching)[left]);
241 left = (*_blossom_rep)[_blossom_set->
242 find(_graph.target((*_ear)[left]))];
243 if (right_set.find(left) != right_set.end()) {
247 left_set.insert(left);
249 if ((*_matching)[right] == INVALID) break;
250 right = _graph.target((*_matching)[right]);
251 right = (*_blossom_rep)[_blossom_set->
252 find(_graph.target((*_ear)[right]))];
253 if (left_set.find(right) != left_set.end()) {
257 right_set.insert(right);
260 if (nca == INVALID) {
261 if ((*_matching)[left] == INVALID) {
263 while (left_set.find(nca) == left_set.end()) {
264 nca = _graph.target((*_matching)[nca]);
265 nca =(*_blossom_rep)[_blossom_set->
266 find(_graph.target((*_ear)[nca]))];
270 while (right_set.find(nca) == right_set.end()) {
271 nca = _graph.target((*_matching)[nca]);
272 nca = (*_blossom_rep)[_blossom_set->
273 find(_graph.target((*_ear)[nca]))];
281 Node node = _graph.u(e);
282 Arc arc = _graph.direct(e, true);
283 Node base = (*_blossom_rep)[_blossom_set->find(node)];
285 while (base != nca) {
286 _ear->set(node, arc);
290 n = _graph.target((*_matching)[n]);
292 n = _graph.target(a);
293 _ear->set(n, _graph.oppositeArc(a));
295 node = _graph.target((*_matching)[base]);
296 _tree_set->erase(base);
297 _tree_set->erase(node);
298 _blossom_set->insert(node, _blossom_set->find(base));
299 _status->set(node, EVEN);
300 _node_queue[_last++] = node;
301 arc = _graph.oppositeArc((*_ear)[node]);
302 node = _graph.target((*_ear)[node]);
303 base = (*_blossom_rep)[_blossom_set->find(node)];
304 _blossom_set->join(_graph.target(arc), base);
308 _blossom_rep->set(_blossom_set->find(nca), nca);
312 Node node = _graph.v(e);
313 Arc arc = _graph.direct(e, false);
314 Node base = (*_blossom_rep)[_blossom_set->find(node)];
316 while (base != nca) {
317 _ear->set(node, arc);
321 n = _graph.target((*_matching)[n]);
323 n = _graph.target(a);
324 _ear->set(n, _graph.oppositeArc(a));
326 node = _graph.target((*_matching)[base]);
327 _tree_set->erase(base);
328 _tree_set->erase(node);
329 _blossom_set->insert(node, _blossom_set->find(base));
330 _status->set(node, EVEN);
331 _node_queue[_last++] = node;
332 arc = _graph.oppositeArc((*_ear)[node]);
333 node = _graph.target((*_ear)[node]);
334 base = (*_blossom_rep)[_blossom_set->find(node)];
335 _blossom_set->join(_graph.target(arc), base);
339 _blossom_rep->set(_blossom_set->find(nca), nca);
344 void extendOnArc(const Arc& a) {
345 Node base = _graph.source(a);
346 Node odd = _graph.target(a);
348 _ear->set(odd, _graph.oppositeArc(a));
349 Node even = _graph.target((*_matching)[odd]);
350 _blossom_rep->set(_blossom_set->insert(even), even);
351 _status->set(odd, ODD);
352 _status->set(even, EVEN);
353 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
354 _tree_set->insert(odd, tree);
355 _tree_set->insert(even, tree);
356 _node_queue[_last++] = even;
360 void augmentOnArc(const Arc& a) {
361 Node even = _graph.source(a);
362 Node odd = _graph.target(a);
364 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
366 _matching->set(odd, _graph.oppositeArc(a));
367 _status->set(odd, MATCHED);
369 Arc arc = (*_matching)[even];
370 _matching->set(even, a);
372 while (arc != INVALID) {
373 odd = _graph.target(arc);
375 even = _graph.target(arc);
376 _matching->set(odd, arc);
377 arc = (*_matching)[even];
378 _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
381 for (typename TreeSet::ItemIt it(*_tree_set, tree);
382 it != INVALID; ++it) {
383 if ((*_status)[it] == ODD) {
384 _status->set(it, MATCHED);
386 int blossom = _blossom_set->find(it);
387 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
388 jt != INVALID; ++jt) {
389 _status->set(jt, MATCHED);
391 _blossom_set->eraseClass(blossom);
394 _tree_set->eraseClass(tree);
400 /// \brief Constructor
403 MaxMatching(const Graph& graph)
404 : _graph(graph), _matching(0), _status(0), _ear(0),
405 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
406 _tree_set_index(0), _tree_set(0) {}
412 /// \name Execution control
413 /// The simplest way to execute the algorithm is to use the member
414 /// \c run() member function.
417 /// If you need more control on the execution, first you must call
418 /// \ref init(), \ref greedyInit() or \ref matchingInit()
419 /// functions, then you can start the algorithm with the \ref
420 /// startParse() or startDense() functions.
424 /// \brief Sets the actual matching to the empty matching.
426 /// Sets the actual matching to the empty matching.
430 for(NodeIt n(_graph); n != INVALID; ++n) {
431 _matching->set(n, INVALID);
432 _status->set(n, UNMATCHED);
436 ///\brief Finds a greedy matching for initial matching.
438 ///For initial matchig it finds a maximal greedy matching.
441 for (NodeIt n(_graph); n != INVALID; ++n) {
442 _matching->set(n, INVALID);
443 _status->set(n, UNMATCHED);
445 for (NodeIt n(_graph); n != INVALID; ++n) {
446 if ((*_matching)[n] == INVALID) {
447 for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
448 Node v = _graph.target(a);
449 if ((*_matching)[v] == INVALID && v != n) {
450 _matching->set(n, a);
451 _status->set(n, MATCHED);
452 _matching->set(v, _graph.oppositeArc(a));
453 _status->set(v, MATCHED);
462 /// \brief Initialize the matching from the map containing a matching.
464 /// Initialize the matching from a \c bool valued \c Edge map. This
465 /// map must have the property that there are no two incident edges
466 /// with true value, ie. it contains a matching.
467 /// \return %True if the map contains a matching.
468 template <typename MatchingMap>
469 bool matchingInit(const MatchingMap& matching) {
472 for (NodeIt n(_graph); n != INVALID; ++n) {
473 _matching->set(n, INVALID);
474 _status->set(n, UNMATCHED);
476 for(EdgeIt e(_graph); e!=INVALID; ++e) {
479 Node u = _graph.u(e);
480 if ((*_matching)[u] != INVALID) return false;
481 _matching->set(u, _graph.direct(e, true));
482 _status->set(u, MATCHED);
484 Node v = _graph.v(e);
485 if ((*_matching)[v] != INVALID) return false;
486 _matching->set(v, _graph.direct(e, false));
487 _status->set(v, MATCHED);
493 /// \brief Starts Edmonds' algorithm
495 /// If runs the original Edmonds' algorithm.
497 for(NodeIt n(_graph); n != INVALID; ++n) {
498 if ((*_status)[n] == UNMATCHED) {
499 (*_blossom_rep)[_blossom_set->insert(n)] = n;
500 _tree_set->insert(n);
501 _status->set(n, EVEN);
507 /// \brief Starts Edmonds' algorithm.
509 /// It runs Edmonds' algorithm with a heuristic of postponing
510 /// shrinks, giving a faster algorithm for dense graphs.
512 for(NodeIt n(_graph); n != INVALID; ++n) {
513 if ((*_status)[n] == UNMATCHED) {
514 (*_blossom_rep)[_blossom_set->insert(n)] = n;
515 _tree_set->insert(n);
516 _status->set(n, EVEN);
523 /// \brief Runs Edmonds' algorithm
525 /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
526 /// or Edmonds' algorithm with a heuristic of
527 /// postponing shrinks for dense graphs.
529 if (countEdges(_graph) < 2 * countNodes(_graph)) {
540 /// \name Primal solution
541 /// Functions for get the primal solution, ie. the matching.
545 ///\brief Returns the size of the actual matching stored.
547 ///Returns the size of the actual matching stored. After \ref
548 ///run() it returns the size of the maximum matching in the graph.
549 int matchingSize() const {
551 for (NodeIt n(_graph); n != INVALID; ++n) {
552 if ((*_matching)[n] != INVALID) {
559 /// \brief Returns true when the edge is in the matching.
561 /// Returns true when the edge is in the matching.
562 bool matching(const Edge& edge) const {
563 return edge == (*_matching)[_graph.u(edge)];
566 /// \brief Returns the matching edge incident to the given node.
568 /// Returns the matching edge of a \c node in the actual matching or
569 /// INVALID if the \c node is not covered by the actual matching.
570 Arc matching(const Node& n) const {
571 return (*_matching)[n];
574 ///\brief Returns the mate of a node in the actual matching.
576 ///Returns the mate of a \c node in the actual matching or
577 ///INVALID if the \c node is not covered by the actual matching.
578 Node mate(const Node& n) const {
579 return (*_matching)[n] != INVALID ?
580 _graph.target((*_matching)[n]) : INVALID;
585 /// \name Dual solution
586 /// Functions for get the dual solution, ie. the decomposition.
590 /// \brief Returns the class of the node in the Edmonds-Gallai
593 /// Returns the class of the node in the Edmonds-Gallai
595 Status decomposition(const Node& n) const {
596 return (*_status)[n];
599 /// \brief Returns true when the node is in the barrier.
601 /// Returns true when the node is in the barrier.
602 bool barrier(const Node& n) const {
603 return (*_status)[n] == ODD;
610 /// \ingroup matching
612 /// \brief Weighted matching in general graphs
614 /// This class provides an efficient implementation of Edmond's
615 /// maximum weighted matching algorithm. The implementation is based
616 /// on extensive use of priority queues and provides
617 /// \f$O(nm\log(n))\f$ time complexity.
619 /// The maximum weighted matching problem is to find undirected
620 /// edges in the graph with maximum overall weight and no two of
621 /// them shares their ends. The problem can be formulated with the
622 /// following linear program.
623 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
624 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
625 \quad \forall B\in\mathcal{O}\f] */
626 /// \f[x_e \ge 0\quad \forall e\in E\f]
627 /// \f[\max \sum_{e\in E}x_ew_e\f]
628 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
629 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
630 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
631 /// subsets of the nodes.
633 /// The algorithm calculates an optimal matching and a proof of the
634 /// optimality. The solution of the dual problem can be used to check
635 /// the result of the algorithm. The dual linear problem is the
636 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
637 z_B \ge w_{uv} \quad \forall uv\in E\f] */
638 /// \f[y_u \ge 0 \quad \forall u \in V\f]
639 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
640 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
641 \frac{\vert B \vert - 1}{2}z_B\f] */
643 /// The algorithm can be executed with \c run() or the \c init() and
644 /// then the \c start() member functions. After it the matching can
645 /// be asked with \c matching() or mate() functions. The dual
646 /// solution can be get with \c nodeValue(), \c blossomNum() and \c
647 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
648 /// "BlossomIt" nested class which is able to iterate on the nodes
649 /// of a blossom. If the value type is integral then the dual
650 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
651 template <typename _Graph,
652 typename _WeightMap = typename _Graph::template EdgeMap<int> >
653 class MaxWeightedMatching {
656 typedef _Graph Graph;
657 typedef _WeightMap WeightMap;
658 typedef typename WeightMap::Value Value;
660 /// \brief Scaling factor for dual solution
662 /// Scaling factor for dual solution, it is equal to 4 or 1
663 /// according to the value type.
664 static const int dualScale =
665 std::numeric_limits<Value>::is_integer ? 4 : 1;
667 typedef typename Graph::template NodeMap<typename Graph::Arc>
672 TEMPLATE_GRAPH_TYPEDEFS(Graph);
674 typedef typename Graph::template NodeMap<Value> NodePotential;
675 typedef std::vector<Node> BlossomNodeList;
677 struct BlossomVariable {
681 BlossomVariable(int _begin, int _end, Value _value)
682 : begin(_begin), end(_end), value(_value) {}
686 typedef std::vector<BlossomVariable> BlossomPotential;
689 const WeightMap& _weight;
691 MatchingMap* _matching;
693 NodePotential* _node_potential;
695 BlossomPotential _blossom_potential;
696 BlossomNodeList _blossom_node_list;
701 typedef RangeMap<int> IntIntMap;
704 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
707 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
716 IntNodeMap *_blossom_index;
717 BlossomSet *_blossom_set;
718 RangeMap<BlossomData>* _blossom_data;
720 IntNodeMap *_node_index;
721 IntArcMap *_node_heap_index;
725 NodeData(IntArcMap& node_heap_index)
726 : heap(node_heap_index) {}
730 BinHeap<Value, IntArcMap> heap;
731 std::map<int, Arc> heap_index;
736 RangeMap<NodeData>* _node_data;
738 typedef ExtendFindEnum<IntIntMap> TreeSet;
740 IntIntMap *_tree_set_index;
743 IntNodeMap *_delta1_index;
744 BinHeap<Value, IntNodeMap> *_delta1;
746 IntIntMap *_delta2_index;
747 BinHeap<Value, IntIntMap> *_delta2;
749 IntEdgeMap *_delta3_index;
750 BinHeap<Value, IntEdgeMap> *_delta3;
752 IntIntMap *_delta4_index;
753 BinHeap<Value, IntIntMap> *_delta4;
757 void createStructures() {
758 _node_num = countNodes(_graph);
759 _blossom_num = _node_num * 3 / 2;
762 _matching = new MatchingMap(_graph);
764 if (!_node_potential) {
765 _node_potential = new NodePotential(_graph);
768 _blossom_index = new IntNodeMap(_graph);
769 _blossom_set = new BlossomSet(*_blossom_index);
770 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
774 _node_index = new IntNodeMap(_graph);
775 _node_heap_index = new IntArcMap(_graph);
776 _node_data = new RangeMap<NodeData>(_node_num,
777 NodeData(*_node_heap_index));
781 _tree_set_index = new IntIntMap(_blossom_num);
782 _tree_set = new TreeSet(*_tree_set_index);
785 _delta1_index = new IntNodeMap(_graph);
786 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
789 _delta2_index = new IntIntMap(_blossom_num);
790 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
793 _delta3_index = new IntEdgeMap(_graph);
794 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
797 _delta4_index = new IntIntMap(_blossom_num);
798 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
802 void destroyStructures() {
803 _node_num = countNodes(_graph);
804 _blossom_num = _node_num * 3 / 2;
809 if (_node_potential) {
810 delete _node_potential;
813 delete _blossom_index;
815 delete _blossom_data;
820 delete _node_heap_index;
825 delete _tree_set_index;
829 delete _delta1_index;
833 delete _delta2_index;
837 delete _delta3_index;
841 delete _delta4_index;
846 void matchedToEven(int blossom, int tree) {
847 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
848 _delta2->erase(blossom);
851 if (!_blossom_set->trivial(blossom)) {
852 (*_blossom_data)[blossom].pot -=
853 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
856 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
859 _blossom_set->increase(n, std::numeric_limits<Value>::max());
860 int ni = (*_node_index)[n];
862 (*_node_data)[ni].heap.clear();
863 (*_node_data)[ni].heap_index.clear();
865 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
867 _delta1->push(n, (*_node_data)[ni].pot);
869 for (InArcIt e(_graph, n); e != INVALID; ++e) {
870 Node v = _graph.source(e);
871 int vb = _blossom_set->find(v);
872 int vi = (*_node_index)[v];
874 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
875 dualScale * _weight[e];
877 if ((*_blossom_data)[vb].status == EVEN) {
878 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
879 _delta3->push(e, rw / 2);
881 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
882 if (_delta3->state(e) != _delta3->IN_HEAP) {
883 _delta3->push(e, rw);
886 typename std::map<int, Arc>::iterator it =
887 (*_node_data)[vi].heap_index.find(tree);
889 if (it != (*_node_data)[vi].heap_index.end()) {
890 if ((*_node_data)[vi].heap[it->second] > rw) {
891 (*_node_data)[vi].heap.replace(it->second, e);
892 (*_node_data)[vi].heap.decrease(e, rw);
896 (*_node_data)[vi].heap.push(e, rw);
897 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
900 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
901 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
903 if ((*_blossom_data)[vb].status == MATCHED) {
904 if (_delta2->state(vb) != _delta2->IN_HEAP) {
905 _delta2->push(vb, _blossom_set->classPrio(vb) -
906 (*_blossom_data)[vb].offset);
907 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
908 (*_blossom_data)[vb].offset){
909 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
910 (*_blossom_data)[vb].offset);
917 (*_blossom_data)[blossom].offset = 0;
920 void matchedToOdd(int blossom) {
921 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
922 _delta2->erase(blossom);
924 (*_blossom_data)[blossom].offset += _delta_sum;
925 if (!_blossom_set->trivial(blossom)) {
926 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
927 (*_blossom_data)[blossom].offset);
931 void evenToMatched(int blossom, int tree) {
932 if (!_blossom_set->trivial(blossom)) {
933 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
936 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
938 int ni = (*_node_index)[n];
939 (*_node_data)[ni].pot -= _delta_sum;
943 for (InArcIt e(_graph, n); e != INVALID; ++e) {
944 Node v = _graph.source(e);
945 int vb = _blossom_set->find(v);
946 int vi = (*_node_index)[v];
948 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
949 dualScale * _weight[e];
952 if (_delta3->state(e) == _delta3->IN_HEAP) {
955 } else if ((*_blossom_data)[vb].status == EVEN) {
957 if (_delta3->state(e) == _delta3->IN_HEAP) {
961 int vt = _tree_set->find(vb);
965 Arc r = _graph.oppositeArc(e);
967 typename std::map<int, Arc>::iterator it =
968 (*_node_data)[ni].heap_index.find(vt);
970 if (it != (*_node_data)[ni].heap_index.end()) {
971 if ((*_node_data)[ni].heap[it->second] > rw) {
972 (*_node_data)[ni].heap.replace(it->second, r);
973 (*_node_data)[ni].heap.decrease(r, rw);
977 (*_node_data)[ni].heap.push(r, rw);
978 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
981 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
982 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
984 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
985 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
986 (*_blossom_data)[blossom].offset);
987 } else if ((*_delta2)[blossom] >
988 _blossom_set->classPrio(blossom) -
989 (*_blossom_data)[blossom].offset){
990 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
991 (*_blossom_data)[blossom].offset);
996 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
997 if (_delta3->state(e) == _delta3->IN_HEAP) {
1002 typename std::map<int, Arc>::iterator it =
1003 (*_node_data)[vi].heap_index.find(tree);
1005 if (it != (*_node_data)[vi].heap_index.end()) {
1006 (*_node_data)[vi].heap.erase(it->second);
1007 (*_node_data)[vi].heap_index.erase(it);
1008 if ((*_node_data)[vi].heap.empty()) {
1009 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1010 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1011 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1014 if ((*_blossom_data)[vb].status == MATCHED) {
1015 if (_blossom_set->classPrio(vb) ==
1016 std::numeric_limits<Value>::max()) {
1018 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1019 (*_blossom_data)[vb].offset) {
1020 _delta2->increase(vb, _blossom_set->classPrio(vb) -
1021 (*_blossom_data)[vb].offset);
1030 void oddToMatched(int blossom) {
1031 (*_blossom_data)[blossom].offset -= _delta_sum;
1033 if (_blossom_set->classPrio(blossom) !=
1034 std::numeric_limits<Value>::max()) {
1035 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1036 (*_blossom_data)[blossom].offset);
1039 if (!_blossom_set->trivial(blossom)) {
1040 _delta4->erase(blossom);
1044 void oddToEven(int blossom, int tree) {
1045 if (!_blossom_set->trivial(blossom)) {
1046 _delta4->erase(blossom);
1047 (*_blossom_data)[blossom].pot -=
1048 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1051 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1052 n != INVALID; ++n) {
1053 int ni = (*_node_index)[n];
1055 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1057 (*_node_data)[ni].heap.clear();
1058 (*_node_data)[ni].heap_index.clear();
1059 (*_node_data)[ni].pot +=
1060 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1062 _delta1->push(n, (*_node_data)[ni].pot);
1064 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1065 Node v = _graph.source(e);
1066 int vb = _blossom_set->find(v);
1067 int vi = (*_node_index)[v];
1069 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1070 dualScale * _weight[e];
1072 if ((*_blossom_data)[vb].status == EVEN) {
1073 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1074 _delta3->push(e, rw / 2);
1076 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1077 if (_delta3->state(e) != _delta3->IN_HEAP) {
1078 _delta3->push(e, rw);
1082 typename std::map<int, Arc>::iterator it =
1083 (*_node_data)[vi].heap_index.find(tree);
1085 if (it != (*_node_data)[vi].heap_index.end()) {
1086 if ((*_node_data)[vi].heap[it->second] > rw) {
1087 (*_node_data)[vi].heap.replace(it->second, e);
1088 (*_node_data)[vi].heap.decrease(e, rw);
1092 (*_node_data)[vi].heap.push(e, rw);
1093 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1096 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1097 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1099 if ((*_blossom_data)[vb].status == MATCHED) {
1100 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1101 _delta2->push(vb, _blossom_set->classPrio(vb) -
1102 (*_blossom_data)[vb].offset);
1103 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1104 (*_blossom_data)[vb].offset) {
1105 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1106 (*_blossom_data)[vb].offset);
1113 (*_blossom_data)[blossom].offset = 0;
1117 void matchedToUnmatched(int blossom) {
1118 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1119 _delta2->erase(blossom);
1122 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1123 n != INVALID; ++n) {
1124 int ni = (*_node_index)[n];
1126 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1128 (*_node_data)[ni].heap.clear();
1129 (*_node_data)[ni].heap_index.clear();
1131 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1132 Node v = _graph.target(e);
1133 int vb = _blossom_set->find(v);
1134 int vi = (*_node_index)[v];
1136 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1137 dualScale * _weight[e];
1139 if ((*_blossom_data)[vb].status == EVEN) {
1140 if (_delta3->state(e) != _delta3->IN_HEAP) {
1141 _delta3->push(e, rw);
1148 void unmatchedToMatched(int blossom) {
1149 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1150 n != INVALID; ++n) {
1151 int ni = (*_node_index)[n];
1153 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1154 Node v = _graph.source(e);
1155 int vb = _blossom_set->find(v);
1156 int vi = (*_node_index)[v];
1158 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1159 dualScale * _weight[e];
1161 if (vb == blossom) {
1162 if (_delta3->state(e) == _delta3->IN_HEAP) {
1165 } else if ((*_blossom_data)[vb].status == EVEN) {
1167 if (_delta3->state(e) == _delta3->IN_HEAP) {
1171 int vt = _tree_set->find(vb);
1173 Arc r = _graph.oppositeArc(e);
1175 typename std::map<int, Arc>::iterator it =
1176 (*_node_data)[ni].heap_index.find(vt);
1178 if (it != (*_node_data)[ni].heap_index.end()) {
1179 if ((*_node_data)[ni].heap[it->second] > rw) {
1180 (*_node_data)[ni].heap.replace(it->second, r);
1181 (*_node_data)[ni].heap.decrease(r, rw);
1185 (*_node_data)[ni].heap.push(r, rw);
1186 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1189 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1190 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1192 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1193 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1194 (*_blossom_data)[blossom].offset);
1195 } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1196 (*_blossom_data)[blossom].offset){
1197 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1198 (*_blossom_data)[blossom].offset);
1202 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1203 if (_delta3->state(e) == _delta3->IN_HEAP) {
1211 void alternatePath(int even, int tree) {
1214 evenToMatched(even, tree);
1215 (*_blossom_data)[even].status = MATCHED;
1217 while ((*_blossom_data)[even].pred != INVALID) {
1218 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1219 (*_blossom_data)[odd].status = MATCHED;
1221 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1223 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1224 (*_blossom_data)[even].status = MATCHED;
1225 evenToMatched(even, tree);
1226 (*_blossom_data)[even].next =
1227 _graph.oppositeArc((*_blossom_data)[odd].pred);
1232 void destroyTree(int tree) {
1233 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1234 if ((*_blossom_data)[b].status == EVEN) {
1235 (*_blossom_data)[b].status = MATCHED;
1236 evenToMatched(b, tree);
1237 } else if ((*_blossom_data)[b].status == ODD) {
1238 (*_blossom_data)[b].status = MATCHED;
1242 _tree_set->eraseClass(tree);
1246 void unmatchNode(const Node& node) {
1247 int blossom = _blossom_set->find(node);
1248 int tree = _tree_set->find(blossom);
1250 alternatePath(blossom, tree);
1253 (*_blossom_data)[blossom].status = UNMATCHED;
1254 (*_blossom_data)[blossom].base = node;
1255 matchedToUnmatched(blossom);
1259 void augmentOnEdge(const Edge& edge) {
1261 int left = _blossom_set->find(_graph.u(edge));
1262 int right = _blossom_set->find(_graph.v(edge));
1264 if ((*_blossom_data)[left].status == EVEN) {
1265 int left_tree = _tree_set->find(left);
1266 alternatePath(left, left_tree);
1267 destroyTree(left_tree);
1269 (*_blossom_data)[left].status = MATCHED;
1270 unmatchedToMatched(left);
1273 if ((*_blossom_data)[right].status == EVEN) {
1274 int right_tree = _tree_set->find(right);
1275 alternatePath(right, right_tree);
1276 destroyTree(right_tree);
1278 (*_blossom_data)[right].status = MATCHED;
1279 unmatchedToMatched(right);
1282 (*_blossom_data)[left].next = _graph.direct(edge, true);
1283 (*_blossom_data)[right].next = _graph.direct(edge, false);
1286 void extendOnArc(const Arc& arc) {
1287 int base = _blossom_set->find(_graph.target(arc));
1288 int tree = _tree_set->find(base);
1290 int odd = _blossom_set->find(_graph.source(arc));
1291 _tree_set->insert(odd, tree);
1292 (*_blossom_data)[odd].status = ODD;
1294 (*_blossom_data)[odd].pred = arc;
1296 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1297 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1298 _tree_set->insert(even, tree);
1299 (*_blossom_data)[even].status = EVEN;
1300 matchedToEven(even, tree);
1303 void shrinkOnEdge(const Edge& edge, int tree) {
1305 std::vector<int> left_path, right_path;
1308 std::set<int> left_set, right_set;
1309 int left = _blossom_set->find(_graph.u(edge));
1310 left_path.push_back(left);
1311 left_set.insert(left);
1313 int right = _blossom_set->find(_graph.v(edge));
1314 right_path.push_back(right);
1315 right_set.insert(right);
1319 if ((*_blossom_data)[left].pred == INVALID) break;
1322 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1323 left_path.push_back(left);
1325 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1326 left_path.push_back(left);
1328 left_set.insert(left);
1330 if (right_set.find(left) != right_set.end()) {
1335 if ((*_blossom_data)[right].pred == INVALID) break;
1338 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1339 right_path.push_back(right);
1341 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1342 right_path.push_back(right);
1344 right_set.insert(right);
1346 if (left_set.find(right) != left_set.end()) {
1354 if ((*_blossom_data)[left].pred == INVALID) {
1356 while (left_set.find(nca) == left_set.end()) {
1358 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1359 right_path.push_back(nca);
1361 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1362 right_path.push_back(nca);
1366 while (right_set.find(nca) == right_set.end()) {
1368 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1369 left_path.push_back(nca);
1371 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1372 left_path.push_back(nca);
1378 std::vector<int> subblossoms;
1381 prev = _graph.direct(edge, true);
1382 for (int i = 0; left_path[i] != nca; i += 2) {
1383 subblossoms.push_back(left_path[i]);
1384 (*_blossom_data)[left_path[i]].next = prev;
1385 _tree_set->erase(left_path[i]);
1387 subblossoms.push_back(left_path[i + 1]);
1388 (*_blossom_data)[left_path[i + 1]].status = EVEN;
1389 oddToEven(left_path[i + 1], tree);
1390 _tree_set->erase(left_path[i + 1]);
1391 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1395 while (right_path[k] != nca) ++k;
1397 subblossoms.push_back(nca);
1398 (*_blossom_data)[nca].next = prev;
1400 for (int i = k - 2; i >= 0; i -= 2) {
1401 subblossoms.push_back(right_path[i + 1]);
1402 (*_blossom_data)[right_path[i + 1]].status = EVEN;
1403 oddToEven(right_path[i + 1], tree);
1404 _tree_set->erase(right_path[i + 1]);
1406 (*_blossom_data)[right_path[i + 1]].next =
1407 (*_blossom_data)[right_path[i + 1]].pred;
1409 subblossoms.push_back(right_path[i]);
1410 _tree_set->erase(right_path[i]);
1414 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1416 for (int i = 0; i < int(subblossoms.size()); ++i) {
1417 if (!_blossom_set->trivial(subblossoms[i])) {
1418 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1420 (*_blossom_data)[subblossoms[i]].status = MATCHED;
1423 (*_blossom_data)[surface].pot = -2 * _delta_sum;
1424 (*_blossom_data)[surface].offset = 0;
1425 (*_blossom_data)[surface].status = EVEN;
1426 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1427 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1429 _tree_set->insert(surface, tree);
1430 _tree_set->erase(nca);
1433 void splitBlossom(int blossom) {
1434 Arc next = (*_blossom_data)[blossom].next;
1435 Arc pred = (*_blossom_data)[blossom].pred;
1437 int tree = _tree_set->find(blossom);
1439 (*_blossom_data)[blossom].status = MATCHED;
1440 oddToMatched(blossom);
1441 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1442 _delta2->erase(blossom);
1445 std::vector<int> subblossoms;
1446 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1448 Value offset = (*_blossom_data)[blossom].offset;
1449 int b = _blossom_set->find(_graph.source(pred));
1450 int d = _blossom_set->find(_graph.source(next));
1452 int ib = -1, id = -1;
1453 for (int i = 0; i < int(subblossoms.size()); ++i) {
1454 if (subblossoms[i] == b) ib = i;
1455 if (subblossoms[i] == d) id = i;
1457 (*_blossom_data)[subblossoms[i]].offset = offset;
1458 if (!_blossom_set->trivial(subblossoms[i])) {
1459 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1461 if (_blossom_set->classPrio(subblossoms[i]) !=
1462 std::numeric_limits<Value>::max()) {
1463 _delta2->push(subblossoms[i],
1464 _blossom_set->classPrio(subblossoms[i]) -
1465 (*_blossom_data)[subblossoms[i]].offset);
1469 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1470 for (int i = (id + 1) % subblossoms.size();
1471 i != ib; i = (i + 2) % subblossoms.size()) {
1472 int sb = subblossoms[i];
1473 int tb = subblossoms[(i + 1) % subblossoms.size()];
1474 (*_blossom_data)[sb].next =
1475 _graph.oppositeArc((*_blossom_data)[tb].next);
1478 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1479 int sb = subblossoms[i];
1480 int tb = subblossoms[(i + 1) % subblossoms.size()];
1481 int ub = subblossoms[(i + 2) % subblossoms.size()];
1483 (*_blossom_data)[sb].status = ODD;
1485 _tree_set->insert(sb, tree);
1486 (*_blossom_data)[sb].pred = pred;
1487 (*_blossom_data)[sb].next =
1488 _graph.oppositeArc((*_blossom_data)[tb].next);
1490 pred = (*_blossom_data)[ub].next;
1492 (*_blossom_data)[tb].status = EVEN;
1493 matchedToEven(tb, tree);
1494 _tree_set->insert(tb, tree);
1495 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1498 (*_blossom_data)[subblossoms[id]].status = ODD;
1499 matchedToOdd(subblossoms[id]);
1500 _tree_set->insert(subblossoms[id], tree);
1501 (*_blossom_data)[subblossoms[id]].next = next;
1502 (*_blossom_data)[subblossoms[id]].pred = pred;
1506 for (int i = (ib + 1) % subblossoms.size();
1507 i != id; i = (i + 2) % subblossoms.size()) {
1508 int sb = subblossoms[i];
1509 int tb = subblossoms[(i + 1) % subblossoms.size()];
1510 (*_blossom_data)[sb].next =
1511 _graph.oppositeArc((*_blossom_data)[tb].next);
1514 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1515 int sb = subblossoms[i];
1516 int tb = subblossoms[(i + 1) % subblossoms.size()];
1517 int ub = subblossoms[(i + 2) % subblossoms.size()];
1519 (*_blossom_data)[sb].status = ODD;
1521 _tree_set->insert(sb, tree);
1522 (*_blossom_data)[sb].next = next;
1523 (*_blossom_data)[sb].pred =
1524 _graph.oppositeArc((*_blossom_data)[tb].next);
1526 (*_blossom_data)[tb].status = EVEN;
1527 matchedToEven(tb, tree);
1528 _tree_set->insert(tb, tree);
1529 (*_blossom_data)[tb].pred =
1530 (*_blossom_data)[tb].next =
1531 _graph.oppositeArc((*_blossom_data)[ub].next);
1532 next = (*_blossom_data)[ub].next;
1535 (*_blossom_data)[subblossoms[ib]].status = ODD;
1536 matchedToOdd(subblossoms[ib]);
1537 _tree_set->insert(subblossoms[ib], tree);
1538 (*_blossom_data)[subblossoms[ib]].next = next;
1539 (*_blossom_data)[subblossoms[ib]].pred = pred;
1541 _tree_set->erase(blossom);
1544 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1545 if (_blossom_set->trivial(blossom)) {
1546 int bi = (*_node_index)[base];
1547 Value pot = (*_node_data)[bi].pot;
1549 _matching->set(base, matching);
1550 _blossom_node_list.push_back(base);
1551 _node_potential->set(base, pot);
1554 Value pot = (*_blossom_data)[blossom].pot;
1555 int bn = _blossom_node_list.size();
1557 std::vector<int> subblossoms;
1558 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1559 int b = _blossom_set->find(base);
1561 for (int i = 0; i < int(subblossoms.size()); ++i) {
1562 if (subblossoms[i] == b) { ib = i; break; }
1565 for (int i = 1; i < int(subblossoms.size()); i += 2) {
1566 int sb = subblossoms[(ib + i) % subblossoms.size()];
1567 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1569 Arc m = (*_blossom_data)[tb].next;
1570 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1571 extractBlossom(tb, _graph.source(m), m);
1573 extractBlossom(subblossoms[ib], base, matching);
1575 int en = _blossom_node_list.size();
1577 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1581 void extractMatching() {
1582 std::vector<int> blossoms;
1583 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1584 blossoms.push_back(c);
1587 for (int i = 0; i < int(blossoms.size()); ++i) {
1588 if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1590 Value offset = (*_blossom_data)[blossoms[i]].offset;
1591 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1592 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1593 n != INVALID; ++n) {
1594 (*_node_data)[(*_node_index)[n]].pot -= offset;
1597 Arc matching = (*_blossom_data)[blossoms[i]].next;
1598 Node base = _graph.source(matching);
1599 extractBlossom(blossoms[i], base, matching);
1601 Node base = (*_blossom_data)[blossoms[i]].base;
1602 extractBlossom(blossoms[i], base, INVALID);
1609 /// \brief Constructor
1612 MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1613 : _graph(graph), _weight(weight), _matching(0),
1614 _node_potential(0), _blossom_potential(), _blossom_node_list(),
1615 _node_num(0), _blossom_num(0),
1617 _blossom_index(0), _blossom_set(0), _blossom_data(0),
1618 _node_index(0), _node_heap_index(0), _node_data(0),
1619 _tree_set_index(0), _tree_set(0),
1621 _delta1_index(0), _delta1(0),
1622 _delta2_index(0), _delta2(0),
1623 _delta3_index(0), _delta3(0),
1624 _delta4_index(0), _delta4(0),
1628 ~MaxWeightedMatching() {
1629 destroyStructures();
1632 /// \name Execution control
1633 /// The simplest way to execute the algorithm is to use the member
1634 /// \c run() member function.
1638 /// \brief Initialize the algorithm
1640 /// Initialize the algorithm
1644 for (ArcIt e(_graph); e != INVALID; ++e) {
1645 _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1647 for (NodeIt n(_graph); n != INVALID; ++n) {
1648 _delta1_index->set(n, _delta1->PRE_HEAP);
1650 for (EdgeIt e(_graph); e != INVALID; ++e) {
1651 _delta3_index->set(e, _delta3->PRE_HEAP);
1653 for (int i = 0; i < _blossom_num; ++i) {
1654 _delta2_index->set(i, _delta2->PRE_HEAP);
1655 _delta4_index->set(i, _delta4->PRE_HEAP);
1659 for (NodeIt n(_graph); n != INVALID; ++n) {
1661 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1662 if (_graph.target(e) == n) continue;
1663 if ((dualScale * _weight[e]) / 2 > max) {
1664 max = (dualScale * _weight[e]) / 2;
1667 _node_index->set(n, index);
1668 (*_node_data)[index].pot = max;
1669 _delta1->push(n, max);
1671 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1673 _tree_set->insert(blossom);
1675 (*_blossom_data)[blossom].status = EVEN;
1676 (*_blossom_data)[blossom].pred = INVALID;
1677 (*_blossom_data)[blossom].next = INVALID;
1678 (*_blossom_data)[blossom].pot = 0;
1679 (*_blossom_data)[blossom].offset = 0;
1682 for (EdgeIt e(_graph); e != INVALID; ++e) {
1683 int si = (*_node_index)[_graph.u(e)];
1684 int ti = (*_node_index)[_graph.v(e)];
1685 if (_graph.u(e) != _graph.v(e)) {
1686 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1687 dualScale * _weight[e]) / 2);
1692 /// \brief Starts the algorithm
1694 /// Starts the algorithm
1700 int unmatched = _node_num;
1701 while (unmatched > 0) {
1702 Value d1 = !_delta1->empty() ?
1703 _delta1->prio() : std::numeric_limits<Value>::max();
1705 Value d2 = !_delta2->empty() ?
1706 _delta2->prio() : std::numeric_limits<Value>::max();
1708 Value d3 = !_delta3->empty() ?
1709 _delta3->prio() : std::numeric_limits<Value>::max();
1711 Value d4 = !_delta4->empty() ?
1712 _delta4->prio() : std::numeric_limits<Value>::max();
1714 _delta_sum = d1; OpType ot = D1;
1715 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1716 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1717 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1723 Node n = _delta1->top();
1730 int blossom = _delta2->top();
1731 Node n = _blossom_set->classTop(blossom);
1732 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1738 Edge e = _delta3->top();
1740 int left_blossom = _blossom_set->find(_graph.u(e));
1741 int right_blossom = _blossom_set->find(_graph.v(e));
1743 if (left_blossom == right_blossom) {
1747 if ((*_blossom_data)[left_blossom].status == EVEN) {
1748 left_tree = _tree_set->find(left_blossom);
1754 if ((*_blossom_data)[right_blossom].status == EVEN) {
1755 right_tree = _tree_set->find(right_blossom);
1761 if (left_tree == right_tree) {
1762 shrinkOnEdge(e, left_tree);
1770 splitBlossom(_delta4->top());
1777 /// \brief Runs %MaxWeightedMatching algorithm.
1779 /// This method runs the %MaxWeightedMatching algorithm.
1781 /// \note mwm.run() is just a shortcut of the following code.
1793 /// \name Primal solution
1794 /// Functions for get the primal solution, ie. the matching.
1798 /// \brief Returns the matching value.
1800 /// Returns the matching value.
1801 Value matchingValue() const {
1803 for (NodeIt n(_graph); n != INVALID; ++n) {
1804 if ((*_matching)[n] != INVALID) {
1805 sum += _weight[(*_matching)[n]];
1811 /// \brief Returns the cardinality of the matching.
1813 /// Returns the cardinality of the matching.
1814 int matchingSize() const {
1816 for (NodeIt n(_graph); n != INVALID; ++n) {
1817 if ((*_matching)[n] != INVALID) {
1824 /// \brief Returns true when the edge is in the matching.
1826 /// Returns true when the edge is in the matching.
1827 bool matching(const Edge& edge) const {
1828 return edge == (*_matching)[_graph.u(edge)];
1831 /// \brief Returns the incident matching arc.
1833 /// Returns the incident matching arc from given node. If the
1834 /// node is not matched then it gives back \c INVALID.
1835 Arc matching(const Node& node) const {
1836 return (*_matching)[node];
1839 /// \brief Returns the mate of the node.
1841 /// Returns the adjancent node in a mathcing arc. If the node is
1842 /// not matched then it gives back \c INVALID.
1843 Node mate(const Node& node) const {
1844 return (*_matching)[node] != INVALID ?
1845 _graph.target((*_matching)[node]) : INVALID;
1850 /// \name Dual solution
1851 /// Functions for get the dual solution.
1855 /// \brief Returns the value of the dual solution.
1857 /// Returns the value of the dual solution. It should be equal to
1858 /// the primal value scaled by \ref dualScale "dual scale".
1859 Value dualValue() const {
1861 for (NodeIt n(_graph); n != INVALID; ++n) {
1862 sum += nodeValue(n);
1864 for (int i = 0; i < blossomNum(); ++i) {
1865 sum += blossomValue(i) * (blossomSize(i) / 2);
1870 /// \brief Returns the value of the node.
1872 /// Returns the the value of the node.
1873 Value nodeValue(const Node& n) const {
1874 return (*_node_potential)[n];
1877 /// \brief Returns the number of the blossoms in the basis.
1879 /// Returns the number of the blossoms in the basis.
1881 int blossomNum() const {
1882 return _blossom_potential.size();
1886 /// \brief Returns the number of the nodes in the blossom.
1888 /// Returns the number of the nodes in the blossom.
1889 int blossomSize(int k) const {
1890 return _blossom_potential[k].end - _blossom_potential[k].begin;
1893 /// \brief Returns the value of the blossom.
1895 /// Returns the the value of the blossom.
1897 Value blossomValue(int k) const {
1898 return _blossom_potential[k].value;
1901 /// \brief Lemon iterator for get the items of the blossom.
1903 /// Lemon iterator for get the nodes of the blossom. This class
1904 /// provides a common style lemon iterator which gives back a
1905 /// subset of the nodes.
1909 /// \brief Constructor.
1911 /// Constructor for get the nodes of the variable.
1912 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1913 : _algorithm(&algorithm)
1915 _index = _algorithm->_blossom_potential[variable].begin;
1916 _last = _algorithm->_blossom_potential[variable].end;
1919 /// \brief Conversion to node.
1921 /// Conversion to node.
1922 operator Node() const {
1923 return _algorithm->_blossom_node_list[_index];
1926 /// \brief Increment operator.
1928 /// Increment operator.
1929 BlossomIt& operator++() {
1934 /// \brief Validity checking
1936 /// Checks whether the iterator is invalid.
1937 bool operator==(Invalid) const { return _index == _last; }
1939 /// \brief Validity checking
1941 /// Checks whether the iterator is valid.
1942 bool operator!=(Invalid) const { return _index != _last; }
1945 const MaxWeightedMatching* _algorithm;
1954 /// \ingroup matching
1956 /// \brief Weighted perfect matching in general graphs
1958 /// This class provides an efficient implementation of Edmond's
1959 /// maximum weighted perfect matching algorithm. The implementation
1960 /// is based on extensive use of priority queues and provides
1961 /// \f$O(nm\log(n))\f$ time complexity.
1963 /// The maximum weighted matching problem is to find undirected
1964 /// edges in the graph with maximum overall weight and no two of
1965 /// them shares their ends and covers all nodes. The problem can be
1966 /// formulated with the following linear program.
1967 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1968 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1969 \quad \forall B\in\mathcal{O}\f] */
1970 /// \f[x_e \ge 0\quad \forall e\in E\f]
1971 /// \f[\max \sum_{e\in E}x_ew_e\f]
1972 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1973 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1974 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1975 /// subsets of the nodes.
1977 /// The algorithm calculates an optimal matching and a proof of the
1978 /// optimality. The solution of the dual problem can be used to check
1979 /// the result of the algorithm. The dual linear problem is the
1980 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1981 w_{uv} \quad \forall uv\in E\f] */
1982 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1983 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1984 \frac{\vert B \vert - 1}{2}z_B\f] */
1986 /// The algorithm can be executed with \c run() or the \c init() and
1987 /// then the \c start() member functions. After it the matching can
1988 /// be asked with \c matching() or mate() functions. The dual
1989 /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1990 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1991 /// "BlossomIt" nested class which is able to iterate on the nodes
1992 /// of a blossom. If the value type is integral then the dual
1993 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1994 template <typename _Graph,
1995 typename _WeightMap = typename _Graph::template EdgeMap<int> >
1996 class MaxWeightedPerfectMatching {
1999 typedef _Graph Graph;
2000 typedef _WeightMap WeightMap;
2001 typedef typename WeightMap::Value Value;
2003 /// \brief Scaling factor for dual solution
2005 /// Scaling factor for dual solution, it is equal to 4 or 1
2006 /// according to the value type.
2007 static const int dualScale =
2008 std::numeric_limits<Value>::is_integer ? 4 : 1;
2010 typedef typename Graph::template NodeMap<typename Graph::Arc>
2015 TEMPLATE_GRAPH_TYPEDEFS(Graph);
2017 typedef typename Graph::template NodeMap<Value> NodePotential;
2018 typedef std::vector<Node> BlossomNodeList;
2020 struct BlossomVariable {
2024 BlossomVariable(int _begin, int _end, Value _value)
2025 : begin(_begin), end(_end), value(_value) {}
2029 typedef std::vector<BlossomVariable> BlossomPotential;
2031 const Graph& _graph;
2032 const WeightMap& _weight;
2034 MatchingMap* _matching;
2036 NodePotential* _node_potential;
2038 BlossomPotential _blossom_potential;
2039 BlossomNodeList _blossom_node_list;
2044 typedef RangeMap<int> IntIntMap;
2047 EVEN = -1, MATCHED = 0, ODD = 1
2050 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2051 struct BlossomData {
2058 IntNodeMap *_blossom_index;
2059 BlossomSet *_blossom_set;
2060 RangeMap<BlossomData>* _blossom_data;
2062 IntNodeMap *_node_index;
2063 IntArcMap *_node_heap_index;
2067 NodeData(IntArcMap& node_heap_index)
2068 : heap(node_heap_index) {}
2072 BinHeap<Value, IntArcMap> heap;
2073 std::map<int, Arc> heap_index;
2078 RangeMap<NodeData>* _node_data;
2080 typedef ExtendFindEnum<IntIntMap> TreeSet;
2082 IntIntMap *_tree_set_index;
2085 IntIntMap *_delta2_index;
2086 BinHeap<Value, IntIntMap> *_delta2;
2088 IntEdgeMap *_delta3_index;
2089 BinHeap<Value, IntEdgeMap> *_delta3;
2091 IntIntMap *_delta4_index;
2092 BinHeap<Value, IntIntMap> *_delta4;
2096 void createStructures() {
2097 _node_num = countNodes(_graph);
2098 _blossom_num = _node_num * 3 / 2;
2101 _matching = new MatchingMap(_graph);
2103 if (!_node_potential) {
2104 _node_potential = new NodePotential(_graph);
2106 if (!_blossom_set) {
2107 _blossom_index = new IntNodeMap(_graph);
2108 _blossom_set = new BlossomSet(*_blossom_index);
2109 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2113 _node_index = new IntNodeMap(_graph);
2114 _node_heap_index = new IntArcMap(_graph);
2115 _node_data = new RangeMap<NodeData>(_node_num,
2116 NodeData(*_node_heap_index));
2120 _tree_set_index = new IntIntMap(_blossom_num);
2121 _tree_set = new TreeSet(*_tree_set_index);
2124 _delta2_index = new IntIntMap(_blossom_num);
2125 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2128 _delta3_index = new IntEdgeMap(_graph);
2129 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2132 _delta4_index = new IntIntMap(_blossom_num);
2133 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2137 void destroyStructures() {
2138 _node_num = countNodes(_graph);
2139 _blossom_num = _node_num * 3 / 2;
2144 if (_node_potential) {
2145 delete _node_potential;
2148 delete _blossom_index;
2149 delete _blossom_set;
2150 delete _blossom_data;
2155 delete _node_heap_index;
2160 delete _tree_set_index;
2164 delete _delta2_index;
2168 delete _delta3_index;
2172 delete _delta4_index;
2177 void matchedToEven(int blossom, int tree) {
2178 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2179 _delta2->erase(blossom);
2182 if (!_blossom_set->trivial(blossom)) {
2183 (*_blossom_data)[blossom].pot -=
2184 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2187 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2188 n != INVALID; ++n) {
2190 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2191 int ni = (*_node_index)[n];
2193 (*_node_data)[ni].heap.clear();
2194 (*_node_data)[ni].heap_index.clear();
2196 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2198 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2199 Node v = _graph.source(e);
2200 int vb = _blossom_set->find(v);
2201 int vi = (*_node_index)[v];
2203 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2204 dualScale * _weight[e];
2206 if ((*_blossom_data)[vb].status == EVEN) {
2207 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2208 _delta3->push(e, rw / 2);
2211 typename std::map<int, Arc>::iterator it =
2212 (*_node_data)[vi].heap_index.find(tree);
2214 if (it != (*_node_data)[vi].heap_index.end()) {
2215 if ((*_node_data)[vi].heap[it->second] > rw) {
2216 (*_node_data)[vi].heap.replace(it->second, e);
2217 (*_node_data)[vi].heap.decrease(e, rw);
2221 (*_node_data)[vi].heap.push(e, rw);
2222 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2225 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2226 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2228 if ((*_blossom_data)[vb].status == MATCHED) {
2229 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2230 _delta2->push(vb, _blossom_set->classPrio(vb) -
2231 (*_blossom_data)[vb].offset);
2232 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2233 (*_blossom_data)[vb].offset){
2234 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2235 (*_blossom_data)[vb].offset);
2242 (*_blossom_data)[blossom].offset = 0;
2245 void matchedToOdd(int blossom) {
2246 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2247 _delta2->erase(blossom);
2249 (*_blossom_data)[blossom].offset += _delta_sum;
2250 if (!_blossom_set->trivial(blossom)) {
2251 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2252 (*_blossom_data)[blossom].offset);
2256 void evenToMatched(int blossom, int tree) {
2257 if (!_blossom_set->trivial(blossom)) {
2258 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2261 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2262 n != INVALID; ++n) {
2263 int ni = (*_node_index)[n];
2264 (*_node_data)[ni].pot -= _delta_sum;
2266 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2267 Node v = _graph.source(e);
2268 int vb = _blossom_set->find(v);
2269 int vi = (*_node_index)[v];
2271 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2272 dualScale * _weight[e];
2274 if (vb == blossom) {
2275 if (_delta3->state(e) == _delta3->IN_HEAP) {
2278 } else if ((*_blossom_data)[vb].status == EVEN) {
2280 if (_delta3->state(e) == _delta3->IN_HEAP) {
2284 int vt = _tree_set->find(vb);
2288 Arc r = _graph.oppositeArc(e);
2290 typename std::map<int, Arc>::iterator it =
2291 (*_node_data)[ni].heap_index.find(vt);
2293 if (it != (*_node_data)[ni].heap_index.end()) {
2294 if ((*_node_data)[ni].heap[it->second] > rw) {
2295 (*_node_data)[ni].heap.replace(it->second, r);
2296 (*_node_data)[ni].heap.decrease(r, rw);
2300 (*_node_data)[ni].heap.push(r, rw);
2301 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2304 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2305 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2307 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2308 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2309 (*_blossom_data)[blossom].offset);
2310 } else if ((*_delta2)[blossom] >
2311 _blossom_set->classPrio(blossom) -
2312 (*_blossom_data)[blossom].offset){
2313 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2314 (*_blossom_data)[blossom].offset);
2320 typename std::map<int, Arc>::iterator it =
2321 (*_node_data)[vi].heap_index.find(tree);
2323 if (it != (*_node_data)[vi].heap_index.end()) {
2324 (*_node_data)[vi].heap.erase(it->second);
2325 (*_node_data)[vi].heap_index.erase(it);
2326 if ((*_node_data)[vi].heap.empty()) {
2327 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2328 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2329 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2332 if ((*_blossom_data)[vb].status == MATCHED) {
2333 if (_blossom_set->classPrio(vb) ==
2334 std::numeric_limits<Value>::max()) {
2336 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2337 (*_blossom_data)[vb].offset) {
2338 _delta2->increase(vb, _blossom_set->classPrio(vb) -
2339 (*_blossom_data)[vb].offset);
2348 void oddToMatched(int blossom) {
2349 (*_blossom_data)[blossom].offset -= _delta_sum;
2351 if (_blossom_set->classPrio(blossom) !=
2352 std::numeric_limits<Value>::max()) {
2353 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2354 (*_blossom_data)[blossom].offset);
2357 if (!_blossom_set->trivial(blossom)) {
2358 _delta4->erase(blossom);
2362 void oddToEven(int blossom, int tree) {
2363 if (!_blossom_set->trivial(blossom)) {
2364 _delta4->erase(blossom);
2365 (*_blossom_data)[blossom].pot -=
2366 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2369 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2370 n != INVALID; ++n) {
2371 int ni = (*_node_index)[n];
2373 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2375 (*_node_data)[ni].heap.clear();
2376 (*_node_data)[ni].heap_index.clear();
2377 (*_node_data)[ni].pot +=
2378 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2380 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2381 Node v = _graph.source(e);
2382 int vb = _blossom_set->find(v);
2383 int vi = (*_node_index)[v];
2385 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2386 dualScale * _weight[e];
2388 if ((*_blossom_data)[vb].status == EVEN) {
2389 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2390 _delta3->push(e, rw / 2);
2394 typename std::map<int, Arc>::iterator it =
2395 (*_node_data)[vi].heap_index.find(tree);
2397 if (it != (*_node_data)[vi].heap_index.end()) {
2398 if ((*_node_data)[vi].heap[it->second] > rw) {
2399 (*_node_data)[vi].heap.replace(it->second, e);
2400 (*_node_data)[vi].heap.decrease(e, rw);
2404 (*_node_data)[vi].heap.push(e, rw);
2405 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2408 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2409 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2411 if ((*_blossom_data)[vb].status == MATCHED) {
2412 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2413 _delta2->push(vb, _blossom_set->classPrio(vb) -
2414 (*_blossom_data)[vb].offset);
2415 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2416 (*_blossom_data)[vb].offset) {
2417 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2418 (*_blossom_data)[vb].offset);
2425 (*_blossom_data)[blossom].offset = 0;
2428 void alternatePath(int even, int tree) {
2431 evenToMatched(even, tree);
2432 (*_blossom_data)[even].status = MATCHED;
2434 while ((*_blossom_data)[even].pred != INVALID) {
2435 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2436 (*_blossom_data)[odd].status = MATCHED;
2438 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2440 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2441 (*_blossom_data)[even].status = MATCHED;
2442 evenToMatched(even, tree);
2443 (*_blossom_data)[even].next =
2444 _graph.oppositeArc((*_blossom_data)[odd].pred);
2449 void destroyTree(int tree) {
2450 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2451 if ((*_blossom_data)[b].status == EVEN) {
2452 (*_blossom_data)[b].status = MATCHED;
2453 evenToMatched(b, tree);
2454 } else if ((*_blossom_data)[b].status == ODD) {
2455 (*_blossom_data)[b].status = MATCHED;
2459 _tree_set->eraseClass(tree);
2462 void augmentOnEdge(const Edge& edge) {
2464 int left = _blossom_set->find(_graph.u(edge));
2465 int right = _blossom_set->find(_graph.v(edge));
2467 int left_tree = _tree_set->find(left);
2468 alternatePath(left, left_tree);
2469 destroyTree(left_tree);
2471 int right_tree = _tree_set->find(right);
2472 alternatePath(right, right_tree);
2473 destroyTree(right_tree);
2475 (*_blossom_data)[left].next = _graph.direct(edge, true);
2476 (*_blossom_data)[right].next = _graph.direct(edge, false);
2479 void extendOnArc(const Arc& arc) {
2480 int base = _blossom_set->find(_graph.target(arc));
2481 int tree = _tree_set->find(base);
2483 int odd = _blossom_set->find(_graph.source(arc));
2484 _tree_set->insert(odd, tree);
2485 (*_blossom_data)[odd].status = ODD;
2487 (*_blossom_data)[odd].pred = arc;
2489 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2490 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2491 _tree_set->insert(even, tree);
2492 (*_blossom_data)[even].status = EVEN;
2493 matchedToEven(even, tree);
2496 void shrinkOnEdge(const Edge& edge, int tree) {
2498 std::vector<int> left_path, right_path;
2501 std::set<int> left_set, right_set;
2502 int left = _blossom_set->find(_graph.u(edge));
2503 left_path.push_back(left);
2504 left_set.insert(left);
2506 int right = _blossom_set->find(_graph.v(edge));
2507 right_path.push_back(right);
2508 right_set.insert(right);
2512 if ((*_blossom_data)[left].pred == INVALID) break;
2515 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2516 left_path.push_back(left);
2518 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2519 left_path.push_back(left);
2521 left_set.insert(left);
2523 if (right_set.find(left) != right_set.end()) {
2528 if ((*_blossom_data)[right].pred == INVALID) break;
2531 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2532 right_path.push_back(right);
2534 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2535 right_path.push_back(right);
2537 right_set.insert(right);
2539 if (left_set.find(right) != left_set.end()) {
2547 if ((*_blossom_data)[left].pred == INVALID) {
2549 while (left_set.find(nca) == left_set.end()) {
2551 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2552 right_path.push_back(nca);
2554 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2555 right_path.push_back(nca);
2559 while (right_set.find(nca) == right_set.end()) {
2561 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2562 left_path.push_back(nca);
2564 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2565 left_path.push_back(nca);
2571 std::vector<int> subblossoms;
2574 prev = _graph.direct(edge, true);
2575 for (int i = 0; left_path[i] != nca; i += 2) {
2576 subblossoms.push_back(left_path[i]);
2577 (*_blossom_data)[left_path[i]].next = prev;
2578 _tree_set->erase(left_path[i]);
2580 subblossoms.push_back(left_path[i + 1]);
2581 (*_blossom_data)[left_path[i + 1]].status = EVEN;
2582 oddToEven(left_path[i + 1], tree);
2583 _tree_set->erase(left_path[i + 1]);
2584 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2588 while (right_path[k] != nca) ++k;
2590 subblossoms.push_back(nca);
2591 (*_blossom_data)[nca].next = prev;
2593 for (int i = k - 2; i >= 0; i -= 2) {
2594 subblossoms.push_back(right_path[i + 1]);
2595 (*_blossom_data)[right_path[i + 1]].status = EVEN;
2596 oddToEven(right_path[i + 1], tree);
2597 _tree_set->erase(right_path[i + 1]);
2599 (*_blossom_data)[right_path[i + 1]].next =
2600 (*_blossom_data)[right_path[i + 1]].pred;
2602 subblossoms.push_back(right_path[i]);
2603 _tree_set->erase(right_path[i]);
2607 _blossom_set->join(subblossoms.begin(), subblossoms.end());
2609 for (int i = 0; i < int(subblossoms.size()); ++i) {
2610 if (!_blossom_set->trivial(subblossoms[i])) {
2611 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2613 (*_blossom_data)[subblossoms[i]].status = MATCHED;
2616 (*_blossom_data)[surface].pot = -2 * _delta_sum;
2617 (*_blossom_data)[surface].offset = 0;
2618 (*_blossom_data)[surface].status = EVEN;
2619 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2620 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2622 _tree_set->insert(surface, tree);
2623 _tree_set->erase(nca);
2626 void splitBlossom(int blossom) {
2627 Arc next = (*_blossom_data)[blossom].next;
2628 Arc pred = (*_blossom_data)[blossom].pred;
2630 int tree = _tree_set->find(blossom);
2632 (*_blossom_data)[blossom].status = MATCHED;
2633 oddToMatched(blossom);
2634 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2635 _delta2->erase(blossom);
2638 std::vector<int> subblossoms;
2639 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2641 Value offset = (*_blossom_data)[blossom].offset;
2642 int b = _blossom_set->find(_graph.source(pred));
2643 int d = _blossom_set->find(_graph.source(next));
2645 int ib = -1, id = -1;
2646 for (int i = 0; i < int(subblossoms.size()); ++i) {
2647 if (subblossoms[i] == b) ib = i;
2648 if (subblossoms[i] == d) id = i;
2650 (*_blossom_data)[subblossoms[i]].offset = offset;
2651 if (!_blossom_set->trivial(subblossoms[i])) {
2652 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2654 if (_blossom_set->classPrio(subblossoms[i]) !=
2655 std::numeric_limits<Value>::max()) {
2656 _delta2->push(subblossoms[i],
2657 _blossom_set->classPrio(subblossoms[i]) -
2658 (*_blossom_data)[subblossoms[i]].offset);
2662 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2663 for (int i = (id + 1) % subblossoms.size();
2664 i != ib; i = (i + 2) % subblossoms.size()) {
2665 int sb = subblossoms[i];
2666 int tb = subblossoms[(i + 1) % subblossoms.size()];
2667 (*_blossom_data)[sb].next =
2668 _graph.oppositeArc((*_blossom_data)[tb].next);
2671 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2672 int sb = subblossoms[i];
2673 int tb = subblossoms[(i + 1) % subblossoms.size()];
2674 int ub = subblossoms[(i + 2) % subblossoms.size()];
2676 (*_blossom_data)[sb].status = ODD;
2678 _tree_set->insert(sb, tree);
2679 (*_blossom_data)[sb].pred = pred;
2680 (*_blossom_data)[sb].next =
2681 _graph.oppositeArc((*_blossom_data)[tb].next);
2683 pred = (*_blossom_data)[ub].next;
2685 (*_blossom_data)[tb].status = EVEN;
2686 matchedToEven(tb, tree);
2687 _tree_set->insert(tb, tree);
2688 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2691 (*_blossom_data)[subblossoms[id]].status = ODD;
2692 matchedToOdd(subblossoms[id]);
2693 _tree_set->insert(subblossoms[id], tree);
2694 (*_blossom_data)[subblossoms[id]].next = next;
2695 (*_blossom_data)[subblossoms[id]].pred = pred;
2699 for (int i = (ib + 1) % subblossoms.size();
2700 i != id; i = (i + 2) % subblossoms.size()) {
2701 int sb = subblossoms[i];
2702 int tb = subblossoms[(i + 1) % subblossoms.size()];
2703 (*_blossom_data)[sb].next =
2704 _graph.oppositeArc((*_blossom_data)[tb].next);
2707 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2708 int sb = subblossoms[i];
2709 int tb = subblossoms[(i + 1) % subblossoms.size()];
2710 int ub = subblossoms[(i + 2) % subblossoms.size()];
2712 (*_blossom_data)[sb].status = ODD;
2714 _tree_set->insert(sb, tree);
2715 (*_blossom_data)[sb].next = next;
2716 (*_blossom_data)[sb].pred =
2717 _graph.oppositeArc((*_blossom_data)[tb].next);
2719 (*_blossom_data)[tb].status = EVEN;
2720 matchedToEven(tb, tree);
2721 _tree_set->insert(tb, tree);
2722 (*_blossom_data)[tb].pred =
2723 (*_blossom_data)[tb].next =
2724 _graph.oppositeArc((*_blossom_data)[ub].next);
2725 next = (*_blossom_data)[ub].next;
2728 (*_blossom_data)[subblossoms[ib]].status = ODD;
2729 matchedToOdd(subblossoms[ib]);
2730 _tree_set->insert(subblossoms[ib], tree);
2731 (*_blossom_data)[subblossoms[ib]].next = next;
2732 (*_blossom_data)[subblossoms[ib]].pred = pred;
2734 _tree_set->erase(blossom);
2737 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2738 if (_blossom_set->trivial(blossom)) {
2739 int bi = (*_node_index)[base];
2740 Value pot = (*_node_data)[bi].pot;
2742 _matching->set(base, matching);
2743 _blossom_node_list.push_back(base);
2744 _node_potential->set(base, pot);
2747 Value pot = (*_blossom_data)[blossom].pot;
2748 int bn = _blossom_node_list.size();
2750 std::vector<int> subblossoms;
2751 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2752 int b = _blossom_set->find(base);
2754 for (int i = 0; i < int(subblossoms.size()); ++i) {
2755 if (subblossoms[i] == b) { ib = i; break; }
2758 for (int i = 1; i < int(subblossoms.size()); i += 2) {
2759 int sb = subblossoms[(ib + i) % subblossoms.size()];
2760 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2762 Arc m = (*_blossom_data)[tb].next;
2763 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2764 extractBlossom(tb, _graph.source(m), m);
2766 extractBlossom(subblossoms[ib], base, matching);
2768 int en = _blossom_node_list.size();
2770 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2774 void extractMatching() {
2775 std::vector<int> blossoms;
2776 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2777 blossoms.push_back(c);
2780 for (int i = 0; i < int(blossoms.size()); ++i) {
2782 Value offset = (*_blossom_data)[blossoms[i]].offset;
2783 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2784 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2785 n != INVALID; ++n) {
2786 (*_node_data)[(*_node_index)[n]].pot -= offset;
2789 Arc matching = (*_blossom_data)[blossoms[i]].next;
2790 Node base = _graph.source(matching);
2791 extractBlossom(blossoms[i], base, matching);
2797 /// \brief Constructor
2800 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2801 : _graph(graph), _weight(weight), _matching(0),
2802 _node_potential(0), _blossom_potential(), _blossom_node_list(),
2803 _node_num(0), _blossom_num(0),
2805 _blossom_index(0), _blossom_set(0), _blossom_data(0),
2806 _node_index(0), _node_heap_index(0), _node_data(0),
2807 _tree_set_index(0), _tree_set(0),
2809 _delta2_index(0), _delta2(0),
2810 _delta3_index(0), _delta3(0),
2811 _delta4_index(0), _delta4(0),
2815 ~MaxWeightedPerfectMatching() {
2816 destroyStructures();
2819 /// \name Execution control
2820 /// The simplest way to execute the algorithm is to use the member
2821 /// \c run() member function.
2825 /// \brief Initialize the algorithm
2827 /// Initialize the algorithm
2831 for (ArcIt e(_graph); e != INVALID; ++e) {
2832 _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
2834 for (EdgeIt e(_graph); e != INVALID; ++e) {
2835 _delta3_index->set(e, _delta3->PRE_HEAP);
2837 for (int i = 0; i < _blossom_num; ++i) {
2838 _delta2_index->set(i, _delta2->PRE_HEAP);
2839 _delta4_index->set(i, _delta4->PRE_HEAP);
2843 for (NodeIt n(_graph); n != INVALID; ++n) {
2844 Value max = - std::numeric_limits<Value>::max();
2845 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2846 if (_graph.target(e) == n) continue;
2847 if ((dualScale * _weight[e]) / 2 > max) {
2848 max = (dualScale * _weight[e]) / 2;
2851 _node_index->set(n, index);
2852 (*_node_data)[index].pot = max;
2854 _blossom_set->insert(n, std::numeric_limits<Value>::max());
2856 _tree_set->insert(blossom);
2858 (*_blossom_data)[blossom].status = EVEN;
2859 (*_blossom_data)[blossom].pred = INVALID;
2860 (*_blossom_data)[blossom].next = INVALID;
2861 (*_blossom_data)[blossom].pot = 0;
2862 (*_blossom_data)[blossom].offset = 0;
2865 for (EdgeIt e(_graph); e != INVALID; ++e) {
2866 int si = (*_node_index)[_graph.u(e)];
2867 int ti = (*_node_index)[_graph.v(e)];
2868 if (_graph.u(e) != _graph.v(e)) {
2869 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2870 dualScale * _weight[e]) / 2);
2875 /// \brief Starts the algorithm
2877 /// Starts the algorithm
2883 int unmatched = _node_num;
2884 while (unmatched > 0) {
2885 Value d2 = !_delta2->empty() ?
2886 _delta2->prio() : std::numeric_limits<Value>::max();
2888 Value d3 = !_delta3->empty() ?
2889 _delta3->prio() : std::numeric_limits<Value>::max();
2891 Value d4 = !_delta4->empty() ?
2892 _delta4->prio() : std::numeric_limits<Value>::max();
2894 _delta_sum = d2; OpType ot = D2;
2895 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2896 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2898 if (_delta_sum == std::numeric_limits<Value>::max()) {
2905 int blossom = _delta2->top();
2906 Node n = _blossom_set->classTop(blossom);
2907 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2913 Edge e = _delta3->top();
2915 int left_blossom = _blossom_set->find(_graph.u(e));
2916 int right_blossom = _blossom_set->find(_graph.v(e));
2918 if (left_blossom == right_blossom) {
2921 int left_tree = _tree_set->find(left_blossom);
2922 int right_tree = _tree_set->find(right_blossom);
2924 if (left_tree == right_tree) {
2925 shrinkOnEdge(e, left_tree);
2933 splitBlossom(_delta4->top());
2941 /// \brief Runs %MaxWeightedPerfectMatching algorithm.
2943 /// This method runs the %MaxWeightedPerfectMatching algorithm.
2945 /// \note mwm.run() is just a shortcut of the following code.
2957 /// \name Primal solution
2958 /// Functions for get the primal solution, ie. the matching.
2962 /// \brief Returns the matching value.
2964 /// Returns the matching value.
2965 Value matchingValue() const {
2967 for (NodeIt n(_graph); n != INVALID; ++n) {
2968 if ((*_matching)[n] != INVALID) {
2969 sum += _weight[(*_matching)[n]];
2975 /// \brief Returns true when the edge is in the matching.
2977 /// Returns true when the edge is in the matching.
2978 bool matching(const Edge& edge) const {
2979 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
2982 /// \brief Returns the incident matching edge.
2984 /// Returns the incident matching arc from given edge.
2985 Arc matching(const Node& node) const {
2986 return (*_matching)[node];
2989 /// \brief Returns the mate of the node.
2991 /// Returns the adjancent node in a mathcing arc.
2992 Node mate(const Node& node) const {
2993 return _graph.target((*_matching)[node]);
2998 /// \name Dual solution
2999 /// Functions for get the dual solution.
3003 /// \brief Returns the value of the dual solution.
3005 /// Returns the value of the dual solution. It should be equal to
3006 /// the primal value scaled by \ref dualScale "dual scale".
3007 Value dualValue() const {
3009 for (NodeIt n(_graph); n != INVALID; ++n) {
3010 sum += nodeValue(n);
3012 for (int i = 0; i < blossomNum(); ++i) {
3013 sum += blossomValue(i) * (blossomSize(i) / 2);
3018 /// \brief Returns the value of the node.
3020 /// Returns the the value of the node.
3021 Value nodeValue(const Node& n) const {
3022 return (*_node_potential)[n];
3025 /// \brief Returns the number of the blossoms in the basis.
3027 /// Returns the number of the blossoms in the basis.
3029 int blossomNum() const {
3030 return _blossom_potential.size();
3034 /// \brief Returns the number of the nodes in the blossom.
3036 /// Returns the number of the nodes in the blossom.
3037 int blossomSize(int k) const {
3038 return _blossom_potential[k].end - _blossom_potential[k].begin;
3041 /// \brief Returns the value of the blossom.
3043 /// Returns the the value of the blossom.
3045 Value blossomValue(int k) const {
3046 return _blossom_potential[k].value;
3049 /// \brief Lemon iterator for get the items of the blossom.
3051 /// Lemon iterator for get the nodes of the blossom. This class
3052 /// provides a common style lemon iterator which gives back a
3053 /// subset of the nodes.
3057 /// \brief Constructor.
3059 /// Constructor for get the nodes of the variable.
3060 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3061 : _algorithm(&algorithm)
3063 _index = _algorithm->_blossom_potential[variable].begin;
3064 _last = _algorithm->_blossom_potential[variable].end;
3067 /// \brief Conversion to node.
3069 /// Conversion to node.
3070 operator Node() const {
3071 return _algorithm->_blossom_node_list[_index];
3074 /// \brief Increment operator.
3076 /// Increment operator.
3077 BlossomIt& operator++() {
3082 /// \brief Validity checking
3084 /// Checks whether the iterator is invalid.
3085 bool operator==(Invalid) const { return _index == _last; }
3087 /// \brief Validity checking
3089 /// Checks whether the iterator is valid.
3090 bool operator!=(Invalid) const { return _index != _last; }
3093 const MaxWeightedPerfectMatching* _algorithm;
3103 } //END OF NAMESPACE LEMON
3105 #endif //LEMON_MAX_MATCHING_H