Fix the rename script.
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
26 #include <lemon/bits/invalid.h>
27 #include <lemon/unionfind.h>
28 #include <lemon/graph_utils.h>
29 #include <lemon/bin_heap.h>
33 ///\brief Maximum matching algorithms in undirected graph.
39 ///\brief Edmonds' alternating forest maximum matching algorithm.
41 ///This class provides Edmonds' alternating forest matching
42 ///algorithm. The starting matching (if any) can be passed to the
43 ///algorithm using some of init functions.
45 ///The dual side of a matching is a map of the nodes to
46 ///MaxMatching::DecompType, having values \c D, \c A and \c C
47 ///showing the Gallai-Edmonds decomposition of the graph. The nodes
48 ///in \c D induce a graph with factor-critical components, the nodes
49 ///in \c A form the barrier, and the nodes in \c C induce a graph
50 ///having a perfect matching. This decomposition can be attained by
51 ///calling \c decomposition() after running the algorithm.
53 ///\param Graph The undirected graph type the algorithm runs on.
55 ///\author Jacint Szabo
56 template <typename Graph>
61 typedef typename Graph::Node Node;
62 typedef typename Graph::Edge Edge;
63 typedef typename Graph::UEdge UEdge;
64 typedef typename Graph::UEdgeIt UEdgeIt;
65 typedef typename Graph::NodeIt NodeIt;
66 typedef typename Graph::IncEdgeIt IncEdgeIt;
68 typedef typename Graph::template NodeMap<int> UFECrossRef;
69 typedef UnionFindEnum<UFECrossRef> UFE;
70 typedef std::vector<Node> NV;
72 typedef typename Graph::template NodeMap<int> EFECrossRef;
73 typedef ExtendFindEnum<EFECrossRef> EFE;
77 ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
79 ///Indicates the Gallai-Edmonds decomposition of the graph, which
80 ///shows an upper bound on the size of a maximum matching. The
81 ///nodes with DecompType \c D induce a graph with factor-critical
82 ///components, the nodes in \c A form the canonical barrier, and the
83 ///nodes in \c C induce a graph having a perfect matching.
92 static const int HEUR_density=2;
94 typename Graph::template NodeMap<Node> _mate;
95 typename Graph::template NodeMap<DecompType> position;
99 MaxMatching(const Graph& _g)
100 : g(_g), _mate(_g), position(_g) {}
102 ///\brief Sets the actual matching to the empty matching.
104 ///Sets the actual matching to the empty matching.
107 for(NodeIt v(g); v!=INVALID; ++v) {
108 _mate.set(v,INVALID);
113 ///\brief Finds a greedy matching for initial matching.
115 ///For initial matchig it finds a maximal greedy matching.
117 for(NodeIt v(g); v!=INVALID; ++v) {
118 _mate.set(v,INVALID);
121 for(NodeIt v(g); v!=INVALID; ++v)
122 if ( _mate[v]==INVALID ) {
123 for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
124 Node y=g.runningNode(e);
125 if ( _mate[y]==INVALID && y!=v ) {
134 ///\brief Initialize the matching from each nodes' mate.
136 ///Initialize the matching from a \c Node valued \c Node map. This
137 ///map must be \e symmetric, i.e. if \c map[u]==v then \c
138 ///map[v]==u must hold, and \c uv will be an edge of the initial
140 template <typename MateMap>
141 void mateMapInit(MateMap& map) {
142 for(NodeIt v(g); v!=INVALID; ++v) {
148 ///\brief Initialize the matching from a node map with the
149 ///incident matching edges.
151 ///Initialize the matching from an \c UEdge valued \c Node map. \c
152 ///map[v] must be an \c UEdge incident to \c v. This map must have
153 ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
154 ///g.oppositeNode(v,map[v])==u holds, and now some edge joining \c
155 ///u to \c v will be an edge of the matching.
156 template<typename MatchingMap>
157 void matchingMapInit(MatchingMap& map) {
158 for(NodeIt v(g); v!=INVALID; ++v) {
162 _mate.set(v,g.oppositeNode(v,e));
164 _mate.set(v,INVALID);
168 ///\brief Initialize the matching from the map containing the
169 ///undirected matching edges.
171 ///Initialize the matching from a \c bool valued \c UEdge map. This
172 ///map must have the property that there are no two incident edges
173 ///\c e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
174 ///map[e]==true form the matching.
175 template <typename MatchingMap>
176 void matchingInit(MatchingMap& map) {
177 for(NodeIt v(g); v!=INVALID; ++v) {
178 _mate.set(v,INVALID);
181 for(UEdgeIt e(g); e!=INVALID; ++e) {
192 ///\brief Runs Edmonds' algorithm.
194 ///Runs Edmonds' algorithm for sparse graphs (number of edges <
195 ///2*number of nodes), and a heuristical Edmonds' algorithm with a
196 ///heuristic of postponing shrinks for dense graphs.
198 if (countUEdges(g) < HEUR_density * countNodes(g)) {
208 ///\brief Starts Edmonds' algorithm.
210 ///If runs the original Edmonds' algorithm.
213 typename Graph::template NodeMap<Node> ear(g,INVALID);
214 //undefined for the base nodes of the blossoms (i.e. for the
215 //representative elements of UFE blossom) and for the nodes in C
217 UFECrossRef blossom_base(g);
218 UFE blossom(blossom_base);
219 NV rep(countNodes(g));
221 EFECrossRef tree_base(g);
224 //If these UFE's would be members of the class then also
225 //blossom_base and tree_base should be a member.
227 //We build only one tree and the other vertices uncovered by the
228 //matching belong to C. (They can be considered as singleton
229 //trees.) If this tree can be augmented or no more
230 //grow/augmentation/shrink is possible then we return to this
232 for(NodeIt v(g); v!=INVALID; ++v) {
233 if (position[v]==C && _mate[v]==INVALID) {
234 rep[blossom.insert(v)] = v;
237 normShrink(v, ear, blossom, rep, tree);
242 ///\brief Starts Edmonds' algorithm.
244 ///It runs Edmonds' algorithm with a heuristic of postponing
245 ///shrinks, giving a faster algorithm for dense graphs.
248 typename Graph::template NodeMap<Node> ear(g,INVALID);
249 //undefined for the base nodes of the blossoms (i.e. for the
250 //representative elements of UFE blossom) and for the nodes in C
252 UFECrossRef blossom_base(g);
253 UFE blossom(blossom_base);
254 NV rep(countNodes(g));
256 EFECrossRef tree_base(g);
259 //If these UFE's would be members of the class then also
260 //blossom_base and tree_base should be a member.
262 //We build only one tree and the other vertices uncovered by the
263 //matching belong to C. (They can be considered as singleton
264 //trees.) If this tree can be augmented or no more
265 //grow/augmentation/shrink is possible then we return to this
267 for(NodeIt v(g); v!=INVALID; ++v) {
268 if ( position[v]==C && _mate[v]==INVALID ) {
269 rep[blossom.insert(v)] = v;
272 lateShrink(v, ear, blossom, rep, tree);
279 ///\brief Returns the size of the actual matching stored.
281 ///Returns the size of the actual matching stored. After \ref
282 ///run() it returns the size of a maximum matching in the graph.
285 for(NodeIt v(g); v!=INVALID; ++v) {
286 if ( _mate[v]!=INVALID ) {
294 ///\brief Returns the mate of a node in the actual matching.
296 ///Returns the mate of a \c node in the actual matching.
297 ///Returns INVALID if the \c node is not covered by the actual matching.
298 Node mate(const Node& node) const {
302 ///\brief Returns the matching edge incident to the given node.
304 ///Returns the matching edge of a \c node in the actual matching.
305 ///Returns INVALID if the \c node is not covered by the actual matching.
306 UEdge matchingEdge(const Node& node) const {
307 if (_mate[node] == INVALID) return INVALID;
308 Node n = node < _mate[node] ? node : _mate[node];
309 for (IncEdgeIt e(g, n); e != INVALID; ++e) {
310 if (g.oppositeNode(n, e) == _mate[n]) {
317 /// \brief Returns the class of the node in the Edmonds-Gallai
320 /// Returns the class of the node in the Edmonds-Gallai
322 DecompType decomposition(const Node& n) {
323 return position[n] == A;
326 /// \brief Returns true when the node is in the barrier.
328 /// Returns true when the node is in the barrier.
329 bool barrier(const Node& n) {
330 return position[n] == A;
333 ///\brief Gives back the matching in a \c Node of mates.
335 ///Writes the stored matching to a \c Node valued \c Node map. The
336 ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
337 ///map[v]==u will hold, and now \c uv is an edge of the matching.
338 template <typename MateMap>
339 void mateMap(MateMap& map) const {
340 for(NodeIt v(g); v!=INVALID; ++v) {
345 ///\brief Gives back the matching in an \c UEdge valued \c Node
348 ///Writes the stored matching to an \c UEdge valued \c Node
349 ///map. \c map[v] will be an \c UEdge incident to \c v. This
350 ///map will have the property that if \c g.oppositeNode(u,map[u])
351 ///== v then \c map[u]==map[v] holds, and now this edge is an edge
353 template <typename MatchingMap>
354 void matchingMap(MatchingMap& map) const {
355 typename Graph::template NodeMap<bool> todo(g,true);
356 for(NodeIt v(g); v!=INVALID; ++v) {
357 if (_mate[v]!=INVALID && v < _mate[v]) {
359 for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
360 if ( g.runningNode(e) == u ) {
373 ///\brief Gives back the matching in a \c bool valued \c UEdge
376 ///Writes the matching stored to a \c bool valued \c Edge
377 ///map. This map will have the property that there are no two
378 ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
379 ///edges \c e with \c map[e]==true form the matching.
380 template<typename MatchingMap>
381 void matching(MatchingMap& map) const {
382 for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
384 typename Graph::template NodeMap<bool> todo(g,true);
385 for(NodeIt v(g); v!=INVALID; ++v) {
386 if ( todo[v] && _mate[v]!=INVALID ) {
388 for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
389 if ( g.runningNode(e) == u ) {
401 ///\brief Returns the canonical decomposition of the graph after running
404 ///After calling any run methods of the class, it writes the
405 ///Gallai-Edmonds canonical decomposition of the graph. \c map
406 ///must be a node map of \ref DecompType 's.
407 template <typename DecompositionMap>
408 void decomposition(DecompositionMap& map) const {
409 for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
412 ///\brief Returns a barrier on the nodes.
414 ///After calling any run methods of the class, it writes a
415 ///canonical barrier on the nodes. The odd component number of the
416 ///remaining graph minus the barrier size is a lower bound for the
417 ///uncovered nodes in the graph. The \c map must be a node map of
419 template <typename BarrierMap>
420 void barrier(BarrierMap& barrier) {
421 for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
427 void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
428 UFE& blossom, NV& rep, EFE& tree) {
429 //We have one tree which we grow, and also shrink but only if it
430 //cannot be postponed. If we augment then we return to the "for"
431 //cycle of runEdmonds().
433 std::queue<Node> Q; //queue of the totally unscanned nodes
436 //queue of the nodes which must be scanned for a possible shrink
438 while ( !Q.empty() ) {
441 for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
442 Node y=g.runningNode(e);
443 //growOrAugment grows if y is covered by the matching and
444 //augments if not. In this latter case it returns 1.
445 if (position[y]==C &&
446 growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
451 while ( !R.empty() ) {
455 for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
456 Node y=g.runningNode(e);
458 if ( position[y] == D && blossom.find(x) != blossom.find(y) )
459 //Recall that we have only one tree.
460 shrink( x, y, ear, blossom, rep, tree, Q);
462 while ( !Q.empty() ) {
465 for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
466 Node w=g.runningNode(f);
467 //growOrAugment grows if y is covered by the matching and
468 //augments if not. In this latter case it returns 1.
469 if (position[w]==C &&
470 growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
475 } // while ( !R.empty() )
478 void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
479 UFE& blossom, NV& rep, EFE& tree) {
480 //We have one tree, which we grow and shrink. If we augment then we
481 //return to the "for" cycle of runEdmonds().
483 std::queue<Node> Q; //queue of the unscanned nodes
485 while ( !Q.empty() ) {
490 for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
491 Node y=g.runningNode(e);
493 switch ( position[y] ) {
494 case D: //x and y must be in the same tree
495 if ( blossom.find(x) != blossom.find(y))
496 //x and y are in the same tree
497 shrink(x, y, ear, blossom, rep, tree, Q);
500 //growOrAugment grows if y is covered by the matching and
501 //augments if not. In this latter case it returns 1.
502 if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
510 void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
511 UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
512 //x and y are the two adjacent vertices in two blossoms.
514 typename Graph::template NodeMap<bool> path(g,false);
516 Node b=rep[blossom.find(x)];
519 while ( b!=INVALID ) {
520 b=rep[blossom.find(ear[b])];
523 } //we go until the root through bases of blossoms and odd vertices
526 Node middle=rep[blossom.find(top)];
528 while ( !path[middle] )
529 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
530 //Until we arrive to a node on the path, we update blossom, tree
531 //and the positions of the odd nodes.
535 middle=rep[blossom.find(top)];
537 Node blossom_base=rep[blossom.find(base)];
538 while ( middle!=blossom_base )
539 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
540 //Until we arrive to a node on the path, we update blossom, tree
541 //and the positions of the odd nodes.
543 rep[blossom.find(base)] = base;
546 void shrinkStep(Node& top, Node& middle, Node& bottom,
547 typename Graph::template NodeMap<Node>& ear,
548 UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
549 //We traverse a blossom and update everything.
553 while ( t!=middle ) {
558 bottom=_mate[middle];
559 position.set(bottom,D);
562 Node oldmiddle=middle;
563 middle=rep[blossom.find(top)];
565 tree.erase(oldmiddle);
566 blossom.insert(bottom);
567 blossom.join(bottom, oldmiddle);
568 blossom.join(top, oldmiddle);
573 bool growOrAugment(Node& y, Node& x, typename Graph::template
574 NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
575 std::queue<Node>& Q) {
576 //x is in a blossom in the tree, y is outside. If y is covered by
577 //the matching we grow, otherwise we augment. In this case we
580 if ( _mate[y]!=INVALID ) { //grow
583 rep[blossom.insert(w)] = w;
586 int t = tree.find(rep[blossom.find(x)]);
591 augment(x, ear, blossom, rep, tree);
599 void augment(Node x, typename Graph::template NodeMap<Node>& ear,
600 UFE& blossom, NV& rep, EFE& tree) {
602 while ( v!=INVALID ) {
610 int y = tree.find(rep[blossom.find(x)]);
611 for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
612 if ( position[tit] == D ) {
613 int b = blossom.find(tit);
614 for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
615 position.set(bit, C);
617 blossom.eraseClass(b);
618 } else position.set(tit, C);
626 /// \ingroup matching
628 /// \brief Weighted matching in general undirected graphs
630 /// This class provides an efficient implementation of Edmond's
631 /// maximum weighted matching algorithm. The implementation is based
632 /// on extensive use of priority queues and provides
633 /// \f$O(nm\log(n))\f$ time complexity.
635 /// The maximum weighted matching problem is to find undirected
636 /// edges in the graph with maximum overall weight and no two of
637 /// them shares their endpoints. The problem can be formulated with
638 /// the next linear program:
639 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
640 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
641 /// \f[x_e \ge 0\quad \forall e\in E\f]
642 /// \f[\max \sum_{e\in E}x_ew_e\f]
643 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
644 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
645 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
648 /// The algorithm calculates an optimal matching and a proof of the
649 /// optimality. The solution of the dual problem can be used to check
650 /// the result of the algorithm. The dual linear problem is the next:
651 /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
652 /// \f[y_u \ge 0 \quad \forall u \in V\f]
653 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
654 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
656 /// The algorithm can be executed with \c run() or the \c init() and
657 /// then the \c start() member functions. After it the matching can
658 /// be asked with \c matching() or mate() functions. The dual
659 /// solution can be get with \c nodeValue(), \c blossomNum() and \c
660 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
661 /// "BlossomIt" nested class which is able to iterate on the nodes
662 /// of a blossom. If the value type is integral then the dual
663 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
665 /// \author Balazs Dezso
666 template <typename _UGraph,
667 typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
668 class MaxWeightedMatching {
671 typedef _UGraph UGraph;
672 typedef _WeightMap WeightMap;
673 typedef typename WeightMap::Value Value;
675 /// \brief Scaling factor for dual solution
677 /// Scaling factor for dual solution, it is equal to 4 or 1
678 /// according to the value type.
679 static const int dualScale =
680 std::numeric_limits<Value>::is_integer ? 4 : 1;
682 typedef typename UGraph::template NodeMap<typename UGraph::Edge>
687 UGRAPH_TYPEDEFS(typename UGraph);
689 typedef typename UGraph::template NodeMap<Value> NodePotential;
690 typedef std::vector<Node> BlossomNodeList;
692 struct BlossomVariable {
696 BlossomVariable(int _begin, int _end, Value _value)
697 : begin(_begin), end(_end), value(_value) {}
701 typedef std::vector<BlossomVariable> BlossomPotential;
703 const UGraph& _ugraph;
704 const WeightMap& _weight;
706 MatchingMap* _matching;
708 NodePotential* _node_potential;
710 BlossomPotential _blossom_potential;
711 BlossomNodeList _blossom_node_list;
716 typedef typename UGraph::template NodeMap<int> NodeIntMap;
717 typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
718 typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
719 typedef IntegerMap<int> IntIntMap;
722 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
725 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
734 NodeIntMap *_blossom_index;
735 BlossomSet *_blossom_set;
736 IntegerMap<BlossomData>* _blossom_data;
738 NodeIntMap *_node_index;
739 EdgeIntMap *_node_heap_index;
743 NodeData(EdgeIntMap& node_heap_index)
744 : heap(node_heap_index) {}
748 BinHeap<Value, EdgeIntMap> heap;
749 std::map<int, Edge> heap_index;
754 IntegerMap<NodeData>* _node_data;
756 typedef ExtendFindEnum<IntIntMap> TreeSet;
758 IntIntMap *_tree_set_index;
761 NodeIntMap *_delta1_index;
762 BinHeap<Value, NodeIntMap> *_delta1;
764 IntIntMap *_delta2_index;
765 BinHeap<Value, IntIntMap> *_delta2;
767 UEdgeIntMap *_delta3_index;
768 BinHeap<Value, UEdgeIntMap> *_delta3;
770 IntIntMap *_delta4_index;
771 BinHeap<Value, IntIntMap> *_delta4;
775 void createStructures() {
776 _node_num = countNodes(_ugraph);
777 _blossom_num = _node_num * 3 / 2;
780 _matching = new MatchingMap(_ugraph);
782 if (!_node_potential) {
783 _node_potential = new NodePotential(_ugraph);
786 _blossom_index = new NodeIntMap(_ugraph);
787 _blossom_set = new BlossomSet(*_blossom_index);
788 _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
792 _node_index = new NodeIntMap(_ugraph);
793 _node_heap_index = new EdgeIntMap(_ugraph);
794 _node_data = new IntegerMap<NodeData>(_node_num,
795 NodeData(*_node_heap_index));
799 _tree_set_index = new IntIntMap(_blossom_num);
800 _tree_set = new TreeSet(*_tree_set_index);
803 _delta1_index = new NodeIntMap(_ugraph);
804 _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
807 _delta2_index = new IntIntMap(_blossom_num);
808 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
811 _delta3_index = new UEdgeIntMap(_ugraph);
812 _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
815 _delta4_index = new IntIntMap(_blossom_num);
816 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
820 void destroyStructures() {
821 _node_num = countNodes(_ugraph);
822 _blossom_num = _node_num * 3 / 2;
827 if (_node_potential) {
828 delete _node_potential;
831 delete _blossom_index;
833 delete _blossom_data;
838 delete _node_heap_index;
843 delete _tree_set_index;
847 delete _delta1_index;
851 delete _delta2_index;
855 delete _delta3_index;
859 delete _delta4_index;
864 void matchedToEven(int blossom, int tree) {
865 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
866 _delta2->erase(blossom);
869 if (!_blossom_set->trivial(blossom)) {
870 (*_blossom_data)[blossom].pot -=
871 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
874 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
877 _blossom_set->increase(n, std::numeric_limits<Value>::max());
878 int ni = (*_node_index)[n];
880 (*_node_data)[ni].heap.clear();
881 (*_node_data)[ni].heap_index.clear();
883 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
885 _delta1->push(n, (*_node_data)[ni].pot);
887 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
888 Node v = _ugraph.source(e);
889 int vb = _blossom_set->find(v);
890 int vi = (*_node_index)[v];
892 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
893 dualScale * _weight[e];
895 if ((*_blossom_data)[vb].status == EVEN) {
896 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
897 _delta3->push(e, rw / 2);
899 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
900 if (_delta3->state(e) != _delta3->IN_HEAP) {
901 _delta3->push(e, rw);
904 typename std::map<int, Edge>::iterator it =
905 (*_node_data)[vi].heap_index.find(tree);
907 if (it != (*_node_data)[vi].heap_index.end()) {
908 if ((*_node_data)[vi].heap[it->second] > rw) {
909 (*_node_data)[vi].heap.replace(it->second, e);
910 (*_node_data)[vi].heap.decrease(e, rw);
914 (*_node_data)[vi].heap.push(e, rw);
915 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
918 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
919 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
921 if ((*_blossom_data)[vb].status == MATCHED) {
922 if (_delta2->state(vb) != _delta2->IN_HEAP) {
923 _delta2->push(vb, _blossom_set->classPrio(vb) -
924 (*_blossom_data)[vb].offset);
925 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
926 (*_blossom_data)[vb].offset){
927 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
928 (*_blossom_data)[vb].offset);
935 (*_blossom_data)[blossom].offset = 0;
938 void matchedToOdd(int blossom) {
939 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
940 _delta2->erase(blossom);
942 (*_blossom_data)[blossom].offset += _delta_sum;
943 if (!_blossom_set->trivial(blossom)) {
944 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
945 (*_blossom_data)[blossom].offset);
949 void evenToMatched(int blossom, int tree) {
950 if (!_blossom_set->trivial(blossom)) {
951 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
954 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
956 int ni = (*_node_index)[n];
957 (*_node_data)[ni].pot -= _delta_sum;
961 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
962 Node v = _ugraph.source(e);
963 int vb = _blossom_set->find(v);
964 int vi = (*_node_index)[v];
966 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
967 dualScale * _weight[e];
970 if (_delta3->state(e) == _delta3->IN_HEAP) {
973 } else if ((*_blossom_data)[vb].status == EVEN) {
975 if (_delta3->state(e) == _delta3->IN_HEAP) {
979 int vt = _tree_set->find(vb);
983 Edge r = _ugraph.oppositeEdge(e);
985 typename std::map<int, Edge>::iterator it =
986 (*_node_data)[ni].heap_index.find(vt);
988 if (it != (*_node_data)[ni].heap_index.end()) {
989 if ((*_node_data)[ni].heap[it->second] > rw) {
990 (*_node_data)[ni].heap.replace(it->second, r);
991 (*_node_data)[ni].heap.decrease(r, rw);
995 (*_node_data)[ni].heap.push(r, rw);
996 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
999 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1000 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1002 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1003 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1004 (*_blossom_data)[blossom].offset);
1005 } else if ((*_delta2)[blossom] >
1006 _blossom_set->classPrio(blossom) -
1007 (*_blossom_data)[blossom].offset){
1008 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1009 (*_blossom_data)[blossom].offset);
1014 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1015 if (_delta3->state(e) == _delta3->IN_HEAP) {
1020 typename std::map<int, Edge>::iterator it =
1021 (*_node_data)[vi].heap_index.find(tree);
1023 if (it != (*_node_data)[vi].heap_index.end()) {
1024 (*_node_data)[vi].heap.erase(it->second);
1025 (*_node_data)[vi].heap_index.erase(it);
1026 if ((*_node_data)[vi].heap.empty()) {
1027 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1028 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1029 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1032 if ((*_blossom_data)[vb].status == MATCHED) {
1033 if (_blossom_set->classPrio(vb) ==
1034 std::numeric_limits<Value>::max()) {
1036 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1037 (*_blossom_data)[vb].offset) {
1038 _delta2->increase(vb, _blossom_set->classPrio(vb) -
1039 (*_blossom_data)[vb].offset);
1048 void oddToMatched(int blossom) {
1049 (*_blossom_data)[blossom].offset -= _delta_sum;
1051 if (_blossom_set->classPrio(blossom) !=
1052 std::numeric_limits<Value>::max()) {
1053 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1054 (*_blossom_data)[blossom].offset);
1057 if (!_blossom_set->trivial(blossom)) {
1058 _delta4->erase(blossom);
1062 void oddToEven(int blossom, int tree) {
1063 if (!_blossom_set->trivial(blossom)) {
1064 _delta4->erase(blossom);
1065 (*_blossom_data)[blossom].pot -=
1066 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1069 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1070 n != INVALID; ++n) {
1071 int ni = (*_node_index)[n];
1073 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1075 (*_node_data)[ni].heap.clear();
1076 (*_node_data)[ni].heap_index.clear();
1077 (*_node_data)[ni].pot +=
1078 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1080 _delta1->push(n, (*_node_data)[ni].pot);
1082 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1083 Node v = _ugraph.source(e);
1084 int vb = _blossom_set->find(v);
1085 int vi = (*_node_index)[v];
1087 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1088 dualScale * _weight[e];
1090 if ((*_blossom_data)[vb].status == EVEN) {
1091 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1092 _delta3->push(e, rw / 2);
1094 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1095 if (_delta3->state(e) != _delta3->IN_HEAP) {
1096 _delta3->push(e, rw);
1100 typename std::map<int, Edge>::iterator it =
1101 (*_node_data)[vi].heap_index.find(tree);
1103 if (it != (*_node_data)[vi].heap_index.end()) {
1104 if ((*_node_data)[vi].heap[it->second] > rw) {
1105 (*_node_data)[vi].heap.replace(it->second, e);
1106 (*_node_data)[vi].heap.decrease(e, rw);
1110 (*_node_data)[vi].heap.push(e, rw);
1111 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1114 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1115 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1117 if ((*_blossom_data)[vb].status == MATCHED) {
1118 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1119 _delta2->push(vb, _blossom_set->classPrio(vb) -
1120 (*_blossom_data)[vb].offset);
1121 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1122 (*_blossom_data)[vb].offset) {
1123 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1124 (*_blossom_data)[vb].offset);
1131 (*_blossom_data)[blossom].offset = 0;
1135 void matchedToUnmatched(int blossom) {
1136 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1137 _delta2->erase(blossom);
1140 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1141 n != INVALID; ++n) {
1142 int ni = (*_node_index)[n];
1144 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1146 (*_node_data)[ni].heap.clear();
1147 (*_node_data)[ni].heap_index.clear();
1149 for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1150 Node v = _ugraph.target(e);
1151 int vb = _blossom_set->find(v);
1152 int vi = (*_node_index)[v];
1154 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1155 dualScale * _weight[e];
1157 if ((*_blossom_data)[vb].status == EVEN) {
1158 if (_delta3->state(e) != _delta3->IN_HEAP) {
1159 _delta3->push(e, rw);
1166 void unmatchedToMatched(int blossom) {
1167 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1168 n != INVALID; ++n) {
1169 int ni = (*_node_index)[n];
1171 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1172 Node v = _ugraph.source(e);
1173 int vb = _blossom_set->find(v);
1174 int vi = (*_node_index)[v];
1176 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1177 dualScale * _weight[e];
1179 if (vb == blossom) {
1180 if (_delta3->state(e) == _delta3->IN_HEAP) {
1183 } else if ((*_blossom_data)[vb].status == EVEN) {
1185 if (_delta3->state(e) == _delta3->IN_HEAP) {
1189 int vt = _tree_set->find(vb);
1191 Edge r = _ugraph.oppositeEdge(e);
1193 typename std::map<int, Edge>::iterator it =
1194 (*_node_data)[ni].heap_index.find(vt);
1196 if (it != (*_node_data)[ni].heap_index.end()) {
1197 if ((*_node_data)[ni].heap[it->second] > rw) {
1198 (*_node_data)[ni].heap.replace(it->second, r);
1199 (*_node_data)[ni].heap.decrease(r, rw);
1203 (*_node_data)[ni].heap.push(r, rw);
1204 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1207 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1208 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1210 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1211 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1212 (*_blossom_data)[blossom].offset);
1213 } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1214 (*_blossom_data)[blossom].offset){
1215 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1216 (*_blossom_data)[blossom].offset);
1220 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1221 if (_delta3->state(e) == _delta3->IN_HEAP) {
1229 void alternatePath(int even, int tree) {
1232 evenToMatched(even, tree);
1233 (*_blossom_data)[even].status = MATCHED;
1235 while ((*_blossom_data)[even].pred != INVALID) {
1236 odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
1237 (*_blossom_data)[odd].status = MATCHED;
1239 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1241 even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
1242 (*_blossom_data)[even].status = MATCHED;
1243 evenToMatched(even, tree);
1244 (*_blossom_data)[even].next =
1245 _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
1250 void destroyTree(int tree) {
1251 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1252 if ((*_blossom_data)[b].status == EVEN) {
1253 (*_blossom_data)[b].status = MATCHED;
1254 evenToMatched(b, tree);
1255 } else if ((*_blossom_data)[b].status == ODD) {
1256 (*_blossom_data)[b].status = MATCHED;
1260 _tree_set->eraseClass(tree);
1264 void unmatchNode(const Node& node) {
1265 int blossom = _blossom_set->find(node);
1266 int tree = _tree_set->find(blossom);
1268 alternatePath(blossom, tree);
1271 (*_blossom_data)[blossom].status = UNMATCHED;
1272 (*_blossom_data)[blossom].base = node;
1273 matchedToUnmatched(blossom);
1277 void augmentOnEdge(const UEdge& edge) {
1279 int left = _blossom_set->find(_ugraph.source(edge));
1280 int right = _blossom_set->find(_ugraph.target(edge));
1282 if ((*_blossom_data)[left].status == EVEN) {
1283 int left_tree = _tree_set->find(left);
1284 alternatePath(left, left_tree);
1285 destroyTree(left_tree);
1287 (*_blossom_data)[left].status = MATCHED;
1288 unmatchedToMatched(left);
1291 if ((*_blossom_data)[right].status == EVEN) {
1292 int right_tree = _tree_set->find(right);
1293 alternatePath(right, right_tree);
1294 destroyTree(right_tree);
1296 (*_blossom_data)[right].status = MATCHED;
1297 unmatchedToMatched(right);
1300 (*_blossom_data)[left].next = _ugraph.direct(edge, true);
1301 (*_blossom_data)[right].next = _ugraph.direct(edge, false);
1304 void extendOnEdge(const Edge& edge) {
1305 int base = _blossom_set->find(_ugraph.target(edge));
1306 int tree = _tree_set->find(base);
1308 int odd = _blossom_set->find(_ugraph.source(edge));
1309 _tree_set->insert(odd, tree);
1310 (*_blossom_data)[odd].status = ODD;
1312 (*_blossom_data)[odd].pred = edge;
1314 int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
1315 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1316 _tree_set->insert(even, tree);
1317 (*_blossom_data)[even].status = EVEN;
1318 matchedToEven(even, tree);
1321 void shrinkOnEdge(const UEdge& uedge, int tree) {
1323 std::vector<int> left_path, right_path;
1326 std::set<int> left_set, right_set;
1327 int left = _blossom_set->find(_ugraph.source(uedge));
1328 left_path.push_back(left);
1329 left_set.insert(left);
1331 int right = _blossom_set->find(_ugraph.target(uedge));
1332 right_path.push_back(right);
1333 right_set.insert(right);
1337 if ((*_blossom_data)[left].pred == INVALID) break;
1340 _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
1341 left_path.push_back(left);
1343 _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
1344 left_path.push_back(left);
1346 left_set.insert(left);
1348 if (right_set.find(left) != right_set.end()) {
1353 if ((*_blossom_data)[right].pred == INVALID) break;
1356 _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
1357 right_path.push_back(right);
1359 _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
1360 right_path.push_back(right);
1362 right_set.insert(right);
1364 if (left_set.find(right) != left_set.end()) {
1372 if ((*_blossom_data)[left].pred == INVALID) {
1374 while (left_set.find(nca) == left_set.end()) {
1376 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1377 right_path.push_back(nca);
1379 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1380 right_path.push_back(nca);
1384 while (right_set.find(nca) == right_set.end()) {
1386 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1387 left_path.push_back(nca);
1389 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1390 left_path.push_back(nca);
1396 std::vector<int> subblossoms;
1399 prev = _ugraph.direct(uedge, true);
1400 for (int i = 0; left_path[i] != nca; i += 2) {
1401 subblossoms.push_back(left_path[i]);
1402 (*_blossom_data)[left_path[i]].next = prev;
1403 _tree_set->erase(left_path[i]);
1405 subblossoms.push_back(left_path[i + 1]);
1406 (*_blossom_data)[left_path[i + 1]].status = EVEN;
1407 oddToEven(left_path[i + 1], tree);
1408 _tree_set->erase(left_path[i + 1]);
1409 prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
1413 while (right_path[k] != nca) ++k;
1415 subblossoms.push_back(nca);
1416 (*_blossom_data)[nca].next = prev;
1418 for (int i = k - 2; i >= 0; i -= 2) {
1419 subblossoms.push_back(right_path[i + 1]);
1420 (*_blossom_data)[right_path[i + 1]].status = EVEN;
1421 oddToEven(right_path[i + 1], tree);
1422 _tree_set->erase(right_path[i + 1]);
1424 (*_blossom_data)[right_path[i + 1]].next =
1425 (*_blossom_data)[right_path[i + 1]].pred;
1427 subblossoms.push_back(right_path[i]);
1428 _tree_set->erase(right_path[i]);
1432 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1434 for (int i = 0; i < int(subblossoms.size()); ++i) {
1435 if (!_blossom_set->trivial(subblossoms[i])) {
1436 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1438 (*_blossom_data)[subblossoms[i]].status = MATCHED;
1441 (*_blossom_data)[surface].pot = -2 * _delta_sum;
1442 (*_blossom_data)[surface].offset = 0;
1443 (*_blossom_data)[surface].status = EVEN;
1444 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1445 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1447 _tree_set->insert(surface, tree);
1448 _tree_set->erase(nca);
1451 void splitBlossom(int blossom) {
1452 Edge next = (*_blossom_data)[blossom].next;
1453 Edge pred = (*_blossom_data)[blossom].pred;
1455 int tree = _tree_set->find(blossom);
1457 (*_blossom_data)[blossom].status = MATCHED;
1458 oddToMatched(blossom);
1459 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1460 _delta2->erase(blossom);
1463 std::vector<int> subblossoms;
1464 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1466 Value offset = (*_blossom_data)[blossom].offset;
1467 int b = _blossom_set->find(_ugraph.source(pred));
1468 int d = _blossom_set->find(_ugraph.source(next));
1470 int ib = -1, id = -1;
1471 for (int i = 0; i < int(subblossoms.size()); ++i) {
1472 if (subblossoms[i] == b) ib = i;
1473 if (subblossoms[i] == d) id = i;
1475 (*_blossom_data)[subblossoms[i]].offset = offset;
1476 if (!_blossom_set->trivial(subblossoms[i])) {
1477 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1479 if (_blossom_set->classPrio(subblossoms[i]) !=
1480 std::numeric_limits<Value>::max()) {
1481 _delta2->push(subblossoms[i],
1482 _blossom_set->classPrio(subblossoms[i]) -
1483 (*_blossom_data)[subblossoms[i]].offset);
1487 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1488 for (int i = (id + 1) % subblossoms.size();
1489 i != ib; i = (i + 2) % subblossoms.size()) {
1490 int sb = subblossoms[i];
1491 int tb = subblossoms[(i + 1) % subblossoms.size()];
1492 (*_blossom_data)[sb].next =
1493 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1496 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1497 int sb = subblossoms[i];
1498 int tb = subblossoms[(i + 1) % subblossoms.size()];
1499 int ub = subblossoms[(i + 2) % subblossoms.size()];
1501 (*_blossom_data)[sb].status = ODD;
1503 _tree_set->insert(sb, tree);
1504 (*_blossom_data)[sb].pred = pred;
1505 (*_blossom_data)[sb].next =
1506 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1508 pred = (*_blossom_data)[ub].next;
1510 (*_blossom_data)[tb].status = EVEN;
1511 matchedToEven(tb, tree);
1512 _tree_set->insert(tb, tree);
1513 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1516 (*_blossom_data)[subblossoms[id]].status = ODD;
1517 matchedToOdd(subblossoms[id]);
1518 _tree_set->insert(subblossoms[id], tree);
1519 (*_blossom_data)[subblossoms[id]].next = next;
1520 (*_blossom_data)[subblossoms[id]].pred = pred;
1524 for (int i = (ib + 1) % subblossoms.size();
1525 i != id; i = (i + 2) % subblossoms.size()) {
1526 int sb = subblossoms[i];
1527 int tb = subblossoms[(i + 1) % subblossoms.size()];
1528 (*_blossom_data)[sb].next =
1529 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1532 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1533 int sb = subblossoms[i];
1534 int tb = subblossoms[(i + 1) % subblossoms.size()];
1535 int ub = subblossoms[(i + 2) % subblossoms.size()];
1537 (*_blossom_data)[sb].status = ODD;
1539 _tree_set->insert(sb, tree);
1540 (*_blossom_data)[sb].next = next;
1541 (*_blossom_data)[sb].pred =
1542 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1544 (*_blossom_data)[tb].status = EVEN;
1545 matchedToEven(tb, tree);
1546 _tree_set->insert(tb, tree);
1547 (*_blossom_data)[tb].pred =
1548 (*_blossom_data)[tb].next =
1549 _ugraph.oppositeEdge((*_blossom_data)[ub].next);
1550 next = (*_blossom_data)[ub].next;
1553 (*_blossom_data)[subblossoms[ib]].status = ODD;
1554 matchedToOdd(subblossoms[ib]);
1555 _tree_set->insert(subblossoms[ib], tree);
1556 (*_blossom_data)[subblossoms[ib]].next = next;
1557 (*_blossom_data)[subblossoms[ib]].pred = pred;
1559 _tree_set->erase(blossom);
1562 void extractBlossom(int blossom, const Node& base, const Edge& matching) {
1563 if (_blossom_set->trivial(blossom)) {
1564 int bi = (*_node_index)[base];
1565 Value pot = (*_node_data)[bi].pot;
1567 _matching->set(base, matching);
1568 _blossom_node_list.push_back(base);
1569 _node_potential->set(base, pot);
1572 Value pot = (*_blossom_data)[blossom].pot;
1573 int bn = _blossom_node_list.size();
1575 std::vector<int> subblossoms;
1576 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1577 int b = _blossom_set->find(base);
1579 for (int i = 0; i < int(subblossoms.size()); ++i) {
1580 if (subblossoms[i] == b) { ib = i; break; }
1583 for (int i = 1; i < int(subblossoms.size()); i += 2) {
1584 int sb = subblossoms[(ib + i) % subblossoms.size()];
1585 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1587 Edge m = (*_blossom_data)[tb].next;
1588 extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
1589 extractBlossom(tb, _ugraph.source(m), m);
1591 extractBlossom(subblossoms[ib], base, matching);
1593 int en = _blossom_node_list.size();
1595 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1599 void extractMatching() {
1600 std::vector<int> blossoms;
1601 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1602 blossoms.push_back(c);
1605 for (int i = 0; i < int(blossoms.size()); ++i) {
1606 if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1608 Value offset = (*_blossom_data)[blossoms[i]].offset;
1609 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1610 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1611 n != INVALID; ++n) {
1612 (*_node_data)[(*_node_index)[n]].pot -= offset;
1615 Edge matching = (*_blossom_data)[blossoms[i]].next;
1616 Node base = _ugraph.source(matching);
1617 extractBlossom(blossoms[i], base, matching);
1619 Node base = (*_blossom_data)[blossoms[i]].base;
1620 extractBlossom(blossoms[i], base, INVALID);
1627 /// \brief Constructor
1630 MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight)
1631 : _ugraph(ugraph), _weight(weight), _matching(0),
1632 _node_potential(0), _blossom_potential(), _blossom_node_list(),
1633 _node_num(0), _blossom_num(0),
1635 _blossom_index(0), _blossom_set(0), _blossom_data(0),
1636 _node_index(0), _node_heap_index(0), _node_data(0),
1637 _tree_set_index(0), _tree_set(0),
1639 _delta1_index(0), _delta1(0),
1640 _delta2_index(0), _delta2(0),
1641 _delta3_index(0), _delta3(0),
1642 _delta4_index(0), _delta4(0),
1646 ~MaxWeightedMatching() {
1647 destroyStructures();
1650 /// \name Execution control
1651 /// The simplest way to execute the algorithm is to use the member
1652 /// \c run() member function.
1656 /// \brief Initialize the algorithm
1658 /// Initialize the algorithm
1662 for (EdgeIt e(_ugraph); e != INVALID; ++e) {
1663 _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
1665 for (NodeIt n(_ugraph); n != INVALID; ++n) {
1666 _delta1_index->set(n, _delta1->PRE_HEAP);
1668 for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
1669 _delta3_index->set(e, _delta3->PRE_HEAP);
1671 for (int i = 0; i < _blossom_num; ++i) {
1672 _delta2_index->set(i, _delta2->PRE_HEAP);
1673 _delta4_index->set(i, _delta4->PRE_HEAP);
1677 for (NodeIt n(_ugraph); n != INVALID; ++n) {
1679 for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1680 if (_ugraph.target(e) == n) continue;
1681 if ((dualScale * _weight[e]) / 2 > max) {
1682 max = (dualScale * _weight[e]) / 2;
1685 _node_index->set(n, index);
1686 (*_node_data)[index].pot = max;
1687 _delta1->push(n, max);
1689 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1691 _tree_set->insert(blossom);
1693 (*_blossom_data)[blossom].status = EVEN;
1694 (*_blossom_data)[blossom].pred = INVALID;
1695 (*_blossom_data)[blossom].next = INVALID;
1696 (*_blossom_data)[blossom].pot = 0;
1697 (*_blossom_data)[blossom].offset = 0;
1700 for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
1701 int si = (*_node_index)[_ugraph.source(e)];
1702 int ti = (*_node_index)[_ugraph.target(e)];
1703 if (_ugraph.source(e) != _ugraph.target(e)) {
1704 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1705 dualScale * _weight[e]) / 2);
1710 /// \brief Starts the algorithm
1712 /// Starts the algorithm
1718 int unmatched = _node_num;
1719 while (unmatched > 0) {
1720 Value d1 = !_delta1->empty() ?
1721 _delta1->prio() : std::numeric_limits<Value>::max();
1723 Value d2 = !_delta2->empty() ?
1724 _delta2->prio() : std::numeric_limits<Value>::max();
1726 Value d3 = !_delta3->empty() ?
1727 _delta3->prio() : std::numeric_limits<Value>::max();
1729 Value d4 = !_delta4->empty() ?
1730 _delta4->prio() : std::numeric_limits<Value>::max();
1732 _delta_sum = d1; OpType ot = D1;
1733 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1734 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1735 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1741 Node n = _delta1->top();
1748 int blossom = _delta2->top();
1749 Node n = _blossom_set->classTop(blossom);
1750 Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
1756 UEdge e = _delta3->top();
1758 int left_blossom = _blossom_set->find(_ugraph.source(e));
1759 int right_blossom = _blossom_set->find(_ugraph.target(e));
1761 if (left_blossom == right_blossom) {
1765 if ((*_blossom_data)[left_blossom].status == EVEN) {
1766 left_tree = _tree_set->find(left_blossom);
1772 if ((*_blossom_data)[right_blossom].status == EVEN) {
1773 right_tree = _tree_set->find(right_blossom);
1779 if (left_tree == right_tree) {
1780 shrinkOnEdge(e, left_tree);
1788 splitBlossom(_delta4->top());
1795 /// \brief Runs %MaxWeightedMatching algorithm.
1797 /// This method runs the %MaxWeightedMatching algorithm.
1799 /// \note mwm.run() is just a shortcut of the following code.
1811 /// \name Primal solution
1812 /// Functions for get the primal solution, ie. the matching.
1816 /// \brief Returns the matching value.
1818 /// Returns the matching value.
1819 Value matchingValue() const {
1821 for (NodeIt n(_ugraph); n != INVALID; ++n) {
1822 if ((*_matching)[n] != INVALID) {
1823 sum += _weight[(*_matching)[n]];
1829 /// \brief Returns true when the edge is in the matching.
1831 /// Returns true when the edge is in the matching.
1832 bool matching(const UEdge& edge) const {
1833 return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
1836 /// \brief Returns the incident matching edge.
1838 /// Returns the incident matching edge from given node. If the
1839 /// node is not matched then it gives back \c INVALID.
1840 Edge matching(const Node& node) const {
1841 return (*_matching)[node];
1844 /// \brief Returns the mate of the node.
1846 /// Returns the adjancent node in a mathcing edge. If the node is
1847 /// not matched then it gives back \c INVALID.
1848 Node mate(const Node& node) const {
1849 return (*_matching)[node] != INVALID ?
1850 _ugraph.target((*_matching)[node]) : INVALID;
1855 /// \name Dual solution
1856 /// Functions for get the dual solution.
1860 /// \brief Returns the value of the dual solution.
1862 /// Returns the value of the dual solution. It should be equal to
1863 /// the primal value scaled by \ref dualScale "dual scale".
1864 Value dualValue() const {
1866 for (NodeIt n(_ugraph); n != INVALID; ++n) {
1867 sum += nodeValue(n);
1869 for (int i = 0; i < blossomNum(); ++i) {
1870 sum += blossomValue(i) * (blossomSize(i) / 2);
1875 /// \brief Returns the value of the node.
1877 /// Returns the the value of the node.
1878 Value nodeValue(const Node& n) const {
1879 return (*_node_potential)[n];
1882 /// \brief Returns the number of the blossoms in the basis.
1884 /// Returns the number of the blossoms in the basis.
1886 int blossomNum() const {
1887 return _blossom_potential.size();
1891 /// \brief Returns the number of the nodes in the blossom.
1893 /// Returns the number of the nodes in the blossom.
1894 int blossomSize(int k) const {
1895 return _blossom_potential[k].end - _blossom_potential[k].begin;
1898 /// \brief Returns the value of the blossom.
1900 /// Returns the the value of the blossom.
1902 Value blossomValue(int k) const {
1903 return _blossom_potential[k].value;
1906 /// \brief Lemon iterator for get the items of the blossom.
1908 /// Lemon iterator for get the nodes of the blossom. This class
1909 /// provides a common style lemon iterator which gives back a
1910 /// subset of the nodes.
1914 /// \brief Constructor.
1916 /// Constructor for get the nodes of the variable.
1917 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1918 : _algorithm(&algorithm)
1920 _index = _algorithm->_blossom_potential[variable].begin;
1921 _last = _algorithm->_blossom_potential[variable].end;
1924 /// \brief Invalid constructor.
1926 /// Invalid constructor.
1927 BlossomIt(Invalid) : _index(-1) {}
1929 /// \brief Conversion to node.
1931 /// Conversion to node.
1932 operator Node() const {
1933 return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1936 /// \brief Increment operator.
1938 /// Increment operator.
1939 BlossomIt& operator++() {
1941 if (_index == _last) {
1947 bool operator==(const BlossomIt& it) const {
1948 return _index == it._index;
1950 bool operator!=(const BlossomIt& it) const {
1951 return _index != it._index;
1955 const MaxWeightedMatching* _algorithm;
1964 /// \ingroup matching
1966 /// \brief Weighted perfect matching in general undirected graphs
1968 /// This class provides an efficient implementation of Edmond's
1969 /// maximum weighted perfecr matching algorithm. The implementation
1970 /// is based on extensive use of priority queues and provides
1971 /// \f$O(nm\log(n))\f$ time complexity.
1973 /// The maximum weighted matching problem is to find undirected
1974 /// edges in the graph with maximum overall weight and no two of
1975 /// them shares their endpoints and covers all nodes. The problem
1976 /// can be formulated with the next linear program:
1977 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1978 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1979 /// \f[x_e \ge 0\quad \forall e\in E\f]
1980 /// \f[\max \sum_{e\in E}x_ew_e\f]
1981 /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1982 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
1983 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1986 /// The algorithm calculates an optimal matching and a proof of the
1987 /// optimality. The solution of the dual problem can be used to check
1988 /// the result of the algorithm. The dual linear problem is the next:
1989 /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
1990 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1991 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1993 /// The algorithm can be executed with \c run() or the \c init() and
1994 /// then the \c start() member functions. After it the matching can
1995 /// be asked with \c matching() or mate() functions. The dual
1996 /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1997 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1998 /// "BlossomIt" nested class which is able to iterate on the nodes
1999 /// of a blossom. If the value type is integral then the dual
2000 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
2002 /// \author Balazs Dezso
2003 template <typename _UGraph,
2004 typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
2005 class MaxWeightedPerfectMatching {
2008 typedef _UGraph UGraph;
2009 typedef _WeightMap WeightMap;
2010 typedef typename WeightMap::Value Value;
2012 /// \brief Scaling factor for dual solution
2014 /// Scaling factor for dual solution, it is equal to 4 or 1
2015 /// according to the value type.
2016 static const int dualScale =
2017 std::numeric_limits<Value>::is_integer ? 4 : 1;
2019 typedef typename UGraph::template NodeMap<typename UGraph::Edge>
2024 UGRAPH_TYPEDEFS(typename UGraph);
2026 typedef typename UGraph::template NodeMap<Value> NodePotential;
2027 typedef std::vector<Node> BlossomNodeList;
2029 struct BlossomVariable {
2033 BlossomVariable(int _begin, int _end, Value _value)
2034 : begin(_begin), end(_end), value(_value) {}
2038 typedef std::vector<BlossomVariable> BlossomPotential;
2040 const UGraph& _ugraph;
2041 const WeightMap& _weight;
2043 MatchingMap* _matching;
2045 NodePotential* _node_potential;
2047 BlossomPotential _blossom_potential;
2048 BlossomNodeList _blossom_node_list;
2053 typedef typename UGraph::template NodeMap<int> NodeIntMap;
2054 typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
2055 typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
2056 typedef IntegerMap<int> IntIntMap;
2059 EVEN = -1, MATCHED = 0, ODD = 1
2062 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
2063 struct BlossomData {
2070 NodeIntMap *_blossom_index;
2071 BlossomSet *_blossom_set;
2072 IntegerMap<BlossomData>* _blossom_data;
2074 NodeIntMap *_node_index;
2075 EdgeIntMap *_node_heap_index;
2079 NodeData(EdgeIntMap& node_heap_index)
2080 : heap(node_heap_index) {}
2084 BinHeap<Value, EdgeIntMap> heap;
2085 std::map<int, Edge> heap_index;
2090 IntegerMap<NodeData>* _node_data;
2092 typedef ExtendFindEnum<IntIntMap> TreeSet;
2094 IntIntMap *_tree_set_index;
2097 IntIntMap *_delta2_index;
2098 BinHeap<Value, IntIntMap> *_delta2;
2100 UEdgeIntMap *_delta3_index;
2101 BinHeap<Value, UEdgeIntMap> *_delta3;
2103 IntIntMap *_delta4_index;
2104 BinHeap<Value, IntIntMap> *_delta4;
2108 void createStructures() {
2109 _node_num = countNodes(_ugraph);
2110 _blossom_num = _node_num * 3 / 2;
2113 _matching = new MatchingMap(_ugraph);
2115 if (!_node_potential) {
2116 _node_potential = new NodePotential(_ugraph);
2118 if (!_blossom_set) {
2119 _blossom_index = new NodeIntMap(_ugraph);
2120 _blossom_set = new BlossomSet(*_blossom_index);
2121 _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
2125 _node_index = new NodeIntMap(_ugraph);
2126 _node_heap_index = new EdgeIntMap(_ugraph);
2127 _node_data = new IntegerMap<NodeData>(_node_num,
2128 NodeData(*_node_heap_index));
2132 _tree_set_index = new IntIntMap(_blossom_num);
2133 _tree_set = new TreeSet(*_tree_set_index);
2136 _delta2_index = new IntIntMap(_blossom_num);
2137 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2140 _delta3_index = new UEdgeIntMap(_ugraph);
2141 _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
2144 _delta4_index = new IntIntMap(_blossom_num);
2145 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2149 void destroyStructures() {
2150 _node_num = countNodes(_ugraph);
2151 _blossom_num = _node_num * 3 / 2;
2156 if (_node_potential) {
2157 delete _node_potential;
2160 delete _blossom_index;
2161 delete _blossom_set;
2162 delete _blossom_data;
2167 delete _node_heap_index;
2172 delete _tree_set_index;
2176 delete _delta2_index;
2180 delete _delta3_index;
2184 delete _delta4_index;
2189 void matchedToEven(int blossom, int tree) {
2190 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2191 _delta2->erase(blossom);
2194 if (!_blossom_set->trivial(blossom)) {
2195 (*_blossom_data)[blossom].pot -=
2196 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2199 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2200 n != INVALID; ++n) {
2202 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2203 int ni = (*_node_index)[n];
2205 (*_node_data)[ni].heap.clear();
2206 (*_node_data)[ni].heap_index.clear();
2208 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2210 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
2211 Node v = _ugraph.source(e);
2212 int vb = _blossom_set->find(v);
2213 int vi = (*_node_index)[v];
2215 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2216 dualScale * _weight[e];
2218 if ((*_blossom_data)[vb].status == EVEN) {
2219 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2220 _delta3->push(e, rw / 2);
2223 typename std::map<int, Edge>::iterator it =
2224 (*_node_data)[vi].heap_index.find(tree);
2226 if (it != (*_node_data)[vi].heap_index.end()) {
2227 if ((*_node_data)[vi].heap[it->second] > rw) {
2228 (*_node_data)[vi].heap.replace(it->second, e);
2229 (*_node_data)[vi].heap.decrease(e, rw);
2233 (*_node_data)[vi].heap.push(e, rw);
2234 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2237 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2238 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2240 if ((*_blossom_data)[vb].status == MATCHED) {
2241 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2242 _delta2->push(vb, _blossom_set->classPrio(vb) -
2243 (*_blossom_data)[vb].offset);
2244 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2245 (*_blossom_data)[vb].offset){
2246 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2247 (*_blossom_data)[vb].offset);
2254 (*_blossom_data)[blossom].offset = 0;
2257 void matchedToOdd(int blossom) {
2258 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2259 _delta2->erase(blossom);
2261 (*_blossom_data)[blossom].offset += _delta_sum;
2262 if (!_blossom_set->trivial(blossom)) {
2263 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2264 (*_blossom_data)[blossom].offset);
2268 void evenToMatched(int blossom, int tree) {
2269 if (!_blossom_set->trivial(blossom)) {
2270 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2273 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2274 n != INVALID; ++n) {
2275 int ni = (*_node_index)[n];
2276 (*_node_data)[ni].pot -= _delta_sum;
2278 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
2279 Node v = _ugraph.source(e);
2280 int vb = _blossom_set->find(v);
2281 int vi = (*_node_index)[v];
2283 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2284 dualScale * _weight[e];
2286 if (vb == blossom) {
2287 if (_delta3->state(e) == _delta3->IN_HEAP) {
2290 } else if ((*_blossom_data)[vb].status == EVEN) {
2292 if (_delta3->state(e) == _delta3->IN_HEAP) {
2296 int vt = _tree_set->find(vb);
2300 Edge r = _ugraph.oppositeEdge(e);
2302 typename std::map<int, Edge>::iterator it =
2303 (*_node_data)[ni].heap_index.find(vt);
2305 if (it != (*_node_data)[ni].heap_index.end()) {
2306 if ((*_node_data)[ni].heap[it->second] > rw) {
2307 (*_node_data)[ni].heap.replace(it->second, r);
2308 (*_node_data)[ni].heap.decrease(r, rw);
2312 (*_node_data)[ni].heap.push(r, rw);
2313 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2316 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2317 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2319 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2320 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2321 (*_blossom_data)[blossom].offset);
2322 } else if ((*_delta2)[blossom] >
2323 _blossom_set->classPrio(blossom) -
2324 (*_blossom_data)[blossom].offset){
2325 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2326 (*_blossom_data)[blossom].offset);
2332 typename std::map<int, Edge>::iterator it =
2333 (*_node_data)[vi].heap_index.find(tree);
2335 if (it != (*_node_data)[vi].heap_index.end()) {
2336 (*_node_data)[vi].heap.erase(it->second);
2337 (*_node_data)[vi].heap_index.erase(it);
2338 if ((*_node_data)[vi].heap.empty()) {
2339 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2340 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2341 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2344 if ((*_blossom_data)[vb].status == MATCHED) {
2345 if (_blossom_set->classPrio(vb) ==
2346 std::numeric_limits<Value>::max()) {
2348 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2349 (*_blossom_data)[vb].offset) {
2350 _delta2->increase(vb, _blossom_set->classPrio(vb) -
2351 (*_blossom_data)[vb].offset);
2360 void oddToMatched(int blossom) {
2361 (*_blossom_data)[blossom].offset -= _delta_sum;
2363 if (_blossom_set->classPrio(blossom) !=
2364 std::numeric_limits<Value>::max()) {
2365 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2366 (*_blossom_data)[blossom].offset);
2369 if (!_blossom_set->trivial(blossom)) {
2370 _delta4->erase(blossom);
2374 void oddToEven(int blossom, int tree) {
2375 if (!_blossom_set->trivial(blossom)) {
2376 _delta4->erase(blossom);
2377 (*_blossom_data)[blossom].pot -=
2378 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2381 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2382 n != INVALID; ++n) {
2383 int ni = (*_node_index)[n];
2385 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2387 (*_node_data)[ni].heap.clear();
2388 (*_node_data)[ni].heap_index.clear();
2389 (*_node_data)[ni].pot +=
2390 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2392 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
2393 Node v = _ugraph.source(e);
2394 int vb = _blossom_set->find(v);
2395 int vi = (*_node_index)[v];
2397 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2398 dualScale * _weight[e];
2400 if ((*_blossom_data)[vb].status == EVEN) {
2401 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2402 _delta3->push(e, rw / 2);
2406 typename std::map<int, Edge>::iterator it =
2407 (*_node_data)[vi].heap_index.find(tree);
2409 if (it != (*_node_data)[vi].heap_index.end()) {
2410 if ((*_node_data)[vi].heap[it->second] > rw) {
2411 (*_node_data)[vi].heap.replace(it->second, e);
2412 (*_node_data)[vi].heap.decrease(e, rw);
2416 (*_node_data)[vi].heap.push(e, rw);
2417 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2420 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2421 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2423 if ((*_blossom_data)[vb].status == MATCHED) {
2424 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2425 _delta2->push(vb, _blossom_set->classPrio(vb) -
2426 (*_blossom_data)[vb].offset);
2427 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2428 (*_blossom_data)[vb].offset) {
2429 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2430 (*_blossom_data)[vb].offset);
2437 (*_blossom_data)[blossom].offset = 0;
2440 void alternatePath(int even, int tree) {
2443 evenToMatched(even, tree);
2444 (*_blossom_data)[even].status = MATCHED;
2446 while ((*_blossom_data)[even].pred != INVALID) {
2447 odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
2448 (*_blossom_data)[odd].status = MATCHED;
2450 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2452 even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
2453 (*_blossom_data)[even].status = MATCHED;
2454 evenToMatched(even, tree);
2455 (*_blossom_data)[even].next =
2456 _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
2461 void destroyTree(int tree) {
2462 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2463 if ((*_blossom_data)[b].status == EVEN) {
2464 (*_blossom_data)[b].status = MATCHED;
2465 evenToMatched(b, tree);
2466 } else if ((*_blossom_data)[b].status == ODD) {
2467 (*_blossom_data)[b].status = MATCHED;
2471 _tree_set->eraseClass(tree);
2474 void augmentOnEdge(const UEdge& edge) {
2476 int left = _blossom_set->find(_ugraph.source(edge));
2477 int right = _blossom_set->find(_ugraph.target(edge));
2479 int left_tree = _tree_set->find(left);
2480 alternatePath(left, left_tree);
2481 destroyTree(left_tree);
2483 int right_tree = _tree_set->find(right);
2484 alternatePath(right, right_tree);
2485 destroyTree(right_tree);
2487 (*_blossom_data)[left].next = _ugraph.direct(edge, true);
2488 (*_blossom_data)[right].next = _ugraph.direct(edge, false);
2491 void extendOnEdge(const Edge& edge) {
2492 int base = _blossom_set->find(_ugraph.target(edge));
2493 int tree = _tree_set->find(base);
2495 int odd = _blossom_set->find(_ugraph.source(edge));
2496 _tree_set->insert(odd, tree);
2497 (*_blossom_data)[odd].status = ODD;
2499 (*_blossom_data)[odd].pred = edge;
2501 int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
2502 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2503 _tree_set->insert(even, tree);
2504 (*_blossom_data)[even].status = EVEN;
2505 matchedToEven(even, tree);
2508 void shrinkOnEdge(const UEdge& uedge, int tree) {
2510 std::vector<int> left_path, right_path;
2513 std::set<int> left_set, right_set;
2514 int left = _blossom_set->find(_ugraph.source(uedge));
2515 left_path.push_back(left);
2516 left_set.insert(left);
2518 int right = _blossom_set->find(_ugraph.target(uedge));
2519 right_path.push_back(right);
2520 right_set.insert(right);
2524 if ((*_blossom_data)[left].pred == INVALID) break;
2527 _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
2528 left_path.push_back(left);
2530 _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
2531 left_path.push_back(left);
2533 left_set.insert(left);
2535 if (right_set.find(left) != right_set.end()) {
2540 if ((*_blossom_data)[right].pred == INVALID) break;
2543 _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
2544 right_path.push_back(right);
2546 _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
2547 right_path.push_back(right);
2549 right_set.insert(right);
2551 if (left_set.find(right) != left_set.end()) {
2559 if ((*_blossom_data)[left].pred == INVALID) {
2561 while (left_set.find(nca) == left_set.end()) {
2563 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
2564 right_path.push_back(nca);
2566 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
2567 right_path.push_back(nca);
2571 while (right_set.find(nca) == right_set.end()) {
2573 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
2574 left_path.push_back(nca);
2576 _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
2577 left_path.push_back(nca);
2583 std::vector<int> subblossoms;
2586 prev = _ugraph.direct(uedge, true);
2587 for (int i = 0; left_path[i] != nca; i += 2) {
2588 subblossoms.push_back(left_path[i]);
2589 (*_blossom_data)[left_path[i]].next = prev;
2590 _tree_set->erase(left_path[i]);
2592 subblossoms.push_back(left_path[i + 1]);
2593 (*_blossom_data)[left_path[i + 1]].status = EVEN;
2594 oddToEven(left_path[i + 1], tree);
2595 _tree_set->erase(left_path[i + 1]);
2596 prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
2600 while (right_path[k] != nca) ++k;
2602 subblossoms.push_back(nca);
2603 (*_blossom_data)[nca].next = prev;
2605 for (int i = k - 2; i >= 0; i -= 2) {
2606 subblossoms.push_back(right_path[i + 1]);
2607 (*_blossom_data)[right_path[i + 1]].status = EVEN;
2608 oddToEven(right_path[i + 1], tree);
2609 _tree_set->erase(right_path[i + 1]);
2611 (*_blossom_data)[right_path[i + 1]].next =
2612 (*_blossom_data)[right_path[i + 1]].pred;
2614 subblossoms.push_back(right_path[i]);
2615 _tree_set->erase(right_path[i]);
2619 _blossom_set->join(subblossoms.begin(), subblossoms.end());
2621 for (int i = 0; i < int(subblossoms.size()); ++i) {
2622 if (!_blossom_set->trivial(subblossoms[i])) {
2623 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2625 (*_blossom_data)[subblossoms[i]].status = MATCHED;
2628 (*_blossom_data)[surface].pot = -2 * _delta_sum;
2629 (*_blossom_data)[surface].offset = 0;
2630 (*_blossom_data)[surface].status = EVEN;
2631 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2632 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2634 _tree_set->insert(surface, tree);
2635 _tree_set->erase(nca);
2638 void splitBlossom(int blossom) {
2639 Edge next = (*_blossom_data)[blossom].next;
2640 Edge pred = (*_blossom_data)[blossom].pred;
2642 int tree = _tree_set->find(blossom);
2644 (*_blossom_data)[blossom].status = MATCHED;
2645 oddToMatched(blossom);
2646 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2647 _delta2->erase(blossom);
2650 std::vector<int> subblossoms;
2651 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2653 Value offset = (*_blossom_data)[blossom].offset;
2654 int b = _blossom_set->find(_ugraph.source(pred));
2655 int d = _blossom_set->find(_ugraph.source(next));
2657 int ib = -1, id = -1;
2658 for (int i = 0; i < int(subblossoms.size()); ++i) {
2659 if (subblossoms[i] == b) ib = i;
2660 if (subblossoms[i] == d) id = i;
2662 (*_blossom_data)[subblossoms[i]].offset = offset;
2663 if (!_blossom_set->trivial(subblossoms[i])) {
2664 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2666 if (_blossom_set->classPrio(subblossoms[i]) !=
2667 std::numeric_limits<Value>::max()) {
2668 _delta2->push(subblossoms[i],
2669 _blossom_set->classPrio(subblossoms[i]) -
2670 (*_blossom_data)[subblossoms[i]].offset);
2674 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2675 for (int i = (id + 1) % subblossoms.size();
2676 i != ib; i = (i + 2) % subblossoms.size()) {
2677 int sb = subblossoms[i];
2678 int tb = subblossoms[(i + 1) % subblossoms.size()];
2679 (*_blossom_data)[sb].next =
2680 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
2683 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2684 int sb = subblossoms[i];
2685 int tb = subblossoms[(i + 1) % subblossoms.size()];
2686 int ub = subblossoms[(i + 2) % subblossoms.size()];
2688 (*_blossom_data)[sb].status = ODD;
2690 _tree_set->insert(sb, tree);
2691 (*_blossom_data)[sb].pred = pred;
2692 (*_blossom_data)[sb].next =
2693 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
2695 pred = (*_blossom_data)[ub].next;
2697 (*_blossom_data)[tb].status = EVEN;
2698 matchedToEven(tb, tree);
2699 _tree_set->insert(tb, tree);
2700 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2703 (*_blossom_data)[subblossoms[id]].status = ODD;
2704 matchedToOdd(subblossoms[id]);
2705 _tree_set->insert(subblossoms[id], tree);
2706 (*_blossom_data)[subblossoms[id]].next = next;
2707 (*_blossom_data)[subblossoms[id]].pred = pred;
2711 for (int i = (ib + 1) % subblossoms.size();
2712 i != id; i = (i + 2) % subblossoms.size()) {
2713 int sb = subblossoms[i];
2714 int tb = subblossoms[(i + 1) % subblossoms.size()];
2715 (*_blossom_data)[sb].next =
2716 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
2719 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2720 int sb = subblossoms[i];
2721 int tb = subblossoms[(i + 1) % subblossoms.size()];
2722 int ub = subblossoms[(i + 2) % subblossoms.size()];
2724 (*_blossom_data)[sb].status = ODD;
2726 _tree_set->insert(sb, tree);
2727 (*_blossom_data)[sb].next = next;
2728 (*_blossom_data)[sb].pred =
2729 _ugraph.oppositeEdge((*_blossom_data)[tb].next);
2731 (*_blossom_data)[tb].status = EVEN;
2732 matchedToEven(tb, tree);
2733 _tree_set->insert(tb, tree);
2734 (*_blossom_data)[tb].pred =
2735 (*_blossom_data)[tb].next =
2736 _ugraph.oppositeEdge((*_blossom_data)[ub].next);
2737 next = (*_blossom_data)[ub].next;
2740 (*_blossom_data)[subblossoms[ib]].status = ODD;
2741 matchedToOdd(subblossoms[ib]);
2742 _tree_set->insert(subblossoms[ib], tree);
2743 (*_blossom_data)[subblossoms[ib]].next = next;
2744 (*_blossom_data)[subblossoms[ib]].pred = pred;
2746 _tree_set->erase(blossom);
2749 void extractBlossom(int blossom, const Node& base, const Edge& matching) {
2750 if (_blossom_set->trivial(blossom)) {
2751 int bi = (*_node_index)[base];
2752 Value pot = (*_node_data)[bi].pot;
2754 _matching->set(base, matching);
2755 _blossom_node_list.push_back(base);
2756 _node_potential->set(base, pot);
2759 Value pot = (*_blossom_data)[blossom].pot;
2760 int bn = _blossom_node_list.size();
2762 std::vector<int> subblossoms;
2763 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2764 int b = _blossom_set->find(base);
2766 for (int i = 0; i < int(subblossoms.size()); ++i) {
2767 if (subblossoms[i] == b) { ib = i; break; }
2770 for (int i = 1; i < int(subblossoms.size()); i += 2) {
2771 int sb = subblossoms[(ib + i) % subblossoms.size()];
2772 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2774 Edge m = (*_blossom_data)[tb].next;
2775 extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
2776 extractBlossom(tb, _ugraph.source(m), m);
2778 extractBlossom(subblossoms[ib], base, matching);
2780 int en = _blossom_node_list.size();
2782 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2786 void extractMatching() {
2787 std::vector<int> blossoms;
2788 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2789 blossoms.push_back(c);
2792 for (int i = 0; i < int(blossoms.size()); ++i) {
2794 Value offset = (*_blossom_data)[blossoms[i]].offset;
2795 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2796 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2797 n != INVALID; ++n) {
2798 (*_node_data)[(*_node_index)[n]].pot -= offset;
2801 Edge matching = (*_blossom_data)[blossoms[i]].next;
2802 Node base = _ugraph.source(matching);
2803 extractBlossom(blossoms[i], base, matching);
2809 /// \brief Constructor
2812 MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight)
2813 : _ugraph(ugraph), _weight(weight), _matching(0),
2814 _node_potential(0), _blossom_potential(), _blossom_node_list(),
2815 _node_num(0), _blossom_num(0),
2817 _blossom_index(0), _blossom_set(0), _blossom_data(0),
2818 _node_index(0), _node_heap_index(0), _node_data(0),
2819 _tree_set_index(0), _tree_set(0),
2821 _delta2_index(0), _delta2(0),
2822 _delta3_index(0), _delta3(0),
2823 _delta4_index(0), _delta4(0),
2827 ~MaxWeightedPerfectMatching() {
2828 destroyStructures();
2831 /// \name Execution control
2832 /// The simplest way to execute the algorithm is to use the member
2833 /// \c run() member function.
2837 /// \brief Initialize the algorithm
2839 /// Initialize the algorithm
2843 for (EdgeIt e(_ugraph); e != INVALID; ++e) {
2844 _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
2846 for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
2847 _delta3_index->set(e, _delta3->PRE_HEAP);
2849 for (int i = 0; i < _blossom_num; ++i) {
2850 _delta2_index->set(i, _delta2->PRE_HEAP);
2851 _delta4_index->set(i, _delta4->PRE_HEAP);
2855 for (NodeIt n(_ugraph); n != INVALID; ++n) {
2856 Value max = std::numeric_limits<Value>::min();
2857 for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
2858 if (_ugraph.target(e) == n) continue;
2859 if ((dualScale * _weight[e]) / 2 > max) {
2860 max = (dualScale * _weight[e]) / 2;
2863 _node_index->set(n, index);
2864 (*_node_data)[index].pot = max;
2866 _blossom_set->insert(n, std::numeric_limits<Value>::max());
2868 _tree_set->insert(blossom);
2870 (*_blossom_data)[blossom].status = EVEN;
2871 (*_blossom_data)[blossom].pred = INVALID;
2872 (*_blossom_data)[blossom].next = INVALID;
2873 (*_blossom_data)[blossom].pot = 0;
2874 (*_blossom_data)[blossom].offset = 0;
2877 for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
2878 int si = (*_node_index)[_ugraph.source(e)];
2879 int ti = (*_node_index)[_ugraph.target(e)];
2880 if (_ugraph.source(e) != _ugraph.target(e)) {
2881 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2882 dualScale * _weight[e]) / 2);
2887 /// \brief Starts the algorithm
2889 /// Starts the algorithm
2895 int unmatched = _node_num;
2896 while (unmatched > 0) {
2897 Value d2 = !_delta2->empty() ?
2898 _delta2->prio() : std::numeric_limits<Value>::max();
2900 Value d3 = !_delta3->empty() ?
2901 _delta3->prio() : std::numeric_limits<Value>::max();
2903 Value d4 = !_delta4->empty() ?
2904 _delta4->prio() : std::numeric_limits<Value>::max();
2906 _delta_sum = d2; OpType ot = D2;
2907 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2908 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2910 if (_delta_sum == std::numeric_limits<Value>::max()) {
2917 int blossom = _delta2->top();
2918 Node n = _blossom_set->classTop(blossom);
2919 Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
2925 UEdge e = _delta3->top();
2927 int left_blossom = _blossom_set->find(_ugraph.source(e));
2928 int right_blossom = _blossom_set->find(_ugraph.target(e));
2930 if (left_blossom == right_blossom) {
2933 int left_tree = _tree_set->find(left_blossom);
2934 int right_tree = _tree_set->find(right_blossom);
2936 if (left_tree == right_tree) {
2937 shrinkOnEdge(e, left_tree);
2945 splitBlossom(_delta4->top());
2953 /// \brief Runs %MaxWeightedPerfectMatching algorithm.
2955 /// This method runs the %MaxWeightedPerfectMatching algorithm.
2957 /// \note mwm.run() is just a shortcut of the following code.
2969 /// \name Primal solution
2970 /// Functions for get the primal solution, ie. the matching.
2974 /// \brief Returns the matching value.
2976 /// Returns the matching value.
2977 Value matchingValue() const {
2979 for (NodeIt n(_ugraph); n != INVALID; ++n) {
2980 if ((*_matching)[n] != INVALID) {
2981 sum += _weight[(*_matching)[n]];
2987 /// \brief Returns true when the edge is in the matching.
2989 /// Returns true when the edge is in the matching.
2990 bool matching(const UEdge& edge) const {
2991 return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
2994 /// \brief Returns the incident matching edge.
2996 /// Returns the incident matching edge from given node.
2997 Edge matching(const Node& node) const {
2998 return (*_matching)[node];
3001 /// \brief Returns the mate of the node.
3003 /// Returns the adjancent node in a mathcing edge.
3004 Node mate(const Node& node) const {
3005 return _ugraph.target((*_matching)[node]);
3010 /// \name Dual solution
3011 /// Functions for get the dual solution.
3015 /// \brief Returns the value of the dual solution.
3017 /// Returns the value of the dual solution. It should be equal to
3018 /// the primal value scaled by \ref dualScale "dual scale".
3019 Value dualValue() const {
3021 for (NodeIt n(_ugraph); n != INVALID; ++n) {
3022 sum += nodeValue(n);
3024 for (int i = 0; i < blossomNum(); ++i) {
3025 sum += blossomValue(i) * (blossomSize(i) / 2);
3030 /// \brief Returns the value of the node.
3032 /// Returns the the value of the node.
3033 Value nodeValue(const Node& n) const {
3034 return (*_node_potential)[n];
3037 /// \brief Returns the number of the blossoms in the basis.
3039 /// Returns the number of the blossoms in the basis.
3041 int blossomNum() const {
3042 return _blossom_potential.size();
3046 /// \brief Returns the number of the nodes in the blossom.
3048 /// Returns the number of the nodes in the blossom.
3049 int blossomSize(int k) const {
3050 return _blossom_potential[k].end - _blossom_potential[k].begin;
3053 /// \brief Returns the value of the blossom.
3055 /// Returns the the value of the blossom.
3057 Value blossomValue(int k) const {
3058 return _blossom_potential[k].value;
3061 /// \brief Lemon iterator for get the items of the blossom.
3063 /// Lemon iterator for get the nodes of the blossom. This class
3064 /// provides a common style lemon iterator which gives back a
3065 /// subset of the nodes.
3069 /// \brief Constructor.
3071 /// Constructor for get the nodes of the variable.
3072 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3073 : _algorithm(&algorithm)
3075 _index = _algorithm->_blossom_potential[variable].begin;
3076 _last = _algorithm->_blossom_potential[variable].end;
3079 /// \brief Invalid constructor.
3081 /// Invalid constructor.
3082 BlossomIt(Invalid) : _index(-1) {}
3084 /// \brief Conversion to node.
3086 /// Conversion to node.
3087 operator Node() const {
3088 return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
3091 /// \brief Increment operator.
3093 /// Increment operator.
3094 BlossomIt& operator++() {
3096 if (_index == _last) {
3102 bool operator==(const BlossomIt& it) const {
3103 return _index == it._index;
3105 bool operator!=(const BlossomIt& it) const {
3106 return _index != it._index;
3110 const MaxWeightedPerfectMatching* _algorithm;
3120 } //END OF NAMESPACE LEMON
3122 #endif //LEMON_MAX_MATCHING_H