1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2009
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_MATCHING_H
20 #define LEMON_MATCHING_H
27 #include <lemon/core.h>
28 #include <lemon/unionfind.h>
29 #include <lemon/bin_heap.h>
30 #include <lemon/maps.h>
31 #include <lemon/fractional_matching.h>
35 ///\brief Maximum matching algorithms in general graphs.
41 /// \brief Maximum cardinality matching in general graphs
43 /// This class implements Edmonds' alternating forest matching algorithm
44 /// for finding a maximum cardinality matching in a general undirected graph.
45 /// It can be started from an arbitrary initial matching
46 /// (the default is the empty one).
48 /// The dual solution of the problem is a map of the nodes to
49 /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
50 /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
51 /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
52 /// with factor-critical components, the nodes in \c ODD/A form the
53 /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
54 /// a perfect matching. The number of the factor-critical components
55 /// minus the number of barrier nodes is a lower bound on the
56 /// unmatched nodes, and the matching is optimal if and only if this bound is
57 /// tight. This decomposition can be obtained using \ref status() or
58 /// \ref statusMap() after running the algorithm.
60 /// \tparam GR The undirected graph type the algorithm runs on.
61 template <typename GR>
65 /// The graph type of the algorithm
67 /// The type of the matching map
68 typedef typename Graph::template NodeMap<typename Graph::Arc>
71 ///\brief Status constants for Gallai-Edmonds decomposition.
73 ///These constants are used for indicating the Gallai-Edmonds
74 ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
75 ///induce a subgraph with factor-critical components, the nodes with
76 ///status \c ODD (or \c A) form the canonical barrier, and the nodes
77 ///with status \c MATCHED (or \c C) induce a subgraph having a
80 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.)
82 MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.)
84 ODD = -1, ///< = -1. (\c A is an alias for \c ODD.)
86 UNMATCHED = -2 ///< = -2.
89 /// The type of the status map
90 typedef typename Graph::template NodeMap<Status> StatusMap;
94 TEMPLATE_GRAPH_TYPEDEFS(Graph);
96 typedef UnionFindEnum<IntNodeMap> BlossomSet;
97 typedef ExtendFindEnum<IntNodeMap> TreeSet;
98 typedef RangeMap<Node> NodeIntMap;
99 typedef MatchingMap EarMap;
100 typedef std::vector<Node> NodeQueue;
103 MatchingMap* _matching;
108 IntNodeMap* _blossom_set_index;
109 BlossomSet* _blossom_set;
110 NodeIntMap* _blossom_rep;
112 IntNodeMap* _tree_set_index;
115 NodeQueue _node_queue;
116 int _process, _postpone, _last;
122 void createStructures() {
123 _node_num = countNodes(_graph);
125 _matching = new MatchingMap(_graph);
128 _status = new StatusMap(_graph);
131 _ear = new EarMap(_graph);
134 _blossom_set_index = new IntNodeMap(_graph);
135 _blossom_set = new BlossomSet(*_blossom_set_index);
138 _blossom_rep = new NodeIntMap(_node_num);
141 _tree_set_index = new IntNodeMap(_graph);
142 _tree_set = new TreeSet(*_tree_set_index);
144 _node_queue.resize(_node_num);
147 void destroyStructures() {
159 delete _blossom_set_index;
165 delete _tree_set_index;
170 void processDense(const Node& n) {
171 _process = _postpone = _last = 0;
172 _node_queue[_last++] = n;
174 while (_process != _last) {
175 Node u = _node_queue[_process++];
176 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
177 Node v = _graph.target(a);
178 if ((*_status)[v] == MATCHED) {
180 } else if ((*_status)[v] == UNMATCHED) {
187 while (_postpone != _last) {
188 Node u = _node_queue[_postpone++];
190 for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
191 Node v = _graph.target(a);
193 if ((*_status)[v] == EVEN) {
194 if (_blossom_set->find(u) != _blossom_set->find(v)) {
199 while (_process != _last) {
200 Node w = _node_queue[_process++];
201 for (OutArcIt b(_graph, w); b != INVALID; ++b) {
202 Node x = _graph.target(b);
203 if ((*_status)[x] == MATCHED) {
205 } else if ((*_status)[x] == UNMATCHED) {
215 void processSparse(const Node& n) {
216 _process = _last = 0;
217 _node_queue[_last++] = n;
218 while (_process != _last) {
219 Node u = _node_queue[_process++];
220 for (OutArcIt a(_graph, u); a != INVALID; ++a) {
221 Node v = _graph.target(a);
223 if ((*_status)[v] == EVEN) {
224 if (_blossom_set->find(u) != _blossom_set->find(v)) {
227 } else if ((*_status)[v] == MATCHED) {
229 } else if ((*_status)[v] == UNMATCHED) {
237 void shrinkOnEdge(const Edge& e) {
241 std::set<Node> left_set, right_set;
243 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
244 left_set.insert(left);
246 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
247 right_set.insert(right);
250 if ((*_matching)[left] == INVALID) break;
251 left = _graph.target((*_matching)[left]);
252 left = (*_blossom_rep)[_blossom_set->
253 find(_graph.target((*_ear)[left]))];
254 if (right_set.find(left) != right_set.end()) {
258 left_set.insert(left);
260 if ((*_matching)[right] == INVALID) break;
261 right = _graph.target((*_matching)[right]);
262 right = (*_blossom_rep)[_blossom_set->
263 find(_graph.target((*_ear)[right]))];
264 if (left_set.find(right) != left_set.end()) {
268 right_set.insert(right);
271 if (nca == INVALID) {
272 if ((*_matching)[left] == INVALID) {
274 while (left_set.find(nca) == left_set.end()) {
275 nca = _graph.target((*_matching)[nca]);
276 nca =(*_blossom_rep)[_blossom_set->
277 find(_graph.target((*_ear)[nca]))];
281 while (right_set.find(nca) == right_set.end()) {
282 nca = _graph.target((*_matching)[nca]);
283 nca = (*_blossom_rep)[_blossom_set->
284 find(_graph.target((*_ear)[nca]))];
292 Node node = _graph.u(e);
293 Arc arc = _graph.direct(e, true);
294 Node base = (*_blossom_rep)[_blossom_set->find(node)];
296 while (base != nca) {
301 n = _graph.target((*_matching)[n]);
303 n = _graph.target(a);
304 (*_ear)[n] = _graph.oppositeArc(a);
306 node = _graph.target((*_matching)[base]);
307 _tree_set->erase(base);
308 _tree_set->erase(node);
309 _blossom_set->insert(node, _blossom_set->find(base));
310 (*_status)[node] = EVEN;
311 _node_queue[_last++] = node;
312 arc = _graph.oppositeArc((*_ear)[node]);
313 node = _graph.target((*_ear)[node]);
314 base = (*_blossom_rep)[_blossom_set->find(node)];
315 _blossom_set->join(_graph.target(arc), base);
319 (*_blossom_rep)[_blossom_set->find(nca)] = nca;
323 Node node = _graph.v(e);
324 Arc arc = _graph.direct(e, false);
325 Node base = (*_blossom_rep)[_blossom_set->find(node)];
327 while (base != nca) {
332 n = _graph.target((*_matching)[n]);
334 n = _graph.target(a);
335 (*_ear)[n] = _graph.oppositeArc(a);
337 node = _graph.target((*_matching)[base]);
338 _tree_set->erase(base);
339 _tree_set->erase(node);
340 _blossom_set->insert(node, _blossom_set->find(base));
341 (*_status)[node] = EVEN;
342 _node_queue[_last++] = node;
343 arc = _graph.oppositeArc((*_ear)[node]);
344 node = _graph.target((*_ear)[node]);
345 base = (*_blossom_rep)[_blossom_set->find(node)];
346 _blossom_set->join(_graph.target(arc), base);
350 (*_blossom_rep)[_blossom_set->find(nca)] = nca;
353 void extendOnArc(const Arc& a) {
354 Node base = _graph.source(a);
355 Node odd = _graph.target(a);
357 (*_ear)[odd] = _graph.oppositeArc(a);
358 Node even = _graph.target((*_matching)[odd]);
359 (*_blossom_rep)[_blossom_set->insert(even)] = even;
360 (*_status)[odd] = ODD;
361 (*_status)[even] = EVEN;
362 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
363 _tree_set->insert(odd, tree);
364 _tree_set->insert(even, tree);
365 _node_queue[_last++] = even;
369 void augmentOnArc(const Arc& a) {
370 Node even = _graph.source(a);
371 Node odd = _graph.target(a);
373 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
375 (*_matching)[odd] = _graph.oppositeArc(a);
376 (*_status)[odd] = MATCHED;
378 Arc arc = (*_matching)[even];
379 (*_matching)[even] = a;
381 while (arc != INVALID) {
382 odd = _graph.target(arc);
384 even = _graph.target(arc);
385 (*_matching)[odd] = arc;
386 arc = (*_matching)[even];
387 (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
390 for (typename TreeSet::ItemIt it(*_tree_set, tree);
391 it != INVALID; ++it) {
392 if ((*_status)[it] == ODD) {
393 (*_status)[it] = MATCHED;
395 int blossom = _blossom_set->find(it);
396 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
397 jt != INVALID; ++jt) {
398 (*_status)[jt] = MATCHED;
400 _blossom_set->eraseClass(blossom);
403 _tree_set->eraseClass(tree);
409 /// \brief Constructor
412 MaxMatching(const Graph& graph)
413 : _graph(graph), _matching(0), _status(0), _ear(0),
414 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
415 _tree_set_index(0), _tree_set(0) {}
421 /// \name Execution Control
422 /// The simplest way to execute the algorithm is to use the
423 /// \c run() member function.\n
424 /// If you need better control on the execution, you have to call
425 /// one of the functions \ref init(), \ref greedyInit() or
426 /// \ref matchingInit() first, then you can start the algorithm with
427 /// \ref startSparse() or \ref startDense().
431 /// \brief Set the initial matching to the empty matching.
433 /// This function sets the initial matching to the empty matching.
436 for(NodeIt n(_graph); n != INVALID; ++n) {
437 (*_matching)[n] = INVALID;
438 (*_status)[n] = UNMATCHED;
442 /// \brief Find an initial matching in a greedy way.
444 /// This function finds an initial matching in a greedy way.
447 for (NodeIt n(_graph); n != INVALID; ++n) {
448 (*_matching)[n] = INVALID;
449 (*_status)[n] = UNMATCHED;
451 for (NodeIt n(_graph); n != INVALID; ++n) {
452 if ((*_matching)[n] == INVALID) {
453 for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
454 Node v = _graph.target(a);
455 if ((*_matching)[v] == INVALID && v != n) {
457 (*_status)[n] = MATCHED;
458 (*_matching)[v] = _graph.oppositeArc(a);
459 (*_status)[v] = MATCHED;
468 /// \brief Initialize the matching from a map.
470 /// This function initializes the matching from a \c bool valued edge
471 /// map. This map should have the property that there are no two incident
472 /// edges with \c true value, i.e. it really contains a matching.
473 /// \return \c true if the map contains a matching.
474 template <typename MatchingMap>
475 bool matchingInit(const MatchingMap& matching) {
478 for (NodeIt n(_graph); n != INVALID; ++n) {
479 (*_matching)[n] = INVALID;
480 (*_status)[n] = UNMATCHED;
482 for(EdgeIt e(_graph); e!=INVALID; ++e) {
485 Node u = _graph.u(e);
486 if ((*_matching)[u] != INVALID) return false;
487 (*_matching)[u] = _graph.direct(e, true);
488 (*_status)[u] = MATCHED;
490 Node v = _graph.v(e);
491 if ((*_matching)[v] != INVALID) return false;
492 (*_matching)[v] = _graph.direct(e, false);
493 (*_status)[v] = MATCHED;
499 /// \brief Start Edmonds' algorithm
501 /// This function runs the original Edmonds' algorithm.
503 /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
504 /// called before using this function.
506 for(NodeIt n(_graph); n != INVALID; ++n) {
507 if ((*_status)[n] == UNMATCHED) {
508 (*_blossom_rep)[_blossom_set->insert(n)] = n;
509 _tree_set->insert(n);
510 (*_status)[n] = EVEN;
516 /// \brief Start Edmonds' algorithm with a heuristic improvement
519 /// This function runs Edmonds' algorithm with a heuristic of postponing
520 /// shrinks, therefore resulting in a faster algorithm for dense graphs.
522 /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
523 /// called before using this function.
525 for(NodeIt n(_graph); n != INVALID; ++n) {
526 if ((*_status)[n] == UNMATCHED) {
527 (*_blossom_rep)[_blossom_set->insert(n)] = n;
528 _tree_set->insert(n);
529 (*_status)[n] = EVEN;
536 /// \brief Run Edmonds' algorithm
538 /// This function runs Edmonds' algorithm. An additional heuristic of
539 /// postponing shrinks is used for relatively dense graphs
540 /// (for which <tt>m>=2*n</tt> holds).
542 if (countEdges(_graph) < 2 * countNodes(_graph)) {
553 /// \name Primal Solution
554 /// Functions to get the primal solution, i.e. the maximum matching.
558 /// \brief Return the size (cardinality) of the matching.
560 /// This function returns the size (cardinality) of the current matching.
561 /// After run() it returns the size of the maximum matching in the graph.
562 int matchingSize() const {
564 for (NodeIt n(_graph); n != INVALID; ++n) {
565 if ((*_matching)[n] != INVALID) {
572 /// \brief Return \c true if the given edge is in the matching.
574 /// This function returns \c true if the given edge is in the current
576 bool matching(const Edge& edge) const {
577 return edge == (*_matching)[_graph.u(edge)];
580 /// \brief Return the matching arc (or edge) incident to the given node.
582 /// This function returns the matching arc (or edge) incident to the
583 /// given node in the current matching or \c INVALID if the node is
584 /// not covered by the matching.
585 Arc matching(const Node& n) const {
586 return (*_matching)[n];
589 /// \brief Return a const reference to the matching map.
591 /// This function returns a const reference to a node map that stores
592 /// the matching arc (or edge) incident to each node.
593 const MatchingMap& matchingMap() const {
597 /// \brief Return the mate of the given node.
599 /// This function returns the mate of the given node in the current
600 /// matching or \c INVALID if the node is not covered by the matching.
601 Node mate(const Node& n) const {
602 return (*_matching)[n] != INVALID ?
603 _graph.target((*_matching)[n]) : INVALID;
608 /// \name Dual Solution
609 /// Functions to get the dual solution, i.e. the Gallai-Edmonds
614 /// \brief Return the status of the given node in the Edmonds-Gallai
617 /// This function returns the \ref Status "status" of the given node
618 /// in the Edmonds-Gallai decomposition.
619 Status status(const Node& n) const {
620 return (*_status)[n];
623 /// \brief Return a const reference to the status map, which stores
624 /// the Edmonds-Gallai decomposition.
626 /// This function returns a const reference to a node map that stores the
627 /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
628 const StatusMap& statusMap() const {
632 /// \brief Return \c true if the given node is in the barrier.
634 /// This function returns \c true if the given node is in the barrier.
635 bool barrier(const Node& n) const {
636 return (*_status)[n] == ODD;
643 /// \ingroup matching
645 /// \brief Weighted matching in general graphs
647 /// This class provides an efficient implementation of Edmond's
648 /// maximum weighted matching algorithm. The implementation is based
649 /// on extensive use of priority queues and provides
650 /// \f$O(nm\log n)\f$ time complexity.
652 /// The maximum weighted matching problem is to find a subset of the
653 /// edges in an undirected graph with maximum overall weight for which
654 /// each node has at most one incident edge.
655 /// It can be formulated with the following linear program.
656 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
657 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
658 \quad \forall B\in\mathcal{O}\f] */
659 /// \f[x_e \ge 0\quad \forall e\in E\f]
660 /// \f[\max \sum_{e\in E}x_ew_e\f]
661 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
662 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
663 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
664 /// subsets of the nodes.
666 /// The algorithm calculates an optimal matching and a proof of the
667 /// optimality. The solution of the dual problem can be used to check
668 /// the result of the algorithm. The dual linear problem is the
670 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
671 z_B \ge w_{uv} \quad \forall uv\in E\f] */
672 /// \f[y_u \ge 0 \quad \forall u \in V\f]
673 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
674 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
675 \frac{\vert B \vert - 1}{2}z_B\f] */
677 /// The algorithm can be executed with the run() function.
678 /// After it the matching (the primal solution) and the dual solution
679 /// can be obtained using the query functions and the
680 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
681 /// which is able to iterate on the nodes of a blossom.
682 /// If the value type is integer, then the dual solution is multiplied
683 /// by \ref MaxWeightedMatching::dualScale "4".
685 /// \tparam GR The undirected graph type the algorithm runs on.
686 /// \tparam WM The type edge weight map. The default type is
687 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
689 template <typename GR, typename WM>
691 template <typename GR,
692 typename WM = typename GR::template EdgeMap<int> >
694 class MaxWeightedMatching {
697 /// The graph type of the algorithm
699 /// The type of the edge weight map
700 typedef WM WeightMap;
701 /// The value type of the edge weights
702 typedef typename WeightMap::Value Value;
704 /// The type of the matching map
705 typedef typename Graph::template NodeMap<typename Graph::Arc>
708 /// \brief Scaling factor for dual solution
710 /// Scaling factor for dual solution. It is equal to 4 or 1
711 /// according to the value type.
712 static const int dualScale =
713 std::numeric_limits<Value>::is_integer ? 4 : 1;
717 TEMPLATE_GRAPH_TYPEDEFS(Graph);
719 typedef typename Graph::template NodeMap<Value> NodePotential;
720 typedef std::vector<Node> BlossomNodeList;
722 struct BlossomVariable {
726 BlossomVariable(int _begin, int _end, Value _value)
727 : begin(_begin), end(_end), value(_value) {}
731 typedef std::vector<BlossomVariable> BlossomPotential;
734 const WeightMap& _weight;
736 MatchingMap* _matching;
738 NodePotential* _node_potential;
740 BlossomPotential _blossom_potential;
741 BlossomNodeList _blossom_node_list;
746 typedef RangeMap<int> IntIntMap;
749 EVEN = -1, MATCHED = 0, ODD = 1
752 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
761 IntNodeMap *_blossom_index;
762 BlossomSet *_blossom_set;
763 RangeMap<BlossomData>* _blossom_data;
765 IntNodeMap *_node_index;
766 IntArcMap *_node_heap_index;
770 NodeData(IntArcMap& node_heap_index)
771 : heap(node_heap_index) {}
775 BinHeap<Value, IntArcMap> heap;
776 std::map<int, Arc> heap_index;
781 RangeMap<NodeData>* _node_data;
783 typedef ExtendFindEnum<IntIntMap> TreeSet;
785 IntIntMap *_tree_set_index;
788 IntNodeMap *_delta1_index;
789 BinHeap<Value, IntNodeMap> *_delta1;
791 IntIntMap *_delta2_index;
792 BinHeap<Value, IntIntMap> *_delta2;
794 IntEdgeMap *_delta3_index;
795 BinHeap<Value, IntEdgeMap> *_delta3;
797 IntIntMap *_delta4_index;
798 BinHeap<Value, IntIntMap> *_delta4;
803 typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
804 FractionalMatching *_fractional;
806 void createStructures() {
807 _node_num = countNodes(_graph);
808 _blossom_num = _node_num * 3 / 2;
811 _matching = new MatchingMap(_graph);
814 if (!_node_potential) {
815 _node_potential = new NodePotential(_graph);
819 _blossom_index = new IntNodeMap(_graph);
820 _blossom_set = new BlossomSet(*_blossom_index);
821 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
822 } else if (_blossom_data->size() != _blossom_num) {
823 delete _blossom_data;
824 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
828 _node_index = new IntNodeMap(_graph);
829 _node_heap_index = new IntArcMap(_graph);
830 _node_data = new RangeMap<NodeData>(_node_num,
831 NodeData(*_node_heap_index));
834 _node_data = new RangeMap<NodeData>(_node_num,
835 NodeData(*_node_heap_index));
839 _tree_set_index = new IntIntMap(_blossom_num);
840 _tree_set = new TreeSet(*_tree_set_index);
842 _tree_set_index->resize(_blossom_num);
846 _delta1_index = new IntNodeMap(_graph);
847 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
851 _delta2_index = new IntIntMap(_blossom_num);
852 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
854 _delta2_index->resize(_blossom_num);
858 _delta3_index = new IntEdgeMap(_graph);
859 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
863 _delta4_index = new IntIntMap(_blossom_num);
864 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
866 _delta4_index->resize(_blossom_num);
870 void destroyStructures() {
874 if (_node_potential) {
875 delete _node_potential;
878 delete _blossom_index;
880 delete _blossom_data;
885 delete _node_heap_index;
890 delete _tree_set_index;
894 delete _delta1_index;
898 delete _delta2_index;
902 delete _delta3_index;
906 delete _delta4_index;
911 void matchedToEven(int blossom, int tree) {
912 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
913 _delta2->erase(blossom);
916 if (!_blossom_set->trivial(blossom)) {
917 (*_blossom_data)[blossom].pot -=
918 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
921 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
924 _blossom_set->increase(n, std::numeric_limits<Value>::max());
925 int ni = (*_node_index)[n];
927 (*_node_data)[ni].heap.clear();
928 (*_node_data)[ni].heap_index.clear();
930 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
932 _delta1->push(n, (*_node_data)[ni].pot);
934 for (InArcIt e(_graph, n); e != INVALID; ++e) {
935 Node v = _graph.source(e);
936 int vb = _blossom_set->find(v);
937 int vi = (*_node_index)[v];
939 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
940 dualScale * _weight[e];
942 if ((*_blossom_data)[vb].status == EVEN) {
943 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
944 _delta3->push(e, rw / 2);
947 typename std::map<int, Arc>::iterator it =
948 (*_node_data)[vi].heap_index.find(tree);
950 if (it != (*_node_data)[vi].heap_index.end()) {
951 if ((*_node_data)[vi].heap[it->second] > rw) {
952 (*_node_data)[vi].heap.replace(it->second, e);
953 (*_node_data)[vi].heap.decrease(e, rw);
957 (*_node_data)[vi].heap.push(e, rw);
958 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
961 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
962 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
964 if ((*_blossom_data)[vb].status == MATCHED) {
965 if (_delta2->state(vb) != _delta2->IN_HEAP) {
966 _delta2->push(vb, _blossom_set->classPrio(vb) -
967 (*_blossom_data)[vb].offset);
968 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
969 (*_blossom_data)[vb].offset) {
970 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
971 (*_blossom_data)[vb].offset);
978 (*_blossom_data)[blossom].offset = 0;
981 void matchedToOdd(int blossom) {
982 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
983 _delta2->erase(blossom);
985 (*_blossom_data)[blossom].offset += _delta_sum;
986 if (!_blossom_set->trivial(blossom)) {
987 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
988 (*_blossom_data)[blossom].offset);
992 void evenToMatched(int blossom, int tree) {
993 if (!_blossom_set->trivial(blossom)) {
994 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
997 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
999 int ni = (*_node_index)[n];
1000 (*_node_data)[ni].pot -= _delta_sum;
1004 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1005 Node v = _graph.source(e);
1006 int vb = _blossom_set->find(v);
1007 int vi = (*_node_index)[v];
1009 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1010 dualScale * _weight[e];
1012 if (vb == blossom) {
1013 if (_delta3->state(e) == _delta3->IN_HEAP) {
1016 } else if ((*_blossom_data)[vb].status == EVEN) {
1018 if (_delta3->state(e) == _delta3->IN_HEAP) {
1022 int vt = _tree_set->find(vb);
1026 Arc r = _graph.oppositeArc(e);
1028 typename std::map<int, Arc>::iterator it =
1029 (*_node_data)[ni].heap_index.find(vt);
1031 if (it != (*_node_data)[ni].heap_index.end()) {
1032 if ((*_node_data)[ni].heap[it->second] > rw) {
1033 (*_node_data)[ni].heap.replace(it->second, r);
1034 (*_node_data)[ni].heap.decrease(r, rw);
1038 (*_node_data)[ni].heap.push(r, rw);
1039 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1042 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1043 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1045 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1046 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1047 (*_blossom_data)[blossom].offset);
1048 } else if ((*_delta2)[blossom] >
1049 _blossom_set->classPrio(blossom) -
1050 (*_blossom_data)[blossom].offset){
1051 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1052 (*_blossom_data)[blossom].offset);
1058 typename std::map<int, Arc>::iterator it =
1059 (*_node_data)[vi].heap_index.find(tree);
1061 if (it != (*_node_data)[vi].heap_index.end()) {
1062 (*_node_data)[vi].heap.erase(it->second);
1063 (*_node_data)[vi].heap_index.erase(it);
1064 if ((*_node_data)[vi].heap.empty()) {
1065 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1066 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1067 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1070 if ((*_blossom_data)[vb].status == MATCHED) {
1071 if (_blossom_set->classPrio(vb) ==
1072 std::numeric_limits<Value>::max()) {
1074 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1075 (*_blossom_data)[vb].offset) {
1076 _delta2->increase(vb, _blossom_set->classPrio(vb) -
1077 (*_blossom_data)[vb].offset);
1086 void oddToMatched(int blossom) {
1087 (*_blossom_data)[blossom].offset -= _delta_sum;
1089 if (_blossom_set->classPrio(blossom) !=
1090 std::numeric_limits<Value>::max()) {
1091 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1092 (*_blossom_data)[blossom].offset);
1095 if (!_blossom_set->trivial(blossom)) {
1096 _delta4->erase(blossom);
1100 void oddToEven(int blossom, int tree) {
1101 if (!_blossom_set->trivial(blossom)) {
1102 _delta4->erase(blossom);
1103 (*_blossom_data)[blossom].pot -=
1104 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1107 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1108 n != INVALID; ++n) {
1109 int ni = (*_node_index)[n];
1111 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1113 (*_node_data)[ni].heap.clear();
1114 (*_node_data)[ni].heap_index.clear();
1115 (*_node_data)[ni].pot +=
1116 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1118 _delta1->push(n, (*_node_data)[ni].pot);
1120 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1121 Node v = _graph.source(e);
1122 int vb = _blossom_set->find(v);
1123 int vi = (*_node_index)[v];
1125 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1126 dualScale * _weight[e];
1128 if ((*_blossom_data)[vb].status == EVEN) {
1129 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1130 _delta3->push(e, rw / 2);
1134 typename std::map<int, Arc>::iterator it =
1135 (*_node_data)[vi].heap_index.find(tree);
1137 if (it != (*_node_data)[vi].heap_index.end()) {
1138 if ((*_node_data)[vi].heap[it->second] > rw) {
1139 (*_node_data)[vi].heap.replace(it->second, e);
1140 (*_node_data)[vi].heap.decrease(e, rw);
1144 (*_node_data)[vi].heap.push(e, rw);
1145 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1148 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1149 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1151 if ((*_blossom_data)[vb].status == MATCHED) {
1152 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1153 _delta2->push(vb, _blossom_set->classPrio(vb) -
1154 (*_blossom_data)[vb].offset);
1155 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1156 (*_blossom_data)[vb].offset) {
1157 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1158 (*_blossom_data)[vb].offset);
1165 (*_blossom_data)[blossom].offset = 0;
1168 void alternatePath(int even, int tree) {
1171 evenToMatched(even, tree);
1172 (*_blossom_data)[even].status = MATCHED;
1174 while ((*_blossom_data)[even].pred != INVALID) {
1175 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1176 (*_blossom_data)[odd].status = MATCHED;
1178 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1180 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1181 (*_blossom_data)[even].status = MATCHED;
1182 evenToMatched(even, tree);
1183 (*_blossom_data)[even].next =
1184 _graph.oppositeArc((*_blossom_data)[odd].pred);
1189 void destroyTree(int tree) {
1190 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1191 if ((*_blossom_data)[b].status == EVEN) {
1192 (*_blossom_data)[b].status = MATCHED;
1193 evenToMatched(b, tree);
1194 } else if ((*_blossom_data)[b].status == ODD) {
1195 (*_blossom_data)[b].status = MATCHED;
1199 _tree_set->eraseClass(tree);
1203 void unmatchNode(const Node& node) {
1204 int blossom = _blossom_set->find(node);
1205 int tree = _tree_set->find(blossom);
1207 alternatePath(blossom, tree);
1210 (*_blossom_data)[blossom].base = node;
1211 (*_blossom_data)[blossom].next = INVALID;
1214 void augmentOnEdge(const Edge& edge) {
1216 int left = _blossom_set->find(_graph.u(edge));
1217 int right = _blossom_set->find(_graph.v(edge));
1219 int left_tree = _tree_set->find(left);
1220 alternatePath(left, left_tree);
1221 destroyTree(left_tree);
1223 int right_tree = _tree_set->find(right);
1224 alternatePath(right, right_tree);
1225 destroyTree(right_tree);
1227 (*_blossom_data)[left].next = _graph.direct(edge, true);
1228 (*_blossom_data)[right].next = _graph.direct(edge, false);
1231 void augmentOnArc(const Arc& arc) {
1233 int left = _blossom_set->find(_graph.source(arc));
1234 int right = _blossom_set->find(_graph.target(arc));
1236 (*_blossom_data)[left].status = MATCHED;
1238 int right_tree = _tree_set->find(right);
1239 alternatePath(right, right_tree);
1240 destroyTree(right_tree);
1242 (*_blossom_data)[left].next = arc;
1243 (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1246 void extendOnArc(const Arc& arc) {
1247 int base = _blossom_set->find(_graph.target(arc));
1248 int tree = _tree_set->find(base);
1250 int odd = _blossom_set->find(_graph.source(arc));
1251 _tree_set->insert(odd, tree);
1252 (*_blossom_data)[odd].status = ODD;
1254 (*_blossom_data)[odd].pred = arc;
1256 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1257 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1258 _tree_set->insert(even, tree);
1259 (*_blossom_data)[even].status = EVEN;
1260 matchedToEven(even, tree);
1263 void shrinkOnEdge(const Edge& edge, int tree) {
1265 std::vector<int> left_path, right_path;
1268 std::set<int> left_set, right_set;
1269 int left = _blossom_set->find(_graph.u(edge));
1270 left_path.push_back(left);
1271 left_set.insert(left);
1273 int right = _blossom_set->find(_graph.v(edge));
1274 right_path.push_back(right);
1275 right_set.insert(right);
1279 if ((*_blossom_data)[left].pred == INVALID) break;
1282 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1283 left_path.push_back(left);
1285 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1286 left_path.push_back(left);
1288 left_set.insert(left);
1290 if (right_set.find(left) != right_set.end()) {
1295 if ((*_blossom_data)[right].pred == INVALID) break;
1298 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1299 right_path.push_back(right);
1301 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1302 right_path.push_back(right);
1304 right_set.insert(right);
1306 if (left_set.find(right) != left_set.end()) {
1314 if ((*_blossom_data)[left].pred == INVALID) {
1316 while (left_set.find(nca) == left_set.end()) {
1318 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1319 right_path.push_back(nca);
1321 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1322 right_path.push_back(nca);
1326 while (right_set.find(nca) == right_set.end()) {
1328 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1329 left_path.push_back(nca);
1331 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1332 left_path.push_back(nca);
1338 std::vector<int> subblossoms;
1341 prev = _graph.direct(edge, true);
1342 for (int i = 0; left_path[i] != nca; i += 2) {
1343 subblossoms.push_back(left_path[i]);
1344 (*_blossom_data)[left_path[i]].next = prev;
1345 _tree_set->erase(left_path[i]);
1347 subblossoms.push_back(left_path[i + 1]);
1348 (*_blossom_data)[left_path[i + 1]].status = EVEN;
1349 oddToEven(left_path[i + 1], tree);
1350 _tree_set->erase(left_path[i + 1]);
1351 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1355 while (right_path[k] != nca) ++k;
1357 subblossoms.push_back(nca);
1358 (*_blossom_data)[nca].next = prev;
1360 for (int i = k - 2; i >= 0; i -= 2) {
1361 subblossoms.push_back(right_path[i + 1]);
1362 (*_blossom_data)[right_path[i + 1]].status = EVEN;
1363 oddToEven(right_path[i + 1], tree);
1364 _tree_set->erase(right_path[i + 1]);
1366 (*_blossom_data)[right_path[i + 1]].next =
1367 (*_blossom_data)[right_path[i + 1]].pred;
1369 subblossoms.push_back(right_path[i]);
1370 _tree_set->erase(right_path[i]);
1374 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1376 for (int i = 0; i < int(subblossoms.size()); ++i) {
1377 if (!_blossom_set->trivial(subblossoms[i])) {
1378 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1380 (*_blossom_data)[subblossoms[i]].status = MATCHED;
1383 (*_blossom_data)[surface].pot = -2 * _delta_sum;
1384 (*_blossom_data)[surface].offset = 0;
1385 (*_blossom_data)[surface].status = EVEN;
1386 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1387 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1389 _tree_set->insert(surface, tree);
1390 _tree_set->erase(nca);
1393 void splitBlossom(int blossom) {
1394 Arc next = (*_blossom_data)[blossom].next;
1395 Arc pred = (*_blossom_data)[blossom].pred;
1397 int tree = _tree_set->find(blossom);
1399 (*_blossom_data)[blossom].status = MATCHED;
1400 oddToMatched(blossom);
1401 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1402 _delta2->erase(blossom);
1405 std::vector<int> subblossoms;
1406 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1408 Value offset = (*_blossom_data)[blossom].offset;
1409 int b = _blossom_set->find(_graph.source(pred));
1410 int d = _blossom_set->find(_graph.source(next));
1412 int ib = -1, id = -1;
1413 for (int i = 0; i < int(subblossoms.size()); ++i) {
1414 if (subblossoms[i] == b) ib = i;
1415 if (subblossoms[i] == d) id = i;
1417 (*_blossom_data)[subblossoms[i]].offset = offset;
1418 if (!_blossom_set->trivial(subblossoms[i])) {
1419 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1421 if (_blossom_set->classPrio(subblossoms[i]) !=
1422 std::numeric_limits<Value>::max()) {
1423 _delta2->push(subblossoms[i],
1424 _blossom_set->classPrio(subblossoms[i]) -
1425 (*_blossom_data)[subblossoms[i]].offset);
1429 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1430 for (int i = (id + 1) % subblossoms.size();
1431 i != ib; i = (i + 2) % subblossoms.size()) {
1432 int sb = subblossoms[i];
1433 int tb = subblossoms[(i + 1) % subblossoms.size()];
1434 (*_blossom_data)[sb].next =
1435 _graph.oppositeArc((*_blossom_data)[tb].next);
1438 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1439 int sb = subblossoms[i];
1440 int tb = subblossoms[(i + 1) % subblossoms.size()];
1441 int ub = subblossoms[(i + 2) % subblossoms.size()];
1443 (*_blossom_data)[sb].status = ODD;
1445 _tree_set->insert(sb, tree);
1446 (*_blossom_data)[sb].pred = pred;
1447 (*_blossom_data)[sb].next =
1448 _graph.oppositeArc((*_blossom_data)[tb].next);
1450 pred = (*_blossom_data)[ub].next;
1452 (*_blossom_data)[tb].status = EVEN;
1453 matchedToEven(tb, tree);
1454 _tree_set->insert(tb, tree);
1455 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1458 (*_blossom_data)[subblossoms[id]].status = ODD;
1459 matchedToOdd(subblossoms[id]);
1460 _tree_set->insert(subblossoms[id], tree);
1461 (*_blossom_data)[subblossoms[id]].next = next;
1462 (*_blossom_data)[subblossoms[id]].pred = pred;
1466 for (int i = (ib + 1) % subblossoms.size();
1467 i != id; i = (i + 2) % subblossoms.size()) {
1468 int sb = subblossoms[i];
1469 int tb = subblossoms[(i + 1) % subblossoms.size()];
1470 (*_blossom_data)[sb].next =
1471 _graph.oppositeArc((*_blossom_data)[tb].next);
1474 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1475 int sb = subblossoms[i];
1476 int tb = subblossoms[(i + 1) % subblossoms.size()];
1477 int ub = subblossoms[(i + 2) % subblossoms.size()];
1479 (*_blossom_data)[sb].status = ODD;
1481 _tree_set->insert(sb, tree);
1482 (*_blossom_data)[sb].next = next;
1483 (*_blossom_data)[sb].pred =
1484 _graph.oppositeArc((*_blossom_data)[tb].next);
1486 (*_blossom_data)[tb].status = EVEN;
1487 matchedToEven(tb, tree);
1488 _tree_set->insert(tb, tree);
1489 (*_blossom_data)[tb].pred =
1490 (*_blossom_data)[tb].next =
1491 _graph.oppositeArc((*_blossom_data)[ub].next);
1492 next = (*_blossom_data)[ub].next;
1495 (*_blossom_data)[subblossoms[ib]].status = ODD;
1496 matchedToOdd(subblossoms[ib]);
1497 _tree_set->insert(subblossoms[ib], tree);
1498 (*_blossom_data)[subblossoms[ib]].next = next;
1499 (*_blossom_data)[subblossoms[ib]].pred = pred;
1501 _tree_set->erase(blossom);
1504 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1505 if (_blossom_set->trivial(blossom)) {
1506 int bi = (*_node_index)[base];
1507 Value pot = (*_node_data)[bi].pot;
1509 (*_matching)[base] = matching;
1510 _blossom_node_list.push_back(base);
1511 (*_node_potential)[base] = pot;
1514 Value pot = (*_blossom_data)[blossom].pot;
1515 int bn = _blossom_node_list.size();
1517 std::vector<int> subblossoms;
1518 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1519 int b = _blossom_set->find(base);
1521 for (int i = 0; i < int(subblossoms.size()); ++i) {
1522 if (subblossoms[i] == b) { ib = i; break; }
1525 for (int i = 1; i < int(subblossoms.size()); i += 2) {
1526 int sb = subblossoms[(ib + i) % subblossoms.size()];
1527 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1529 Arc m = (*_blossom_data)[tb].next;
1530 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1531 extractBlossom(tb, _graph.source(m), m);
1533 extractBlossom(subblossoms[ib], base, matching);
1535 int en = _blossom_node_list.size();
1537 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1541 void extractMatching() {
1542 std::vector<int> blossoms;
1543 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1544 blossoms.push_back(c);
1547 for (int i = 0; i < int(blossoms.size()); ++i) {
1548 if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1550 Value offset = (*_blossom_data)[blossoms[i]].offset;
1551 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1552 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1553 n != INVALID; ++n) {
1554 (*_node_data)[(*_node_index)[n]].pot -= offset;
1557 Arc matching = (*_blossom_data)[blossoms[i]].next;
1558 Node base = _graph.source(matching);
1559 extractBlossom(blossoms[i], base, matching);
1561 Node base = (*_blossom_data)[blossoms[i]].base;
1562 extractBlossom(blossoms[i], base, INVALID);
1569 /// \brief Constructor
1572 MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1573 : _graph(graph), _weight(weight), _matching(0),
1574 _node_potential(0), _blossom_potential(), _blossom_node_list(),
1575 _node_num(0), _blossom_num(0),
1577 _blossom_index(0), _blossom_set(0), _blossom_data(0),
1578 _node_index(0), _node_heap_index(0), _node_data(0),
1579 _tree_set_index(0), _tree_set(0),
1581 _delta1_index(0), _delta1(0),
1582 _delta2_index(0), _delta2(0),
1583 _delta3_index(0), _delta3(0),
1584 _delta4_index(0), _delta4(0),
1586 _delta_sum(), _unmatched(0),
1591 ~MaxWeightedMatching() {
1592 destroyStructures();
1598 /// \name Execution Control
1599 /// The simplest way to execute the algorithm is to use the
1600 /// \ref run() member function.
1604 /// \brief Initialize the algorithm
1606 /// This function initializes the algorithm.
1610 _blossom_node_list.clear();
1611 _blossom_potential.clear();
1613 for (ArcIt e(_graph); e != INVALID; ++e) {
1614 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1616 for (NodeIt n(_graph); n != INVALID; ++n) {
1617 (*_delta1_index)[n] = _delta1->PRE_HEAP;
1619 for (EdgeIt e(_graph); e != INVALID; ++e) {
1620 (*_delta3_index)[e] = _delta3->PRE_HEAP;
1622 for (int i = 0; i < _blossom_num; ++i) {
1623 (*_delta2_index)[i] = _delta2->PRE_HEAP;
1624 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1627 _unmatched = _node_num;
1633 _blossom_set->clear();
1637 for (NodeIt n(_graph); n != INVALID; ++n) {
1639 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1640 if (_graph.target(e) == n) continue;
1641 if ((dualScale * _weight[e]) / 2 > max) {
1642 max = (dualScale * _weight[e]) / 2;
1645 (*_node_index)[n] = index;
1646 (*_node_data)[index].heap_index.clear();
1647 (*_node_data)[index].heap.clear();
1648 (*_node_data)[index].pot = max;
1649 _delta1->push(n, max);
1651 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1653 _tree_set->insert(blossom);
1655 (*_blossom_data)[blossom].status = EVEN;
1656 (*_blossom_data)[blossom].pred = INVALID;
1657 (*_blossom_data)[blossom].next = INVALID;
1658 (*_blossom_data)[blossom].pot = 0;
1659 (*_blossom_data)[blossom].offset = 0;
1662 for (EdgeIt e(_graph); e != INVALID; ++e) {
1663 int si = (*_node_index)[_graph.u(e)];
1664 int ti = (*_node_index)[_graph.v(e)];
1665 if (_graph.u(e) != _graph.v(e)) {
1666 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1667 dualScale * _weight[e]) / 2);
1672 /// \brief Initialize the algorithm with fractional matching
1674 /// This function initializes the algorithm with a fractional
1675 /// matching. This initialization is also called jumpstart heuristic.
1676 void fractionalInit() {
1679 if (_fractional == 0) {
1680 _fractional = new FractionalMatching(_graph, _weight, false);
1684 for (ArcIt e(_graph); e != INVALID; ++e) {
1685 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1687 for (NodeIt n(_graph); n != INVALID; ++n) {
1688 (*_delta1_index)[n] = _delta1->PRE_HEAP;
1690 for (EdgeIt e(_graph); e != INVALID; ++e) {
1691 (*_delta3_index)[e] = _delta3->PRE_HEAP;
1693 for (int i = 0; i < _blossom_num; ++i) {
1694 (*_delta2_index)[i] = _delta2->PRE_HEAP;
1695 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1701 for (NodeIt n(_graph); n != INVALID; ++n) {
1702 Value pot = _fractional->nodeValue(n);
1703 (*_node_index)[n] = index;
1704 (*_node_data)[index].pot = pot;
1706 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1708 (*_blossom_data)[blossom].status = MATCHED;
1709 (*_blossom_data)[blossom].pred = INVALID;
1710 (*_blossom_data)[blossom].next = _fractional->matching(n);
1711 if (_fractional->matching(n) == INVALID) {
1712 (*_blossom_data)[blossom].base = n;
1714 (*_blossom_data)[blossom].pot = 0;
1715 (*_blossom_data)[blossom].offset = 0;
1719 typename Graph::template NodeMap<bool> processed(_graph, false);
1720 for (NodeIt n(_graph); n != INVALID; ++n) {
1721 if (processed[n]) continue;
1722 processed[n] = true;
1723 if (_fractional->matching(n) == INVALID) continue;
1725 Node v = _graph.target(_fractional->matching(n));
1727 processed[v] = true;
1728 v = _graph.target(_fractional->matching(v));
1733 std::vector<int> subblossoms(num);
1735 subblossoms[--num] = _blossom_set->find(n);
1736 _delta1->push(n, _fractional->nodeValue(n));
1737 v = _graph.target(_fractional->matching(n));
1739 subblossoms[--num] = _blossom_set->find(v);
1740 _delta1->push(v, _fractional->nodeValue(v));
1741 v = _graph.target(_fractional->matching(v));
1745 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1746 (*_blossom_data)[surface].status = EVEN;
1747 (*_blossom_data)[surface].pred = INVALID;
1748 (*_blossom_data)[surface].next = INVALID;
1749 (*_blossom_data)[surface].pot = 0;
1750 (*_blossom_data)[surface].offset = 0;
1752 _tree_set->insert(surface);
1757 for (EdgeIt e(_graph); e != INVALID; ++e) {
1758 int si = (*_node_index)[_graph.u(e)];
1759 int sb = _blossom_set->find(_graph.u(e));
1760 int ti = (*_node_index)[_graph.v(e)];
1761 int tb = _blossom_set->find(_graph.v(e));
1762 if ((*_blossom_data)[sb].status == EVEN &&
1763 (*_blossom_data)[tb].status == EVEN && sb != tb) {
1764 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1765 dualScale * _weight[e]) / 2);
1769 for (NodeIt n(_graph); n != INVALID; ++n) {
1770 int nb = _blossom_set->find(n);
1771 if ((*_blossom_data)[nb].status != MATCHED) continue;
1772 int ni = (*_node_index)[n];
1774 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1775 Node v = _graph.target(e);
1776 int vb = _blossom_set->find(v);
1777 int vi = (*_node_index)[v];
1779 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1780 dualScale * _weight[e];
1782 if ((*_blossom_data)[vb].status == EVEN) {
1784 int vt = _tree_set->find(vb);
1786 typename std::map<int, Arc>::iterator it =
1787 (*_node_data)[ni].heap_index.find(vt);
1789 if (it != (*_node_data)[ni].heap_index.end()) {
1790 if ((*_node_data)[ni].heap[it->second] > rw) {
1791 (*_node_data)[ni].heap.replace(it->second, e);
1792 (*_node_data)[ni].heap.decrease(e, rw);
1796 (*_node_data)[ni].heap.push(e, rw);
1797 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1802 if (!(*_node_data)[ni].heap.empty()) {
1803 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1804 _delta2->push(nb, _blossom_set->classPrio(nb));
1809 /// \brief Start the algorithm
1811 /// This function starts the algorithm.
1813 /// \pre \ref init() or \ref fractionalInit() must be called
1814 /// before using this function.
1820 while (_unmatched > 0) {
1821 Value d1 = !_delta1->empty() ?
1822 _delta1->prio() : std::numeric_limits<Value>::max();
1824 Value d2 = !_delta2->empty() ?
1825 _delta2->prio() : std::numeric_limits<Value>::max();
1827 Value d3 = !_delta3->empty() ?
1828 _delta3->prio() : std::numeric_limits<Value>::max();
1830 Value d4 = !_delta4->empty() ?
1831 _delta4->prio() : std::numeric_limits<Value>::max();
1833 _delta_sum = d3; OpType ot = D3;
1834 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1835 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1836 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1841 Node n = _delta1->top();
1848 int blossom = _delta2->top();
1849 Node n = _blossom_set->classTop(blossom);
1850 Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1851 if ((*_blossom_data)[blossom].next == INVALID) {
1861 Edge e = _delta3->top();
1863 int left_blossom = _blossom_set->find(_graph.u(e));
1864 int right_blossom = _blossom_set->find(_graph.v(e));
1866 if (left_blossom == right_blossom) {
1869 int left_tree = _tree_set->find(left_blossom);
1870 int right_tree = _tree_set->find(right_blossom);
1872 if (left_tree == right_tree) {
1873 shrinkOnEdge(e, left_tree);
1881 splitBlossom(_delta4->top());
1888 /// \brief Run the algorithm.
1890 /// This method runs the \c %MaxWeightedMatching algorithm.
1892 /// \note mwm.run() is just a shortcut of the following code.
1894 /// mwm.fractionalInit();
1904 /// \name Primal Solution
1905 /// Functions to get the primal solution, i.e. the maximum weighted
1907 /// Either \ref run() or \ref start() function should be called before
1912 /// \brief Return the weight of the matching.
1914 /// This function returns the weight of the found matching.
1916 /// \pre Either run() or start() must be called before using this function.
1917 Value matchingWeight() const {
1919 for (NodeIt n(_graph); n != INVALID; ++n) {
1920 if ((*_matching)[n] != INVALID) {
1921 sum += _weight[(*_matching)[n]];
1927 /// \brief Return the size (cardinality) of the matching.
1929 /// This function returns the size (cardinality) of the found matching.
1931 /// \pre Either run() or start() must be called before using this function.
1932 int matchingSize() const {
1934 for (NodeIt n(_graph); n != INVALID; ++n) {
1935 if ((*_matching)[n] != INVALID) {
1942 /// \brief Return \c true if the given edge is in the matching.
1944 /// This function returns \c true if the given edge is in the found
1947 /// \pre Either run() or start() must be called before using this function.
1948 bool matching(const Edge& edge) const {
1949 return edge == (*_matching)[_graph.u(edge)];
1952 /// \brief Return the matching arc (or edge) incident to the given node.
1954 /// This function returns the matching arc (or edge) incident to the
1955 /// given node in the found matching or \c INVALID if the node is
1956 /// not covered by the matching.
1958 /// \pre Either run() or start() must be called before using this function.
1959 Arc matching(const Node& node) const {
1960 return (*_matching)[node];
1963 /// \brief Return a const reference to the matching map.
1965 /// This function returns a const reference to a node map that stores
1966 /// the matching arc (or edge) incident to each node.
1967 const MatchingMap& matchingMap() const {
1971 /// \brief Return the mate of the given node.
1973 /// This function returns the mate of the given node in the found
1974 /// matching or \c INVALID if the node is not covered by the matching.
1976 /// \pre Either run() or start() must be called before using this function.
1977 Node mate(const Node& node) const {
1978 return (*_matching)[node] != INVALID ?
1979 _graph.target((*_matching)[node]) : INVALID;
1984 /// \name Dual Solution
1985 /// Functions to get the dual solution.\n
1986 /// Either \ref run() or \ref start() function should be called before
1991 /// \brief Return the value of the dual solution.
1993 /// This function returns the value of the dual solution.
1994 /// It should be equal to the primal value scaled by \ref dualScale
1997 /// \pre Either run() or start() must be called before using this function.
1998 Value dualValue() const {
2000 for (NodeIt n(_graph); n != INVALID; ++n) {
2001 sum += nodeValue(n);
2003 for (int i = 0; i < blossomNum(); ++i) {
2004 sum += blossomValue(i) * (blossomSize(i) / 2);
2009 /// \brief Return the dual value (potential) of the given node.
2011 /// This function returns the dual value (potential) of the given node.
2013 /// \pre Either run() or start() must be called before using this function.
2014 Value nodeValue(const Node& n) const {
2015 return (*_node_potential)[n];
2018 /// \brief Return the number of the blossoms in the basis.
2020 /// This function returns the number of the blossoms in the basis.
2022 /// \pre Either run() or start() must be called before using this function.
2024 int blossomNum() const {
2025 return _blossom_potential.size();
2028 /// \brief Return the number of the nodes in the given blossom.
2030 /// This function returns the number of the nodes in the given blossom.
2032 /// \pre Either run() or start() must be called before using this function.
2034 int blossomSize(int k) const {
2035 return _blossom_potential[k].end - _blossom_potential[k].begin;
2038 /// \brief Return the dual value (ptential) of the given blossom.
2040 /// This function returns the dual value (ptential) of the given blossom.
2042 /// \pre Either run() or start() must be called before using this function.
2043 Value blossomValue(int k) const {
2044 return _blossom_potential[k].value;
2047 /// \brief Iterator for obtaining the nodes of a blossom.
2049 /// This class provides an iterator for obtaining the nodes of the
2050 /// given blossom. It lists a subset of the nodes.
2051 /// Before using this iterator, you must allocate a
2052 /// MaxWeightedMatching class and execute it.
2056 /// \brief Constructor.
2058 /// Constructor to get the nodes of the given variable.
2060 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
2061 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
2062 /// called before initializing this iterator.
2063 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
2064 : _algorithm(&algorithm)
2066 _index = _algorithm->_blossom_potential[variable].begin;
2067 _last = _algorithm->_blossom_potential[variable].end;
2070 /// \brief Conversion to \c Node.
2072 /// Conversion to \c Node.
2073 operator Node() const {
2074 return _algorithm->_blossom_node_list[_index];
2077 /// \brief Increment operator.
2079 /// Increment operator.
2080 BlossomIt& operator++() {
2085 /// \brief Validity checking
2087 /// Checks whether the iterator is invalid.
2088 bool operator==(Invalid) const { return _index == _last; }
2090 /// \brief Validity checking
2092 /// Checks whether the iterator is valid.
2093 bool operator!=(Invalid) const { return _index != _last; }
2096 const MaxWeightedMatching* _algorithm;
2105 /// \ingroup matching
2107 /// \brief Weighted perfect matching in general graphs
2109 /// This class provides an efficient implementation of Edmond's
2110 /// maximum weighted perfect matching algorithm. The implementation
2111 /// is based on extensive use of priority queues and provides
2112 /// \f$O(nm\log n)\f$ time complexity.
2114 /// The maximum weighted perfect matching problem is to find a subset of
2115 /// the edges in an undirected graph with maximum overall weight for which
2116 /// each node has exactly one incident edge.
2117 /// It can be formulated with the following linear program.
2118 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2119 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2120 \quad \forall B\in\mathcal{O}\f] */
2121 /// \f[x_e \ge 0\quad \forall e\in E\f]
2122 /// \f[\max \sum_{e\in E}x_ew_e\f]
2123 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2124 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2125 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2126 /// subsets of the nodes.
2128 /// The algorithm calculates an optimal matching and a proof of the
2129 /// optimality. The solution of the dual problem can be used to check
2130 /// the result of the algorithm. The dual linear problem is the
2132 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2133 w_{uv} \quad \forall uv\in E\f] */
2134 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2135 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2136 \frac{\vert B \vert - 1}{2}z_B\f] */
2138 /// The algorithm can be executed with the run() function.
2139 /// After it the matching (the primal solution) and the dual solution
2140 /// can be obtained using the query functions and the
2141 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2142 /// which is able to iterate on the nodes of a blossom.
2143 /// If the value type is integer, then the dual solution is multiplied
2144 /// by \ref MaxWeightedMatching::dualScale "4".
2146 /// \tparam GR The undirected graph type the algorithm runs on.
2147 /// \tparam WM The type edge weight map. The default type is
2148 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2150 template <typename GR, typename WM>
2152 template <typename GR,
2153 typename WM = typename GR::template EdgeMap<int> >
2155 class MaxWeightedPerfectMatching {
2158 /// The graph type of the algorithm
2160 /// The type of the edge weight map
2161 typedef WM WeightMap;
2162 /// The value type of the edge weights
2163 typedef typename WeightMap::Value Value;
2165 /// \brief Scaling factor for dual solution
2167 /// Scaling factor for dual solution, it is equal to 4 or 1
2168 /// according to the value type.
2169 static const int dualScale =
2170 std::numeric_limits<Value>::is_integer ? 4 : 1;
2172 /// The type of the matching map
2173 typedef typename Graph::template NodeMap<typename Graph::Arc>
2178 TEMPLATE_GRAPH_TYPEDEFS(Graph);
2180 typedef typename Graph::template NodeMap<Value> NodePotential;
2181 typedef std::vector<Node> BlossomNodeList;
2183 struct BlossomVariable {
2187 BlossomVariable(int _begin, int _end, Value _value)
2188 : begin(_begin), end(_end), value(_value) {}
2192 typedef std::vector<BlossomVariable> BlossomPotential;
2194 const Graph& _graph;
2195 const WeightMap& _weight;
2197 MatchingMap* _matching;
2199 NodePotential* _node_potential;
2201 BlossomPotential _blossom_potential;
2202 BlossomNodeList _blossom_node_list;
2207 typedef RangeMap<int> IntIntMap;
2210 EVEN = -1, MATCHED = 0, ODD = 1
2213 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2214 struct BlossomData {
2221 IntNodeMap *_blossom_index;
2222 BlossomSet *_blossom_set;
2223 RangeMap<BlossomData>* _blossom_data;
2225 IntNodeMap *_node_index;
2226 IntArcMap *_node_heap_index;
2230 NodeData(IntArcMap& node_heap_index)
2231 : heap(node_heap_index) {}
2235 BinHeap<Value, IntArcMap> heap;
2236 std::map<int, Arc> heap_index;
2241 RangeMap<NodeData>* _node_data;
2243 typedef ExtendFindEnum<IntIntMap> TreeSet;
2245 IntIntMap *_tree_set_index;
2248 IntIntMap *_delta2_index;
2249 BinHeap<Value, IntIntMap> *_delta2;
2251 IntEdgeMap *_delta3_index;
2252 BinHeap<Value, IntEdgeMap> *_delta3;
2254 IntIntMap *_delta4_index;
2255 BinHeap<Value, IntIntMap> *_delta4;
2260 typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
2262 FractionalMatching *_fractional;
2264 void createStructures() {
2265 _node_num = countNodes(_graph);
2266 _blossom_num = _node_num * 3 / 2;
2269 _matching = new MatchingMap(_graph);
2272 if (!_node_potential) {
2273 _node_potential = new NodePotential(_graph);
2276 if (!_blossom_set) {
2277 _blossom_index = new IntNodeMap(_graph);
2278 _blossom_set = new BlossomSet(*_blossom_index);
2279 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2280 } else if (_blossom_data->size() != _blossom_num) {
2281 delete _blossom_data;
2282 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2286 _node_index = new IntNodeMap(_graph);
2287 _node_heap_index = new IntArcMap(_graph);
2288 _node_data = new RangeMap<NodeData>(_node_num,
2289 NodeData(*_node_heap_index));
2290 } else if (_node_data->size() != _node_num) {
2292 _node_data = new RangeMap<NodeData>(_node_num,
2293 NodeData(*_node_heap_index));
2297 _tree_set_index = new IntIntMap(_blossom_num);
2298 _tree_set = new TreeSet(*_tree_set_index);
2300 _tree_set_index->resize(_blossom_num);
2304 _delta2_index = new IntIntMap(_blossom_num);
2305 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2307 _delta2_index->resize(_blossom_num);
2311 _delta3_index = new IntEdgeMap(_graph);
2312 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2316 _delta4_index = new IntIntMap(_blossom_num);
2317 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2319 _delta4_index->resize(_blossom_num);
2323 void destroyStructures() {
2327 if (_node_potential) {
2328 delete _node_potential;
2331 delete _blossom_index;
2332 delete _blossom_set;
2333 delete _blossom_data;
2338 delete _node_heap_index;
2343 delete _tree_set_index;
2347 delete _delta2_index;
2351 delete _delta3_index;
2355 delete _delta4_index;
2360 void matchedToEven(int blossom, int tree) {
2361 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2362 _delta2->erase(blossom);
2365 if (!_blossom_set->trivial(blossom)) {
2366 (*_blossom_data)[blossom].pot -=
2367 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2370 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2371 n != INVALID; ++n) {
2373 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2374 int ni = (*_node_index)[n];
2376 (*_node_data)[ni].heap.clear();
2377 (*_node_data)[ni].heap_index.clear();
2379 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2381 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2382 Node v = _graph.source(e);
2383 int vb = _blossom_set->find(v);
2384 int vi = (*_node_index)[v];
2386 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2387 dualScale * _weight[e];
2389 if ((*_blossom_data)[vb].status == EVEN) {
2390 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2391 _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 matchedToOdd(int blossom) {
2429 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2430 _delta2->erase(blossom);
2432 (*_blossom_data)[blossom].offset += _delta_sum;
2433 if (!_blossom_set->trivial(blossom)) {
2434 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2435 (*_blossom_data)[blossom].offset);
2439 void evenToMatched(int blossom, int tree) {
2440 if (!_blossom_set->trivial(blossom)) {
2441 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2444 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2445 n != INVALID; ++n) {
2446 int ni = (*_node_index)[n];
2447 (*_node_data)[ni].pot -= _delta_sum;
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 (vb == blossom) {
2458 if (_delta3->state(e) == _delta3->IN_HEAP) {
2461 } else if ((*_blossom_data)[vb].status == EVEN) {
2463 if (_delta3->state(e) == _delta3->IN_HEAP) {
2467 int vt = _tree_set->find(vb);
2471 Arc r = _graph.oppositeArc(e);
2473 typename std::map<int, Arc>::iterator it =
2474 (*_node_data)[ni].heap_index.find(vt);
2476 if (it != (*_node_data)[ni].heap_index.end()) {
2477 if ((*_node_data)[ni].heap[it->second] > rw) {
2478 (*_node_data)[ni].heap.replace(it->second, r);
2479 (*_node_data)[ni].heap.decrease(r, rw);
2483 (*_node_data)[ni].heap.push(r, rw);
2484 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2487 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2488 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2490 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2491 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2492 (*_blossom_data)[blossom].offset);
2493 } else if ((*_delta2)[blossom] >
2494 _blossom_set->classPrio(blossom) -
2495 (*_blossom_data)[blossom].offset){
2496 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2497 (*_blossom_data)[blossom].offset);
2503 typename std::map<int, Arc>::iterator it =
2504 (*_node_data)[vi].heap_index.find(tree);
2506 if (it != (*_node_data)[vi].heap_index.end()) {
2507 (*_node_data)[vi].heap.erase(it->second);
2508 (*_node_data)[vi].heap_index.erase(it);
2509 if ((*_node_data)[vi].heap.empty()) {
2510 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2511 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2512 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2515 if ((*_blossom_data)[vb].status == MATCHED) {
2516 if (_blossom_set->classPrio(vb) ==
2517 std::numeric_limits<Value>::max()) {
2519 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2520 (*_blossom_data)[vb].offset) {
2521 _delta2->increase(vb, _blossom_set->classPrio(vb) -
2522 (*_blossom_data)[vb].offset);
2531 void oddToMatched(int blossom) {
2532 (*_blossom_data)[blossom].offset -= _delta_sum;
2534 if (_blossom_set->classPrio(blossom) !=
2535 std::numeric_limits<Value>::max()) {
2536 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2537 (*_blossom_data)[blossom].offset);
2540 if (!_blossom_set->trivial(blossom)) {
2541 _delta4->erase(blossom);
2545 void oddToEven(int blossom, int tree) {
2546 if (!_blossom_set->trivial(blossom)) {
2547 _delta4->erase(blossom);
2548 (*_blossom_data)[blossom].pot -=
2549 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2552 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2553 n != INVALID; ++n) {
2554 int ni = (*_node_index)[n];
2556 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2558 (*_node_data)[ni].heap.clear();
2559 (*_node_data)[ni].heap_index.clear();
2560 (*_node_data)[ni].pot +=
2561 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2563 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2564 Node v = _graph.source(e);
2565 int vb = _blossom_set->find(v);
2566 int vi = (*_node_index)[v];
2568 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2569 dualScale * _weight[e];
2571 if ((*_blossom_data)[vb].status == EVEN) {
2572 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2573 _delta3->push(e, rw / 2);
2577 typename std::map<int, Arc>::iterator it =
2578 (*_node_data)[vi].heap_index.find(tree);
2580 if (it != (*_node_data)[vi].heap_index.end()) {
2581 if ((*_node_data)[vi].heap[it->second] > rw) {
2582 (*_node_data)[vi].heap.replace(it->second, e);
2583 (*_node_data)[vi].heap.decrease(e, rw);
2587 (*_node_data)[vi].heap.push(e, rw);
2588 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2591 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2592 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2594 if ((*_blossom_data)[vb].status == MATCHED) {
2595 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2596 _delta2->push(vb, _blossom_set->classPrio(vb) -
2597 (*_blossom_data)[vb].offset);
2598 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2599 (*_blossom_data)[vb].offset) {
2600 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2601 (*_blossom_data)[vb].offset);
2608 (*_blossom_data)[blossom].offset = 0;
2611 void alternatePath(int even, int tree) {
2614 evenToMatched(even, tree);
2615 (*_blossom_data)[even].status = MATCHED;
2617 while ((*_blossom_data)[even].pred != INVALID) {
2618 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2619 (*_blossom_data)[odd].status = MATCHED;
2621 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2623 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2624 (*_blossom_data)[even].status = MATCHED;
2625 evenToMatched(even, tree);
2626 (*_blossom_data)[even].next =
2627 _graph.oppositeArc((*_blossom_data)[odd].pred);
2632 void destroyTree(int tree) {
2633 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2634 if ((*_blossom_data)[b].status == EVEN) {
2635 (*_blossom_data)[b].status = MATCHED;
2636 evenToMatched(b, tree);
2637 } else if ((*_blossom_data)[b].status == ODD) {
2638 (*_blossom_data)[b].status = MATCHED;
2642 _tree_set->eraseClass(tree);
2645 void augmentOnEdge(const Edge& edge) {
2647 int left = _blossom_set->find(_graph.u(edge));
2648 int right = _blossom_set->find(_graph.v(edge));
2650 int left_tree = _tree_set->find(left);
2651 alternatePath(left, left_tree);
2652 destroyTree(left_tree);
2654 int right_tree = _tree_set->find(right);
2655 alternatePath(right, right_tree);
2656 destroyTree(right_tree);
2658 (*_blossom_data)[left].next = _graph.direct(edge, true);
2659 (*_blossom_data)[right].next = _graph.direct(edge, false);
2662 void extendOnArc(const Arc& arc) {
2663 int base = _blossom_set->find(_graph.target(arc));
2664 int tree = _tree_set->find(base);
2666 int odd = _blossom_set->find(_graph.source(arc));
2667 _tree_set->insert(odd, tree);
2668 (*_blossom_data)[odd].status = ODD;
2670 (*_blossom_data)[odd].pred = arc;
2672 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2673 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2674 _tree_set->insert(even, tree);
2675 (*_blossom_data)[even].status = EVEN;
2676 matchedToEven(even, tree);
2679 void shrinkOnEdge(const Edge& edge, int tree) {
2681 std::vector<int> left_path, right_path;
2684 std::set<int> left_set, right_set;
2685 int left = _blossom_set->find(_graph.u(edge));
2686 left_path.push_back(left);
2687 left_set.insert(left);
2689 int right = _blossom_set->find(_graph.v(edge));
2690 right_path.push_back(right);
2691 right_set.insert(right);
2695 if ((*_blossom_data)[left].pred == INVALID) break;
2698 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2699 left_path.push_back(left);
2701 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2702 left_path.push_back(left);
2704 left_set.insert(left);
2706 if (right_set.find(left) != right_set.end()) {
2711 if ((*_blossom_data)[right].pred == INVALID) break;
2714 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2715 right_path.push_back(right);
2717 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2718 right_path.push_back(right);
2720 right_set.insert(right);
2722 if (left_set.find(right) != left_set.end()) {
2730 if ((*_blossom_data)[left].pred == INVALID) {
2732 while (left_set.find(nca) == left_set.end()) {
2734 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2735 right_path.push_back(nca);
2737 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2738 right_path.push_back(nca);
2742 while (right_set.find(nca) == right_set.end()) {
2744 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2745 left_path.push_back(nca);
2747 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2748 left_path.push_back(nca);
2754 std::vector<int> subblossoms;
2757 prev = _graph.direct(edge, true);
2758 for (int i = 0; left_path[i] != nca; i += 2) {
2759 subblossoms.push_back(left_path[i]);
2760 (*_blossom_data)[left_path[i]].next = prev;
2761 _tree_set->erase(left_path[i]);
2763 subblossoms.push_back(left_path[i + 1]);
2764 (*_blossom_data)[left_path[i + 1]].status = EVEN;
2765 oddToEven(left_path[i + 1], tree);
2766 _tree_set->erase(left_path[i + 1]);
2767 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2771 while (right_path[k] != nca) ++k;
2773 subblossoms.push_back(nca);
2774 (*_blossom_data)[nca].next = prev;
2776 for (int i = k - 2; i >= 0; i -= 2) {
2777 subblossoms.push_back(right_path[i + 1]);
2778 (*_blossom_data)[right_path[i + 1]].status = EVEN;
2779 oddToEven(right_path[i + 1], tree);
2780 _tree_set->erase(right_path[i + 1]);
2782 (*_blossom_data)[right_path[i + 1]].next =
2783 (*_blossom_data)[right_path[i + 1]].pred;
2785 subblossoms.push_back(right_path[i]);
2786 _tree_set->erase(right_path[i]);
2790 _blossom_set->join(subblossoms.begin(), subblossoms.end());
2792 for (int i = 0; i < int(subblossoms.size()); ++i) {
2793 if (!_blossom_set->trivial(subblossoms[i])) {
2794 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2796 (*_blossom_data)[subblossoms[i]].status = MATCHED;
2799 (*_blossom_data)[surface].pot = -2 * _delta_sum;
2800 (*_blossom_data)[surface].offset = 0;
2801 (*_blossom_data)[surface].status = EVEN;
2802 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2803 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2805 _tree_set->insert(surface, tree);
2806 _tree_set->erase(nca);
2809 void splitBlossom(int blossom) {
2810 Arc next = (*_blossom_data)[blossom].next;
2811 Arc pred = (*_blossom_data)[blossom].pred;
2813 int tree = _tree_set->find(blossom);
2815 (*_blossom_data)[blossom].status = MATCHED;
2816 oddToMatched(blossom);
2817 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2818 _delta2->erase(blossom);
2821 std::vector<int> subblossoms;
2822 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2824 Value offset = (*_blossom_data)[blossom].offset;
2825 int b = _blossom_set->find(_graph.source(pred));
2826 int d = _blossom_set->find(_graph.source(next));
2828 int ib = -1, id = -1;
2829 for (int i = 0; i < int(subblossoms.size()); ++i) {
2830 if (subblossoms[i] == b) ib = i;
2831 if (subblossoms[i] == d) id = i;
2833 (*_blossom_data)[subblossoms[i]].offset = offset;
2834 if (!_blossom_set->trivial(subblossoms[i])) {
2835 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2837 if (_blossom_set->classPrio(subblossoms[i]) !=
2838 std::numeric_limits<Value>::max()) {
2839 _delta2->push(subblossoms[i],
2840 _blossom_set->classPrio(subblossoms[i]) -
2841 (*_blossom_data)[subblossoms[i]].offset);
2845 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2846 for (int i = (id + 1) % subblossoms.size();
2847 i != ib; i = (i + 2) % subblossoms.size()) {
2848 int sb = subblossoms[i];
2849 int tb = subblossoms[(i + 1) % subblossoms.size()];
2850 (*_blossom_data)[sb].next =
2851 _graph.oppositeArc((*_blossom_data)[tb].next);
2854 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2855 int sb = subblossoms[i];
2856 int tb = subblossoms[(i + 1) % subblossoms.size()];
2857 int ub = subblossoms[(i + 2) % subblossoms.size()];
2859 (*_blossom_data)[sb].status = ODD;
2861 _tree_set->insert(sb, tree);
2862 (*_blossom_data)[sb].pred = pred;
2863 (*_blossom_data)[sb].next =
2864 _graph.oppositeArc((*_blossom_data)[tb].next);
2866 pred = (*_blossom_data)[ub].next;
2868 (*_blossom_data)[tb].status = EVEN;
2869 matchedToEven(tb, tree);
2870 _tree_set->insert(tb, tree);
2871 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2874 (*_blossom_data)[subblossoms[id]].status = ODD;
2875 matchedToOdd(subblossoms[id]);
2876 _tree_set->insert(subblossoms[id], tree);
2877 (*_blossom_data)[subblossoms[id]].next = next;
2878 (*_blossom_data)[subblossoms[id]].pred = pred;
2882 for (int i = (ib + 1) % subblossoms.size();
2883 i != id; i = (i + 2) % subblossoms.size()) {
2884 int sb = subblossoms[i];
2885 int tb = subblossoms[(i + 1) % subblossoms.size()];
2886 (*_blossom_data)[sb].next =
2887 _graph.oppositeArc((*_blossom_data)[tb].next);
2890 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2891 int sb = subblossoms[i];
2892 int tb = subblossoms[(i + 1) % subblossoms.size()];
2893 int ub = subblossoms[(i + 2) % subblossoms.size()];
2895 (*_blossom_data)[sb].status = ODD;
2897 _tree_set->insert(sb, tree);
2898 (*_blossom_data)[sb].next = next;
2899 (*_blossom_data)[sb].pred =
2900 _graph.oppositeArc((*_blossom_data)[tb].next);
2902 (*_blossom_data)[tb].status = EVEN;
2903 matchedToEven(tb, tree);
2904 _tree_set->insert(tb, tree);
2905 (*_blossom_data)[tb].pred =
2906 (*_blossom_data)[tb].next =
2907 _graph.oppositeArc((*_blossom_data)[ub].next);
2908 next = (*_blossom_data)[ub].next;
2911 (*_blossom_data)[subblossoms[ib]].status = ODD;
2912 matchedToOdd(subblossoms[ib]);
2913 _tree_set->insert(subblossoms[ib], tree);
2914 (*_blossom_data)[subblossoms[ib]].next = next;
2915 (*_blossom_data)[subblossoms[ib]].pred = pred;
2917 _tree_set->erase(blossom);
2920 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2921 if (_blossom_set->trivial(blossom)) {
2922 int bi = (*_node_index)[base];
2923 Value pot = (*_node_data)[bi].pot;
2925 (*_matching)[base] = matching;
2926 _blossom_node_list.push_back(base);
2927 (*_node_potential)[base] = pot;
2930 Value pot = (*_blossom_data)[blossom].pot;
2931 int bn = _blossom_node_list.size();
2933 std::vector<int> subblossoms;
2934 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2935 int b = _blossom_set->find(base);
2937 for (int i = 0; i < int(subblossoms.size()); ++i) {
2938 if (subblossoms[i] == b) { ib = i; break; }
2941 for (int i = 1; i < int(subblossoms.size()); i += 2) {
2942 int sb = subblossoms[(ib + i) % subblossoms.size()];
2943 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2945 Arc m = (*_blossom_data)[tb].next;
2946 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2947 extractBlossom(tb, _graph.source(m), m);
2949 extractBlossom(subblossoms[ib], base, matching);
2951 int en = _blossom_node_list.size();
2953 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2957 void extractMatching() {
2958 std::vector<int> blossoms;
2959 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2960 blossoms.push_back(c);
2963 for (int i = 0; i < int(blossoms.size()); ++i) {
2965 Value offset = (*_blossom_data)[blossoms[i]].offset;
2966 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2967 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2968 n != INVALID; ++n) {
2969 (*_node_data)[(*_node_index)[n]].pot -= offset;
2972 Arc matching = (*_blossom_data)[blossoms[i]].next;
2973 Node base = _graph.source(matching);
2974 extractBlossom(blossoms[i], base, matching);
2980 /// \brief Constructor
2983 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2984 : _graph(graph), _weight(weight), _matching(0),
2985 _node_potential(0), _blossom_potential(), _blossom_node_list(),
2986 _node_num(0), _blossom_num(0),
2988 _blossom_index(0), _blossom_set(0), _blossom_data(0),
2989 _node_index(0), _node_heap_index(0), _node_data(0),
2990 _tree_set_index(0), _tree_set(0),
2992 _delta2_index(0), _delta2(0),
2993 _delta3_index(0), _delta3(0),
2994 _delta4_index(0), _delta4(0),
2996 _delta_sum(), _unmatched(0),
3001 ~MaxWeightedPerfectMatching() {
3002 destroyStructures();
3008 /// \name Execution Control
3009 /// The simplest way to execute the algorithm is to use the
3010 /// \ref run() member function.
3014 /// \brief Initialize the algorithm
3016 /// This function initializes the algorithm.
3020 _blossom_node_list.clear();
3021 _blossom_potential.clear();
3023 for (ArcIt e(_graph); e != INVALID; ++e) {
3024 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
3026 for (EdgeIt e(_graph); e != INVALID; ++e) {
3027 (*_delta3_index)[e] = _delta3->PRE_HEAP;
3029 for (int i = 0; i < _blossom_num; ++i) {
3030 (*_delta2_index)[i] = _delta2->PRE_HEAP;
3031 (*_delta4_index)[i] = _delta4->PRE_HEAP;
3034 _unmatched = _node_num;
3039 _blossom_set->clear();
3043 for (NodeIt n(_graph); n != INVALID; ++n) {
3044 Value max = - std::numeric_limits<Value>::max();
3045 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3046 if (_graph.target(e) == n) continue;
3047 if ((dualScale * _weight[e]) / 2 > max) {
3048 max = (dualScale * _weight[e]) / 2;
3051 (*_node_index)[n] = index;
3052 (*_node_data)[index].heap_index.clear();
3053 (*_node_data)[index].heap.clear();
3054 (*_node_data)[index].pot = max;
3056 _blossom_set->insert(n, std::numeric_limits<Value>::max());
3058 _tree_set->insert(blossom);
3060 (*_blossom_data)[blossom].status = EVEN;
3061 (*_blossom_data)[blossom].pred = INVALID;
3062 (*_blossom_data)[blossom].next = INVALID;
3063 (*_blossom_data)[blossom].pot = 0;
3064 (*_blossom_data)[blossom].offset = 0;
3067 for (EdgeIt e(_graph); e != INVALID; ++e) {
3068 int si = (*_node_index)[_graph.u(e)];
3069 int ti = (*_node_index)[_graph.v(e)];
3070 if (_graph.u(e) != _graph.v(e)) {
3071 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3072 dualScale * _weight[e]) / 2);
3077 /// \brief Initialize the algorithm with fractional matching
3079 /// This function initializes the algorithm with a fractional
3080 /// matching. This initialization is also called jumpstart heuristic.
3081 void fractionalInit() {
3084 if (_fractional == 0) {
3085 _fractional = new FractionalMatching(_graph, _weight, false);
3087 if (!_fractional->run()) {
3092 for (ArcIt e(_graph); e != INVALID; ++e) {
3093 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
3095 for (EdgeIt e(_graph); e != INVALID; ++e) {
3096 (*_delta3_index)[e] = _delta3->PRE_HEAP;
3098 for (int i = 0; i < _blossom_num; ++i) {
3099 (*_delta2_index)[i] = _delta2->PRE_HEAP;
3100 (*_delta4_index)[i] = _delta4->PRE_HEAP;
3106 for (NodeIt n(_graph); n != INVALID; ++n) {
3107 Value pot = _fractional->nodeValue(n);
3108 (*_node_index)[n] = index;
3109 (*_node_data)[index].pot = pot;
3111 _blossom_set->insert(n, std::numeric_limits<Value>::max());
3113 (*_blossom_data)[blossom].status = MATCHED;
3114 (*_blossom_data)[blossom].pred = INVALID;
3115 (*_blossom_data)[blossom].next = _fractional->matching(n);
3116 (*_blossom_data)[blossom].pot = 0;
3117 (*_blossom_data)[blossom].offset = 0;
3121 typename Graph::template NodeMap<bool> processed(_graph, false);
3122 for (NodeIt n(_graph); n != INVALID; ++n) {
3123 if (processed[n]) continue;
3124 processed[n] = true;
3125 if (_fractional->matching(n) == INVALID) continue;
3127 Node v = _graph.target(_fractional->matching(n));
3129 processed[v] = true;
3130 v = _graph.target(_fractional->matching(v));
3135 std::vector<int> subblossoms(num);
3137 subblossoms[--num] = _blossom_set->find(n);
3138 v = _graph.target(_fractional->matching(n));
3140 subblossoms[--num] = _blossom_set->find(v);
3141 v = _graph.target(_fractional->matching(v));
3145 _blossom_set->join(subblossoms.begin(), subblossoms.end());
3146 (*_blossom_data)[surface].status = EVEN;
3147 (*_blossom_data)[surface].pred = INVALID;
3148 (*_blossom_data)[surface].next = INVALID;
3149 (*_blossom_data)[surface].pot = 0;
3150 (*_blossom_data)[surface].offset = 0;
3152 _tree_set->insert(surface);
3157 for (EdgeIt e(_graph); e != INVALID; ++e) {
3158 int si = (*_node_index)[_graph.u(e)];
3159 int sb = _blossom_set->find(_graph.u(e));
3160 int ti = (*_node_index)[_graph.v(e)];
3161 int tb = _blossom_set->find(_graph.v(e));
3162 if ((*_blossom_data)[sb].status == EVEN &&
3163 (*_blossom_data)[tb].status == EVEN && sb != tb) {
3164 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3165 dualScale * _weight[e]) / 2);
3169 for (NodeIt n(_graph); n != INVALID; ++n) {
3170 int nb = _blossom_set->find(n);
3171 if ((*_blossom_data)[nb].status != MATCHED) continue;
3172 int ni = (*_node_index)[n];
3174 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3175 Node v = _graph.target(e);
3176 int vb = _blossom_set->find(v);
3177 int vi = (*_node_index)[v];
3179 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3180 dualScale * _weight[e];
3182 if ((*_blossom_data)[vb].status == EVEN) {
3184 int vt = _tree_set->find(vb);
3186 typename std::map<int, Arc>::iterator it =
3187 (*_node_data)[ni].heap_index.find(vt);
3189 if (it != (*_node_data)[ni].heap_index.end()) {
3190 if ((*_node_data)[ni].heap[it->second] > rw) {
3191 (*_node_data)[ni].heap.replace(it->second, e);
3192 (*_node_data)[ni].heap.decrease(e, rw);
3196 (*_node_data)[ni].heap.push(e, rw);
3197 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
3202 if (!(*_node_data)[ni].heap.empty()) {
3203 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3204 _delta2->push(nb, _blossom_set->classPrio(nb));
3209 /// \brief Start the algorithm
3211 /// This function starts the algorithm.
3213 /// \pre \ref init() or \ref fractionalInit() must be called before
3214 /// using this function.
3220 if (_unmatched == -1) return false;
3222 while (_unmatched > 0) {
3223 Value d2 = !_delta2->empty() ?
3224 _delta2->prio() : std::numeric_limits<Value>::max();
3226 Value d3 = !_delta3->empty() ?
3227 _delta3->prio() : std::numeric_limits<Value>::max();
3229 Value d4 = !_delta4->empty() ?
3230 _delta4->prio() : std::numeric_limits<Value>::max();
3232 _delta_sum = d3; OpType ot = D3;
3233 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3234 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3236 if (_delta_sum == std::numeric_limits<Value>::max()) {
3243 int blossom = _delta2->top();
3244 Node n = _blossom_set->classTop(blossom);
3245 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
3251 Edge e = _delta3->top();
3253 int left_blossom = _blossom_set->find(_graph.u(e));
3254 int right_blossom = _blossom_set->find(_graph.v(e));
3256 if (left_blossom == right_blossom) {
3259 int left_tree = _tree_set->find(left_blossom);
3260 int right_tree = _tree_set->find(right_blossom);
3262 if (left_tree == right_tree) {
3263 shrinkOnEdge(e, left_tree);
3271 splitBlossom(_delta4->top());
3279 /// \brief Run the algorithm.
3281 /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3283 /// \note mwpm.run() is just a shortcut of the following code.
3285 /// mwpm.fractionalInit();
3295 /// \name Primal Solution
3296 /// Functions to get the primal solution, i.e. the maximum weighted
3297 /// perfect matching.\n
3298 /// Either \ref run() or \ref start() function should be called before
3303 /// \brief Return the weight of the matching.
3305 /// This function returns the weight of the found matching.
3307 /// \pre Either run() or start() must be called before using this function.
3308 Value matchingWeight() const {
3310 for (NodeIt n(_graph); n != INVALID; ++n) {
3311 if ((*_matching)[n] != INVALID) {
3312 sum += _weight[(*_matching)[n]];
3318 /// \brief Return \c true if the given edge is in the matching.
3320 /// This function returns \c true if the given edge is in the found
3323 /// \pre Either run() or start() must be called before using this function.
3324 bool matching(const Edge& edge) const {
3325 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3328 /// \brief Return the matching arc (or edge) incident to the given node.
3330 /// This function returns the matching arc (or edge) incident to the
3331 /// given node in the found matching or \c INVALID if the node is
3332 /// not covered by the matching.
3334 /// \pre Either run() or start() must be called before using this function.
3335 Arc matching(const Node& node) const {
3336 return (*_matching)[node];
3339 /// \brief Return a const reference to the matching map.
3341 /// This function returns a const reference to a node map that stores
3342 /// the matching arc (or edge) incident to each node.
3343 const MatchingMap& matchingMap() const {
3347 /// \brief Return the mate of the given node.
3349 /// This function returns the mate of the given node in the found
3350 /// matching or \c INVALID if the node is not covered by the matching.
3352 /// \pre Either run() or start() must be called before using this function.
3353 Node mate(const Node& node) const {
3354 return _graph.target((*_matching)[node]);
3359 /// \name Dual Solution
3360 /// Functions to get the dual solution.\n
3361 /// Either \ref run() or \ref start() function should be called before
3366 /// \brief Return the value of the dual solution.
3368 /// This function returns the value of the dual solution.
3369 /// It should be equal to the primal value scaled by \ref dualScale
3372 /// \pre Either run() or start() must be called before using this function.
3373 Value dualValue() const {
3375 for (NodeIt n(_graph); n != INVALID; ++n) {
3376 sum += nodeValue(n);
3378 for (int i = 0; i < blossomNum(); ++i) {
3379 sum += blossomValue(i) * (blossomSize(i) / 2);
3384 /// \brief Return the dual value (potential) of the given node.
3386 /// This function returns the dual value (potential) of the given node.
3388 /// \pre Either run() or start() must be called before using this function.
3389 Value nodeValue(const Node& n) const {
3390 return (*_node_potential)[n];
3393 /// \brief Return the number of the blossoms in the basis.
3395 /// This function returns the number of the blossoms in the basis.
3397 /// \pre Either run() or start() must be called before using this function.
3399 int blossomNum() const {
3400 return _blossom_potential.size();
3403 /// \brief Return the number of the nodes in the given blossom.
3405 /// This function returns the number of the nodes in the given blossom.
3407 /// \pre Either run() or start() must be called before using this function.
3409 int blossomSize(int k) const {
3410 return _blossom_potential[k].end - _blossom_potential[k].begin;
3413 /// \brief Return the dual value (ptential) of the given blossom.
3415 /// This function returns the dual value (ptential) of the given blossom.
3417 /// \pre Either run() or start() must be called before using this function.
3418 Value blossomValue(int k) const {
3419 return _blossom_potential[k].value;
3422 /// \brief Iterator for obtaining the nodes of a blossom.
3424 /// This class provides an iterator for obtaining the nodes of the
3425 /// given blossom. It lists a subset of the nodes.
3426 /// Before using this iterator, you must allocate a
3427 /// MaxWeightedPerfectMatching class and execute it.
3431 /// \brief Constructor.
3433 /// Constructor to get the nodes of the given variable.
3435 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3436 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3437 /// must be called before initializing this iterator.
3438 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3439 : _algorithm(&algorithm)
3441 _index = _algorithm->_blossom_potential[variable].begin;
3442 _last = _algorithm->_blossom_potential[variable].end;
3445 /// \brief Conversion to \c Node.
3447 /// Conversion to \c Node.
3448 operator Node() const {
3449 return _algorithm->_blossom_node_list[_index];
3452 /// \brief Increment operator.
3454 /// Increment operator.
3455 BlossomIt& operator++() {
3460 /// \brief Validity checking
3462 /// This function checks whether the iterator is invalid.
3463 bool operator==(Invalid) const { return _index == _last; }
3465 /// \brief Validity checking
3467 /// This function checks whether the iterator is valid.
3468 bool operator!=(Invalid) const { return _index != _last; }
3471 const MaxWeightedPerfectMatching* _algorithm;
3480 } //END OF NAMESPACE LEMON
3482 #endif //LEMON_MATCHING_H