COIN-OR::LEMON - Graph Library

source: lemon/lemon/matching.h @ 954:07ec2b52e53d

Last change on this file since 954:07ec2b52e53d was 954:07ec2b52e53d, checked in by Alpar Juttner <alpar@…>, 9 years ago

Merge #356

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