1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
2  *
3  * This file is a part of LEMON, a generic C++ optimization library.
4  *
5  * Copyright (C) 2003-2013
6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
8  *
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.
12  *
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
15  * purpose.
16  *
17  */
18 
19 #ifndef LEMON_MATCHING_H
20 #define LEMON_MATCHING_H
21 
22 #include <vector>
23 #include <queue>
24 #include <set>
25 #include <limits>
26 
27 #include <lemon/core.h>
28 #include <lemon/unionfind.h>
29 #include <lemon/bin_heap.h>
30 #include <lemon/maps.h>
31 #include <lemon/fractional_matching.h>
32 
33 ///\ingroup matching
34 ///\file
35 ///\brief Maximum matching algorithms in general graphs.
36 
37 namespace lemon {
38 
39   /// \ingroup matching
40   ///
41   /// \brief Maximum cardinality matching in general graphs
42   ///
43   /// This class implements Edmonds' alternating forest matching algorithm
44   /// for finding a maximum cardinality matching in a general undirected graph.
45   /// It can be started from an arbitrary initial matching
46   /// (the default is the empty one).
47   ///
48   /// The dual solution of the problem is a map of the nodes to
49   /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
50   /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
51   /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
52   /// with factor-critical components, the nodes in \c ODD/A form the
53   /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
54   /// a perfect matching. The number of the factor-critical components
55   /// minus the number of barrier nodes is a lower bound on the
56   /// unmatched nodes, and the matching is optimal if and only if this bound is
57   /// tight. This decomposition can be obtained using \ref status() or
58   /// \ref statusMap() after running the algorithm.
59   ///
60   /// \tparam GR The undirected graph type the algorithm runs on.
61   template <typename GR>
62   class MaxMatching {
63   public:
64 
65     /// The graph type of the algorithm
66     typedef GR Graph;
67     /// The type of the matching map
68     typedef typename Graph::template NodeMap<typename Graph::Arc>
69     MatchingMap;
70 
71     ///\brief Status constants for Gallai-Edmonds decomposition.
72     ///
73     ///These constants are used for indicating the Gallai-Edmonds
74     ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
75     ///induce a subgraph with factor-critical components, the nodes with
76     ///status \c ODD (or \c A) form the canonical barrier, and the nodes
77     ///with status \c MATCHED (or \c C) induce a subgraph having a
78     ///perfect matching.
79     enum Status {
80       EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
81       D = 1,
82       MATCHED = 0,    ///< = 0. (\c C is an alias for \c MATCHED.)
83       C = 0,
84       ODD = -1,       ///< = -1. (\c A is an alias for \c ODD.)
85       A = -1,
86       UNMATCHED = -2  ///< = -2.
87     };
88 
89     /// The type of the status map
90     typedef typename Graph::template NodeMap<Status> StatusMap;
91 
92   private:
93 
94     TEMPLATE_GRAPH_TYPEDEFS(Graph);
95 
96     typedef UnionFindEnum<IntNodeMap> BlossomSet;
97     typedef ExtendFindEnum<IntNodeMap> TreeSet;
98     typedef RangeMap<Node> NodeIntMap;
99     typedef MatchingMap EarMap;
100     typedef std::vector<Node> NodeQueue;
101 
102     const Graph& _graph;
103     MatchingMap* _matching;
104     StatusMap* _status;
105 
106     EarMap* _ear;
107 
108     IntNodeMap* _blossom_set_index;
109     BlossomSet* _blossom_set;
110     NodeIntMap* _blossom_rep;
111 
112     IntNodeMap* _tree_set_index;
113     TreeSet* _tree_set;
114 
115     NodeQueue _node_queue;
116     int _process, _postpone, _last;
117 
118     int _node_num;
119 
120   private:
121 
createStructures()122     void createStructures() {
123       _node_num = countNodes(_graph);
124       if (!_matching) {
125         _matching = new MatchingMap(_graph);
126       }
127       if (!_status) {
128         _status = new StatusMap(_graph);
129       }
130       if (!_ear) {
131         _ear = new EarMap(_graph);
132       }
133       if (!_blossom_set) {
134         _blossom_set_index = new IntNodeMap(_graph);
135         _blossom_set = new BlossomSet(*_blossom_set_index);
136       }
137       if (!_blossom_rep) {
138         _blossom_rep = new NodeIntMap(_node_num);
139       }
140       if (!_tree_set) {
141         _tree_set_index = new IntNodeMap(_graph);
142         _tree_set = new TreeSet(*_tree_set_index);
143       }
144       _node_queue.resize(_node_num);
145     }
146 
destroyStructures()147     void destroyStructures() {
148       if (_matching) {
149         delete _matching;
150       }
151       if (_status) {
152         delete _status;
153       }
154       if (_ear) {
155         delete _ear;
156       }
157       if (_blossom_set) {
158         delete _blossom_set;
159         delete _blossom_set_index;
160       }
161       if (_blossom_rep) {
162         delete _blossom_rep;
163       }
164       if (_tree_set) {
165         delete _tree_set_index;
166         delete _tree_set;
167       }
168     }
169 
processDense(const Node & n)170     void processDense(const Node& n) {
171       _process = _postpone = _last = 0;
172       _node_queue[_last++] = n;
173 
174       while (_process != _last) {
175         Node u = _node_queue[_process++];
176         for (OutArcIt a(_graph, u); a != INVALID; ++a) {
177           Node v = _graph.target(a);
178           if ((*_status)[v] == MATCHED) {
179             extendOnArc(a);
180           } else if ((*_status)[v] == UNMATCHED) {
181             augmentOnArc(a);
182             return;
183           }
184         }
185       }
186 
187       while (_postpone != _last) {
188         Node u = _node_queue[_postpone++];
189 
190         for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
191           Node v = _graph.target(a);
192 
193           if ((*_status)[v] == EVEN) {
194             if (_blossom_set->find(u) != _blossom_set->find(v)) {
195               shrinkOnEdge(a);
196             }
197           }
198 
199           while (_process != _last) {
200             Node w = _node_queue[_process++];
201             for (OutArcIt b(_graph, w); b != INVALID; ++b) {
202               Node x = _graph.target(b);
203               if ((*_status)[x] == MATCHED) {
204                 extendOnArc(b);
205               } else if ((*_status)[x] == UNMATCHED) {
206                 augmentOnArc(b);
207                 return;
208               }
209             }
210           }
211         }
212       }
213     }
214 
processSparse(const Node & n)215     void processSparse(const Node& n) {
216       _process = _last = 0;
217       _node_queue[_last++] = n;
218       while (_process != _last) {
219         Node u = _node_queue[_process++];
220         for (OutArcIt a(_graph, u); a != INVALID; ++a) {
221           Node v = _graph.target(a);
222 
223           if ((*_status)[v] == EVEN) {
224             if (_blossom_set->find(u) != _blossom_set->find(v)) {
225               shrinkOnEdge(a);
226             }
227           } else if ((*_status)[v] == MATCHED) {
228             extendOnArc(a);
229           } else if ((*_status)[v] == UNMATCHED) {
230             augmentOnArc(a);
231             return;
232           }
233         }
234       }
235     }
236 
shrinkOnEdge(const Edge & e)237     void shrinkOnEdge(const Edge& e) {
238       Node nca = INVALID;
239 
240       {
241         std::set<Node> left_set, right_set;
242 
243         Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
244         left_set.insert(left);
245 
246         Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
247         right_set.insert(right);
248 
249         while (true) {
250           if ((*_matching)[left] == INVALID) break;
251           left = _graph.target((*_matching)[left]);
252           left = (*_blossom_rep)[_blossom_set->
253                                  find(_graph.target((*_ear)[left]))];
254           if (right_set.find(left) != right_set.end()) {
255             nca = left;
256             break;
257           }
258           left_set.insert(left);
259 
260           if ((*_matching)[right] == INVALID) break;
261           right = _graph.target((*_matching)[right]);
262           right = (*_blossom_rep)[_blossom_set->
263                                   find(_graph.target((*_ear)[right]))];
264           if (left_set.find(right) != left_set.end()) {
265             nca = right;
266             break;
267           }
268           right_set.insert(right);
269         }
270 
271         if (nca == INVALID) {
272           if ((*_matching)[left] == INVALID) {
273             nca = right;
274             while (left_set.find(nca) == left_set.end()) {
275               nca = _graph.target((*_matching)[nca]);
276               nca =(*_blossom_rep)[_blossom_set->
277                                    find(_graph.target((*_ear)[nca]))];
278             }
279           } else {
280             nca = left;
281             while (right_set.find(nca) == right_set.end()) {
282               nca = _graph.target((*_matching)[nca]);
283               nca = (*_blossom_rep)[_blossom_set->
284                                    find(_graph.target((*_ear)[nca]))];
285             }
286           }
287         }
288       }
289 
290       {
291 
292         Node node = _graph.u(e);
293         Arc arc = _graph.direct(e, true);
294         Node base = (*_blossom_rep)[_blossom_set->find(node)];
295 
296         while (base != nca) {
297           (*_ear)[node] = arc;
298 
299           Node n = node;
300           while (n != base) {
301             n = _graph.target((*_matching)[n]);
302             Arc a = (*_ear)[n];
303             n = _graph.target(a);
304             (*_ear)[n] = _graph.oppositeArc(a);
305           }
306           node = _graph.target((*_matching)[base]);
307           _tree_set->erase(base);
308           _tree_set->erase(node);
309           _blossom_set->insert(node, _blossom_set->find(base));
310           (*_status)[node] = EVEN;
311           _node_queue[_last++] = node;
312           arc = _graph.oppositeArc((*_ear)[node]);
313           node = _graph.target((*_ear)[node]);
314           base = (*_blossom_rep)[_blossom_set->find(node)];
315           _blossom_set->join(_graph.target(arc), base);
316         }
317       }
318 
319       (*_blossom_rep)[_blossom_set->find(nca)] = nca;
320 
321       {
322 
323         Node node = _graph.v(e);
324         Arc arc = _graph.direct(e, false);
325         Node base = (*_blossom_rep)[_blossom_set->find(node)];
326 
327         while (base != nca) {
328           (*_ear)[node] = arc;
329 
330           Node n = node;
331           while (n != base) {
332             n = _graph.target((*_matching)[n]);
333             Arc a = (*_ear)[n];
334             n = _graph.target(a);
335             (*_ear)[n] = _graph.oppositeArc(a);
336           }
337           node = _graph.target((*_matching)[base]);
338           _tree_set->erase(base);
339           _tree_set->erase(node);
340           _blossom_set->insert(node, _blossom_set->find(base));
341           (*_status)[node] = EVEN;
342           _node_queue[_last++] = node;
343           arc = _graph.oppositeArc((*_ear)[node]);
344           node = _graph.target((*_ear)[node]);
345           base = (*_blossom_rep)[_blossom_set->find(node)];
346           _blossom_set->join(_graph.target(arc), base);
347         }
348       }
349 
350       (*_blossom_rep)[_blossom_set->find(nca)] = nca;
351     }
352 
extendOnArc(const Arc & a)353     void extendOnArc(const Arc& a) {
354       Node base = _graph.source(a);
355       Node odd = _graph.target(a);
356 
357       (*_ear)[odd] = _graph.oppositeArc(a);
358       Node even = _graph.target((*_matching)[odd]);
359       (*_blossom_rep)[_blossom_set->insert(even)] = even;
360       (*_status)[odd] = ODD;
361       (*_status)[even] = EVEN;
362       int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
363       _tree_set->insert(odd, tree);
364       _tree_set->insert(even, tree);
365       _node_queue[_last++] = even;
366 
367     }
368 
augmentOnArc(const Arc & a)369     void augmentOnArc(const Arc& a) {
370       Node even = _graph.source(a);
371       Node odd = _graph.target(a);
372 
373       int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
374 
375       (*_matching)[odd] = _graph.oppositeArc(a);
376       (*_status)[odd] = MATCHED;
377 
378       Arc arc = (*_matching)[even];
379       (*_matching)[even] = a;
380 
381       while (arc != INVALID) {
382         odd = _graph.target(arc);
383         arc = (*_ear)[odd];
384         even = _graph.target(arc);
385         (*_matching)[odd] = arc;
386         arc = (*_matching)[even];
387         (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
388       }
389 
390       for (typename TreeSet::ItemIt it(*_tree_set, tree);
391            it != INVALID; ++it) {
392         if ((*_status)[it] == ODD) {
393           (*_status)[it] = MATCHED;
394         } else {
395           int blossom = _blossom_set->find(it);
396           for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
397                jt != INVALID; ++jt) {
398             (*_status)[jt] = MATCHED;
399           }
400           _blossom_set->eraseClass(blossom);
401         }
402       }
403       _tree_set->eraseClass(tree);
404 
405     }
406 
407   public:
408 
409     /// \brief Constructor
410     ///
411     /// Constructor.
MaxMatching(const Graph & graph)412     MaxMatching(const Graph& graph)
413       : _graph(graph), _matching(0), _status(0), _ear(0),
414         _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
415         _tree_set_index(0), _tree_set(0) {}
416 
~MaxMatching()417     ~MaxMatching() {
418       destroyStructures();
419     }
420 
421     /// \name Execution Control
422     /// The simplest way to execute the algorithm is to use the
423     /// \c run() member function.\n
424     /// If you need better control on the execution, you have to call
425     /// one of the functions \ref init(), \ref greedyInit() or
426     /// \ref matchingInit() first, then you can start the algorithm with
427     /// \ref startSparse() or \ref startDense().
428 
429     ///@{
430 
431     /// \brief Set the initial matching to the empty matching.
432     ///
433     /// This function sets the initial matching to the empty matching.
init()434     void init() {
435       createStructures();
436       for(NodeIt n(_graph); n != INVALID; ++n) {
437         (*_matching)[n] = INVALID;
438         (*_status)[n] = UNMATCHED;
439       }
440     }
441 
442     /// \brief Find an initial matching in a greedy way.
443     ///
444     /// This function finds an initial matching in a greedy way.
greedyInit()445     void greedyInit() {
446       createStructures();
447       for (NodeIt n(_graph); n != INVALID; ++n) {
448         (*_matching)[n] = INVALID;
449         (*_status)[n] = UNMATCHED;
450       }
451       for (NodeIt n(_graph); n != INVALID; ++n) {
452         if ((*_matching)[n] == INVALID) {
453           for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
454             Node v = _graph.target(a);
455             if ((*_matching)[v] == INVALID && v != n) {
456               (*_matching)[n] = a;
457               (*_status)[n] = MATCHED;
458               (*_matching)[v] = _graph.oppositeArc(a);
459               (*_status)[v] = MATCHED;
460               break;
461             }
462           }
463         }
464       }
465     }
466 
467 
468     /// \brief Initialize the matching from a map.
469     ///
470     /// This function initializes the matching from a \c bool valued edge
471     /// map. This map should have the property that there are no two incident
472     /// edges with \c true value, i.e. it really contains a matching.
473     /// \return \c true if the map contains a matching.
474     template <typename MatchingMap>
matchingInit(const MatchingMap & matching)475     bool matchingInit(const MatchingMap& matching) {
476       createStructures();
477 
478       for (NodeIt n(_graph); n != INVALID; ++n) {
479         (*_matching)[n] = INVALID;
480         (*_status)[n] = UNMATCHED;
481       }
482       for(EdgeIt e(_graph); e!=INVALID; ++e) {
483         if (matching[e]) {
484 
485           Node u = _graph.u(e);
486           if ((*_matching)[u] != INVALID) return false;
487           (*_matching)[u] = _graph.direct(e, true);
488           (*_status)[u] = MATCHED;
489 
490           Node v = _graph.v(e);
491           if ((*_matching)[v] != INVALID) return false;
492           (*_matching)[v] = _graph.direct(e, false);
493           (*_status)[v] = MATCHED;
494         }
495       }
496       return true;
497     }
498 
499     /// \brief Start Edmonds' algorithm
500     ///
501     /// This function runs the original Edmonds' algorithm.
502     ///
503     /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
504     /// called before using this function.
startSparse()505     void startSparse() {
506       for(NodeIt n(_graph); n != INVALID; ++n) {
507         if ((*_status)[n] == UNMATCHED) {
508           (*_blossom_rep)[_blossom_set->insert(n)] = n;
509           _tree_set->insert(n);
510           (*_status)[n] = EVEN;
511           processSparse(n);
512         }
513       }
514     }
515 
516     /// \brief Start Edmonds' algorithm with a heuristic improvement
517     /// for dense graphs
518     ///
519     /// This function runs Edmonds' algorithm with a heuristic of postponing
520     /// shrinks, therefore resulting in a faster algorithm for dense graphs.
521     ///
522     /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
523     /// called before using this function.
startDense()524     void startDense() {
525       for(NodeIt n(_graph); n != INVALID; ++n) {
526         if ((*_status)[n] == UNMATCHED) {
527           (*_blossom_rep)[_blossom_set->insert(n)] = n;
528           _tree_set->insert(n);
529           (*_status)[n] = EVEN;
530           processDense(n);
531         }
532       }
533     }
534 
535 
536     /// \brief Run Edmonds' algorithm
537     ///
538     /// This function runs Edmonds' algorithm. An additional heuristic of
539     /// postponing shrinks is used for relatively dense graphs
540     /// (for which <tt>m>=2*n</tt> holds).
run()541     void run() {
542       if (countEdges(_graph) < 2 * countNodes(_graph)) {
543         greedyInit();
544         startSparse();
545       } else {
546         init();
547         startDense();
548       }
549     }
550 
551     /// @}
552 
553     /// \name Primal Solution
554     /// Functions to get the primal solution, i.e. the maximum matching.
555 
556     /// @{
557 
558     /// \brief Return the size (cardinality) of the matching.
559     ///
560     /// This function returns the size (cardinality) of the current matching.
561     /// After run() it returns the size of the maximum matching in the graph.
matchingSize()562     int matchingSize() const {
563       int size = 0;
564       for (NodeIt n(_graph); n != INVALID; ++n) {
565         if ((*_matching)[n] != INVALID) {
566           ++size;
567         }
568       }
569       return size / 2;
570     }
571 
572     /// \brief Return \c true if the given edge is in the matching.
573     ///
574     /// This function returns \c true if the given edge is in the current
575     /// matching.
matching(const Edge & edge)576     bool matching(const Edge& edge) const {
577       return edge == (*_matching)[_graph.u(edge)];
578     }
579 
580     /// \brief Return the matching arc (or edge) incident to the given node.
581     ///
582     /// This function returns the matching arc (or edge) incident to the
583     /// given node in the current matching or \c INVALID if the node is
584     /// not covered by the matching.
matching(const Node & n)585     Arc matching(const Node& n) const {
586       return (*_matching)[n];
587     }
588 
589     /// \brief Return a const reference to the matching map.
590     ///
591     /// This function returns a const reference to a node map that stores
592     /// the matching arc (or edge) incident to each node.
matchingMap()593     const MatchingMap& matchingMap() const {
594       return *_matching;
595     }
596 
597     /// \brief Return the mate of the given node.
598     ///
599     /// This function returns the mate of the given node in the current
600     /// matching or \c INVALID if the node is not covered by the matching.
mate(const Node & n)601     Node mate(const Node& n) const {
602       return (*_matching)[n] != INVALID ?
603         _graph.target((*_matching)[n]) : INVALID;
604     }
605 
606     /// @}
607 
608     /// \name Dual Solution
609     /// Functions to get the dual solution, i.e. the Gallai-Edmonds
610     /// decomposition.
611 
612     /// @{
613 
614     /// \brief Return the status of the given node in the Edmonds-Gallai
615     /// decomposition.
616     ///
617     /// This function returns the \ref Status "status" of the given node
618     /// in the Edmonds-Gallai decomposition.
status(const Node & n)619     Status status(const Node& n) const {
620       return (*_status)[n];
621     }
622 
623     /// \brief Return a const reference to the status map, which stores
624     /// the Edmonds-Gallai decomposition.
625     ///
626     /// This function returns a const reference to a node map that stores the
627     /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
statusMap()628     const StatusMap& statusMap() const {
629       return *_status;
630     }
631 
632     /// \brief Return \c true if the given node is in the barrier.
633     ///
634     /// This function returns \c true if the given node is in the barrier.
barrier(const Node & n)635     bool barrier(const Node& n) const {
636       return (*_status)[n] == ODD;
637     }
638 
639     /// @}
640 
641   };
642 
643   /// \ingroup matching
644   ///
645   /// \brief Weighted matching in general graphs
646   ///
647   /// This class provides an efficient implementation of Edmond's
648   /// maximum weighted matching algorithm. The implementation is based
649   /// on extensive use of priority queues and provides
650   /// \f$O(nm\log n)\f$ time complexity.
651   ///
652   /// The maximum weighted matching problem is to find a subset of the
653   /// edges in an undirected graph with maximum overall weight for which
654   /// each node has at most one incident edge.
655   /// It can be formulated with the following linear program.
656   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
657   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
658       \quad \forall B\in\mathcal{O}\f] */
659   /// \f[x_e \ge 0\quad \forall e\in E\f]
660   /// \f[\max \sum_{e\in E}x_ew_e\f]
661   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
662   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
663   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
664   /// subsets of the nodes.
665   ///
666   /// The algorithm calculates an optimal matching and a proof of the
667   /// optimality. The solution of the dual problem can be used to check
668   /// the result of the algorithm. The dual linear problem is the
669   /// following.
670   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
671       z_B \ge w_{uv} \quad \forall uv\in E\f] */
672   /// \f[y_u \ge 0 \quad \forall u \in V\f]
673   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
674   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
675       \frac{\vert B \vert - 1}{2}z_B\f] */
676   ///
677   /// The algorithm can be executed with the run() function.
678   /// After it the matching (the primal solution) and the dual solution
679   /// can be obtained using the query functions and the
680   /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
681   /// which is able to iterate on the nodes of a blossom.
682   /// If the value type is integer, then the dual solution is multiplied
683   /// by \ref MaxWeightedMatching::dualScale "4".
684   ///
685   /// \tparam GR The undirected graph type the algorithm runs on.
686   /// \tparam WM The type edge weight map. The default type is
687   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
688 #ifdef DOXYGEN
689   template <typename GR, typename WM>
690 #else
691   template <typename GR,
692             typename WM = typename GR::template EdgeMap<int> >
693 #endif
694   class MaxWeightedMatching {
695   public:
696 
697     /// The graph type of the algorithm
698     typedef GR Graph;
699     /// The type of the edge weight map
700     typedef WM WeightMap;
701     /// The value type of the edge weights
702     typedef typename WeightMap::Value Value;
703 
704     /// The type of the matching map
705     typedef typename Graph::template NodeMap<typename Graph::Arc>
706     MatchingMap;
707 
708     /// \brief Scaling factor for dual solution
709     ///
710     /// Scaling factor for dual solution. It is equal to 4 or 1
711     /// according to the value type.
712     static const int dualScale =
713       std::numeric_limits<Value>::is_integer ? 4 : 1;
714 
715   private:
716 
717     TEMPLATE_GRAPH_TYPEDEFS(Graph);
718 
719     typedef typename Graph::template NodeMap<Value> NodePotential;
720     typedef std::vector<Node> BlossomNodeList;
721 
722     struct BlossomVariable {
723       int begin, end;
724       Value value;
725 
BlossomVariableBlossomVariable726       BlossomVariable(int _begin, int _end, Value _value)
727         : begin(_begin), end(_end), value(_value) {}
728 
729     };
730 
731     typedef std::vector<BlossomVariable> BlossomPotential;
732 
733     const Graph& _graph;
734     const WeightMap& _weight;
735 
736     MatchingMap* _matching;
737 
738     NodePotential* _node_potential;
739 
740     BlossomPotential _blossom_potential;
741     BlossomNodeList _blossom_node_list;
742 
743     int _node_num;
744     int _blossom_num;
745 
746     typedef RangeMap<int> IntIntMap;
747 
748     enum Status {
749       EVEN = -1, MATCHED = 0, ODD = 1
750     };
751 
752     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
753     struct BlossomData {
754       int tree;
755       Status status;
756       Arc pred, next;
757       Value pot, offset;
758       Node base;
759     };
760 
761     IntNodeMap *_blossom_index;
762     BlossomSet *_blossom_set;
763     RangeMap<BlossomData>* _blossom_data;
764 
765     IntNodeMap *_node_index;
766     IntArcMap *_node_heap_index;
767 
768     struct NodeData {
769 
NodeDataNodeData770       NodeData(IntArcMap& node_heap_index)
771         : heap(node_heap_index) {}
772 
773       int blossom;
774       Value pot;
775       BinHeap<Value, IntArcMap> heap;
776       std::map<int, Arc> heap_index;
777 
778       int tree;
779     };
780 
781     RangeMap<NodeData>* _node_data;
782 
783     typedef ExtendFindEnum<IntIntMap> TreeSet;
784 
785     IntIntMap *_tree_set_index;
786     TreeSet *_tree_set;
787 
788     IntNodeMap *_delta1_index;
789     BinHeap<Value, IntNodeMap> *_delta1;
790 
791     IntIntMap *_delta2_index;
792     BinHeap<Value, IntIntMap> *_delta2;
793 
794     IntEdgeMap *_delta3_index;
795     BinHeap<Value, IntEdgeMap> *_delta3;
796 
797     IntIntMap *_delta4_index;
798     BinHeap<Value, IntIntMap> *_delta4;
799 
800     Value _delta_sum;
801     int _unmatched;
802 
803     typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
804     FractionalMatching *_fractional;
805 
createStructures()806     void createStructures() {
807       _node_num = countNodes(_graph);
808       _blossom_num = _node_num * 3 / 2;
809 
810       if (!_matching) {
811         _matching = new MatchingMap(_graph);
812       }
813 
814       if (!_node_potential) {
815         _node_potential = new NodePotential(_graph);
816       }
817 
818       if (!_blossom_set) {
819         _blossom_index = new IntNodeMap(_graph);
820         _blossom_set = new BlossomSet(*_blossom_index);
821         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
822       } else if (_blossom_data->size() != _blossom_num) {
823         delete _blossom_data;
824         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
825       }
826 
827       if (!_node_index) {
828         _node_index = new IntNodeMap(_graph);
829         _node_heap_index = new IntArcMap(_graph);
830         _node_data = new RangeMap<NodeData>(_node_num,
831                                             NodeData(*_node_heap_index));
832       } else {
833         delete _node_data;
834         _node_data = new RangeMap<NodeData>(_node_num,
835                                             NodeData(*_node_heap_index));
836       }
837 
838       if (!_tree_set) {
839         _tree_set_index = new IntIntMap(_blossom_num);
840         _tree_set = new TreeSet(*_tree_set_index);
841       } else {
842         _tree_set_index->resize(_blossom_num);
843       }
844 
845       if (!_delta1) {
846         _delta1_index = new IntNodeMap(_graph);
847         _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
848       }
849 
850       if (!_delta2) {
851         _delta2_index = new IntIntMap(_blossom_num);
852         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
853       } else {
854         _delta2_index->resize(_blossom_num);
855       }
856 
857       if (!_delta3) {
858         _delta3_index = new IntEdgeMap(_graph);
859         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
860       }
861 
862       if (!_delta4) {
863         _delta4_index = new IntIntMap(_blossom_num);
864         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
865       } else {
866         _delta4_index->resize(_blossom_num);
867       }
868     }
869 
destroyStructures()870     void destroyStructures() {
871       if (_matching) {
872         delete _matching;
873       }
874       if (_node_potential) {
875         delete _node_potential;
876       }
877       if (_blossom_set) {
878         delete _blossom_index;
879         delete _blossom_set;
880         delete _blossom_data;
881       }
882 
883       if (_node_index) {
884         delete _node_index;
885         delete _node_heap_index;
886         delete _node_data;
887       }
888 
889       if (_tree_set) {
890         delete _tree_set_index;
891         delete _tree_set;
892       }
893       if (_delta1) {
894         delete _delta1_index;
895         delete _delta1;
896       }
897       if (_delta2) {
898         delete _delta2_index;
899         delete _delta2;
900       }
901       if (_delta3) {
902         delete _delta3_index;
903         delete _delta3;
904       }
905       if (_delta4) {
906         delete _delta4_index;
907         delete _delta4;
908       }
909     }
910 
matchedToEven(int blossom,int tree)911     void matchedToEven(int blossom, int tree) {
912       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
913         _delta2->erase(blossom);
914       }
915 
916       if (!_blossom_set->trivial(blossom)) {
917         (*_blossom_data)[blossom].pot -=
918           2 * (_delta_sum - (*_blossom_data)[blossom].offset);
919       }
920 
921       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
922            n != INVALID; ++n) {
923 
924         _blossom_set->increase(n, std::numeric_limits<Value>::max());
925         int ni = (*_node_index)[n];
926 
927         (*_node_data)[ni].heap.clear();
928         (*_node_data)[ni].heap_index.clear();
929 
930         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
931 
932         _delta1->push(n, (*_node_data)[ni].pot);
933 
934         for (InArcIt e(_graph, n); e != INVALID; ++e) {
935           Node v = _graph.source(e);
936           int vb = _blossom_set->find(v);
937           int vi = (*_node_index)[v];
938 
939           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
940             dualScale * _weight[e];
941 
942           if ((*_blossom_data)[vb].status == EVEN) {
943             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
944               _delta3->push(e, rw / 2);
945             }
946           } else {
947             typename std::map<int, Arc>::iterator it =
948               (*_node_data)[vi].heap_index.find(tree);
949 
950             if (it != (*_node_data)[vi].heap_index.end()) {
951               if ((*_node_data)[vi].heap[it->second] > rw) {
952                 (*_node_data)[vi].heap.replace(it->second, e);
953                 (*_node_data)[vi].heap.decrease(e, rw);
954                 it->second = e;
955               }
956             } else {
957               (*_node_data)[vi].heap.push(e, rw);
958               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
959             }
960 
961             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
962               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
963 
964               if ((*_blossom_data)[vb].status == MATCHED) {
965                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
966                   _delta2->push(vb, _blossom_set->classPrio(vb) -
967                                (*_blossom_data)[vb].offset);
968                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
969                            (*_blossom_data)[vb].offset) {
970                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
971                                    (*_blossom_data)[vb].offset);
972                 }
973               }
974             }
975           }
976         }
977       }
978       (*_blossom_data)[blossom].offset = 0;
979     }
980 
matchedToOdd(int blossom)981     void matchedToOdd(int blossom) {
982       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
983         _delta2->erase(blossom);
984       }
985       (*_blossom_data)[blossom].offset += _delta_sum;
986       if (!_blossom_set->trivial(blossom)) {
987         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
988                       (*_blossom_data)[blossom].offset);
989       }
990     }
991 
evenToMatched(int blossom,int tree)992     void evenToMatched(int blossom, int tree) {
993       if (!_blossom_set->trivial(blossom)) {
994         (*_blossom_data)[blossom].pot += 2 * _delta_sum;
995       }
996 
997       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
998            n != INVALID; ++n) {
999         int ni = (*_node_index)[n];
1000         (*_node_data)[ni].pot -= _delta_sum;
1001 
1002         _delta1->erase(n);
1003 
1004         for (InArcIt e(_graph, n); e != INVALID; ++e) {
1005           Node v = _graph.source(e);
1006           int vb = _blossom_set->find(v);
1007           int vi = (*_node_index)[v];
1008 
1009           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1010             dualScale * _weight[e];
1011 
1012           if (vb == blossom) {
1013             if (_delta3->state(e) == _delta3->IN_HEAP) {
1014               _delta3->erase(e);
1015             }
1016           } else if ((*_blossom_data)[vb].status == EVEN) {
1017 
1018             if (_delta3->state(e) == _delta3->IN_HEAP) {
1019               _delta3->erase(e);
1020             }
1021 
1022             int vt = _tree_set->find(vb);
1023 
1024             if (vt != tree) {
1025 
1026               Arc r = _graph.oppositeArc(e);
1027 
1028               typename std::map<int, Arc>::iterator it =
1029                 (*_node_data)[ni].heap_index.find(vt);
1030 
1031               if (it != (*_node_data)[ni].heap_index.end()) {
1032                 if ((*_node_data)[ni].heap[it->second] > rw) {
1033                   (*_node_data)[ni].heap.replace(it->second, r);
1034                   (*_node_data)[ni].heap.decrease(r, rw);
1035                   it->second = r;
1036                 }
1037               } else {
1038                 (*_node_data)[ni].heap.push(r, rw);
1039                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1040               }
1041 
1042               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1043                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1044 
1045                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1046                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1047                                (*_blossom_data)[blossom].offset);
1048                 } else if ((*_delta2)[blossom] >
1049                            _blossom_set->classPrio(blossom) -
1050                            (*_blossom_data)[blossom].offset){
1051                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1052                                    (*_blossom_data)[blossom].offset);
1053                 }
1054               }
1055             }
1056           } else {
1057 
1058             typename std::map<int, Arc>::iterator it =
1059               (*_node_data)[vi].heap_index.find(tree);
1060 
1061             if (it != (*_node_data)[vi].heap_index.end()) {
1062               (*_node_data)[vi].heap.erase(it->second);
1063               (*_node_data)[vi].heap_index.erase(it);
1064               if ((*_node_data)[vi].heap.empty()) {
1065                 _blossom_set->increase(v, std::numeric_limits<Value>::max());
1066               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1067                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1068               }
1069 
1070               if ((*_blossom_data)[vb].status == MATCHED) {
1071                 if (_blossom_set->classPrio(vb) ==
1072                     std::numeric_limits<Value>::max()) {
1073                   _delta2->erase(vb);
1074                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1075                            (*_blossom_data)[vb].offset) {
1076                   _delta2->increase(vb, _blossom_set->classPrio(vb) -
1077                                    (*_blossom_data)[vb].offset);
1078                 }
1079               }
1080             }
1081           }
1082         }
1083       }
1084     }
1085 
oddToMatched(int blossom)1086     void oddToMatched(int blossom) {
1087       (*_blossom_data)[blossom].offset -= _delta_sum;
1088 
1089       if (_blossom_set->classPrio(blossom) !=
1090           std::numeric_limits<Value>::max()) {
1091         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1092                       (*_blossom_data)[blossom].offset);
1093       }
1094 
1095       if (!_blossom_set->trivial(blossom)) {
1096         _delta4->erase(blossom);
1097       }
1098     }
1099 
oddToEven(int blossom,int tree)1100     void oddToEven(int blossom, int tree) {
1101       if (!_blossom_set->trivial(blossom)) {
1102         _delta4->erase(blossom);
1103         (*_blossom_data)[blossom].pot -=
1104           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1105       }
1106 
1107       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1108            n != INVALID; ++n) {
1109         int ni = (*_node_index)[n];
1110 
1111         _blossom_set->increase(n, std::numeric_limits<Value>::max());
1112 
1113         (*_node_data)[ni].heap.clear();
1114         (*_node_data)[ni].heap_index.clear();
1115         (*_node_data)[ni].pot +=
1116           2 * _delta_sum - (*_blossom_data)[blossom].offset;
1117 
1118         _delta1->push(n, (*_node_data)[ni].pot);
1119 
1120         for (InArcIt e(_graph, n); e != INVALID; ++e) {
1121           Node v = _graph.source(e);
1122           int vb = _blossom_set->find(v);
1123           int vi = (*_node_index)[v];
1124 
1125           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1126             dualScale * _weight[e];
1127 
1128           if ((*_blossom_data)[vb].status == EVEN) {
1129             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1130               _delta3->push(e, rw / 2);
1131             }
1132           } else {
1133 
1134             typename std::map<int, Arc>::iterator it =
1135               (*_node_data)[vi].heap_index.find(tree);
1136 
1137             if (it != (*_node_data)[vi].heap_index.end()) {
1138               if ((*_node_data)[vi].heap[it->second] > rw) {
1139                 (*_node_data)[vi].heap.replace(it->second, e);
1140                 (*_node_data)[vi].heap.decrease(e, rw);
1141                 it->second = e;
1142               }
1143             } else {
1144               (*_node_data)[vi].heap.push(e, rw);
1145               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1146             }
1147 
1148             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1149               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1150 
1151               if ((*_blossom_data)[vb].status == MATCHED) {
1152                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
1153                   _delta2->push(vb, _blossom_set->classPrio(vb) -
1154                                (*_blossom_data)[vb].offset);
1155                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1156                            (*_blossom_data)[vb].offset) {
1157                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1158                                    (*_blossom_data)[vb].offset);
1159                 }
1160               }
1161             }
1162           }
1163         }
1164       }
1165       (*_blossom_data)[blossom].offset = 0;
1166     }
1167 
alternatePath(int even,int tree)1168     void alternatePath(int even, int tree) {
1169       int odd;
1170 
1171       evenToMatched(even, tree);
1172       (*_blossom_data)[even].status = MATCHED;
1173 
1174       while ((*_blossom_data)[even].pred != INVALID) {
1175         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1176         (*_blossom_data)[odd].status = MATCHED;
1177         oddToMatched(odd);
1178         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1179 
1180         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1181         (*_blossom_data)[even].status = MATCHED;
1182         evenToMatched(even, tree);
1183         (*_blossom_data)[even].next =
1184           _graph.oppositeArc((*_blossom_data)[odd].pred);
1185       }
1186 
1187     }
1188 
destroyTree(int tree)1189     void destroyTree(int tree) {
1190       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1191         if ((*_blossom_data)[b].status == EVEN) {
1192           (*_blossom_data)[b].status = MATCHED;
1193           evenToMatched(b, tree);
1194         } else if ((*_blossom_data)[b].status == ODD) {
1195           (*_blossom_data)[b].status = MATCHED;
1196           oddToMatched(b);
1197         }
1198       }
1199       _tree_set->eraseClass(tree);
1200     }
1201 
1202 
unmatchNode(const Node & node)1203     void unmatchNode(const Node& node) {
1204       int blossom = _blossom_set->find(node);
1205       int tree = _tree_set->find(blossom);
1206 
1207       alternatePath(blossom, tree);
1208       destroyTree(tree);
1209 
1210       (*_blossom_data)[blossom].base = node;
1211       (*_blossom_data)[blossom].next = INVALID;
1212     }
1213 
augmentOnEdge(const Edge & edge)1214     void augmentOnEdge(const Edge& edge) {
1215 
1216       int left = _blossom_set->find(_graph.u(edge));
1217       int right = _blossom_set->find(_graph.v(edge));
1218 
1219       int left_tree = _tree_set->find(left);
1220       alternatePath(left, left_tree);
1221       destroyTree(left_tree);
1222 
1223       int right_tree = _tree_set->find(right);
1224       alternatePath(right, right_tree);
1225       destroyTree(right_tree);
1226 
1227       (*_blossom_data)[left].next = _graph.direct(edge, true);
1228       (*_blossom_data)[right].next = _graph.direct(edge, false);
1229     }
1230 
augmentOnArc(const Arc & arc)1231     void augmentOnArc(const Arc& arc) {
1232 
1233       int left = _blossom_set->find(_graph.source(arc));
1234       int right = _blossom_set->find(_graph.target(arc));
1235 
1236       (*_blossom_data)[left].status = MATCHED;
1237 
1238       int right_tree = _tree_set->find(right);
1239       alternatePath(right, right_tree);
1240       destroyTree(right_tree);
1241 
1242       (*_blossom_data)[left].next = arc;
1243       (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1244     }
1245 
extendOnArc(const Arc & arc)1246     void extendOnArc(const Arc& arc) {
1247       int base = _blossom_set->find(_graph.target(arc));
1248       int tree = _tree_set->find(base);
1249 
1250       int odd = _blossom_set->find(_graph.source(arc));
1251       _tree_set->insert(odd, tree);
1252       (*_blossom_data)[odd].status = ODD;
1253       matchedToOdd(odd);
1254       (*_blossom_data)[odd].pred = arc;
1255 
1256       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1257       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1258       _tree_set->insert(even, tree);
1259       (*_blossom_data)[even].status = EVEN;
1260       matchedToEven(even, tree);
1261     }
1262 
shrinkOnEdge(const Edge & edge,int tree)1263     void shrinkOnEdge(const Edge& edge, int tree) {
1264       int nca = -1;
1265       std::vector<int> left_path, right_path;
1266 
1267       {
1268         std::set<int> left_set, right_set;
1269         int left = _blossom_set->find(_graph.u(edge));
1270         left_path.push_back(left);
1271         left_set.insert(left);
1272 
1273         int right = _blossom_set->find(_graph.v(edge));
1274         right_path.push_back(right);
1275         right_set.insert(right);
1276 
1277         while (true) {
1278 
1279           if ((*_blossom_data)[left].pred == INVALID) break;
1280 
1281           left =
1282             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1283           left_path.push_back(left);
1284           left =
1285             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1286           left_path.push_back(left);
1287 
1288           left_set.insert(left);
1289 
1290           if (right_set.find(left) != right_set.end()) {
1291             nca = left;
1292             break;
1293           }
1294 
1295           if ((*_blossom_data)[right].pred == INVALID) break;
1296 
1297           right =
1298             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1299           right_path.push_back(right);
1300           right =
1301             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1302           right_path.push_back(right);
1303 
1304           right_set.insert(right);
1305 
1306           if (left_set.find(right) != left_set.end()) {
1307             nca = right;
1308             break;
1309           }
1310 
1311         }
1312 
1313         if (nca == -1) {
1314           if ((*_blossom_data)[left].pred == INVALID) {
1315             nca = right;
1316             while (left_set.find(nca) == left_set.end()) {
1317               nca =
1318                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1319               right_path.push_back(nca);
1320               nca =
1321                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1322               right_path.push_back(nca);
1323             }
1324           } else {
1325             nca = left;
1326             while (right_set.find(nca) == right_set.end()) {
1327               nca =
1328                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1329               left_path.push_back(nca);
1330               nca =
1331                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1332               left_path.push_back(nca);
1333             }
1334           }
1335         }
1336       }
1337 
1338       std::vector<int> subblossoms;
1339       Arc prev;
1340 
1341       prev = _graph.direct(edge, true);
1342       for (int i = 0; left_path[i] != nca; i += 2) {
1343         subblossoms.push_back(left_path[i]);
1344         (*_blossom_data)[left_path[i]].next = prev;
1345         _tree_set->erase(left_path[i]);
1346 
1347         subblossoms.push_back(left_path[i + 1]);
1348         (*_blossom_data)[left_path[i + 1]].status = EVEN;
1349         oddToEven(left_path[i + 1], tree);
1350         _tree_set->erase(left_path[i + 1]);
1351         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1352       }
1353 
1354       int k = 0;
1355       while (right_path[k] != nca) ++k;
1356 
1357       subblossoms.push_back(nca);
1358       (*_blossom_data)[nca].next = prev;
1359 
1360       for (int i = k - 2; i >= 0; i -= 2) {
1361         subblossoms.push_back(right_path[i + 1]);
1362         (*_blossom_data)[right_path[i + 1]].status = EVEN;
1363         oddToEven(right_path[i + 1], tree);
1364         _tree_set->erase(right_path[i + 1]);
1365 
1366         (*_blossom_data)[right_path[i + 1]].next =
1367           (*_blossom_data)[right_path[i + 1]].pred;
1368 
1369         subblossoms.push_back(right_path[i]);
1370         _tree_set->erase(right_path[i]);
1371       }
1372 
1373       int surface =
1374         _blossom_set->join(subblossoms.begin(), subblossoms.end());
1375 
1376       for (int i = 0; i < int(subblossoms.size()); ++i) {
1377         if (!_blossom_set->trivial(subblossoms[i])) {
1378           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1379         }
1380         (*_blossom_data)[subblossoms[i]].status = MATCHED;
1381       }
1382 
1383       (*_blossom_data)[surface].pot = -2 * _delta_sum;
1384       (*_blossom_data)[surface].offset = 0;
1385       (*_blossom_data)[surface].status = EVEN;
1386       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1387       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1388 
1389       _tree_set->insert(surface, tree);
1390       _tree_set->erase(nca);
1391     }
1392 
splitBlossom(int blossom)1393     void splitBlossom(int blossom) {
1394       Arc next = (*_blossom_data)[blossom].next;
1395       Arc pred = (*_blossom_data)[blossom].pred;
1396 
1397       int tree = _tree_set->find(blossom);
1398 
1399       (*_blossom_data)[blossom].status = MATCHED;
1400       oddToMatched(blossom);
1401       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1402         _delta2->erase(blossom);
1403       }
1404 
1405       std::vector<int> subblossoms;
1406       _blossom_set->split(blossom, std::back_inserter(subblossoms));
1407 
1408       Value offset = (*_blossom_data)[blossom].offset;
1409       int b = _blossom_set->find(_graph.source(pred));
1410       int d = _blossom_set->find(_graph.source(next));
1411 
1412       int ib = -1, id = -1;
1413       for (int i = 0; i < int(subblossoms.size()); ++i) {
1414         if (subblossoms[i] == b) ib = i;
1415         if (subblossoms[i] == d) id = i;
1416 
1417         (*_blossom_data)[subblossoms[i]].offset = offset;
1418         if (!_blossom_set->trivial(subblossoms[i])) {
1419           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1420         }
1421         if (_blossom_set->classPrio(subblossoms[i]) !=
1422             std::numeric_limits<Value>::max()) {
1423           _delta2->push(subblossoms[i],
1424                         _blossom_set->classPrio(subblossoms[i]) -
1425                         (*_blossom_data)[subblossoms[i]].offset);
1426         }
1427       }
1428 
1429       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1430         for (int i = (id + 1) % subblossoms.size();
1431              i != ib; i = (i + 2) % subblossoms.size()) {
1432           int sb = subblossoms[i];
1433           int tb = subblossoms[(i + 1) % subblossoms.size()];
1434           (*_blossom_data)[sb].next =
1435             _graph.oppositeArc((*_blossom_data)[tb].next);
1436         }
1437 
1438         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1439           int sb = subblossoms[i];
1440           int tb = subblossoms[(i + 1) % subblossoms.size()];
1441           int ub = subblossoms[(i + 2) % subblossoms.size()];
1442 
1443           (*_blossom_data)[sb].status = ODD;
1444           matchedToOdd(sb);
1445           _tree_set->insert(sb, tree);
1446           (*_blossom_data)[sb].pred = pred;
1447           (*_blossom_data)[sb].next =
1448             _graph.oppositeArc((*_blossom_data)[tb].next);
1449 
1450           pred = (*_blossom_data)[ub].next;
1451 
1452           (*_blossom_data)[tb].status = EVEN;
1453           matchedToEven(tb, tree);
1454           _tree_set->insert(tb, tree);
1455           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1456         }
1457 
1458         (*_blossom_data)[subblossoms[id]].status = ODD;
1459         matchedToOdd(subblossoms[id]);
1460         _tree_set->insert(subblossoms[id], tree);
1461         (*_blossom_data)[subblossoms[id]].next = next;
1462         (*_blossom_data)[subblossoms[id]].pred = pred;
1463 
1464       } else {
1465 
1466         for (int i = (ib + 1) % subblossoms.size();
1467              i != id; i = (i + 2) % subblossoms.size()) {
1468           int sb = subblossoms[i];
1469           int tb = subblossoms[(i + 1) % subblossoms.size()];
1470           (*_blossom_data)[sb].next =
1471             _graph.oppositeArc((*_blossom_data)[tb].next);
1472         }
1473 
1474         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1475           int sb = subblossoms[i];
1476           int tb = subblossoms[(i + 1) % subblossoms.size()];
1477           int ub = subblossoms[(i + 2) % subblossoms.size()];
1478 
1479           (*_blossom_data)[sb].status = ODD;
1480           matchedToOdd(sb);
1481           _tree_set->insert(sb, tree);
1482           (*_blossom_data)[sb].next = next;
1483           (*_blossom_data)[sb].pred =
1484             _graph.oppositeArc((*_blossom_data)[tb].next);
1485 
1486           (*_blossom_data)[tb].status = EVEN;
1487           matchedToEven(tb, tree);
1488           _tree_set->insert(tb, tree);
1489           (*_blossom_data)[tb].pred =
1490             (*_blossom_data)[tb].next =
1491             _graph.oppositeArc((*_blossom_data)[ub].next);
1492           next = (*_blossom_data)[ub].next;
1493         }
1494 
1495         (*_blossom_data)[subblossoms[ib]].status = ODD;
1496         matchedToOdd(subblossoms[ib]);
1497         _tree_set->insert(subblossoms[ib], tree);
1498         (*_blossom_data)[subblossoms[ib]].next = next;
1499         (*_blossom_data)[subblossoms[ib]].pred = pred;
1500       }
1501       _tree_set->erase(blossom);
1502     }
1503 
extractBlossom(int blossom,const Node & base,const Arc & matching)1504     void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1505       if (_blossom_set->trivial(blossom)) {
1506         int bi = (*_node_index)[base];
1507         Value pot = (*_node_data)[bi].pot;
1508 
1509         (*_matching)[base] = matching;
1510         _blossom_node_list.push_back(base);
1511         (*_node_potential)[base] = pot;
1512       } else {
1513 
1514         Value pot = (*_blossom_data)[blossom].pot;
1515         int bn = _blossom_node_list.size();
1516 
1517         std::vector<int> subblossoms;
1518         _blossom_set->split(blossom, std::back_inserter(subblossoms));
1519         int b = _blossom_set->find(base);
1520         int ib = -1;
1521         for (int i = 0; i < int(subblossoms.size()); ++i) {
1522           if (subblossoms[i] == b) { ib = i; break; }
1523         }
1524 
1525         for (int i = 1; i < int(subblossoms.size()); i += 2) {
1526           int sb = subblossoms[(ib + i) % subblossoms.size()];
1527           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1528 
1529           Arc m = (*_blossom_data)[tb].next;
1530           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1531           extractBlossom(tb, _graph.source(m), m);
1532         }
1533         extractBlossom(subblossoms[ib], base, matching);
1534 
1535         int en = _blossom_node_list.size();
1536 
1537         _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1538       }
1539     }
1540 
extractMatching()1541     void extractMatching() {
1542       std::vector<int> blossoms;
1543       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1544         blossoms.push_back(c);
1545       }
1546 
1547       for (int i = 0; i < int(blossoms.size()); ++i) {
1548         if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1549 
1550           Value offset = (*_blossom_data)[blossoms[i]].offset;
1551           (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1552           for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1553                n != INVALID; ++n) {
1554             (*_node_data)[(*_node_index)[n]].pot -= offset;
1555           }
1556 
1557           Arc matching = (*_blossom_data)[blossoms[i]].next;
1558           Node base = _graph.source(matching);
1559           extractBlossom(blossoms[i], base, matching);
1560         } else {
1561           Node base = (*_blossom_data)[blossoms[i]].base;
1562           extractBlossom(blossoms[i], base, INVALID);
1563         }
1564       }
1565     }
1566 
1567   public:
1568 
1569     /// \brief Constructor
1570     ///
1571     /// Constructor.
MaxWeightedMatching(const Graph & graph,const WeightMap & weight)1572     MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1573       : _graph(graph), _weight(weight), _matching(0),
1574         _node_potential(0), _blossom_potential(), _blossom_node_list(),
1575         _node_num(0), _blossom_num(0),
1576 
1577         _blossom_index(0), _blossom_set(0), _blossom_data(0),
1578         _node_index(0), _node_heap_index(0), _node_data(0),
1579         _tree_set_index(0), _tree_set(0),
1580 
1581         _delta1_index(0), _delta1(0),
1582         _delta2_index(0), _delta2(0),
1583         _delta3_index(0), _delta3(0),
1584         _delta4_index(0), _delta4(0),
1585 
1586         _delta_sum(), _unmatched(0),
1587 
1588         _fractional(0)
1589     {}
1590 
~MaxWeightedMatching()1591     ~MaxWeightedMatching() {
1592       destroyStructures();
1593       if (_fractional) {
1594         delete _fractional;
1595       }
1596     }
1597 
1598     /// \name Execution Control
1599     /// The simplest way to execute the algorithm is to use the
1600     /// \ref run() member function.
1601 
1602     ///@{
1603 
1604     /// \brief Initialize the algorithm
1605     ///
1606     /// This function initializes the algorithm.
init()1607     void init() {
1608       createStructures();
1609 
1610       _blossom_node_list.clear();
1611       _blossom_potential.clear();
1612 
1613       for (ArcIt e(_graph); e != INVALID; ++e) {
1614         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1615       }
1616       for (NodeIt n(_graph); n != INVALID; ++n) {
1617         (*_delta1_index)[n] = _delta1->PRE_HEAP;
1618       }
1619       for (EdgeIt e(_graph); e != INVALID; ++e) {
1620         (*_delta3_index)[e] = _delta3->PRE_HEAP;
1621       }
1622       for (int i = 0; i < _blossom_num; ++i) {
1623         (*_delta2_index)[i] = _delta2->PRE_HEAP;
1624         (*_delta4_index)[i] = _delta4->PRE_HEAP;
1625       }
1626 
1627       _unmatched = _node_num;
1628 
1629       _delta1->clear();
1630       _delta2->clear();
1631       _delta3->clear();
1632       _delta4->clear();
1633       _blossom_set->clear();
1634       _tree_set->clear();
1635 
1636       int index = 0;
1637       for (NodeIt n(_graph); n != INVALID; ++n) {
1638         Value max = 0;
1639         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1640           if (_graph.target(e) == n) continue;
1641           if ((dualScale * _weight[e]) / 2 > max) {
1642             max = (dualScale * _weight[e]) / 2;
1643           }
1644         }
1645         (*_node_index)[n] = index;
1646         (*_node_data)[index].heap_index.clear();
1647         (*_node_data)[index].heap.clear();
1648         (*_node_data)[index].pot = max;
1649         _delta1->push(n, max);
1650         int blossom =
1651           _blossom_set->insert(n, std::numeric_limits<Value>::max());
1652 
1653         _tree_set->insert(blossom);
1654 
1655         (*_blossom_data)[blossom].status = EVEN;
1656         (*_blossom_data)[blossom].pred = INVALID;
1657         (*_blossom_data)[blossom].next = INVALID;
1658         (*_blossom_data)[blossom].pot = 0;
1659         (*_blossom_data)[blossom].offset = 0;
1660         ++index;
1661       }
1662       for (EdgeIt e(_graph); e != INVALID; ++e) {
1663         int si = (*_node_index)[_graph.u(e)];
1664         int ti = (*_node_index)[_graph.v(e)];
1665         if (_graph.u(e) != _graph.v(e)) {
1666           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1667                             dualScale * _weight[e]) / 2);
1668         }
1669       }
1670     }
1671 
1672     /// \brief Initialize the algorithm with fractional matching
1673     ///
1674     /// This function initializes the algorithm with a fractional
1675     /// matching. This initialization is also called jumpstart heuristic.
fractionalInit()1676     void fractionalInit() {
1677       createStructures();
1678 
1679       _blossom_node_list.clear();
1680       _blossom_potential.clear();
1681 
1682       if (_fractional == 0) {
1683         _fractional = new FractionalMatching(_graph, _weight, false);
1684       }
1685       _fractional->run();
1686 
1687       for (ArcIt e(_graph); e != INVALID; ++e) {
1688         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1689       }
1690       for (NodeIt n(_graph); n != INVALID; ++n) {
1691         (*_delta1_index)[n] = _delta1->PRE_HEAP;
1692       }
1693       for (EdgeIt e(_graph); e != INVALID; ++e) {
1694         (*_delta3_index)[e] = _delta3->PRE_HEAP;
1695       }
1696       for (int i = 0; i < _blossom_num; ++i) {
1697         (*_delta2_index)[i] = _delta2->PRE_HEAP;
1698         (*_delta4_index)[i] = _delta4->PRE_HEAP;
1699       }
1700 
1701       _unmatched = 0;
1702 
1703       _delta1->clear();
1704       _delta2->clear();
1705       _delta3->clear();
1706       _delta4->clear();
1707       _blossom_set->clear();
1708       _tree_set->clear();
1709 
1710       int index = 0;
1711       for (NodeIt n(_graph); n != INVALID; ++n) {
1712         Value pot = _fractional->nodeValue(n);
1713         (*_node_index)[n] = index;
1714         (*_node_data)[index].pot = pot;
1715         (*_node_data)[index].heap_index.clear();
1716         (*_node_data)[index].heap.clear();
1717         int blossom =
1718           _blossom_set->insert(n, std::numeric_limits<Value>::max());
1719 
1720         (*_blossom_data)[blossom].status = MATCHED;
1721         (*_blossom_data)[blossom].pred = INVALID;
1722         (*_blossom_data)[blossom].next = _fractional->matching(n);
1723         if (_fractional->matching(n) == INVALID) {
1724           (*_blossom_data)[blossom].base = n;
1725         }
1726         (*_blossom_data)[blossom].pot = 0;
1727         (*_blossom_data)[blossom].offset = 0;
1728         ++index;
1729       }
1730 
1731       typename Graph::template NodeMap<bool> processed(_graph, false);
1732       for (NodeIt n(_graph); n != INVALID; ++n) {
1733         if (processed[n]) continue;
1734         processed[n] = true;
1735         if (_fractional->matching(n) == INVALID) continue;
1736         int num = 1;
1737         Node v = _graph.target(_fractional->matching(n));
1738         while (n != v) {
1739           processed[v] = true;
1740           v = _graph.target(_fractional->matching(v));
1741           ++num;
1742         }
1743 
1744         if (num % 2 == 1) {
1745           std::vector<int> subblossoms(num);
1746 
1747           subblossoms[--num] = _blossom_set->find(n);
1748           _delta1->push(n, _fractional->nodeValue(n));
1749           v = _graph.target(_fractional->matching(n));
1750           while (n != v) {
1751             subblossoms[--num] = _blossom_set->find(v);
1752             _delta1->push(v, _fractional->nodeValue(v));
1753             v = _graph.target(_fractional->matching(v));
1754           }
1755 
1756           int surface =
1757             _blossom_set->join(subblossoms.begin(), subblossoms.end());
1758           (*_blossom_data)[surface].status = EVEN;
1759           (*_blossom_data)[surface].pred = INVALID;
1760           (*_blossom_data)[surface].next = INVALID;
1761           (*_blossom_data)[surface].pot = 0;
1762           (*_blossom_data)[surface].offset = 0;
1763 
1764           _tree_set->insert(surface);
1765           ++_unmatched;
1766         }
1767       }
1768 
1769       for (EdgeIt e(_graph); e != INVALID; ++e) {
1770         int si = (*_node_index)[_graph.u(e)];
1771         int sb = _blossom_set->find(_graph.u(e));
1772         int ti = (*_node_index)[_graph.v(e)];
1773         int tb = _blossom_set->find(_graph.v(e));
1774         if ((*_blossom_data)[sb].status == EVEN &&
1775             (*_blossom_data)[tb].status == EVEN && sb != tb) {
1776           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1777                             dualScale * _weight[e]) / 2);
1778         }
1779       }
1780 
1781       for (NodeIt n(_graph); n != INVALID; ++n) {
1782         int nb = _blossom_set->find(n);
1783         if ((*_blossom_data)[nb].status != MATCHED) continue;
1784         int ni = (*_node_index)[n];
1785 
1786         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1787           Node v = _graph.target(e);
1788           int vb = _blossom_set->find(v);
1789           int vi = (*_node_index)[v];
1790 
1791           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1792             dualScale * _weight[e];
1793 
1794           if ((*_blossom_data)[vb].status == EVEN) {
1795 
1796             int vt = _tree_set->find(vb);
1797 
1798             typename std::map<int, Arc>::iterator it =
1799               (*_node_data)[ni].heap_index.find(vt);
1800 
1801             if (it != (*_node_data)[ni].heap_index.end()) {
1802               if ((*_node_data)[ni].heap[it->second] > rw) {
1803                 (*_node_data)[ni].heap.replace(it->second, e);
1804                 (*_node_data)[ni].heap.decrease(e, rw);
1805                 it->second = e;
1806               }
1807             } else {
1808               (*_node_data)[ni].heap.push(e, rw);
1809               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1810             }
1811           }
1812         }
1813 
1814         if (!(*_node_data)[ni].heap.empty()) {
1815           _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1816           _delta2->push(nb, _blossom_set->classPrio(nb));
1817         }
1818       }
1819     }
1820 
1821     /// \brief Start the algorithm
1822     ///
1823     /// This function starts the algorithm.
1824     ///
1825     /// \pre \ref init() or \ref fractionalInit() must be called
1826     /// before using this function.
start()1827     void start() {
1828       enum OpType {
1829         D1, D2, D3, D4
1830       };
1831 
1832       while (_unmatched > 0) {
1833         Value d1 = !_delta1->empty() ?
1834           _delta1->prio() : std::numeric_limits<Value>::max();
1835 
1836         Value d2 = !_delta2->empty() ?
1837           _delta2->prio() : std::numeric_limits<Value>::max();
1838 
1839         Value d3 = !_delta3->empty() ?
1840           _delta3->prio() : std::numeric_limits<Value>::max();
1841 
1842         Value d4 = !_delta4->empty() ?
1843           _delta4->prio() : std::numeric_limits<Value>::max();
1844 
1845         _delta_sum = d3; OpType ot = D3;
1846         if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1847         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1848         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1849 
1850         switch (ot) {
1851         case D1:
1852           {
1853             Node n = _delta1->top();
1854             unmatchNode(n);
1855             --_unmatched;
1856           }
1857           break;
1858         case D2:
1859           {
1860             int blossom = _delta2->top();
1861             Node n = _blossom_set->classTop(blossom);
1862             Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1863             if ((*_blossom_data)[blossom].next == INVALID) {
1864               augmentOnArc(a);
1865               --_unmatched;
1866             } else {
1867               extendOnArc(a);
1868             }
1869           }
1870           break;
1871         case D3:
1872           {
1873             Edge e = _delta3->top();
1874 
1875             int left_blossom = _blossom_set->find(_graph.u(e));
1876             int right_blossom = _blossom_set->find(_graph.v(e));
1877 
1878             if (left_blossom == right_blossom) {
1879               _delta3->pop();
1880             } else {
1881               int left_tree = _tree_set->find(left_blossom);
1882               int right_tree = _tree_set->find(right_blossom);
1883 
1884               if (left_tree == right_tree) {
1885                 shrinkOnEdge(e, left_tree);
1886               } else {
1887                 augmentOnEdge(e);
1888                 _unmatched -= 2;
1889               }
1890             }
1891           } break;
1892         case D4:
1893           splitBlossom(_delta4->top());
1894           break;
1895         }
1896       }
1897       extractMatching();
1898     }
1899 
1900     /// \brief Run the algorithm.
1901     ///
1902     /// This method runs the \c %MaxWeightedMatching algorithm.
1903     ///
1904     /// \note mwm.run() is just a shortcut of the following code.
1905     /// \code
1906     ///   mwm.fractionalInit();
1907     ///   mwm.start();
1908     /// \endcode
run()1909     void run() {
1910       fractionalInit();
1911       start();
1912     }
1913 
1914     /// @}
1915 
1916     /// \name Primal Solution
1917     /// Functions to get the primal solution, i.e. the maximum weighted
1918     /// matching.\n
1919     /// Either \ref run() or \ref start() function should be called before
1920     /// using them.
1921 
1922     /// @{
1923 
1924     /// \brief Return the weight of the matching.
1925     ///
1926     /// This function returns the weight of the found matching.
1927     ///
1928     /// \pre Either run() or start() must be called before using this function.
matchingWeight()1929     Value matchingWeight() const {
1930       Value sum = 0;
1931       for (NodeIt n(_graph); n != INVALID; ++n) {
1932         if ((*_matching)[n] != INVALID) {
1933           sum += _weight[(*_matching)[n]];
1934         }
1935       }
1936       return sum / 2;
1937     }
1938 
1939     /// \brief Return the size (cardinality) of the matching.
1940     ///
1941     /// This function returns the size (cardinality) of the found matching.
1942     ///
1943     /// \pre Either run() or start() must be called before using this function.
matchingSize()1944     int matchingSize() const {
1945       int num = 0;
1946       for (NodeIt n(_graph); n != INVALID; ++n) {
1947         if ((*_matching)[n] != INVALID) {
1948           ++num;
1949         }
1950       }
1951       return num /= 2;
1952     }
1953 
1954     /// \brief Return \c true if the given edge is in the matching.
1955     ///
1956     /// This function returns \c true if the given edge is in the found
1957     /// matching.
1958     ///
1959     /// \pre Either run() or start() must be called before using this function.
matching(const Edge & edge)1960     bool matching(const Edge& edge) const {
1961       return edge == (*_matching)[_graph.u(edge)];
1962     }
1963 
1964     /// \brief Return the matching arc (or edge) incident to the given node.
1965     ///
1966     /// This function returns the matching arc (or edge) incident to the
1967     /// given node in the found matching or \c INVALID if the node is
1968     /// not covered by the matching.
1969     ///
1970     /// \pre Either run() or start() must be called before using this function.
matching(const Node & node)1971     Arc matching(const Node& node) const {
1972       return (*_matching)[node];
1973     }
1974 
1975     /// \brief Return a const reference to the matching map.
1976     ///
1977     /// This function returns a const reference to a node map that stores
1978     /// the matching arc (or edge) incident to each node.
matchingMap()1979     const MatchingMap& matchingMap() const {
1980       return *_matching;
1981     }
1982 
1983     /// \brief Return the mate of the given node.
1984     ///
1985     /// This function returns the mate of the given node in the found
1986     /// matching or \c INVALID if the node is not covered by the matching.
1987     ///
1988     /// \pre Either run() or start() must be called before using this function.
mate(const Node & node)1989     Node mate(const Node& node) const {
1990       return (*_matching)[node] != INVALID ?
1991         _graph.target((*_matching)[node]) : INVALID;
1992     }
1993 
1994     /// @}
1995 
1996     /// \name Dual Solution
1997     /// Functions to get the dual solution.\n
1998     /// Either \ref run() or \ref start() function should be called before
1999     /// using them.
2000 
2001     /// @{
2002 
2003     /// \brief Return the value of the dual solution.
2004     ///
2005     /// This function returns the value of the dual solution.
2006     /// It should be equal to the primal value scaled by \ref dualScale
2007     /// "dual scale".
2008     ///
2009     /// \pre Either run() or start() must be called before using this function.
dualValue()2010     Value dualValue() const {
2011       Value sum = 0;
2012       for (NodeIt n(_graph); n != INVALID; ++n) {
2013         sum += nodeValue(n);
2014       }
2015       for (int i = 0; i < blossomNum(); ++i) {
2016         sum += blossomValue(i) * (blossomSize(i) / 2);
2017       }
2018       return sum;
2019     }
2020 
2021     /// \brief Return the dual value (potential) of the given node.
2022     ///
2023     /// This function returns the dual value (potential) of the given node.
2024     ///
2025     /// \pre Either run() or start() must be called before using this function.
nodeValue(const Node & n)2026     Value nodeValue(const Node& n) const {
2027       return (*_node_potential)[n];
2028     }
2029 
2030     /// \brief Return the number of the blossoms in the basis.
2031     ///
2032     /// This function returns the number of the blossoms in the basis.
2033     ///
2034     /// \pre Either run() or start() must be called before using this function.
2035     /// \see BlossomIt
blossomNum()2036     int blossomNum() const {
2037       return _blossom_potential.size();
2038     }
2039 
2040     /// \brief Return the number of the nodes in the given blossom.
2041     ///
2042     /// This function returns the number of the nodes in the given blossom.
2043     ///
2044     /// \pre Either run() or start() must be called before using this function.
2045     /// \see BlossomIt
blossomSize(int k)2046     int blossomSize(int k) const {
2047       return _blossom_potential[k].end - _blossom_potential[k].begin;
2048     }
2049 
2050     /// \brief Return the dual value (ptential) of the given blossom.
2051     ///
2052     /// This function returns the dual value (ptential) of the given blossom.
2053     ///
2054     /// \pre Either run() or start() must be called before using this function.
blossomValue(int k)2055     Value blossomValue(int k) const {
2056       return _blossom_potential[k].value;
2057     }
2058 
2059     /// \brief Iterator for obtaining the nodes of a blossom.
2060     ///
2061     /// This class provides an iterator for obtaining the nodes of the
2062     /// given blossom. It lists a subset of the nodes.
2063     /// Before using this iterator, you must allocate a
2064     /// MaxWeightedMatching class and execute it.
2065     class BlossomIt {
2066     public:
2067 
2068       /// \brief Constructor.
2069       ///
2070       /// Constructor to get the nodes of the given variable.
2071       ///
2072       /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
2073       /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
2074       /// called before initializing this iterator.
BlossomIt(const MaxWeightedMatching & algorithm,int variable)2075       BlossomIt(const MaxWeightedMatching& algorithm, int variable)
2076         : _algorithm(&algorithm)
2077       {
2078         _index = _algorithm->_blossom_potential[variable].begin;
2079         _last = _algorithm->_blossom_potential[variable].end;
2080       }
2081 
2082       /// \brief Conversion to \c Node.
2083       ///
2084       /// Conversion to \c Node.
Node()2085       operator Node() const {
2086         return _algorithm->_blossom_node_list[_index];
2087       }
2088 
2089       /// \brief Increment operator.
2090       ///
2091       /// Increment operator.
2092       BlossomIt& operator++() {
2093         ++_index;
2094         return *this;
2095       }
2096 
2097       /// \brief Validity checking
2098       ///
2099       /// Checks whether the iterator is invalid.
2100       bool operator==(Invalid) const { return _index == _last; }
2101 
2102       /// \brief Validity checking
2103       ///
2104       /// Checks whether the iterator is valid.
2105       bool operator!=(Invalid) const { return _index != _last; }
2106 
2107     private:
2108       const MaxWeightedMatching* _algorithm;
2109       int _last;
2110       int _index;
2111     };
2112 
2113     /// @}
2114 
2115   };
2116 
2117   /// \ingroup matching
2118   ///
2119   /// \brief Weighted perfect matching in general graphs
2120   ///
2121   /// This class provides an efficient implementation of Edmond's
2122   /// maximum weighted perfect matching algorithm. The implementation
2123   /// is based on extensive use of priority queues and provides
2124   /// \f$O(nm\log n)\f$ time complexity.
2125   ///
2126   /// The maximum weighted perfect matching problem is to find a subset of
2127   /// the edges in an undirected graph with maximum overall weight for which
2128   /// each node has exactly one incident edge.
2129   /// It can be formulated with the following linear program.
2130   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2131   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2132       \quad \forall B\in\mathcal{O}\f] */
2133   /// \f[x_e \ge 0\quad \forall e\in E\f]
2134   /// \f[\max \sum_{e\in E}x_ew_e\f]
2135   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2136   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2137   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2138   /// subsets of the nodes.
2139   ///
2140   /// The algorithm calculates an optimal matching and a proof of the
2141   /// optimality. The solution of the dual problem can be used to check
2142   /// the result of the algorithm. The dual linear problem is the
2143   /// following.
2144   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2145       w_{uv} \quad \forall uv\in E\f] */
2146   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2147   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2148       \frac{\vert B \vert - 1}{2}z_B\f] */
2149   ///
2150   /// The algorithm can be executed with the run() function.
2151   /// After it the matching (the primal solution) and the dual solution
2152   /// can be obtained using the query functions and the
2153   /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2154   /// which is able to iterate on the nodes of a blossom.
2155   /// If the value type is integer, then the dual solution is multiplied
2156   /// by \ref MaxWeightedMatching::dualScale "4".
2157   ///
2158   /// \tparam GR The undirected graph type the algorithm runs on.
2159   /// \tparam WM The type edge weight map. The default type is
2160   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2161 #ifdef DOXYGEN
2162   template <typename GR, typename WM>
2163 #else
2164   template <typename GR,
2165             typename WM = typename GR::template EdgeMap<int> >
2166 #endif
2167   class MaxWeightedPerfectMatching {
2168   public:
2169 
2170     /// The graph type of the algorithm
2171     typedef GR Graph;
2172     /// The type of the edge weight map
2173     typedef WM WeightMap;
2174     /// The value type of the edge weights
2175     typedef typename WeightMap::Value Value;
2176 
2177     /// \brief Scaling factor for dual solution
2178     ///
2179     /// Scaling factor for dual solution, it is equal to 4 or 1
2180     /// according to the value type.
2181     static const int dualScale =
2182       std::numeric_limits<Value>::is_integer ? 4 : 1;
2183 
2184     /// The type of the matching map
2185     typedef typename Graph::template NodeMap<typename Graph::Arc>
2186     MatchingMap;
2187 
2188   private:
2189 
2190     TEMPLATE_GRAPH_TYPEDEFS(Graph);
2191 
2192     typedef typename Graph::template NodeMap<Value> NodePotential;
2193     typedef std::vector<Node> BlossomNodeList;
2194 
2195     struct BlossomVariable {
2196       int begin, end;
2197       Value value;
2198 
BlossomVariableBlossomVariable2199       BlossomVariable(int _begin, int _end, Value _value)
2200         : begin(_begin), end(_end), value(_value) {}
2201 
2202     };
2203 
2204     typedef std::vector<BlossomVariable> BlossomPotential;
2205 
2206     const Graph& _graph;
2207     const WeightMap& _weight;
2208 
2209     MatchingMap* _matching;
2210 
2211     NodePotential* _node_potential;
2212 
2213     BlossomPotential _blossom_potential;
2214     BlossomNodeList _blossom_node_list;
2215 
2216     int _node_num;
2217     int _blossom_num;
2218 
2219     typedef RangeMap<int> IntIntMap;
2220 
2221     enum Status {
2222       EVEN = -1, MATCHED = 0, ODD = 1
2223     };
2224 
2225     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2226     struct BlossomData {
2227       int tree;
2228       Status status;
2229       Arc pred, next;
2230       Value pot, offset;
2231     };
2232 
2233     IntNodeMap *_blossom_index;
2234     BlossomSet *_blossom_set;
2235     RangeMap<BlossomData>* _blossom_data;
2236 
2237     IntNodeMap *_node_index;
2238     IntArcMap *_node_heap_index;
2239 
2240     struct NodeData {
2241 
NodeDataNodeData2242       NodeData(IntArcMap& node_heap_index)
2243         : heap(node_heap_index) {}
2244 
2245       int blossom;
2246       Value pot;
2247       BinHeap<Value, IntArcMap> heap;
2248       std::map<int, Arc> heap_index;
2249 
2250       int tree;
2251     };
2252 
2253     RangeMap<NodeData>* _node_data;
2254 
2255     typedef ExtendFindEnum<IntIntMap> TreeSet;
2256 
2257     IntIntMap *_tree_set_index;
2258     TreeSet *_tree_set;
2259 
2260     IntIntMap *_delta2_index;
2261     BinHeap<Value, IntIntMap> *_delta2;
2262 
2263     IntEdgeMap *_delta3_index;
2264     BinHeap<Value, IntEdgeMap> *_delta3;
2265 
2266     IntIntMap *_delta4_index;
2267     BinHeap<Value, IntIntMap> *_delta4;
2268 
2269     Value _delta_sum;
2270     int _unmatched;
2271 
2272     typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
2273     FractionalMatching;
2274     FractionalMatching *_fractional;
2275 
createStructures()2276     void createStructures() {
2277       _node_num = countNodes(_graph);
2278       _blossom_num = _node_num * 3 / 2;
2279 
2280       if (!_matching) {
2281         _matching = new MatchingMap(_graph);
2282       }
2283 
2284       if (!_node_potential) {
2285         _node_potential = new NodePotential(_graph);
2286       }
2287 
2288       if (!_blossom_set) {
2289         _blossom_index = new IntNodeMap(_graph);
2290         _blossom_set = new BlossomSet(*_blossom_index);
2291         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2292       } else if (_blossom_data->size() != _blossom_num) {
2293         delete _blossom_data;
2294         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2295       }
2296 
2297       if (!_node_index) {
2298         _node_index = new IntNodeMap(_graph);
2299         _node_heap_index = new IntArcMap(_graph);
2300         _node_data = new RangeMap<NodeData>(_node_num,
2301                                             NodeData(*_node_heap_index));
2302       } else if (_node_data->size() != _node_num) {
2303         delete _node_data;
2304         _node_data = new RangeMap<NodeData>(_node_num,
2305                                             NodeData(*_node_heap_index));
2306       }
2307 
2308       if (!_tree_set) {
2309         _tree_set_index = new IntIntMap(_blossom_num);
2310         _tree_set = new TreeSet(*_tree_set_index);
2311       } else {
2312         _tree_set_index->resize(_blossom_num);
2313       }
2314 
2315       if (!_delta2) {
2316         _delta2_index = new IntIntMap(_blossom_num);
2317         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2318       } else {
2319         _delta2_index->resize(_blossom_num);
2320       }
2321 
2322       if (!_delta3) {
2323         _delta3_index = new IntEdgeMap(_graph);
2324         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2325       }
2326 
2327       if (!_delta4) {
2328         _delta4_index = new IntIntMap(_blossom_num);
2329         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2330       } else {
2331         _delta4_index->resize(_blossom_num);
2332       }
2333     }
2334 
destroyStructures()2335     void destroyStructures() {
2336       if (_matching) {
2337         delete _matching;
2338       }
2339       if (_node_potential) {
2340         delete _node_potential;
2341       }
2342       if (_blossom_set) {
2343         delete _blossom_index;
2344         delete _blossom_set;
2345         delete _blossom_data;
2346       }
2347 
2348       if (_node_index) {
2349         delete _node_index;
2350         delete _node_heap_index;
2351         delete _node_data;
2352       }
2353 
2354       if (_tree_set) {
2355         delete _tree_set_index;
2356         delete _tree_set;
2357       }
2358       if (_delta2) {
2359         delete _delta2_index;
2360         delete _delta2;
2361       }
2362       if (_delta3) {
2363         delete _delta3_index;
2364         delete _delta3;
2365       }
2366       if (_delta4) {
2367         delete _delta4_index;
2368         delete _delta4;
2369       }
2370     }
2371 
matchedToEven(int blossom,int tree)2372     void matchedToEven(int blossom, int tree) {
2373       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2374         _delta2->erase(blossom);
2375       }
2376 
2377       if (!_blossom_set->trivial(blossom)) {
2378         (*_blossom_data)[blossom].pot -=
2379           2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2380       }
2381 
2382       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2383            n != INVALID; ++n) {
2384 
2385         _blossom_set->increase(n, std::numeric_limits<Value>::max());
2386         int ni = (*_node_index)[n];
2387 
2388         (*_node_data)[ni].heap.clear();
2389         (*_node_data)[ni].heap_index.clear();
2390 
2391         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2392 
2393         for (InArcIt e(_graph, n); e != INVALID; ++e) {
2394           Node v = _graph.source(e);
2395           int vb = _blossom_set->find(v);
2396           int vi = (*_node_index)[v];
2397 
2398           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2399             dualScale * _weight[e];
2400 
2401           if ((*_blossom_data)[vb].status == EVEN) {
2402             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2403               _delta3->push(e, rw / 2);
2404             }
2405           } else {
2406             typename std::map<int, Arc>::iterator it =
2407               (*_node_data)[vi].heap_index.find(tree);
2408 
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);
2413                 it->second = e;
2414               }
2415             } else {
2416               (*_node_data)[vi].heap.push(e, rw);
2417               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2418             }
2419 
2420             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2421               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2422 
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);
2431                 }
2432               }
2433             }
2434           }
2435         }
2436       }
2437       (*_blossom_data)[blossom].offset = 0;
2438     }
2439 
matchedToOdd(int blossom)2440     void matchedToOdd(int blossom) {
2441       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2442         _delta2->erase(blossom);
2443       }
2444       (*_blossom_data)[blossom].offset += _delta_sum;
2445       if (!_blossom_set->trivial(blossom)) {
2446         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2447                      (*_blossom_data)[blossom].offset);
2448       }
2449     }
2450 
evenToMatched(int blossom,int tree)2451     void evenToMatched(int blossom, int tree) {
2452       if (!_blossom_set->trivial(blossom)) {
2453         (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2454       }
2455 
2456       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2457            n != INVALID; ++n) {
2458         int ni = (*_node_index)[n];
2459         (*_node_data)[ni].pot -= _delta_sum;
2460 
2461         for (InArcIt e(_graph, n); e != INVALID; ++e) {
2462           Node v = _graph.source(e);
2463           int vb = _blossom_set->find(v);
2464           int vi = (*_node_index)[v];
2465 
2466           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2467             dualScale * _weight[e];
2468 
2469           if (vb == blossom) {
2470             if (_delta3->state(e) == _delta3->IN_HEAP) {
2471               _delta3->erase(e);
2472             }
2473           } else if ((*_blossom_data)[vb].status == EVEN) {
2474 
2475             if (_delta3->state(e) == _delta3->IN_HEAP) {
2476               _delta3->erase(e);
2477             }
2478 
2479             int vt = _tree_set->find(vb);
2480 
2481             if (vt != tree) {
2482 
2483               Arc r = _graph.oppositeArc(e);
2484 
2485               typename std::map<int, Arc>::iterator it =
2486                 (*_node_data)[ni].heap_index.find(vt);
2487 
2488               if (it != (*_node_data)[ni].heap_index.end()) {
2489                 if ((*_node_data)[ni].heap[it->second] > rw) {
2490                   (*_node_data)[ni].heap.replace(it->second, r);
2491                   (*_node_data)[ni].heap.decrease(r, rw);
2492                   it->second = r;
2493                 }
2494               } else {
2495                 (*_node_data)[ni].heap.push(r, rw);
2496                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2497               }
2498 
2499               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2500                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2501 
2502                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2503                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2504                                (*_blossom_data)[blossom].offset);
2505                 } else if ((*_delta2)[blossom] >
2506                            _blossom_set->classPrio(blossom) -
2507                            (*_blossom_data)[blossom].offset){
2508                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2509                                    (*_blossom_data)[blossom].offset);
2510                 }
2511               }
2512             }
2513           } else {
2514 
2515             typename std::map<int, Arc>::iterator it =
2516               (*_node_data)[vi].heap_index.find(tree);
2517 
2518             if (it != (*_node_data)[vi].heap_index.end()) {
2519               (*_node_data)[vi].heap.erase(it->second);
2520               (*_node_data)[vi].heap_index.erase(it);
2521               if ((*_node_data)[vi].heap.empty()) {
2522                 _blossom_set->increase(v, std::numeric_limits<Value>::max());
2523               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2524                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2525               }
2526 
2527               if ((*_blossom_data)[vb].status == MATCHED) {
2528                 if (_blossom_set->classPrio(vb) ==
2529                     std::numeric_limits<Value>::max()) {
2530                   _delta2->erase(vb);
2531                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2532                            (*_blossom_data)[vb].offset) {
2533                   _delta2->increase(vb, _blossom_set->classPrio(vb) -
2534                                    (*_blossom_data)[vb].offset);
2535                 }
2536               }
2537             }
2538           }
2539         }
2540       }
2541     }
2542 
oddToMatched(int blossom)2543     void oddToMatched(int blossom) {
2544       (*_blossom_data)[blossom].offset -= _delta_sum;
2545 
2546       if (_blossom_set->classPrio(blossom) !=
2547           std::numeric_limits<Value>::max()) {
2548         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2549                        (*_blossom_data)[blossom].offset);
2550       }
2551 
2552       if (!_blossom_set->trivial(blossom)) {
2553         _delta4->erase(blossom);
2554       }
2555     }
2556 
oddToEven(int blossom,int tree)2557     void oddToEven(int blossom, int tree) {
2558       if (!_blossom_set->trivial(blossom)) {
2559         _delta4->erase(blossom);
2560         (*_blossom_data)[blossom].pot -=
2561           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2562       }
2563 
2564       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2565            n != INVALID; ++n) {
2566         int ni = (*_node_index)[n];
2567 
2568         _blossom_set->increase(n, std::numeric_limits<Value>::max());
2569 
2570         (*_node_data)[ni].heap.clear();
2571         (*_node_data)[ni].heap_index.clear();
2572         (*_node_data)[ni].pot +=
2573           2 * _delta_sum - (*_blossom_data)[blossom].offset;
2574 
2575         for (InArcIt e(_graph, n); e != INVALID; ++e) {
2576           Node v = _graph.source(e);
2577           int vb = _blossom_set->find(v);
2578           int vi = (*_node_index)[v];
2579 
2580           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2581             dualScale * _weight[e];
2582 
2583           if ((*_blossom_data)[vb].status == EVEN) {
2584             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2585               _delta3->push(e, rw / 2);
2586             }
2587           } else {
2588 
2589             typename std::map<int, Arc>::iterator it =
2590               (*_node_data)[vi].heap_index.find(tree);
2591 
2592             if (it != (*_node_data)[vi].heap_index.end()) {
2593               if ((*_node_data)[vi].heap[it->second] > rw) {
2594                 (*_node_data)[vi].heap.replace(it->second, e);
2595                 (*_node_data)[vi].heap.decrease(e, rw);
2596                 it->second = e;
2597               }
2598             } else {
2599               (*_node_data)[vi].heap.push(e, rw);
2600               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2601             }
2602 
2603             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2604               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2605 
2606               if ((*_blossom_data)[vb].status == MATCHED) {
2607                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
2608                   _delta2->push(vb, _blossom_set->classPrio(vb) -
2609                                (*_blossom_data)[vb].offset);
2610                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2611                            (*_blossom_data)[vb].offset) {
2612                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2613                                    (*_blossom_data)[vb].offset);
2614                 }
2615               }
2616             }
2617           }
2618         }
2619       }
2620       (*_blossom_data)[blossom].offset = 0;
2621     }
2622 
alternatePath(int even,int tree)2623     void alternatePath(int even, int tree) {
2624       int odd;
2625 
2626       evenToMatched(even, tree);
2627       (*_blossom_data)[even].status = MATCHED;
2628 
2629       while ((*_blossom_data)[even].pred != INVALID) {
2630         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2631         (*_blossom_data)[odd].status = MATCHED;
2632         oddToMatched(odd);
2633         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2634 
2635         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2636         (*_blossom_data)[even].status = MATCHED;
2637         evenToMatched(even, tree);
2638         (*_blossom_data)[even].next =
2639           _graph.oppositeArc((*_blossom_data)[odd].pred);
2640       }
2641 
2642     }
2643 
destroyTree(int tree)2644     void destroyTree(int tree) {
2645       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2646         if ((*_blossom_data)[b].status == EVEN) {
2647           (*_blossom_data)[b].status = MATCHED;
2648           evenToMatched(b, tree);
2649         } else if ((*_blossom_data)[b].status == ODD) {
2650           (*_blossom_data)[b].status = MATCHED;
2651           oddToMatched(b);
2652         }
2653       }
2654       _tree_set->eraseClass(tree);
2655     }
2656 
augmentOnEdge(const Edge & edge)2657     void augmentOnEdge(const Edge& edge) {
2658 
2659       int left = _blossom_set->find(_graph.u(edge));
2660       int right = _blossom_set->find(_graph.v(edge));
2661 
2662       int left_tree = _tree_set->find(left);
2663       alternatePath(left, left_tree);
2664       destroyTree(left_tree);
2665 
2666       int right_tree = _tree_set->find(right);
2667       alternatePath(right, right_tree);
2668       destroyTree(right_tree);
2669 
2670       (*_blossom_data)[left].next = _graph.direct(edge, true);
2671       (*_blossom_data)[right].next = _graph.direct(edge, false);
2672     }
2673 
extendOnArc(const Arc & arc)2674     void extendOnArc(const Arc& arc) {
2675       int base = _blossom_set->find(_graph.target(arc));
2676       int tree = _tree_set->find(base);
2677 
2678       int odd = _blossom_set->find(_graph.source(arc));
2679       _tree_set->insert(odd, tree);
2680       (*_blossom_data)[odd].status = ODD;
2681       matchedToOdd(odd);
2682       (*_blossom_data)[odd].pred = arc;
2683 
2684       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2685       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2686       _tree_set->insert(even, tree);
2687       (*_blossom_data)[even].status = EVEN;
2688       matchedToEven(even, tree);
2689     }
2690 
shrinkOnEdge(const Edge & edge,int tree)2691     void shrinkOnEdge(const Edge& edge, int tree) {
2692       int nca = -1;
2693       std::vector<int> left_path, right_path;
2694 
2695       {
2696         std::set<int> left_set, right_set;
2697         int left = _blossom_set->find(_graph.u(edge));
2698         left_path.push_back(left);
2699         left_set.insert(left);
2700 
2701         int right = _blossom_set->find(_graph.v(edge));
2702         right_path.push_back(right);
2703         right_set.insert(right);
2704 
2705         while (true) {
2706 
2707           if ((*_blossom_data)[left].pred == INVALID) break;
2708 
2709           left =
2710             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2711           left_path.push_back(left);
2712           left =
2713             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2714           left_path.push_back(left);
2715 
2716           left_set.insert(left);
2717 
2718           if (right_set.find(left) != right_set.end()) {
2719             nca = left;
2720             break;
2721           }
2722 
2723           if ((*_blossom_data)[right].pred == INVALID) break;
2724 
2725           right =
2726             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2727           right_path.push_back(right);
2728           right =
2729             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2730           right_path.push_back(right);
2731 
2732           right_set.insert(right);
2733 
2734           if (left_set.find(right) != left_set.end()) {
2735             nca = right;
2736             break;
2737           }
2738 
2739         }
2740 
2741         if (nca == -1) {
2742           if ((*_blossom_data)[left].pred == INVALID) {
2743             nca = right;
2744             while (left_set.find(nca) == left_set.end()) {
2745               nca =
2746                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2747               right_path.push_back(nca);
2748               nca =
2749                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2750               right_path.push_back(nca);
2751             }
2752           } else {
2753             nca = left;
2754             while (right_set.find(nca) == right_set.end()) {
2755               nca =
2756                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2757               left_path.push_back(nca);
2758               nca =
2759                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2760               left_path.push_back(nca);
2761             }
2762           }
2763         }
2764       }
2765 
2766       std::vector<int> subblossoms;
2767       Arc prev;
2768 
2769       prev = _graph.direct(edge, true);
2770       for (int i = 0; left_path[i] != nca; i += 2) {
2771         subblossoms.push_back(left_path[i]);
2772         (*_blossom_data)[left_path[i]].next = prev;
2773         _tree_set->erase(left_path[i]);
2774 
2775         subblossoms.push_back(left_path[i + 1]);
2776         (*_blossom_data)[left_path[i + 1]].status = EVEN;
2777         oddToEven(left_path[i + 1], tree);
2778         _tree_set->erase(left_path[i + 1]);
2779         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2780       }
2781 
2782       int k = 0;
2783       while (right_path[k] != nca) ++k;
2784 
2785       subblossoms.push_back(nca);
2786       (*_blossom_data)[nca].next = prev;
2787 
2788       for (int i = k - 2; i >= 0; i -= 2) {
2789         subblossoms.push_back(right_path[i + 1]);
2790         (*_blossom_data)[right_path[i + 1]].status = EVEN;
2791         oddToEven(right_path[i + 1], tree);
2792         _tree_set->erase(right_path[i + 1]);
2793 
2794         (*_blossom_data)[right_path[i + 1]].next =
2795           (*_blossom_data)[right_path[i + 1]].pred;
2796 
2797         subblossoms.push_back(right_path[i]);
2798         _tree_set->erase(right_path[i]);
2799       }
2800 
2801       int surface =
2802         _blossom_set->join(subblossoms.begin(), subblossoms.end());
2803 
2804       for (int i = 0; i < int(subblossoms.size()); ++i) {
2805         if (!_blossom_set->trivial(subblossoms[i])) {
2806           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2807         }
2808         (*_blossom_data)[subblossoms[i]].status = MATCHED;
2809       }
2810 
2811       (*_blossom_data)[surface].pot = -2 * _delta_sum;
2812       (*_blossom_data)[surface].offset = 0;
2813       (*_blossom_data)[surface].status = EVEN;
2814       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2815       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2816 
2817       _tree_set->insert(surface, tree);
2818       _tree_set->erase(nca);
2819     }
2820 
splitBlossom(int blossom)2821     void splitBlossom(int blossom) {
2822       Arc next = (*_blossom_data)[blossom].next;
2823       Arc pred = (*_blossom_data)[blossom].pred;
2824 
2825       int tree = _tree_set->find(blossom);
2826 
2827       (*_blossom_data)[blossom].status = MATCHED;
2828       oddToMatched(blossom);
2829       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2830         _delta2->erase(blossom);
2831       }
2832 
2833       std::vector<int> subblossoms;
2834       _blossom_set->split(blossom, std::back_inserter(subblossoms));
2835 
2836       Value offset = (*_blossom_data)[blossom].offset;
2837       int b = _blossom_set->find(_graph.source(pred));
2838       int d = _blossom_set->find(_graph.source(next));
2839 
2840       int ib = -1, id = -1;
2841       for (int i = 0; i < int(subblossoms.size()); ++i) {
2842         if (subblossoms[i] == b) ib = i;
2843         if (subblossoms[i] == d) id = i;
2844 
2845         (*_blossom_data)[subblossoms[i]].offset = offset;
2846         if (!_blossom_set->trivial(subblossoms[i])) {
2847           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2848         }
2849         if (_blossom_set->classPrio(subblossoms[i]) !=
2850             std::numeric_limits<Value>::max()) {
2851           _delta2->push(subblossoms[i],
2852                         _blossom_set->classPrio(subblossoms[i]) -
2853                         (*_blossom_data)[subblossoms[i]].offset);
2854         }
2855       }
2856 
2857       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2858         for (int i = (id + 1) % subblossoms.size();
2859              i != ib; i = (i + 2) % subblossoms.size()) {
2860           int sb = subblossoms[i];
2861           int tb = subblossoms[(i + 1) % subblossoms.size()];
2862           (*_blossom_data)[sb].next =
2863             _graph.oppositeArc((*_blossom_data)[tb].next);
2864         }
2865 
2866         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2867           int sb = subblossoms[i];
2868           int tb = subblossoms[(i + 1) % subblossoms.size()];
2869           int ub = subblossoms[(i + 2) % subblossoms.size()];
2870 
2871           (*_blossom_data)[sb].status = ODD;
2872           matchedToOdd(sb);
2873           _tree_set->insert(sb, tree);
2874           (*_blossom_data)[sb].pred = pred;
2875           (*_blossom_data)[sb].next =
2876                            _graph.oppositeArc((*_blossom_data)[tb].next);
2877 
2878           pred = (*_blossom_data)[ub].next;
2879 
2880           (*_blossom_data)[tb].status = EVEN;
2881           matchedToEven(tb, tree);
2882           _tree_set->insert(tb, tree);
2883           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2884         }
2885 
2886         (*_blossom_data)[subblossoms[id]].status = ODD;
2887         matchedToOdd(subblossoms[id]);
2888         _tree_set->insert(subblossoms[id], tree);
2889         (*_blossom_data)[subblossoms[id]].next = next;
2890         (*_blossom_data)[subblossoms[id]].pred = pred;
2891 
2892       } else {
2893 
2894         for (int i = (ib + 1) % subblossoms.size();
2895              i != id; i = (i + 2) % subblossoms.size()) {
2896           int sb = subblossoms[i];
2897           int tb = subblossoms[(i + 1) % subblossoms.size()];
2898           (*_blossom_data)[sb].next =
2899             _graph.oppositeArc((*_blossom_data)[tb].next);
2900         }
2901 
2902         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2903           int sb = subblossoms[i];
2904           int tb = subblossoms[(i + 1) % subblossoms.size()];
2905           int ub = subblossoms[(i + 2) % subblossoms.size()];
2906 
2907           (*_blossom_data)[sb].status = ODD;
2908           matchedToOdd(sb);
2909           _tree_set->insert(sb, tree);
2910           (*_blossom_data)[sb].next = next;
2911           (*_blossom_data)[sb].pred =
2912             _graph.oppositeArc((*_blossom_data)[tb].next);
2913 
2914           (*_blossom_data)[tb].status = EVEN;
2915           matchedToEven(tb, tree);
2916           _tree_set->insert(tb, tree);
2917           (*_blossom_data)[tb].pred =
2918             (*_blossom_data)[tb].next =
2919             _graph.oppositeArc((*_blossom_data)[ub].next);
2920           next = (*_blossom_data)[ub].next;
2921         }
2922 
2923         (*_blossom_data)[subblossoms[ib]].status = ODD;
2924         matchedToOdd(subblossoms[ib]);
2925         _tree_set->insert(subblossoms[ib], tree);
2926         (*_blossom_data)[subblossoms[ib]].next = next;
2927         (*_blossom_data)[subblossoms[ib]].pred = pred;
2928       }
2929       _tree_set->erase(blossom);
2930     }
2931 
extractBlossom(int blossom,const Node & base,const Arc & matching)2932     void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2933       if (_blossom_set->trivial(blossom)) {
2934         int bi = (*_node_index)[base];
2935         Value pot = (*_node_data)[bi].pot;
2936 
2937         (*_matching)[base] = matching;
2938         _blossom_node_list.push_back(base);
2939         (*_node_potential)[base] = pot;
2940       } else {
2941 
2942         Value pot = (*_blossom_data)[blossom].pot;
2943         int bn = _blossom_node_list.size();
2944 
2945         std::vector<int> subblossoms;
2946         _blossom_set->split(blossom, std::back_inserter(subblossoms));
2947         int b = _blossom_set->find(base);
2948         int ib = -1;
2949         for (int i = 0; i < int(subblossoms.size()); ++i) {
2950           if (subblossoms[i] == b) { ib = i; break; }
2951         }
2952 
2953         for (int i = 1; i < int(subblossoms.size()); i += 2) {
2954           int sb = subblossoms[(ib + i) % subblossoms.size()];
2955           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2956 
2957           Arc m = (*_blossom_data)[tb].next;
2958           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2959           extractBlossom(tb, _graph.source(m), m);
2960         }
2961         extractBlossom(subblossoms[ib], base, matching);
2962 
2963         int en = _blossom_node_list.size();
2964 
2965         _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2966       }
2967     }
2968 
extractMatching()2969     void extractMatching() {
2970       std::vector<int> blossoms;
2971       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2972         blossoms.push_back(c);
2973       }
2974 
2975       for (int i = 0; i < int(blossoms.size()); ++i) {
2976 
2977         Value offset = (*_blossom_data)[blossoms[i]].offset;
2978         (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2979         for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2980              n != INVALID; ++n) {
2981           (*_node_data)[(*_node_index)[n]].pot -= offset;
2982         }
2983 
2984         Arc matching = (*_blossom_data)[blossoms[i]].next;
2985         Node base = _graph.source(matching);
2986         extractBlossom(blossoms[i], base, matching);
2987       }
2988     }
2989 
2990   public:
2991 
2992     /// \brief Constructor
2993     ///
2994     /// Constructor.
MaxWeightedPerfectMatching(const Graph & graph,const WeightMap & weight)2995     MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2996       : _graph(graph), _weight(weight), _matching(0),
2997         _node_potential(0), _blossom_potential(), _blossom_node_list(),
2998         _node_num(0), _blossom_num(0),
2999 
3000         _blossom_index(0), _blossom_set(0), _blossom_data(0),
3001         _node_index(0), _node_heap_index(0), _node_data(0),
3002         _tree_set_index(0), _tree_set(0),
3003 
3004         _delta2_index(0), _delta2(0),
3005         _delta3_index(0), _delta3(0),
3006         _delta4_index(0), _delta4(0),
3007 
3008         _delta_sum(), _unmatched(0),
3009 
3010         _fractional(0)
3011     {}
3012 
~MaxWeightedPerfectMatching()3013     ~MaxWeightedPerfectMatching() {
3014       destroyStructures();
3015       if (_fractional) {
3016         delete _fractional;
3017       }
3018     }
3019 
3020     /// \name Execution Control
3021     /// The simplest way to execute the algorithm is to use the
3022     /// \ref run() member function.
3023 
3024     ///@{
3025 
3026     /// \brief Initialize the algorithm
3027     ///
3028     /// This function initializes the algorithm.
init()3029     void init() {
3030       createStructures();
3031 
3032       _blossom_node_list.clear();
3033       _blossom_potential.clear();
3034 
3035       for (ArcIt e(_graph); e != INVALID; ++e) {
3036         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
3037       }
3038       for (EdgeIt e(_graph); e != INVALID; ++e) {
3039         (*_delta3_index)[e] = _delta3->PRE_HEAP;
3040       }
3041       for (int i = 0; i < _blossom_num; ++i) {
3042         (*_delta2_index)[i] = _delta2->PRE_HEAP;
3043         (*_delta4_index)[i] = _delta4->PRE_HEAP;
3044       }
3045 
3046       _unmatched = _node_num;
3047 
3048       _delta2->clear();
3049       _delta3->clear();
3050       _delta4->clear();
3051       _blossom_set->clear();
3052       _tree_set->clear();
3053 
3054       int index = 0;
3055       for (NodeIt n(_graph); n != INVALID; ++n) {
3056         Value max = - std::numeric_limits<Value>::max();
3057         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3058           if (_graph.target(e) == n) continue;
3059           if ((dualScale * _weight[e]) / 2 > max) {
3060             max = (dualScale * _weight[e]) / 2;
3061           }
3062         }
3063         (*_node_index)[n] = index;
3064         (*_node_data)[index].heap_index.clear();
3065         (*_node_data)[index].heap.clear();
3066         (*_node_data)[index].pot = max;
3067         int blossom =
3068           _blossom_set->insert(n, std::numeric_limits<Value>::max());
3069 
3070         _tree_set->insert(blossom);
3071 
3072         (*_blossom_data)[blossom].status = EVEN;
3073         (*_blossom_data)[blossom].pred = INVALID;
3074         (*_blossom_data)[blossom].next = INVALID;
3075         (*_blossom_data)[blossom].pot = 0;
3076         (*_blossom_data)[blossom].offset = 0;
3077         ++index;
3078       }
3079       for (EdgeIt e(_graph); e != INVALID; ++e) {
3080         int si = (*_node_index)[_graph.u(e)];
3081         int ti = (*_node_index)[_graph.v(e)];
3082         if (_graph.u(e) != _graph.v(e)) {
3083           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3084                             dualScale * _weight[e]) / 2);
3085         }
3086       }
3087     }
3088 
3089     /// \brief Initialize the algorithm with fractional matching
3090     ///
3091     /// This function initializes the algorithm with a fractional
3092     /// matching. This initialization is also called jumpstart heuristic.
fractionalInit()3093     void fractionalInit() {
3094       createStructures();
3095 
3096       _blossom_node_list.clear();
3097       _blossom_potential.clear();
3098 
3099       if (_fractional == 0) {
3100         _fractional = new FractionalMatching(_graph, _weight, false);
3101       }
3102       if (!_fractional->run()) {
3103         _unmatched = -1;
3104         return;
3105       }
3106 
3107       for (ArcIt e(_graph); e != INVALID; ++e) {
3108         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
3109       }
3110       for (EdgeIt e(_graph); e != INVALID; ++e) {
3111         (*_delta3_index)[e] = _delta3->PRE_HEAP;
3112       }
3113       for (int i = 0; i < _blossom_num; ++i) {
3114         (*_delta2_index)[i] = _delta2->PRE_HEAP;
3115         (*_delta4_index)[i] = _delta4->PRE_HEAP;
3116       }
3117 
3118       _unmatched = 0;
3119 
3120       _delta2->clear();
3121       _delta3->clear();
3122       _delta4->clear();
3123       _blossom_set->clear();
3124       _tree_set->clear();
3125 
3126       int index = 0;
3127       for (NodeIt n(_graph); n != INVALID; ++n) {
3128         Value pot = _fractional->nodeValue(n);
3129         (*_node_index)[n] = index;
3130         (*_node_data)[index].pot = pot;
3131         (*_node_data)[index].heap_index.clear();
3132         (*_node_data)[index].heap.clear();
3133         int blossom =
3134           _blossom_set->insert(n, std::numeric_limits<Value>::max());
3135 
3136         (*_blossom_data)[blossom].status = MATCHED;
3137         (*_blossom_data)[blossom].pred = INVALID;
3138         (*_blossom_data)[blossom].next = _fractional->matching(n);
3139         (*_blossom_data)[blossom].pot = 0;
3140         (*_blossom_data)[blossom].offset = 0;
3141         ++index;
3142       }
3143 
3144       typename Graph::template NodeMap<bool> processed(_graph, false);
3145       for (NodeIt n(_graph); n != INVALID; ++n) {
3146         if (processed[n]) continue;
3147         processed[n] = true;
3148         if (_fractional->matching(n) == INVALID) continue;
3149         int num = 1;
3150         Node v = _graph.target(_fractional->matching(n));
3151         while (n != v) {
3152           processed[v] = true;
3153           v = _graph.target(_fractional->matching(v));
3154           ++num;
3155         }
3156 
3157         if (num % 2 == 1) {
3158           std::vector<int> subblossoms(num);
3159 
3160           subblossoms[--num] = _blossom_set->find(n);
3161           v = _graph.target(_fractional->matching(n));
3162           while (n != v) {
3163             subblossoms[--num] = _blossom_set->find(v);
3164             v = _graph.target(_fractional->matching(v));
3165           }
3166 
3167           int surface =
3168             _blossom_set->join(subblossoms.begin(), subblossoms.end());
3169           (*_blossom_data)[surface].status = EVEN;
3170           (*_blossom_data)[surface].pred = INVALID;
3171           (*_blossom_data)[surface].next = INVALID;
3172           (*_blossom_data)[surface].pot = 0;
3173           (*_blossom_data)[surface].offset = 0;
3174 
3175           _tree_set->insert(surface);
3176           ++_unmatched;
3177         }
3178       }
3179 
3180       for (EdgeIt e(_graph); e != INVALID; ++e) {
3181         int si = (*_node_index)[_graph.u(e)];
3182         int sb = _blossom_set->find(_graph.u(e));
3183         int ti = (*_node_index)[_graph.v(e)];
3184         int tb = _blossom_set->find(_graph.v(e));
3185         if ((*_blossom_data)[sb].status == EVEN &&
3186             (*_blossom_data)[tb].status == EVEN && sb != tb) {
3187           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3188                             dualScale * _weight[e]) / 2);
3189         }
3190       }
3191 
3192       for (NodeIt n(_graph); n != INVALID; ++n) {
3193         int nb = _blossom_set->find(n);
3194         if ((*_blossom_data)[nb].status != MATCHED) continue;
3195         int ni = (*_node_index)[n];
3196 
3197         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3198           Node v = _graph.target(e);
3199           int vb = _blossom_set->find(v);
3200           int vi = (*_node_index)[v];
3201 
3202           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3203             dualScale * _weight[e];
3204 
3205           if ((*_blossom_data)[vb].status == EVEN) {
3206 
3207             int vt = _tree_set->find(vb);
3208 
3209             typename std::map<int, Arc>::iterator it =
3210               (*_node_data)[ni].heap_index.find(vt);
3211 
3212             if (it != (*_node_data)[ni].heap_index.end()) {
3213               if ((*_node_data)[ni].heap[it->second] > rw) {
3214                 (*_node_data)[ni].heap.replace(it->second, e);
3215                 (*_node_data)[ni].heap.decrease(e, rw);
3216                 it->second = e;
3217               }
3218             } else {
3219               (*_node_data)[ni].heap.push(e, rw);
3220               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
3221             }
3222           }
3223         }
3224 
3225         if (!(*_node_data)[ni].heap.empty()) {
3226           _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3227           _delta2->push(nb, _blossom_set->classPrio(nb));
3228         }
3229       }
3230     }
3231 
3232     /// \brief Start the algorithm
3233     ///
3234     /// This function starts the algorithm.
3235     ///
3236     /// \pre \ref init() or \ref fractionalInit() must be called before
3237     /// using this function.
start()3238     bool start() {
3239       enum OpType {
3240         D2, D3, D4
3241       };
3242 
3243       if (_unmatched == -1) return false;
3244 
3245       while (_unmatched > 0) {
3246         Value d2 = !_delta2->empty() ?
3247           _delta2->prio() : std::numeric_limits<Value>::max();
3248 
3249         Value d3 = !_delta3->empty() ?
3250           _delta3->prio() : std::numeric_limits<Value>::max();
3251 
3252         Value d4 = !_delta4->empty() ?
3253           _delta4->prio() : std::numeric_limits<Value>::max();
3254 
3255         _delta_sum = d3; OpType ot = D3;
3256         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3257         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3258 
3259         if (_delta_sum == std::numeric_limits<Value>::max()) {
3260           return false;
3261         }
3262 
3263         switch (ot) {
3264         case D2:
3265           {
3266             int blossom = _delta2->top();
3267             Node n = _blossom_set->classTop(blossom);
3268             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
3269             extendOnArc(e);
3270           }
3271           break;
3272         case D3:
3273           {
3274             Edge e = _delta3->top();
3275 
3276             int left_blossom = _blossom_set->find(_graph.u(e));
3277             int right_blossom = _blossom_set->find(_graph.v(e));
3278 
3279             if (left_blossom == right_blossom) {
3280               _delta3->pop();
3281             } else {
3282               int left_tree = _tree_set->find(left_blossom);
3283               int right_tree = _tree_set->find(right_blossom);
3284 
3285               if (left_tree == right_tree) {
3286                 shrinkOnEdge(e, left_tree);
3287               } else {
3288                 augmentOnEdge(e);
3289                 _unmatched -= 2;
3290               }
3291             }
3292           } break;
3293         case D4:
3294           splitBlossom(_delta4->top());
3295           break;
3296         }
3297       }
3298       extractMatching();
3299       return true;
3300     }
3301 
3302     /// \brief Run the algorithm.
3303     ///
3304     /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3305     ///
3306     /// \note mwpm.run() is just a shortcut of the following code.
3307     /// \code
3308     ///   mwpm.fractionalInit();
3309     ///   mwpm.start();
3310     /// \endcode
run()3311     bool run() {
3312       fractionalInit();
3313       return start();
3314     }
3315 
3316     /// @}
3317 
3318     /// \name Primal Solution
3319     /// Functions to get the primal solution, i.e. the maximum weighted
3320     /// perfect matching.\n
3321     /// Either \ref run() or \ref start() function should be called before
3322     /// using them.
3323 
3324     /// @{
3325 
3326     /// \brief Return the weight of the matching.
3327     ///
3328     /// This function returns the weight of the found matching.
3329     ///
3330     /// \pre Either run() or start() must be called before using this function.
matchingWeight()3331     Value matchingWeight() const {
3332       Value sum = 0;
3333       for (NodeIt n(_graph); n != INVALID; ++n) {
3334         if ((*_matching)[n] != INVALID) {
3335           sum += _weight[(*_matching)[n]];
3336         }
3337       }
3338       return sum / 2;
3339     }
3340 
3341     /// \brief Return \c true if the given edge is in the matching.
3342     ///
3343     /// This function returns \c true if the given edge is in the found
3344     /// matching.
3345     ///
3346     /// \pre Either run() or start() must be called before using this function.
matching(const Edge & edge)3347     bool matching(const Edge& edge) const {
3348       return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3349     }
3350 
3351     /// \brief Return the matching arc (or edge) incident to the given node.
3352     ///
3353     /// This function returns the matching arc (or edge) incident to the
3354     /// given node in the found matching or \c INVALID if the node is
3355     /// not covered by the matching.
3356     ///
3357     /// \pre Either run() or start() must be called before using this function.
matching(const Node & node)3358     Arc matching(const Node& node) const {
3359       return (*_matching)[node];
3360     }
3361 
3362     /// \brief Return a const reference to the matching map.
3363     ///
3364     /// This function returns a const reference to a node map that stores
3365     /// the matching arc (or edge) incident to each node.
matchingMap()3366     const MatchingMap& matchingMap() const {
3367       return *_matching;
3368     }
3369 
3370     /// \brief Return the mate of the given node.
3371     ///
3372     /// This function returns the mate of the given node in the found
3373     /// matching or \c INVALID if the node is not covered by the matching.
3374     ///
3375     /// \pre Either run() or start() must be called before using this function.
mate(const Node & node)3376     Node mate(const Node& node) const {
3377       return _graph.target((*_matching)[node]);
3378     }
3379 
3380     /// @}
3381 
3382     /// \name Dual Solution
3383     /// Functions to get the dual solution.\n
3384     /// Either \ref run() or \ref start() function should be called before
3385     /// using them.
3386 
3387     /// @{
3388 
3389     /// \brief Return the value of the dual solution.
3390     ///
3391     /// This function returns the value of the dual solution.
3392     /// It should be equal to the primal value scaled by \ref dualScale
3393     /// "dual scale".
3394     ///
3395     /// \pre Either run() or start() must be called before using this function.
dualValue()3396     Value dualValue() const {
3397       Value sum = 0;
3398       for (NodeIt n(_graph); n != INVALID; ++n) {
3399         sum += nodeValue(n);
3400       }
3401       for (int i = 0; i < blossomNum(); ++i) {
3402         sum += blossomValue(i) * (blossomSize(i) / 2);
3403       }
3404       return sum;
3405     }
3406 
3407     /// \brief Return the dual value (potential) of the given node.
3408     ///
3409     /// This function returns the dual value (potential) of the given node.
3410     ///
3411     /// \pre Either run() or start() must be called before using this function.
nodeValue(const Node & n)3412     Value nodeValue(const Node& n) const {
3413       return (*_node_potential)[n];
3414     }
3415 
3416     /// \brief Return the number of the blossoms in the basis.
3417     ///
3418     /// This function returns the number of the blossoms in the basis.
3419     ///
3420     /// \pre Either run() or start() must be called before using this function.
3421     /// \see BlossomIt
blossomNum()3422     int blossomNum() const {
3423       return _blossom_potential.size();
3424     }
3425 
3426     /// \brief Return the number of the nodes in the given blossom.
3427     ///
3428     /// This function returns the number of the nodes in the given blossom.
3429     ///
3430     /// \pre Either run() or start() must be called before using this function.
3431     /// \see BlossomIt
blossomSize(int k)3432     int blossomSize(int k) const {
3433       return _blossom_potential[k].end - _blossom_potential[k].begin;
3434     }
3435 
3436     /// \brief Return the dual value (ptential) of the given blossom.
3437     ///
3438     /// This function returns the dual value (ptential) of the given blossom.
3439     ///
3440     /// \pre Either run() or start() must be called before using this function.
blossomValue(int k)3441     Value blossomValue(int k) const {
3442       return _blossom_potential[k].value;
3443     }
3444 
3445     /// \brief Iterator for obtaining the nodes of a blossom.
3446     ///
3447     /// This class provides an iterator for obtaining the nodes of the
3448     /// given blossom. It lists a subset of the nodes.
3449     /// Before using this iterator, you must allocate a
3450     /// MaxWeightedPerfectMatching class and execute it.
3451     class BlossomIt {
3452     public:
3453 
3454       /// \brief Constructor.
3455       ///
3456       /// Constructor to get the nodes of the given variable.
3457       ///
3458       /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3459       /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3460       /// must be called before initializing this iterator.
BlossomIt(const MaxWeightedPerfectMatching & algorithm,int variable)3461       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3462         : _algorithm(&algorithm)
3463       {
3464         _index = _algorithm->_blossom_potential[variable].begin;
3465         _last = _algorithm->_blossom_potential[variable].end;
3466       }
3467 
3468       /// \brief Conversion to \c Node.
3469       ///
3470       /// Conversion to \c Node.
Node()3471       operator Node() const {
3472         return _algorithm->_blossom_node_list[_index];
3473       }
3474 
3475       /// \brief Increment operator.
3476       ///
3477       /// Increment operator.
3478       BlossomIt& operator++() {
3479         ++_index;
3480         return *this;
3481       }
3482 
3483       /// \brief Validity checking
3484       ///
3485       /// This function checks whether the iterator is invalid.
3486       bool operator==(Invalid) const { return _index == _last; }
3487 
3488       /// \brief Validity checking
3489       ///
3490       /// This function checks whether the iterator is valid.
3491       bool operator!=(Invalid) const { return _index != _last; }
3492 
3493     private:
3494       const MaxWeightedPerfectMatching* _algorithm;
3495       int _last;
3496       int _index;
3497     };
3498 
3499     /// @}
3500 
3501   };
3502 
3503 } //END OF NAMESPACE LEMON
3504 
3505 #endif //LEMON_MATCHING_H
3506