lemon/max_matching.h
changeset 594 d657c71db7db
parent 593 7ac52d6a268e
child 595 16d7255a6849
     1.1 --- a/lemon/max_matching.h	Fri Apr 17 09:54:14 2009 +0200
     1.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.3 @@ -1,3244 +0,0 @@
     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