1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_MAX_MATCHING_H
20 #define LEMON_MAX_MATCHING_H
27 #include <lemon/core.h>
28 #include <lemon/unionfind.h>
29 #include <lemon/bin_heap.h>
30 #include <lemon/maps.h>
34 ///\brief Maximum matching algorithms in graph.
40 ///\brief Edmonds' alternating forest maximum matching algorithm.
42 ///This class provides Edmonds' alternating forest matching
43 ///algorithm. The starting matching (if any) can be passed to the
44 ///algorithm using some of init functions.
46 ///The dual side of a matching is a map of the nodes to
47 ///MaxMatching::DecompType, having values \c D, \c A and \c C
48 ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
49 ///in \c D induce a digraph with factor-critical components, the nodes
50 ///in \c A form the barrier, and the nodes in \c C induce a digraph
51 ///having a perfect matching. This decomposition can be attained by
52 ///calling \c decomposition() after running the algorithm.
54 ///\param Digraph The graph type the algorithm runs on.
55 template <typename Graph>
60 TEMPLATE_GRAPH_TYPEDEFS(Graph);
62 typedef typename Graph::template NodeMap<int> UFECrossRef;
63 typedef UnionFindEnum<UFECrossRef> UFE;
64 typedef std::vector<Node> NV;
66 typedef typename Graph::template NodeMap<int> EFECrossRef;
67 typedef ExtendFindEnum<EFECrossRef> EFE;
71 ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
73 ///Indicates the Gallai-Edmonds decomposition of the digraph, which
74 ///shows an upper bound on the size of a maximum matching. The
75 ///nodes with DecompType \c D induce a digraph with factor-critical
76 ///components, the nodes in \c A form the canonical barrier, and the
77 ///nodes in \c C induce a digraph having a perfect matching.
86 static const int HEUR_density=2;
88 typename Graph::template NodeMap<Node> _mate;
89 typename Graph::template NodeMap<DecompType> position;
93 MaxMatching(const Graph& _g)
94 : g(_g), _mate(_g), position(_g) {}
96 ///\brief Sets the actual matching to the empty matching.
98 ///Sets the actual matching to the empty matching.
101 for(NodeIt v(g); v!=INVALID; ++v) {
102 _mate.set(v,INVALID);
107 ///\brief Finds a greedy matching for initial matching.
109 ///For initial matchig it finds a maximal greedy matching.
111 for(NodeIt v(g); v!=INVALID; ++v) {
112 _mate.set(v,INVALID);
115 for(NodeIt v(g); v!=INVALID; ++v)
116 if ( _mate[v]==INVALID ) {
117 for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
118 Node y=g.runningNode(e);
119 if ( _mate[y]==INVALID && y!=v ) {
128 ///\brief Initialize the matching from each nodes' mate.
130 ///Initialize the matching from a \c Node valued \c Node map. This
131 ///map must be \e symmetric, i.e. if \c map[u]==v then \c
132 ///map[v]==u must hold, and \c uv will be an arc of the initial
134 template <typename MateMap>
135 void mateMapInit(MateMap& map) {
136 for(NodeIt v(g); v!=INVALID; ++v) {
142 ///\brief Initialize the matching from a node map with the
143 ///incident matching arcs.
145 ///Initialize the matching from an \c Edge valued \c Node map. \c
146 ///map[v] must be an \c Edge incident to \c v. This map must have
147 ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
148 ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
149 ///u to \c v will be an arc of the matching.
150 template<typename MatchingMap>
151 void matchingMapInit(MatchingMap& map) {
152 for(NodeIt v(g); v!=INVALID; ++v) {
156 _mate.set(v,g.oppositeNode(v,e));
158 _mate.set(v,INVALID);
162 ///\brief Initialize the matching from the map containing the
163 ///undirected matching arcs.
165 ///Initialize the matching from a \c bool valued \c Edge map. This
166 ///map must have the property that there are no two incident arcs
167 ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
168 ///map[e]==true form the matching.
169 template <typename MatchingMap>
170 void matchingInit(MatchingMap& map) {
171 for(NodeIt v(g); v!=INVALID; ++v) {
172 _mate.set(v,INVALID);
175 for(EdgeIt e(g); e!=INVALID; ++e) {
186 ///\brief Runs Edmonds' algorithm.
188 ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
189 ///2*number of nodes), and a heuristical Edmonds' algorithm with a
190 ///heuristic of postponing shrinks for dense digraphs.
192 if (countEdges(g) < HEUR_density * countNodes(g)) {
202 ///\brief Starts Edmonds' algorithm.
204 ///If runs the original Edmonds' algorithm.
207 typename Graph::template NodeMap<Node> ear(g,INVALID);
208 //undefined for the base nodes of the blossoms (i.e. for the
209 //representative elements of UFE blossom) and for the nodes in C
211 UFECrossRef blossom_base(g);
212 UFE blossom(blossom_base);
213 NV rep(countNodes(g));
215 EFECrossRef tree_base(g);
218 //If these UFE's would be members of the class then also
219 //blossom_base and tree_base should be a member.
221 //We build only one tree and the other vertices uncovered by the
222 //matching belong to C. (They can be considered as singleton
223 //trees.) If this tree can be augmented or no more
224 //grow/augmentation/shrink is possible then we return to this
226 for(NodeIt v(g); v!=INVALID; ++v) {
227 if (position[v]==C && _mate[v]==INVALID) {
228 rep[blossom.insert(v)] = v;
231 normShrink(v, ear, blossom, rep, tree);
236 ///\brief Starts Edmonds' algorithm.
238 ///It runs Edmonds' algorithm with a heuristic of postponing
239 ///shrinks, giving a faster algorithm for dense digraphs.
242 typename Graph::template NodeMap<Node> ear(g,INVALID);
243 //undefined for the base nodes of the blossoms (i.e. for the
244 //representative elements of UFE blossom) and for the nodes in C
246 UFECrossRef blossom_base(g);
247 UFE blossom(blossom_base);
248 NV rep(countNodes(g));
250 EFECrossRef tree_base(g);
253 //If these UFE's would be members of the class then also
254 //blossom_base and tree_base should be a member.
256 //We build only one tree and the other vertices uncovered by the
257 //matching belong to C. (They can be considered as singleton
258 //trees.) If this tree can be augmented or no more
259 //grow/augmentation/shrink is possible then we return to this
261 for(NodeIt v(g); v!=INVALID; ++v) {
262 if ( position[v]==C && _mate[v]==INVALID ) {
263 rep[blossom.insert(v)] = v;
266 lateShrink(v, ear, blossom, rep, tree);
273 ///\brief Returns the size of the actual matching stored.
275 ///Returns the size of the actual matching stored. After \ref
276 ///run() it returns the size of a maximum matching in the digraph.
279 for(NodeIt v(g); v!=INVALID; ++v) {
280 if ( _mate[v]!=INVALID ) {
288 ///\brief Returns the mate of a node in the actual matching.
290 ///Returns the mate of a \c node in the actual matching.
291 ///Returns INVALID if the \c node is not covered by the actual matching.
292 Node mate(const Node& node) const {
296 ///\brief Returns the matching arc incident to the given node.
298 ///Returns the matching arc of a \c node in the actual matching.
299 ///Returns INVALID if the \c node is not covered by the actual matching.
300 Edge matchingArc(const Node& node) const {
301 if (_mate[node] == INVALID) return INVALID;
302 Node n = node < _mate[node] ? node : _mate[node];
303 for (IncEdgeIt e(g, n); e != INVALID; ++e) {
304 if (g.oppositeNode(n, e) == _mate[n]) {
311 /// \brief Returns the class of the node in the Edmonds-Gallai
314 /// Returns the class of the node in the Edmonds-Gallai
316 DecompType decomposition(const Node& n) {
317 return position[n] == A;
320 /// \brief Returns true when the node is in the barrier.
322 /// Returns true when the node is in the barrier.
323 bool barrier(const Node& n) {
324 return position[n] == A;
327 ///\brief Gives back the matching in a \c Node of mates.
329 ///Writes the stored matching to a \c Node valued \c Node map. The
330 ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
331 ///map[v]==u will hold, and now \c uv is an arc of the matching.
332 template <typename MateMap>
333 void mateMap(MateMap& map) const {
334 for(NodeIt v(g); v!=INVALID; ++v) {
339 ///\brief Gives back the matching in an \c Edge valued \c Node
342 ///Writes the stored matching to an \c Edge valued \c Node
343 ///map. \c map[v] will be an \c Edge incident to \c v. This
344 ///map will have the property that if \c g.oppositeNode(u,map[u])
345 ///== v then \c map[u]==map[v] holds, and now this arc is an arc
347 template <typename MatchingMap>
348 void matchingMap(MatchingMap& map) const {
349 typename Graph::template NodeMap<bool> todo(g,true);
350 for(NodeIt v(g); v!=INVALID; ++v) {
351 if (_mate[v]!=INVALID && v < _mate[v]) {
353 for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
354 if ( g.runningNode(e) == u ) {
367 ///\brief Gives back the matching in a \c bool valued \c Edge
370 ///Writes the matching stored to a \c bool valued \c Arc
371 ///map. This map will have the property that there are no two
372 ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
373 ///arcs \c e with \c map[e]==true form the matching.
374 template<typename MatchingMap>
375 void matching(MatchingMap& map) const {
376 for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
378 typename Graph::template NodeMap<bool> todo(g,true);
379 for(NodeIt v(g); v!=INVALID; ++v) {
380 if ( todo[v] && _mate[v]!=INVALID ) {
382 for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
383 if ( g.runningNode(e) == u ) {
395 ///\brief Returns the canonical decomposition of the digraph after running
398 ///After calling any run methods of the class, it writes the
399 ///Gallai-Edmonds canonical decomposition of the digraph. \c map
400 ///must be a node map of \ref DecompType 's.
401 template <typename DecompositionMap>
402 void decomposition(DecompositionMap& map) const {
403 for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
406 ///\brief Returns a barrier on the nodes.
408 ///After calling any run methods of the class, it writes a
409 ///canonical barrier on the nodes. The odd component number of the
410 ///remaining digraph minus the barrier size is a lower bound for the
411 ///uncovered nodes in the digraph. The \c map must be a node map of
413 template <typename BarrierMap>
414 void barrier(BarrierMap& barrier) {
415 for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
421 void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
422 UFE& blossom, NV& rep, EFE& tree) {
423 //We have one tree which we grow, and also shrink but only if it
424 //cannot be postponed. If we augment then we return to the "for"
425 //cycle of runEdmonds().
427 std::queue<Node> Q; //queue of the totally unscanned nodes
430 //queue of the nodes which must be scanned for a possible shrink
432 while ( !Q.empty() ) {
435 for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
436 Node y=g.runningNode(e);
437 //growOrAugment grows if y is covered by the matching and
438 //augments if not. In this latter case it returns 1.
439 if (position[y]==C &&
440 growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
445 while ( !R.empty() ) {
449 for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
450 Node y=g.runningNode(e);
452 if ( position[y] == D && blossom.find(x) != blossom.find(y) )
453 //Recall that we have only one tree.
454 shrink( x, y, ear, blossom, rep, tree, Q);
456 while ( !Q.empty() ) {
459 for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
460 Node w=g.runningNode(f);
461 //growOrAugment grows if y is covered by the matching and
462 //augments if not. In this latter case it returns 1.
463 if (position[w]==C &&
464 growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
469 } // while ( !R.empty() )
472 void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
473 UFE& blossom, NV& rep, EFE& tree) {
474 //We have one tree, which we grow and shrink. If we augment then we
475 //return to the "for" cycle of runEdmonds().
477 std::queue<Node> Q; //queue of the unscanned nodes
479 while ( !Q.empty() ) {
484 for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
485 Node y=g.runningNode(e);
487 switch ( position[y] ) {
488 case D: //x and y must be in the same tree
489 if ( blossom.find(x) != blossom.find(y))
490 //x and y are in the same tree
491 shrink(x, y, ear, blossom, rep, tree, Q);
494 //growOrAugment grows if y is covered by the matching and
495 //augments if not. In this latter case it returns 1.
496 if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
504 void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
505 UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
506 //x and y are the two adjacent vertices in two blossoms.
508 typename Graph::template NodeMap<bool> path(g,false);
510 Node b=rep[blossom.find(x)];
513 while ( b!=INVALID ) {
514 b=rep[blossom.find(ear[b])];
517 } //we go until the root through bases of blossoms and odd vertices
520 Node middle=rep[blossom.find(top)];
522 while ( !path[middle] )
523 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
524 //Until we arrive to a node on the path, we update blossom, tree
525 //and the positions of the odd nodes.
529 middle=rep[blossom.find(top)];
531 Node blossom_base=rep[blossom.find(base)];
532 while ( middle!=blossom_base )
533 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
534 //Until we arrive to a node on the path, we update blossom, tree
535 //and the positions of the odd nodes.
537 rep[blossom.find(base)] = base;
540 void shrinkStep(Node& top, Node& middle, Node& bottom,
541 typename Graph::template NodeMap<Node>& ear,
542 UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
543 //We traverse a blossom and update everything.
547 while ( t!=middle ) {
552 bottom=_mate[middle];
553 position.set(bottom,D);
556 Node oldmiddle=middle;
557 middle=rep[blossom.find(top)];
559 tree.erase(oldmiddle);
560 blossom.insert(bottom);
561 blossom.join(bottom, oldmiddle);
562 blossom.join(top, oldmiddle);
567 bool growOrAugment(Node& y, Node& x, typename Graph::template
568 NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
569 std::queue<Node>& Q) {
570 //x is in a blossom in the tree, y is outside. If y is covered by
571 //the matching we grow, otherwise we augment. In this case we
574 if ( _mate[y]!=INVALID ) { //grow
577 rep[blossom.insert(w)] = w;
580 int t = tree.find(rep[blossom.find(x)]);
585 augment(x, ear, blossom, rep, tree);
593 void augment(Node x, typename Graph::template NodeMap<Node>& ear,
594 UFE& blossom, NV& rep, EFE& tree) {
596 while ( v!=INVALID ) {
604 int y = tree.find(rep[blossom.find(x)]);
605 for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
606 if ( position[tit] == D ) {
607 int b = blossom.find(tit);
608 for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
609 position.set(bit, C);
611 blossom.eraseClass(b);
612 } else position.set(tit, C);
620 /// \ingroup matching
622 /// \brief Weighted matching in general graphs
624 /// This class provides an efficient implementation of Edmond's
625 /// maximum weighted matching algorithm. The implementation is based
626 /// on extensive use of priority queues and provides
627 /// \f$O(nm\log(n))\f$ time complexity.
629 /// The maximum weighted matching problem is to find undirected
630 /// arcs in the digraph with maximum overall weight and no two of
631 /// them shares their endpoints. The problem can be formulated with
632 /// the next linear program:
633 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
634 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
635 /// \f[x_e \ge 0\quad \forall e\in E\f]
636 /// \f[\max \sum_{e\in E}x_ew_e\f]
637 /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
638 /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
639 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
642 /// The algorithm calculates an optimal matching and a proof of the
643 /// optimality. The solution of the dual problem can be used to check
644 /// the result of the algorithm. The dual linear problem is the next:
645 /// \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]
646 /// \f[y_u \ge 0 \quad \forall u \in V\f]
647 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
648 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
650 /// The algorithm can be executed with \c run() or the \c init() and
651 /// then the \c start() member functions. After it the matching can
652 /// be asked with \c matching() or mate() functions. The dual
653 /// solution can be get with \c nodeValue(), \c blossomNum() and \c
654 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
655 /// "BlossomIt" nested class which is able to iterate on the nodes
656 /// of a blossom. If the value type is integral then the dual
657 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
658 template <typename _Graph,
659 typename _WeightMap = typename _Graph::template EdgeMap<int> >
660 class MaxWeightedMatching {
663 typedef _Graph Graph;
664 typedef _WeightMap WeightMap;
665 typedef typename WeightMap::Value Value;
667 /// \brief Scaling factor for dual solution
669 /// Scaling factor for dual solution, it is equal to 4 or 1
670 /// according to the value type.
671 static const int dualScale =
672 std::numeric_limits<Value>::is_integer ? 4 : 1;
674 typedef typename Graph::template NodeMap<typename Graph::Arc>
679 TEMPLATE_GRAPH_TYPEDEFS(Graph);
681 typedef typename Graph::template NodeMap<Value> NodePotential;
682 typedef std::vector<Node> BlossomNodeList;
684 struct BlossomVariable {
688 BlossomVariable(int _begin, int _end, Value _value)
689 : begin(_begin), end(_end), value(_value) {}
693 typedef std::vector<BlossomVariable> BlossomPotential;
696 const WeightMap& _weight;
698 MatchingMap* _matching;
700 NodePotential* _node_potential;
702 BlossomPotential _blossom_potential;
703 BlossomNodeList _blossom_node_list;
708 typedef typename Graph::template NodeMap<int> NodeIntMap;
709 typedef typename Graph::template ArcMap<int> ArcIntMap;
710 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
711 typedef RangeMap<int> IntIntMap;
714 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
717 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
726 NodeIntMap *_blossom_index;
727 BlossomSet *_blossom_set;
728 RangeMap<BlossomData>* _blossom_data;
730 NodeIntMap *_node_index;
731 ArcIntMap *_node_heap_index;
735 NodeData(ArcIntMap& node_heap_index)
736 : heap(node_heap_index) {}
740 BinHeap<Value, ArcIntMap> heap;
741 std::map<int, Arc> heap_index;
746 RangeMap<NodeData>* _node_data;
748 typedef ExtendFindEnum<IntIntMap> TreeSet;
750 IntIntMap *_tree_set_index;
753 NodeIntMap *_delta1_index;
754 BinHeap<Value, NodeIntMap> *_delta1;
756 IntIntMap *_delta2_index;
757 BinHeap<Value, IntIntMap> *_delta2;
759 EdgeIntMap *_delta3_index;
760 BinHeap<Value, EdgeIntMap> *_delta3;
762 IntIntMap *_delta4_index;
763 BinHeap<Value, IntIntMap> *_delta4;
767 void createStructures() {
768 _node_num = countNodes(_graph);
769 _blossom_num = _node_num * 3 / 2;
772 _matching = new MatchingMap(_graph);
774 if (!_node_potential) {
775 _node_potential = new NodePotential(_graph);
778 _blossom_index = new NodeIntMap(_graph);
779 _blossom_set = new BlossomSet(*_blossom_index);
780 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
784 _node_index = new NodeIntMap(_graph);
785 _node_heap_index = new ArcIntMap(_graph);
786 _node_data = new RangeMap<NodeData>(_node_num,
787 NodeData(*_node_heap_index));
791 _tree_set_index = new IntIntMap(_blossom_num);
792 _tree_set = new TreeSet(*_tree_set_index);
795 _delta1_index = new NodeIntMap(_graph);
796 _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
799 _delta2_index = new IntIntMap(_blossom_num);
800 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
803 _delta3_index = new EdgeIntMap(_graph);
804 _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
807 _delta4_index = new IntIntMap(_blossom_num);
808 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
812 void destroyStructures() {
813 _node_num = countNodes(_graph);
814 _blossom_num = _node_num * 3 / 2;
819 if (_node_potential) {
820 delete _node_potential;
823 delete _blossom_index;
825 delete _blossom_data;
830 delete _node_heap_index;
835 delete _tree_set_index;
839 delete _delta1_index;
843 delete _delta2_index;
847 delete _delta3_index;
851 delete _delta4_index;
856 void matchedToEven(int blossom, int tree) {
857 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
858 _delta2->erase(blossom);
861 if (!_blossom_set->trivial(blossom)) {
862 (*_blossom_data)[blossom].pot -=
863 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
866 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
869 _blossom_set->increase(n, std::numeric_limits<Value>::max());
870 int ni = (*_node_index)[n];
872 (*_node_data)[ni].heap.clear();
873 (*_node_data)[ni].heap_index.clear();
875 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
877 _delta1->push(n, (*_node_data)[ni].pot);
879 for (InArcIt e(_graph, n); e != INVALID; ++e) {
880 Node v = _graph.source(e);
881 int vb = _blossom_set->find(v);
882 int vi = (*_node_index)[v];
884 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
885 dualScale * _weight[e];
887 if ((*_blossom_data)[vb].status == EVEN) {
888 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
889 _delta3->push(e, rw / 2);
891 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
892 if (_delta3->state(e) != _delta3->IN_HEAP) {
893 _delta3->push(e, rw);
896 typename std::map<int, Arc>::iterator it =
897 (*_node_data)[vi].heap_index.find(tree);
899 if (it != (*_node_data)[vi].heap_index.end()) {
900 if ((*_node_data)[vi].heap[it->second] > rw) {
901 (*_node_data)[vi].heap.replace(it->second, e);
902 (*_node_data)[vi].heap.decrease(e, rw);
906 (*_node_data)[vi].heap.push(e, rw);
907 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
910 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
911 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
913 if ((*_blossom_data)[vb].status == MATCHED) {
914 if (_delta2->state(vb) != _delta2->IN_HEAP) {
915 _delta2->push(vb, _blossom_set->classPrio(vb) -
916 (*_blossom_data)[vb].offset);
917 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
918 (*_blossom_data)[vb].offset){
919 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
920 (*_blossom_data)[vb].offset);
927 (*_blossom_data)[blossom].offset = 0;
930 void matchedToOdd(int blossom) {
931 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
932 _delta2->erase(blossom);
934 (*_blossom_data)[blossom].offset += _delta_sum;
935 if (!_blossom_set->trivial(blossom)) {
936 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
937 (*_blossom_data)[blossom].offset);
941 void evenToMatched(int blossom, int tree) {
942 if (!_blossom_set->trivial(blossom)) {
943 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
946 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
948 int ni = (*_node_index)[n];
949 (*_node_data)[ni].pot -= _delta_sum;
953 for (InArcIt e(_graph, n); e != INVALID; ++e) {
954 Node v = _graph.source(e);
955 int vb = _blossom_set->find(v);
956 int vi = (*_node_index)[v];
958 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
959 dualScale * _weight[e];
962 if (_delta3->state(e) == _delta3->IN_HEAP) {
965 } else if ((*_blossom_data)[vb].status == EVEN) {
967 if (_delta3->state(e) == _delta3->IN_HEAP) {
971 int vt = _tree_set->find(vb);
975 Arc r = _graph.oppositeArc(e);
977 typename std::map<int, Arc>::iterator it =
978 (*_node_data)[ni].heap_index.find(vt);
980 if (it != (*_node_data)[ni].heap_index.end()) {
981 if ((*_node_data)[ni].heap[it->second] > rw) {
982 (*_node_data)[ni].heap.replace(it->second, r);
983 (*_node_data)[ni].heap.decrease(r, rw);
987 (*_node_data)[ni].heap.push(r, rw);
988 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
991 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
992 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
994 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
995 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
996 (*_blossom_data)[blossom].offset);
997 } else if ((*_delta2)[blossom] >
998 _blossom_set->classPrio(blossom) -
999 (*_blossom_data)[blossom].offset){
1000 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1001 (*_blossom_data)[blossom].offset);
1006 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1007 if (_delta3->state(e) == _delta3->IN_HEAP) {
1012 typename std::map<int, Arc>::iterator it =
1013 (*_node_data)[vi].heap_index.find(tree);
1015 if (it != (*_node_data)[vi].heap_index.end()) {
1016 (*_node_data)[vi].heap.erase(it->second);
1017 (*_node_data)[vi].heap_index.erase(it);
1018 if ((*_node_data)[vi].heap.empty()) {
1019 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1020 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1021 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1024 if ((*_blossom_data)[vb].status == MATCHED) {
1025 if (_blossom_set->classPrio(vb) ==
1026 std::numeric_limits<Value>::max()) {
1028 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1029 (*_blossom_data)[vb].offset) {
1030 _delta2->increase(vb, _blossom_set->classPrio(vb) -
1031 (*_blossom_data)[vb].offset);
1040 void oddToMatched(int blossom) {
1041 (*_blossom_data)[blossom].offset -= _delta_sum;
1043 if (_blossom_set->classPrio(blossom) !=
1044 std::numeric_limits<Value>::max()) {
1045 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1046 (*_blossom_data)[blossom].offset);
1049 if (!_blossom_set->trivial(blossom)) {
1050 _delta4->erase(blossom);
1054 void oddToEven(int blossom, int tree) {
1055 if (!_blossom_set->trivial(blossom)) {
1056 _delta4->erase(blossom);
1057 (*_blossom_data)[blossom].pot -=
1058 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1061 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1062 n != INVALID; ++n) {
1063 int ni = (*_node_index)[n];
1065 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1067 (*_node_data)[ni].heap.clear();
1068 (*_node_data)[ni].heap_index.clear();
1069 (*_node_data)[ni].pot +=
1070 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1072 _delta1->push(n, (*_node_data)[ni].pot);
1074 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1075 Node v = _graph.source(e);
1076 int vb = _blossom_set->find(v);
1077 int vi = (*_node_index)[v];
1079 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1080 dualScale * _weight[e];
1082 if ((*_blossom_data)[vb].status == EVEN) {
1083 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1084 _delta3->push(e, rw / 2);
1086 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1087 if (_delta3->state(e) != _delta3->IN_HEAP) {
1088 _delta3->push(e, rw);
1092 typename std::map<int, Arc>::iterator it =
1093 (*_node_data)[vi].heap_index.find(tree);
1095 if (it != (*_node_data)[vi].heap_index.end()) {
1096 if ((*_node_data)[vi].heap[it->second] > rw) {
1097 (*_node_data)[vi].heap.replace(it->second, e);
1098 (*_node_data)[vi].heap.decrease(e, rw);
1102 (*_node_data)[vi].heap.push(e, rw);
1103 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1106 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1107 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1109 if ((*_blossom_data)[vb].status == MATCHED) {
1110 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1111 _delta2->push(vb, _blossom_set->classPrio(vb) -
1112 (*_blossom_data)[vb].offset);
1113 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1114 (*_blossom_data)[vb].offset) {
1115 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1116 (*_blossom_data)[vb].offset);
1123 (*_blossom_data)[blossom].offset = 0;
1127 void matchedToUnmatched(int blossom) {
1128 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1129 _delta2->erase(blossom);
1132 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1133 n != INVALID; ++n) {
1134 int ni = (*_node_index)[n];
1136 _blossom_set->increase(n, std::numeric_limits<Value>::max());
1138 (*_node_data)[ni].heap.clear();
1139 (*_node_data)[ni].heap_index.clear();
1141 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1142 Node v = _graph.target(e);
1143 int vb = _blossom_set->find(v);
1144 int vi = (*_node_index)[v];
1146 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1147 dualScale * _weight[e];
1149 if ((*_blossom_data)[vb].status == EVEN) {
1150 if (_delta3->state(e) != _delta3->IN_HEAP) {
1151 _delta3->push(e, rw);
1158 void unmatchedToMatched(int blossom) {
1159 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1160 n != INVALID; ++n) {
1161 int ni = (*_node_index)[n];
1163 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1164 Node v = _graph.source(e);
1165 int vb = _blossom_set->find(v);
1166 int vi = (*_node_index)[v];
1168 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1169 dualScale * _weight[e];
1171 if (vb == blossom) {
1172 if (_delta3->state(e) == _delta3->IN_HEAP) {
1175 } else if ((*_blossom_data)[vb].status == EVEN) {
1177 if (_delta3->state(e) == _delta3->IN_HEAP) {
1181 int vt = _tree_set->find(vb);
1183 Arc r = _graph.oppositeArc(e);
1185 typename std::map<int, Arc>::iterator it =
1186 (*_node_data)[ni].heap_index.find(vt);
1188 if (it != (*_node_data)[ni].heap_index.end()) {
1189 if ((*_node_data)[ni].heap[it->second] > rw) {
1190 (*_node_data)[ni].heap.replace(it->second, r);
1191 (*_node_data)[ni].heap.decrease(r, rw);
1195 (*_node_data)[ni].heap.push(r, rw);
1196 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1199 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1200 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1202 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1203 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1204 (*_blossom_data)[blossom].offset);
1205 } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1206 (*_blossom_data)[blossom].offset){
1207 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1208 (*_blossom_data)[blossom].offset);
1212 } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1213 if (_delta3->state(e) == _delta3->IN_HEAP) {
1221 void alternatePath(int even, int tree) {
1224 evenToMatched(even, tree);
1225 (*_blossom_data)[even].status = MATCHED;
1227 while ((*_blossom_data)[even].pred != INVALID) {
1228 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1229 (*_blossom_data)[odd].status = MATCHED;
1231 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1233 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1234 (*_blossom_data)[even].status = MATCHED;
1235 evenToMatched(even, tree);
1236 (*_blossom_data)[even].next =
1237 _graph.oppositeArc((*_blossom_data)[odd].pred);
1242 void destroyTree(int tree) {
1243 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1244 if ((*_blossom_data)[b].status == EVEN) {
1245 (*_blossom_data)[b].status = MATCHED;
1246 evenToMatched(b, tree);
1247 } else if ((*_blossom_data)[b].status == ODD) {
1248 (*_blossom_data)[b].status = MATCHED;
1252 _tree_set->eraseClass(tree);
1256 void unmatchNode(const Node& node) {
1257 int blossom = _blossom_set->find(node);
1258 int tree = _tree_set->find(blossom);
1260 alternatePath(blossom, tree);
1263 (*_blossom_data)[blossom].status = UNMATCHED;
1264 (*_blossom_data)[blossom].base = node;
1265 matchedToUnmatched(blossom);
1269 void augmentOnArc(const Edge& arc) {
1271 int left = _blossom_set->find(_graph.u(arc));
1272 int right = _blossom_set->find(_graph.v(arc));
1274 if ((*_blossom_data)[left].status == EVEN) {
1275 int left_tree = _tree_set->find(left);
1276 alternatePath(left, left_tree);
1277 destroyTree(left_tree);
1279 (*_blossom_data)[left].status = MATCHED;
1280 unmatchedToMatched(left);
1283 if ((*_blossom_data)[right].status == EVEN) {
1284 int right_tree = _tree_set->find(right);
1285 alternatePath(right, right_tree);
1286 destroyTree(right_tree);
1288 (*_blossom_data)[right].status = MATCHED;
1289 unmatchedToMatched(right);
1292 (*_blossom_data)[left].next = _graph.direct(arc, true);
1293 (*_blossom_data)[right].next = _graph.direct(arc, false);
1296 void extendOnArc(const Arc& arc) {
1297 int base = _blossom_set->find(_graph.target(arc));
1298 int tree = _tree_set->find(base);
1300 int odd = _blossom_set->find(_graph.source(arc));
1301 _tree_set->insert(odd, tree);
1302 (*_blossom_data)[odd].status = ODD;
1304 (*_blossom_data)[odd].pred = arc;
1306 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1307 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1308 _tree_set->insert(even, tree);
1309 (*_blossom_data)[even].status = EVEN;
1310 matchedToEven(even, tree);
1313 void shrinkOnArc(const Edge& edge, int tree) {
1315 std::vector<int> left_path, right_path;
1318 std::set<int> left_set, right_set;
1319 int left = _blossom_set->find(_graph.u(edge));
1320 left_path.push_back(left);
1321 left_set.insert(left);
1323 int right = _blossom_set->find(_graph.v(edge));
1324 right_path.push_back(right);
1325 right_set.insert(right);
1329 if ((*_blossom_data)[left].pred == INVALID) break;
1332 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1333 left_path.push_back(left);
1335 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1336 left_path.push_back(left);
1338 left_set.insert(left);
1340 if (right_set.find(left) != right_set.end()) {
1345 if ((*_blossom_data)[right].pred == INVALID) break;
1348 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1349 right_path.push_back(right);
1351 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1352 right_path.push_back(right);
1354 right_set.insert(right);
1356 if (left_set.find(right) != left_set.end()) {
1364 if ((*_blossom_data)[left].pred == INVALID) {
1366 while (left_set.find(nca) == left_set.end()) {
1368 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1369 right_path.push_back(nca);
1371 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1372 right_path.push_back(nca);
1376 while (right_set.find(nca) == right_set.end()) {
1378 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1379 left_path.push_back(nca);
1381 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1382 left_path.push_back(nca);
1388 std::vector<int> subblossoms;
1391 prev = _graph.direct(edge, true);
1392 for (int i = 0; left_path[i] != nca; i += 2) {
1393 subblossoms.push_back(left_path[i]);
1394 (*_blossom_data)[left_path[i]].next = prev;
1395 _tree_set->erase(left_path[i]);
1397 subblossoms.push_back(left_path[i + 1]);
1398 (*_blossom_data)[left_path[i + 1]].status = EVEN;
1399 oddToEven(left_path[i + 1], tree);
1400 _tree_set->erase(left_path[i + 1]);
1401 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1405 while (right_path[k] != nca) ++k;
1407 subblossoms.push_back(nca);
1408 (*_blossom_data)[nca].next = prev;
1410 for (int i = k - 2; i >= 0; i -= 2) {
1411 subblossoms.push_back(right_path[i + 1]);
1412 (*_blossom_data)[right_path[i + 1]].status = EVEN;
1413 oddToEven(right_path[i + 1], tree);
1414 _tree_set->erase(right_path[i + 1]);
1416 (*_blossom_data)[right_path[i + 1]].next =
1417 (*_blossom_data)[right_path[i + 1]].pred;
1419 subblossoms.push_back(right_path[i]);
1420 _tree_set->erase(right_path[i]);
1424 _blossom_set->join(subblossoms.begin(), subblossoms.end());
1426 for (int i = 0; i < int(subblossoms.size()); ++i) {
1427 if (!_blossom_set->trivial(subblossoms[i])) {
1428 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1430 (*_blossom_data)[subblossoms[i]].status = MATCHED;
1433 (*_blossom_data)[surface].pot = -2 * _delta_sum;
1434 (*_blossom_data)[surface].offset = 0;
1435 (*_blossom_data)[surface].status = EVEN;
1436 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1437 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1439 _tree_set->insert(surface, tree);
1440 _tree_set->erase(nca);
1443 void splitBlossom(int blossom) {
1444 Arc next = (*_blossom_data)[blossom].next;
1445 Arc pred = (*_blossom_data)[blossom].pred;
1447 int tree = _tree_set->find(blossom);
1449 (*_blossom_data)[blossom].status = MATCHED;
1450 oddToMatched(blossom);
1451 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1452 _delta2->erase(blossom);
1455 std::vector<int> subblossoms;
1456 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1458 Value offset = (*_blossom_data)[blossom].offset;
1459 int b = _blossom_set->find(_graph.source(pred));
1460 int d = _blossom_set->find(_graph.source(next));
1462 int ib = -1, id = -1;
1463 for (int i = 0; i < int(subblossoms.size()); ++i) {
1464 if (subblossoms[i] == b) ib = i;
1465 if (subblossoms[i] == d) id = i;
1467 (*_blossom_data)[subblossoms[i]].offset = offset;
1468 if (!_blossom_set->trivial(subblossoms[i])) {
1469 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1471 if (_blossom_set->classPrio(subblossoms[i]) !=
1472 std::numeric_limits<Value>::max()) {
1473 _delta2->push(subblossoms[i],
1474 _blossom_set->classPrio(subblossoms[i]) -
1475 (*_blossom_data)[subblossoms[i]].offset);
1479 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1480 for (int i = (id + 1) % subblossoms.size();
1481 i != ib; i = (i + 2) % subblossoms.size()) {
1482 int sb = subblossoms[i];
1483 int tb = subblossoms[(i + 1) % subblossoms.size()];
1484 (*_blossom_data)[sb].next =
1485 _graph.oppositeArc((*_blossom_data)[tb].next);
1488 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1489 int sb = subblossoms[i];
1490 int tb = subblossoms[(i + 1) % subblossoms.size()];
1491 int ub = subblossoms[(i + 2) % subblossoms.size()];
1493 (*_blossom_data)[sb].status = ODD;
1495 _tree_set->insert(sb, tree);
1496 (*_blossom_data)[sb].pred = pred;
1497 (*_blossom_data)[sb].next =
1498 _graph.oppositeArc((*_blossom_data)[tb].next);
1500 pred = (*_blossom_data)[ub].next;
1502 (*_blossom_data)[tb].status = EVEN;
1503 matchedToEven(tb, tree);
1504 _tree_set->insert(tb, tree);
1505 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1508 (*_blossom_data)[subblossoms[id]].status = ODD;
1509 matchedToOdd(subblossoms[id]);
1510 _tree_set->insert(subblossoms[id], tree);
1511 (*_blossom_data)[subblossoms[id]].next = next;
1512 (*_blossom_data)[subblossoms[id]].pred = pred;
1516 for (int i = (ib + 1) % subblossoms.size();
1517 i != id; i = (i + 2) % subblossoms.size()) {
1518 int sb = subblossoms[i];
1519 int tb = subblossoms[(i + 1) % subblossoms.size()];
1520 (*_blossom_data)[sb].next =
1521 _graph.oppositeArc((*_blossom_data)[tb].next);
1524 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1525 int sb = subblossoms[i];
1526 int tb = subblossoms[(i + 1) % subblossoms.size()];
1527 int ub = subblossoms[(i + 2) % subblossoms.size()];
1529 (*_blossom_data)[sb].status = ODD;
1531 _tree_set->insert(sb, tree);
1532 (*_blossom_data)[sb].next = next;
1533 (*_blossom_data)[sb].pred =
1534 _graph.oppositeArc((*_blossom_data)[tb].next);
1536 (*_blossom_data)[tb].status = EVEN;
1537 matchedToEven(tb, tree);
1538 _tree_set->insert(tb, tree);
1539 (*_blossom_data)[tb].pred =
1540 (*_blossom_data)[tb].next =
1541 _graph.oppositeArc((*_blossom_data)[ub].next);
1542 next = (*_blossom_data)[ub].next;
1545 (*_blossom_data)[subblossoms[ib]].status = ODD;
1546 matchedToOdd(subblossoms[ib]);
1547 _tree_set->insert(subblossoms[ib], tree);
1548 (*_blossom_data)[subblossoms[ib]].next = next;
1549 (*_blossom_data)[subblossoms[ib]].pred = pred;
1551 _tree_set->erase(blossom);
1554 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1555 if (_blossom_set->trivial(blossom)) {
1556 int bi = (*_node_index)[base];
1557 Value pot = (*_node_data)[bi].pot;
1559 _matching->set(base, matching);
1560 _blossom_node_list.push_back(base);
1561 _node_potential->set(base, pot);
1564 Value pot = (*_blossom_data)[blossom].pot;
1565 int bn = _blossom_node_list.size();
1567 std::vector<int> subblossoms;
1568 _blossom_set->split(blossom, std::back_inserter(subblossoms));
1569 int b = _blossom_set->find(base);
1571 for (int i = 0; i < int(subblossoms.size()); ++i) {
1572 if (subblossoms[i] == b) { ib = i; break; }
1575 for (int i = 1; i < int(subblossoms.size()); i += 2) {
1576 int sb = subblossoms[(ib + i) % subblossoms.size()];
1577 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1579 Arc m = (*_blossom_data)[tb].next;
1580 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1581 extractBlossom(tb, _graph.source(m), m);
1583 extractBlossom(subblossoms[ib], base, matching);
1585 int en = _blossom_node_list.size();
1587 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1591 void extractMatching() {
1592 std::vector<int> blossoms;
1593 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1594 blossoms.push_back(c);
1597 for (int i = 0; i < int(blossoms.size()); ++i) {
1598 if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1600 Value offset = (*_blossom_data)[blossoms[i]].offset;
1601 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1602 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1603 n != INVALID; ++n) {
1604 (*_node_data)[(*_node_index)[n]].pot -= offset;
1607 Arc matching = (*_blossom_data)[blossoms[i]].next;
1608 Node base = _graph.source(matching);
1609 extractBlossom(blossoms[i], base, matching);
1611 Node base = (*_blossom_data)[blossoms[i]].base;
1612 extractBlossom(blossoms[i], base, INVALID);
1619 /// \brief Constructor
1622 MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1623 : _graph(graph), _weight(weight), _matching(0),
1624 _node_potential(0), _blossom_potential(), _blossom_node_list(),
1625 _node_num(0), _blossom_num(0),
1627 _blossom_index(0), _blossom_set(0), _blossom_data(0),
1628 _node_index(0), _node_heap_index(0), _node_data(0),
1629 _tree_set_index(0), _tree_set(0),
1631 _delta1_index(0), _delta1(0),
1632 _delta2_index(0), _delta2(0),
1633 _delta3_index(0), _delta3(0),
1634 _delta4_index(0), _delta4(0),
1638 ~MaxWeightedMatching() {
1639 destroyStructures();
1642 /// \name Execution control
1643 /// The simplest way to execute the algorithm is to use the member
1644 /// \c run() member function.
1648 /// \brief Initialize the algorithm
1650 /// Initialize the algorithm
1654 for (ArcIt e(_graph); e != INVALID; ++e) {
1655 _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
1657 for (NodeIt n(_graph); n != INVALID; ++n) {
1658 _delta1_index->set(n, _delta1->PRE_HEAP);
1660 for (EdgeIt e(_graph); e != INVALID; ++e) {
1661 _delta3_index->set(e, _delta3->PRE_HEAP);
1663 for (int i = 0; i < _blossom_num; ++i) {
1664 _delta2_index->set(i, _delta2->PRE_HEAP);
1665 _delta4_index->set(i, _delta4->PRE_HEAP);
1669 for (NodeIt n(_graph); n != INVALID; ++n) {
1671 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1672 if (_graph.target(e) == n) continue;
1673 if ((dualScale * _weight[e]) / 2 > max) {
1674 max = (dualScale * _weight[e]) / 2;
1677 _node_index->set(n, index);
1678 (*_node_data)[index].pot = max;
1679 _delta1->push(n, max);
1681 _blossom_set->insert(n, std::numeric_limits<Value>::max());
1683 _tree_set->insert(blossom);
1685 (*_blossom_data)[blossom].status = EVEN;
1686 (*_blossom_data)[blossom].pred = INVALID;
1687 (*_blossom_data)[blossom].next = INVALID;
1688 (*_blossom_data)[blossom].pot = 0;
1689 (*_blossom_data)[blossom].offset = 0;
1692 for (EdgeIt e(_graph); e != INVALID; ++e) {
1693 int si = (*_node_index)[_graph.u(e)];
1694 int ti = (*_node_index)[_graph.v(e)];
1695 if (_graph.u(e) != _graph.v(e)) {
1696 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1697 dualScale * _weight[e]) / 2);
1702 /// \brief Starts the algorithm
1704 /// Starts the algorithm
1710 int unmatched = _node_num;
1711 while (unmatched > 0) {
1712 Value d1 = !_delta1->empty() ?
1713 _delta1->prio() : std::numeric_limits<Value>::max();
1715 Value d2 = !_delta2->empty() ?
1716 _delta2->prio() : std::numeric_limits<Value>::max();
1718 Value d3 = !_delta3->empty() ?
1719 _delta3->prio() : std::numeric_limits<Value>::max();
1721 Value d4 = !_delta4->empty() ?
1722 _delta4->prio() : std::numeric_limits<Value>::max();
1724 _delta_sum = d1; OpType ot = D1;
1725 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1726 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1727 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1733 Node n = _delta1->top();
1740 int blossom = _delta2->top();
1741 Node n = _blossom_set->classTop(blossom);
1742 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1748 Edge e = _delta3->top();
1750 int left_blossom = _blossom_set->find(_graph.u(e));
1751 int right_blossom = _blossom_set->find(_graph.v(e));
1753 if (left_blossom == right_blossom) {
1757 if ((*_blossom_data)[left_blossom].status == EVEN) {
1758 left_tree = _tree_set->find(left_blossom);
1764 if ((*_blossom_data)[right_blossom].status == EVEN) {
1765 right_tree = _tree_set->find(right_blossom);
1771 if (left_tree == right_tree) {
1772 shrinkOnArc(e, left_tree);
1780 splitBlossom(_delta4->top());
1787 /// \brief Runs %MaxWeightedMatching algorithm.
1789 /// This method runs the %MaxWeightedMatching algorithm.
1791 /// \note mwm.run() is just a shortcut of the following code.
1803 /// \name Primal solution
1804 /// Functions for get the primal solution, ie. the matching.
1808 /// \brief Returns the matching value.
1810 /// Returns the matching value.
1811 Value matchingValue() const {
1813 for (NodeIt n(_graph); n != INVALID; ++n) {
1814 if ((*_matching)[n] != INVALID) {
1815 sum += _weight[(*_matching)[n]];
1821 /// \brief Returns true when the arc is in the matching.
1823 /// Returns true when the arc is in the matching.
1824 bool matching(const Edge& arc) const {
1825 return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
1828 /// \brief Returns the incident matching arc.
1830 /// Returns the incident matching arc from given node. If the
1831 /// node is not matched then it gives back \c INVALID.
1832 Arc matching(const Node& node) const {
1833 return (*_matching)[node];
1836 /// \brief Returns the mate of the node.
1838 /// Returns the adjancent node in a mathcing arc. If the node is
1839 /// not matched then it gives back \c INVALID.
1840 Node mate(const Node& node) const {
1841 return (*_matching)[node] != INVALID ?
1842 _graph.target((*_matching)[node]) : INVALID;
1847 /// \name Dual solution
1848 /// Functions for get the dual solution.
1852 /// \brief Returns the value of the dual solution.
1854 /// Returns the value of the dual solution. It should be equal to
1855 /// the primal value scaled by \ref dualScale "dual scale".
1856 Value dualValue() const {
1858 for (NodeIt n(_graph); n != INVALID; ++n) {
1859 sum += nodeValue(n);
1861 for (int i = 0; i < blossomNum(); ++i) {
1862 sum += blossomValue(i) * (blossomSize(i) / 2);
1867 /// \brief Returns the value of the node.
1869 /// Returns the the value of the node.
1870 Value nodeValue(const Node& n) const {
1871 return (*_node_potential)[n];
1874 /// \brief Returns the number of the blossoms in the basis.
1876 /// Returns the number of the blossoms in the basis.
1878 int blossomNum() const {
1879 return _blossom_potential.size();
1883 /// \brief Returns the number of the nodes in the blossom.
1885 /// Returns the number of the nodes in the blossom.
1886 int blossomSize(int k) const {
1887 return _blossom_potential[k].end - _blossom_potential[k].begin;
1890 /// \brief Returns the value of the blossom.
1892 /// Returns the the value of the blossom.
1894 Value blossomValue(int k) const {
1895 return _blossom_potential[k].value;
1898 /// \brief Lemon iterator for get the items of the blossom.
1900 /// Lemon iterator for get the nodes of the blossom. This class
1901 /// provides a common style lemon iterator which gives back a
1902 /// subset of the nodes.
1906 /// \brief Constructor.
1908 /// Constructor for get the nodes of the variable.
1909 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1910 : _algorithm(&algorithm)
1912 _index = _algorithm->_blossom_potential[variable].begin;
1913 _last = _algorithm->_blossom_potential[variable].end;
1916 /// \brief Invalid constructor.
1918 /// Invalid constructor.
1919 BlossomIt(Invalid) : _index(-1) {}
1921 /// \brief Conversion to node.
1923 /// Conversion to node.
1924 operator Node() const {
1925 return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1928 /// \brief Increment operator.
1930 /// Increment operator.
1931 BlossomIt& operator++() {
1933 if (_index == _last) {
1939 bool operator==(const BlossomIt& it) const {
1940 return _index == it._index;
1942 bool operator!=(const BlossomIt& it) const {
1943 return _index != it._index;
1947 const MaxWeightedMatching* _algorithm;
1956 /// \ingroup matching
1958 /// \brief Weighted perfect matching in general graphs
1960 /// This class provides an efficient implementation of Edmond's
1961 /// maximum weighted perfecr matching algorithm. The implementation
1962 /// is based on extensive use of priority queues and provides
1963 /// \f$O(nm\log(n))\f$ time complexity.
1965 /// The maximum weighted matching problem is to find undirected
1966 /// arcs in the digraph with maximum overall weight and no two of
1967 /// them shares their endpoints and covers all nodes. The problem
1968 /// can be formulated with the next linear program:
1969 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1970 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1971 /// \f[x_e \ge 0\quad \forall e\in E\f]
1972 /// \f[\max \sum_{e\in E}x_ew_e\f]
1973 /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
1974 /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
1975 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1978 /// The algorithm calculates an optimal matching and a proof of the
1979 /// optimality. The solution of the dual problem can be used to check
1980 /// the result of the algorithm. The dual linear problem is the next:
1981 /// \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]
1982 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1983 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1985 /// The algorithm can be executed with \c run() or the \c init() and
1986 /// then the \c start() member functions. After it the matching can
1987 /// be asked with \c matching() or mate() functions. The dual
1988 /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1989 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1990 /// "BlossomIt" nested class which is able to iterate on the nodes
1991 /// of a blossom. If the value type is integral then the dual
1992 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1993 template <typename _Graph,
1994 typename _WeightMap = typename _Graph::template EdgeMap<int> >
1995 class MaxWeightedPerfectMatching {
1998 typedef _Graph Graph;
1999 typedef _WeightMap WeightMap;
2000 typedef typename WeightMap::Value Value;
2002 /// \brief Scaling factor for dual solution
2004 /// Scaling factor for dual solution, it is equal to 4 or 1
2005 /// according to the value type.
2006 static const int dualScale =
2007 std::numeric_limits<Value>::is_integer ? 4 : 1;
2009 typedef typename Graph::template NodeMap<typename Graph::Arc>
2014 TEMPLATE_GRAPH_TYPEDEFS(Graph);
2016 typedef typename Graph::template NodeMap<Value> NodePotential;
2017 typedef std::vector<Node> BlossomNodeList;
2019 struct BlossomVariable {
2023 BlossomVariable(int _begin, int _end, Value _value)
2024 : begin(_begin), end(_end), value(_value) {}
2028 typedef std::vector<BlossomVariable> BlossomPotential;
2030 const Graph& _graph;
2031 const WeightMap& _weight;
2033 MatchingMap* _matching;
2035 NodePotential* _node_potential;
2037 BlossomPotential _blossom_potential;
2038 BlossomNodeList _blossom_node_list;
2043 typedef typename Graph::template NodeMap<int> NodeIntMap;
2044 typedef typename Graph::template ArcMap<int> ArcIntMap;
2045 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
2046 typedef RangeMap<int> IntIntMap;
2049 EVEN = -1, MATCHED = 0, ODD = 1
2052 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
2053 struct BlossomData {
2060 NodeIntMap *_blossom_index;
2061 BlossomSet *_blossom_set;
2062 RangeMap<BlossomData>* _blossom_data;
2064 NodeIntMap *_node_index;
2065 ArcIntMap *_node_heap_index;
2069 NodeData(ArcIntMap& node_heap_index)
2070 : heap(node_heap_index) {}
2074 BinHeap<Value, ArcIntMap> heap;
2075 std::map<int, Arc> heap_index;
2080 RangeMap<NodeData>* _node_data;
2082 typedef ExtendFindEnum<IntIntMap> TreeSet;
2084 IntIntMap *_tree_set_index;
2087 IntIntMap *_delta2_index;
2088 BinHeap<Value, IntIntMap> *_delta2;
2090 EdgeIntMap *_delta3_index;
2091 BinHeap<Value, EdgeIntMap> *_delta3;
2093 IntIntMap *_delta4_index;
2094 BinHeap<Value, IntIntMap> *_delta4;
2098 void createStructures() {
2099 _node_num = countNodes(_graph);
2100 _blossom_num = _node_num * 3 / 2;
2103 _matching = new MatchingMap(_graph);
2105 if (!_node_potential) {
2106 _node_potential = new NodePotential(_graph);
2108 if (!_blossom_set) {
2109 _blossom_index = new NodeIntMap(_graph);
2110 _blossom_set = new BlossomSet(*_blossom_index);
2111 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2115 _node_index = new NodeIntMap(_graph);
2116 _node_heap_index = new ArcIntMap(_graph);
2117 _node_data = new RangeMap<NodeData>(_node_num,
2118 NodeData(*_node_heap_index));
2122 _tree_set_index = new IntIntMap(_blossom_num);
2123 _tree_set = new TreeSet(*_tree_set_index);
2126 _delta2_index = new IntIntMap(_blossom_num);
2127 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2130 _delta3_index = new EdgeIntMap(_graph);
2131 _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
2134 _delta4_index = new IntIntMap(_blossom_num);
2135 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2139 void destroyStructures() {
2140 _node_num = countNodes(_graph);
2141 _blossom_num = _node_num * 3 / 2;
2146 if (_node_potential) {
2147 delete _node_potential;
2150 delete _blossom_index;
2151 delete _blossom_set;
2152 delete _blossom_data;
2157 delete _node_heap_index;
2162 delete _tree_set_index;
2166 delete _delta2_index;
2170 delete _delta3_index;
2174 delete _delta4_index;
2179 void matchedToEven(int blossom, int tree) {
2180 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2181 _delta2->erase(blossom);
2184 if (!_blossom_set->trivial(blossom)) {
2185 (*_blossom_data)[blossom].pot -=
2186 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2189 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2190 n != INVALID; ++n) {
2192 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2193 int ni = (*_node_index)[n];
2195 (*_node_data)[ni].heap.clear();
2196 (*_node_data)[ni].heap_index.clear();
2198 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2200 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2201 Node v = _graph.source(e);
2202 int vb = _blossom_set->find(v);
2203 int vi = (*_node_index)[v];
2205 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2206 dualScale * _weight[e];
2208 if ((*_blossom_data)[vb].status == EVEN) {
2209 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2210 _delta3->push(e, rw / 2);
2213 typename std::map<int, Arc>::iterator it =
2214 (*_node_data)[vi].heap_index.find(tree);
2216 if (it != (*_node_data)[vi].heap_index.end()) {
2217 if ((*_node_data)[vi].heap[it->second] > rw) {
2218 (*_node_data)[vi].heap.replace(it->second, e);
2219 (*_node_data)[vi].heap.decrease(e, rw);
2223 (*_node_data)[vi].heap.push(e, rw);
2224 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2227 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2228 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2230 if ((*_blossom_data)[vb].status == MATCHED) {
2231 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2232 _delta2->push(vb, _blossom_set->classPrio(vb) -
2233 (*_blossom_data)[vb].offset);
2234 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2235 (*_blossom_data)[vb].offset){
2236 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2237 (*_blossom_data)[vb].offset);
2244 (*_blossom_data)[blossom].offset = 0;
2247 void matchedToOdd(int blossom) {
2248 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2249 _delta2->erase(blossom);
2251 (*_blossom_data)[blossom].offset += _delta_sum;
2252 if (!_blossom_set->trivial(blossom)) {
2253 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2254 (*_blossom_data)[blossom].offset);
2258 void evenToMatched(int blossom, int tree) {
2259 if (!_blossom_set->trivial(blossom)) {
2260 (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2263 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2264 n != INVALID; ++n) {
2265 int ni = (*_node_index)[n];
2266 (*_node_data)[ni].pot -= _delta_sum;
2268 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2269 Node v = _graph.source(e);
2270 int vb = _blossom_set->find(v);
2271 int vi = (*_node_index)[v];
2273 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2274 dualScale * _weight[e];
2276 if (vb == blossom) {
2277 if (_delta3->state(e) == _delta3->IN_HEAP) {
2280 } else if ((*_blossom_data)[vb].status == EVEN) {
2282 if (_delta3->state(e) == _delta3->IN_HEAP) {
2286 int vt = _tree_set->find(vb);
2290 Arc r = _graph.oppositeArc(e);
2292 typename std::map<int, Arc>::iterator it =
2293 (*_node_data)[ni].heap_index.find(vt);
2295 if (it != (*_node_data)[ni].heap_index.end()) {
2296 if ((*_node_data)[ni].heap[it->second] > rw) {
2297 (*_node_data)[ni].heap.replace(it->second, r);
2298 (*_node_data)[ni].heap.decrease(r, rw);
2302 (*_node_data)[ni].heap.push(r, rw);
2303 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2306 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2307 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2309 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2310 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2311 (*_blossom_data)[blossom].offset);
2312 } else if ((*_delta2)[blossom] >
2313 _blossom_set->classPrio(blossom) -
2314 (*_blossom_data)[blossom].offset){
2315 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2316 (*_blossom_data)[blossom].offset);
2322 typename std::map<int, Arc>::iterator it =
2323 (*_node_data)[vi].heap_index.find(tree);
2325 if (it != (*_node_data)[vi].heap_index.end()) {
2326 (*_node_data)[vi].heap.erase(it->second);
2327 (*_node_data)[vi].heap_index.erase(it);
2328 if ((*_node_data)[vi].heap.empty()) {
2329 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2330 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2331 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2334 if ((*_blossom_data)[vb].status == MATCHED) {
2335 if (_blossom_set->classPrio(vb) ==
2336 std::numeric_limits<Value>::max()) {
2338 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2339 (*_blossom_data)[vb].offset) {
2340 _delta2->increase(vb, _blossom_set->classPrio(vb) -
2341 (*_blossom_data)[vb].offset);
2350 void oddToMatched(int blossom) {
2351 (*_blossom_data)[blossom].offset -= _delta_sum;
2353 if (_blossom_set->classPrio(blossom) !=
2354 std::numeric_limits<Value>::max()) {
2355 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2356 (*_blossom_data)[blossom].offset);
2359 if (!_blossom_set->trivial(blossom)) {
2360 _delta4->erase(blossom);
2364 void oddToEven(int blossom, int tree) {
2365 if (!_blossom_set->trivial(blossom)) {
2366 _delta4->erase(blossom);
2367 (*_blossom_data)[blossom].pot -=
2368 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2371 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2372 n != INVALID; ++n) {
2373 int ni = (*_node_index)[n];
2375 _blossom_set->increase(n, std::numeric_limits<Value>::max());
2377 (*_node_data)[ni].heap.clear();
2378 (*_node_data)[ni].heap_index.clear();
2379 (*_node_data)[ni].pot +=
2380 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2382 for (InArcIt e(_graph, n); e != INVALID; ++e) {
2383 Node v = _graph.source(e);
2384 int vb = _blossom_set->find(v);
2385 int vi = (*_node_index)[v];
2387 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2388 dualScale * _weight[e];
2390 if ((*_blossom_data)[vb].status == EVEN) {
2391 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2392 _delta3->push(e, rw / 2);
2396 typename std::map<int, Arc>::iterator it =
2397 (*_node_data)[vi].heap_index.find(tree);
2399 if (it != (*_node_data)[vi].heap_index.end()) {
2400 if ((*_node_data)[vi].heap[it->second] > rw) {
2401 (*_node_data)[vi].heap.replace(it->second, e);
2402 (*_node_data)[vi].heap.decrease(e, rw);
2406 (*_node_data)[vi].heap.push(e, rw);
2407 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2410 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2411 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2413 if ((*_blossom_data)[vb].status == MATCHED) {
2414 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2415 _delta2->push(vb, _blossom_set->classPrio(vb) -
2416 (*_blossom_data)[vb].offset);
2417 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2418 (*_blossom_data)[vb].offset) {
2419 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2420 (*_blossom_data)[vb].offset);
2427 (*_blossom_data)[blossom].offset = 0;
2430 void alternatePath(int even, int tree) {
2433 evenToMatched(even, tree);
2434 (*_blossom_data)[even].status = MATCHED;
2436 while ((*_blossom_data)[even].pred != INVALID) {
2437 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2438 (*_blossom_data)[odd].status = MATCHED;
2440 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2442 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2443 (*_blossom_data)[even].status = MATCHED;
2444 evenToMatched(even, tree);
2445 (*_blossom_data)[even].next =
2446 _graph.oppositeArc((*_blossom_data)[odd].pred);
2451 void destroyTree(int tree) {
2452 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2453 if ((*_blossom_data)[b].status == EVEN) {
2454 (*_blossom_data)[b].status = MATCHED;
2455 evenToMatched(b, tree);
2456 } else if ((*_blossom_data)[b].status == ODD) {
2457 (*_blossom_data)[b].status = MATCHED;
2461 _tree_set->eraseClass(tree);
2464 void augmentOnArc(const Edge& arc) {
2466 int left = _blossom_set->find(_graph.u(arc));
2467 int right = _blossom_set->find(_graph.v(arc));
2469 int left_tree = _tree_set->find(left);
2470 alternatePath(left, left_tree);
2471 destroyTree(left_tree);
2473 int right_tree = _tree_set->find(right);
2474 alternatePath(right, right_tree);
2475 destroyTree(right_tree);
2477 (*_blossom_data)[left].next = _graph.direct(arc, true);
2478 (*_blossom_data)[right].next = _graph.direct(arc, false);
2481 void extendOnArc(const Arc& arc) {
2482 int base = _blossom_set->find(_graph.target(arc));
2483 int tree = _tree_set->find(base);
2485 int odd = _blossom_set->find(_graph.source(arc));
2486 _tree_set->insert(odd, tree);
2487 (*_blossom_data)[odd].status = ODD;
2489 (*_blossom_data)[odd].pred = arc;
2491 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2492 (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2493 _tree_set->insert(even, tree);
2494 (*_blossom_data)[even].status = EVEN;
2495 matchedToEven(even, tree);
2498 void shrinkOnArc(const Edge& edge, int tree) {
2500 std::vector<int> left_path, right_path;
2503 std::set<int> left_set, right_set;
2504 int left = _blossom_set->find(_graph.u(edge));
2505 left_path.push_back(left);
2506 left_set.insert(left);
2508 int right = _blossom_set->find(_graph.v(edge));
2509 right_path.push_back(right);
2510 right_set.insert(right);
2514 if ((*_blossom_data)[left].pred == INVALID) break;
2517 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2518 left_path.push_back(left);
2520 _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2521 left_path.push_back(left);
2523 left_set.insert(left);
2525 if (right_set.find(left) != right_set.end()) {
2530 if ((*_blossom_data)[right].pred == INVALID) break;
2533 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2534 right_path.push_back(right);
2536 _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2537 right_path.push_back(right);
2539 right_set.insert(right);
2541 if (left_set.find(right) != left_set.end()) {
2549 if ((*_blossom_data)[left].pred == INVALID) {
2551 while (left_set.find(nca) == left_set.end()) {
2553 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2554 right_path.push_back(nca);
2556 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2557 right_path.push_back(nca);
2561 while (right_set.find(nca) == right_set.end()) {
2563 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2564 left_path.push_back(nca);
2566 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2567 left_path.push_back(nca);
2573 std::vector<int> subblossoms;
2576 prev = _graph.direct(edge, true);
2577 for (int i = 0; left_path[i] != nca; i += 2) {
2578 subblossoms.push_back(left_path[i]);
2579 (*_blossom_data)[left_path[i]].next = prev;
2580 _tree_set->erase(left_path[i]);
2582 subblossoms.push_back(left_path[i + 1]);
2583 (*_blossom_data)[left_path[i + 1]].status = EVEN;
2584 oddToEven(left_path[i + 1], tree);
2585 _tree_set->erase(left_path[i + 1]);
2586 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2590 while (right_path[k] != nca) ++k;
2592 subblossoms.push_back(nca);
2593 (*_blossom_data)[nca].next = prev;
2595 for (int i = k - 2; i >= 0; i -= 2) {
2596 subblossoms.push_back(right_path[i + 1]);
2597 (*_blossom_data)[right_path[i + 1]].status = EVEN;
2598 oddToEven(right_path[i + 1], tree);
2599 _tree_set->erase(right_path[i + 1]);
2601 (*_blossom_data)[right_path[i + 1]].next =
2602 (*_blossom_data)[right_path[i + 1]].pred;
2604 subblossoms.push_back(right_path[i]);
2605 _tree_set->erase(right_path[i]);
2609 _blossom_set->join(subblossoms.begin(), subblossoms.end());
2611 for (int i = 0; i < int(subblossoms.size()); ++i) {
2612 if (!_blossom_set->trivial(subblossoms[i])) {
2613 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2615 (*_blossom_data)[subblossoms[i]].status = MATCHED;
2618 (*_blossom_data)[surface].pot = -2 * _delta_sum;
2619 (*_blossom_data)[surface].offset = 0;
2620 (*_blossom_data)[surface].status = EVEN;
2621 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2622 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2624 _tree_set->insert(surface, tree);
2625 _tree_set->erase(nca);
2628 void splitBlossom(int blossom) {
2629 Arc next = (*_blossom_data)[blossom].next;
2630 Arc pred = (*_blossom_data)[blossom].pred;
2632 int tree = _tree_set->find(blossom);
2634 (*_blossom_data)[blossom].status = MATCHED;
2635 oddToMatched(blossom);
2636 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2637 _delta2->erase(blossom);
2640 std::vector<int> subblossoms;
2641 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2643 Value offset = (*_blossom_data)[blossom].offset;
2644 int b = _blossom_set->find(_graph.source(pred));
2645 int d = _blossom_set->find(_graph.source(next));
2647 int ib = -1, id = -1;
2648 for (int i = 0; i < int(subblossoms.size()); ++i) {
2649 if (subblossoms[i] == b) ib = i;
2650 if (subblossoms[i] == d) id = i;
2652 (*_blossom_data)[subblossoms[i]].offset = offset;
2653 if (!_blossom_set->trivial(subblossoms[i])) {
2654 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2656 if (_blossom_set->classPrio(subblossoms[i]) !=
2657 std::numeric_limits<Value>::max()) {
2658 _delta2->push(subblossoms[i],
2659 _blossom_set->classPrio(subblossoms[i]) -
2660 (*_blossom_data)[subblossoms[i]].offset);
2664 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2665 for (int i = (id + 1) % subblossoms.size();
2666 i != ib; i = (i + 2) % subblossoms.size()) {
2667 int sb = subblossoms[i];
2668 int tb = subblossoms[(i + 1) % subblossoms.size()];
2669 (*_blossom_data)[sb].next =
2670 _graph.oppositeArc((*_blossom_data)[tb].next);
2673 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2674 int sb = subblossoms[i];
2675 int tb = subblossoms[(i + 1) % subblossoms.size()];
2676 int ub = subblossoms[(i + 2) % subblossoms.size()];
2678 (*_blossom_data)[sb].status = ODD;
2680 _tree_set->insert(sb, tree);
2681 (*_blossom_data)[sb].pred = pred;
2682 (*_blossom_data)[sb].next =
2683 _graph.oppositeArc((*_blossom_data)[tb].next);
2685 pred = (*_blossom_data)[ub].next;
2687 (*_blossom_data)[tb].status = EVEN;
2688 matchedToEven(tb, tree);
2689 _tree_set->insert(tb, tree);
2690 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2693 (*_blossom_data)[subblossoms[id]].status = ODD;
2694 matchedToOdd(subblossoms[id]);
2695 _tree_set->insert(subblossoms[id], tree);
2696 (*_blossom_data)[subblossoms[id]].next = next;
2697 (*_blossom_data)[subblossoms[id]].pred = pred;
2701 for (int i = (ib + 1) % subblossoms.size();
2702 i != id; i = (i + 2) % subblossoms.size()) {
2703 int sb = subblossoms[i];
2704 int tb = subblossoms[(i + 1) % subblossoms.size()];
2705 (*_blossom_data)[sb].next =
2706 _graph.oppositeArc((*_blossom_data)[tb].next);
2709 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2710 int sb = subblossoms[i];
2711 int tb = subblossoms[(i + 1) % subblossoms.size()];
2712 int ub = subblossoms[(i + 2) % subblossoms.size()];
2714 (*_blossom_data)[sb].status = ODD;
2716 _tree_set->insert(sb, tree);
2717 (*_blossom_data)[sb].next = next;
2718 (*_blossom_data)[sb].pred =
2719 _graph.oppositeArc((*_blossom_data)[tb].next);
2721 (*_blossom_data)[tb].status = EVEN;
2722 matchedToEven(tb, tree);
2723 _tree_set->insert(tb, tree);
2724 (*_blossom_data)[tb].pred =
2725 (*_blossom_data)[tb].next =
2726 _graph.oppositeArc((*_blossom_data)[ub].next);
2727 next = (*_blossom_data)[ub].next;
2730 (*_blossom_data)[subblossoms[ib]].status = ODD;
2731 matchedToOdd(subblossoms[ib]);
2732 _tree_set->insert(subblossoms[ib], tree);
2733 (*_blossom_data)[subblossoms[ib]].next = next;
2734 (*_blossom_data)[subblossoms[ib]].pred = pred;
2736 _tree_set->erase(blossom);
2739 void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2740 if (_blossom_set->trivial(blossom)) {
2741 int bi = (*_node_index)[base];
2742 Value pot = (*_node_data)[bi].pot;
2744 _matching->set(base, matching);
2745 _blossom_node_list.push_back(base);
2746 _node_potential->set(base, pot);
2749 Value pot = (*_blossom_data)[blossom].pot;
2750 int bn = _blossom_node_list.size();
2752 std::vector<int> subblossoms;
2753 _blossom_set->split(blossom, std::back_inserter(subblossoms));
2754 int b = _blossom_set->find(base);
2756 for (int i = 0; i < int(subblossoms.size()); ++i) {
2757 if (subblossoms[i] == b) { ib = i; break; }
2760 for (int i = 1; i < int(subblossoms.size()); i += 2) {
2761 int sb = subblossoms[(ib + i) % subblossoms.size()];
2762 int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2764 Arc m = (*_blossom_data)[tb].next;
2765 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2766 extractBlossom(tb, _graph.source(m), m);
2768 extractBlossom(subblossoms[ib], base, matching);
2770 int en = _blossom_node_list.size();
2772 _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2776 void extractMatching() {
2777 std::vector<int> blossoms;
2778 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2779 blossoms.push_back(c);
2782 for (int i = 0; i < int(blossoms.size()); ++i) {
2784 Value offset = (*_blossom_data)[blossoms[i]].offset;
2785 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2786 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2787 n != INVALID; ++n) {
2788 (*_node_data)[(*_node_index)[n]].pot -= offset;
2791 Arc matching = (*_blossom_data)[blossoms[i]].next;
2792 Node base = _graph.source(matching);
2793 extractBlossom(blossoms[i], base, matching);
2799 /// \brief Constructor
2802 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2803 : _graph(graph), _weight(weight), _matching(0),
2804 _node_potential(0), _blossom_potential(), _blossom_node_list(),
2805 _node_num(0), _blossom_num(0),
2807 _blossom_index(0), _blossom_set(0), _blossom_data(0),
2808 _node_index(0), _node_heap_index(0), _node_data(0),
2809 _tree_set_index(0), _tree_set(0),
2811 _delta2_index(0), _delta2(0),
2812 _delta3_index(0), _delta3(0),
2813 _delta4_index(0), _delta4(0),
2817 ~MaxWeightedPerfectMatching() {
2818 destroyStructures();
2821 /// \name Execution control
2822 /// The simplest way to execute the algorithm is to use the member
2823 /// \c run() member function.
2827 /// \brief Initialize the algorithm
2829 /// Initialize the algorithm
2833 for (ArcIt e(_graph); e != INVALID; ++e) {
2834 _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
2836 for (EdgeIt e(_graph); e != INVALID; ++e) {
2837 _delta3_index->set(e, _delta3->PRE_HEAP);
2839 for (int i = 0; i < _blossom_num; ++i) {
2840 _delta2_index->set(i, _delta2->PRE_HEAP);
2841 _delta4_index->set(i, _delta4->PRE_HEAP);
2845 for (NodeIt n(_graph); n != INVALID; ++n) {
2846 Value max = - std::numeric_limits<Value>::max();
2847 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2848 if (_graph.target(e) == n) continue;
2849 if ((dualScale * _weight[e]) / 2 > max) {
2850 max = (dualScale * _weight[e]) / 2;
2853 _node_index->set(n, index);
2854 (*_node_data)[index].pot = max;
2856 _blossom_set->insert(n, std::numeric_limits<Value>::max());
2858 _tree_set->insert(blossom);
2860 (*_blossom_data)[blossom].status = EVEN;
2861 (*_blossom_data)[blossom].pred = INVALID;
2862 (*_blossom_data)[blossom].next = INVALID;
2863 (*_blossom_data)[blossom].pot = 0;
2864 (*_blossom_data)[blossom].offset = 0;
2867 for (EdgeIt e(_graph); e != INVALID; ++e) {
2868 int si = (*_node_index)[_graph.u(e)];
2869 int ti = (*_node_index)[_graph.v(e)];
2870 if (_graph.u(e) != _graph.v(e)) {
2871 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2872 dualScale * _weight[e]) / 2);
2877 /// \brief Starts the algorithm
2879 /// Starts the algorithm
2885 int unmatched = _node_num;
2886 while (unmatched > 0) {
2887 Value d2 = !_delta2->empty() ?
2888 _delta2->prio() : std::numeric_limits<Value>::max();
2890 Value d3 = !_delta3->empty() ?
2891 _delta3->prio() : std::numeric_limits<Value>::max();
2893 Value d4 = !_delta4->empty() ?
2894 _delta4->prio() : std::numeric_limits<Value>::max();
2896 _delta_sum = d2; OpType ot = D2;
2897 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2898 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2900 if (_delta_sum == std::numeric_limits<Value>::max()) {
2907 int blossom = _delta2->top();
2908 Node n = _blossom_set->classTop(blossom);
2909 Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2915 Edge e = _delta3->top();
2917 int left_blossom = _blossom_set->find(_graph.u(e));
2918 int right_blossom = _blossom_set->find(_graph.v(e));
2920 if (left_blossom == right_blossom) {
2923 int left_tree = _tree_set->find(left_blossom);
2924 int right_tree = _tree_set->find(right_blossom);
2926 if (left_tree == right_tree) {
2927 shrinkOnArc(e, left_tree);
2935 splitBlossom(_delta4->top());
2943 /// \brief Runs %MaxWeightedPerfectMatching algorithm.
2945 /// This method runs the %MaxWeightedPerfectMatching algorithm.
2947 /// \note mwm.run() is just a shortcut of the following code.
2959 /// \name Primal solution
2960 /// Functions for get the primal solution, ie. the matching.
2964 /// \brief Returns the matching value.
2966 /// Returns the matching value.
2967 Value matchingValue() const {
2969 for (NodeIt n(_graph); n != INVALID; ++n) {
2970 if ((*_matching)[n] != INVALID) {
2971 sum += _weight[(*_matching)[n]];
2977 /// \brief Returns true when the arc is in the matching.
2979 /// Returns true when the arc is in the matching.
2980 bool matching(const Edge& arc) const {
2981 return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
2984 /// \brief Returns the incident matching arc.
2986 /// Returns the incident matching arc from given node.
2987 Arc matching(const Node& node) const {
2988 return (*_matching)[node];
2991 /// \brief Returns the mate of the node.
2993 /// Returns the adjancent node in a mathcing arc.
2994 Node mate(const Node& node) const {
2995 return _graph.target((*_matching)[node]);
3000 /// \name Dual solution
3001 /// Functions for get the dual solution.
3005 /// \brief Returns the value of the dual solution.
3007 /// Returns the value of the dual solution. It should be equal to
3008 /// the primal value scaled by \ref dualScale "dual scale".
3009 Value dualValue() const {
3011 for (NodeIt n(_graph); n != INVALID; ++n) {
3012 sum += nodeValue(n);
3014 for (int i = 0; i < blossomNum(); ++i) {
3015 sum += blossomValue(i) * (blossomSize(i) / 2);
3020 /// \brief Returns the value of the node.
3022 /// Returns the the value of the node.
3023 Value nodeValue(const Node& n) const {
3024 return (*_node_potential)[n];
3027 /// \brief Returns the number of the blossoms in the basis.
3029 /// Returns the number of the blossoms in the basis.
3031 int blossomNum() const {
3032 return _blossom_potential.size();
3036 /// \brief Returns the number of the nodes in the blossom.
3038 /// Returns the number of the nodes in the blossom.
3039 int blossomSize(int k) const {
3040 return _blossom_potential[k].end - _blossom_potential[k].begin;
3043 /// \brief Returns the value of the blossom.
3045 /// Returns the the value of the blossom.
3047 Value blossomValue(int k) const {
3048 return _blossom_potential[k].value;
3051 /// \brief Lemon iterator for get the items of the blossom.
3053 /// Lemon iterator for get the nodes of the blossom. This class
3054 /// provides a common style lemon iterator which gives back a
3055 /// subset of the nodes.
3059 /// \brief Constructor.
3061 /// Constructor for get the nodes of the variable.
3062 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3063 : _algorithm(&algorithm)
3065 _index = _algorithm->_blossom_potential[variable].begin;
3066 _last = _algorithm->_blossom_potential[variable].end;
3069 /// \brief Invalid constructor.
3071 /// Invalid constructor.
3072 BlossomIt(Invalid) : _index(-1) {}
3074 /// \brief Conversion to node.
3076 /// Conversion to node.
3077 operator Node() const {
3078 return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
3081 /// \brief Increment operator.
3083 /// Increment operator.
3084 BlossomIt& operator++() {
3086 if (_index == _last) {
3092 bool operator==(const BlossomIt& it) const {
3093 return _index == it._index;
3095 bool operator!=(const BlossomIt& it) const {
3096 return _index != it._index;
3100 const MaxWeightedPerfectMatching* _algorithm;
3110 } //END OF NAMESPACE LEMON
3112 #endif //LEMON_MAX_MATCHING_H