Merge #314
authorAlpar Juttner <alpar@cs.elte.hu>
Tue, 16 Mar 2010 21:27:35 +0100
changeset 953d8ea85825e02
parent 952 23da1b807280
parent 951 41d7ac528c3a
child 954 07ec2b52e53d
Merge #314
doc/groups.dox
lemon/Makefile.am
test/CMakeLists.txt
test/Makefile.am
     1.1 --- a/doc/groups.dox	Tue Mar 16 21:18:39 2010 +0100
     1.2 +++ b/doc/groups.dox	Tue Mar 16 21:27:35 2010 +0100
     1.3 @@ -386,7 +386,7 @@
     1.4  also provide functions to query the minimum cut, which is the dual
     1.5  problem of maximum flow.
     1.6  
     1.7 -\ref Circulation is a preflow push-relabel algorithm implemented directly 
     1.8 +\ref Circulation is a preflow push-relabel algorithm implemented directly
     1.9  for finding feasible circulations, which is a somewhat different problem,
    1.10  but it is strongly related to maximum flow.
    1.11  For more information, see \ref Circulation.
    1.12 @@ -522,6 +522,13 @@
    1.13  - \ref MaxWeightedPerfectMatching
    1.14    Edmond's blossom shrinking algorithm for calculating maximum weighted
    1.15    perfect matching in general graphs.
    1.16 +- \ref MaxFractionalMatching Push-relabel algorithm for calculating
    1.17 +  maximum cardinality fractional matching in general graphs.
    1.18 +- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
    1.19 +  maximum weighted fractional matching in general graphs.
    1.20 +- \ref MaxWeightedPerfectFractionalMatching
    1.21 +  Augmenting path algorithm for calculating maximum weighted
    1.22 +  perfect fractional matching in general graphs.
    1.23  
    1.24  \image html matching.png
    1.25  \image latex matching.eps "Min Cost Perfect Matching" width=\textwidth
     2.1 --- a/lemon/Makefile.am	Tue Mar 16 21:18:39 2010 +0100
     2.2 +++ b/lemon/Makefile.am	Tue Mar 16 21:27:35 2010 +0100
     2.3 @@ -84,6 +84,7 @@
     2.4  	lemon/error.h \
     2.5  	lemon/euler.h \
     2.6  	lemon/fib_heap.h \
     2.7 +	lemon/fractional_matching.h \
     2.8  	lemon/full_graph.h \
     2.9  	lemon/glpk.h \
    2.10  	lemon/gomory_hu.h \
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/lemon/fractional_matching.h	Tue Mar 16 21:27:35 2010 +0100
     3.3 @@ -0,0 +1,2130 @@
     3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library.
     3.7 + *
     3.8 + * Copyright (C) 2003-2009
     3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    3.11 + *
    3.12 + * Permission to use, modify and distribute this software is granted
    3.13 + * provided that this copyright notice appears in all copies. For
    3.14 + * precise terms see the accompanying LICENSE file.
    3.15 + *
    3.16 + * This software is provided "AS IS" with no warranty of any kind,
    3.17 + * express or implied, and with no claim as to its suitability for any
    3.18 + * purpose.
    3.19 + *
    3.20 + */
    3.21 +
    3.22 +#ifndef LEMON_FRACTIONAL_MATCHING_H
    3.23 +#define LEMON_FRACTIONAL_MATCHING_H
    3.24 +
    3.25 +#include <vector>
    3.26 +#include <queue>
    3.27 +#include <set>
    3.28 +#include <limits>
    3.29 +
    3.30 +#include <lemon/core.h>
    3.31 +#include <lemon/unionfind.h>
    3.32 +#include <lemon/bin_heap.h>
    3.33 +#include <lemon/maps.h>
    3.34 +#include <lemon/assert.h>
    3.35 +#include <lemon/elevator.h>
    3.36 +
    3.37 +///\ingroup matching
    3.38 +///\file
    3.39 +///\brief Fractional matching algorithms in general graphs.
    3.40 +
    3.41 +namespace lemon {
    3.42 +
    3.43 +  /// \brief Default traits class of MaxFractionalMatching class.
    3.44 +  ///
    3.45 +  /// Default traits class of MaxFractionalMatching class.
    3.46 +  /// \tparam GR Graph type.
    3.47 +  template <typename GR>
    3.48 +  struct MaxFractionalMatchingDefaultTraits {
    3.49 +
    3.50 +    /// \brief The type of the graph the algorithm runs on.
    3.51 +    typedef GR Graph;
    3.52 +
    3.53 +    /// \brief The type of the map that stores the matching.
    3.54 +    ///
    3.55 +    /// The type of the map that stores the matching arcs.
    3.56 +    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
    3.57 +    typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
    3.58 +
    3.59 +    /// \brief Instantiates a MatchingMap.
    3.60 +    ///
    3.61 +    /// This function instantiates a \ref MatchingMap.
    3.62 +    /// \param graph The graph for which we would like to define
    3.63 +    /// the matching map.
    3.64 +    static MatchingMap* createMatchingMap(const Graph& graph) {
    3.65 +      return new MatchingMap(graph);
    3.66 +    }
    3.67 +
    3.68 +    /// \brief The elevator type used by MaxFractionalMatching algorithm.
    3.69 +    ///
    3.70 +    /// The elevator type used by MaxFractionalMatching algorithm.
    3.71 +    ///
    3.72 +    /// \sa Elevator
    3.73 +    /// \sa LinkedElevator
    3.74 +    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
    3.75 +
    3.76 +    /// \brief Instantiates an Elevator.
    3.77 +    ///
    3.78 +    /// This function instantiates an \ref Elevator.
    3.79 +    /// \param graph The graph for which we would like to define
    3.80 +    /// the elevator.
    3.81 +    /// \param max_level The maximum level of the elevator.
    3.82 +    static Elevator* createElevator(const Graph& graph, int max_level) {
    3.83 +      return new Elevator(graph, max_level);
    3.84 +    }
    3.85 +  };
    3.86 +
    3.87 +  /// \ingroup matching
    3.88 +  ///
    3.89 +  /// \brief Max cardinality fractional matching
    3.90 +  ///
    3.91 +  /// This class provides an implementation of fractional matching
    3.92 +  /// algorithm based on push-relabel principle.
    3.93 +  ///
    3.94 +  /// The maximum cardinality fractional matching is a relaxation of the
    3.95 +  /// maximum cardinality matching problem where the odd set constraints
    3.96 +  /// are omitted.
    3.97 +  /// It can be formulated with the following linear program.
    3.98 +  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
    3.99 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
   3.100 +  /// \f[\max \sum_{e\in E}x_e\f]
   3.101 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   3.102 +  /// \f$X\f$. The result can be represented as the union of a
   3.103 +  /// matching with one value edges and a set of odd length cycles
   3.104 +  /// with half value edges.
   3.105 +  ///
   3.106 +  /// The algorithm calculates an optimal fractional matching and a
   3.107 +  /// barrier. The number of adjacents of any node set minus the size
   3.108 +  /// of node set is a lower bound on the uncovered nodes in the
   3.109 +  /// graph. For maximum matching a barrier is computed which
   3.110 +  /// maximizes this difference.
   3.111 +  ///
   3.112 +  /// The algorithm can be executed with the run() function.  After it
   3.113 +  /// the matching (the primal solution) and the barrier (the dual
   3.114 +  /// solution) can be obtained using the query functions.
   3.115 +  ///
   3.116 +  /// The primal solution is multiplied by
   3.117 +  /// \ref MaxFractionalMatching::primalScale "2".
   3.118 +  ///
   3.119 +  /// \tparam GR The undirected graph type the algorithm runs on.
   3.120 +#ifdef DOXYGEN
   3.121 +  template <typename GR, typename TR>
   3.122 +#else
   3.123 +  template <typename GR,
   3.124 +            typename TR = MaxFractionalMatchingDefaultTraits<GR> >
   3.125 +#endif
   3.126 +  class MaxFractionalMatching {
   3.127 +  public:
   3.128 +
   3.129 +    /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
   3.130 +    /// class" of the algorithm.
   3.131 +    typedef TR Traits;
   3.132 +    /// The type of the graph the algorithm runs on.
   3.133 +    typedef typename TR::Graph Graph;
   3.134 +    /// The type of the matching map.
   3.135 +    typedef typename TR::MatchingMap MatchingMap;
   3.136 +    /// The type of the elevator.
   3.137 +    typedef typename TR::Elevator Elevator;
   3.138 +
   3.139 +    /// \brief Scaling factor for primal solution
   3.140 +    ///
   3.141 +    /// Scaling factor for primal solution.
   3.142 +    static const int primalScale = 2;
   3.143 +
   3.144 +  private:
   3.145 +
   3.146 +    const Graph &_graph;
   3.147 +    int _node_num;
   3.148 +    bool _allow_loops;
   3.149 +    int _empty_level;
   3.150 +
   3.151 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
   3.152 +
   3.153 +    bool _local_matching;
   3.154 +    MatchingMap *_matching;
   3.155 +
   3.156 +    bool _local_level;
   3.157 +    Elevator *_level;
   3.158 +
   3.159 +    typedef typename Graph::template NodeMap<int> InDegMap;
   3.160 +    InDegMap *_indeg;
   3.161 +
   3.162 +    void createStructures() {
   3.163 +      _node_num = countNodes(_graph);
   3.164 +
   3.165 +      if (!_matching) {
   3.166 +        _local_matching = true;
   3.167 +        _matching = Traits::createMatchingMap(_graph);
   3.168 +      }
   3.169 +      if (!_level) {
   3.170 +        _local_level = true;
   3.171 +        _level = Traits::createElevator(_graph, _node_num);
   3.172 +      }
   3.173 +      if (!_indeg) {
   3.174 +        _indeg = new InDegMap(_graph);
   3.175 +      }
   3.176 +    }
   3.177 +
   3.178 +    void destroyStructures() {
   3.179 +      if (_local_matching) {
   3.180 +        delete _matching;
   3.181 +      }
   3.182 +      if (_local_level) {
   3.183 +        delete _level;
   3.184 +      }
   3.185 +      if (_indeg) {
   3.186 +        delete _indeg;
   3.187 +      }
   3.188 +    }
   3.189 +
   3.190 +    void postprocessing() {
   3.191 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.192 +        if ((*_indeg)[n] != 0) continue;
   3.193 +        _indeg->set(n, -1);
   3.194 +        Node u = n;
   3.195 +        while ((*_matching)[u] != INVALID) {
   3.196 +          Node v = _graph.target((*_matching)[u]);
   3.197 +          _indeg->set(v, -1);
   3.198 +          Arc a = _graph.oppositeArc((*_matching)[u]);
   3.199 +          u = _graph.target((*_matching)[v]);
   3.200 +          _indeg->set(u, -1);
   3.201 +          _matching->set(v, a);
   3.202 +        }
   3.203 +      }
   3.204 +
   3.205 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.206 +        if ((*_indeg)[n] != 1) continue;
   3.207 +        _indeg->set(n, -1);
   3.208 +
   3.209 +        int num = 1;
   3.210 +        Node u = _graph.target((*_matching)[n]);
   3.211 +        while (u != n) {
   3.212 +          _indeg->set(u, -1);
   3.213 +          u = _graph.target((*_matching)[u]);
   3.214 +          ++num;
   3.215 +        }
   3.216 +        if (num % 2 == 0 && num > 2) {
   3.217 +          Arc prev = _graph.oppositeArc((*_matching)[n]);
   3.218 +          Node v = _graph.target((*_matching)[n]);
   3.219 +          u = _graph.target((*_matching)[v]);
   3.220 +          _matching->set(v, prev);
   3.221 +          while (u != n) {
   3.222 +            prev = _graph.oppositeArc((*_matching)[u]);
   3.223 +            v = _graph.target((*_matching)[u]);
   3.224 +            u = _graph.target((*_matching)[v]);
   3.225 +            _matching->set(v, prev);
   3.226 +          }
   3.227 +        }
   3.228 +      }
   3.229 +    }
   3.230 +
   3.231 +  public:
   3.232 +
   3.233 +    typedef MaxFractionalMatching Create;
   3.234 +
   3.235 +    ///\name Named Template Parameters
   3.236 +
   3.237 +    ///@{
   3.238 +
   3.239 +    template <typename T>
   3.240 +    struct SetMatchingMapTraits : public Traits {
   3.241 +      typedef T MatchingMap;
   3.242 +      static MatchingMap *createMatchingMap(const Graph&) {
   3.243 +        LEMON_ASSERT(false, "MatchingMap is not initialized");
   3.244 +        return 0; // ignore warnings
   3.245 +      }
   3.246 +    };
   3.247 +
   3.248 +    /// \brief \ref named-templ-param "Named parameter" for setting
   3.249 +    /// MatchingMap type
   3.250 +    ///
   3.251 +    /// \ref named-templ-param "Named parameter" for setting MatchingMap
   3.252 +    /// type.
   3.253 +    template <typename T>
   3.254 +    struct SetMatchingMap
   3.255 +      : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
   3.256 +      typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
   3.257 +    };
   3.258 +
   3.259 +    template <typename T>
   3.260 +    struct SetElevatorTraits : public Traits {
   3.261 +      typedef T Elevator;
   3.262 +      static Elevator *createElevator(const Graph&, int) {
   3.263 +        LEMON_ASSERT(false, "Elevator is not initialized");
   3.264 +        return 0; // ignore warnings
   3.265 +      }
   3.266 +    };
   3.267 +
   3.268 +    /// \brief \ref named-templ-param "Named parameter" for setting
   3.269 +    /// Elevator type
   3.270 +    ///
   3.271 +    /// \ref named-templ-param "Named parameter" for setting Elevator
   3.272 +    /// type. If this named parameter is used, then an external
   3.273 +    /// elevator object must be passed to the algorithm using the
   3.274 +    /// \ref elevator(Elevator&) "elevator()" function before calling
   3.275 +    /// \ref run() or \ref init().
   3.276 +    /// \sa SetStandardElevator
   3.277 +    template <typename T>
   3.278 +    struct SetElevator
   3.279 +      : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
   3.280 +      typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
   3.281 +    };
   3.282 +
   3.283 +    template <typename T>
   3.284 +    struct SetStandardElevatorTraits : public Traits {
   3.285 +      typedef T Elevator;
   3.286 +      static Elevator *createElevator(const Graph& graph, int max_level) {
   3.287 +        return new Elevator(graph, max_level);
   3.288 +      }
   3.289 +    };
   3.290 +
   3.291 +    /// \brief \ref named-templ-param "Named parameter" for setting
   3.292 +    /// Elevator type with automatic allocation
   3.293 +    ///
   3.294 +    /// \ref named-templ-param "Named parameter" for setting Elevator
   3.295 +    /// type with automatic allocation.
   3.296 +    /// The Elevator should have standard constructor interface to be
   3.297 +    /// able to automatically created by the algorithm (i.e. the
   3.298 +    /// graph and the maximum level should be passed to it).
   3.299 +    /// However an external elevator object could also be passed to the
   3.300 +    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
   3.301 +    /// before calling \ref run() or \ref init().
   3.302 +    /// \sa SetElevator
   3.303 +    template <typename T>
   3.304 +    struct SetStandardElevator
   3.305 +      : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
   3.306 +      typedef MaxFractionalMatching<Graph,
   3.307 +                                    SetStandardElevatorTraits<T> > Create;
   3.308 +    };
   3.309 +
   3.310 +    /// @}
   3.311 +
   3.312 +  protected:
   3.313 +
   3.314 +    MaxFractionalMatching() {}
   3.315 +
   3.316 +  public:
   3.317 +
   3.318 +    /// \brief Constructor
   3.319 +    ///
   3.320 +    /// Constructor.
   3.321 +    ///
   3.322 +    MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
   3.323 +      : _graph(graph), _allow_loops(allow_loops),
   3.324 +        _local_matching(false), _matching(0),
   3.325 +        _local_level(false), _level(0),  _indeg(0)
   3.326 +    {}
   3.327 +
   3.328 +    ~MaxFractionalMatching() {
   3.329 +      destroyStructures();
   3.330 +    }
   3.331 +
   3.332 +    /// \brief Sets the matching map.
   3.333 +    ///
   3.334 +    /// Sets the matching map.
   3.335 +    /// If you don't use this function before calling \ref run() or
   3.336 +    /// \ref init(), an instance will be allocated automatically.
   3.337 +    /// The destructor deallocates this automatically allocated map,
   3.338 +    /// of course.
   3.339 +    /// \return <tt>(*this)</tt>
   3.340 +    MaxFractionalMatching& matchingMap(MatchingMap& map) {
   3.341 +      if (_local_matching) {
   3.342 +        delete _matching;
   3.343 +        _local_matching = false;
   3.344 +      }
   3.345 +      _matching = &map;
   3.346 +      return *this;
   3.347 +    }
   3.348 +
   3.349 +    /// \brief Sets the elevator used by algorithm.
   3.350 +    ///
   3.351 +    /// Sets the elevator used by algorithm.
   3.352 +    /// If you don't use this function before calling \ref run() or
   3.353 +    /// \ref init(), an instance will be allocated automatically.
   3.354 +    /// The destructor deallocates this automatically allocated elevator,
   3.355 +    /// of course.
   3.356 +    /// \return <tt>(*this)</tt>
   3.357 +    MaxFractionalMatching& elevator(Elevator& elevator) {
   3.358 +      if (_local_level) {
   3.359 +        delete _level;
   3.360 +        _local_level = false;
   3.361 +      }
   3.362 +      _level = &elevator;
   3.363 +      return *this;
   3.364 +    }
   3.365 +
   3.366 +    /// \brief Returns a const reference to the elevator.
   3.367 +    ///
   3.368 +    /// Returns a const reference to the elevator.
   3.369 +    ///
   3.370 +    /// \pre Either \ref run() or \ref init() must be called before
   3.371 +    /// using this function.
   3.372 +    const Elevator& elevator() const {
   3.373 +      return *_level;
   3.374 +    }
   3.375 +
   3.376 +    /// \name Execution control
   3.377 +    /// The simplest way to execute the algorithm is to use one of the
   3.378 +    /// member functions called \c run(). \n
   3.379 +    /// If you need more control on the execution, first
   3.380 +    /// you must call \ref init() and then one variant of the start()
   3.381 +    /// member.
   3.382 +
   3.383 +    /// @{
   3.384 +
   3.385 +    /// \brief Initializes the internal data structures.
   3.386 +    ///
   3.387 +    /// Initializes the internal data structures and sets the initial
   3.388 +    /// matching.
   3.389 +    void init() {
   3.390 +      createStructures();
   3.391 +
   3.392 +      _level->initStart();
   3.393 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.394 +        _indeg->set(n, 0);
   3.395 +        _matching->set(n, INVALID);
   3.396 +        _level->initAddItem(n);
   3.397 +      }
   3.398 +      _level->initFinish();
   3.399 +
   3.400 +      _empty_level = _node_num;
   3.401 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.402 +        for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   3.403 +          if (_graph.target(a) == n && !_allow_loops) continue;
   3.404 +          _matching->set(n, a);
   3.405 +          Node v = _graph.target((*_matching)[n]);
   3.406 +          _indeg->set(v, (*_indeg)[v] + 1);
   3.407 +          break;
   3.408 +        }
   3.409 +      }
   3.410 +
   3.411 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.412 +        if ((*_indeg)[n] == 0) {
   3.413 +          _level->activate(n);
   3.414 +        }
   3.415 +      }
   3.416 +    }
   3.417 +
   3.418 +    /// \brief Starts the algorithm and computes a fractional matching
   3.419 +    ///
   3.420 +    /// The algorithm computes a maximum fractional matching.
   3.421 +    ///
   3.422 +    /// \param postprocess The algorithm computes first a matching
   3.423 +    /// which is a union of a matching with one value edges, cycles
   3.424 +    /// with half value edges and even length paths with half value
   3.425 +    /// edges. If the parameter is true, then after the push-relabel
   3.426 +    /// algorithm it postprocesses the matching to contain only
   3.427 +    /// matching edges and half value odd cycles.
   3.428 +    void start(bool postprocess = true) {
   3.429 +      Node n;
   3.430 +      while ((n = _level->highestActive()) != INVALID) {
   3.431 +        int level = _level->highestActiveLevel();
   3.432 +        int new_level = _level->maxLevel();
   3.433 +        for (InArcIt a(_graph, n); a != INVALID; ++a) {
   3.434 +          Node u = _graph.source(a);
   3.435 +          if (n == u && !_allow_loops) continue;
   3.436 +          Node v = _graph.target((*_matching)[u]);
   3.437 +          if ((*_level)[v] < level) {
   3.438 +            _indeg->set(v, (*_indeg)[v] - 1);
   3.439 +            if ((*_indeg)[v] == 0) {
   3.440 +              _level->activate(v);
   3.441 +            }
   3.442 +            _matching->set(u, a);
   3.443 +            _indeg->set(n, (*_indeg)[n] + 1);
   3.444 +            _level->deactivate(n);
   3.445 +            goto no_more_push;
   3.446 +          } else if (new_level > (*_level)[v]) {
   3.447 +            new_level = (*_level)[v];
   3.448 +          }
   3.449 +        }
   3.450 +
   3.451 +        if (new_level + 1 < _level->maxLevel()) {
   3.452 +          _level->liftHighestActive(new_level + 1);
   3.453 +        } else {
   3.454 +          _level->liftHighestActiveToTop();
   3.455 +        }
   3.456 +        if (_level->emptyLevel(level)) {
   3.457 +          _level->liftToTop(level);
   3.458 +        }
   3.459 +      no_more_push:
   3.460 +        ;
   3.461 +      }
   3.462 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.463 +        if ((*_matching)[n] == INVALID) continue;
   3.464 +        Node u = _graph.target((*_matching)[n]);
   3.465 +        if ((*_indeg)[u] > 1) {
   3.466 +          _indeg->set(u, (*_indeg)[u] - 1);
   3.467 +          _matching->set(n, INVALID);
   3.468 +        }
   3.469 +      }
   3.470 +      if (postprocess) {
   3.471 +        postprocessing();
   3.472 +      }
   3.473 +    }
   3.474 +
   3.475 +    /// \brief Starts the algorithm and computes a perfect fractional
   3.476 +    /// matching
   3.477 +    ///
   3.478 +    /// The algorithm computes a perfect fractional matching. If it
   3.479 +    /// does not exists, then the algorithm returns false and the
   3.480 +    /// matching is undefined and the barrier.
   3.481 +    ///
   3.482 +    /// \param postprocess The algorithm computes first a matching
   3.483 +    /// which is a union of a matching with one value edges, cycles
   3.484 +    /// with half value edges and even length paths with half value
   3.485 +    /// edges. If the parameter is true, then after the push-relabel
   3.486 +    /// algorithm it postprocesses the matching to contain only
   3.487 +    /// matching edges and half value odd cycles.
   3.488 +    bool startPerfect(bool postprocess = true) {
   3.489 +      Node n;
   3.490 +      while ((n = _level->highestActive()) != INVALID) {
   3.491 +        int level = _level->highestActiveLevel();
   3.492 +        int new_level = _level->maxLevel();
   3.493 +        for (InArcIt a(_graph, n); a != INVALID; ++a) {
   3.494 +          Node u = _graph.source(a);
   3.495 +          if (n == u && !_allow_loops) continue;
   3.496 +          Node v = _graph.target((*_matching)[u]);
   3.497 +          if ((*_level)[v] < level) {
   3.498 +            _indeg->set(v, (*_indeg)[v] - 1);
   3.499 +            if ((*_indeg)[v] == 0) {
   3.500 +              _level->activate(v);
   3.501 +            }
   3.502 +            _matching->set(u, a);
   3.503 +            _indeg->set(n, (*_indeg)[n] + 1);
   3.504 +            _level->deactivate(n);
   3.505 +            goto no_more_push;
   3.506 +          } else if (new_level > (*_level)[v]) {
   3.507 +            new_level = (*_level)[v];
   3.508 +          }
   3.509 +        }
   3.510 +
   3.511 +        if (new_level + 1 < _level->maxLevel()) {
   3.512 +          _level->liftHighestActive(new_level + 1);
   3.513 +        } else {
   3.514 +          _level->liftHighestActiveToTop();
   3.515 +          _empty_level = _level->maxLevel() - 1;
   3.516 +          return false;
   3.517 +        }
   3.518 +        if (_level->emptyLevel(level)) {
   3.519 +          _level->liftToTop(level);
   3.520 +          _empty_level = level;
   3.521 +          return false;
   3.522 +        }
   3.523 +      no_more_push:
   3.524 +        ;
   3.525 +      }
   3.526 +      if (postprocess) {
   3.527 +        postprocessing();
   3.528 +      }
   3.529 +      return true;
   3.530 +    }
   3.531 +
   3.532 +    /// \brief Runs the algorithm
   3.533 +    ///
   3.534 +    /// Just a shortcut for the next code:
   3.535 +    ///\code
   3.536 +    /// init();
   3.537 +    /// start();
   3.538 +    ///\endcode
   3.539 +    void run(bool postprocess = true) {
   3.540 +      init();
   3.541 +      start(postprocess);
   3.542 +    }
   3.543 +
   3.544 +    /// \brief Runs the algorithm to find a perfect fractional matching
   3.545 +    ///
   3.546 +    /// Just a shortcut for the next code:
   3.547 +    ///\code
   3.548 +    /// init();
   3.549 +    /// startPerfect();
   3.550 +    ///\endcode
   3.551 +    bool runPerfect(bool postprocess = true) {
   3.552 +      init();
   3.553 +      return startPerfect(postprocess);
   3.554 +    }
   3.555 +
   3.556 +    ///@}
   3.557 +
   3.558 +    /// \name Query Functions
   3.559 +    /// The result of the %Matching algorithm can be obtained using these
   3.560 +    /// functions.\n
   3.561 +    /// Before the use of these functions,
   3.562 +    /// either run() or start() must be called.
   3.563 +    ///@{
   3.564 +
   3.565 +
   3.566 +    /// \brief Return the number of covered nodes in the matching.
   3.567 +    ///
   3.568 +    /// This function returns the number of covered nodes in the matching.
   3.569 +    ///
   3.570 +    /// \pre Either run() or start() must be called before using this function.
   3.571 +    int matchingSize() const {
   3.572 +      int num = 0;
   3.573 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.574 +        if ((*_matching)[n] != INVALID) {
   3.575 +          ++num;
   3.576 +        }
   3.577 +      }
   3.578 +      return num;
   3.579 +    }
   3.580 +
   3.581 +    /// \brief Returns a const reference to the matching map.
   3.582 +    ///
   3.583 +    /// Returns a const reference to the node map storing the found
   3.584 +    /// fractional matching. This method can be called after
   3.585 +    /// running the algorithm.
   3.586 +    ///
   3.587 +    /// \pre Either \ref run() or \ref init() must be called before
   3.588 +    /// using this function.
   3.589 +    const MatchingMap& matchingMap() const {
   3.590 +      return *_matching;
   3.591 +    }
   3.592 +
   3.593 +    /// \brief Return \c true if the given edge is in the matching.
   3.594 +    ///
   3.595 +    /// This function returns \c true if the given edge is in the
   3.596 +    /// found matching. The result is scaled by \ref primalScale
   3.597 +    /// "primal scale".
   3.598 +    ///
   3.599 +    /// \pre Either run() or start() must be called before using this function.
   3.600 +    int matching(const Edge& edge) const {
   3.601 +      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
   3.602 +        (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
   3.603 +    }
   3.604 +
   3.605 +    /// \brief Return the fractional matching arc (or edge) incident
   3.606 +    /// to the given node.
   3.607 +    ///
   3.608 +    /// This function returns one of the fractional matching arc (or
   3.609 +    /// edge) incident to the given node in the found matching or \c
   3.610 +    /// INVALID if the node is not covered by the matching or if the
   3.611 +    /// node is on an odd length cycle then it is the successor edge
   3.612 +    /// on the cycle.
   3.613 +    ///
   3.614 +    /// \pre Either run() or start() must be called before using this function.
   3.615 +    Arc matching(const Node& node) const {
   3.616 +      return (*_matching)[node];
   3.617 +    }
   3.618 +
   3.619 +    /// \brief Returns true if the node is in the barrier
   3.620 +    ///
   3.621 +    /// The barrier is a subset of the nodes. If the nodes in the
   3.622 +    /// barrier have less adjacent nodes than the size of the barrier,
   3.623 +    /// then at least as much nodes cannot be covered as the
   3.624 +    /// difference of the two subsets.
   3.625 +    bool barrier(const Node& node) const {
   3.626 +      return (*_level)[node] >= _empty_level;
   3.627 +    }
   3.628 +
   3.629 +    /// @}
   3.630 +
   3.631 +  };
   3.632 +
   3.633 +  /// \ingroup matching
   3.634 +  ///
   3.635 +  /// \brief Weighted fractional matching in general graphs
   3.636 +  ///
   3.637 +  /// This class provides an efficient implementation of fractional
   3.638 +  /// matching algorithm. The implementation uses priority queues and
   3.639 +  /// provides \f$O(nm\log n)\f$ time complexity.
   3.640 +  ///
   3.641 +  /// The maximum weighted fractional matching is a relaxation of the
   3.642 +  /// maximum weighted matching problem where the odd set constraints
   3.643 +  /// are omitted.
   3.644 +  /// It can be formulated with the following linear program.
   3.645 +  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   3.646 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
   3.647 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
   3.648 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   3.649 +  /// \f$X\f$. The result must be the union of a matching with one
   3.650 +  /// value edges and a set of odd length cycles with half value edges.
   3.651 +  ///
   3.652 +  /// The algorithm calculates an optimal fractional matching and a
   3.653 +  /// proof of the optimality. The solution of the dual problem can be
   3.654 +  /// used to check the result of the algorithm. The dual linear
   3.655 +  /// problem is the following.
   3.656 +  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
   3.657 +  /// \f[y_u \ge 0 \quad \forall u \in V\f]
   3.658 +  /// \f[\min \sum_{u \in V}y_u \f]
   3.659 +  ///
   3.660 +  /// The algorithm can be executed with the run() function.
   3.661 +  /// After it the matching (the primal solution) and the dual solution
   3.662 +  /// can be obtained using the query functions.
   3.663 +  ///
   3.664 +  /// The primal solution is multiplied by
   3.665 +  /// \ref MaxWeightedFractionalMatching::primalScale "2".
   3.666 +  /// If the value type is integer, then the dual
   3.667 +  /// solution is scaled by
   3.668 +  /// \ref MaxWeightedFractionalMatching::dualScale "4".
   3.669 +  ///
   3.670 +  /// \tparam GR The undirected graph type the algorithm runs on.
   3.671 +  /// \tparam WM The type edge weight map. The default type is
   3.672 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   3.673 +#ifdef DOXYGEN
   3.674 +  template <typename GR, typename WM>
   3.675 +#else
   3.676 +  template <typename GR,
   3.677 +            typename WM = typename GR::template EdgeMap<int> >
   3.678 +#endif
   3.679 +  class MaxWeightedFractionalMatching {
   3.680 +  public:
   3.681 +
   3.682 +    /// The graph type of the algorithm
   3.683 +    typedef GR Graph;
   3.684 +    /// The type of the edge weight map
   3.685 +    typedef WM WeightMap;
   3.686 +    /// The value type of the edge weights
   3.687 +    typedef typename WeightMap::Value Value;
   3.688 +
   3.689 +    /// The type of the matching map
   3.690 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
   3.691 +    MatchingMap;
   3.692 +
   3.693 +    /// \brief Scaling factor for primal solution
   3.694 +    ///
   3.695 +    /// Scaling factor for primal solution.
   3.696 +    static const int primalScale = 2;
   3.697 +
   3.698 +    /// \brief Scaling factor for dual solution
   3.699 +    ///
   3.700 +    /// Scaling factor for dual solution. It is equal to 4 or 1
   3.701 +    /// according to the value type.
   3.702 +    static const int dualScale =
   3.703 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
   3.704 +
   3.705 +  private:
   3.706 +
   3.707 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
   3.708 +
   3.709 +    typedef typename Graph::template NodeMap<Value> NodePotential;
   3.710 +
   3.711 +    const Graph& _graph;
   3.712 +    const WeightMap& _weight;
   3.713 +
   3.714 +    MatchingMap* _matching;
   3.715 +    NodePotential* _node_potential;
   3.716 +
   3.717 +    int _node_num;
   3.718 +    bool _allow_loops;
   3.719 +
   3.720 +    enum Status {
   3.721 +      EVEN = -1, MATCHED = 0, ODD = 1
   3.722 +    };
   3.723 +
   3.724 +    typedef typename Graph::template NodeMap<Status> StatusMap;
   3.725 +    StatusMap* _status;
   3.726 +
   3.727 +    typedef typename Graph::template NodeMap<Arc> PredMap;
   3.728 +    PredMap* _pred;
   3.729 +
   3.730 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
   3.731 +
   3.732 +    IntNodeMap *_tree_set_index;
   3.733 +    TreeSet *_tree_set;
   3.734 +
   3.735 +    IntNodeMap *_delta1_index;
   3.736 +    BinHeap<Value, IntNodeMap> *_delta1;
   3.737 +
   3.738 +    IntNodeMap *_delta2_index;
   3.739 +    BinHeap<Value, IntNodeMap> *_delta2;
   3.740 +
   3.741 +    IntEdgeMap *_delta3_index;
   3.742 +    BinHeap<Value, IntEdgeMap> *_delta3;
   3.743 +
   3.744 +    Value _delta_sum;
   3.745 +
   3.746 +    void createStructures() {
   3.747 +      _node_num = countNodes(_graph);
   3.748 +
   3.749 +      if (!_matching) {
   3.750 +        _matching = new MatchingMap(_graph);
   3.751 +      }
   3.752 +      if (!_node_potential) {
   3.753 +        _node_potential = new NodePotential(_graph);
   3.754 +      }
   3.755 +      if (!_status) {
   3.756 +        _status = new StatusMap(_graph);
   3.757 +      }
   3.758 +      if (!_pred) {
   3.759 +        _pred = new PredMap(_graph);
   3.760 +      }
   3.761 +      if (!_tree_set) {
   3.762 +        _tree_set_index = new IntNodeMap(_graph);
   3.763 +        _tree_set = new TreeSet(*_tree_set_index);
   3.764 +      }
   3.765 +      if (!_delta1) {
   3.766 +        _delta1_index = new IntNodeMap(_graph);
   3.767 +        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   3.768 +      }
   3.769 +      if (!_delta2) {
   3.770 +        _delta2_index = new IntNodeMap(_graph);
   3.771 +        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
   3.772 +      }
   3.773 +      if (!_delta3) {
   3.774 +        _delta3_index = new IntEdgeMap(_graph);
   3.775 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   3.776 +      }
   3.777 +    }
   3.778 +
   3.779 +    void destroyStructures() {
   3.780 +      if (_matching) {
   3.781 +        delete _matching;
   3.782 +      }
   3.783 +      if (_node_potential) {
   3.784 +        delete _node_potential;
   3.785 +      }
   3.786 +      if (_status) {
   3.787 +        delete _status;
   3.788 +      }
   3.789 +      if (_pred) {
   3.790 +        delete _pred;
   3.791 +      }
   3.792 +      if (_tree_set) {
   3.793 +        delete _tree_set_index;
   3.794 +        delete _tree_set;
   3.795 +      }
   3.796 +      if (_delta1) {
   3.797 +        delete _delta1_index;
   3.798 +        delete _delta1;
   3.799 +      }
   3.800 +      if (_delta2) {
   3.801 +        delete _delta2_index;
   3.802 +        delete _delta2;
   3.803 +      }
   3.804 +      if (_delta3) {
   3.805 +        delete _delta3_index;
   3.806 +        delete _delta3;
   3.807 +      }
   3.808 +    }
   3.809 +
   3.810 +    void matchedToEven(Node node, int tree) {
   3.811 +      _tree_set->insert(node, tree);
   3.812 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   3.813 +      _delta1->push(node, (*_node_potential)[node]);
   3.814 +
   3.815 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
   3.816 +        _delta2->erase(node);
   3.817 +      }
   3.818 +
   3.819 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
   3.820 +        Node v = _graph.source(a);
   3.821 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   3.822 +          dualScale * _weight[a];
   3.823 +        if (node == v) {
   3.824 +          if (_allow_loops && _graph.direction(a)) {
   3.825 +            _delta3->push(a, rw / 2);
   3.826 +          }
   3.827 +        } else if ((*_status)[v] == EVEN) {
   3.828 +          _delta3->push(a, rw / 2);
   3.829 +        } else if ((*_status)[v] == MATCHED) {
   3.830 +          if (_delta2->state(v) != _delta2->IN_HEAP) {
   3.831 +            _pred->set(v, a);
   3.832 +            _delta2->push(v, rw);
   3.833 +          } else if ((*_delta2)[v] > rw) {
   3.834 +            _pred->set(v, a);
   3.835 +            _delta2->decrease(v, rw);
   3.836 +          }
   3.837 +        }
   3.838 +      }
   3.839 +    }
   3.840 +
   3.841 +    void matchedToOdd(Node node, int tree) {
   3.842 +      _tree_set->insert(node, tree);
   3.843 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   3.844 +
   3.845 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
   3.846 +        _delta2->erase(node);
   3.847 +      }
   3.848 +    }
   3.849 +
   3.850 +    void evenToMatched(Node node, int tree) {
   3.851 +      _delta1->erase(node);
   3.852 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   3.853 +      Arc min = INVALID;
   3.854 +      Value minrw = std::numeric_limits<Value>::max();
   3.855 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
   3.856 +        Node v = _graph.source(a);
   3.857 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   3.858 +          dualScale * _weight[a];
   3.859 +
   3.860 +        if (node == v) {
   3.861 +          if (_allow_loops && _graph.direction(a)) {
   3.862 +            _delta3->erase(a);
   3.863 +          }
   3.864 +        } else if ((*_status)[v] == EVEN) {
   3.865 +          _delta3->erase(a);
   3.866 +          if (minrw > rw) {
   3.867 +            min = _graph.oppositeArc(a);
   3.868 +            minrw = rw;
   3.869 +          }
   3.870 +        } else if ((*_status)[v]  == MATCHED) {
   3.871 +          if ((*_pred)[v] == a) {
   3.872 +            Arc mina = INVALID;
   3.873 +            Value minrwa = std::numeric_limits<Value>::max();
   3.874 +            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
   3.875 +              Node va = _graph.target(aa);
   3.876 +              if ((*_status)[va] != EVEN ||
   3.877 +                  _tree_set->find(va) == tree) continue;
   3.878 +              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
   3.879 +                dualScale * _weight[aa];
   3.880 +              if (minrwa > rwa) {
   3.881 +                minrwa = rwa;
   3.882 +                mina = aa;
   3.883 +              }
   3.884 +            }
   3.885 +            if (mina != INVALID) {
   3.886 +              _pred->set(v, mina);
   3.887 +              _delta2->increase(v, minrwa);
   3.888 +            } else {
   3.889 +              _pred->set(v, INVALID);
   3.890 +              _delta2->erase(v);
   3.891 +            }
   3.892 +          }
   3.893 +        }
   3.894 +      }
   3.895 +      if (min != INVALID) {
   3.896 +        _pred->set(node, min);
   3.897 +        _delta2->push(node, minrw);
   3.898 +      } else {
   3.899 +        _pred->set(node, INVALID);
   3.900 +      }
   3.901 +    }
   3.902 +
   3.903 +    void oddToMatched(Node node) {
   3.904 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   3.905 +      Arc min = INVALID;
   3.906 +      Value minrw = std::numeric_limits<Value>::max();
   3.907 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
   3.908 +        Node v = _graph.source(a);
   3.909 +        if ((*_status)[v] != EVEN) continue;
   3.910 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   3.911 +          dualScale * _weight[a];
   3.912 +
   3.913 +        if (minrw > rw) {
   3.914 +          min = _graph.oppositeArc(a);
   3.915 +          minrw = rw;
   3.916 +        }
   3.917 +      }
   3.918 +      if (min != INVALID) {
   3.919 +        _pred->set(node, min);
   3.920 +        _delta2->push(node, minrw);
   3.921 +      } else {
   3.922 +        _pred->set(node, INVALID);
   3.923 +      }
   3.924 +    }
   3.925 +
   3.926 +    void alternatePath(Node even, int tree) {
   3.927 +      Node odd;
   3.928 +
   3.929 +      _status->set(even, MATCHED);
   3.930 +      evenToMatched(even, tree);
   3.931 +
   3.932 +      Arc prev = (*_matching)[even];
   3.933 +      while (prev != INVALID) {
   3.934 +        odd = _graph.target(prev);
   3.935 +        even = _graph.target((*_pred)[odd]);
   3.936 +        _matching->set(odd, (*_pred)[odd]);
   3.937 +        _status->set(odd, MATCHED);
   3.938 +        oddToMatched(odd);
   3.939 +
   3.940 +        prev = (*_matching)[even];
   3.941 +        _status->set(even, MATCHED);
   3.942 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
   3.943 +        evenToMatched(even, tree);
   3.944 +      }
   3.945 +    }
   3.946 +
   3.947 +    void destroyTree(int tree) {
   3.948 +      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
   3.949 +        if ((*_status)[n] == EVEN) {
   3.950 +          _status->set(n, MATCHED);
   3.951 +          evenToMatched(n, tree);
   3.952 +        } else if ((*_status)[n] == ODD) {
   3.953 +          _status->set(n, MATCHED);
   3.954 +          oddToMatched(n);
   3.955 +        }
   3.956 +      }
   3.957 +      _tree_set->eraseClass(tree);
   3.958 +    }
   3.959 +
   3.960 +
   3.961 +    void unmatchNode(const Node& node) {
   3.962 +      int tree = _tree_set->find(node);
   3.963 +
   3.964 +      alternatePath(node, tree);
   3.965 +      destroyTree(tree);
   3.966 +
   3.967 +      _matching->set(node, INVALID);
   3.968 +    }
   3.969 +
   3.970 +
   3.971 +    void augmentOnEdge(const Edge& edge) {
   3.972 +      Node left = _graph.u(edge);
   3.973 +      int left_tree = _tree_set->find(left);
   3.974 +
   3.975 +      alternatePath(left, left_tree);
   3.976 +      destroyTree(left_tree);
   3.977 +      _matching->set(left, _graph.direct(edge, true));
   3.978 +
   3.979 +      Node right = _graph.v(edge);
   3.980 +      int right_tree = _tree_set->find(right);
   3.981 +
   3.982 +      alternatePath(right, right_tree);
   3.983 +      destroyTree(right_tree);
   3.984 +      _matching->set(right, _graph.direct(edge, false));
   3.985 +    }
   3.986 +
   3.987 +    void augmentOnArc(const Arc& arc) {
   3.988 +      Node left = _graph.source(arc);
   3.989 +      _status->set(left, MATCHED);
   3.990 +      _matching->set(left, arc);
   3.991 +      _pred->set(left, arc);
   3.992 +
   3.993 +      Node right = _graph.target(arc);
   3.994 +      int right_tree = _tree_set->find(right);
   3.995 +
   3.996 +      alternatePath(right, right_tree);
   3.997 +      destroyTree(right_tree);
   3.998 +      _matching->set(right, _graph.oppositeArc(arc));
   3.999 +    }
  3.1000 +
  3.1001 +    void extendOnArc(const Arc& arc) {
  3.1002 +      Node base = _graph.target(arc);
  3.1003 +      int tree = _tree_set->find(base);
  3.1004 +
  3.1005 +      Node odd = _graph.source(arc);
  3.1006 +      _tree_set->insert(odd, tree);
  3.1007 +      _status->set(odd, ODD);
  3.1008 +      matchedToOdd(odd, tree);
  3.1009 +      _pred->set(odd, arc);
  3.1010 +
  3.1011 +      Node even = _graph.target((*_matching)[odd]);
  3.1012 +      _tree_set->insert(even, tree);
  3.1013 +      _status->set(even, EVEN);
  3.1014 +      matchedToEven(even, tree);
  3.1015 +    }
  3.1016 +
  3.1017 +    void cycleOnEdge(const Edge& edge, int tree) {
  3.1018 +      Node nca = INVALID;
  3.1019 +      std::vector<Node> left_path, right_path;
  3.1020 +
  3.1021 +      {
  3.1022 +        std::set<Node> left_set, right_set;
  3.1023 +        Node left = _graph.u(edge);
  3.1024 +        left_path.push_back(left);
  3.1025 +        left_set.insert(left);
  3.1026 +
  3.1027 +        Node right = _graph.v(edge);
  3.1028 +        right_path.push_back(right);
  3.1029 +        right_set.insert(right);
  3.1030 +
  3.1031 +        while (true) {
  3.1032 +
  3.1033 +          if (left_set.find(right) != left_set.end()) {
  3.1034 +            nca = right;
  3.1035 +            break;
  3.1036 +          }
  3.1037 +
  3.1038 +          if ((*_matching)[left] == INVALID) break;
  3.1039 +
  3.1040 +          left = _graph.target((*_matching)[left]);
  3.1041 +          left_path.push_back(left);
  3.1042 +          left = _graph.target((*_pred)[left]);
  3.1043 +          left_path.push_back(left);
  3.1044 +
  3.1045 +          left_set.insert(left);
  3.1046 +
  3.1047 +          if (right_set.find(left) != right_set.end()) {
  3.1048 +            nca = left;
  3.1049 +            break;
  3.1050 +          }
  3.1051 +
  3.1052 +          if ((*_matching)[right] == INVALID) break;
  3.1053 +
  3.1054 +          right = _graph.target((*_matching)[right]);
  3.1055 +          right_path.push_back(right);
  3.1056 +          right = _graph.target((*_pred)[right]);
  3.1057 +          right_path.push_back(right);
  3.1058 +
  3.1059 +          right_set.insert(right);
  3.1060 +
  3.1061 +        }
  3.1062 +
  3.1063 +        if (nca == INVALID) {
  3.1064 +          if ((*_matching)[left] == INVALID) {
  3.1065 +            nca = right;
  3.1066 +            while (left_set.find(nca) == left_set.end()) {
  3.1067 +              nca = _graph.target((*_matching)[nca]);
  3.1068 +              right_path.push_back(nca);
  3.1069 +              nca = _graph.target((*_pred)[nca]);
  3.1070 +              right_path.push_back(nca);
  3.1071 +            }
  3.1072 +          } else {
  3.1073 +            nca = left;
  3.1074 +            while (right_set.find(nca) == right_set.end()) {
  3.1075 +              nca = _graph.target((*_matching)[nca]);
  3.1076 +              left_path.push_back(nca);
  3.1077 +              nca = _graph.target((*_pred)[nca]);
  3.1078 +              left_path.push_back(nca);
  3.1079 +            }
  3.1080 +          }
  3.1081 +        }
  3.1082 +      }
  3.1083 +
  3.1084 +      alternatePath(nca, tree);
  3.1085 +      Arc prev;
  3.1086 +
  3.1087 +      prev = _graph.direct(edge, true);
  3.1088 +      for (int i = 0; left_path[i] != nca; i += 2) {
  3.1089 +        _matching->set(left_path[i], prev);
  3.1090 +        _status->set(left_path[i], MATCHED);
  3.1091 +        evenToMatched(left_path[i], tree);
  3.1092 +
  3.1093 +        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  3.1094 +        _status->set(left_path[i + 1], MATCHED);
  3.1095 +        oddToMatched(left_path[i + 1]);
  3.1096 +      }
  3.1097 +      _matching->set(nca, prev);
  3.1098 +
  3.1099 +      for (int i = 0; right_path[i] != nca; i += 2) {
  3.1100 +        _status->set(right_path[i], MATCHED);
  3.1101 +        evenToMatched(right_path[i], tree);
  3.1102 +
  3.1103 +        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  3.1104 +        _status->set(right_path[i + 1], MATCHED);
  3.1105 +        oddToMatched(right_path[i + 1]);
  3.1106 +      }
  3.1107 +
  3.1108 +      destroyTree(tree);
  3.1109 +    }
  3.1110 +
  3.1111 +    void extractCycle(const Arc &arc) {
  3.1112 +      Node left = _graph.source(arc);
  3.1113 +      Node odd = _graph.target((*_matching)[left]);
  3.1114 +      Arc prev;
  3.1115 +      while (odd != left) {
  3.1116 +        Node even = _graph.target((*_matching)[odd]);
  3.1117 +        prev = (*_matching)[odd];
  3.1118 +        odd = _graph.target((*_matching)[even]);
  3.1119 +        _matching->set(even, _graph.oppositeArc(prev));
  3.1120 +      }
  3.1121 +      _matching->set(left, arc);
  3.1122 +
  3.1123 +      Node right = _graph.target(arc);
  3.1124 +      int right_tree = _tree_set->find(right);
  3.1125 +      alternatePath(right, right_tree);
  3.1126 +      destroyTree(right_tree);
  3.1127 +      _matching->set(right, _graph.oppositeArc(arc));
  3.1128 +    }
  3.1129 +
  3.1130 +  public:
  3.1131 +
  3.1132 +    /// \brief Constructor
  3.1133 +    ///
  3.1134 +    /// Constructor.
  3.1135 +    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
  3.1136 +                                  bool allow_loops = true)
  3.1137 +      : _graph(graph), _weight(weight), _matching(0),
  3.1138 +      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  3.1139 +      _status(0),  _pred(0),
  3.1140 +      _tree_set_index(0), _tree_set(0),
  3.1141 +
  3.1142 +      _delta1_index(0), _delta1(0),
  3.1143 +      _delta2_index(0), _delta2(0),
  3.1144 +      _delta3_index(0), _delta3(0),
  3.1145 +
  3.1146 +      _delta_sum() {}
  3.1147 +
  3.1148 +    ~MaxWeightedFractionalMatching() {
  3.1149 +      destroyStructures();
  3.1150 +    }
  3.1151 +
  3.1152 +    /// \name Execution Control
  3.1153 +    /// The simplest way to execute the algorithm is to use the
  3.1154 +    /// \ref run() member function.
  3.1155 +
  3.1156 +    ///@{
  3.1157 +
  3.1158 +    /// \brief Initialize the algorithm
  3.1159 +    ///
  3.1160 +    /// This function initializes the algorithm.
  3.1161 +    void init() {
  3.1162 +      createStructures();
  3.1163 +
  3.1164 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1165 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
  3.1166 +        (*_delta2_index)[n] = _delta2->PRE_HEAP;
  3.1167 +      }
  3.1168 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  3.1169 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  3.1170 +      }
  3.1171 +
  3.1172 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1173 +        Value max = 0;
  3.1174 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  3.1175 +          if (_graph.target(e) == n && !_allow_loops) continue;
  3.1176 +          if ((dualScale * _weight[e]) / 2 > max) {
  3.1177 +            max = (dualScale * _weight[e]) / 2;
  3.1178 +          }
  3.1179 +        }
  3.1180 +        _node_potential->set(n, max);
  3.1181 +        _delta1->push(n, max);
  3.1182 +
  3.1183 +        _tree_set->insert(n);
  3.1184 +
  3.1185 +        _matching->set(n, INVALID);
  3.1186 +        _status->set(n, EVEN);
  3.1187 +      }
  3.1188 +
  3.1189 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  3.1190 +        Node left = _graph.u(e);
  3.1191 +        Node right = _graph.v(e);
  3.1192 +        if (left == right && !_allow_loops) continue;
  3.1193 +        _delta3->push(e, ((*_node_potential)[left] +
  3.1194 +                          (*_node_potential)[right] -
  3.1195 +                          dualScale * _weight[e]) / 2);
  3.1196 +      }
  3.1197 +    }
  3.1198 +
  3.1199 +    /// \brief Start the algorithm
  3.1200 +    ///
  3.1201 +    /// This function starts the algorithm.
  3.1202 +    ///
  3.1203 +    /// \pre \ref init() must be called before using this function.
  3.1204 +    void start() {
  3.1205 +      enum OpType {
  3.1206 +        D1, D2, D3
  3.1207 +      };
  3.1208 +
  3.1209 +      int unmatched = _node_num;
  3.1210 +      while (unmatched > 0) {
  3.1211 +        Value d1 = !_delta1->empty() ?
  3.1212 +          _delta1->prio() : std::numeric_limits<Value>::max();
  3.1213 +
  3.1214 +        Value d2 = !_delta2->empty() ?
  3.1215 +          _delta2->prio() : std::numeric_limits<Value>::max();
  3.1216 +
  3.1217 +        Value d3 = !_delta3->empty() ?
  3.1218 +          _delta3->prio() : std::numeric_limits<Value>::max();
  3.1219 +
  3.1220 +        _delta_sum = d3; OpType ot = D3;
  3.1221 +        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
  3.1222 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  3.1223 +
  3.1224 +        switch (ot) {
  3.1225 +        case D1:
  3.1226 +          {
  3.1227 +            Node n = _delta1->top();
  3.1228 +            unmatchNode(n);
  3.1229 +            --unmatched;
  3.1230 +          }
  3.1231 +          break;
  3.1232 +        case D2:
  3.1233 +          {
  3.1234 +            Node n = _delta2->top();
  3.1235 +            Arc a = (*_pred)[n];
  3.1236 +            if ((*_matching)[n] == INVALID) {
  3.1237 +              augmentOnArc(a);
  3.1238 +              --unmatched;
  3.1239 +            } else {
  3.1240 +              Node v = _graph.target((*_matching)[n]);
  3.1241 +              if ((*_matching)[n] !=
  3.1242 +                  _graph.oppositeArc((*_matching)[v])) {
  3.1243 +                extractCycle(a);
  3.1244 +                --unmatched;
  3.1245 +              } else {
  3.1246 +                extendOnArc(a);
  3.1247 +              }
  3.1248 +            }
  3.1249 +          } break;
  3.1250 +        case D3:
  3.1251 +          {
  3.1252 +            Edge e = _delta3->top();
  3.1253 +
  3.1254 +            Node left = _graph.u(e);
  3.1255 +            Node right = _graph.v(e);
  3.1256 +
  3.1257 +            int left_tree = _tree_set->find(left);
  3.1258 +            int right_tree = _tree_set->find(right);
  3.1259 +
  3.1260 +            if (left_tree == right_tree) {
  3.1261 +              cycleOnEdge(e, left_tree);
  3.1262 +              --unmatched;
  3.1263 +            } else {
  3.1264 +              augmentOnEdge(e);
  3.1265 +              unmatched -= 2;
  3.1266 +            }
  3.1267 +          } break;
  3.1268 +        }
  3.1269 +      }
  3.1270 +    }
  3.1271 +
  3.1272 +    /// \brief Run the algorithm.
  3.1273 +    ///
  3.1274 +    /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
  3.1275 +    ///
  3.1276 +    /// \note mwfm.run() is just a shortcut of the following code.
  3.1277 +    /// \code
  3.1278 +    ///   mwfm.init();
  3.1279 +    ///   mwfm.start();
  3.1280 +    /// \endcode
  3.1281 +    void run() {
  3.1282 +      init();
  3.1283 +      start();
  3.1284 +    }
  3.1285 +
  3.1286 +    /// @}
  3.1287 +
  3.1288 +    /// \name Primal Solution
  3.1289 +    /// Functions to get the primal solution, i.e. the maximum weighted
  3.1290 +    /// matching.\n
  3.1291 +    /// Either \ref run() or \ref start() function should be called before
  3.1292 +    /// using them.
  3.1293 +
  3.1294 +    /// @{
  3.1295 +
  3.1296 +    /// \brief Return the weight of the matching.
  3.1297 +    ///
  3.1298 +    /// This function returns the weight of the found matching. This
  3.1299 +    /// value is scaled by \ref primalScale "primal scale".
  3.1300 +    ///
  3.1301 +    /// \pre Either run() or start() must be called before using this function.
  3.1302 +    Value matchingWeight() const {
  3.1303 +      Value sum = 0;
  3.1304 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1305 +        if ((*_matching)[n] != INVALID) {
  3.1306 +          sum += _weight[(*_matching)[n]];
  3.1307 +        }
  3.1308 +      }
  3.1309 +      return sum * primalScale / 2;
  3.1310 +    }
  3.1311 +
  3.1312 +    /// \brief Return the number of covered nodes in the matching.
  3.1313 +    ///
  3.1314 +    /// This function returns the number of covered nodes in the matching.
  3.1315 +    ///
  3.1316 +    /// \pre Either run() or start() must be called before using this function.
  3.1317 +    int matchingSize() const {
  3.1318 +      int num = 0;
  3.1319 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1320 +        if ((*_matching)[n] != INVALID) {
  3.1321 +          ++num;
  3.1322 +        }
  3.1323 +      }
  3.1324 +      return num;
  3.1325 +    }
  3.1326 +
  3.1327 +    /// \brief Return \c true if the given edge is in the matching.
  3.1328 +    ///
  3.1329 +    /// This function returns \c true if the given edge is in the
  3.1330 +    /// found matching. The result is scaled by \ref primalScale
  3.1331 +    /// "primal scale".
  3.1332 +    ///
  3.1333 +    /// \pre Either run() or start() must be called before using this function.
  3.1334 +    int matching(const Edge& edge) const {
  3.1335 +      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  3.1336 +        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
  3.1337 +    }
  3.1338 +
  3.1339 +    /// \brief Return the fractional matching arc (or edge) incident
  3.1340 +    /// to the given node.
  3.1341 +    ///
  3.1342 +    /// This function returns one of the fractional matching arc (or
  3.1343 +    /// edge) incident to the given node in the found matching or \c
  3.1344 +    /// INVALID if the node is not covered by the matching or if the
  3.1345 +    /// node is on an odd length cycle then it is the successor edge
  3.1346 +    /// on the cycle.
  3.1347 +    ///
  3.1348 +    /// \pre Either run() or start() must be called before using this function.
  3.1349 +    Arc matching(const Node& node) const {
  3.1350 +      return (*_matching)[node];
  3.1351 +    }
  3.1352 +
  3.1353 +    /// \brief Return a const reference to the matching map.
  3.1354 +    ///
  3.1355 +    /// This function returns a const reference to a node map that stores
  3.1356 +    /// the matching arc (or edge) incident to each node.
  3.1357 +    const MatchingMap& matchingMap() const {
  3.1358 +      return *_matching;
  3.1359 +    }
  3.1360 +
  3.1361 +    /// @}
  3.1362 +
  3.1363 +    /// \name Dual Solution
  3.1364 +    /// Functions to get the dual solution.\n
  3.1365 +    /// Either \ref run() or \ref start() function should be called before
  3.1366 +    /// using them.
  3.1367 +
  3.1368 +    /// @{
  3.1369 +
  3.1370 +    /// \brief Return the value of the dual solution.
  3.1371 +    ///
  3.1372 +    /// This function returns the value of the dual solution.
  3.1373 +    /// It should be equal to the primal value scaled by \ref dualScale
  3.1374 +    /// "dual scale".
  3.1375 +    ///
  3.1376 +    /// \pre Either run() or start() must be called before using this function.
  3.1377 +    Value dualValue() const {
  3.1378 +      Value sum = 0;
  3.1379 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1380 +        sum += nodeValue(n);
  3.1381 +      }
  3.1382 +      return sum;
  3.1383 +    }
  3.1384 +
  3.1385 +    /// \brief Return the dual value (potential) of the given node.
  3.1386 +    ///
  3.1387 +    /// This function returns the dual value (potential) of the given node.
  3.1388 +    ///
  3.1389 +    /// \pre Either run() or start() must be called before using this function.
  3.1390 +    Value nodeValue(const Node& n) const {
  3.1391 +      return (*_node_potential)[n];
  3.1392 +    }
  3.1393 +
  3.1394 +    /// @}
  3.1395 +
  3.1396 +  };
  3.1397 +
  3.1398 +  /// \ingroup matching
  3.1399 +  ///
  3.1400 +  /// \brief Weighted fractional perfect matching in general graphs
  3.1401 +  ///
  3.1402 +  /// This class provides an efficient implementation of fractional
  3.1403 +  /// matching algorithm. The implementation uses priority queues and
  3.1404 +  /// provides \f$O(nm\log n)\f$ time complexity.
  3.1405 +  ///
  3.1406 +  /// The maximum weighted fractional perfect matching is a relaxation
  3.1407 +  /// of the maximum weighted perfect matching problem where the odd
  3.1408 +  /// set constraints are omitted.
  3.1409 +  /// It can be formulated with the following linear program.
  3.1410 +  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  3.1411 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
  3.1412 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
  3.1413 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  3.1414 +  /// \f$X\f$. The result must be the union of a matching with one
  3.1415 +  /// value edges and a set of odd length cycles with half value edges.
  3.1416 +  ///
  3.1417 +  /// The algorithm calculates an optimal fractional matching and a
  3.1418 +  /// proof of the optimality. The solution of the dual problem can be
  3.1419 +  /// used to check the result of the algorithm. The dual linear
  3.1420 +  /// problem is the following.
  3.1421 +  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
  3.1422 +  /// \f[\min \sum_{u \in V}y_u \f]
  3.1423 +  ///
  3.1424 +  /// The algorithm can be executed with the run() function.
  3.1425 +  /// After it the matching (the primal solution) and the dual solution
  3.1426 +  /// can be obtained using the query functions.
  3.1427 +  ///
  3.1428 +  /// The primal solution is multiplied by
  3.1429 +  /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
  3.1430 +  /// If the value type is integer, then the dual
  3.1431 +  /// solution is scaled by
  3.1432 +  /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
  3.1433 +  ///
  3.1434 +  /// \tparam GR The undirected graph type the algorithm runs on.
  3.1435 +  /// \tparam WM The type edge weight map. The default type is
  3.1436 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  3.1437 +#ifdef DOXYGEN
  3.1438 +  template <typename GR, typename WM>
  3.1439 +#else
  3.1440 +  template <typename GR,
  3.1441 +            typename WM = typename GR::template EdgeMap<int> >
  3.1442 +#endif
  3.1443 +  class MaxWeightedPerfectFractionalMatching {
  3.1444 +  public:
  3.1445 +
  3.1446 +    /// The graph type of the algorithm
  3.1447 +    typedef GR Graph;
  3.1448 +    /// The type of the edge weight map
  3.1449 +    typedef WM WeightMap;
  3.1450 +    /// The value type of the edge weights
  3.1451 +    typedef typename WeightMap::Value Value;
  3.1452 +
  3.1453 +    /// The type of the matching map
  3.1454 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
  3.1455 +    MatchingMap;
  3.1456 +
  3.1457 +    /// \brief Scaling factor for primal solution
  3.1458 +    ///
  3.1459 +    /// Scaling factor for primal solution.
  3.1460 +    static const int primalScale = 2;
  3.1461 +
  3.1462 +    /// \brief Scaling factor for dual solution
  3.1463 +    ///
  3.1464 +    /// Scaling factor for dual solution. It is equal to 4 or 1
  3.1465 +    /// according to the value type.
  3.1466 +    static const int dualScale =
  3.1467 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
  3.1468 +
  3.1469 +  private:
  3.1470 +
  3.1471 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
  3.1472 +
  3.1473 +    typedef typename Graph::template NodeMap<Value> NodePotential;
  3.1474 +
  3.1475 +    const Graph& _graph;
  3.1476 +    const WeightMap& _weight;
  3.1477 +
  3.1478 +    MatchingMap* _matching;
  3.1479 +    NodePotential* _node_potential;
  3.1480 +
  3.1481 +    int _node_num;
  3.1482 +    bool _allow_loops;
  3.1483 +
  3.1484 +    enum Status {
  3.1485 +      EVEN = -1, MATCHED = 0, ODD = 1
  3.1486 +    };
  3.1487 +
  3.1488 +    typedef typename Graph::template NodeMap<Status> StatusMap;
  3.1489 +    StatusMap* _status;
  3.1490 +
  3.1491 +    typedef typename Graph::template NodeMap<Arc> PredMap;
  3.1492 +    PredMap* _pred;
  3.1493 +
  3.1494 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
  3.1495 +
  3.1496 +    IntNodeMap *_tree_set_index;
  3.1497 +    TreeSet *_tree_set;
  3.1498 +
  3.1499 +    IntNodeMap *_delta2_index;
  3.1500 +    BinHeap<Value, IntNodeMap> *_delta2;
  3.1501 +
  3.1502 +    IntEdgeMap *_delta3_index;
  3.1503 +    BinHeap<Value, IntEdgeMap> *_delta3;
  3.1504 +
  3.1505 +    Value _delta_sum;
  3.1506 +
  3.1507 +    void createStructures() {
  3.1508 +      _node_num = countNodes(_graph);
  3.1509 +
  3.1510 +      if (!_matching) {
  3.1511 +        _matching = new MatchingMap(_graph);
  3.1512 +      }
  3.1513 +      if (!_node_potential) {
  3.1514 +        _node_potential = new NodePotential(_graph);
  3.1515 +      }
  3.1516 +      if (!_status) {
  3.1517 +        _status = new StatusMap(_graph);
  3.1518 +      }
  3.1519 +      if (!_pred) {
  3.1520 +        _pred = new PredMap(_graph);
  3.1521 +      }
  3.1522 +      if (!_tree_set) {
  3.1523 +        _tree_set_index = new IntNodeMap(_graph);
  3.1524 +        _tree_set = new TreeSet(*_tree_set_index);
  3.1525 +      }
  3.1526 +      if (!_delta2) {
  3.1527 +        _delta2_index = new IntNodeMap(_graph);
  3.1528 +        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
  3.1529 +      }
  3.1530 +      if (!_delta3) {
  3.1531 +        _delta3_index = new IntEdgeMap(_graph);
  3.1532 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  3.1533 +      }
  3.1534 +    }
  3.1535 +
  3.1536 +    void destroyStructures() {
  3.1537 +      if (_matching) {
  3.1538 +        delete _matching;
  3.1539 +      }
  3.1540 +      if (_node_potential) {
  3.1541 +        delete _node_potential;
  3.1542 +      }
  3.1543 +      if (_status) {
  3.1544 +        delete _status;
  3.1545 +      }
  3.1546 +      if (_pred) {
  3.1547 +        delete _pred;
  3.1548 +      }
  3.1549 +      if (_tree_set) {
  3.1550 +        delete _tree_set_index;
  3.1551 +        delete _tree_set;
  3.1552 +      }
  3.1553 +      if (_delta2) {
  3.1554 +        delete _delta2_index;
  3.1555 +        delete _delta2;
  3.1556 +      }
  3.1557 +      if (_delta3) {
  3.1558 +        delete _delta3_index;
  3.1559 +        delete _delta3;
  3.1560 +      }
  3.1561 +    }
  3.1562 +
  3.1563 +    void matchedToEven(Node node, int tree) {
  3.1564 +      _tree_set->insert(node, tree);
  3.1565 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  3.1566 +
  3.1567 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
  3.1568 +        _delta2->erase(node);
  3.1569 +      }
  3.1570 +
  3.1571 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
  3.1572 +        Node v = _graph.source(a);
  3.1573 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  3.1574 +          dualScale * _weight[a];
  3.1575 +        if (node == v) {
  3.1576 +          if (_allow_loops && _graph.direction(a)) {
  3.1577 +            _delta3->push(a, rw / 2);
  3.1578 +          }
  3.1579 +        } else if ((*_status)[v] == EVEN) {
  3.1580 +          _delta3->push(a, rw / 2);
  3.1581 +        } else if ((*_status)[v] == MATCHED) {
  3.1582 +          if (_delta2->state(v) != _delta2->IN_HEAP) {
  3.1583 +            _pred->set(v, a);
  3.1584 +            _delta2->push(v, rw);
  3.1585 +          } else if ((*_delta2)[v] > rw) {
  3.1586 +            _pred->set(v, a);
  3.1587 +            _delta2->decrease(v, rw);
  3.1588 +          }
  3.1589 +        }
  3.1590 +      }
  3.1591 +    }
  3.1592 +
  3.1593 +    void matchedToOdd(Node node, int tree) {
  3.1594 +      _tree_set->insert(node, tree);
  3.1595 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  3.1596 +
  3.1597 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
  3.1598 +        _delta2->erase(node);
  3.1599 +      }
  3.1600 +    }
  3.1601 +
  3.1602 +    void evenToMatched(Node node, int tree) {
  3.1603 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  3.1604 +      Arc min = INVALID;
  3.1605 +      Value minrw = std::numeric_limits<Value>::max();
  3.1606 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
  3.1607 +        Node v = _graph.source(a);
  3.1608 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  3.1609 +          dualScale * _weight[a];
  3.1610 +
  3.1611 +        if (node == v) {
  3.1612 +          if (_allow_loops && _graph.direction(a)) {
  3.1613 +            _delta3->erase(a);
  3.1614 +          }
  3.1615 +        } else if ((*_status)[v] == EVEN) {
  3.1616 +          _delta3->erase(a);
  3.1617 +          if (minrw > rw) {
  3.1618 +            min = _graph.oppositeArc(a);
  3.1619 +            minrw = rw;
  3.1620 +          }
  3.1621 +        } else if ((*_status)[v]  == MATCHED) {
  3.1622 +          if ((*_pred)[v] == a) {
  3.1623 +            Arc mina = INVALID;
  3.1624 +            Value minrwa = std::numeric_limits<Value>::max();
  3.1625 +            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
  3.1626 +              Node va = _graph.target(aa);
  3.1627 +              if ((*_status)[va] != EVEN ||
  3.1628 +                  _tree_set->find(va) == tree) continue;
  3.1629 +              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
  3.1630 +                dualScale * _weight[aa];
  3.1631 +              if (minrwa > rwa) {
  3.1632 +                minrwa = rwa;
  3.1633 +                mina = aa;
  3.1634 +              }
  3.1635 +            }
  3.1636 +            if (mina != INVALID) {
  3.1637 +              _pred->set(v, mina);
  3.1638 +              _delta2->increase(v, minrwa);
  3.1639 +            } else {
  3.1640 +              _pred->set(v, INVALID);
  3.1641 +              _delta2->erase(v);
  3.1642 +            }
  3.1643 +          }
  3.1644 +        }
  3.1645 +      }
  3.1646 +      if (min != INVALID) {
  3.1647 +        _pred->set(node, min);
  3.1648 +        _delta2->push(node, minrw);
  3.1649 +      } else {
  3.1650 +        _pred->set(node, INVALID);
  3.1651 +      }
  3.1652 +    }
  3.1653 +
  3.1654 +    void oddToMatched(Node node) {
  3.1655 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  3.1656 +      Arc min = INVALID;
  3.1657 +      Value minrw = std::numeric_limits<Value>::max();
  3.1658 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
  3.1659 +        Node v = _graph.source(a);
  3.1660 +        if ((*_status)[v] != EVEN) continue;
  3.1661 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  3.1662 +          dualScale * _weight[a];
  3.1663 +
  3.1664 +        if (minrw > rw) {
  3.1665 +          min = _graph.oppositeArc(a);
  3.1666 +          minrw = rw;
  3.1667 +        }
  3.1668 +      }
  3.1669 +      if (min != INVALID) {
  3.1670 +        _pred->set(node, min);
  3.1671 +        _delta2->push(node, minrw);
  3.1672 +      } else {
  3.1673 +        _pred->set(node, INVALID);
  3.1674 +      }
  3.1675 +    }
  3.1676 +
  3.1677 +    void alternatePath(Node even, int tree) {
  3.1678 +      Node odd;
  3.1679 +
  3.1680 +      _status->set(even, MATCHED);
  3.1681 +      evenToMatched(even, tree);
  3.1682 +
  3.1683 +      Arc prev = (*_matching)[even];
  3.1684 +      while (prev != INVALID) {
  3.1685 +        odd = _graph.target(prev);
  3.1686 +        even = _graph.target((*_pred)[odd]);
  3.1687 +        _matching->set(odd, (*_pred)[odd]);
  3.1688 +        _status->set(odd, MATCHED);
  3.1689 +        oddToMatched(odd);
  3.1690 +
  3.1691 +        prev = (*_matching)[even];
  3.1692 +        _status->set(even, MATCHED);
  3.1693 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
  3.1694 +        evenToMatched(even, tree);
  3.1695 +      }
  3.1696 +    }
  3.1697 +
  3.1698 +    void destroyTree(int tree) {
  3.1699 +      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
  3.1700 +        if ((*_status)[n] == EVEN) {
  3.1701 +          _status->set(n, MATCHED);
  3.1702 +          evenToMatched(n, tree);
  3.1703 +        } else if ((*_status)[n] == ODD) {
  3.1704 +          _status->set(n, MATCHED);
  3.1705 +          oddToMatched(n);
  3.1706 +        }
  3.1707 +      }
  3.1708 +      _tree_set->eraseClass(tree);
  3.1709 +    }
  3.1710 +
  3.1711 +    void augmentOnEdge(const Edge& edge) {
  3.1712 +      Node left = _graph.u(edge);
  3.1713 +      int left_tree = _tree_set->find(left);
  3.1714 +
  3.1715 +      alternatePath(left, left_tree);
  3.1716 +      destroyTree(left_tree);
  3.1717 +      _matching->set(left, _graph.direct(edge, true));
  3.1718 +
  3.1719 +      Node right = _graph.v(edge);
  3.1720 +      int right_tree = _tree_set->find(right);
  3.1721 +
  3.1722 +      alternatePath(right, right_tree);
  3.1723 +      destroyTree(right_tree);
  3.1724 +      _matching->set(right, _graph.direct(edge, false));
  3.1725 +    }
  3.1726 +
  3.1727 +    void augmentOnArc(const Arc& arc) {
  3.1728 +      Node left = _graph.source(arc);
  3.1729 +      _status->set(left, MATCHED);
  3.1730 +      _matching->set(left, arc);
  3.1731 +      _pred->set(left, arc);
  3.1732 +
  3.1733 +      Node right = _graph.target(arc);
  3.1734 +      int right_tree = _tree_set->find(right);
  3.1735 +
  3.1736 +      alternatePath(right, right_tree);
  3.1737 +      destroyTree(right_tree);
  3.1738 +      _matching->set(right, _graph.oppositeArc(arc));
  3.1739 +    }
  3.1740 +
  3.1741 +    void extendOnArc(const Arc& arc) {
  3.1742 +      Node base = _graph.target(arc);
  3.1743 +      int tree = _tree_set->find(base);
  3.1744 +
  3.1745 +      Node odd = _graph.source(arc);
  3.1746 +      _tree_set->insert(odd, tree);
  3.1747 +      _status->set(odd, ODD);
  3.1748 +      matchedToOdd(odd, tree);
  3.1749 +      _pred->set(odd, arc);
  3.1750 +
  3.1751 +      Node even = _graph.target((*_matching)[odd]);
  3.1752 +      _tree_set->insert(even, tree);
  3.1753 +      _status->set(even, EVEN);
  3.1754 +      matchedToEven(even, tree);
  3.1755 +    }
  3.1756 +
  3.1757 +    void cycleOnEdge(const Edge& edge, int tree) {
  3.1758 +      Node nca = INVALID;
  3.1759 +      std::vector<Node> left_path, right_path;
  3.1760 +
  3.1761 +      {
  3.1762 +        std::set<Node> left_set, right_set;
  3.1763 +        Node left = _graph.u(edge);
  3.1764 +        left_path.push_back(left);
  3.1765 +        left_set.insert(left);
  3.1766 +
  3.1767 +        Node right = _graph.v(edge);
  3.1768 +        right_path.push_back(right);
  3.1769 +        right_set.insert(right);
  3.1770 +
  3.1771 +        while (true) {
  3.1772 +
  3.1773 +          if (left_set.find(right) != left_set.end()) {
  3.1774 +            nca = right;
  3.1775 +            break;
  3.1776 +          }
  3.1777 +
  3.1778 +          if ((*_matching)[left] == INVALID) break;
  3.1779 +
  3.1780 +          left = _graph.target((*_matching)[left]);
  3.1781 +          left_path.push_back(left);
  3.1782 +          left = _graph.target((*_pred)[left]);
  3.1783 +          left_path.push_back(left);
  3.1784 +
  3.1785 +          left_set.insert(left);
  3.1786 +
  3.1787 +          if (right_set.find(left) != right_set.end()) {
  3.1788 +            nca = left;
  3.1789 +            break;
  3.1790 +          }
  3.1791 +
  3.1792 +          if ((*_matching)[right] == INVALID) break;
  3.1793 +
  3.1794 +          right = _graph.target((*_matching)[right]);
  3.1795 +          right_path.push_back(right);
  3.1796 +          right = _graph.target((*_pred)[right]);
  3.1797 +          right_path.push_back(right);
  3.1798 +
  3.1799 +          right_set.insert(right);
  3.1800 +
  3.1801 +        }
  3.1802 +
  3.1803 +        if (nca == INVALID) {
  3.1804 +          if ((*_matching)[left] == INVALID) {
  3.1805 +            nca = right;
  3.1806 +            while (left_set.find(nca) == left_set.end()) {
  3.1807 +              nca = _graph.target((*_matching)[nca]);
  3.1808 +              right_path.push_back(nca);
  3.1809 +              nca = _graph.target((*_pred)[nca]);
  3.1810 +              right_path.push_back(nca);
  3.1811 +            }
  3.1812 +          } else {
  3.1813 +            nca = left;
  3.1814 +            while (right_set.find(nca) == right_set.end()) {
  3.1815 +              nca = _graph.target((*_matching)[nca]);
  3.1816 +              left_path.push_back(nca);
  3.1817 +              nca = _graph.target((*_pred)[nca]);
  3.1818 +              left_path.push_back(nca);
  3.1819 +            }
  3.1820 +          }
  3.1821 +        }
  3.1822 +      }
  3.1823 +
  3.1824 +      alternatePath(nca, tree);
  3.1825 +      Arc prev;
  3.1826 +
  3.1827 +      prev = _graph.direct(edge, true);
  3.1828 +      for (int i = 0; left_path[i] != nca; i += 2) {
  3.1829 +        _matching->set(left_path[i], prev);
  3.1830 +        _status->set(left_path[i], MATCHED);
  3.1831 +        evenToMatched(left_path[i], tree);
  3.1832 +
  3.1833 +        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  3.1834 +        _status->set(left_path[i + 1], MATCHED);
  3.1835 +        oddToMatched(left_path[i + 1]);
  3.1836 +      }
  3.1837 +      _matching->set(nca, prev);
  3.1838 +
  3.1839 +      for (int i = 0; right_path[i] != nca; i += 2) {
  3.1840 +        _status->set(right_path[i], MATCHED);
  3.1841 +        evenToMatched(right_path[i], tree);
  3.1842 +
  3.1843 +        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  3.1844 +        _status->set(right_path[i + 1], MATCHED);
  3.1845 +        oddToMatched(right_path[i + 1]);
  3.1846 +      }
  3.1847 +
  3.1848 +      destroyTree(tree);
  3.1849 +    }
  3.1850 +
  3.1851 +    void extractCycle(const Arc &arc) {
  3.1852 +      Node left = _graph.source(arc);
  3.1853 +      Node odd = _graph.target((*_matching)[left]);
  3.1854 +      Arc prev;
  3.1855 +      while (odd != left) {
  3.1856 +        Node even = _graph.target((*_matching)[odd]);
  3.1857 +        prev = (*_matching)[odd];
  3.1858 +        odd = _graph.target((*_matching)[even]);
  3.1859 +        _matching->set(even, _graph.oppositeArc(prev));
  3.1860 +      }
  3.1861 +      _matching->set(left, arc);
  3.1862 +
  3.1863 +      Node right = _graph.target(arc);
  3.1864 +      int right_tree = _tree_set->find(right);
  3.1865 +      alternatePath(right, right_tree);
  3.1866 +      destroyTree(right_tree);
  3.1867 +      _matching->set(right, _graph.oppositeArc(arc));
  3.1868 +    }
  3.1869 +
  3.1870 +  public:
  3.1871 +
  3.1872 +    /// \brief Constructor
  3.1873 +    ///
  3.1874 +    /// Constructor.
  3.1875 +    MaxWeightedPerfectFractionalMatching(const Graph& graph,
  3.1876 +                                         const WeightMap& weight,
  3.1877 +                                         bool allow_loops = true)
  3.1878 +      : _graph(graph), _weight(weight), _matching(0),
  3.1879 +      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  3.1880 +      _status(0),  _pred(0),
  3.1881 +      _tree_set_index(0), _tree_set(0),
  3.1882 +
  3.1883 +      _delta2_index(0), _delta2(0),
  3.1884 +      _delta3_index(0), _delta3(0),
  3.1885 +
  3.1886 +      _delta_sum() {}
  3.1887 +
  3.1888 +    ~MaxWeightedPerfectFractionalMatching() {
  3.1889 +      destroyStructures();
  3.1890 +    }
  3.1891 +
  3.1892 +    /// \name Execution Control
  3.1893 +    /// The simplest way to execute the algorithm is to use the
  3.1894 +    /// \ref run() member function.
  3.1895 +
  3.1896 +    ///@{
  3.1897 +
  3.1898 +    /// \brief Initialize the algorithm
  3.1899 +    ///
  3.1900 +    /// This function initializes the algorithm.
  3.1901 +    void init() {
  3.1902 +      createStructures();
  3.1903 +
  3.1904 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1905 +        (*_delta2_index)[n] = _delta2->PRE_HEAP;
  3.1906 +      }
  3.1907 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  3.1908 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  3.1909 +      }
  3.1910 +
  3.1911 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.1912 +        Value max = - std::numeric_limits<Value>::max();
  3.1913 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  3.1914 +          if (_graph.target(e) == n && !_allow_loops) continue;
  3.1915 +          if ((dualScale * _weight[e]) / 2 > max) {
  3.1916 +            max = (dualScale * _weight[e]) / 2;
  3.1917 +          }
  3.1918 +        }
  3.1919 +        _node_potential->set(n, max);
  3.1920 +
  3.1921 +        _tree_set->insert(n);
  3.1922 +
  3.1923 +        _matching->set(n, INVALID);
  3.1924 +        _status->set(n, EVEN);
  3.1925 +      }
  3.1926 +
  3.1927 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  3.1928 +        Node left = _graph.u(e);
  3.1929 +        Node right = _graph.v(e);
  3.1930 +        if (left == right && !_allow_loops) continue;
  3.1931 +        _delta3->push(e, ((*_node_potential)[left] +
  3.1932 +                          (*_node_potential)[right] -
  3.1933 +                          dualScale * _weight[e]) / 2);
  3.1934 +      }
  3.1935 +    }
  3.1936 +
  3.1937 +    /// \brief Start the algorithm
  3.1938 +    ///
  3.1939 +    /// This function starts the algorithm.
  3.1940 +    ///
  3.1941 +    /// \pre \ref init() must be called before using this function.
  3.1942 +    bool start() {
  3.1943 +      enum OpType {
  3.1944 +        D2, D3
  3.1945 +      };
  3.1946 +
  3.1947 +      int unmatched = _node_num;
  3.1948 +      while (unmatched > 0) {
  3.1949 +        Value d2 = !_delta2->empty() ?
  3.1950 +          _delta2->prio() : std::numeric_limits<Value>::max();
  3.1951 +
  3.1952 +        Value d3 = !_delta3->empty() ?
  3.1953 +          _delta3->prio() : std::numeric_limits<Value>::max();
  3.1954 +
  3.1955 +        _delta_sum = d3; OpType ot = D3;
  3.1956 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  3.1957 +
  3.1958 +        if (_delta_sum == std::numeric_limits<Value>::max()) {
  3.1959 +          return false;
  3.1960 +        }
  3.1961 +
  3.1962 +        switch (ot) {
  3.1963 +        case D2:
  3.1964 +          {
  3.1965 +            Node n = _delta2->top();
  3.1966 +            Arc a = (*_pred)[n];
  3.1967 +            if ((*_matching)[n] == INVALID) {
  3.1968 +              augmentOnArc(a);
  3.1969 +              --unmatched;
  3.1970 +            } else {
  3.1971 +              Node v = _graph.target((*_matching)[n]);
  3.1972 +              if ((*_matching)[n] !=
  3.1973 +                  _graph.oppositeArc((*_matching)[v])) {
  3.1974 +                extractCycle(a);
  3.1975 +                --unmatched;
  3.1976 +              } else {
  3.1977 +                extendOnArc(a);
  3.1978 +              }
  3.1979 +            }
  3.1980 +          } break;
  3.1981 +        case D3:
  3.1982 +          {
  3.1983 +            Edge e = _delta3->top();
  3.1984 +
  3.1985 +            Node left = _graph.u(e);
  3.1986 +            Node right = _graph.v(e);
  3.1987 +
  3.1988 +            int left_tree = _tree_set->find(left);
  3.1989 +            int right_tree = _tree_set->find(right);
  3.1990 +
  3.1991 +            if (left_tree == right_tree) {
  3.1992 +              cycleOnEdge(e, left_tree);
  3.1993 +              --unmatched;
  3.1994 +            } else {
  3.1995 +              augmentOnEdge(e);
  3.1996 +              unmatched -= 2;
  3.1997 +            }
  3.1998 +          } break;
  3.1999 +        }
  3.2000 +      }
  3.2001 +      return true;
  3.2002 +    }
  3.2003 +
  3.2004 +    /// \brief Run the algorithm.
  3.2005 +    ///
  3.2006 +    /// This method runs the \c %MaxWeightedPerfectFractionalMatching 
  3.2007 +    /// algorithm.
  3.2008 +    ///
  3.2009 +    /// \note mwfm.run() is just a shortcut of the following code.
  3.2010 +    /// \code
  3.2011 +    ///   mwpfm.init();
  3.2012 +    ///   mwpfm.start();
  3.2013 +    /// \endcode
  3.2014 +    bool run() {
  3.2015 +      init();
  3.2016 +      return start();
  3.2017 +    }
  3.2018 +
  3.2019 +    /// @}
  3.2020 +
  3.2021 +    /// \name Primal Solution
  3.2022 +    /// Functions to get the primal solution, i.e. the maximum weighted
  3.2023 +    /// matching.\n
  3.2024 +    /// Either \ref run() or \ref start() function should be called before
  3.2025 +    /// using them.
  3.2026 +
  3.2027 +    /// @{
  3.2028 +
  3.2029 +    /// \brief Return the weight of the matching.
  3.2030 +    ///
  3.2031 +    /// This function returns the weight of the found matching. This
  3.2032 +    /// value is scaled by \ref primalScale "primal scale".
  3.2033 +    ///
  3.2034 +    /// \pre Either run() or start() must be called before using this function.
  3.2035 +    Value matchingWeight() const {
  3.2036 +      Value sum = 0;
  3.2037 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.2038 +        if ((*_matching)[n] != INVALID) {
  3.2039 +          sum += _weight[(*_matching)[n]];
  3.2040 +        }
  3.2041 +      }
  3.2042 +      return sum * primalScale / 2;
  3.2043 +    }
  3.2044 +
  3.2045 +    /// \brief Return the number of covered nodes in the matching.
  3.2046 +    ///
  3.2047 +    /// This function returns the number of covered nodes in the matching.
  3.2048 +    ///
  3.2049 +    /// \pre Either run() or start() must be called before using this function.
  3.2050 +    int matchingSize() const {
  3.2051 +      int num = 0;
  3.2052 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.2053 +        if ((*_matching)[n] != INVALID) {
  3.2054 +          ++num;
  3.2055 +        }
  3.2056 +      }
  3.2057 +      return num;
  3.2058 +    }
  3.2059 +
  3.2060 +    /// \brief Return \c true if the given edge is in the matching.
  3.2061 +    ///
  3.2062 +    /// This function returns \c true if the given edge is in the
  3.2063 +    /// found matching. The result is scaled by \ref primalScale
  3.2064 +    /// "primal scale".
  3.2065 +    ///
  3.2066 +    /// \pre Either run() or start() must be called before using this function.
  3.2067 +    int matching(const Edge& edge) const {
  3.2068 +      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  3.2069 +        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
  3.2070 +    }
  3.2071 +
  3.2072 +    /// \brief Return the fractional matching arc (or edge) incident
  3.2073 +    /// to the given node.
  3.2074 +    ///
  3.2075 +    /// This function returns one of the fractional matching arc (or
  3.2076 +    /// edge) incident to the given node in the found matching or \c
  3.2077 +    /// INVALID if the node is not covered by the matching or if the
  3.2078 +    /// node is on an odd length cycle then it is the successor edge
  3.2079 +    /// on the cycle.
  3.2080 +    ///
  3.2081 +    /// \pre Either run() or start() must be called before using this function.
  3.2082 +    Arc matching(const Node& node) const {
  3.2083 +      return (*_matching)[node];
  3.2084 +    }
  3.2085 +
  3.2086 +    /// \brief Return a const reference to the matching map.
  3.2087 +    ///
  3.2088 +    /// This function returns a const reference to a node map that stores
  3.2089 +    /// the matching arc (or edge) incident to each node.
  3.2090 +    const MatchingMap& matchingMap() const {
  3.2091 +      return *_matching;
  3.2092 +    }
  3.2093 +
  3.2094 +    /// @}
  3.2095 +
  3.2096 +    /// \name Dual Solution
  3.2097 +    /// Functions to get the dual solution.\n
  3.2098 +    /// Either \ref run() or \ref start() function should be called before
  3.2099 +    /// using them.
  3.2100 +
  3.2101 +    /// @{
  3.2102 +
  3.2103 +    /// \brief Return the value of the dual solution.
  3.2104 +    ///
  3.2105 +    /// This function returns the value of the dual solution.
  3.2106 +    /// It should be equal to the primal value scaled by \ref dualScale
  3.2107 +    /// "dual scale".
  3.2108 +    ///
  3.2109 +    /// \pre Either run() or start() must be called before using this function.
  3.2110 +    Value dualValue() const {
  3.2111 +      Value sum = 0;
  3.2112 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  3.2113 +        sum += nodeValue(n);
  3.2114 +      }
  3.2115 +      return sum;
  3.2116 +    }
  3.2117 +
  3.2118 +    /// \brief Return the dual value (potential) of the given node.
  3.2119 +    ///
  3.2120 +    /// This function returns the dual value (potential) of the given node.
  3.2121 +    ///
  3.2122 +    /// \pre Either run() or start() must be called before using this function.
  3.2123 +    Value nodeValue(const Node& n) const {
  3.2124 +      return (*_node_potential)[n];
  3.2125 +    }
  3.2126 +
  3.2127 +    /// @}
  3.2128 +
  3.2129 +  };
  3.2130 +
  3.2131 +} //END OF NAMESPACE LEMON
  3.2132 +
  3.2133 +#endif //LEMON_FRACTIONAL_MATCHING_H
     4.1 --- a/lemon/matching.h	Tue Mar 16 21:18:39 2010 +0100
     4.2 +++ b/lemon/matching.h	Tue Mar 16 21:27:35 2010 +0100
     4.3 @@ -16,8 +16,8 @@
     4.4   *
     4.5   */
     4.6  
     4.7 -#ifndef LEMON_MAX_MATCHING_H
     4.8 -#define LEMON_MAX_MATCHING_H
     4.9 +#ifndef LEMON_MATCHING_H
    4.10 +#define LEMON_MATCHING_H
    4.11  
    4.12  #include <vector>
    4.13  #include <queue>
    4.14 @@ -28,6 +28,7 @@
    4.15  #include <lemon/unionfind.h>
    4.16  #include <lemon/bin_heap.h>
    4.17  #include <lemon/maps.h>
    4.18 +#include <lemon/fractional_matching.h>
    4.19  
    4.20  ///\ingroup matching
    4.21  ///\file
    4.22 @@ -41,7 +42,7 @@
    4.23    ///
    4.24    /// This class implements Edmonds' alternating forest matching algorithm
    4.25    /// for finding a maximum cardinality matching in a general undirected graph.
    4.26 -  /// It can be started from an arbitrary initial matching 
    4.27 +  /// It can be started from an arbitrary initial matching
    4.28    /// (the default is the empty one).
    4.29    ///
    4.30    /// The dual solution of the problem is a map of the nodes to
    4.31 @@ -69,11 +70,11 @@
    4.32  
    4.33      ///\brief Status constants for Gallai-Edmonds decomposition.
    4.34      ///
    4.35 -    ///These constants are used for indicating the Gallai-Edmonds 
    4.36 +    ///These constants are used for indicating the Gallai-Edmonds
    4.37      ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    4.38      ///induce a subgraph with factor-critical components, the nodes with
    4.39      ///status \c ODD (or \c A) form the canonical barrier, and the nodes
    4.40 -    ///with status \c MATCHED (or \c C) induce a subgraph having a 
    4.41 +    ///with status \c MATCHED (or \c C) induce a subgraph having a
    4.42      ///perfect matching.
    4.43      enum Status {
    4.44        EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
    4.45 @@ -512,7 +513,7 @@
    4.46        }
    4.47      }
    4.48  
    4.49 -    /// \brief Start Edmonds' algorithm with a heuristic improvement 
    4.50 +    /// \brief Start Edmonds' algorithm with a heuristic improvement
    4.51      /// for dense graphs
    4.52      ///
    4.53      /// This function runs Edmonds' algorithm with a heuristic of postponing
    4.54 @@ -534,8 +535,8 @@
    4.55  
    4.56      /// \brief Run Edmonds' algorithm
    4.57      ///
    4.58 -    /// This function runs Edmonds' algorithm. An additional heuristic of 
    4.59 -    /// postponing shrinks is used for relatively dense graphs 
    4.60 +    /// This function runs Edmonds' algorithm. An additional heuristic of
    4.61 +    /// postponing shrinks is used for relatively dense graphs
    4.62      /// (for which <tt>m>=2*n</tt> holds).
    4.63      void run() {
    4.64        if (countEdges(_graph) < 2 * countNodes(_graph)) {
    4.65 @@ -556,7 +557,7 @@
    4.66  
    4.67      /// \brief Return the size (cardinality) of the matching.
    4.68      ///
    4.69 -    /// This function returns the size (cardinality) of the current matching. 
    4.70 +    /// This function returns the size (cardinality) of the current matching.
    4.71      /// After run() it returns the size of the maximum matching in the graph.
    4.72      int matchingSize() const {
    4.73        int size = 0;
    4.74 @@ -570,7 +571,7 @@
    4.75  
    4.76      /// \brief Return \c true if the given edge is in the matching.
    4.77      ///
    4.78 -    /// This function returns \c true if the given edge is in the current 
    4.79 +    /// This function returns \c true if the given edge is in the current
    4.80      /// matching.
    4.81      bool matching(const Edge& edge) const {
    4.82        return edge == (*_matching)[_graph.u(edge)];
    4.83 @@ -579,7 +580,7 @@
    4.84      /// \brief Return the matching arc (or edge) incident to the given node.
    4.85      ///
    4.86      /// This function returns the matching arc (or edge) incident to the
    4.87 -    /// given node in the current matching or \c INVALID if the node is 
    4.88 +    /// given node in the current matching or \c INVALID if the node is
    4.89      /// not covered by the matching.
    4.90      Arc matching(const Node& n) const {
    4.91        return (*_matching)[n];
    4.92 @@ -595,7 +596,7 @@
    4.93  
    4.94      /// \brief Return the mate of the given node.
    4.95      ///
    4.96 -    /// This function returns the mate of the given node in the current 
    4.97 +    /// This function returns the mate of the given node in the current
    4.98      /// matching or \c INVALID if the node is not covered by the matching.
    4.99      Node mate(const Node& n) const {
   4.100        return (*_matching)[n] != INVALID ?
   4.101 @@ -605,7 +606,7 @@
   4.102      /// @}
   4.103  
   4.104      /// \name Dual Solution
   4.105 -    /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
   4.106 +    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
   4.107      /// decomposition.
   4.108  
   4.109      /// @{
   4.110 @@ -648,8 +649,8 @@
   4.111    /// on extensive use of priority queues and provides
   4.112    /// \f$O(nm\log n)\f$ time complexity.
   4.113    ///
   4.114 -  /// The maximum weighted matching problem is to find a subset of the 
   4.115 -  /// edges in an undirected graph with maximum overall weight for which 
   4.116 +  /// The maximum weighted matching problem is to find a subset of the
   4.117 +  /// edges in an undirected graph with maximum overall weight for which
   4.118    /// each node has at most one incident edge.
   4.119    /// It can be formulated with the following linear program.
   4.120    /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   4.121 @@ -673,16 +674,16 @@
   4.122    /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
   4.123        \frac{\vert B \vert - 1}{2}z_B\f] */
   4.124    ///
   4.125 -  /// The algorithm can be executed with the run() function. 
   4.126 +  /// The algorithm can be executed with the run() function.
   4.127    /// After it the matching (the primal solution) and the dual solution
   4.128 -  /// can be obtained using the query functions and the 
   4.129 -  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
   4.130 -  /// which is able to iterate on the nodes of a blossom. 
   4.131 +  /// can be obtained using the query functions and the
   4.132 +  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
   4.133 +  /// which is able to iterate on the nodes of a blossom.
   4.134    /// If the value type is integer, then the dual solution is multiplied
   4.135    /// by \ref MaxWeightedMatching::dualScale "4".
   4.136    ///
   4.137    /// \tparam GR The undirected graph type the algorithm runs on.
   4.138 -  /// \tparam WM The type edge weight map. The default type is 
   4.139 +  /// \tparam WM The type edge weight map. The default type is
   4.140    /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   4.141  #ifdef DOXYGEN
   4.142    template <typename GR, typename WM>
   4.143 @@ -745,7 +746,7 @@
   4.144      typedef RangeMap<int> IntIntMap;
   4.145  
   4.146      enum Status {
   4.147 -      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
   4.148 +      EVEN = -1, MATCHED = 0, ODD = 1
   4.149      };
   4.150  
   4.151      typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
   4.152 @@ -797,6 +798,10 @@
   4.153      BinHeap<Value, IntIntMap> *_delta4;
   4.154  
   4.155      Value _delta_sum;
   4.156 +    int _unmatched;
   4.157 +
   4.158 +    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
   4.159 +    FractionalMatching *_fractional;
   4.160  
   4.161      void createStructures() {
   4.162        _node_num = countNodes(_graph);
   4.163 @@ -844,9 +849,6 @@
   4.164      }
   4.165  
   4.166      void destroyStructures() {
   4.167 -      _node_num = countNodes(_graph);
   4.168 -      _blossom_num = _node_num * 3 / 2;
   4.169 -
   4.170        if (_matching) {
   4.171          delete _matching;
   4.172        }
   4.173 @@ -922,10 +924,6 @@
   4.174              if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   4.175                _delta3->push(e, rw / 2);
   4.176              }
   4.177 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   4.178 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   4.179 -              _delta3->push(e, rw);
   4.180 -            }
   4.181            } else {
   4.182              typename std::map<int, Arc>::iterator it =
   4.183                (*_node_data)[vi].heap_index.find(tree);
   4.184 @@ -949,202 +947,6 @@
   4.185                    _delta2->push(vb, _blossom_set->classPrio(vb) -
   4.186                                 (*_blossom_data)[vb].offset);
   4.187                  } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   4.188 -                           (*_blossom_data)[vb].offset){
   4.189 -                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   4.190 -                                   (*_blossom_data)[vb].offset);
   4.191 -                }
   4.192 -              }
   4.193 -            }
   4.194 -          }
   4.195 -        }
   4.196 -      }
   4.197 -      (*_blossom_data)[blossom].offset = 0;
   4.198 -    }
   4.199 -
   4.200 -    void matchedToOdd(int blossom) {
   4.201 -      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   4.202 -        _delta2->erase(blossom);
   4.203 -      }
   4.204 -      (*_blossom_data)[blossom].offset += _delta_sum;
   4.205 -      if (!_blossom_set->trivial(blossom)) {
   4.206 -        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   4.207 -                     (*_blossom_data)[blossom].offset);
   4.208 -      }
   4.209 -    }
   4.210 -
   4.211 -    void evenToMatched(int blossom, int tree) {
   4.212 -      if (!_blossom_set->trivial(blossom)) {
   4.213 -        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   4.214 -      }
   4.215 -
   4.216 -      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   4.217 -           n != INVALID; ++n) {
   4.218 -        int ni = (*_node_index)[n];
   4.219 -        (*_node_data)[ni].pot -= _delta_sum;
   4.220 -
   4.221 -        _delta1->erase(n);
   4.222 -
   4.223 -        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   4.224 -          Node v = _graph.source(e);
   4.225 -          int vb = _blossom_set->find(v);
   4.226 -          int vi = (*_node_index)[v];
   4.227 -
   4.228 -          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   4.229 -            dualScale * _weight[e];
   4.230 -
   4.231 -          if (vb == blossom) {
   4.232 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.233 -              _delta3->erase(e);
   4.234 -            }
   4.235 -          } else if ((*_blossom_data)[vb].status == EVEN) {
   4.236 -
   4.237 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.238 -              _delta3->erase(e);
   4.239 -            }
   4.240 -
   4.241 -            int vt = _tree_set->find(vb);
   4.242 -
   4.243 -            if (vt != tree) {
   4.244 -
   4.245 -              Arc r = _graph.oppositeArc(e);
   4.246 -
   4.247 -              typename std::map<int, Arc>::iterator it =
   4.248 -                (*_node_data)[ni].heap_index.find(vt);
   4.249 -
   4.250 -              if (it != (*_node_data)[ni].heap_index.end()) {
   4.251 -                if ((*_node_data)[ni].heap[it->second] > rw) {
   4.252 -                  (*_node_data)[ni].heap.replace(it->second, r);
   4.253 -                  (*_node_data)[ni].heap.decrease(r, rw);
   4.254 -                  it->second = r;
   4.255 -                }
   4.256 -              } else {
   4.257 -                (*_node_data)[ni].heap.push(r, rw);
   4.258 -                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   4.259 -              }
   4.260 -
   4.261 -              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   4.262 -                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   4.263 -
   4.264 -                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   4.265 -                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   4.266 -                               (*_blossom_data)[blossom].offset);
   4.267 -                } else if ((*_delta2)[blossom] >
   4.268 -                           _blossom_set->classPrio(blossom) -
   4.269 -                           (*_blossom_data)[blossom].offset){
   4.270 -                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   4.271 -                                   (*_blossom_data)[blossom].offset);
   4.272 -                }
   4.273 -              }
   4.274 -            }
   4.275 -
   4.276 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   4.277 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.278 -              _delta3->erase(e);
   4.279 -            }
   4.280 -          } else {
   4.281 -
   4.282 -            typename std::map<int, Arc>::iterator it =
   4.283 -              (*_node_data)[vi].heap_index.find(tree);
   4.284 -
   4.285 -            if (it != (*_node_data)[vi].heap_index.end()) {
   4.286 -              (*_node_data)[vi].heap.erase(it->second);
   4.287 -              (*_node_data)[vi].heap_index.erase(it);
   4.288 -              if ((*_node_data)[vi].heap.empty()) {
   4.289 -                _blossom_set->increase(v, std::numeric_limits<Value>::max());
   4.290 -              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
   4.291 -                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
   4.292 -              }
   4.293 -
   4.294 -              if ((*_blossom_data)[vb].status == MATCHED) {
   4.295 -                if (_blossom_set->classPrio(vb) ==
   4.296 -                    std::numeric_limits<Value>::max()) {
   4.297 -                  _delta2->erase(vb);
   4.298 -                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
   4.299 -                           (*_blossom_data)[vb].offset) {
   4.300 -                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
   4.301 -                                   (*_blossom_data)[vb].offset);
   4.302 -                }
   4.303 -              }
   4.304 -            }
   4.305 -          }
   4.306 -        }
   4.307 -      }
   4.308 -    }
   4.309 -
   4.310 -    void oddToMatched(int blossom) {
   4.311 -      (*_blossom_data)[blossom].offset -= _delta_sum;
   4.312 -
   4.313 -      if (_blossom_set->classPrio(blossom) !=
   4.314 -          std::numeric_limits<Value>::max()) {
   4.315 -        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   4.316 -                       (*_blossom_data)[blossom].offset);
   4.317 -      }
   4.318 -
   4.319 -      if (!_blossom_set->trivial(blossom)) {
   4.320 -        _delta4->erase(blossom);
   4.321 -      }
   4.322 -    }
   4.323 -
   4.324 -    void oddToEven(int blossom, int tree) {
   4.325 -      if (!_blossom_set->trivial(blossom)) {
   4.326 -        _delta4->erase(blossom);
   4.327 -        (*_blossom_data)[blossom].pot -=
   4.328 -          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
   4.329 -      }
   4.330 -
   4.331 -      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   4.332 -           n != INVALID; ++n) {
   4.333 -        int ni = (*_node_index)[n];
   4.334 -
   4.335 -        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   4.336 -
   4.337 -        (*_node_data)[ni].heap.clear();
   4.338 -        (*_node_data)[ni].heap_index.clear();
   4.339 -        (*_node_data)[ni].pot +=
   4.340 -          2 * _delta_sum - (*_blossom_data)[blossom].offset;
   4.341 -
   4.342 -        _delta1->push(n, (*_node_data)[ni].pot);
   4.343 -
   4.344 -        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   4.345 -          Node v = _graph.source(e);
   4.346 -          int vb = _blossom_set->find(v);
   4.347 -          int vi = (*_node_index)[v];
   4.348 -
   4.349 -          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   4.350 -            dualScale * _weight[e];
   4.351 -
   4.352 -          if ((*_blossom_data)[vb].status == EVEN) {
   4.353 -            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   4.354 -              _delta3->push(e, rw / 2);
   4.355 -            }
   4.356 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   4.357 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   4.358 -              _delta3->push(e, rw);
   4.359 -            }
   4.360 -          } else {
   4.361 -
   4.362 -            typename std::map<int, Arc>::iterator it =
   4.363 -              (*_node_data)[vi].heap_index.find(tree);
   4.364 -
   4.365 -            if (it != (*_node_data)[vi].heap_index.end()) {
   4.366 -              if ((*_node_data)[vi].heap[it->second] > rw) {
   4.367 -                (*_node_data)[vi].heap.replace(it->second, e);
   4.368 -                (*_node_data)[vi].heap.decrease(e, rw);
   4.369 -                it->second = e;
   4.370 -              }
   4.371 -            } else {
   4.372 -              (*_node_data)[vi].heap.push(e, rw);
   4.373 -              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   4.374 -            }
   4.375 -
   4.376 -            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   4.377 -              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   4.378 -
   4.379 -              if ((*_blossom_data)[vb].status == MATCHED) {
   4.380 -                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   4.381 -                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   4.382 -                               (*_blossom_data)[vb].offset);
   4.383 -                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   4.384                             (*_blossom_data)[vb].offset) {
   4.385                    _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   4.386                                     (*_blossom_data)[vb].offset);
   4.387 @@ -1157,43 +959,145 @@
   4.388        (*_blossom_data)[blossom].offset = 0;
   4.389      }
   4.390  
   4.391 -
   4.392 -    void matchedToUnmatched(int blossom) {
   4.393 +    void matchedToOdd(int blossom) {
   4.394        if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   4.395          _delta2->erase(blossom);
   4.396        }
   4.397 +      (*_blossom_data)[blossom].offset += _delta_sum;
   4.398 +      if (!_blossom_set->trivial(blossom)) {
   4.399 +        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   4.400 +                      (*_blossom_data)[blossom].offset);
   4.401 +      }
   4.402 +    }
   4.403 +
   4.404 +    void evenToMatched(int blossom, int tree) {
   4.405 +      if (!_blossom_set->trivial(blossom)) {
   4.406 +        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   4.407 +      }
   4.408  
   4.409        for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   4.410             n != INVALID; ++n) {
   4.411          int ni = (*_node_index)[n];
   4.412 -
   4.413 -        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   4.414 -
   4.415 -        (*_node_data)[ni].heap.clear();
   4.416 -        (*_node_data)[ni].heap_index.clear();
   4.417 -
   4.418 -        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   4.419 -          Node v = _graph.target(e);
   4.420 +        (*_node_data)[ni].pot -= _delta_sum;
   4.421 +
   4.422 +        _delta1->erase(n);
   4.423 +
   4.424 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   4.425 +          Node v = _graph.source(e);
   4.426            int vb = _blossom_set->find(v);
   4.427            int vi = (*_node_index)[v];
   4.428  
   4.429            Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   4.430              dualScale * _weight[e];
   4.431  
   4.432 -          if ((*_blossom_data)[vb].status == EVEN) {
   4.433 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   4.434 -              _delta3->push(e, rw);
   4.435 +          if (vb == blossom) {
   4.436 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.437 +              _delta3->erase(e);
   4.438 +            }
   4.439 +          } else if ((*_blossom_data)[vb].status == EVEN) {
   4.440 +
   4.441 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.442 +              _delta3->erase(e);
   4.443 +            }
   4.444 +
   4.445 +            int vt = _tree_set->find(vb);
   4.446 +
   4.447 +            if (vt != tree) {
   4.448 +
   4.449 +              Arc r = _graph.oppositeArc(e);
   4.450 +
   4.451 +              typename std::map<int, Arc>::iterator it =
   4.452 +                (*_node_data)[ni].heap_index.find(vt);
   4.453 +
   4.454 +              if (it != (*_node_data)[ni].heap_index.end()) {
   4.455 +                if ((*_node_data)[ni].heap[it->second] > rw) {
   4.456 +                  (*_node_data)[ni].heap.replace(it->second, r);
   4.457 +                  (*_node_data)[ni].heap.decrease(r, rw);
   4.458 +                  it->second = r;
   4.459 +                }
   4.460 +              } else {
   4.461 +                (*_node_data)[ni].heap.push(r, rw);
   4.462 +                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   4.463 +              }
   4.464 +
   4.465 +              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   4.466 +                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   4.467 +
   4.468 +                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   4.469 +                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   4.470 +                               (*_blossom_data)[blossom].offset);
   4.471 +                } else if ((*_delta2)[blossom] >
   4.472 +                           _blossom_set->classPrio(blossom) -
   4.473 +                           (*_blossom_data)[blossom].offset){
   4.474 +                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   4.475 +                                   (*_blossom_data)[blossom].offset);
   4.476 +                }
   4.477 +              }
   4.478 +            }
   4.479 +          } else {
   4.480 +
   4.481 +            typename std::map<int, Arc>::iterator it =
   4.482 +              (*_node_data)[vi].heap_index.find(tree);
   4.483 +
   4.484 +            if (it != (*_node_data)[vi].heap_index.end()) {
   4.485 +              (*_node_data)[vi].heap.erase(it->second);
   4.486 +              (*_node_data)[vi].heap_index.erase(it);
   4.487 +              if ((*_node_data)[vi].heap.empty()) {
   4.488 +                _blossom_set->increase(v, std::numeric_limits<Value>::max());
   4.489 +              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
   4.490 +                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
   4.491 +              }
   4.492 +
   4.493 +              if ((*_blossom_data)[vb].status == MATCHED) {
   4.494 +                if (_blossom_set->classPrio(vb) ==
   4.495 +                    std::numeric_limits<Value>::max()) {
   4.496 +                  _delta2->erase(vb);
   4.497 +                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
   4.498 +                           (*_blossom_data)[vb].offset) {
   4.499 +                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
   4.500 +                                   (*_blossom_data)[vb].offset);
   4.501 +                }
   4.502 +              }
   4.503              }
   4.504            }
   4.505          }
   4.506        }
   4.507      }
   4.508  
   4.509 -    void unmatchedToMatched(int blossom) {
   4.510 +    void oddToMatched(int blossom) {
   4.511 +      (*_blossom_data)[blossom].offset -= _delta_sum;
   4.512 +
   4.513 +      if (_blossom_set->classPrio(blossom) !=
   4.514 +          std::numeric_limits<Value>::max()) {
   4.515 +        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   4.516 +                      (*_blossom_data)[blossom].offset);
   4.517 +      }
   4.518 +
   4.519 +      if (!_blossom_set->trivial(blossom)) {
   4.520 +        _delta4->erase(blossom);
   4.521 +      }
   4.522 +    }
   4.523 +
   4.524 +    void oddToEven(int blossom, int tree) {
   4.525 +      if (!_blossom_set->trivial(blossom)) {
   4.526 +        _delta4->erase(blossom);
   4.527 +        (*_blossom_data)[blossom].pot -=
   4.528 +          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
   4.529 +      }
   4.530 +
   4.531        for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   4.532             n != INVALID; ++n) {
   4.533          int ni = (*_node_index)[n];
   4.534  
   4.535 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   4.536 +
   4.537 +        (*_node_data)[ni].heap.clear();
   4.538 +        (*_node_data)[ni].heap_index.clear();
   4.539 +        (*_node_data)[ni].pot +=
   4.540 +          2 * _delta_sum - (*_blossom_data)[blossom].offset;
   4.541 +
   4.542 +        _delta1->push(n, (*_node_data)[ni].pot);
   4.543 +
   4.544          for (InArcIt e(_graph, n); e != INVALID; ++e) {
   4.545            Node v = _graph.source(e);
   4.546            int vb = _blossom_set->find(v);
   4.547 @@ -1202,54 +1106,44 @@
   4.548            Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   4.549              dualScale * _weight[e];
   4.550  
   4.551 -          if (vb == blossom) {
   4.552 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.553 -              _delta3->erase(e);
   4.554 +          if ((*_blossom_data)[vb].status == EVEN) {
   4.555 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   4.556 +              _delta3->push(e, rw / 2);
   4.557              }
   4.558 -          } else if ((*_blossom_data)[vb].status == EVEN) {
   4.559 -
   4.560 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.561 -              _delta3->erase(e);
   4.562 -            }
   4.563 -
   4.564 -            int vt = _tree_set->find(vb);
   4.565 -
   4.566 -            Arc r = _graph.oppositeArc(e);
   4.567 +          } else {
   4.568  
   4.569              typename std::map<int, Arc>::iterator it =
   4.570 -              (*_node_data)[ni].heap_index.find(vt);
   4.571 -
   4.572 -            if (it != (*_node_data)[ni].heap_index.end()) {
   4.573 -              if ((*_node_data)[ni].heap[it->second] > rw) {
   4.574 -                (*_node_data)[ni].heap.replace(it->second, r);
   4.575 -                (*_node_data)[ni].heap.decrease(r, rw);
   4.576 -                it->second = r;
   4.577 +              (*_node_data)[vi].heap_index.find(tree);
   4.578 +
   4.579 +            if (it != (*_node_data)[vi].heap_index.end()) {
   4.580 +              if ((*_node_data)[vi].heap[it->second] > rw) {
   4.581 +                (*_node_data)[vi].heap.replace(it->second, e);
   4.582 +                (*_node_data)[vi].heap.decrease(e, rw);
   4.583 +                it->second = e;
   4.584                }
   4.585              } else {
   4.586 -              (*_node_data)[ni].heap.push(r, rw);
   4.587 -              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   4.588 +              (*_node_data)[vi].heap.push(e, rw);
   4.589 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   4.590              }
   4.591  
   4.592 -            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   4.593 -              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   4.594 -
   4.595 -              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   4.596 -                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   4.597 -                             (*_blossom_data)[blossom].offset);
   4.598 -              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
   4.599 -                         (*_blossom_data)[blossom].offset){
   4.600 -                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   4.601 -                                 (*_blossom_data)[blossom].offset);
   4.602 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   4.603 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   4.604 +
   4.605 +              if ((*_blossom_data)[vb].status == MATCHED) {
   4.606 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   4.607 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   4.608 +                               (*_blossom_data)[vb].offset);
   4.609 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   4.610 +                           (*_blossom_data)[vb].offset) {
   4.611 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   4.612 +                                   (*_blossom_data)[vb].offset);
   4.613 +                }
   4.614                }
   4.615              }
   4.616 -
   4.617 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   4.618 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   4.619 -              _delta3->erase(e);
   4.620 -            }
   4.621            }
   4.622          }
   4.623        }
   4.624 +      (*_blossom_data)[blossom].offset = 0;
   4.625      }
   4.626  
   4.627      void alternatePath(int even, int tree) {
   4.628 @@ -1294,39 +1188,42 @@
   4.629        alternatePath(blossom, tree);
   4.630        destroyTree(tree);
   4.631  
   4.632 -      (*_blossom_data)[blossom].status = UNMATCHED;
   4.633        (*_blossom_data)[blossom].base = node;
   4.634 -      matchedToUnmatched(blossom);
   4.635 +      (*_blossom_data)[blossom].next = INVALID;
   4.636      }
   4.637  
   4.638 -
   4.639      void augmentOnEdge(const Edge& edge) {
   4.640  
   4.641        int left = _blossom_set->find(_graph.u(edge));
   4.642        int right = _blossom_set->find(_graph.v(edge));
   4.643  
   4.644 -      if ((*_blossom_data)[left].status == EVEN) {
   4.645 -        int left_tree = _tree_set->find(left);
   4.646 -        alternatePath(left, left_tree);
   4.647 -        destroyTree(left_tree);
   4.648 -      } else {
   4.649 -        (*_blossom_data)[left].status = MATCHED;
   4.650 -        unmatchedToMatched(left);
   4.651 -      }
   4.652 -
   4.653 -      if ((*_blossom_data)[right].status == EVEN) {
   4.654 -        int right_tree = _tree_set->find(right);
   4.655 -        alternatePath(right, right_tree);
   4.656 -        destroyTree(right_tree);
   4.657 -      } else {
   4.658 -        (*_blossom_data)[right].status = MATCHED;
   4.659 -        unmatchedToMatched(right);
   4.660 -      }
   4.661 +      int left_tree = _tree_set->find(left);
   4.662 +      alternatePath(left, left_tree);
   4.663 +      destroyTree(left_tree);
   4.664 +
   4.665 +      int right_tree = _tree_set->find(right);
   4.666 +      alternatePath(right, right_tree);
   4.667 +      destroyTree(right_tree);
   4.668  
   4.669        (*_blossom_data)[left].next = _graph.direct(edge, true);
   4.670        (*_blossom_data)[right].next = _graph.direct(edge, false);
   4.671      }
   4.672  
   4.673 +    void augmentOnArc(const Arc& arc) {
   4.674 +
   4.675 +      int left = _blossom_set->find(_graph.source(arc));
   4.676 +      int right = _blossom_set->find(_graph.target(arc));
   4.677 +
   4.678 +      (*_blossom_data)[left].status = MATCHED;
   4.679 +
   4.680 +      int right_tree = _tree_set->find(right);
   4.681 +      alternatePath(right, right_tree);
   4.682 +      destroyTree(right_tree);
   4.683 +
   4.684 +      (*_blossom_data)[left].next = arc;
   4.685 +      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
   4.686 +    }
   4.687 +
   4.688      void extendOnArc(const Arc& arc) {
   4.689        int base = _blossom_set->find(_graph.target(arc));
   4.690        int tree = _tree_set->find(base);
   4.691 @@ -1529,7 +1426,7 @@
   4.692            _tree_set->insert(sb, tree);
   4.693            (*_blossom_data)[sb].pred = pred;
   4.694            (*_blossom_data)[sb].next =
   4.695 -                           _graph.oppositeArc((*_blossom_data)[tb].next);
   4.696 +            _graph.oppositeArc((*_blossom_data)[tb].next);
   4.697  
   4.698            pred = (*_blossom_data)[ub].next;
   4.699  
   4.700 @@ -1629,7 +1526,7 @@
   4.701        }
   4.702  
   4.703        for (int i = 0; i < int(blossoms.size()); ++i) {
   4.704 -        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
   4.705 +        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
   4.706  
   4.707            Value offset = (*_blossom_data)[blossoms[i]].offset;
   4.708            (*_blossom_data)[blossoms[i]].pot += 2 * offset;
   4.709 @@ -1667,10 +1564,16 @@
   4.710          _delta3_index(0), _delta3(0),
   4.711          _delta4_index(0), _delta4(0),
   4.712  
   4.713 -        _delta_sum() {}
   4.714 +        _delta_sum(), _unmatched(0),
   4.715 +
   4.716 +        _fractional(0)
   4.717 +    {}
   4.718  
   4.719      ~MaxWeightedMatching() {
   4.720        destroyStructures();
   4.721 +      if (_fractional) {
   4.722 +        delete _fractional;
   4.723 +      }
   4.724      }
   4.725  
   4.726      /// \name Execution Control
   4.727 @@ -1699,6 +1602,8 @@
   4.728          (*_delta4_index)[i] = _delta4->PRE_HEAP;
   4.729        }
   4.730  
   4.731 +      _unmatched = _node_num;
   4.732 +
   4.733        int index = 0;
   4.734        for (NodeIt n(_graph); n != INVALID; ++n) {
   4.735          Value max = 0;
   4.736 @@ -1733,18 +1638,155 @@
   4.737        }
   4.738      }
   4.739  
   4.740 +    /// \brief Initialize the algorithm with fractional matching
   4.741 +    ///
   4.742 +    /// This function initializes the algorithm with a fractional
   4.743 +    /// matching. This initialization is also called jumpstart heuristic.
   4.744 +    void fractionalInit() {
   4.745 +      createStructures();
   4.746 +      
   4.747 +      if (_fractional == 0) {
   4.748 +        _fractional = new FractionalMatching(_graph, _weight, false);
   4.749 +      }
   4.750 +      _fractional->run();
   4.751 +
   4.752 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   4.753 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
   4.754 +      }
   4.755 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   4.756 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
   4.757 +      }
   4.758 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   4.759 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
   4.760 +      }
   4.761 +      for (int i = 0; i < _blossom_num; ++i) {
   4.762 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
   4.763 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
   4.764 +      }
   4.765 +
   4.766 +      _unmatched = 0;
   4.767 +
   4.768 +      int index = 0;
   4.769 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   4.770 +        Value pot = _fractional->nodeValue(n);
   4.771 +        (*_node_index)[n] = index;
   4.772 +        (*_node_data)[index].pot = pot;
   4.773 +        int blossom =
   4.774 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
   4.775 +
   4.776 +        (*_blossom_data)[blossom].status = MATCHED;
   4.777 +        (*_blossom_data)[blossom].pred = INVALID;
   4.778 +        (*_blossom_data)[blossom].next = _fractional->matching(n);
   4.779 +        if (_fractional->matching(n) == INVALID) {
   4.780 +          (*_blossom_data)[blossom].base = n;
   4.781 +        }
   4.782 +        (*_blossom_data)[blossom].pot = 0;
   4.783 +        (*_blossom_data)[blossom].offset = 0;
   4.784 +        ++index;
   4.785 +      }
   4.786 +
   4.787 +      typename Graph::template NodeMap<bool> processed(_graph, false);
   4.788 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   4.789 +        if (processed[n]) continue;
   4.790 +        processed[n] = true;
   4.791 +        if (_fractional->matching(n) == INVALID) continue;
   4.792 +        int num = 1;
   4.793 +        Node v = _graph.target(_fractional->matching(n));
   4.794 +        while (n != v) {
   4.795 +          processed[v] = true;
   4.796 +          v = _graph.target(_fractional->matching(v));
   4.797 +          ++num;
   4.798 +        }
   4.799 +
   4.800 +        if (num % 2 == 1) {
   4.801 +          std::vector<int> subblossoms(num);
   4.802 +
   4.803 +          subblossoms[--num] = _blossom_set->find(n);
   4.804 +          _delta1->push(n, _fractional->nodeValue(n));
   4.805 +          v = _graph.target(_fractional->matching(n));
   4.806 +          while (n != v) {
   4.807 +            subblossoms[--num] = _blossom_set->find(v);
   4.808 +            _delta1->push(v, _fractional->nodeValue(v));
   4.809 +            v = _graph.target(_fractional->matching(v));            
   4.810 +          }
   4.811 +          
   4.812 +          int surface = 
   4.813 +            _blossom_set->join(subblossoms.begin(), subblossoms.end());
   4.814 +          (*_blossom_data)[surface].status = EVEN;
   4.815 +          (*_blossom_data)[surface].pred = INVALID;
   4.816 +          (*_blossom_data)[surface].next = INVALID;
   4.817 +          (*_blossom_data)[surface].pot = 0;
   4.818 +          (*_blossom_data)[surface].offset = 0;
   4.819 +          
   4.820 +          _tree_set->insert(surface);
   4.821 +          ++_unmatched;
   4.822 +        }
   4.823 +      }
   4.824 +
   4.825 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   4.826 +        int si = (*_node_index)[_graph.u(e)];
   4.827 +        int sb = _blossom_set->find(_graph.u(e));
   4.828 +        int ti = (*_node_index)[_graph.v(e)];
   4.829 +        int tb = _blossom_set->find(_graph.v(e));
   4.830 +        if ((*_blossom_data)[sb].status == EVEN &&
   4.831 +            (*_blossom_data)[tb].status == EVEN && sb != tb) {
   4.832 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
   4.833 +                            dualScale * _weight[e]) / 2);
   4.834 +        }
   4.835 +      }
   4.836 +
   4.837 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   4.838 +        int nb = _blossom_set->find(n);
   4.839 +        if ((*_blossom_data)[nb].status != MATCHED) continue;
   4.840 +        int ni = (*_node_index)[n];
   4.841 +
   4.842 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   4.843 +          Node v = _graph.target(e);
   4.844 +          int vb = _blossom_set->find(v);
   4.845 +          int vi = (*_node_index)[v];
   4.846 +
   4.847 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   4.848 +            dualScale * _weight[e];
   4.849 +
   4.850 +          if ((*_blossom_data)[vb].status == EVEN) {
   4.851 +
   4.852 +            int vt = _tree_set->find(vb);
   4.853 +
   4.854 +            typename std::map<int, Arc>::iterator it =
   4.855 +              (*_node_data)[ni].heap_index.find(vt);
   4.856 +
   4.857 +            if (it != (*_node_data)[ni].heap_index.end()) {
   4.858 +              if ((*_node_data)[ni].heap[it->second] > rw) {
   4.859 +                (*_node_data)[ni].heap.replace(it->second, e);
   4.860 +                (*_node_data)[ni].heap.decrease(e, rw);
   4.861 +                it->second = e;
   4.862 +              }
   4.863 +            } else {
   4.864 +              (*_node_data)[ni].heap.push(e, rw);
   4.865 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
   4.866 +            }
   4.867 +          }
   4.868 +        }
   4.869 +            
   4.870 +        if (!(*_node_data)[ni].heap.empty()) {
   4.871 +          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   4.872 +          _delta2->push(nb, _blossom_set->classPrio(nb));
   4.873 +        }
   4.874 +      }
   4.875 +    }
   4.876 +
   4.877      /// \brief Start the algorithm
   4.878      ///
   4.879      /// This function starts the algorithm.
   4.880      ///
   4.881 -    /// \pre \ref init() must be called before using this function.
   4.882 +    /// \pre \ref init() or \ref fractionalInit() must be called
   4.883 +    /// before using this function.
   4.884      void start() {
   4.885        enum OpType {
   4.886          D1, D2, D3, D4
   4.887        };
   4.888  
   4.889 -      int unmatched = _node_num;
   4.890 -      while (unmatched > 0) {
   4.891 +      while (_unmatched > 0) {
   4.892          Value d1 = !_delta1->empty() ?
   4.893            _delta1->prio() : std::numeric_limits<Value>::max();
   4.894  
   4.895 @@ -1757,26 +1799,30 @@
   4.896          Value d4 = !_delta4->empty() ?
   4.897            _delta4->prio() : std::numeric_limits<Value>::max();
   4.898  
   4.899 -        _delta_sum = d1; OpType ot = D1;
   4.900 +        _delta_sum = d3; OpType ot = D3;
   4.901 +        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
   4.902          if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
   4.903 -        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
   4.904          if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
   4.905  
   4.906 -
   4.907          switch (ot) {
   4.908          case D1:
   4.909            {
   4.910              Node n = _delta1->top();
   4.911              unmatchNode(n);
   4.912 -            --unmatched;
   4.913 +            --_unmatched;
   4.914            }
   4.915            break;
   4.916          case D2:
   4.917            {
   4.918              int blossom = _delta2->top();
   4.919              Node n = _blossom_set->classTop(blossom);
   4.920 -            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
   4.921 -            extendOnArc(e);
   4.922 +            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
   4.923 +            if ((*_blossom_data)[blossom].next == INVALID) {
   4.924 +              augmentOnArc(a);
   4.925 +              --_unmatched;
   4.926 +            } else {
   4.927 +              extendOnArc(a);
   4.928 +            }
   4.929            }
   4.930            break;
   4.931          case D3:
   4.932 @@ -1789,26 +1835,14 @@
   4.933              if (left_blossom == right_blossom) {
   4.934                _delta3->pop();
   4.935              } else {
   4.936 -              int left_tree;
   4.937 -              if ((*_blossom_data)[left_blossom].status == EVEN) {
   4.938 -                left_tree = _tree_set->find(left_blossom);
   4.939 -              } else {
   4.940 -                left_tree = -1;
   4.941 -                ++unmatched;
   4.942 -              }
   4.943 -              int right_tree;
   4.944 -              if ((*_blossom_data)[right_blossom].status == EVEN) {
   4.945 -                right_tree = _tree_set->find(right_blossom);
   4.946 -              } else {
   4.947 -                right_tree = -1;
   4.948 -                ++unmatched;
   4.949 -              }
   4.950 +              int left_tree = _tree_set->find(left_blossom);
   4.951 +              int right_tree = _tree_set->find(right_blossom);
   4.952  
   4.953                if (left_tree == right_tree) {
   4.954                  shrinkOnEdge(e, left_tree);
   4.955                } else {
   4.956                  augmentOnEdge(e);
   4.957 -                unmatched -= 2;
   4.958 +                _unmatched -= 2;
   4.959                }
   4.960              }
   4.961            } break;
   4.962 @@ -1826,18 +1860,18 @@
   4.963      ///
   4.964      /// \note mwm.run() is just a shortcut of the following code.
   4.965      /// \code
   4.966 -    ///   mwm.init();
   4.967 +    ///   mwm.fractionalInit();
   4.968      ///   mwm.start();
   4.969      /// \endcode
   4.970      void run() {
   4.971 -      init();
   4.972 +      fractionalInit();
   4.973        start();
   4.974      }
   4.975  
   4.976      /// @}
   4.977  
   4.978      /// \name Primal Solution
   4.979 -    /// Functions to get the primal solution, i.e. the maximum weighted 
   4.980 +    /// Functions to get the primal solution, i.e. the maximum weighted
   4.981      /// matching.\n
   4.982      /// Either \ref run() or \ref start() function should be called before
   4.983      /// using them.
   4.984 @@ -1856,7 +1890,7 @@
   4.985            sum += _weight[(*_matching)[n]];
   4.986          }
   4.987        }
   4.988 -      return sum /= 2;
   4.989 +      return sum / 2;
   4.990      }
   4.991  
   4.992      /// \brief Return the size (cardinality) of the matching.
   4.993 @@ -1876,7 +1910,7 @@
   4.994  
   4.995      /// \brief Return \c true if the given edge is in the matching.
   4.996      ///
   4.997 -    /// This function returns \c true if the given edge is in the found 
   4.998 +    /// This function returns \c true if the given edge is in the found
   4.999      /// matching.
  4.1000      ///
  4.1001      /// \pre Either run() or start() must be called before using this function.
  4.1002 @@ -1887,7 +1921,7 @@
  4.1003      /// \brief Return the matching arc (or edge) incident to the given node.
  4.1004      ///
  4.1005      /// This function returns the matching arc (or edge) incident to the
  4.1006 -    /// given node in the found matching or \c INVALID if the node is 
  4.1007 +    /// given node in the found matching or \c INVALID if the node is
  4.1008      /// not covered by the matching.
  4.1009      ///
  4.1010      /// \pre Either run() or start() must be called before using this function.
  4.1011 @@ -1905,7 +1939,7 @@
  4.1012  
  4.1013      /// \brief Return the mate of the given node.
  4.1014      ///
  4.1015 -    /// This function returns the mate of the given node in the found 
  4.1016 +    /// This function returns the mate of the given node in the found
  4.1017      /// matching or \c INVALID if the node is not covered by the matching.
  4.1018      ///
  4.1019      /// \pre Either run() or start() must be called before using this function.
  4.1020 @@ -1925,8 +1959,8 @@
  4.1021  
  4.1022      /// \brief Return the value of the dual solution.
  4.1023      ///
  4.1024 -    /// This function returns the value of the dual solution. 
  4.1025 -    /// It should be equal to the primal value scaled by \ref dualScale 
  4.1026 +    /// This function returns the value of the dual solution.
  4.1027 +    /// It should be equal to the primal value scaled by \ref dualScale
  4.1028      /// "dual scale".
  4.1029      ///
  4.1030      /// \pre Either run() or start() must be called before using this function.
  4.1031 @@ -1981,9 +2015,9 @@
  4.1032  
  4.1033      /// \brief Iterator for obtaining the nodes of a blossom.
  4.1034      ///
  4.1035 -    /// This class provides an iterator for obtaining the nodes of the 
  4.1036 +    /// This class provides an iterator for obtaining the nodes of the
  4.1037      /// given blossom. It lists a subset of the nodes.
  4.1038 -    /// Before using this iterator, you must allocate a 
  4.1039 +    /// Before using this iterator, you must allocate a
  4.1040      /// MaxWeightedMatching class and execute it.
  4.1041      class BlossomIt {
  4.1042      public:
  4.1043 @@ -1992,8 +2026,8 @@
  4.1044        ///
  4.1045        /// Constructor to get the nodes of the given variable.
  4.1046        ///
  4.1047 -      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
  4.1048 -      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
  4.1049 +      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
  4.1050 +      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
  4.1051        /// called before initializing this iterator.
  4.1052        BlossomIt(const MaxWeightedMatching& algorithm, int variable)
  4.1053          : _algorithm(&algorithm)
  4.1054 @@ -2046,8 +2080,8 @@
  4.1055    /// is based on extensive use of priority queues and provides
  4.1056    /// \f$O(nm\log n)\f$ time complexity.
  4.1057    ///
  4.1058 -  /// The maximum weighted perfect matching problem is to find a subset of 
  4.1059 -  /// the edges in an undirected graph with maximum overall weight for which 
  4.1060 +  /// The maximum weighted perfect matching problem is to find a subset of
  4.1061 +  /// the edges in an undirected graph with maximum overall weight for which
  4.1062    /// each node has exactly one incident edge.
  4.1063    /// It can be formulated with the following linear program.
  4.1064    /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  4.1065 @@ -2070,16 +2104,16 @@
  4.1066    /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  4.1067        \frac{\vert B \vert - 1}{2}z_B\f] */
  4.1068    ///
  4.1069 -  /// The algorithm can be executed with the run() function. 
  4.1070 +  /// The algorithm can be executed with the run() function.
  4.1071    /// After it the matching (the primal solution) and the dual solution
  4.1072 -  /// can be obtained using the query functions and the 
  4.1073 -  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
  4.1074 -  /// which is able to iterate on the nodes of a blossom. 
  4.1075 +  /// can be obtained using the query functions and the
  4.1076 +  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
  4.1077 +  /// which is able to iterate on the nodes of a blossom.
  4.1078    /// If the value type is integer, then the dual solution is multiplied
  4.1079    /// by \ref MaxWeightedMatching::dualScale "4".
  4.1080    ///
  4.1081    /// \tparam GR The undirected graph type the algorithm runs on.
  4.1082 -  /// \tparam WM The type edge weight map. The default type is 
  4.1083 +  /// \tparam WM The type edge weight map. The default type is
  4.1084    /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  4.1085  #ifdef DOXYGEN
  4.1086    template <typename GR, typename WM>
  4.1087 @@ -2190,6 +2224,11 @@
  4.1088      BinHeap<Value, IntIntMap> *_delta4;
  4.1089  
  4.1090      Value _delta_sum;
  4.1091 +    int _unmatched;
  4.1092 +
  4.1093 +    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 
  4.1094 +    FractionalMatching;
  4.1095 +    FractionalMatching *_fractional;
  4.1096  
  4.1097      void createStructures() {
  4.1098        _node_num = countNodes(_graph);
  4.1099 @@ -2233,9 +2272,6 @@
  4.1100      }
  4.1101  
  4.1102      void destroyStructures() {
  4.1103 -      _node_num = countNodes(_graph);
  4.1104 -      _blossom_num = _node_num * 3 / 2;
  4.1105 -
  4.1106        if (_matching) {
  4.1107          delete _matching;
  4.1108        }
  4.1109 @@ -2908,10 +2944,16 @@
  4.1110          _delta3_index(0), _delta3(0),
  4.1111          _delta4_index(0), _delta4(0),
  4.1112  
  4.1113 -        _delta_sum() {}
  4.1114 +        _delta_sum(), _unmatched(0),
  4.1115 +
  4.1116 +        _fractional(0)
  4.1117 +    {}
  4.1118  
  4.1119      ~MaxWeightedPerfectMatching() {
  4.1120        destroyStructures();
  4.1121 +      if (_fractional) {
  4.1122 +        delete _fractional;
  4.1123 +      }
  4.1124      }
  4.1125  
  4.1126      /// \name Execution Control
  4.1127 @@ -2937,6 +2979,8 @@
  4.1128          (*_delta4_index)[i] = _delta4->PRE_HEAP;
  4.1129        }
  4.1130  
  4.1131 +      _unmatched = _node_num;
  4.1132 +
  4.1133        int index = 0;
  4.1134        for (NodeIt n(_graph); n != INVALID; ++n) {
  4.1135          Value max = - std::numeric_limits<Value>::max();
  4.1136 @@ -2970,18 +3014,152 @@
  4.1137        }
  4.1138      }
  4.1139  
  4.1140 +    /// \brief Initialize the algorithm with fractional matching
  4.1141 +    ///
  4.1142 +    /// This function initializes the algorithm with a fractional
  4.1143 +    /// matching. This initialization is also called jumpstart heuristic.
  4.1144 +    void fractionalInit() {
  4.1145 +      createStructures();
  4.1146 +      
  4.1147 +      if (_fractional == 0) {
  4.1148 +        _fractional = new FractionalMatching(_graph, _weight, false);
  4.1149 +      }
  4.1150 +      if (!_fractional->run()) {
  4.1151 +        _unmatched = -1;
  4.1152 +        return;
  4.1153 +      }
  4.1154 +
  4.1155 +      for (ArcIt e(_graph); e != INVALID; ++e) {
  4.1156 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
  4.1157 +      }
  4.1158 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  4.1159 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  4.1160 +      }
  4.1161 +      for (int i = 0; i < _blossom_num; ++i) {
  4.1162 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
  4.1163 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
  4.1164 +      }
  4.1165 +
  4.1166 +      _unmatched = 0;
  4.1167 +
  4.1168 +      int index = 0;
  4.1169 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  4.1170 +        Value pot = _fractional->nodeValue(n);
  4.1171 +        (*_node_index)[n] = index;
  4.1172 +        (*_node_data)[index].pot = pot;
  4.1173 +        int blossom =
  4.1174 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
  4.1175 +
  4.1176 +        (*_blossom_data)[blossom].status = MATCHED;
  4.1177 +        (*_blossom_data)[blossom].pred = INVALID;
  4.1178 +        (*_blossom_data)[blossom].next = _fractional->matching(n);
  4.1179 +        (*_blossom_data)[blossom].pot = 0;
  4.1180 +        (*_blossom_data)[blossom].offset = 0;
  4.1181 +        ++index;
  4.1182 +      }
  4.1183 +
  4.1184 +      typename Graph::template NodeMap<bool> processed(_graph, false);
  4.1185 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  4.1186 +        if (processed[n]) continue;
  4.1187 +        processed[n] = true;
  4.1188 +        if (_fractional->matching(n) == INVALID) continue;
  4.1189 +        int num = 1;
  4.1190 +        Node v = _graph.target(_fractional->matching(n));
  4.1191 +        while (n != v) {
  4.1192 +          processed[v] = true;
  4.1193 +          v = _graph.target(_fractional->matching(v));
  4.1194 +          ++num;
  4.1195 +        }
  4.1196 +
  4.1197 +        if (num % 2 == 1) {
  4.1198 +          std::vector<int> subblossoms(num);
  4.1199 +
  4.1200 +          subblossoms[--num] = _blossom_set->find(n);
  4.1201 +          v = _graph.target(_fractional->matching(n));
  4.1202 +          while (n != v) {
  4.1203 +            subblossoms[--num] = _blossom_set->find(v);
  4.1204 +            v = _graph.target(_fractional->matching(v));            
  4.1205 +          }
  4.1206 +          
  4.1207 +          int surface = 
  4.1208 +            _blossom_set->join(subblossoms.begin(), subblossoms.end());
  4.1209 +          (*_blossom_data)[surface].status = EVEN;
  4.1210 +          (*_blossom_data)[surface].pred = INVALID;
  4.1211 +          (*_blossom_data)[surface].next = INVALID;
  4.1212 +          (*_blossom_data)[surface].pot = 0;
  4.1213 +          (*_blossom_data)[surface].offset = 0;
  4.1214 +          
  4.1215 +          _tree_set->insert(surface);
  4.1216 +          ++_unmatched;
  4.1217 +        }
  4.1218 +      }
  4.1219 +
  4.1220 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  4.1221 +        int si = (*_node_index)[_graph.u(e)];
  4.1222 +        int sb = _blossom_set->find(_graph.u(e));
  4.1223 +        int ti = (*_node_index)[_graph.v(e)];
  4.1224 +        int tb = _blossom_set->find(_graph.v(e));
  4.1225 +        if ((*_blossom_data)[sb].status == EVEN &&
  4.1226 +            (*_blossom_data)[tb].status == EVEN && sb != tb) {
  4.1227 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  4.1228 +                            dualScale * _weight[e]) / 2);
  4.1229 +        }
  4.1230 +      }
  4.1231 +
  4.1232 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  4.1233 +        int nb = _blossom_set->find(n);
  4.1234 +        if ((*_blossom_data)[nb].status != MATCHED) continue;
  4.1235 +        int ni = (*_node_index)[n];
  4.1236 +
  4.1237 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  4.1238 +          Node v = _graph.target(e);
  4.1239 +          int vb = _blossom_set->find(v);
  4.1240 +          int vi = (*_node_index)[v];
  4.1241 +
  4.1242 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  4.1243 +            dualScale * _weight[e];
  4.1244 +
  4.1245 +          if ((*_blossom_data)[vb].status == EVEN) {
  4.1246 +
  4.1247 +            int vt = _tree_set->find(vb);
  4.1248 +
  4.1249 +            typename std::map<int, Arc>::iterator it =
  4.1250 +              (*_node_data)[ni].heap_index.find(vt);
  4.1251 +
  4.1252 +            if (it != (*_node_data)[ni].heap_index.end()) {
  4.1253 +              if ((*_node_data)[ni].heap[it->second] > rw) {
  4.1254 +                (*_node_data)[ni].heap.replace(it->second, e);
  4.1255 +                (*_node_data)[ni].heap.decrease(e, rw);
  4.1256 +                it->second = e;
  4.1257 +              }
  4.1258 +            } else {
  4.1259 +              (*_node_data)[ni].heap.push(e, rw);
  4.1260 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
  4.1261 +            }
  4.1262 +          }
  4.1263 +        }
  4.1264 +            
  4.1265 +        if (!(*_node_data)[ni].heap.empty()) {
  4.1266 +          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  4.1267 +          _delta2->push(nb, _blossom_set->classPrio(nb));
  4.1268 +        }
  4.1269 +      }
  4.1270 +    }
  4.1271 +
  4.1272      /// \brief Start the algorithm
  4.1273      ///
  4.1274      /// This function starts the algorithm.
  4.1275      ///
  4.1276 -    /// \pre \ref init() must be called before using this function.
  4.1277 +    /// \pre \ref init() or \ref fractionalInit() must be called before
  4.1278 +    /// using this function.
  4.1279      bool start() {
  4.1280        enum OpType {
  4.1281          D2, D3, D4
  4.1282        };
  4.1283  
  4.1284 -      int unmatched = _node_num;
  4.1285 -      while (unmatched > 0) {
  4.1286 +      if (_unmatched == -1) return false;
  4.1287 +
  4.1288 +      while (_unmatched > 0) {
  4.1289          Value d2 = !_delta2->empty() ?
  4.1290            _delta2->prio() : std::numeric_limits<Value>::max();
  4.1291  
  4.1292 @@ -2991,8 +3169,8 @@
  4.1293          Value d4 = !_delta4->empty() ?
  4.1294            _delta4->prio() : std::numeric_limits<Value>::max();
  4.1295  
  4.1296 -        _delta_sum = d2; OpType ot = D2;
  4.1297 -        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  4.1298 +        _delta_sum = d3; OpType ot = D3;
  4.1299 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  4.1300          if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  4.1301  
  4.1302          if (_delta_sum == std::numeric_limits<Value>::max()) {
  4.1303 @@ -3025,7 +3203,7 @@
  4.1304                  shrinkOnEdge(e, left_tree);
  4.1305                } else {
  4.1306                  augmentOnEdge(e);
  4.1307 -                unmatched -= 2;
  4.1308 +                _unmatched -= 2;
  4.1309                }
  4.1310              }
  4.1311            } break;
  4.1312 @@ -3044,18 +3222,18 @@
  4.1313      ///
  4.1314      /// \note mwpm.run() is just a shortcut of the following code.
  4.1315      /// \code
  4.1316 -    ///   mwpm.init();
  4.1317 +    ///   mwpm.fractionalInit();
  4.1318      ///   mwpm.start();
  4.1319      /// \endcode
  4.1320      bool run() {
  4.1321 -      init();
  4.1322 +      fractionalInit();
  4.1323        return start();
  4.1324      }
  4.1325  
  4.1326      /// @}
  4.1327  
  4.1328      /// \name Primal Solution
  4.1329 -    /// Functions to get the primal solution, i.e. the maximum weighted 
  4.1330 +    /// Functions to get the primal solution, i.e. the maximum weighted
  4.1331      /// perfect matching.\n
  4.1332      /// Either \ref run() or \ref start() function should be called before
  4.1333      /// using them.
  4.1334 @@ -3074,12 +3252,12 @@
  4.1335            sum += _weight[(*_matching)[n]];
  4.1336          }
  4.1337        }
  4.1338 -      return sum /= 2;
  4.1339 +      return sum / 2;
  4.1340      }
  4.1341  
  4.1342      /// \brief Return \c true if the given edge is in the matching.
  4.1343      ///
  4.1344 -    /// This function returns \c true if the given edge is in the found 
  4.1345 +    /// This function returns \c true if the given edge is in the found
  4.1346      /// matching.
  4.1347      ///
  4.1348      /// \pre Either run() or start() must be called before using this function.
  4.1349 @@ -3090,7 +3268,7 @@
  4.1350      /// \brief Return the matching arc (or edge) incident to the given node.
  4.1351      ///
  4.1352      /// This function returns the matching arc (or edge) incident to the
  4.1353 -    /// given node in the found matching or \c INVALID if the node is 
  4.1354 +    /// given node in the found matching or \c INVALID if the node is
  4.1355      /// not covered by the matching.
  4.1356      ///
  4.1357      /// \pre Either run() or start() must be called before using this function.
  4.1358 @@ -3108,7 +3286,7 @@
  4.1359  
  4.1360      /// \brief Return the mate of the given node.
  4.1361      ///
  4.1362 -    /// This function returns the mate of the given node in the found 
  4.1363 +    /// This function returns the mate of the given node in the found
  4.1364      /// matching or \c INVALID if the node is not covered by the matching.
  4.1365      ///
  4.1366      /// \pre Either run() or start() must be called before using this function.
  4.1367 @@ -3127,8 +3305,8 @@
  4.1368  
  4.1369      /// \brief Return the value of the dual solution.
  4.1370      ///
  4.1371 -    /// This function returns the value of the dual solution. 
  4.1372 -    /// It should be equal to the primal value scaled by \ref dualScale 
  4.1373 +    /// This function returns the value of the dual solution.
  4.1374 +    /// It should be equal to the primal value scaled by \ref dualScale
  4.1375      /// "dual scale".
  4.1376      ///
  4.1377      /// \pre Either run() or start() must be called before using this function.
  4.1378 @@ -3183,9 +3361,9 @@
  4.1379  
  4.1380      /// \brief Iterator for obtaining the nodes of a blossom.
  4.1381      ///
  4.1382 -    /// This class provides an iterator for obtaining the nodes of the 
  4.1383 +    /// This class provides an iterator for obtaining the nodes of the
  4.1384      /// given blossom. It lists a subset of the nodes.
  4.1385 -    /// Before using this iterator, you must allocate a 
  4.1386 +    /// Before using this iterator, you must allocate a
  4.1387      /// MaxWeightedPerfectMatching class and execute it.
  4.1388      class BlossomIt {
  4.1389      public:
  4.1390 @@ -3194,8 +3372,8 @@
  4.1391        ///
  4.1392        /// Constructor to get the nodes of the given variable.
  4.1393        ///
  4.1394 -      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
  4.1395 -      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
  4.1396 +      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
  4.1397 +      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
  4.1398        /// must be called before initializing this iterator.
  4.1399        BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
  4.1400          : _algorithm(&algorithm)
  4.1401 @@ -3241,4 +3419,4 @@
  4.1402  
  4.1403  } //END OF NAMESPACE LEMON
  4.1404  
  4.1405 -#endif //LEMON_MAX_MATCHING_H
  4.1406 +#endif //LEMON_MATCHING_H
     5.1 --- a/test/CMakeLists.txt	Tue Mar 16 21:18:39 2010 +0100
     5.2 +++ b/test/CMakeLists.txt	Tue Mar 16 21:27:35 2010 +0100
     5.3 @@ -21,6 +21,7 @@
     5.4    edge_set_test
     5.5    error_test
     5.6    euler_test
     5.7 +  fractional_matching_test
     5.8    gomory_hu_test
     5.9    graph_copy_test
    5.10    graph_test
     6.1 --- a/test/Makefile.am	Tue Mar 16 21:18:39 2010 +0100
     6.2 +++ b/test/Makefile.am	Tue Mar 16 21:27:35 2010 +0100
     6.3 @@ -23,6 +23,7 @@
     6.4  	test/edge_set_test \
     6.5  	test/error_test \
     6.6  	test/euler_test \
     6.7 +	test/fractional_matching_test \
     6.8  	test/gomory_hu_test \
     6.9  	test/graph_copy_test \
    6.10  	test/graph_test \
    6.11 @@ -71,6 +72,7 @@
    6.12  test_edge_set_test_SOURCES = test/edge_set_test.cc
    6.13  test_error_test_SOURCES = test/error_test.cc
    6.14  test_euler_test_SOURCES = test/euler_test.cc
    6.15 +test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
    6.16  test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
    6.17  test_graph_copy_test_SOURCES = test/graph_copy_test.cc
    6.18  test_graph_test_SOURCES = test/graph_test.cc
     7.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     7.2 +++ b/test/fractional_matching_test.cc	Tue Mar 16 21:27:35 2010 +0100
     7.3 @@ -0,0 +1,525 @@
     7.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     7.5 + *
     7.6 + * This file is a part of LEMON, a generic C++ optimization library.
     7.7 + *
     7.8 + * Copyright (C) 2003-2009
     7.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    7.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    7.11 + *
    7.12 + * Permission to use, modify and distribute this software is granted
    7.13 + * provided that this copyright notice appears in all copies. For
    7.14 + * precise terms see the accompanying LICENSE file.
    7.15 + *
    7.16 + * This software is provided "AS IS" with no warranty of any kind,
    7.17 + * express or implied, and with no claim as to its suitability for any
    7.18 + * purpose.
    7.19 + *
    7.20 + */
    7.21 +
    7.22 +#include <iostream>
    7.23 +#include <sstream>
    7.24 +#include <vector>
    7.25 +#include <queue>
    7.26 +#include <cstdlib>
    7.27 +
    7.28 +#include <lemon/fractional_matching.h>
    7.29 +#include <lemon/smart_graph.h>
    7.30 +#include <lemon/concepts/graph.h>
    7.31 +#include <lemon/concepts/maps.h>
    7.32 +#include <lemon/lgf_reader.h>
    7.33 +#include <lemon/math.h>
    7.34 +
    7.35 +#include "test_tools.h"
    7.36 +
    7.37 +using namespace std;
    7.38 +using namespace lemon;
    7.39 +
    7.40 +GRAPH_TYPEDEFS(SmartGraph);
    7.41 +
    7.42 +
    7.43 +const int lgfn = 4;
    7.44 +const std::string lgf[lgfn] = {
    7.45 +  "@nodes\n"
    7.46 +  "label\n"
    7.47 +  "0\n"
    7.48 +  "1\n"
    7.49 +  "2\n"
    7.50 +  "3\n"
    7.51 +  "4\n"
    7.52 +  "5\n"
    7.53 +  "6\n"
    7.54 +  "7\n"
    7.55 +  "@edges\n"
    7.56 +  "     label  weight\n"
    7.57 +  "7 4  0      984\n"
    7.58 +  "0 7  1      73\n"
    7.59 +  "7 1  2      204\n"
    7.60 +  "2 3  3      583\n"
    7.61 +  "2 7  4      565\n"
    7.62 +  "2 1  5      582\n"
    7.63 +  "0 4  6      551\n"
    7.64 +  "2 5  7      385\n"
    7.65 +  "1 5  8      561\n"
    7.66 +  "5 3  9      484\n"
    7.67 +  "7 5  10     904\n"
    7.68 +  "3 6  11     47\n"
    7.69 +  "7 6  12     888\n"
    7.70 +  "3 0  13     747\n"
    7.71 +  "6 1  14     310\n",
    7.72 +
    7.73 +  "@nodes\n"
    7.74 +  "label\n"
    7.75 +  "0\n"
    7.76 +  "1\n"
    7.77 +  "2\n"
    7.78 +  "3\n"
    7.79 +  "4\n"
    7.80 +  "5\n"
    7.81 +  "6\n"
    7.82 +  "7\n"
    7.83 +  "@edges\n"
    7.84 +  "     label  weight\n"
    7.85 +  "2 5  0      710\n"
    7.86 +  "0 5  1      241\n"
    7.87 +  "2 4  2      856\n"
    7.88 +  "2 6  3      762\n"
    7.89 +  "4 1  4      747\n"
    7.90 +  "6 1  5      962\n"
    7.91 +  "4 7  6      723\n"
    7.92 +  "1 7  7      661\n"
    7.93 +  "2 3  8      376\n"
    7.94 +  "1 0  9      416\n"
    7.95 +  "6 7  10     391\n",
    7.96 +
    7.97 +  "@nodes\n"
    7.98 +  "label\n"
    7.99 +  "0\n"
   7.100 +  "1\n"
   7.101 +  "2\n"
   7.102 +  "3\n"
   7.103 +  "4\n"
   7.104 +  "5\n"
   7.105 +  "6\n"
   7.106 +  "7\n"
   7.107 +  "@edges\n"
   7.108 +  "     label  weight\n"
   7.109 +  "6 2  0      553\n"
   7.110 +  "0 7  1      653\n"
   7.111 +  "6 3  2      22\n"
   7.112 +  "4 7  3      846\n"
   7.113 +  "7 2  4      981\n"
   7.114 +  "7 6  5      250\n"
   7.115 +  "5 2  6      539\n",
   7.116 +
   7.117 +  "@nodes\n"
   7.118 +  "label\n"
   7.119 +  "0\n"
   7.120 +  "@edges\n"
   7.121 +  "     label  weight\n"
   7.122 +  "0 0  0      100\n"
   7.123 +};
   7.124 +
   7.125 +void checkMaxFractionalMatchingCompile()
   7.126 +{
   7.127 +  typedef concepts::Graph Graph;
   7.128 +  typedef Graph::Node Node;
   7.129 +  typedef Graph::Edge Edge;
   7.130 +
   7.131 +  Graph g;
   7.132 +  Node n;
   7.133 +  Edge e;
   7.134 +
   7.135 +  MaxFractionalMatching<Graph> mat_test(g);
   7.136 +  const MaxFractionalMatching<Graph>&
   7.137 +    const_mat_test = mat_test;
   7.138 +
   7.139 +  mat_test.init();
   7.140 +  mat_test.start();
   7.141 +  mat_test.start(true);
   7.142 +  mat_test.startPerfect();
   7.143 +  mat_test.startPerfect(true);
   7.144 +  mat_test.run();
   7.145 +  mat_test.run(true);
   7.146 +  mat_test.runPerfect();
   7.147 +  mat_test.runPerfect(true);
   7.148 +
   7.149 +  const_mat_test.matchingSize();
   7.150 +  const_mat_test.matching(e);
   7.151 +  const_mat_test.matching(n);
   7.152 +  const MaxFractionalMatching<Graph>::MatchingMap& mmap =
   7.153 +    const_mat_test.matchingMap();
   7.154 +  e = mmap[n];
   7.155 +
   7.156 +  const_mat_test.barrier(n);
   7.157 +}
   7.158 +
   7.159 +void checkMaxWeightedFractionalMatchingCompile()
   7.160 +{
   7.161 +  typedef concepts::Graph Graph;
   7.162 +  typedef Graph::Node Node;
   7.163 +  typedef Graph::Edge Edge;
   7.164 +  typedef Graph::EdgeMap<int> WeightMap;
   7.165 +
   7.166 +  Graph g;
   7.167 +  Node n;
   7.168 +  Edge e;
   7.169 +  WeightMap w(g);
   7.170 +
   7.171 +  MaxWeightedFractionalMatching<Graph> mat_test(g, w);
   7.172 +  const MaxWeightedFractionalMatching<Graph>&
   7.173 +    const_mat_test = mat_test;
   7.174 +
   7.175 +  mat_test.init();
   7.176 +  mat_test.start();
   7.177 +  mat_test.run();
   7.178 +
   7.179 +  const_mat_test.matchingWeight();
   7.180 +  const_mat_test.matchingSize();
   7.181 +  const_mat_test.matching(e);
   7.182 +  const_mat_test.matching(n);
   7.183 +  const MaxWeightedFractionalMatching<Graph>::MatchingMap& mmap =
   7.184 +    const_mat_test.matchingMap();
   7.185 +  e = mmap[n];
   7.186 +
   7.187 +  const_mat_test.dualValue();
   7.188 +  const_mat_test.nodeValue(n);
   7.189 +}
   7.190 +
   7.191 +void checkMaxWeightedPerfectFractionalMatchingCompile()
   7.192 +{
   7.193 +  typedef concepts::Graph Graph;
   7.194 +  typedef Graph::Node Node;
   7.195 +  typedef Graph::Edge Edge;
   7.196 +  typedef Graph::EdgeMap<int> WeightMap;
   7.197 +
   7.198 +  Graph g;
   7.199 +  Node n;
   7.200 +  Edge e;
   7.201 +  WeightMap w(g);
   7.202 +
   7.203 +  MaxWeightedPerfectFractionalMatching<Graph> mat_test(g, w);
   7.204 +  const MaxWeightedPerfectFractionalMatching<Graph>&
   7.205 +    const_mat_test = mat_test;
   7.206 +
   7.207 +  mat_test.init();
   7.208 +  mat_test.start();
   7.209 +  mat_test.run();
   7.210 +
   7.211 +  const_mat_test.matchingWeight();
   7.212 +  const_mat_test.matching(e);
   7.213 +  const_mat_test.matching(n);
   7.214 +  const MaxWeightedPerfectFractionalMatching<Graph>::MatchingMap& mmap =
   7.215 +    const_mat_test.matchingMap();
   7.216 +  e = mmap[n];
   7.217 +
   7.218 +  const_mat_test.dualValue();
   7.219 +  const_mat_test.nodeValue(n);
   7.220 +}
   7.221 +
   7.222 +void checkFractionalMatching(const SmartGraph& graph,
   7.223 +                             const MaxFractionalMatching<SmartGraph>& mfm,
   7.224 +                             bool allow_loops = true) {
   7.225 +  int pv = 0;
   7.226 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.227 +    int indeg = 0;
   7.228 +    for (InArcIt a(graph, n); a != INVALID; ++a) {
   7.229 +      if (mfm.matching(graph.source(a)) == a) {
   7.230 +        ++indeg;
   7.231 +      }
   7.232 +    }
   7.233 +    if (mfm.matching(n) != INVALID) {
   7.234 +      check(indeg == 1, "Invalid matching");
   7.235 +      ++pv;
   7.236 +    } else {
   7.237 +      check(indeg == 0, "Invalid matching");
   7.238 +    }
   7.239 +  }
   7.240 +  check(pv == mfm.matchingSize(), "Wrong matching size");
   7.241 +
   7.242 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   7.243 +    check((e == mfm.matching(graph.u(e)) ? 1 : 0) +
   7.244 +          (e == mfm.matching(graph.v(e)) ? 1 : 0) == 
   7.245 +          mfm.matching(e), "Invalid matching");
   7.246 +  }
   7.247 +
   7.248 +  SmartGraph::NodeMap<bool> processed(graph, false);
   7.249 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.250 +    if (processed[n]) continue;
   7.251 +    processed[n] = true;
   7.252 +    if (mfm.matching(n) == INVALID) continue;
   7.253 +    int num = 1;
   7.254 +    Node v = graph.target(mfm.matching(n));
   7.255 +    while (v != n) {
   7.256 +      processed[v] = true;
   7.257 +      ++num;
   7.258 +      v = graph.target(mfm.matching(v));
   7.259 +    }
   7.260 +    check(num == 2 || num % 2 == 1, "Wrong cycle size");
   7.261 +    check(allow_loops || num != 1, "Wrong cycle size");
   7.262 +  }
   7.263 +
   7.264 +  int anum = 0, bnum = 0;
   7.265 +  SmartGraph::NodeMap<bool> neighbours(graph, false);
   7.266 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.267 +    if (!mfm.barrier(n)) continue;
   7.268 +    ++anum;
   7.269 +    for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
   7.270 +      Node u = graph.source(a);
   7.271 +      if (!allow_loops && u == n) continue;
   7.272 +      if (!neighbours[u]) {
   7.273 +        neighbours[u] = true;
   7.274 +        ++bnum;
   7.275 +      }
   7.276 +    }
   7.277 +  }
   7.278 +  check(anum - bnum + mfm.matchingSize() == countNodes(graph),
   7.279 +        "Wrong barrier");
   7.280 +}
   7.281 +
   7.282 +void checkPerfectFractionalMatching(const SmartGraph& graph,
   7.283 +                             const MaxFractionalMatching<SmartGraph>& mfm,
   7.284 +                             bool perfect, bool allow_loops = true) {
   7.285 +  if (perfect) {
   7.286 +    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.287 +      int indeg = 0;
   7.288 +      for (InArcIt a(graph, n); a != INVALID; ++a) {
   7.289 +        if (mfm.matching(graph.source(a)) == a) {
   7.290 +          ++indeg;
   7.291 +        }
   7.292 +      }
   7.293 +      check(mfm.matching(n) != INVALID, "Invalid matching");
   7.294 +      check(indeg == 1, "Invalid matching");
   7.295 +    }
   7.296 +    for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   7.297 +      check((e == mfm.matching(graph.u(e)) ? 1 : 0) +
   7.298 +            (e == mfm.matching(graph.v(e)) ? 1 : 0) == 
   7.299 +            mfm.matching(e), "Invalid matching");
   7.300 +    }
   7.301 +  } else {
   7.302 +    int anum = 0, bnum = 0;
   7.303 +    SmartGraph::NodeMap<bool> neighbours(graph, false);
   7.304 +    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.305 +      if (!mfm.barrier(n)) continue;
   7.306 +      ++anum;
   7.307 +      for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
   7.308 +        Node u = graph.source(a);
   7.309 +        if (!allow_loops && u == n) continue;
   7.310 +        if (!neighbours[u]) {
   7.311 +          neighbours[u] = true;
   7.312 +          ++bnum;
   7.313 +        }
   7.314 +      }
   7.315 +    }
   7.316 +    check(anum - bnum > 0, "Wrong barrier");
   7.317 +  }
   7.318 +}
   7.319 +
   7.320 +void checkWeightedFractionalMatching(const SmartGraph& graph,
   7.321 +                   const SmartGraph::EdgeMap<int>& weight,
   7.322 +                   const MaxWeightedFractionalMatching<SmartGraph>& mwfm,
   7.323 +                   bool allow_loops = true) {
   7.324 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   7.325 +    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
   7.326 +    int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e))
   7.327 +      - weight[e] * mwfm.dualScale;
   7.328 +
   7.329 +    check(rw >= 0, "Negative reduced weight");
   7.330 +    check(rw == 0 || !mwfm.matching(e),
   7.331 +          "Non-zero reduced weight on matching edge");
   7.332 +  }
   7.333 +
   7.334 +  int pv = 0;
   7.335 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.336 +    int indeg = 0;
   7.337 +    for (InArcIt a(graph, n); a != INVALID; ++a) {
   7.338 +      if (mwfm.matching(graph.source(a)) == a) {
   7.339 +        ++indeg;
   7.340 +      }
   7.341 +    }
   7.342 +    check(indeg <= 1, "Invalid matching");
   7.343 +    if (mwfm.matching(n) != INVALID) {
   7.344 +      check(mwfm.nodeValue(n) >= 0, "Invalid node value");
   7.345 +      check(indeg == 1, "Invalid matching");
   7.346 +      pv += weight[mwfm.matching(n)];
   7.347 +      SmartGraph::Node o = graph.target(mwfm.matching(n));
   7.348 +    } else {
   7.349 +      check(mwfm.nodeValue(n) == 0, "Invalid matching");
   7.350 +      check(indeg == 0, "Invalid matching");
   7.351 +    }
   7.352 +  }
   7.353 +
   7.354 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   7.355 +    check((e == mwfm.matching(graph.u(e)) ? 1 : 0) +
   7.356 +          (e == mwfm.matching(graph.v(e)) ? 1 : 0) == 
   7.357 +          mwfm.matching(e), "Invalid matching");
   7.358 +  }
   7.359 +
   7.360 +  int dv = 0;
   7.361 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.362 +    dv += mwfm.nodeValue(n);
   7.363 +  }
   7.364 +
   7.365 +  check(pv * mwfm.dualScale == dv * 2, "Wrong duality");
   7.366 +
   7.367 +  SmartGraph::NodeMap<bool> processed(graph, false);
   7.368 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.369 +    if (processed[n]) continue;
   7.370 +    processed[n] = true;
   7.371 +    if (mwfm.matching(n) == INVALID) continue;
   7.372 +    int num = 1;
   7.373 +    Node v = graph.target(mwfm.matching(n));
   7.374 +    while (v != n) {
   7.375 +      processed[v] = true;
   7.376 +      ++num;
   7.377 +      v = graph.target(mwfm.matching(v));
   7.378 +    }
   7.379 +    check(num == 2 || num % 2 == 1, "Wrong cycle size");
   7.380 +    check(allow_loops || num != 1, "Wrong cycle size");
   7.381 +  }
   7.382 +
   7.383 +  return;
   7.384 +}
   7.385 +
   7.386 +void checkWeightedPerfectFractionalMatching(const SmartGraph& graph,
   7.387 +                const SmartGraph::EdgeMap<int>& weight,
   7.388 +                const MaxWeightedPerfectFractionalMatching<SmartGraph>& mwpfm,
   7.389 +                bool allow_loops = true) {
   7.390 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   7.391 +    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
   7.392 +    int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e))
   7.393 +      - weight[e] * mwpfm.dualScale;
   7.394 +
   7.395 +    check(rw >= 0, "Negative reduced weight");
   7.396 +    check(rw == 0 || !mwpfm.matching(e),
   7.397 +          "Non-zero reduced weight on matching edge");
   7.398 +  }
   7.399 +
   7.400 +  int pv = 0;
   7.401 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.402 +    int indeg = 0;
   7.403 +    for (InArcIt a(graph, n); a != INVALID; ++a) {
   7.404 +      if (mwpfm.matching(graph.source(a)) == a) {
   7.405 +        ++indeg;
   7.406 +      }
   7.407 +    }
   7.408 +    check(mwpfm.matching(n) != INVALID, "Invalid perfect matching");
   7.409 +    check(indeg == 1, "Invalid perfect matching");
   7.410 +    pv += weight[mwpfm.matching(n)];
   7.411 +    SmartGraph::Node o = graph.target(mwpfm.matching(n));
   7.412 +  }
   7.413 +
   7.414 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   7.415 +    check((e == mwpfm.matching(graph.u(e)) ? 1 : 0) +
   7.416 +          (e == mwpfm.matching(graph.v(e)) ? 1 : 0) == 
   7.417 +          mwpfm.matching(e), "Invalid matching");
   7.418 +  }
   7.419 +
   7.420 +  int dv = 0;
   7.421 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.422 +    dv += mwpfm.nodeValue(n);
   7.423 +  }
   7.424 +
   7.425 +  check(pv * mwpfm.dualScale == dv * 2, "Wrong duality");
   7.426 +
   7.427 +  SmartGraph::NodeMap<bool> processed(graph, false);
   7.428 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   7.429 +    if (processed[n]) continue;
   7.430 +    processed[n] = true;
   7.431 +    if (mwpfm.matching(n) == INVALID) continue;
   7.432 +    int num = 1;
   7.433 +    Node v = graph.target(mwpfm.matching(n));
   7.434 +    while (v != n) {
   7.435 +      processed[v] = true;
   7.436 +      ++num;
   7.437 +      v = graph.target(mwpfm.matching(v));
   7.438 +    }
   7.439 +    check(num == 2 || num % 2 == 1, "Wrong cycle size");
   7.440 +    check(allow_loops || num != 1, "Wrong cycle size");
   7.441 +  }
   7.442 +
   7.443 +  return;
   7.444 +}
   7.445 +
   7.446 +
   7.447 +int main() {
   7.448 +
   7.449 +  for (int i = 0; i < lgfn; ++i) {
   7.450 +    SmartGraph graph;
   7.451 +    SmartGraph::EdgeMap<int> weight(graph);
   7.452 +
   7.453 +    istringstream lgfs(lgf[i]);
   7.454 +    graphReader(graph, lgfs).
   7.455 +      edgeMap("weight", weight).run();
   7.456 +
   7.457 +    bool perfect_with_loops;
   7.458 +    {
   7.459 +      MaxFractionalMatching<SmartGraph> mfm(graph, true);
   7.460 +      mfm.run();
   7.461 +      checkFractionalMatching(graph, mfm, true);
   7.462 +      perfect_with_loops = mfm.matchingSize() == countNodes(graph);
   7.463 +    }
   7.464 +
   7.465 +    bool perfect_without_loops;
   7.466 +    {
   7.467 +      MaxFractionalMatching<SmartGraph> mfm(graph, false);
   7.468 +      mfm.run();
   7.469 +      checkFractionalMatching(graph, mfm, false);
   7.470 +      perfect_without_loops = mfm.matchingSize() == countNodes(graph);
   7.471 +    }
   7.472 +
   7.473 +    {
   7.474 +      MaxFractionalMatching<SmartGraph> mfm(graph, true);
   7.475 +      bool result = mfm.runPerfect();
   7.476 +      checkPerfectFractionalMatching(graph, mfm, result, true);
   7.477 +      check(result == perfect_with_loops, "Wrong perfect matching");
   7.478 +    }
   7.479 +
   7.480 +    {
   7.481 +      MaxFractionalMatching<SmartGraph> mfm(graph, false);
   7.482 +      bool result = mfm.runPerfect();
   7.483 +      checkPerfectFractionalMatching(graph, mfm, result, false);
   7.484 +      check(result == perfect_without_loops, "Wrong perfect matching");
   7.485 +    }
   7.486 +
   7.487 +    {
   7.488 +      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, true);
   7.489 +      mwfm.run();
   7.490 +      checkWeightedFractionalMatching(graph, weight, mwfm, true);
   7.491 +    }
   7.492 +
   7.493 +    {
   7.494 +      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, false);
   7.495 +      mwfm.run();
   7.496 +      checkWeightedFractionalMatching(graph, weight, mwfm, false);
   7.497 +    }
   7.498 +
   7.499 +    {
   7.500 +      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
   7.501 +                                                             true);
   7.502 +      bool perfect = mwpfm.run();
   7.503 +      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
   7.504 +            "Perfect matching found");
   7.505 +      check(perfect == perfect_with_loops, "Wrong perfect matching");
   7.506 +
   7.507 +      if (perfect) {
   7.508 +        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true);
   7.509 +      }
   7.510 +    }
   7.511 +
   7.512 +    {
   7.513 +      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
   7.514 +                                                             false);
   7.515 +      bool perfect = mwpfm.run();
   7.516 +      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
   7.517 +            "Perfect matching found");
   7.518 +      check(perfect == perfect_without_loops, "Wrong perfect matching");
   7.519 +
   7.520 +      if (perfect) {
   7.521 +        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false);
   7.522 +      }
   7.523 +    }
   7.524 +
   7.525 +  }
   7.526 +
   7.527 +  return 0;
   7.528 +}
     8.1 --- a/test/matching_test.cc	Tue Mar 16 21:18:39 2010 +0100
     8.2 +++ b/test/matching_test.cc	Tue Mar 16 21:27:35 2010 +0100
     8.3 @@ -401,22 +401,46 @@
     8.4      graphReader(graph, lgfs).
     8.5        edgeMap("weight", weight).run();
     8.6  
     8.7 -    MaxMatching<SmartGraph> mm(graph);
     8.8 -    mm.run();
     8.9 -    checkMatching(graph, mm);
    8.10 +    bool perfect;
    8.11 +    {
    8.12 +      MaxMatching<SmartGraph> mm(graph);
    8.13 +      mm.run();
    8.14 +      checkMatching(graph, mm);
    8.15 +      perfect = 2 * mm.matchingSize() == countNodes(graph);
    8.16 +    }
    8.17  
    8.18 -    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    8.19 -    mwm.run();
    8.20 -    checkWeightedMatching(graph, weight, mwm);
    8.21 +    {
    8.22 +      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    8.23 +      mwm.run();
    8.24 +      checkWeightedMatching(graph, weight, mwm);
    8.25 +    }
    8.26  
    8.27 -    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    8.28 -    bool perfect = mwpm.run();
    8.29 +    {
    8.30 +      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    8.31 +      mwm.init();
    8.32 +      mwm.start();
    8.33 +      checkWeightedMatching(graph, weight, mwm);
    8.34 +    }
    8.35  
    8.36 -    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
    8.37 -          "Perfect matching found");
    8.38 +    {
    8.39 +      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    8.40 +      bool result = mwpm.run();
    8.41 +      
    8.42 +      check(result == perfect, "Perfect matching found");
    8.43 +      if (perfect) {
    8.44 +        checkWeightedPerfectMatching(graph, weight, mwpm);
    8.45 +      }
    8.46 +    }
    8.47  
    8.48 -    if (perfect) {
    8.49 -      checkWeightedPerfectMatching(graph, weight, mwpm);
    8.50 +    {
    8.51 +      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    8.52 +      mwpm.init();
    8.53 +      bool result = mwpm.start();
    8.54 +      
    8.55 +      check(result == perfect, "Perfect matching found");
    8.56 +      if (perfect) {
    8.57 +        checkWeightedPerfectMatching(graph, weight, mwpm);
    8.58 +      }
    8.59      }
    8.60    }
    8.61