lemon/matching.h
changeset 784 1a7fe3bef514
parent 594 d657c71db7db
child 867 5b926cc36a4b
child 868 0513ccfea967
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/matching.h	Thu Nov 05 15:50:01 2009 +0100
     1.3 @@ -0,0 +1,3244 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     1.8 + * Copyright (C) 2003-2009
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +#ifndef LEMON_MAX_MATCHING_H
    1.23 +#define LEMON_MAX_MATCHING_H
    1.24 +
    1.25 +#include <vector>
    1.26 +#include <queue>
    1.27 +#include <set>
    1.28 +#include <limits>
    1.29 +
    1.30 +#include <lemon/core.h>
    1.31 +#include <lemon/unionfind.h>
    1.32 +#include <lemon/bin_heap.h>
    1.33 +#include <lemon/maps.h>
    1.34 +
    1.35 +///\ingroup matching
    1.36 +///\file
    1.37 +///\brief Maximum matching algorithms in general graphs.
    1.38 +
    1.39 +namespace lemon {
    1.40 +
    1.41 +  /// \ingroup matching
    1.42 +  ///
    1.43 +  /// \brief Maximum cardinality matching in general graphs
    1.44 +  ///
    1.45 +  /// This class implements Edmonds' alternating forest matching algorithm
    1.46 +  /// for finding a maximum cardinality matching in a general undirected graph.
    1.47 +  /// It can be started from an arbitrary initial matching 
    1.48 +  /// (the default is the empty one).
    1.49 +  ///
    1.50 +  /// The dual solution of the problem is a map of the nodes to
    1.51 +  /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
    1.52 +  /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
    1.53 +  /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
    1.54 +  /// with factor-critical components, the nodes in \c ODD/A form the
    1.55 +  /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
    1.56 +  /// a perfect matching. The number of the factor-critical components
    1.57 +  /// minus the number of barrier nodes is a lower bound on the
    1.58 +  /// unmatched nodes, and the matching is optimal if and only if this bound is
    1.59 +  /// tight. This decomposition can be obtained using \ref status() or
    1.60 +  /// \ref statusMap() after running the algorithm.
    1.61 +  ///
    1.62 +  /// \tparam GR The undirected graph type the algorithm runs on.
    1.63 +  template <typename GR>
    1.64 +  class MaxMatching {
    1.65 +  public:
    1.66 +
    1.67 +    /// The graph type of the algorithm
    1.68 +    typedef GR Graph;
    1.69 +    /// The type of the matching map
    1.70 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
    1.71 +    MatchingMap;
    1.72 +
    1.73 +    ///\brief Status constants for Gallai-Edmonds decomposition.
    1.74 +    ///
    1.75 +    ///These constants are used for indicating the Gallai-Edmonds 
    1.76 +    ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    1.77 +    ///induce a subgraph with factor-critical components, the nodes with
    1.78 +    ///status \c ODD (or \c A) form the canonical barrier, and the nodes
    1.79 +    ///with status \c MATCHED (or \c C) induce a subgraph having a 
    1.80 +    ///perfect matching.
    1.81 +    enum Status {
    1.82 +      EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
    1.83 +      D = 1,
    1.84 +      MATCHED = 0,    ///< = 0. (\c C is an alias for \c MATCHED.)
    1.85 +      C = 0,
    1.86 +      ODD = -1,       ///< = -1. (\c A is an alias for \c ODD.)
    1.87 +      A = -1,
    1.88 +      UNMATCHED = -2  ///< = -2.
    1.89 +    };
    1.90 +
    1.91 +    /// The type of the status map
    1.92 +    typedef typename Graph::template NodeMap<Status> StatusMap;
    1.93 +
    1.94 +  private:
    1.95 +
    1.96 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
    1.97 +
    1.98 +    typedef UnionFindEnum<IntNodeMap> BlossomSet;
    1.99 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
   1.100 +    typedef RangeMap<Node> NodeIntMap;
   1.101 +    typedef MatchingMap EarMap;
   1.102 +    typedef std::vector<Node> NodeQueue;
   1.103 +
   1.104 +    const Graph& _graph;
   1.105 +    MatchingMap* _matching;
   1.106 +    StatusMap* _status;
   1.107 +
   1.108 +    EarMap* _ear;
   1.109 +
   1.110 +    IntNodeMap* _blossom_set_index;
   1.111 +    BlossomSet* _blossom_set;
   1.112 +    NodeIntMap* _blossom_rep;
   1.113 +
   1.114 +    IntNodeMap* _tree_set_index;
   1.115 +    TreeSet* _tree_set;
   1.116 +
   1.117 +    NodeQueue _node_queue;
   1.118 +    int _process, _postpone, _last;
   1.119 +
   1.120 +    int _node_num;
   1.121 +
   1.122 +  private:
   1.123 +
   1.124 +    void createStructures() {
   1.125 +      _node_num = countNodes(_graph);
   1.126 +      if (!_matching) {
   1.127 +        _matching = new MatchingMap(_graph);
   1.128 +      }
   1.129 +      if (!_status) {
   1.130 +        _status = new StatusMap(_graph);
   1.131 +      }
   1.132 +      if (!_ear) {
   1.133 +        _ear = new EarMap(_graph);
   1.134 +      }
   1.135 +      if (!_blossom_set) {
   1.136 +        _blossom_set_index = new IntNodeMap(_graph);
   1.137 +        _blossom_set = new BlossomSet(*_blossom_set_index);
   1.138 +      }
   1.139 +      if (!_blossom_rep) {
   1.140 +        _blossom_rep = new NodeIntMap(_node_num);
   1.141 +      }
   1.142 +      if (!_tree_set) {
   1.143 +        _tree_set_index = new IntNodeMap(_graph);
   1.144 +        _tree_set = new TreeSet(*_tree_set_index);
   1.145 +      }
   1.146 +      _node_queue.resize(_node_num);
   1.147 +    }
   1.148 +
   1.149 +    void destroyStructures() {
   1.150 +      if (_matching) {
   1.151 +        delete _matching;
   1.152 +      }
   1.153 +      if (_status) {
   1.154 +        delete _status;
   1.155 +      }
   1.156 +      if (_ear) {
   1.157 +        delete _ear;
   1.158 +      }
   1.159 +      if (_blossom_set) {
   1.160 +        delete _blossom_set;
   1.161 +        delete _blossom_set_index;
   1.162 +      }
   1.163 +      if (_blossom_rep) {
   1.164 +        delete _blossom_rep;
   1.165 +      }
   1.166 +      if (_tree_set) {
   1.167 +        delete _tree_set_index;
   1.168 +        delete _tree_set;
   1.169 +      }
   1.170 +    }
   1.171 +
   1.172 +    void processDense(const Node& n) {
   1.173 +      _process = _postpone = _last = 0;
   1.174 +      _node_queue[_last++] = n;
   1.175 +
   1.176 +      while (_process != _last) {
   1.177 +        Node u = _node_queue[_process++];
   1.178 +        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   1.179 +          Node v = _graph.target(a);
   1.180 +          if ((*_status)[v] == MATCHED) {
   1.181 +            extendOnArc(a);
   1.182 +          } else if ((*_status)[v] == UNMATCHED) {
   1.183 +            augmentOnArc(a);
   1.184 +            return;
   1.185 +          }
   1.186 +        }
   1.187 +      }
   1.188 +
   1.189 +      while (_postpone != _last) {
   1.190 +        Node u = _node_queue[_postpone++];
   1.191 +
   1.192 +        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
   1.193 +          Node v = _graph.target(a);
   1.194 +
   1.195 +          if ((*_status)[v] == EVEN) {
   1.196 +            if (_blossom_set->find(u) != _blossom_set->find(v)) {
   1.197 +              shrinkOnEdge(a);
   1.198 +            }
   1.199 +          }
   1.200 +
   1.201 +          while (_process != _last) {
   1.202 +            Node w = _node_queue[_process++];
   1.203 +            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
   1.204 +              Node x = _graph.target(b);
   1.205 +              if ((*_status)[x] == MATCHED) {
   1.206 +                extendOnArc(b);
   1.207 +              } else if ((*_status)[x] == UNMATCHED) {
   1.208 +                augmentOnArc(b);
   1.209 +                return;
   1.210 +              }
   1.211 +            }
   1.212 +          }
   1.213 +        }
   1.214 +      }
   1.215 +    }
   1.216 +
   1.217 +    void processSparse(const Node& n) {
   1.218 +      _process = _last = 0;
   1.219 +      _node_queue[_last++] = n;
   1.220 +      while (_process != _last) {
   1.221 +        Node u = _node_queue[_process++];
   1.222 +        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   1.223 +          Node v = _graph.target(a);
   1.224 +
   1.225 +          if ((*_status)[v] == EVEN) {
   1.226 +            if (_blossom_set->find(u) != _blossom_set->find(v)) {
   1.227 +              shrinkOnEdge(a);
   1.228 +            }
   1.229 +          } else if ((*_status)[v] == MATCHED) {
   1.230 +            extendOnArc(a);
   1.231 +          } else if ((*_status)[v] == UNMATCHED) {
   1.232 +            augmentOnArc(a);
   1.233 +            return;
   1.234 +          }
   1.235 +        }
   1.236 +      }
   1.237 +    }
   1.238 +
   1.239 +    void shrinkOnEdge(const Edge& e) {
   1.240 +      Node nca = INVALID;
   1.241 +
   1.242 +      {
   1.243 +        std::set<Node> left_set, right_set;
   1.244 +
   1.245 +        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
   1.246 +        left_set.insert(left);
   1.247 +
   1.248 +        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
   1.249 +        right_set.insert(right);
   1.250 +
   1.251 +        while (true) {
   1.252 +          if ((*_matching)[left] == INVALID) break;
   1.253 +          left = _graph.target((*_matching)[left]);
   1.254 +          left = (*_blossom_rep)[_blossom_set->
   1.255 +                                 find(_graph.target((*_ear)[left]))];
   1.256 +          if (right_set.find(left) != right_set.end()) {
   1.257 +            nca = left;
   1.258 +            break;
   1.259 +          }
   1.260 +          left_set.insert(left);
   1.261 +
   1.262 +          if ((*_matching)[right] == INVALID) break;
   1.263 +          right = _graph.target((*_matching)[right]);
   1.264 +          right = (*_blossom_rep)[_blossom_set->
   1.265 +                                  find(_graph.target((*_ear)[right]))];
   1.266 +          if (left_set.find(right) != left_set.end()) {
   1.267 +            nca = right;
   1.268 +            break;
   1.269 +          }
   1.270 +          right_set.insert(right);
   1.271 +        }
   1.272 +
   1.273 +        if (nca == INVALID) {
   1.274 +          if ((*_matching)[left] == INVALID) {
   1.275 +            nca = right;
   1.276 +            while (left_set.find(nca) == left_set.end()) {
   1.277 +              nca = _graph.target((*_matching)[nca]);
   1.278 +              nca =(*_blossom_rep)[_blossom_set->
   1.279 +                                   find(_graph.target((*_ear)[nca]))];
   1.280 +            }
   1.281 +          } else {
   1.282 +            nca = left;
   1.283 +            while (right_set.find(nca) == right_set.end()) {
   1.284 +              nca = _graph.target((*_matching)[nca]);
   1.285 +              nca = (*_blossom_rep)[_blossom_set->
   1.286 +                                   find(_graph.target((*_ear)[nca]))];
   1.287 +            }
   1.288 +          }
   1.289 +        }
   1.290 +      }
   1.291 +
   1.292 +      {
   1.293 +
   1.294 +        Node node = _graph.u(e);
   1.295 +        Arc arc = _graph.direct(e, true);
   1.296 +        Node base = (*_blossom_rep)[_blossom_set->find(node)];
   1.297 +
   1.298 +        while (base != nca) {
   1.299 +          (*_ear)[node] = arc;
   1.300 +
   1.301 +          Node n = node;
   1.302 +          while (n != base) {
   1.303 +            n = _graph.target((*_matching)[n]);
   1.304 +            Arc a = (*_ear)[n];
   1.305 +            n = _graph.target(a);
   1.306 +            (*_ear)[n] = _graph.oppositeArc(a);
   1.307 +          }
   1.308 +          node = _graph.target((*_matching)[base]);
   1.309 +          _tree_set->erase(base);
   1.310 +          _tree_set->erase(node);
   1.311 +          _blossom_set->insert(node, _blossom_set->find(base));
   1.312 +          (*_status)[node] = EVEN;
   1.313 +          _node_queue[_last++] = node;
   1.314 +          arc = _graph.oppositeArc((*_ear)[node]);
   1.315 +          node = _graph.target((*_ear)[node]);
   1.316 +          base = (*_blossom_rep)[_blossom_set->find(node)];
   1.317 +          _blossom_set->join(_graph.target(arc), base);
   1.318 +        }
   1.319 +      }
   1.320 +
   1.321 +      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
   1.322 +
   1.323 +      {
   1.324 +
   1.325 +        Node node = _graph.v(e);
   1.326 +        Arc arc = _graph.direct(e, false);
   1.327 +        Node base = (*_blossom_rep)[_blossom_set->find(node)];
   1.328 +
   1.329 +        while (base != nca) {
   1.330 +          (*_ear)[node] = arc;
   1.331 +
   1.332 +          Node n = node;
   1.333 +          while (n != base) {
   1.334 +            n = _graph.target((*_matching)[n]);
   1.335 +            Arc a = (*_ear)[n];
   1.336 +            n = _graph.target(a);
   1.337 +            (*_ear)[n] = _graph.oppositeArc(a);
   1.338 +          }
   1.339 +          node = _graph.target((*_matching)[base]);
   1.340 +          _tree_set->erase(base);
   1.341 +          _tree_set->erase(node);
   1.342 +          _blossom_set->insert(node, _blossom_set->find(base));
   1.343 +          (*_status)[node] = EVEN;
   1.344 +          _node_queue[_last++] = node;
   1.345 +          arc = _graph.oppositeArc((*_ear)[node]);
   1.346 +          node = _graph.target((*_ear)[node]);
   1.347 +          base = (*_blossom_rep)[_blossom_set->find(node)];
   1.348 +          _blossom_set->join(_graph.target(arc), base);
   1.349 +        }
   1.350 +      }
   1.351 +
   1.352 +      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
   1.353 +    }
   1.354 +
   1.355 +    void extendOnArc(const Arc& a) {
   1.356 +      Node base = _graph.source(a);
   1.357 +      Node odd = _graph.target(a);
   1.358 +
   1.359 +      (*_ear)[odd] = _graph.oppositeArc(a);
   1.360 +      Node even = _graph.target((*_matching)[odd]);
   1.361 +      (*_blossom_rep)[_blossom_set->insert(even)] = even;
   1.362 +      (*_status)[odd] = ODD;
   1.363 +      (*_status)[even] = EVEN;
   1.364 +      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
   1.365 +      _tree_set->insert(odd, tree);
   1.366 +      _tree_set->insert(even, tree);
   1.367 +      _node_queue[_last++] = even;
   1.368 +
   1.369 +    }
   1.370 +
   1.371 +    void augmentOnArc(const Arc& a) {
   1.372 +      Node even = _graph.source(a);
   1.373 +      Node odd = _graph.target(a);
   1.374 +
   1.375 +      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
   1.376 +
   1.377 +      (*_matching)[odd] = _graph.oppositeArc(a);
   1.378 +      (*_status)[odd] = MATCHED;
   1.379 +
   1.380 +      Arc arc = (*_matching)[even];
   1.381 +      (*_matching)[even] = a;
   1.382 +
   1.383 +      while (arc != INVALID) {
   1.384 +        odd = _graph.target(arc);
   1.385 +        arc = (*_ear)[odd];
   1.386 +        even = _graph.target(arc);
   1.387 +        (*_matching)[odd] = arc;
   1.388 +        arc = (*_matching)[even];
   1.389 +        (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
   1.390 +      }
   1.391 +
   1.392 +      for (typename TreeSet::ItemIt it(*_tree_set, tree);
   1.393 +           it != INVALID; ++it) {
   1.394 +        if ((*_status)[it] == ODD) {
   1.395 +          (*_status)[it] = MATCHED;
   1.396 +        } else {
   1.397 +          int blossom = _blossom_set->find(it);
   1.398 +          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
   1.399 +               jt != INVALID; ++jt) {
   1.400 +            (*_status)[jt] = MATCHED;
   1.401 +          }
   1.402 +          _blossom_set->eraseClass(blossom);
   1.403 +        }
   1.404 +      }
   1.405 +      _tree_set->eraseClass(tree);
   1.406 +
   1.407 +    }
   1.408 +
   1.409 +  public:
   1.410 +
   1.411 +    /// \brief Constructor
   1.412 +    ///
   1.413 +    /// Constructor.
   1.414 +    MaxMatching(const Graph& graph)
   1.415 +      : _graph(graph), _matching(0), _status(0), _ear(0),
   1.416 +        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
   1.417 +        _tree_set_index(0), _tree_set(0) {}
   1.418 +
   1.419 +    ~MaxMatching() {
   1.420 +      destroyStructures();
   1.421 +    }
   1.422 +
   1.423 +    /// \name Execution Control
   1.424 +    /// The simplest way to execute the algorithm is to use the
   1.425 +    /// \c run() member function.\n
   1.426 +    /// If you need better control on the execution, you have to call
   1.427 +    /// one of the functions \ref init(), \ref greedyInit() or
   1.428 +    /// \ref matchingInit() first, then you can start the algorithm with
   1.429 +    /// \ref startSparse() or \ref startDense().
   1.430 +
   1.431 +    ///@{
   1.432 +
   1.433 +    /// \brief Set the initial matching to the empty matching.
   1.434 +    ///
   1.435 +    /// This function sets the initial matching to the empty matching.
   1.436 +    void init() {
   1.437 +      createStructures();
   1.438 +      for(NodeIt n(_graph); n != INVALID; ++n) {
   1.439 +        (*_matching)[n] = INVALID;
   1.440 +        (*_status)[n] = UNMATCHED;
   1.441 +      }
   1.442 +    }
   1.443 +
   1.444 +    /// \brief Find an initial matching in a greedy way.
   1.445 +    ///
   1.446 +    /// This function finds an initial matching in a greedy way.
   1.447 +    void greedyInit() {
   1.448 +      createStructures();
   1.449 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.450 +        (*_matching)[n] = INVALID;
   1.451 +        (*_status)[n] = UNMATCHED;
   1.452 +      }
   1.453 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.454 +        if ((*_matching)[n] == INVALID) {
   1.455 +          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
   1.456 +            Node v = _graph.target(a);
   1.457 +            if ((*_matching)[v] == INVALID && v != n) {
   1.458 +              (*_matching)[n] = a;
   1.459 +              (*_status)[n] = MATCHED;
   1.460 +              (*_matching)[v] = _graph.oppositeArc(a);
   1.461 +              (*_status)[v] = MATCHED;
   1.462 +              break;
   1.463 +            }
   1.464 +          }
   1.465 +        }
   1.466 +      }
   1.467 +    }
   1.468 +
   1.469 +
   1.470 +    /// \brief Initialize the matching from a map.
   1.471 +    ///
   1.472 +    /// This function initializes the matching from a \c bool valued edge
   1.473 +    /// map. This map should have the property that there are no two incident
   1.474 +    /// edges with \c true value, i.e. it really contains a matching.
   1.475 +    /// \return \c true if the map contains a matching.
   1.476 +    template <typename MatchingMap>
   1.477 +    bool matchingInit(const MatchingMap& matching) {
   1.478 +      createStructures();
   1.479 +
   1.480 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.481 +        (*_matching)[n] = INVALID;
   1.482 +        (*_status)[n] = UNMATCHED;
   1.483 +      }
   1.484 +      for(EdgeIt e(_graph); e!=INVALID; ++e) {
   1.485 +        if (matching[e]) {
   1.486 +
   1.487 +          Node u = _graph.u(e);
   1.488 +          if ((*_matching)[u] != INVALID) return false;
   1.489 +          (*_matching)[u] = _graph.direct(e, true);
   1.490 +          (*_status)[u] = MATCHED;
   1.491 +
   1.492 +          Node v = _graph.v(e);
   1.493 +          if ((*_matching)[v] != INVALID) return false;
   1.494 +          (*_matching)[v] = _graph.direct(e, false);
   1.495 +          (*_status)[v] = MATCHED;
   1.496 +        }
   1.497 +      }
   1.498 +      return true;
   1.499 +    }
   1.500 +
   1.501 +    /// \brief Start Edmonds' algorithm
   1.502 +    ///
   1.503 +    /// This function runs the original Edmonds' algorithm.
   1.504 +    ///
   1.505 +    /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
   1.506 +    /// called before using this function.
   1.507 +    void startSparse() {
   1.508 +      for(NodeIt n(_graph); n != INVALID; ++n) {
   1.509 +        if ((*_status)[n] == UNMATCHED) {
   1.510 +          (*_blossom_rep)[_blossom_set->insert(n)] = n;
   1.511 +          _tree_set->insert(n);
   1.512 +          (*_status)[n] = EVEN;
   1.513 +          processSparse(n);
   1.514 +        }
   1.515 +      }
   1.516 +    }
   1.517 +
   1.518 +    /// \brief Start Edmonds' algorithm with a heuristic improvement 
   1.519 +    /// for dense graphs
   1.520 +    ///
   1.521 +    /// This function runs Edmonds' algorithm with a heuristic of postponing
   1.522 +    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
   1.523 +    ///
   1.524 +    /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
   1.525 +    /// called before using this function.
   1.526 +    void startDense() {
   1.527 +      for(NodeIt n(_graph); n != INVALID; ++n) {
   1.528 +        if ((*_status)[n] == UNMATCHED) {
   1.529 +          (*_blossom_rep)[_blossom_set->insert(n)] = n;
   1.530 +          _tree_set->insert(n);
   1.531 +          (*_status)[n] = EVEN;
   1.532 +          processDense(n);
   1.533 +        }
   1.534 +      }
   1.535 +    }
   1.536 +
   1.537 +
   1.538 +    /// \brief Run Edmonds' algorithm
   1.539 +    ///
   1.540 +    /// This function runs Edmonds' algorithm. An additional heuristic of 
   1.541 +    /// postponing shrinks is used for relatively dense graphs 
   1.542 +    /// (for which <tt>m>=2*n</tt> holds).
   1.543 +    void run() {
   1.544 +      if (countEdges(_graph) < 2 * countNodes(_graph)) {
   1.545 +        greedyInit();
   1.546 +        startSparse();
   1.547 +      } else {
   1.548 +        init();
   1.549 +        startDense();
   1.550 +      }
   1.551 +    }
   1.552 +
   1.553 +    /// @}
   1.554 +
   1.555 +    /// \name Primal Solution
   1.556 +    /// Functions to get the primal solution, i.e. the maximum matching.
   1.557 +
   1.558 +    /// @{
   1.559 +
   1.560 +    /// \brief Return the size (cardinality) of the matching.
   1.561 +    ///
   1.562 +    /// This function returns the size (cardinality) of the current matching. 
   1.563 +    /// After run() it returns the size of the maximum matching in the graph.
   1.564 +    int matchingSize() const {
   1.565 +      int size = 0;
   1.566 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.567 +        if ((*_matching)[n] != INVALID) {
   1.568 +          ++size;
   1.569 +        }
   1.570 +      }
   1.571 +      return size / 2;
   1.572 +    }
   1.573 +
   1.574 +    /// \brief Return \c true if the given edge is in the matching.
   1.575 +    ///
   1.576 +    /// This function returns \c true if the given edge is in the current 
   1.577 +    /// matching.
   1.578 +    bool matching(const Edge& edge) const {
   1.579 +      return edge == (*_matching)[_graph.u(edge)];
   1.580 +    }
   1.581 +
   1.582 +    /// \brief Return the matching arc (or edge) incident to the given node.
   1.583 +    ///
   1.584 +    /// This function returns the matching arc (or edge) incident to the
   1.585 +    /// given node in the current matching or \c INVALID if the node is 
   1.586 +    /// not covered by the matching.
   1.587 +    Arc matching(const Node& n) const {
   1.588 +      return (*_matching)[n];
   1.589 +    }
   1.590 +
   1.591 +    /// \brief Return a const reference to the matching map.
   1.592 +    ///
   1.593 +    /// This function returns a const reference to a node map that stores
   1.594 +    /// the matching arc (or edge) incident to each node.
   1.595 +    const MatchingMap& matchingMap() const {
   1.596 +      return *_matching;
   1.597 +    }
   1.598 +
   1.599 +    /// \brief Return the mate of the given node.
   1.600 +    ///
   1.601 +    /// This function returns the mate of the given node in the current 
   1.602 +    /// matching or \c INVALID if the node is not covered by the matching.
   1.603 +    Node mate(const Node& n) const {
   1.604 +      return (*_matching)[n] != INVALID ?
   1.605 +        _graph.target((*_matching)[n]) : INVALID;
   1.606 +    }
   1.607 +
   1.608 +    /// @}
   1.609 +
   1.610 +    /// \name Dual Solution
   1.611 +    /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
   1.612 +    /// decomposition.
   1.613 +
   1.614 +    /// @{
   1.615 +
   1.616 +    /// \brief Return the status of the given node in the Edmonds-Gallai
   1.617 +    /// decomposition.
   1.618 +    ///
   1.619 +    /// This function returns the \ref Status "status" of the given node
   1.620 +    /// in the Edmonds-Gallai decomposition.
   1.621 +    Status status(const Node& n) const {
   1.622 +      return (*_status)[n];
   1.623 +    }
   1.624 +
   1.625 +    /// \brief Return a const reference to the status map, which stores
   1.626 +    /// the Edmonds-Gallai decomposition.
   1.627 +    ///
   1.628 +    /// This function returns a const reference to a node map that stores the
   1.629 +    /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
   1.630 +    const StatusMap& statusMap() const {
   1.631 +      return *_status;
   1.632 +    }
   1.633 +
   1.634 +    /// \brief Return \c true if the given node is in the barrier.
   1.635 +    ///
   1.636 +    /// This function returns \c true if the given node is in the barrier.
   1.637 +    bool barrier(const Node& n) const {
   1.638 +      return (*_status)[n] == ODD;
   1.639 +    }
   1.640 +
   1.641 +    /// @}
   1.642 +
   1.643 +  };
   1.644 +
   1.645 +  /// \ingroup matching
   1.646 +  ///
   1.647 +  /// \brief Weighted matching in general graphs
   1.648 +  ///
   1.649 +  /// This class provides an efficient implementation of Edmond's
   1.650 +  /// maximum weighted matching algorithm. The implementation is based
   1.651 +  /// on extensive use of priority queues and provides
   1.652 +  /// \f$O(nm\log n)\f$ time complexity.
   1.653 +  ///
   1.654 +  /// The maximum weighted matching problem is to find a subset of the 
   1.655 +  /// edges in an undirected graph with maximum overall weight for which 
   1.656 +  /// each node has at most one incident edge.
   1.657 +  /// It can be formulated with the following linear program.
   1.658 +  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   1.659 +  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
   1.660 +      \quad \forall B\in\mathcal{O}\f] */
   1.661 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
   1.662 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
   1.663 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   1.664 +  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
   1.665 +  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
   1.666 +  /// subsets of the nodes.
   1.667 +  ///
   1.668 +  /// The algorithm calculates an optimal matching and a proof of the
   1.669 +  /// optimality. The solution of the dual problem can be used to check
   1.670 +  /// the result of the algorithm. The dual linear problem is the
   1.671 +  /// following.
   1.672 +  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
   1.673 +      z_B \ge w_{uv} \quad \forall uv\in E\f] */
   1.674 +  /// \f[y_u \ge 0 \quad \forall u \in V\f]
   1.675 +  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
   1.676 +  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
   1.677 +      \frac{\vert B \vert - 1}{2}z_B\f] */
   1.678 +  ///
   1.679 +  /// The algorithm can be executed with the run() function. 
   1.680 +  /// After it the matching (the primal solution) and the dual solution
   1.681 +  /// can be obtained using the query functions and the 
   1.682 +  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
   1.683 +  /// which is able to iterate on the nodes of a blossom. 
   1.684 +  /// If the value type is integer, then the dual solution is multiplied
   1.685 +  /// by \ref MaxWeightedMatching::dualScale "4".
   1.686 +  ///
   1.687 +  /// \tparam GR The undirected graph type the algorithm runs on.
   1.688 +  /// \tparam WM The type edge weight map. The default type is 
   1.689 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   1.690 +#ifdef DOXYGEN
   1.691 +  template <typename GR, typename WM>
   1.692 +#else
   1.693 +  template <typename GR,
   1.694 +            typename WM = typename GR::template EdgeMap<int> >
   1.695 +#endif
   1.696 +  class MaxWeightedMatching {
   1.697 +  public:
   1.698 +
   1.699 +    /// The graph type of the algorithm
   1.700 +    typedef GR Graph;
   1.701 +    /// The type of the edge weight map
   1.702 +    typedef WM WeightMap;
   1.703 +    /// The value type of the edge weights
   1.704 +    typedef typename WeightMap::Value Value;
   1.705 +
   1.706 +    /// The type of the matching map
   1.707 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
   1.708 +    MatchingMap;
   1.709 +
   1.710 +    /// \brief Scaling factor for dual solution
   1.711 +    ///
   1.712 +    /// Scaling factor for dual solution. It is equal to 4 or 1
   1.713 +    /// according to the value type.
   1.714 +    static const int dualScale =
   1.715 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
   1.716 +
   1.717 +  private:
   1.718 +
   1.719 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
   1.720 +
   1.721 +    typedef typename Graph::template NodeMap<Value> NodePotential;
   1.722 +    typedef std::vector<Node> BlossomNodeList;
   1.723 +
   1.724 +    struct BlossomVariable {
   1.725 +      int begin, end;
   1.726 +      Value value;
   1.727 +
   1.728 +      BlossomVariable(int _begin, int _end, Value _value)
   1.729 +        : begin(_begin), end(_end), value(_value) {}
   1.730 +
   1.731 +    };
   1.732 +
   1.733 +    typedef std::vector<BlossomVariable> BlossomPotential;
   1.734 +
   1.735 +    const Graph& _graph;
   1.736 +    const WeightMap& _weight;
   1.737 +
   1.738 +    MatchingMap* _matching;
   1.739 +
   1.740 +    NodePotential* _node_potential;
   1.741 +
   1.742 +    BlossomPotential _blossom_potential;
   1.743 +    BlossomNodeList _blossom_node_list;
   1.744 +
   1.745 +    int _node_num;
   1.746 +    int _blossom_num;
   1.747 +
   1.748 +    typedef RangeMap<int> IntIntMap;
   1.749 +
   1.750 +    enum Status {
   1.751 +      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
   1.752 +    };
   1.753 +
   1.754 +    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
   1.755 +    struct BlossomData {
   1.756 +      int tree;
   1.757 +      Status status;
   1.758 +      Arc pred, next;
   1.759 +      Value pot, offset;
   1.760 +      Node base;
   1.761 +    };
   1.762 +
   1.763 +    IntNodeMap *_blossom_index;
   1.764 +    BlossomSet *_blossom_set;
   1.765 +    RangeMap<BlossomData>* _blossom_data;
   1.766 +
   1.767 +    IntNodeMap *_node_index;
   1.768 +    IntArcMap *_node_heap_index;
   1.769 +
   1.770 +    struct NodeData {
   1.771 +
   1.772 +      NodeData(IntArcMap& node_heap_index)
   1.773 +        : heap(node_heap_index) {}
   1.774 +
   1.775 +      int blossom;
   1.776 +      Value pot;
   1.777 +      BinHeap<Value, IntArcMap> heap;
   1.778 +      std::map<int, Arc> heap_index;
   1.779 +
   1.780 +      int tree;
   1.781 +    };
   1.782 +
   1.783 +    RangeMap<NodeData>* _node_data;
   1.784 +
   1.785 +    typedef ExtendFindEnum<IntIntMap> TreeSet;
   1.786 +
   1.787 +    IntIntMap *_tree_set_index;
   1.788 +    TreeSet *_tree_set;
   1.789 +
   1.790 +    IntNodeMap *_delta1_index;
   1.791 +    BinHeap<Value, IntNodeMap> *_delta1;
   1.792 +
   1.793 +    IntIntMap *_delta2_index;
   1.794 +    BinHeap<Value, IntIntMap> *_delta2;
   1.795 +
   1.796 +    IntEdgeMap *_delta3_index;
   1.797 +    BinHeap<Value, IntEdgeMap> *_delta3;
   1.798 +
   1.799 +    IntIntMap *_delta4_index;
   1.800 +    BinHeap<Value, IntIntMap> *_delta4;
   1.801 +
   1.802 +    Value _delta_sum;
   1.803 +
   1.804 +    void createStructures() {
   1.805 +      _node_num = countNodes(_graph);
   1.806 +      _blossom_num = _node_num * 3 / 2;
   1.807 +
   1.808 +      if (!_matching) {
   1.809 +        _matching = new MatchingMap(_graph);
   1.810 +      }
   1.811 +      if (!_node_potential) {
   1.812 +        _node_potential = new NodePotential(_graph);
   1.813 +      }
   1.814 +      if (!_blossom_set) {
   1.815 +        _blossom_index = new IntNodeMap(_graph);
   1.816 +        _blossom_set = new BlossomSet(*_blossom_index);
   1.817 +        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
   1.818 +      }
   1.819 +
   1.820 +      if (!_node_index) {
   1.821 +        _node_index = new IntNodeMap(_graph);
   1.822 +        _node_heap_index = new IntArcMap(_graph);
   1.823 +        _node_data = new RangeMap<NodeData>(_node_num,
   1.824 +                                              NodeData(*_node_heap_index));
   1.825 +      }
   1.826 +
   1.827 +      if (!_tree_set) {
   1.828 +        _tree_set_index = new IntIntMap(_blossom_num);
   1.829 +        _tree_set = new TreeSet(*_tree_set_index);
   1.830 +      }
   1.831 +      if (!_delta1) {
   1.832 +        _delta1_index = new IntNodeMap(_graph);
   1.833 +        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   1.834 +      }
   1.835 +      if (!_delta2) {
   1.836 +        _delta2_index = new IntIntMap(_blossom_num);
   1.837 +        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
   1.838 +      }
   1.839 +      if (!_delta3) {
   1.840 +        _delta3_index = new IntEdgeMap(_graph);
   1.841 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   1.842 +      }
   1.843 +      if (!_delta4) {
   1.844 +        _delta4_index = new IntIntMap(_blossom_num);
   1.845 +        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
   1.846 +      }
   1.847 +    }
   1.848 +
   1.849 +    void destroyStructures() {
   1.850 +      _node_num = countNodes(_graph);
   1.851 +      _blossom_num = _node_num * 3 / 2;
   1.852 +
   1.853 +      if (_matching) {
   1.854 +        delete _matching;
   1.855 +      }
   1.856 +      if (_node_potential) {
   1.857 +        delete _node_potential;
   1.858 +      }
   1.859 +      if (_blossom_set) {
   1.860 +        delete _blossom_index;
   1.861 +        delete _blossom_set;
   1.862 +        delete _blossom_data;
   1.863 +      }
   1.864 +
   1.865 +      if (_node_index) {
   1.866 +        delete _node_index;
   1.867 +        delete _node_heap_index;
   1.868 +        delete _node_data;
   1.869 +      }
   1.870 +
   1.871 +      if (_tree_set) {
   1.872 +        delete _tree_set_index;
   1.873 +        delete _tree_set;
   1.874 +      }
   1.875 +      if (_delta1) {
   1.876 +        delete _delta1_index;
   1.877 +        delete _delta1;
   1.878 +      }
   1.879 +      if (_delta2) {
   1.880 +        delete _delta2_index;
   1.881 +        delete _delta2;
   1.882 +      }
   1.883 +      if (_delta3) {
   1.884 +        delete _delta3_index;
   1.885 +        delete _delta3;
   1.886 +      }
   1.887 +      if (_delta4) {
   1.888 +        delete _delta4_index;
   1.889 +        delete _delta4;
   1.890 +      }
   1.891 +    }
   1.892 +
   1.893 +    void matchedToEven(int blossom, int tree) {
   1.894 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   1.895 +        _delta2->erase(blossom);
   1.896 +      }
   1.897 +
   1.898 +      if (!_blossom_set->trivial(blossom)) {
   1.899 +        (*_blossom_data)[blossom].pot -=
   1.900 +          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
   1.901 +      }
   1.902 +
   1.903 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.904 +           n != INVALID; ++n) {
   1.905 +
   1.906 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.907 +        int ni = (*_node_index)[n];
   1.908 +
   1.909 +        (*_node_data)[ni].heap.clear();
   1.910 +        (*_node_data)[ni].heap_index.clear();
   1.911 +
   1.912 +        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
   1.913 +
   1.914 +        _delta1->push(n, (*_node_data)[ni].pot);
   1.915 +
   1.916 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.917 +          Node v = _graph.source(e);
   1.918 +          int vb = _blossom_set->find(v);
   1.919 +          int vi = (*_node_index)[v];
   1.920 +
   1.921 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.922 +            dualScale * _weight[e];
   1.923 +
   1.924 +          if ((*_blossom_data)[vb].status == EVEN) {
   1.925 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.926 +              _delta3->push(e, rw / 2);
   1.927 +            }
   1.928 +          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.929 +            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.930 +              _delta3->push(e, rw);
   1.931 +            }
   1.932 +          } else {
   1.933 +            typename std::map<int, Arc>::iterator it =
   1.934 +              (*_node_data)[vi].heap_index.find(tree);
   1.935 +
   1.936 +            if (it != (*_node_data)[vi].heap_index.end()) {
   1.937 +              if ((*_node_data)[vi].heap[it->second] > rw) {
   1.938 +                (*_node_data)[vi].heap.replace(it->second, e);
   1.939 +                (*_node_data)[vi].heap.decrease(e, rw);
   1.940 +                it->second = e;
   1.941 +              }
   1.942 +            } else {
   1.943 +              (*_node_data)[vi].heap.push(e, rw);
   1.944 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   1.945 +            }
   1.946 +
   1.947 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   1.948 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   1.949 +
   1.950 +              if ((*_blossom_data)[vb].status == MATCHED) {
   1.951 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   1.952 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.953 +                               (*_blossom_data)[vb].offset);
   1.954 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.955 +                           (*_blossom_data)[vb].offset){
   1.956 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.957 +                                   (*_blossom_data)[vb].offset);
   1.958 +                }
   1.959 +              }
   1.960 +            }
   1.961 +          }
   1.962 +        }
   1.963 +      }
   1.964 +      (*_blossom_data)[blossom].offset = 0;
   1.965 +    }
   1.966 +
   1.967 +    void matchedToOdd(int blossom) {
   1.968 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   1.969 +        _delta2->erase(blossom);
   1.970 +      }
   1.971 +      (*_blossom_data)[blossom].offset += _delta_sum;
   1.972 +      if (!_blossom_set->trivial(blossom)) {
   1.973 +        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   1.974 +                     (*_blossom_data)[blossom].offset);
   1.975 +      }
   1.976 +    }
   1.977 +
   1.978 +    void evenToMatched(int blossom, int tree) {
   1.979 +      if (!_blossom_set->trivial(blossom)) {
   1.980 +        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   1.981 +      }
   1.982 +
   1.983 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.984 +           n != INVALID; ++n) {
   1.985 +        int ni = (*_node_index)[n];
   1.986 +        (*_node_data)[ni].pot -= _delta_sum;
   1.987 +
   1.988 +        _delta1->erase(n);
   1.989 +
   1.990 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.991 +          Node v = _graph.source(e);
   1.992 +          int vb = _blossom_set->find(v);
   1.993 +          int vi = (*_node_index)[v];
   1.994 +
   1.995 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.996 +            dualScale * _weight[e];
   1.997 +
   1.998 +          if (vb == blossom) {
   1.999 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.1000 +              _delta3->erase(e);
  1.1001 +            }
  1.1002 +          } else if ((*_blossom_data)[vb].status == EVEN) {
  1.1003 +
  1.1004 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.1005 +              _delta3->erase(e);
  1.1006 +            }
  1.1007 +
  1.1008 +            int vt = _tree_set->find(vb);
  1.1009 +
  1.1010 +            if (vt != tree) {
  1.1011 +
  1.1012 +              Arc r = _graph.oppositeArc(e);
  1.1013 +
  1.1014 +              typename std::map<int, Arc>::iterator it =
  1.1015 +                (*_node_data)[ni].heap_index.find(vt);
  1.1016 +
  1.1017 +              if (it != (*_node_data)[ni].heap_index.end()) {
  1.1018 +                if ((*_node_data)[ni].heap[it->second] > rw) {
  1.1019 +                  (*_node_data)[ni].heap.replace(it->second, r);
  1.1020 +                  (*_node_data)[ni].heap.decrease(r, rw);
  1.1021 +                  it->second = r;
  1.1022 +                }
  1.1023 +              } else {
  1.1024 +                (*_node_data)[ni].heap.push(r, rw);
  1.1025 +                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  1.1026 +              }
  1.1027 +
  1.1028 +              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  1.1029 +                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1.1030 +
  1.1031 +                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  1.1032 +                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1.1033 +                               (*_blossom_data)[blossom].offset);
  1.1034 +                } else if ((*_delta2)[blossom] >
  1.1035 +                           _blossom_set->classPrio(blossom) -
  1.1036 +                           (*_blossom_data)[blossom].offset){
  1.1037 +                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  1.1038 +                                   (*_blossom_data)[blossom].offset);
  1.1039 +                }
  1.1040 +              }
  1.1041 +            }
  1.1042 +
  1.1043 +          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1.1044 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.1045 +              _delta3->erase(e);
  1.1046 +            }
  1.1047 +          } else {
  1.1048 +
  1.1049 +            typename std::map<int, Arc>::iterator it =
  1.1050 +              (*_node_data)[vi].heap_index.find(tree);
  1.1051 +
  1.1052 +            if (it != (*_node_data)[vi].heap_index.end()) {
  1.1053 +              (*_node_data)[vi].heap.erase(it->second);
  1.1054 +              (*_node_data)[vi].heap_index.erase(it);
  1.1055 +              if ((*_node_data)[vi].heap.empty()) {
  1.1056 +                _blossom_set->increase(v, std::numeric_limits<Value>::max());
  1.1057 +              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
  1.1058 +                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
  1.1059 +              }
  1.1060 +
  1.1061 +              if ((*_blossom_data)[vb].status == MATCHED) {
  1.1062 +                if (_blossom_set->classPrio(vb) ==
  1.1063 +                    std::numeric_limits<Value>::max()) {
  1.1064 +                  _delta2->erase(vb);
  1.1065 +                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
  1.1066 +                           (*_blossom_data)[vb].offset) {
  1.1067 +                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
  1.1068 +                                   (*_blossom_data)[vb].offset);
  1.1069 +                }
  1.1070 +              }
  1.1071 +            }
  1.1072 +          }
  1.1073 +        }
  1.1074 +      }
  1.1075 +    }
  1.1076 +
  1.1077 +    void oddToMatched(int blossom) {
  1.1078 +      (*_blossom_data)[blossom].offset -= _delta_sum;
  1.1079 +
  1.1080 +      if (_blossom_set->classPrio(blossom) !=
  1.1081 +          std::numeric_limits<Value>::max()) {
  1.1082 +        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1.1083 +                       (*_blossom_data)[blossom].offset);
  1.1084 +      }
  1.1085 +
  1.1086 +      if (!_blossom_set->trivial(blossom)) {
  1.1087 +        _delta4->erase(blossom);
  1.1088 +      }
  1.1089 +    }
  1.1090 +
  1.1091 +    void oddToEven(int blossom, int tree) {
  1.1092 +      if (!_blossom_set->trivial(blossom)) {
  1.1093 +        _delta4->erase(blossom);
  1.1094 +        (*_blossom_data)[blossom].pot -=
  1.1095 +          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
  1.1096 +      }
  1.1097 +
  1.1098 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1.1099 +           n != INVALID; ++n) {
  1.1100 +        int ni = (*_node_index)[n];
  1.1101 +
  1.1102 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1.1103 +
  1.1104 +        (*_node_data)[ni].heap.clear();
  1.1105 +        (*_node_data)[ni].heap_index.clear();
  1.1106 +        (*_node_data)[ni].pot +=
  1.1107 +          2 * _delta_sum - (*_blossom_data)[blossom].offset;
  1.1108 +
  1.1109 +        _delta1->push(n, (*_node_data)[ni].pot);
  1.1110 +
  1.1111 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1.1112 +          Node v = _graph.source(e);
  1.1113 +          int vb = _blossom_set->find(v);
  1.1114 +          int vi = (*_node_index)[v];
  1.1115 +
  1.1116 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.1117 +            dualScale * _weight[e];
  1.1118 +
  1.1119 +          if ((*_blossom_data)[vb].status == EVEN) {
  1.1120 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  1.1121 +              _delta3->push(e, rw / 2);
  1.1122 +            }
  1.1123 +          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1.1124 +            if (_delta3->state(e) != _delta3->IN_HEAP) {
  1.1125 +              _delta3->push(e, rw);
  1.1126 +            }
  1.1127 +          } else {
  1.1128 +
  1.1129 +            typename std::map<int, Arc>::iterator it =
  1.1130 +              (*_node_data)[vi].heap_index.find(tree);
  1.1131 +
  1.1132 +            if (it != (*_node_data)[vi].heap_index.end()) {
  1.1133 +              if ((*_node_data)[vi].heap[it->second] > rw) {
  1.1134 +                (*_node_data)[vi].heap.replace(it->second, e);
  1.1135 +                (*_node_data)[vi].heap.decrease(e, rw);
  1.1136 +                it->second = e;
  1.1137 +              }
  1.1138 +            } else {
  1.1139 +              (*_node_data)[vi].heap.push(e, rw);
  1.1140 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  1.1141 +            }
  1.1142 +
  1.1143 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  1.1144 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  1.1145 +
  1.1146 +              if ((*_blossom_data)[vb].status == MATCHED) {
  1.1147 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
  1.1148 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
  1.1149 +                               (*_blossom_data)[vb].offset);
  1.1150 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  1.1151 +                           (*_blossom_data)[vb].offset) {
  1.1152 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  1.1153 +                                   (*_blossom_data)[vb].offset);
  1.1154 +                }
  1.1155 +              }
  1.1156 +            }
  1.1157 +          }
  1.1158 +        }
  1.1159 +      }
  1.1160 +      (*_blossom_data)[blossom].offset = 0;
  1.1161 +    }
  1.1162 +
  1.1163 +
  1.1164 +    void matchedToUnmatched(int blossom) {
  1.1165 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1.1166 +        _delta2->erase(blossom);
  1.1167 +      }
  1.1168 +
  1.1169 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1.1170 +           n != INVALID; ++n) {
  1.1171 +        int ni = (*_node_index)[n];
  1.1172 +
  1.1173 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1.1174 +
  1.1175 +        (*_node_data)[ni].heap.clear();
  1.1176 +        (*_node_data)[ni].heap_index.clear();
  1.1177 +
  1.1178 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1.1179 +          Node v = _graph.target(e);
  1.1180 +          int vb = _blossom_set->find(v);
  1.1181 +          int vi = (*_node_index)[v];
  1.1182 +
  1.1183 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.1184 +            dualScale * _weight[e];
  1.1185 +
  1.1186 +          if ((*_blossom_data)[vb].status == EVEN) {
  1.1187 +            if (_delta3->state(e) != _delta3->IN_HEAP) {
  1.1188 +              _delta3->push(e, rw);
  1.1189 +            }
  1.1190 +          }
  1.1191 +        }
  1.1192 +      }
  1.1193 +    }
  1.1194 +
  1.1195 +    void unmatchedToMatched(int blossom) {
  1.1196 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1.1197 +           n != INVALID; ++n) {
  1.1198 +        int ni = (*_node_index)[n];
  1.1199 +
  1.1200 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1.1201 +          Node v = _graph.source(e);
  1.1202 +          int vb = _blossom_set->find(v);
  1.1203 +          int vi = (*_node_index)[v];
  1.1204 +
  1.1205 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.1206 +            dualScale * _weight[e];
  1.1207 +
  1.1208 +          if (vb == blossom) {
  1.1209 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.1210 +              _delta3->erase(e);
  1.1211 +            }
  1.1212 +          } else if ((*_blossom_data)[vb].status == EVEN) {
  1.1213 +
  1.1214 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.1215 +              _delta3->erase(e);
  1.1216 +            }
  1.1217 +
  1.1218 +            int vt = _tree_set->find(vb);
  1.1219 +
  1.1220 +            Arc r = _graph.oppositeArc(e);
  1.1221 +
  1.1222 +            typename std::map<int, Arc>::iterator it =
  1.1223 +              (*_node_data)[ni].heap_index.find(vt);
  1.1224 +
  1.1225 +            if (it != (*_node_data)[ni].heap_index.end()) {
  1.1226 +              if ((*_node_data)[ni].heap[it->second] > rw) {
  1.1227 +                (*_node_data)[ni].heap.replace(it->second, r);
  1.1228 +                (*_node_data)[ni].heap.decrease(r, rw);
  1.1229 +                it->second = r;
  1.1230 +              }
  1.1231 +            } else {
  1.1232 +              (*_node_data)[ni].heap.push(r, rw);
  1.1233 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  1.1234 +            }
  1.1235 +
  1.1236 +            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  1.1237 +              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1.1238 +
  1.1239 +              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  1.1240 +                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1.1241 +                             (*_blossom_data)[blossom].offset);
  1.1242 +              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
  1.1243 +                         (*_blossom_data)[blossom].offset){
  1.1244 +                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  1.1245 +                                 (*_blossom_data)[blossom].offset);
  1.1246 +              }
  1.1247 +            }
  1.1248 +
  1.1249 +          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1.1250 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.1251 +              _delta3->erase(e);
  1.1252 +            }
  1.1253 +          }
  1.1254 +        }
  1.1255 +      }
  1.1256 +    }
  1.1257 +
  1.1258 +    void alternatePath(int even, int tree) {
  1.1259 +      int odd;
  1.1260 +
  1.1261 +      evenToMatched(even, tree);
  1.1262 +      (*_blossom_data)[even].status = MATCHED;
  1.1263 +
  1.1264 +      while ((*_blossom_data)[even].pred != INVALID) {
  1.1265 +        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
  1.1266 +        (*_blossom_data)[odd].status = MATCHED;
  1.1267 +        oddToMatched(odd);
  1.1268 +        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
  1.1269 +
  1.1270 +        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
  1.1271 +        (*_blossom_data)[even].status = MATCHED;
  1.1272 +        evenToMatched(even, tree);
  1.1273 +        (*_blossom_data)[even].next =
  1.1274 +          _graph.oppositeArc((*_blossom_data)[odd].pred);
  1.1275 +      }
  1.1276 +
  1.1277 +    }
  1.1278 +
  1.1279 +    void destroyTree(int tree) {
  1.1280 +      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
  1.1281 +        if ((*_blossom_data)[b].status == EVEN) {
  1.1282 +          (*_blossom_data)[b].status = MATCHED;
  1.1283 +          evenToMatched(b, tree);
  1.1284 +        } else if ((*_blossom_data)[b].status == ODD) {
  1.1285 +          (*_blossom_data)[b].status = MATCHED;
  1.1286 +          oddToMatched(b);
  1.1287 +        }
  1.1288 +      }
  1.1289 +      _tree_set->eraseClass(tree);
  1.1290 +    }
  1.1291 +
  1.1292 +
  1.1293 +    void unmatchNode(const Node& node) {
  1.1294 +      int blossom = _blossom_set->find(node);
  1.1295 +      int tree = _tree_set->find(blossom);
  1.1296 +
  1.1297 +      alternatePath(blossom, tree);
  1.1298 +      destroyTree(tree);
  1.1299 +
  1.1300 +      (*_blossom_data)[blossom].status = UNMATCHED;
  1.1301 +      (*_blossom_data)[blossom].base = node;
  1.1302 +      matchedToUnmatched(blossom);
  1.1303 +    }
  1.1304 +
  1.1305 +
  1.1306 +    void augmentOnEdge(const Edge& edge) {
  1.1307 +
  1.1308 +      int left = _blossom_set->find(_graph.u(edge));
  1.1309 +      int right = _blossom_set->find(_graph.v(edge));
  1.1310 +
  1.1311 +      if ((*_blossom_data)[left].status == EVEN) {
  1.1312 +        int left_tree = _tree_set->find(left);
  1.1313 +        alternatePath(left, left_tree);
  1.1314 +        destroyTree(left_tree);
  1.1315 +      } else {
  1.1316 +        (*_blossom_data)[left].status = MATCHED;
  1.1317 +        unmatchedToMatched(left);
  1.1318 +      }
  1.1319 +
  1.1320 +      if ((*_blossom_data)[right].status == EVEN) {
  1.1321 +        int right_tree = _tree_set->find(right);
  1.1322 +        alternatePath(right, right_tree);
  1.1323 +        destroyTree(right_tree);
  1.1324 +      } else {
  1.1325 +        (*_blossom_data)[right].status = MATCHED;
  1.1326 +        unmatchedToMatched(right);
  1.1327 +      }
  1.1328 +
  1.1329 +      (*_blossom_data)[left].next = _graph.direct(edge, true);
  1.1330 +      (*_blossom_data)[right].next = _graph.direct(edge, false);
  1.1331 +    }
  1.1332 +
  1.1333 +    void extendOnArc(const Arc& arc) {
  1.1334 +      int base = _blossom_set->find(_graph.target(arc));
  1.1335 +      int tree = _tree_set->find(base);
  1.1336 +
  1.1337 +      int odd = _blossom_set->find(_graph.source(arc));
  1.1338 +      _tree_set->insert(odd, tree);
  1.1339 +      (*_blossom_data)[odd].status = ODD;
  1.1340 +      matchedToOdd(odd);
  1.1341 +      (*_blossom_data)[odd].pred = arc;
  1.1342 +
  1.1343 +      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
  1.1344 +      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
  1.1345 +      _tree_set->insert(even, tree);
  1.1346 +      (*_blossom_data)[even].status = EVEN;
  1.1347 +      matchedToEven(even, tree);
  1.1348 +    }
  1.1349 +
  1.1350 +    void shrinkOnEdge(const Edge& edge, int tree) {
  1.1351 +      int nca = -1;
  1.1352 +      std::vector<int> left_path, right_path;
  1.1353 +
  1.1354 +      {
  1.1355 +        std::set<int> left_set, right_set;
  1.1356 +        int left = _blossom_set->find(_graph.u(edge));
  1.1357 +        left_path.push_back(left);
  1.1358 +        left_set.insert(left);
  1.1359 +
  1.1360 +        int right = _blossom_set->find(_graph.v(edge));
  1.1361 +        right_path.push_back(right);
  1.1362 +        right_set.insert(right);
  1.1363 +
  1.1364 +        while (true) {
  1.1365 +
  1.1366 +          if ((*_blossom_data)[left].pred == INVALID) break;
  1.1367 +
  1.1368 +          left =
  1.1369 +            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1.1370 +          left_path.push_back(left);
  1.1371 +          left =
  1.1372 +            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1.1373 +          left_path.push_back(left);
  1.1374 +
  1.1375 +          left_set.insert(left);
  1.1376 +
  1.1377 +          if (right_set.find(left) != right_set.end()) {
  1.1378 +            nca = left;
  1.1379 +            break;
  1.1380 +          }
  1.1381 +
  1.1382 +          if ((*_blossom_data)[right].pred == INVALID) break;
  1.1383 +
  1.1384 +          right =
  1.1385 +            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1.1386 +          right_path.push_back(right);
  1.1387 +          right =
  1.1388 +            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1.1389 +          right_path.push_back(right);
  1.1390 +
  1.1391 +          right_set.insert(right);
  1.1392 +
  1.1393 +          if (left_set.find(right) != left_set.end()) {
  1.1394 +            nca = right;
  1.1395 +            break;
  1.1396 +          }
  1.1397 +
  1.1398 +        }
  1.1399 +
  1.1400 +        if (nca == -1) {
  1.1401 +          if ((*_blossom_data)[left].pred == INVALID) {
  1.1402 +            nca = right;
  1.1403 +            while (left_set.find(nca) == left_set.end()) {
  1.1404 +              nca =
  1.1405 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.1406 +              right_path.push_back(nca);
  1.1407 +              nca =
  1.1408 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.1409 +              right_path.push_back(nca);
  1.1410 +            }
  1.1411 +          } else {
  1.1412 +            nca = left;
  1.1413 +            while (right_set.find(nca) == right_set.end()) {
  1.1414 +              nca =
  1.1415 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.1416 +              left_path.push_back(nca);
  1.1417 +              nca =
  1.1418 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.1419 +              left_path.push_back(nca);
  1.1420 +            }
  1.1421 +          }
  1.1422 +        }
  1.1423 +      }
  1.1424 +
  1.1425 +      std::vector<int> subblossoms;
  1.1426 +      Arc prev;
  1.1427 +
  1.1428 +      prev = _graph.direct(edge, true);
  1.1429 +      for (int i = 0; left_path[i] != nca; i += 2) {
  1.1430 +        subblossoms.push_back(left_path[i]);
  1.1431 +        (*_blossom_data)[left_path[i]].next = prev;
  1.1432 +        _tree_set->erase(left_path[i]);
  1.1433 +
  1.1434 +        subblossoms.push_back(left_path[i + 1]);
  1.1435 +        (*_blossom_data)[left_path[i + 1]].status = EVEN;
  1.1436 +        oddToEven(left_path[i + 1], tree);
  1.1437 +        _tree_set->erase(left_path[i + 1]);
  1.1438 +        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
  1.1439 +      }
  1.1440 +
  1.1441 +      int k = 0;
  1.1442 +      while (right_path[k] != nca) ++k;
  1.1443 +
  1.1444 +      subblossoms.push_back(nca);
  1.1445 +      (*_blossom_data)[nca].next = prev;
  1.1446 +
  1.1447 +      for (int i = k - 2; i >= 0; i -= 2) {
  1.1448 +        subblossoms.push_back(right_path[i + 1]);
  1.1449 +        (*_blossom_data)[right_path[i + 1]].status = EVEN;
  1.1450 +        oddToEven(right_path[i + 1], tree);
  1.1451 +        _tree_set->erase(right_path[i + 1]);
  1.1452 +
  1.1453 +        (*_blossom_data)[right_path[i + 1]].next =
  1.1454 +          (*_blossom_data)[right_path[i + 1]].pred;
  1.1455 +
  1.1456 +        subblossoms.push_back(right_path[i]);
  1.1457 +        _tree_set->erase(right_path[i]);
  1.1458 +      }
  1.1459 +
  1.1460 +      int surface =
  1.1461 +        _blossom_set->join(subblossoms.begin(), subblossoms.end());
  1.1462 +
  1.1463 +      for (int i = 0; i < int(subblossoms.size()); ++i) {
  1.1464 +        if (!_blossom_set->trivial(subblossoms[i])) {
  1.1465 +          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
  1.1466 +        }
  1.1467 +        (*_blossom_data)[subblossoms[i]].status = MATCHED;
  1.1468 +      }
  1.1469 +
  1.1470 +      (*_blossom_data)[surface].pot = -2 * _delta_sum;
  1.1471 +      (*_blossom_data)[surface].offset = 0;
  1.1472 +      (*_blossom_data)[surface].status = EVEN;
  1.1473 +      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
  1.1474 +      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
  1.1475 +
  1.1476 +      _tree_set->insert(surface, tree);
  1.1477 +      _tree_set->erase(nca);
  1.1478 +    }
  1.1479 +
  1.1480 +    void splitBlossom(int blossom) {
  1.1481 +      Arc next = (*_blossom_data)[blossom].next;
  1.1482 +      Arc pred = (*_blossom_data)[blossom].pred;
  1.1483 +
  1.1484 +      int tree = _tree_set->find(blossom);
  1.1485 +
  1.1486 +      (*_blossom_data)[blossom].status = MATCHED;
  1.1487 +      oddToMatched(blossom);
  1.1488 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1.1489 +        _delta2->erase(blossom);
  1.1490 +      }
  1.1491 +
  1.1492 +      std::vector<int> subblossoms;
  1.1493 +      _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1.1494 +
  1.1495 +      Value offset = (*_blossom_data)[blossom].offset;
  1.1496 +      int b = _blossom_set->find(_graph.source(pred));
  1.1497 +      int d = _blossom_set->find(_graph.source(next));
  1.1498 +
  1.1499 +      int ib = -1, id = -1;
  1.1500 +      for (int i = 0; i < int(subblossoms.size()); ++i) {
  1.1501 +        if (subblossoms[i] == b) ib = i;
  1.1502 +        if (subblossoms[i] == d) id = i;
  1.1503 +
  1.1504 +        (*_blossom_data)[subblossoms[i]].offset = offset;
  1.1505 +        if (!_blossom_set->trivial(subblossoms[i])) {
  1.1506 +          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
  1.1507 +        }
  1.1508 +        if (_blossom_set->classPrio(subblossoms[i]) !=
  1.1509 +            std::numeric_limits<Value>::max()) {
  1.1510 +          _delta2->push(subblossoms[i],
  1.1511 +                        _blossom_set->classPrio(subblossoms[i]) -
  1.1512 +                        (*_blossom_data)[subblossoms[i]].offset);
  1.1513 +        }
  1.1514 +      }
  1.1515 +
  1.1516 +      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
  1.1517 +        for (int i = (id + 1) % subblossoms.size();
  1.1518 +             i != ib; i = (i + 2) % subblossoms.size()) {
  1.1519 +          int sb = subblossoms[i];
  1.1520 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.1521 +          (*_blossom_data)[sb].next =
  1.1522 +            _graph.oppositeArc((*_blossom_data)[tb].next);
  1.1523 +        }
  1.1524 +
  1.1525 +        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
  1.1526 +          int sb = subblossoms[i];
  1.1527 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.1528 +          int ub = subblossoms[(i + 2) % subblossoms.size()];
  1.1529 +
  1.1530 +          (*_blossom_data)[sb].status = ODD;
  1.1531 +          matchedToOdd(sb);
  1.1532 +          _tree_set->insert(sb, tree);
  1.1533 +          (*_blossom_data)[sb].pred = pred;
  1.1534 +          (*_blossom_data)[sb].next =
  1.1535 +                           _graph.oppositeArc((*_blossom_data)[tb].next);
  1.1536 +
  1.1537 +          pred = (*_blossom_data)[ub].next;
  1.1538 +
  1.1539 +          (*_blossom_data)[tb].status = EVEN;
  1.1540 +          matchedToEven(tb, tree);
  1.1541 +          _tree_set->insert(tb, tree);
  1.1542 +          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
  1.1543 +        }
  1.1544 +
  1.1545 +        (*_blossom_data)[subblossoms[id]].status = ODD;
  1.1546 +        matchedToOdd(subblossoms[id]);
  1.1547 +        _tree_set->insert(subblossoms[id], tree);
  1.1548 +        (*_blossom_data)[subblossoms[id]].next = next;
  1.1549 +        (*_blossom_data)[subblossoms[id]].pred = pred;
  1.1550 +
  1.1551 +      } else {
  1.1552 +
  1.1553 +        for (int i = (ib + 1) % subblossoms.size();
  1.1554 +             i != id; i = (i + 2) % subblossoms.size()) {
  1.1555 +          int sb = subblossoms[i];
  1.1556 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.1557 +          (*_blossom_data)[sb].next =
  1.1558 +            _graph.oppositeArc((*_blossom_data)[tb].next);
  1.1559 +        }
  1.1560 +
  1.1561 +        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
  1.1562 +          int sb = subblossoms[i];
  1.1563 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.1564 +          int ub = subblossoms[(i + 2) % subblossoms.size()];
  1.1565 +
  1.1566 +          (*_blossom_data)[sb].status = ODD;
  1.1567 +          matchedToOdd(sb);
  1.1568 +          _tree_set->insert(sb, tree);
  1.1569 +          (*_blossom_data)[sb].next = next;
  1.1570 +          (*_blossom_data)[sb].pred =
  1.1571 +            _graph.oppositeArc((*_blossom_data)[tb].next);
  1.1572 +
  1.1573 +          (*_blossom_data)[tb].status = EVEN;
  1.1574 +          matchedToEven(tb, tree);
  1.1575 +          _tree_set->insert(tb, tree);
  1.1576 +          (*_blossom_data)[tb].pred =
  1.1577 +            (*_blossom_data)[tb].next =
  1.1578 +            _graph.oppositeArc((*_blossom_data)[ub].next);
  1.1579 +          next = (*_blossom_data)[ub].next;
  1.1580 +        }
  1.1581 +
  1.1582 +        (*_blossom_data)[subblossoms[ib]].status = ODD;
  1.1583 +        matchedToOdd(subblossoms[ib]);
  1.1584 +        _tree_set->insert(subblossoms[ib], tree);
  1.1585 +        (*_blossom_data)[subblossoms[ib]].next = next;
  1.1586 +        (*_blossom_data)[subblossoms[ib]].pred = pred;
  1.1587 +      }
  1.1588 +      _tree_set->erase(blossom);
  1.1589 +    }
  1.1590 +
  1.1591 +    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
  1.1592 +      if (_blossom_set->trivial(blossom)) {
  1.1593 +        int bi = (*_node_index)[base];
  1.1594 +        Value pot = (*_node_data)[bi].pot;
  1.1595 +
  1.1596 +        (*_matching)[base] = matching;
  1.1597 +        _blossom_node_list.push_back(base);
  1.1598 +        (*_node_potential)[base] = pot;
  1.1599 +      } else {
  1.1600 +
  1.1601 +        Value pot = (*_blossom_data)[blossom].pot;
  1.1602 +        int bn = _blossom_node_list.size();
  1.1603 +
  1.1604 +        std::vector<int> subblossoms;
  1.1605 +        _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1.1606 +        int b = _blossom_set->find(base);
  1.1607 +        int ib = -1;
  1.1608 +        for (int i = 0; i < int(subblossoms.size()); ++i) {
  1.1609 +          if (subblossoms[i] == b) { ib = i; break; }
  1.1610 +        }
  1.1611 +
  1.1612 +        for (int i = 1; i < int(subblossoms.size()); i += 2) {
  1.1613 +          int sb = subblossoms[(ib + i) % subblossoms.size()];
  1.1614 +          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
  1.1615 +
  1.1616 +          Arc m = (*_blossom_data)[tb].next;
  1.1617 +          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
  1.1618 +          extractBlossom(tb, _graph.source(m), m);
  1.1619 +        }
  1.1620 +        extractBlossom(subblossoms[ib], base, matching);
  1.1621 +
  1.1622 +        int en = _blossom_node_list.size();
  1.1623 +
  1.1624 +        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
  1.1625 +      }
  1.1626 +    }
  1.1627 +
  1.1628 +    void extractMatching() {
  1.1629 +      std::vector<int> blossoms;
  1.1630 +      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
  1.1631 +        blossoms.push_back(c);
  1.1632 +      }
  1.1633 +
  1.1634 +      for (int i = 0; i < int(blossoms.size()); ++i) {
  1.1635 +        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
  1.1636 +
  1.1637 +          Value offset = (*_blossom_data)[blossoms[i]].offset;
  1.1638 +          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
  1.1639 +          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
  1.1640 +               n != INVALID; ++n) {
  1.1641 +            (*_node_data)[(*_node_index)[n]].pot -= offset;
  1.1642 +          }
  1.1643 +
  1.1644 +          Arc matching = (*_blossom_data)[blossoms[i]].next;
  1.1645 +          Node base = _graph.source(matching);
  1.1646 +          extractBlossom(blossoms[i], base, matching);
  1.1647 +        } else {
  1.1648 +          Node base = (*_blossom_data)[blossoms[i]].base;
  1.1649 +          extractBlossom(blossoms[i], base, INVALID);
  1.1650 +        }
  1.1651 +      }
  1.1652 +    }
  1.1653 +
  1.1654 +  public:
  1.1655 +
  1.1656 +    /// \brief Constructor
  1.1657 +    ///
  1.1658 +    /// Constructor.
  1.1659 +    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
  1.1660 +      : _graph(graph), _weight(weight), _matching(0),
  1.1661 +        _node_potential(0), _blossom_potential(), _blossom_node_list(),
  1.1662 +        _node_num(0), _blossom_num(0),
  1.1663 +
  1.1664 +        _blossom_index(0), _blossom_set(0), _blossom_data(0),
  1.1665 +        _node_index(0), _node_heap_index(0), _node_data(0),
  1.1666 +        _tree_set_index(0), _tree_set(0),
  1.1667 +
  1.1668 +        _delta1_index(0), _delta1(0),
  1.1669 +        _delta2_index(0), _delta2(0),
  1.1670 +        _delta3_index(0), _delta3(0),
  1.1671 +        _delta4_index(0), _delta4(0),
  1.1672 +
  1.1673 +        _delta_sum() {}
  1.1674 +
  1.1675 +    ~MaxWeightedMatching() {
  1.1676 +      destroyStructures();
  1.1677 +    }
  1.1678 +
  1.1679 +    /// \name Execution Control
  1.1680 +    /// The simplest way to execute the algorithm is to use the
  1.1681 +    /// \ref run() member function.
  1.1682 +
  1.1683 +    ///@{
  1.1684 +
  1.1685 +    /// \brief Initialize the algorithm
  1.1686 +    ///
  1.1687 +    /// This function initializes the algorithm.
  1.1688 +    void init() {
  1.1689 +      createStructures();
  1.1690 +
  1.1691 +      for (ArcIt e(_graph); e != INVALID; ++e) {
  1.1692 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
  1.1693 +      }
  1.1694 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1695 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
  1.1696 +      }
  1.1697 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1698 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1.1699 +      }
  1.1700 +      for (int i = 0; i < _blossom_num; ++i) {
  1.1701 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
  1.1702 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
  1.1703 +      }
  1.1704 +
  1.1705 +      int index = 0;
  1.1706 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1707 +        Value max = 0;
  1.1708 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1.1709 +          if (_graph.target(e) == n) continue;
  1.1710 +          if ((dualScale * _weight[e]) / 2 > max) {
  1.1711 +            max = (dualScale * _weight[e]) / 2;
  1.1712 +          }
  1.1713 +        }
  1.1714 +        (*_node_index)[n] = index;
  1.1715 +        (*_node_data)[index].pot = max;
  1.1716 +        _delta1->push(n, max);
  1.1717 +        int blossom =
  1.1718 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
  1.1719 +
  1.1720 +        _tree_set->insert(blossom);
  1.1721 +
  1.1722 +        (*_blossom_data)[blossom].status = EVEN;
  1.1723 +        (*_blossom_data)[blossom].pred = INVALID;
  1.1724 +        (*_blossom_data)[blossom].next = INVALID;
  1.1725 +        (*_blossom_data)[blossom].pot = 0;
  1.1726 +        (*_blossom_data)[blossom].offset = 0;
  1.1727 +        ++index;
  1.1728 +      }
  1.1729 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1730 +        int si = (*_node_index)[_graph.u(e)];
  1.1731 +        int ti = (*_node_index)[_graph.v(e)];
  1.1732 +        if (_graph.u(e) != _graph.v(e)) {
  1.1733 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  1.1734 +                            dualScale * _weight[e]) / 2);
  1.1735 +        }
  1.1736 +      }
  1.1737 +    }
  1.1738 +
  1.1739 +    /// \brief Start the algorithm
  1.1740 +    ///
  1.1741 +    /// This function starts the algorithm.
  1.1742 +    ///
  1.1743 +    /// \pre \ref init() must be called before using this function.
  1.1744 +    void start() {
  1.1745 +      enum OpType {
  1.1746 +        D1, D2, D3, D4
  1.1747 +      };
  1.1748 +
  1.1749 +      int unmatched = _node_num;
  1.1750 +      while (unmatched > 0) {
  1.1751 +        Value d1 = !_delta1->empty() ?
  1.1752 +          _delta1->prio() : std::numeric_limits<Value>::max();
  1.1753 +
  1.1754 +        Value d2 = !_delta2->empty() ?
  1.1755 +          _delta2->prio() : std::numeric_limits<Value>::max();
  1.1756 +
  1.1757 +        Value d3 = !_delta3->empty() ?
  1.1758 +          _delta3->prio() : std::numeric_limits<Value>::max();
  1.1759 +
  1.1760 +        Value d4 = !_delta4->empty() ?
  1.1761 +          _delta4->prio() : std::numeric_limits<Value>::max();
  1.1762 +
  1.1763 +        _delta_sum = d1; OpType ot = D1;
  1.1764 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1.1765 +        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  1.1766 +        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  1.1767 +
  1.1768 +
  1.1769 +        switch (ot) {
  1.1770 +        case D1:
  1.1771 +          {
  1.1772 +            Node n = _delta1->top();
  1.1773 +            unmatchNode(n);
  1.1774 +            --unmatched;
  1.1775 +          }
  1.1776 +          break;
  1.1777 +        case D2:
  1.1778 +          {
  1.1779 +            int blossom = _delta2->top();
  1.1780 +            Node n = _blossom_set->classTop(blossom);
  1.1781 +            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
  1.1782 +            extendOnArc(e);
  1.1783 +          }
  1.1784 +          break;
  1.1785 +        case D3:
  1.1786 +          {
  1.1787 +            Edge e = _delta3->top();
  1.1788 +
  1.1789 +            int left_blossom = _blossom_set->find(_graph.u(e));
  1.1790 +            int right_blossom = _blossom_set->find(_graph.v(e));
  1.1791 +
  1.1792 +            if (left_blossom == right_blossom) {
  1.1793 +              _delta3->pop();
  1.1794 +            } else {
  1.1795 +              int left_tree;
  1.1796 +              if ((*_blossom_data)[left_blossom].status == EVEN) {
  1.1797 +                left_tree = _tree_set->find(left_blossom);
  1.1798 +              } else {
  1.1799 +                left_tree = -1;
  1.1800 +                ++unmatched;
  1.1801 +              }
  1.1802 +              int right_tree;
  1.1803 +              if ((*_blossom_data)[right_blossom].status == EVEN) {
  1.1804 +                right_tree = _tree_set->find(right_blossom);
  1.1805 +              } else {
  1.1806 +                right_tree = -1;
  1.1807 +                ++unmatched;
  1.1808 +              }
  1.1809 +
  1.1810 +              if (left_tree == right_tree) {
  1.1811 +                shrinkOnEdge(e, left_tree);
  1.1812 +              } else {
  1.1813 +                augmentOnEdge(e);
  1.1814 +                unmatched -= 2;
  1.1815 +              }
  1.1816 +            }
  1.1817 +          } break;
  1.1818 +        case D4:
  1.1819 +          splitBlossom(_delta4->top());
  1.1820 +          break;
  1.1821 +        }
  1.1822 +      }
  1.1823 +      extractMatching();
  1.1824 +    }
  1.1825 +
  1.1826 +    /// \brief Run the algorithm.
  1.1827 +    ///
  1.1828 +    /// This method runs the \c %MaxWeightedMatching algorithm.
  1.1829 +    ///
  1.1830 +    /// \note mwm.run() is just a shortcut of the following code.
  1.1831 +    /// \code
  1.1832 +    ///   mwm.init();
  1.1833 +    ///   mwm.start();
  1.1834 +    /// \endcode
  1.1835 +    void run() {
  1.1836 +      init();
  1.1837 +      start();
  1.1838 +    }
  1.1839 +
  1.1840 +    /// @}
  1.1841 +
  1.1842 +    /// \name Primal Solution
  1.1843 +    /// Functions to get the primal solution, i.e. the maximum weighted 
  1.1844 +    /// matching.\n
  1.1845 +    /// Either \ref run() or \ref start() function should be called before
  1.1846 +    /// using them.
  1.1847 +
  1.1848 +    /// @{
  1.1849 +
  1.1850 +    /// \brief Return the weight of the matching.
  1.1851 +    ///
  1.1852 +    /// This function returns the weight of the found matching.
  1.1853 +    ///
  1.1854 +    /// \pre Either run() or start() must be called before using this function.
  1.1855 +    Value matchingWeight() const {
  1.1856 +      Value sum = 0;
  1.1857 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1858 +        if ((*_matching)[n] != INVALID) {
  1.1859 +          sum += _weight[(*_matching)[n]];
  1.1860 +        }
  1.1861 +      }
  1.1862 +      return sum /= 2;
  1.1863 +    }
  1.1864 +
  1.1865 +    /// \brief Return the size (cardinality) of the matching.
  1.1866 +    ///
  1.1867 +    /// This function returns the size (cardinality) of the found matching.
  1.1868 +    ///
  1.1869 +    /// \pre Either run() or start() must be called before using this function.
  1.1870 +    int matchingSize() const {
  1.1871 +      int num = 0;
  1.1872 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1873 +        if ((*_matching)[n] != INVALID) {
  1.1874 +          ++num;
  1.1875 +        }
  1.1876 +      }
  1.1877 +      return num /= 2;
  1.1878 +    }
  1.1879 +
  1.1880 +    /// \brief Return \c true if the given edge is in the matching.
  1.1881 +    ///
  1.1882 +    /// This function returns \c true if the given edge is in the found 
  1.1883 +    /// matching.
  1.1884 +    ///
  1.1885 +    /// \pre Either run() or start() must be called before using this function.
  1.1886 +    bool matching(const Edge& edge) const {
  1.1887 +      return edge == (*_matching)[_graph.u(edge)];
  1.1888 +    }
  1.1889 +
  1.1890 +    /// \brief Return the matching arc (or edge) incident to the given node.
  1.1891 +    ///
  1.1892 +    /// This function returns the matching arc (or edge) incident to the
  1.1893 +    /// given node in the found matching or \c INVALID if the node is 
  1.1894 +    /// not covered by the matching.
  1.1895 +    ///
  1.1896 +    /// \pre Either run() or start() must be called before using this function.
  1.1897 +    Arc matching(const Node& node) const {
  1.1898 +      return (*_matching)[node];
  1.1899 +    }
  1.1900 +
  1.1901 +    /// \brief Return a const reference to the matching map.
  1.1902 +    ///
  1.1903 +    /// This function returns a const reference to a node map that stores
  1.1904 +    /// the matching arc (or edge) incident to each node.
  1.1905 +    const MatchingMap& matchingMap() const {
  1.1906 +      return *_matching;
  1.1907 +    }
  1.1908 +
  1.1909 +    /// \brief Return the mate of the given node.
  1.1910 +    ///
  1.1911 +    /// This function returns the mate of the given node in the found 
  1.1912 +    /// matching or \c INVALID if the node is not covered by the matching.
  1.1913 +    ///
  1.1914 +    /// \pre Either run() or start() must be called before using this function.
  1.1915 +    Node mate(const Node& node) const {
  1.1916 +      return (*_matching)[node] != INVALID ?
  1.1917 +        _graph.target((*_matching)[node]) : INVALID;
  1.1918 +    }
  1.1919 +
  1.1920 +    /// @}
  1.1921 +
  1.1922 +    /// \name Dual Solution
  1.1923 +    /// Functions to get the dual solution.\n
  1.1924 +    /// Either \ref run() or \ref start() function should be called before
  1.1925 +    /// using them.
  1.1926 +
  1.1927 +    /// @{
  1.1928 +
  1.1929 +    /// \brief Return the value of the dual solution.
  1.1930 +    ///
  1.1931 +    /// This function returns the value of the dual solution. 
  1.1932 +    /// It should be equal to the primal value scaled by \ref dualScale 
  1.1933 +    /// "dual scale".
  1.1934 +    ///
  1.1935 +    /// \pre Either run() or start() must be called before using this function.
  1.1936 +    Value dualValue() const {
  1.1937 +      Value sum = 0;
  1.1938 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1939 +        sum += nodeValue(n);
  1.1940 +      }
  1.1941 +      for (int i = 0; i < blossomNum(); ++i) {
  1.1942 +        sum += blossomValue(i) * (blossomSize(i) / 2);
  1.1943 +      }
  1.1944 +      return sum;
  1.1945 +    }
  1.1946 +
  1.1947 +    /// \brief Return the dual value (potential) of the given node.
  1.1948 +    ///
  1.1949 +    /// This function returns the dual value (potential) of the given node.
  1.1950 +    ///
  1.1951 +    /// \pre Either run() or start() must be called before using this function.
  1.1952 +    Value nodeValue(const Node& n) const {
  1.1953 +      return (*_node_potential)[n];
  1.1954 +    }
  1.1955 +
  1.1956 +    /// \brief Return the number of the blossoms in the basis.
  1.1957 +    ///
  1.1958 +    /// This function returns the number of the blossoms in the basis.
  1.1959 +    ///
  1.1960 +    /// \pre Either run() or start() must be called before using this function.
  1.1961 +    /// \see BlossomIt
  1.1962 +    int blossomNum() const {
  1.1963 +      return _blossom_potential.size();
  1.1964 +    }
  1.1965 +
  1.1966 +    /// \brief Return the number of the nodes in the given blossom.
  1.1967 +    ///
  1.1968 +    /// This function returns the number of the nodes in the given blossom.
  1.1969 +    ///
  1.1970 +    /// \pre Either run() or start() must be called before using this function.
  1.1971 +    /// \see BlossomIt
  1.1972 +    int blossomSize(int k) const {
  1.1973 +      return _blossom_potential[k].end - _blossom_potential[k].begin;
  1.1974 +    }
  1.1975 +
  1.1976 +    /// \brief Return the dual value (ptential) of the given blossom.
  1.1977 +    ///
  1.1978 +    /// This function returns the dual value (ptential) of the given blossom.
  1.1979 +    ///
  1.1980 +    /// \pre Either run() or start() must be called before using this function.
  1.1981 +    Value blossomValue(int k) const {
  1.1982 +      return _blossom_potential[k].value;
  1.1983 +    }
  1.1984 +
  1.1985 +    /// \brief Iterator for obtaining the nodes of a blossom.
  1.1986 +    ///
  1.1987 +    /// This class provides an iterator for obtaining the nodes of the 
  1.1988 +    /// given blossom. It lists a subset of the nodes.
  1.1989 +    /// Before using this iterator, you must allocate a 
  1.1990 +    /// MaxWeightedMatching class and execute it.
  1.1991 +    class BlossomIt {
  1.1992 +    public:
  1.1993 +
  1.1994 +      /// \brief Constructor.
  1.1995 +      ///
  1.1996 +      /// Constructor to get the nodes of the given variable.
  1.1997 +      ///
  1.1998 +      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
  1.1999 +      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
  1.2000 +      /// called before initializing this iterator.
  1.2001 +      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
  1.2002 +        : _algorithm(&algorithm)
  1.2003 +      {
  1.2004 +        _index = _algorithm->_blossom_potential[variable].begin;
  1.2005 +        _last = _algorithm->_blossom_potential[variable].end;
  1.2006 +      }
  1.2007 +
  1.2008 +      /// \brief Conversion to \c Node.
  1.2009 +      ///
  1.2010 +      /// Conversion to \c Node.
  1.2011 +      operator Node() const {
  1.2012 +        return _algorithm->_blossom_node_list[_index];
  1.2013 +      }
  1.2014 +
  1.2015 +      /// \brief Increment operator.
  1.2016 +      ///
  1.2017 +      /// Increment operator.
  1.2018 +      BlossomIt& operator++() {
  1.2019 +        ++_index;
  1.2020 +        return *this;
  1.2021 +      }
  1.2022 +
  1.2023 +      /// \brief Validity checking
  1.2024 +      ///
  1.2025 +      /// Checks whether the iterator is invalid.
  1.2026 +      bool operator==(Invalid) const { return _index == _last; }
  1.2027 +
  1.2028 +      /// \brief Validity checking
  1.2029 +      ///
  1.2030 +      /// Checks whether the iterator is valid.
  1.2031 +      bool operator!=(Invalid) const { return _index != _last; }
  1.2032 +
  1.2033 +    private:
  1.2034 +      const MaxWeightedMatching* _algorithm;
  1.2035 +      int _last;
  1.2036 +      int _index;
  1.2037 +    };
  1.2038 +
  1.2039 +    /// @}
  1.2040 +
  1.2041 +  };
  1.2042 +
  1.2043 +  /// \ingroup matching
  1.2044 +  ///
  1.2045 +  /// \brief Weighted perfect matching in general graphs
  1.2046 +  ///
  1.2047 +  /// This class provides an efficient implementation of Edmond's
  1.2048 +  /// maximum weighted perfect matching algorithm. The implementation
  1.2049 +  /// is based on extensive use of priority queues and provides
  1.2050 +  /// \f$O(nm\log n)\f$ time complexity.
  1.2051 +  ///
  1.2052 +  /// The maximum weighted perfect matching problem is to find a subset of 
  1.2053 +  /// the edges in an undirected graph with maximum overall weight for which 
  1.2054 +  /// each node has exactly one incident edge.
  1.2055 +  /// It can be formulated with the following linear program.
  1.2056 +  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1.2057 +  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
  1.2058 +      \quad \forall B\in\mathcal{O}\f] */
  1.2059 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
  1.2060 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
  1.2061 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1.2062 +  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
  1.2063 +  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
  1.2064 +  /// subsets of the nodes.
  1.2065 +  ///
  1.2066 +  /// The algorithm calculates an optimal matching and a proof of the
  1.2067 +  /// optimality. The solution of the dual problem can be used to check
  1.2068 +  /// the result of the algorithm. The dual linear problem is the
  1.2069 +  /// following.
  1.2070 +  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
  1.2071 +      w_{uv} \quad \forall uv\in E\f] */
  1.2072 +  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
  1.2073 +  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  1.2074 +      \frac{\vert B \vert - 1}{2}z_B\f] */
  1.2075 +  ///
  1.2076 +  /// The algorithm can be executed with the run() function. 
  1.2077 +  /// After it the matching (the primal solution) and the dual solution
  1.2078 +  /// can be obtained using the query functions and the 
  1.2079 +  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
  1.2080 +  /// which is able to iterate on the nodes of a blossom. 
  1.2081 +  /// If the value type is integer, then the dual solution is multiplied
  1.2082 +  /// by \ref MaxWeightedMatching::dualScale "4".
  1.2083 +  ///
  1.2084 +  /// \tparam GR The undirected graph type the algorithm runs on.
  1.2085 +  /// \tparam WM The type edge weight map. The default type is 
  1.2086 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  1.2087 +#ifdef DOXYGEN
  1.2088 +  template <typename GR, typename WM>
  1.2089 +#else
  1.2090 +  template <typename GR,
  1.2091 +            typename WM = typename GR::template EdgeMap<int> >
  1.2092 +#endif
  1.2093 +  class MaxWeightedPerfectMatching {
  1.2094 +  public:
  1.2095 +
  1.2096 +    /// The graph type of the algorithm
  1.2097 +    typedef GR Graph;
  1.2098 +    /// The type of the edge weight map
  1.2099 +    typedef WM WeightMap;
  1.2100 +    /// The value type of the edge weights
  1.2101 +    typedef typename WeightMap::Value Value;
  1.2102 +
  1.2103 +    /// \brief Scaling factor for dual solution
  1.2104 +    ///
  1.2105 +    /// Scaling factor for dual solution, it is equal to 4 or 1
  1.2106 +    /// according to the value type.
  1.2107 +    static const int dualScale =
  1.2108 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
  1.2109 +
  1.2110 +    /// The type of the matching map
  1.2111 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
  1.2112 +    MatchingMap;
  1.2113 +
  1.2114 +  private:
  1.2115 +
  1.2116 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
  1.2117 +
  1.2118 +    typedef typename Graph::template NodeMap<Value> NodePotential;
  1.2119 +    typedef std::vector<Node> BlossomNodeList;
  1.2120 +
  1.2121 +    struct BlossomVariable {
  1.2122 +      int begin, end;
  1.2123 +      Value value;
  1.2124 +
  1.2125 +      BlossomVariable(int _begin, int _end, Value _value)
  1.2126 +        : begin(_begin), end(_end), value(_value) {}
  1.2127 +
  1.2128 +    };
  1.2129 +
  1.2130 +    typedef std::vector<BlossomVariable> BlossomPotential;
  1.2131 +
  1.2132 +    const Graph& _graph;
  1.2133 +    const WeightMap& _weight;
  1.2134 +
  1.2135 +    MatchingMap* _matching;
  1.2136 +
  1.2137 +    NodePotential* _node_potential;
  1.2138 +
  1.2139 +    BlossomPotential _blossom_potential;
  1.2140 +    BlossomNodeList _blossom_node_list;
  1.2141 +
  1.2142 +    int _node_num;
  1.2143 +    int _blossom_num;
  1.2144 +
  1.2145 +    typedef RangeMap<int> IntIntMap;
  1.2146 +
  1.2147 +    enum Status {
  1.2148 +      EVEN = -1, MATCHED = 0, ODD = 1
  1.2149 +    };
  1.2150 +
  1.2151 +    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
  1.2152 +    struct BlossomData {
  1.2153 +      int tree;
  1.2154 +      Status status;
  1.2155 +      Arc pred, next;
  1.2156 +      Value pot, offset;
  1.2157 +    };
  1.2158 +
  1.2159 +    IntNodeMap *_blossom_index;
  1.2160 +    BlossomSet *_blossom_set;
  1.2161 +    RangeMap<BlossomData>* _blossom_data;
  1.2162 +
  1.2163 +    IntNodeMap *_node_index;
  1.2164 +    IntArcMap *_node_heap_index;
  1.2165 +
  1.2166 +    struct NodeData {
  1.2167 +
  1.2168 +      NodeData(IntArcMap& node_heap_index)
  1.2169 +        : heap(node_heap_index) {}
  1.2170 +
  1.2171 +      int blossom;
  1.2172 +      Value pot;
  1.2173 +      BinHeap<Value, IntArcMap> heap;
  1.2174 +      std::map<int, Arc> heap_index;
  1.2175 +
  1.2176 +      int tree;
  1.2177 +    };
  1.2178 +
  1.2179 +    RangeMap<NodeData>* _node_data;
  1.2180 +
  1.2181 +    typedef ExtendFindEnum<IntIntMap> TreeSet;
  1.2182 +
  1.2183 +    IntIntMap *_tree_set_index;
  1.2184 +    TreeSet *_tree_set;
  1.2185 +
  1.2186 +    IntIntMap *_delta2_index;
  1.2187 +    BinHeap<Value, IntIntMap> *_delta2;
  1.2188 +
  1.2189 +    IntEdgeMap *_delta3_index;
  1.2190 +    BinHeap<Value, IntEdgeMap> *_delta3;
  1.2191 +
  1.2192 +    IntIntMap *_delta4_index;
  1.2193 +    BinHeap<Value, IntIntMap> *_delta4;
  1.2194 +
  1.2195 +    Value _delta_sum;
  1.2196 +
  1.2197 +    void createStructures() {
  1.2198 +      _node_num = countNodes(_graph);
  1.2199 +      _blossom_num = _node_num * 3 / 2;
  1.2200 +
  1.2201 +      if (!_matching) {
  1.2202 +        _matching = new MatchingMap(_graph);
  1.2203 +      }
  1.2204 +      if (!_node_potential) {
  1.2205 +        _node_potential = new NodePotential(_graph);
  1.2206 +      }
  1.2207 +      if (!_blossom_set) {
  1.2208 +        _blossom_index = new IntNodeMap(_graph);
  1.2209 +        _blossom_set = new BlossomSet(*_blossom_index);
  1.2210 +        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
  1.2211 +      }
  1.2212 +
  1.2213 +      if (!_node_index) {
  1.2214 +        _node_index = new IntNodeMap(_graph);
  1.2215 +        _node_heap_index = new IntArcMap(_graph);
  1.2216 +        _node_data = new RangeMap<NodeData>(_node_num,
  1.2217 +                                            NodeData(*_node_heap_index));
  1.2218 +      }
  1.2219 +
  1.2220 +      if (!_tree_set) {
  1.2221 +        _tree_set_index = new IntIntMap(_blossom_num);
  1.2222 +        _tree_set = new TreeSet(*_tree_set_index);
  1.2223 +      }
  1.2224 +      if (!_delta2) {
  1.2225 +        _delta2_index = new IntIntMap(_blossom_num);
  1.2226 +        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
  1.2227 +      }
  1.2228 +      if (!_delta3) {
  1.2229 +        _delta3_index = new IntEdgeMap(_graph);
  1.2230 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  1.2231 +      }
  1.2232 +      if (!_delta4) {
  1.2233 +        _delta4_index = new IntIntMap(_blossom_num);
  1.2234 +        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
  1.2235 +      }
  1.2236 +    }
  1.2237 +
  1.2238 +    void destroyStructures() {
  1.2239 +      _node_num = countNodes(_graph);
  1.2240 +      _blossom_num = _node_num * 3 / 2;
  1.2241 +
  1.2242 +      if (_matching) {
  1.2243 +        delete _matching;
  1.2244 +      }
  1.2245 +      if (_node_potential) {
  1.2246 +        delete _node_potential;
  1.2247 +      }
  1.2248 +      if (_blossom_set) {
  1.2249 +        delete _blossom_index;
  1.2250 +        delete _blossom_set;
  1.2251 +        delete _blossom_data;
  1.2252 +      }
  1.2253 +
  1.2254 +      if (_node_index) {
  1.2255 +        delete _node_index;
  1.2256 +        delete _node_heap_index;
  1.2257 +        delete _node_data;
  1.2258 +      }
  1.2259 +
  1.2260 +      if (_tree_set) {
  1.2261 +        delete _tree_set_index;
  1.2262 +        delete _tree_set;
  1.2263 +      }
  1.2264 +      if (_delta2) {
  1.2265 +        delete _delta2_index;
  1.2266 +        delete _delta2;
  1.2267 +      }
  1.2268 +      if (_delta3) {
  1.2269 +        delete _delta3_index;
  1.2270 +        delete _delta3;
  1.2271 +      }
  1.2272 +      if (_delta4) {
  1.2273 +        delete _delta4_index;
  1.2274 +        delete _delta4;
  1.2275 +      }
  1.2276 +    }
  1.2277 +
  1.2278 +    void matchedToEven(int blossom, int tree) {
  1.2279 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1.2280 +        _delta2->erase(blossom);
  1.2281 +      }
  1.2282 +
  1.2283 +      if (!_blossom_set->trivial(blossom)) {
  1.2284 +        (*_blossom_data)[blossom].pot -=
  1.2285 +          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
  1.2286 +      }
  1.2287 +
  1.2288 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1.2289 +           n != INVALID; ++n) {
  1.2290 +
  1.2291 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1.2292 +        int ni = (*_node_index)[n];
  1.2293 +
  1.2294 +        (*_node_data)[ni].heap.clear();
  1.2295 +        (*_node_data)[ni].heap_index.clear();
  1.2296 +
  1.2297 +        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
  1.2298 +
  1.2299 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1.2300 +          Node v = _graph.source(e);
  1.2301 +          int vb = _blossom_set->find(v);
  1.2302 +          int vi = (*_node_index)[v];
  1.2303 +
  1.2304 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.2305 +            dualScale * _weight[e];
  1.2306 +
  1.2307 +          if ((*_blossom_data)[vb].status == EVEN) {
  1.2308 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  1.2309 +              _delta3->push(e, rw / 2);
  1.2310 +            }
  1.2311 +          } else {
  1.2312 +            typename std::map<int, Arc>::iterator it =
  1.2313 +              (*_node_data)[vi].heap_index.find(tree);
  1.2314 +
  1.2315 +            if (it != (*_node_data)[vi].heap_index.end()) {
  1.2316 +              if ((*_node_data)[vi].heap[it->second] > rw) {
  1.2317 +                (*_node_data)[vi].heap.replace(it->second, e);
  1.2318 +                (*_node_data)[vi].heap.decrease(e, rw);
  1.2319 +                it->second = e;
  1.2320 +              }
  1.2321 +            } else {
  1.2322 +              (*_node_data)[vi].heap.push(e, rw);
  1.2323 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  1.2324 +            }
  1.2325 +
  1.2326 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  1.2327 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  1.2328 +
  1.2329 +              if ((*_blossom_data)[vb].status == MATCHED) {
  1.2330 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
  1.2331 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
  1.2332 +                               (*_blossom_data)[vb].offset);
  1.2333 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  1.2334 +                           (*_blossom_data)[vb].offset){
  1.2335 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  1.2336 +                                   (*_blossom_data)[vb].offset);
  1.2337 +                }
  1.2338 +              }
  1.2339 +            }
  1.2340 +          }
  1.2341 +        }
  1.2342 +      }
  1.2343 +      (*_blossom_data)[blossom].offset = 0;
  1.2344 +    }
  1.2345 +
  1.2346 +    void matchedToOdd(int blossom) {
  1.2347 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1.2348 +        _delta2->erase(blossom);
  1.2349 +      }
  1.2350 +      (*_blossom_data)[blossom].offset += _delta_sum;
  1.2351 +      if (!_blossom_set->trivial(blossom)) {
  1.2352 +        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
  1.2353 +                     (*_blossom_data)[blossom].offset);
  1.2354 +      }
  1.2355 +    }
  1.2356 +
  1.2357 +    void evenToMatched(int blossom, int tree) {
  1.2358 +      if (!_blossom_set->trivial(blossom)) {
  1.2359 +        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
  1.2360 +      }
  1.2361 +
  1.2362 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1.2363 +           n != INVALID; ++n) {
  1.2364 +        int ni = (*_node_index)[n];
  1.2365 +        (*_node_data)[ni].pot -= _delta_sum;
  1.2366 +
  1.2367 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1.2368 +          Node v = _graph.source(e);
  1.2369 +          int vb = _blossom_set->find(v);
  1.2370 +          int vi = (*_node_index)[v];
  1.2371 +
  1.2372 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.2373 +            dualScale * _weight[e];
  1.2374 +
  1.2375 +          if (vb == blossom) {
  1.2376 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.2377 +              _delta3->erase(e);
  1.2378 +            }
  1.2379 +          } else if ((*_blossom_data)[vb].status == EVEN) {
  1.2380 +
  1.2381 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
  1.2382 +              _delta3->erase(e);
  1.2383 +            }
  1.2384 +
  1.2385 +            int vt = _tree_set->find(vb);
  1.2386 +
  1.2387 +            if (vt != tree) {
  1.2388 +
  1.2389 +              Arc r = _graph.oppositeArc(e);
  1.2390 +
  1.2391 +              typename std::map<int, Arc>::iterator it =
  1.2392 +                (*_node_data)[ni].heap_index.find(vt);
  1.2393 +
  1.2394 +              if (it != (*_node_data)[ni].heap_index.end()) {
  1.2395 +                if ((*_node_data)[ni].heap[it->second] > rw) {
  1.2396 +                  (*_node_data)[ni].heap.replace(it->second, r);
  1.2397 +                  (*_node_data)[ni].heap.decrease(r, rw);
  1.2398 +                  it->second = r;
  1.2399 +                }
  1.2400 +              } else {
  1.2401 +                (*_node_data)[ni].heap.push(r, rw);
  1.2402 +                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  1.2403 +              }
  1.2404 +
  1.2405 +              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  1.2406 +                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1.2407 +
  1.2408 +                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  1.2409 +                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1.2410 +                               (*_blossom_data)[blossom].offset);
  1.2411 +                } else if ((*_delta2)[blossom] >
  1.2412 +                           _blossom_set->classPrio(blossom) -
  1.2413 +                           (*_blossom_data)[blossom].offset){
  1.2414 +                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  1.2415 +                                   (*_blossom_data)[blossom].offset);
  1.2416 +                }
  1.2417 +              }
  1.2418 +            }
  1.2419 +          } else {
  1.2420 +
  1.2421 +            typename std::map<int, Arc>::iterator it =
  1.2422 +              (*_node_data)[vi].heap_index.find(tree);
  1.2423 +
  1.2424 +            if (it != (*_node_data)[vi].heap_index.end()) {
  1.2425 +              (*_node_data)[vi].heap.erase(it->second);
  1.2426 +              (*_node_data)[vi].heap_index.erase(it);
  1.2427 +              if ((*_node_data)[vi].heap.empty()) {
  1.2428 +                _blossom_set->increase(v, std::numeric_limits<Value>::max());
  1.2429 +              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
  1.2430 +                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
  1.2431 +              }
  1.2432 +
  1.2433 +              if ((*_blossom_data)[vb].status == MATCHED) {
  1.2434 +                if (_blossom_set->classPrio(vb) ==
  1.2435 +                    std::numeric_limits<Value>::max()) {
  1.2436 +                  _delta2->erase(vb);
  1.2437 +                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
  1.2438 +                           (*_blossom_data)[vb].offset) {
  1.2439 +                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
  1.2440 +                                   (*_blossom_data)[vb].offset);
  1.2441 +                }
  1.2442 +              }
  1.2443 +            }
  1.2444 +          }
  1.2445 +        }
  1.2446 +      }
  1.2447 +    }
  1.2448 +
  1.2449 +    void oddToMatched(int blossom) {
  1.2450 +      (*_blossom_data)[blossom].offset -= _delta_sum;
  1.2451 +
  1.2452 +      if (_blossom_set->classPrio(blossom) !=
  1.2453 +          std::numeric_limits<Value>::max()) {
  1.2454 +        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1.2455 +                       (*_blossom_data)[blossom].offset);
  1.2456 +      }
  1.2457 +
  1.2458 +      if (!_blossom_set->trivial(blossom)) {
  1.2459 +        _delta4->erase(blossom);
  1.2460 +      }
  1.2461 +    }
  1.2462 +
  1.2463 +    void oddToEven(int blossom, int tree) {
  1.2464 +      if (!_blossom_set->trivial(blossom)) {
  1.2465 +        _delta4->erase(blossom);
  1.2466 +        (*_blossom_data)[blossom].pot -=
  1.2467 +          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
  1.2468 +      }
  1.2469 +
  1.2470 +      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1.2471 +           n != INVALID; ++n) {
  1.2472 +        int ni = (*_node_index)[n];
  1.2473 +
  1.2474 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1.2475 +
  1.2476 +        (*_node_data)[ni].heap.clear();
  1.2477 +        (*_node_data)[ni].heap_index.clear();
  1.2478 +        (*_node_data)[ni].pot +=
  1.2479 +          2 * _delta_sum - (*_blossom_data)[blossom].offset;
  1.2480 +
  1.2481 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1.2482 +          Node v = _graph.source(e);
  1.2483 +          int vb = _blossom_set->find(v);
  1.2484 +          int vi = (*_node_index)[v];
  1.2485 +
  1.2486 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.2487 +            dualScale * _weight[e];
  1.2488 +
  1.2489 +          if ((*_blossom_data)[vb].status == EVEN) {
  1.2490 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  1.2491 +              _delta3->push(e, rw / 2);
  1.2492 +            }
  1.2493 +          } else {
  1.2494 +
  1.2495 +            typename std::map<int, Arc>::iterator it =
  1.2496 +              (*_node_data)[vi].heap_index.find(tree);
  1.2497 +
  1.2498 +            if (it != (*_node_data)[vi].heap_index.end()) {
  1.2499 +              if ((*_node_data)[vi].heap[it->second] > rw) {
  1.2500 +                (*_node_data)[vi].heap.replace(it->second, e);
  1.2501 +                (*_node_data)[vi].heap.decrease(e, rw);
  1.2502 +                it->second = e;
  1.2503 +              }
  1.2504 +            } else {
  1.2505 +              (*_node_data)[vi].heap.push(e, rw);
  1.2506 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  1.2507 +            }
  1.2508 +
  1.2509 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  1.2510 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  1.2511 +
  1.2512 +              if ((*_blossom_data)[vb].status == MATCHED) {
  1.2513 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
  1.2514 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
  1.2515 +                               (*_blossom_data)[vb].offset);
  1.2516 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  1.2517 +                           (*_blossom_data)[vb].offset) {
  1.2518 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  1.2519 +                                   (*_blossom_data)[vb].offset);
  1.2520 +                }
  1.2521 +              }
  1.2522 +            }
  1.2523 +          }
  1.2524 +        }
  1.2525 +      }
  1.2526 +      (*_blossom_data)[blossom].offset = 0;
  1.2527 +    }
  1.2528 +
  1.2529 +    void alternatePath(int even, int tree) {
  1.2530 +      int odd;
  1.2531 +
  1.2532 +      evenToMatched(even, tree);
  1.2533 +      (*_blossom_data)[even].status = MATCHED;
  1.2534 +
  1.2535 +      while ((*_blossom_data)[even].pred != INVALID) {
  1.2536 +        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
  1.2537 +        (*_blossom_data)[odd].status = MATCHED;
  1.2538 +        oddToMatched(odd);
  1.2539 +        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
  1.2540 +
  1.2541 +        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
  1.2542 +        (*_blossom_data)[even].status = MATCHED;
  1.2543 +        evenToMatched(even, tree);
  1.2544 +        (*_blossom_data)[even].next =
  1.2545 +          _graph.oppositeArc((*_blossom_data)[odd].pred);
  1.2546 +      }
  1.2547 +
  1.2548 +    }
  1.2549 +
  1.2550 +    void destroyTree(int tree) {
  1.2551 +      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
  1.2552 +        if ((*_blossom_data)[b].status == EVEN) {
  1.2553 +          (*_blossom_data)[b].status = MATCHED;
  1.2554 +          evenToMatched(b, tree);
  1.2555 +        } else if ((*_blossom_data)[b].status == ODD) {
  1.2556 +          (*_blossom_data)[b].status = MATCHED;
  1.2557 +          oddToMatched(b);
  1.2558 +        }
  1.2559 +      }
  1.2560 +      _tree_set->eraseClass(tree);
  1.2561 +    }
  1.2562 +
  1.2563 +    void augmentOnEdge(const Edge& edge) {
  1.2564 +
  1.2565 +      int left = _blossom_set->find(_graph.u(edge));
  1.2566 +      int right = _blossom_set->find(_graph.v(edge));
  1.2567 +
  1.2568 +      int left_tree = _tree_set->find(left);
  1.2569 +      alternatePath(left, left_tree);
  1.2570 +      destroyTree(left_tree);
  1.2571 +
  1.2572 +      int right_tree = _tree_set->find(right);
  1.2573 +      alternatePath(right, right_tree);
  1.2574 +      destroyTree(right_tree);
  1.2575 +
  1.2576 +      (*_blossom_data)[left].next = _graph.direct(edge, true);
  1.2577 +      (*_blossom_data)[right].next = _graph.direct(edge, false);
  1.2578 +    }
  1.2579 +
  1.2580 +    void extendOnArc(const Arc& arc) {
  1.2581 +      int base = _blossom_set->find(_graph.target(arc));
  1.2582 +      int tree = _tree_set->find(base);
  1.2583 +
  1.2584 +      int odd = _blossom_set->find(_graph.source(arc));
  1.2585 +      _tree_set->insert(odd, tree);
  1.2586 +      (*_blossom_data)[odd].status = ODD;
  1.2587 +      matchedToOdd(odd);
  1.2588 +      (*_blossom_data)[odd].pred = arc;
  1.2589 +
  1.2590 +      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
  1.2591 +      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
  1.2592 +      _tree_set->insert(even, tree);
  1.2593 +      (*_blossom_data)[even].status = EVEN;
  1.2594 +      matchedToEven(even, tree);
  1.2595 +    }
  1.2596 +
  1.2597 +    void shrinkOnEdge(const Edge& edge, int tree) {
  1.2598 +      int nca = -1;
  1.2599 +      std::vector<int> left_path, right_path;
  1.2600 +
  1.2601 +      {
  1.2602 +        std::set<int> left_set, right_set;
  1.2603 +        int left = _blossom_set->find(_graph.u(edge));
  1.2604 +        left_path.push_back(left);
  1.2605 +        left_set.insert(left);
  1.2606 +
  1.2607 +        int right = _blossom_set->find(_graph.v(edge));
  1.2608 +        right_path.push_back(right);
  1.2609 +        right_set.insert(right);
  1.2610 +
  1.2611 +        while (true) {
  1.2612 +
  1.2613 +          if ((*_blossom_data)[left].pred == INVALID) break;
  1.2614 +
  1.2615 +          left =
  1.2616 +            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1.2617 +          left_path.push_back(left);
  1.2618 +          left =
  1.2619 +            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1.2620 +          left_path.push_back(left);
  1.2621 +
  1.2622 +          left_set.insert(left);
  1.2623 +
  1.2624 +          if (right_set.find(left) != right_set.end()) {
  1.2625 +            nca = left;
  1.2626 +            break;
  1.2627 +          }
  1.2628 +
  1.2629 +          if ((*_blossom_data)[right].pred == INVALID) break;
  1.2630 +
  1.2631 +          right =
  1.2632 +            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1.2633 +          right_path.push_back(right);
  1.2634 +          right =
  1.2635 +            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1.2636 +          right_path.push_back(right);
  1.2637 +
  1.2638 +          right_set.insert(right);
  1.2639 +
  1.2640 +          if (left_set.find(right) != left_set.end()) {
  1.2641 +            nca = right;
  1.2642 +            break;
  1.2643 +          }
  1.2644 +
  1.2645 +        }
  1.2646 +
  1.2647 +        if (nca == -1) {
  1.2648 +          if ((*_blossom_data)[left].pred == INVALID) {
  1.2649 +            nca = right;
  1.2650 +            while (left_set.find(nca) == left_set.end()) {
  1.2651 +              nca =
  1.2652 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.2653 +              right_path.push_back(nca);
  1.2654 +              nca =
  1.2655 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.2656 +              right_path.push_back(nca);
  1.2657 +            }
  1.2658 +          } else {
  1.2659 +            nca = left;
  1.2660 +            while (right_set.find(nca) == right_set.end()) {
  1.2661 +              nca =
  1.2662 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.2663 +              left_path.push_back(nca);
  1.2664 +              nca =
  1.2665 +                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1.2666 +              left_path.push_back(nca);
  1.2667 +            }
  1.2668 +          }
  1.2669 +        }
  1.2670 +      }
  1.2671 +
  1.2672 +      std::vector<int> subblossoms;
  1.2673 +      Arc prev;
  1.2674 +
  1.2675 +      prev = _graph.direct(edge, true);
  1.2676 +      for (int i = 0; left_path[i] != nca; i += 2) {
  1.2677 +        subblossoms.push_back(left_path[i]);
  1.2678 +        (*_blossom_data)[left_path[i]].next = prev;
  1.2679 +        _tree_set->erase(left_path[i]);
  1.2680 +
  1.2681 +        subblossoms.push_back(left_path[i + 1]);
  1.2682 +        (*_blossom_data)[left_path[i + 1]].status = EVEN;
  1.2683 +        oddToEven(left_path[i + 1], tree);
  1.2684 +        _tree_set->erase(left_path[i + 1]);
  1.2685 +        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
  1.2686 +      }
  1.2687 +
  1.2688 +      int k = 0;
  1.2689 +      while (right_path[k] != nca) ++k;
  1.2690 +
  1.2691 +      subblossoms.push_back(nca);
  1.2692 +      (*_blossom_data)[nca].next = prev;
  1.2693 +
  1.2694 +      for (int i = k - 2; i >= 0; i -= 2) {
  1.2695 +        subblossoms.push_back(right_path[i + 1]);
  1.2696 +        (*_blossom_data)[right_path[i + 1]].status = EVEN;
  1.2697 +        oddToEven(right_path[i + 1], tree);
  1.2698 +        _tree_set->erase(right_path[i + 1]);
  1.2699 +
  1.2700 +        (*_blossom_data)[right_path[i + 1]].next =
  1.2701 +          (*_blossom_data)[right_path[i + 1]].pred;
  1.2702 +
  1.2703 +        subblossoms.push_back(right_path[i]);
  1.2704 +        _tree_set->erase(right_path[i]);
  1.2705 +      }
  1.2706 +
  1.2707 +      int surface =
  1.2708 +        _blossom_set->join(subblossoms.begin(), subblossoms.end());
  1.2709 +
  1.2710 +      for (int i = 0; i < int(subblossoms.size()); ++i) {
  1.2711 +        if (!_blossom_set->trivial(subblossoms[i])) {
  1.2712 +          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
  1.2713 +        }
  1.2714 +        (*_blossom_data)[subblossoms[i]].status = MATCHED;
  1.2715 +      }
  1.2716 +
  1.2717 +      (*_blossom_data)[surface].pot = -2 * _delta_sum;
  1.2718 +      (*_blossom_data)[surface].offset = 0;
  1.2719 +      (*_blossom_data)[surface].status = EVEN;
  1.2720 +      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
  1.2721 +      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
  1.2722 +
  1.2723 +      _tree_set->insert(surface, tree);
  1.2724 +      _tree_set->erase(nca);
  1.2725 +    }
  1.2726 +
  1.2727 +    void splitBlossom(int blossom) {
  1.2728 +      Arc next = (*_blossom_data)[blossom].next;
  1.2729 +      Arc pred = (*_blossom_data)[blossom].pred;
  1.2730 +
  1.2731 +      int tree = _tree_set->find(blossom);
  1.2732 +
  1.2733 +      (*_blossom_data)[blossom].status = MATCHED;
  1.2734 +      oddToMatched(blossom);
  1.2735 +      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1.2736 +        _delta2->erase(blossom);
  1.2737 +      }
  1.2738 +
  1.2739 +      std::vector<int> subblossoms;
  1.2740 +      _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1.2741 +
  1.2742 +      Value offset = (*_blossom_data)[blossom].offset;
  1.2743 +      int b = _blossom_set->find(_graph.source(pred));
  1.2744 +      int d = _blossom_set->find(_graph.source(next));
  1.2745 +
  1.2746 +      int ib = -1, id = -1;
  1.2747 +      for (int i = 0; i < int(subblossoms.size()); ++i) {
  1.2748 +        if (subblossoms[i] == b) ib = i;
  1.2749 +        if (subblossoms[i] == d) id = i;
  1.2750 +
  1.2751 +        (*_blossom_data)[subblossoms[i]].offset = offset;
  1.2752 +        if (!_blossom_set->trivial(subblossoms[i])) {
  1.2753 +          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
  1.2754 +        }
  1.2755 +        if (_blossom_set->classPrio(subblossoms[i]) !=
  1.2756 +            std::numeric_limits<Value>::max()) {
  1.2757 +          _delta2->push(subblossoms[i],
  1.2758 +                        _blossom_set->classPrio(subblossoms[i]) -
  1.2759 +                        (*_blossom_data)[subblossoms[i]].offset);
  1.2760 +        }
  1.2761 +      }
  1.2762 +
  1.2763 +      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
  1.2764 +        for (int i = (id + 1) % subblossoms.size();
  1.2765 +             i != ib; i = (i + 2) % subblossoms.size()) {
  1.2766 +          int sb = subblossoms[i];
  1.2767 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.2768 +          (*_blossom_data)[sb].next =
  1.2769 +            _graph.oppositeArc((*_blossom_data)[tb].next);
  1.2770 +        }
  1.2771 +
  1.2772 +        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
  1.2773 +          int sb = subblossoms[i];
  1.2774 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.2775 +          int ub = subblossoms[(i + 2) % subblossoms.size()];
  1.2776 +
  1.2777 +          (*_blossom_data)[sb].status = ODD;
  1.2778 +          matchedToOdd(sb);
  1.2779 +          _tree_set->insert(sb, tree);
  1.2780 +          (*_blossom_data)[sb].pred = pred;
  1.2781 +          (*_blossom_data)[sb].next =
  1.2782 +                           _graph.oppositeArc((*_blossom_data)[tb].next);
  1.2783 +
  1.2784 +          pred = (*_blossom_data)[ub].next;
  1.2785 +
  1.2786 +          (*_blossom_data)[tb].status = EVEN;
  1.2787 +          matchedToEven(tb, tree);
  1.2788 +          _tree_set->insert(tb, tree);
  1.2789 +          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
  1.2790 +        }
  1.2791 +
  1.2792 +        (*_blossom_data)[subblossoms[id]].status = ODD;
  1.2793 +        matchedToOdd(subblossoms[id]);
  1.2794 +        _tree_set->insert(subblossoms[id], tree);
  1.2795 +        (*_blossom_data)[subblossoms[id]].next = next;
  1.2796 +        (*_blossom_data)[subblossoms[id]].pred = pred;
  1.2797 +
  1.2798 +      } else {
  1.2799 +
  1.2800 +        for (int i = (ib + 1) % subblossoms.size();
  1.2801 +             i != id; i = (i + 2) % subblossoms.size()) {
  1.2802 +          int sb = subblossoms[i];
  1.2803 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.2804 +          (*_blossom_data)[sb].next =
  1.2805 +            _graph.oppositeArc((*_blossom_data)[tb].next);
  1.2806 +        }
  1.2807 +
  1.2808 +        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
  1.2809 +          int sb = subblossoms[i];
  1.2810 +          int tb = subblossoms[(i + 1) % subblossoms.size()];
  1.2811 +          int ub = subblossoms[(i + 2) % subblossoms.size()];
  1.2812 +
  1.2813 +          (*_blossom_data)[sb].status = ODD;
  1.2814 +          matchedToOdd(sb);
  1.2815 +          _tree_set->insert(sb, tree);
  1.2816 +          (*_blossom_data)[sb].next = next;
  1.2817 +          (*_blossom_data)[sb].pred =
  1.2818 +            _graph.oppositeArc((*_blossom_data)[tb].next);
  1.2819 +
  1.2820 +          (*_blossom_data)[tb].status = EVEN;
  1.2821 +          matchedToEven(tb, tree);
  1.2822 +          _tree_set->insert(tb, tree);
  1.2823 +          (*_blossom_data)[tb].pred =
  1.2824 +            (*_blossom_data)[tb].next =
  1.2825 +            _graph.oppositeArc((*_blossom_data)[ub].next);
  1.2826 +          next = (*_blossom_data)[ub].next;
  1.2827 +        }
  1.2828 +
  1.2829 +        (*_blossom_data)[subblossoms[ib]].status = ODD;
  1.2830 +        matchedToOdd(subblossoms[ib]);
  1.2831 +        _tree_set->insert(subblossoms[ib], tree);
  1.2832 +        (*_blossom_data)[subblossoms[ib]].next = next;
  1.2833 +        (*_blossom_data)[subblossoms[ib]].pred = pred;
  1.2834 +      }
  1.2835 +      _tree_set->erase(blossom);
  1.2836 +    }
  1.2837 +
  1.2838 +    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
  1.2839 +      if (_blossom_set->trivial(blossom)) {
  1.2840 +        int bi = (*_node_index)[base];
  1.2841 +        Value pot = (*_node_data)[bi].pot;
  1.2842 +
  1.2843 +        (*_matching)[base] = matching;
  1.2844 +        _blossom_node_list.push_back(base);
  1.2845 +        (*_node_potential)[base] = pot;
  1.2846 +      } else {
  1.2847 +
  1.2848 +        Value pot = (*_blossom_data)[blossom].pot;
  1.2849 +        int bn = _blossom_node_list.size();
  1.2850 +
  1.2851 +        std::vector<int> subblossoms;
  1.2852 +        _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1.2853 +        int b = _blossom_set->find(base);
  1.2854 +        int ib = -1;
  1.2855 +        for (int i = 0; i < int(subblossoms.size()); ++i) {
  1.2856 +          if (subblossoms[i] == b) { ib = i; break; }
  1.2857 +        }
  1.2858 +
  1.2859 +        for (int i = 1; i < int(subblossoms.size()); i += 2) {
  1.2860 +          int sb = subblossoms[(ib + i) % subblossoms.size()];
  1.2861 +          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
  1.2862 +
  1.2863 +          Arc m = (*_blossom_data)[tb].next;
  1.2864 +          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
  1.2865 +          extractBlossom(tb, _graph.source(m), m);
  1.2866 +        }
  1.2867 +        extractBlossom(subblossoms[ib], base, matching);
  1.2868 +
  1.2869 +        int en = _blossom_node_list.size();
  1.2870 +
  1.2871 +        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
  1.2872 +      }
  1.2873 +    }
  1.2874 +
  1.2875 +    void extractMatching() {
  1.2876 +      std::vector<int> blossoms;
  1.2877 +      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
  1.2878 +        blossoms.push_back(c);
  1.2879 +      }
  1.2880 +
  1.2881 +      for (int i = 0; i < int(blossoms.size()); ++i) {
  1.2882 +
  1.2883 +        Value offset = (*_blossom_data)[blossoms[i]].offset;
  1.2884 +        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
  1.2885 +        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
  1.2886 +             n != INVALID; ++n) {
  1.2887 +          (*_node_data)[(*_node_index)[n]].pot -= offset;
  1.2888 +        }
  1.2889 +
  1.2890 +        Arc matching = (*_blossom_data)[blossoms[i]].next;
  1.2891 +        Node base = _graph.source(matching);
  1.2892 +        extractBlossom(blossoms[i], base, matching);
  1.2893 +      }
  1.2894 +    }
  1.2895 +
  1.2896 +  public:
  1.2897 +
  1.2898 +    /// \brief Constructor
  1.2899 +    ///
  1.2900 +    /// Constructor.
  1.2901 +    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
  1.2902 +      : _graph(graph), _weight(weight), _matching(0),
  1.2903 +        _node_potential(0), _blossom_potential(), _blossom_node_list(),
  1.2904 +        _node_num(0), _blossom_num(0),
  1.2905 +
  1.2906 +        _blossom_index(0), _blossom_set(0), _blossom_data(0),
  1.2907 +        _node_index(0), _node_heap_index(0), _node_data(0),
  1.2908 +        _tree_set_index(0), _tree_set(0),
  1.2909 +
  1.2910 +        _delta2_index(0), _delta2(0),
  1.2911 +        _delta3_index(0), _delta3(0),
  1.2912 +        _delta4_index(0), _delta4(0),
  1.2913 +
  1.2914 +        _delta_sum() {}
  1.2915 +
  1.2916 +    ~MaxWeightedPerfectMatching() {
  1.2917 +      destroyStructures();
  1.2918 +    }
  1.2919 +
  1.2920 +    /// \name Execution Control
  1.2921 +    /// The simplest way to execute the algorithm is to use the
  1.2922 +    /// \ref run() member function.
  1.2923 +
  1.2924 +    ///@{
  1.2925 +
  1.2926 +    /// \brief Initialize the algorithm
  1.2927 +    ///
  1.2928 +    /// This function initializes the algorithm.
  1.2929 +    void init() {
  1.2930 +      createStructures();
  1.2931 +
  1.2932 +      for (ArcIt e(_graph); e != INVALID; ++e) {
  1.2933 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
  1.2934 +      }
  1.2935 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.2936 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1.2937 +      }
  1.2938 +      for (int i = 0; i < _blossom_num; ++i) {
  1.2939 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
  1.2940 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
  1.2941 +      }
  1.2942 +
  1.2943 +      int index = 0;
  1.2944 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.2945 +        Value max = - std::numeric_limits<Value>::max();
  1.2946 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1.2947 +          if (_graph.target(e) == n) continue;
  1.2948 +          if ((dualScale * _weight[e]) / 2 > max) {
  1.2949 +            max = (dualScale * _weight[e]) / 2;
  1.2950 +          }
  1.2951 +        }
  1.2952 +        (*_node_index)[n] = index;
  1.2953 +        (*_node_data)[index].pot = max;
  1.2954 +        int blossom =
  1.2955 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
  1.2956 +
  1.2957 +        _tree_set->insert(blossom);
  1.2958 +
  1.2959 +        (*_blossom_data)[blossom].status = EVEN;
  1.2960 +        (*_blossom_data)[blossom].pred = INVALID;
  1.2961 +        (*_blossom_data)[blossom].next = INVALID;
  1.2962 +        (*_blossom_data)[blossom].pot = 0;
  1.2963 +        (*_blossom_data)[blossom].offset = 0;
  1.2964 +        ++index;
  1.2965 +      }
  1.2966 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.2967 +        int si = (*_node_index)[_graph.u(e)];
  1.2968 +        int ti = (*_node_index)[_graph.v(e)];
  1.2969 +        if (_graph.u(e) != _graph.v(e)) {
  1.2970 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  1.2971 +                            dualScale * _weight[e]) / 2);
  1.2972 +        }
  1.2973 +      }
  1.2974 +    }
  1.2975 +
  1.2976 +    /// \brief Start the algorithm
  1.2977 +    ///
  1.2978 +    /// This function starts the algorithm.
  1.2979 +    ///
  1.2980 +    /// \pre \ref init() must be called before using this function.
  1.2981 +    bool start() {
  1.2982 +      enum OpType {
  1.2983 +        D2, D3, D4
  1.2984 +      };
  1.2985 +
  1.2986 +      int unmatched = _node_num;
  1.2987 +      while (unmatched > 0) {
  1.2988 +        Value d2 = !_delta2->empty() ?
  1.2989 +          _delta2->prio() : std::numeric_limits<Value>::max();
  1.2990 +
  1.2991 +        Value d3 = !_delta3->empty() ?
  1.2992 +          _delta3->prio() : std::numeric_limits<Value>::max();
  1.2993 +
  1.2994 +        Value d4 = !_delta4->empty() ?
  1.2995 +          _delta4->prio() : std::numeric_limits<Value>::max();
  1.2996 +
  1.2997 +        _delta_sum = d2; OpType ot = D2;
  1.2998 +        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  1.2999 +        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  1.3000 +
  1.3001 +        if (_delta_sum == std::numeric_limits<Value>::max()) {
  1.3002 +          return false;
  1.3003 +        }
  1.3004 +
  1.3005 +        switch (ot) {
  1.3006 +        case D2:
  1.3007 +          {
  1.3008 +            int blossom = _delta2->top();
  1.3009 +            Node n = _blossom_set->classTop(blossom);
  1.3010 +            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
  1.3011 +            extendOnArc(e);
  1.3012 +          }
  1.3013 +          break;
  1.3014 +        case D3:
  1.3015 +          {
  1.3016 +            Edge e = _delta3->top();
  1.3017 +
  1.3018 +            int left_blossom = _blossom_set->find(_graph.u(e));
  1.3019 +            int right_blossom = _blossom_set->find(_graph.v(e));
  1.3020 +
  1.3021 +            if (left_blossom == right_blossom) {
  1.3022 +              _delta3->pop();
  1.3023 +            } else {
  1.3024 +              int left_tree = _tree_set->find(left_blossom);
  1.3025 +              int right_tree = _tree_set->find(right_blossom);
  1.3026 +
  1.3027 +              if (left_tree == right_tree) {
  1.3028 +                shrinkOnEdge(e, left_tree);
  1.3029 +              } else {
  1.3030 +                augmentOnEdge(e);
  1.3031 +                unmatched -= 2;
  1.3032 +              }
  1.3033 +            }
  1.3034 +          } break;
  1.3035 +        case D4:
  1.3036 +          splitBlossom(_delta4->top());
  1.3037 +          break;
  1.3038 +        }
  1.3039 +      }
  1.3040 +      extractMatching();
  1.3041 +      return true;
  1.3042 +    }
  1.3043 +
  1.3044 +    /// \brief Run the algorithm.
  1.3045 +    ///
  1.3046 +    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
  1.3047 +    ///
  1.3048 +    /// \note mwpm.run() is just a shortcut of the following code.
  1.3049 +    /// \code
  1.3050 +    ///   mwpm.init();
  1.3051 +    ///   mwpm.start();
  1.3052 +    /// \endcode
  1.3053 +    bool run() {
  1.3054 +      init();
  1.3055 +      return start();
  1.3056 +    }
  1.3057 +
  1.3058 +    /// @}
  1.3059 +
  1.3060 +    /// \name Primal Solution
  1.3061 +    /// Functions to get the primal solution, i.e. the maximum weighted 
  1.3062 +    /// perfect matching.\n
  1.3063 +    /// Either \ref run() or \ref start() function should be called before
  1.3064 +    /// using them.
  1.3065 +
  1.3066 +    /// @{
  1.3067 +
  1.3068 +    /// \brief Return the weight of the matching.
  1.3069 +    ///
  1.3070 +    /// This function returns the weight of the found matching.
  1.3071 +    ///
  1.3072 +    /// \pre Either run() or start() must be called before using this function.
  1.3073 +    Value matchingWeight() const {
  1.3074 +      Value sum = 0;
  1.3075 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.3076 +        if ((*_matching)[n] != INVALID) {
  1.3077 +          sum += _weight[(*_matching)[n]];
  1.3078 +        }
  1.3079 +      }
  1.3080 +      return sum /= 2;
  1.3081 +    }
  1.3082 +
  1.3083 +    /// \brief Return \c true if the given edge is in the matching.
  1.3084 +    ///
  1.3085 +    /// This function returns \c true if the given edge is in the found 
  1.3086 +    /// matching.
  1.3087 +    ///
  1.3088 +    /// \pre Either run() or start() must be called before using this function.
  1.3089 +    bool matching(const Edge& edge) const {
  1.3090 +      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
  1.3091 +    }
  1.3092 +
  1.3093 +    /// \brief Return the matching arc (or edge) incident to the given node.
  1.3094 +    ///
  1.3095 +    /// This function returns the matching arc (or edge) incident to the
  1.3096 +    /// given node in the found matching or \c INVALID if the node is 
  1.3097 +    /// not covered by the matching.
  1.3098 +    ///
  1.3099 +    /// \pre Either run() or start() must be called before using this function.
  1.3100 +    Arc matching(const Node& node) const {
  1.3101 +      return (*_matching)[node];
  1.3102 +    }
  1.3103 +
  1.3104 +    /// \brief Return a const reference to the matching map.
  1.3105 +    ///
  1.3106 +    /// This function returns a const reference to a node map that stores
  1.3107 +    /// the matching arc (or edge) incident to each node.
  1.3108 +    const MatchingMap& matchingMap() const {
  1.3109 +      return *_matching;
  1.3110 +    }
  1.3111 +
  1.3112 +    /// \brief Return the mate of the given node.
  1.3113 +    ///
  1.3114 +    /// This function returns the mate of the given node in the found 
  1.3115 +    /// matching or \c INVALID if the node is not covered by the matching.
  1.3116 +    ///
  1.3117 +    /// \pre Either run() or start() must be called before using this function.
  1.3118 +    Node mate(const Node& node) const {
  1.3119 +      return _graph.target((*_matching)[node]);
  1.3120 +    }
  1.3121 +
  1.3122 +    /// @}
  1.3123 +
  1.3124 +    /// \name Dual Solution
  1.3125 +    /// Functions to get the dual solution.\n
  1.3126 +    /// Either \ref run() or \ref start() function should be called before
  1.3127 +    /// using them.
  1.3128 +
  1.3129 +    /// @{
  1.3130 +
  1.3131 +    /// \brief Return the value of the dual solution.
  1.3132 +    ///
  1.3133 +    /// This function returns the value of the dual solution. 
  1.3134 +    /// It should be equal to the primal value scaled by \ref dualScale 
  1.3135 +    /// "dual scale".
  1.3136 +    ///
  1.3137 +    /// \pre Either run() or start() must be called before using this function.
  1.3138 +    Value dualValue() const {
  1.3139 +      Value sum = 0;
  1.3140 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.3141 +        sum += nodeValue(n);
  1.3142 +      }
  1.3143 +      for (int i = 0; i < blossomNum(); ++i) {
  1.3144 +        sum += blossomValue(i) * (blossomSize(i) / 2);
  1.3145 +      }
  1.3146 +      return sum;
  1.3147 +    }
  1.3148 +
  1.3149 +    /// \brief Return the dual value (potential) of the given node.
  1.3150 +    ///
  1.3151 +    /// This function returns the dual value (potential) of the given node.
  1.3152 +    ///
  1.3153 +    /// \pre Either run() or start() must be called before using this function.
  1.3154 +    Value nodeValue(const Node& n) const {
  1.3155 +      return (*_node_potential)[n];
  1.3156 +    }
  1.3157 +
  1.3158 +    /// \brief Return the number of the blossoms in the basis.
  1.3159 +    ///
  1.3160 +    /// This function returns the number of the blossoms in the basis.
  1.3161 +    ///
  1.3162 +    /// \pre Either run() or start() must be called before using this function.
  1.3163 +    /// \see BlossomIt
  1.3164 +    int blossomNum() const {
  1.3165 +      return _blossom_potential.size();
  1.3166 +    }
  1.3167 +
  1.3168 +    /// \brief Return the number of the nodes in the given blossom.
  1.3169 +    ///
  1.3170 +    /// This function returns the number of the nodes in the given blossom.
  1.3171 +    ///
  1.3172 +    /// \pre Either run() or start() must be called before using this function.
  1.3173 +    /// \see BlossomIt
  1.3174 +    int blossomSize(int k) const {
  1.3175 +      return _blossom_potential[k].end - _blossom_potential[k].begin;
  1.3176 +    }
  1.3177 +
  1.3178 +    /// \brief Return the dual value (ptential) of the given blossom.
  1.3179 +    ///
  1.3180 +    /// This function returns the dual value (ptential) of the given blossom.
  1.3181 +    ///
  1.3182 +    /// \pre Either run() or start() must be called before using this function.
  1.3183 +    Value blossomValue(int k) const {
  1.3184 +      return _blossom_potential[k].value;
  1.3185 +    }
  1.3186 +
  1.3187 +    /// \brief Iterator for obtaining the nodes of a blossom.
  1.3188 +    ///
  1.3189 +    /// This class provides an iterator for obtaining the nodes of the 
  1.3190 +    /// given blossom. It lists a subset of the nodes.
  1.3191 +    /// Before using this iterator, you must allocate a 
  1.3192 +    /// MaxWeightedPerfectMatching class and execute it.
  1.3193 +    class BlossomIt {
  1.3194 +    public:
  1.3195 +
  1.3196 +      /// \brief Constructor.
  1.3197 +      ///
  1.3198 +      /// Constructor to get the nodes of the given variable.
  1.3199 +      ///
  1.3200 +      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
  1.3201 +      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
  1.3202 +      /// must be called before initializing this iterator.
  1.3203 +      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
  1.3204 +        : _algorithm(&algorithm)
  1.3205 +      {
  1.3206 +        _index = _algorithm->_blossom_potential[variable].begin;
  1.3207 +        _last = _algorithm->_blossom_potential[variable].end;
  1.3208 +      }
  1.3209 +
  1.3210 +      /// \brief Conversion to \c Node.
  1.3211 +      ///
  1.3212 +      /// Conversion to \c Node.
  1.3213 +      operator Node() const {
  1.3214 +        return _algorithm->_blossom_node_list[_index];
  1.3215 +      }
  1.3216 +
  1.3217 +      /// \brief Increment operator.
  1.3218 +      ///
  1.3219 +      /// Increment operator.
  1.3220 +      BlossomIt& operator++() {
  1.3221 +        ++_index;
  1.3222 +        return *this;
  1.3223 +      }
  1.3224 +
  1.3225 +      /// \brief Validity checking
  1.3226 +      ///
  1.3227 +      /// This function checks whether the iterator is invalid.
  1.3228 +      bool operator==(Invalid) const { return _index == _last; }
  1.3229 +
  1.3230 +      /// \brief Validity checking
  1.3231 +      ///
  1.3232 +      /// This function checks whether the iterator is valid.
  1.3233 +      bool operator!=(Invalid) const { return _index != _last; }
  1.3234 +
  1.3235 +    private:
  1.3236 +      const MaxWeightedPerfectMatching* _algorithm;
  1.3237 +      int _last;
  1.3238 +      int _index;
  1.3239 +    };
  1.3240 +
  1.3241 +    /// @}
  1.3242 +
  1.3243 +  };
  1.3244 +
  1.3245 +} //END OF NAMESPACE LEMON
  1.3246 +
  1.3247 +#endif //LEMON_MAX_MATCHING_H