COIN-OR::LEMON - Graph Library

source: lemon-main/lemon/matching.h @ 1211:a278d16bd2d0

Last change on this file since 1211:a278d16bd2d0 was 1208:c6aa2cc1af04, checked in by Balazs Dezso <deba@…>, 3 years ago

Factor out recursion from weighted matching algorithms (#638)

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