lemon/fractional_matching.h
changeset 869 636dadefe1e6
child 871 86613aa28a0c
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/fractional_matching.h	Fri Sep 25 21:51:36 2009 +0200
     1.3 @@ -0,0 +1,2135 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     1.8 + * Copyright (C) 2003-2009
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +#ifndef LEMON_FRACTIONAL_MATCHING_H
    1.23 +#define LEMON_FRACTIONAL_MATCHING_H
    1.24 +
    1.25 +#include <vector>
    1.26 +#include <queue>
    1.27 +#include <set>
    1.28 +#include <limits>
    1.29 +
    1.30 +#include <lemon/core.h>
    1.31 +#include <lemon/unionfind.h>
    1.32 +#include <lemon/bin_heap.h>
    1.33 +#include <lemon/maps.h>
    1.34 +#include <lemon/assert.h>
    1.35 +#include <lemon/elevator.h>
    1.36 +
    1.37 +///\ingroup matching
    1.38 +///\file
    1.39 +///\brief Fractional matching algorithms in general graphs.
    1.40 +
    1.41 +namespace lemon {
    1.42 +
    1.43 +  /// \brief Default traits class of MaxFractionalMatching class.
    1.44 +  ///
    1.45 +  /// Default traits class of MaxFractionalMatching class.
    1.46 +  /// \tparam GR Graph type.
    1.47 +  template <typename GR>
    1.48 +  struct MaxFractionalMatchingDefaultTraits {
    1.49 +
    1.50 +    /// \brief The type of the graph the algorithm runs on.
    1.51 +    typedef GR Graph;
    1.52 +
    1.53 +    /// \brief The type of the map that stores the matching.
    1.54 +    ///
    1.55 +    /// The type of the map that stores the matching arcs.
    1.56 +    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
    1.57 +    typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
    1.58 +
    1.59 +    /// \brief Instantiates a MatchingMap.
    1.60 +    ///
    1.61 +    /// This function instantiates a \ref MatchingMap.
    1.62 +    /// \param graph The graph for which we would like to define
    1.63 +    /// the matching map.
    1.64 +    static MatchingMap* createMatchingMap(const Graph& graph) {
    1.65 +      return new MatchingMap(graph);
    1.66 +    }
    1.67 +
    1.68 +    /// \brief The elevator type used by MaxFractionalMatching algorithm.
    1.69 +    ///
    1.70 +    /// The elevator type used by MaxFractionalMatching algorithm.
    1.71 +    ///
    1.72 +    /// \sa Elevator
    1.73 +    /// \sa LinkedElevator
    1.74 +    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
    1.75 +
    1.76 +    /// \brief Instantiates an Elevator.
    1.77 +    ///
    1.78 +    /// This function instantiates an \ref Elevator.
    1.79 +    /// \param graph The graph for which we would like to define
    1.80 +    /// the elevator.
    1.81 +    /// \param max_level The maximum level of the elevator.
    1.82 +    static Elevator* createElevator(const Graph& graph, int max_level) {
    1.83 +      return new Elevator(graph, max_level);
    1.84 +    }
    1.85 +  };
    1.86 +
    1.87 +  /// \ingroup matching
    1.88 +  ///
    1.89 +  /// \brief Max cardinality fractional matching
    1.90 +  ///
    1.91 +  /// This class provides an implementation of fractional matching
    1.92 +  /// algorithm based on push-relabel principle.
    1.93 +  ///
    1.94 +  /// The maximum cardinality fractional matching is a relaxation of the
    1.95 +  /// maximum cardinality matching problem where the odd set constraints
    1.96 +  /// are omitted.
    1.97 +  /// It can be formulated with the following linear program.
    1.98 +  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
    1.99 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
   1.100 +  /// \f[\max \sum_{e\in E}x_e\f]
   1.101 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   1.102 +  /// \f$X\f$. The result can be represented as the union of a
   1.103 +  /// matching with one value edges and a set of odd length cycles
   1.104 +  /// with half value edges.
   1.105 +  ///
   1.106 +  /// The algorithm calculates an optimal fractional matching and a
   1.107 +  /// barrier. The number of adjacents of any node set minus the size
   1.108 +  /// of node set is a lower bound on the uncovered nodes in the
   1.109 +  /// graph. For maximum matching a barrier is computed which
   1.110 +  /// maximizes this difference.
   1.111 +  ///
   1.112 +  /// The algorithm can be executed with the run() function.  After it
   1.113 +  /// the matching (the primal solution) and the barrier (the dual
   1.114 +  /// solution) can be obtained using the query functions.
   1.115 +  ///
   1.116 +  /// The primal solution is multiplied by
   1.117 +  /// \ref MaxWeightedMatching::primalScale "2".
   1.118 +  ///
   1.119 +  /// \tparam GR The undirected graph type the algorithm runs on.
   1.120 +#ifdef DOXYGEN
   1.121 +  template <typename GR, typename TR>
   1.122 +#else
   1.123 +  template <typename GR,
   1.124 +            typename TR = MaxFractionalMatchingDefaultTraits<GR> >
   1.125 +#endif
   1.126 +  class MaxFractionalMatching {
   1.127 +  public:
   1.128 +
   1.129 +    /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
   1.130 +    /// class" of the algorithm.
   1.131 +    typedef TR Traits;
   1.132 +    /// The type of the graph the algorithm runs on.
   1.133 +    typedef typename TR::Graph Graph;
   1.134 +    /// The type of the matching map.
   1.135 +    typedef typename TR::MatchingMap MatchingMap;
   1.136 +    /// The type of the elevator.
   1.137 +    typedef typename TR::Elevator Elevator;
   1.138 +
   1.139 +    /// \brief Scaling factor for primal solution
   1.140 +    ///
   1.141 +    /// Scaling factor for primal solution.
   1.142 +    static const int primalScale = 2;
   1.143 +
   1.144 +  private:
   1.145 +
   1.146 +    const Graph &_graph;
   1.147 +    int _node_num;
   1.148 +    bool _allow_loops;
   1.149 +    int _empty_level;
   1.150 +
   1.151 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
   1.152 +
   1.153 +    bool _local_matching;
   1.154 +    MatchingMap *_matching;
   1.155 +
   1.156 +    bool _local_level;
   1.157 +    Elevator *_level;
   1.158 +
   1.159 +    typedef typename Graph::template NodeMap<int> InDegMap;
   1.160 +    InDegMap *_indeg;
   1.161 +
   1.162 +    void createStructures() {
   1.163 +      _node_num = countNodes(_graph);
   1.164 +
   1.165 +      if (!_matching) {
   1.166 +        _local_matching = true;
   1.167 +        _matching = Traits::createMatchingMap(_graph);
   1.168 +      }
   1.169 +      if (!_level) {
   1.170 +        _local_level = true;
   1.171 +        _level = Traits::createElevator(_graph, _node_num);
   1.172 +      }
   1.173 +      if (!_indeg) {
   1.174 +        _indeg = new InDegMap(_graph);
   1.175 +      }
   1.176 +    }
   1.177 +
   1.178 +    void destroyStructures() {
   1.179 +      if (_local_matching) {
   1.180 +        delete _matching;
   1.181 +      }
   1.182 +      if (_local_level) {
   1.183 +        delete _level;
   1.184 +      }
   1.185 +      if (_indeg) {
   1.186 +        delete _indeg;
   1.187 +      }
   1.188 +    }
   1.189 +
   1.190 +    void postprocessing() {
   1.191 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.192 +        if ((*_indeg)[n] != 0) continue;
   1.193 +        _indeg->set(n, -1);
   1.194 +        Node u = n;
   1.195 +        while ((*_matching)[u] != INVALID) {
   1.196 +          Node v = _graph.target((*_matching)[u]);
   1.197 +          _indeg->set(v, -1);
   1.198 +          Arc a = _graph.oppositeArc((*_matching)[u]);
   1.199 +          u = _graph.target((*_matching)[v]);
   1.200 +          _indeg->set(u, -1);
   1.201 +          _matching->set(v, a);
   1.202 +        }
   1.203 +      }
   1.204 +
   1.205 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.206 +        if ((*_indeg)[n] != 1) continue;
   1.207 +        _indeg->set(n, -1);
   1.208 +
   1.209 +        int num = 1;
   1.210 +        Node u = _graph.target((*_matching)[n]);
   1.211 +        while (u != n) {
   1.212 +          _indeg->set(u, -1);
   1.213 +          u = _graph.target((*_matching)[u]);
   1.214 +          ++num;
   1.215 +        }
   1.216 +        if (num % 2 == 0 && num > 2) {
   1.217 +          Arc prev = _graph.oppositeArc((*_matching)[n]);
   1.218 +          Node v = _graph.target((*_matching)[n]);
   1.219 +          u = _graph.target((*_matching)[v]);
   1.220 +          _matching->set(v, prev);
   1.221 +          while (u != n) {
   1.222 +            prev = _graph.oppositeArc((*_matching)[u]);
   1.223 +            v = _graph.target((*_matching)[u]);
   1.224 +            u = _graph.target((*_matching)[v]);
   1.225 +            _matching->set(v, prev);
   1.226 +          }
   1.227 +        }
   1.228 +      }
   1.229 +    }
   1.230 +
   1.231 +  public:
   1.232 +
   1.233 +    typedef MaxFractionalMatching Create;
   1.234 +
   1.235 +    ///\name Named Template Parameters
   1.236 +
   1.237 +    ///@{
   1.238 +
   1.239 +    template <typename T>
   1.240 +    struct SetMatchingMapTraits : public Traits {
   1.241 +      typedef T MatchingMap;
   1.242 +      static MatchingMap *createMatchingMap(const Graph&) {
   1.243 +        LEMON_ASSERT(false, "MatchingMap is not initialized");
   1.244 +        return 0; // ignore warnings
   1.245 +      }
   1.246 +    };
   1.247 +
   1.248 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.249 +    /// MatchingMap type
   1.250 +    ///
   1.251 +    /// \ref named-templ-param "Named parameter" for setting MatchingMap
   1.252 +    /// type.
   1.253 +    template <typename T>
   1.254 +    struct SetMatchingMap
   1.255 +      : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
   1.256 +      typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
   1.257 +    };
   1.258 +
   1.259 +    template <typename T>
   1.260 +    struct SetElevatorTraits : public Traits {
   1.261 +      typedef T Elevator;
   1.262 +      static Elevator *createElevator(const Graph&, int) {
   1.263 +        LEMON_ASSERT(false, "Elevator is not initialized");
   1.264 +        return 0; // ignore warnings
   1.265 +      }
   1.266 +    };
   1.267 +
   1.268 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.269 +    /// Elevator type
   1.270 +    ///
   1.271 +    /// \ref named-templ-param "Named parameter" for setting Elevator
   1.272 +    /// type. If this named parameter is used, then an external
   1.273 +    /// elevator object must be passed to the algorithm using the
   1.274 +    /// \ref elevator(Elevator&) "elevator()" function before calling
   1.275 +    /// \ref run() or \ref init().
   1.276 +    /// \sa SetStandardElevator
   1.277 +    template <typename T>
   1.278 +    struct SetElevator
   1.279 +      : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
   1.280 +      typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
   1.281 +    };
   1.282 +
   1.283 +    template <typename T>
   1.284 +    struct SetStandardElevatorTraits : public Traits {
   1.285 +      typedef T Elevator;
   1.286 +      static Elevator *createElevator(const Graph& graph, int max_level) {
   1.287 +        return new Elevator(graph, max_level);
   1.288 +      }
   1.289 +    };
   1.290 +
   1.291 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.292 +    /// Elevator type with automatic allocation
   1.293 +    ///
   1.294 +    /// \ref named-templ-param "Named parameter" for setting Elevator
   1.295 +    /// type with automatic allocation.
   1.296 +    /// The Elevator should have standard constructor interface to be
   1.297 +    /// able to automatically created by the algorithm (i.e. the
   1.298 +    /// graph and the maximum level should be passed to it).
   1.299 +    /// However an external elevator object could also be passed to the
   1.300 +    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
   1.301 +    /// before calling \ref run() or \ref init().
   1.302 +    /// \sa SetElevator
   1.303 +    template <typename T>
   1.304 +    struct SetStandardElevator
   1.305 +      : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
   1.306 +      typedef MaxFractionalMatching<Graph,
   1.307 +                                    SetStandardElevatorTraits<T> > Create;
   1.308 +    };
   1.309 +
   1.310 +    /// @}
   1.311 +
   1.312 +  protected:
   1.313 +
   1.314 +    MaxFractionalMatching() {}
   1.315 +
   1.316 +  public:
   1.317 +
   1.318 +    /// \brief Constructor
   1.319 +    ///
   1.320 +    /// Constructor.
   1.321 +    ///
   1.322 +    MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
   1.323 +      : _graph(graph), _allow_loops(allow_loops),
   1.324 +        _local_matching(false), _matching(0),
   1.325 +        _local_level(false), _level(0),  _indeg(0)
   1.326 +    {}
   1.327 +
   1.328 +    ~MaxFractionalMatching() {
   1.329 +      destroyStructures();
   1.330 +    }
   1.331 +
   1.332 +    /// \brief Sets the matching map.
   1.333 +    ///
   1.334 +    /// Sets the matching map.
   1.335 +    /// If you don't use this function before calling \ref run() or
   1.336 +    /// \ref init(), an instance will be allocated automatically.
   1.337 +    /// The destructor deallocates this automatically allocated map,
   1.338 +    /// of course.
   1.339 +    /// \return <tt>(*this)</tt>
   1.340 +    MaxFractionalMatching& matchingMap(MatchingMap& map) {
   1.341 +      if (_local_matching) {
   1.342 +        delete _matching;
   1.343 +        _local_matching = false;
   1.344 +      }
   1.345 +      _matching = &map;
   1.346 +      return *this;
   1.347 +    }
   1.348 +
   1.349 +    /// \brief Sets the elevator used by algorithm.
   1.350 +    ///
   1.351 +    /// Sets the elevator used by algorithm.
   1.352 +    /// If you don't use this function before calling \ref run() or
   1.353 +    /// \ref init(), an instance will be allocated automatically.
   1.354 +    /// The destructor deallocates this automatically allocated elevator,
   1.355 +    /// of course.
   1.356 +    /// \return <tt>(*this)</tt>
   1.357 +    MaxFractionalMatching& elevator(Elevator& elevator) {
   1.358 +      if (_local_level) {
   1.359 +        delete _level;
   1.360 +        _local_level = false;
   1.361 +      }
   1.362 +      _level = &elevator;
   1.363 +      return *this;
   1.364 +    }
   1.365 +
   1.366 +    /// \brief Returns a const reference to the elevator.
   1.367 +    ///
   1.368 +    /// Returns a const reference to the elevator.
   1.369 +    ///
   1.370 +    /// \pre Either \ref run() or \ref init() must be called before
   1.371 +    /// using this function.
   1.372 +    const Elevator& elevator() const {
   1.373 +      return *_level;
   1.374 +    }
   1.375 +
   1.376 +    /// \name Execution control
   1.377 +    /// The simplest way to execute the algorithm is to use one of the
   1.378 +    /// member functions called \c run(). \n
   1.379 +    /// If you need more control on the execution, first
   1.380 +    /// you must call \ref init() and then one variant of the start()
   1.381 +    /// member.
   1.382 +
   1.383 +    /// @{
   1.384 +
   1.385 +    /// \brief Initializes the internal data structures.
   1.386 +    ///
   1.387 +    /// Initializes the internal data structures and sets the initial
   1.388 +    /// matching.
   1.389 +    void init() {
   1.390 +      createStructures();
   1.391 +
   1.392 +      _level->initStart();
   1.393 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.394 +        _indeg->set(n, 0);
   1.395 +        _matching->set(n, INVALID);
   1.396 +        _level->initAddItem(n);
   1.397 +      }
   1.398 +      _level->initFinish();
   1.399 +
   1.400 +      _empty_level = _node_num;
   1.401 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.402 +        for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   1.403 +          if (_graph.target(a) == n && !_allow_loops) continue;
   1.404 +          _matching->set(n, a);
   1.405 +          Node v = _graph.target((*_matching)[n]);
   1.406 +          _indeg->set(v, (*_indeg)[v] + 1);
   1.407 +          break;
   1.408 +        }
   1.409 +      }
   1.410 +
   1.411 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.412 +        if ((*_indeg)[n] == 0) {
   1.413 +          _level->activate(n);
   1.414 +        }
   1.415 +      }
   1.416 +    }
   1.417 +
   1.418 +    /// \brief Starts the algorithm and computes a fractional matching
   1.419 +    ///
   1.420 +    /// The algorithm computes a maximum fractional matching.
   1.421 +    ///
   1.422 +    /// \param postprocess The algorithm computes first a matching
   1.423 +    /// which is a union of a matching with one value edges, cycles
   1.424 +    /// with half value edges and even length paths with half value
   1.425 +    /// edges. If the parameter is true, then after the push-relabel
   1.426 +    /// algorithm it postprocesses the matching to contain only
   1.427 +    /// matching edges and half value odd cycles.
   1.428 +    void start(bool postprocess = true) {
   1.429 +      Node n;
   1.430 +      while ((n = _level->highestActive()) != INVALID) {
   1.431 +        int level = _level->highestActiveLevel();
   1.432 +        int new_level = _level->maxLevel();
   1.433 +        for (InArcIt a(_graph, n); a != INVALID; ++a) {
   1.434 +          Node u = _graph.source(a);
   1.435 +          if (n == u && !_allow_loops) continue;
   1.436 +          Node v = _graph.target((*_matching)[u]);
   1.437 +          if ((*_level)[v] < level) {
   1.438 +            _indeg->set(v, (*_indeg)[v] - 1);
   1.439 +            if ((*_indeg)[v] == 0) {
   1.440 +              _level->activate(v);
   1.441 +            }
   1.442 +            _matching->set(u, a);
   1.443 +            _indeg->set(n, (*_indeg)[n] + 1);
   1.444 +            _level->deactivate(n);
   1.445 +            goto no_more_push;
   1.446 +          } else if (new_level > (*_level)[v]) {
   1.447 +            new_level = (*_level)[v];
   1.448 +          }
   1.449 +        }
   1.450 +
   1.451 +        if (new_level + 1 < _level->maxLevel()) {
   1.452 +          _level->liftHighestActive(new_level + 1);
   1.453 +        } else {
   1.454 +          _level->liftHighestActiveToTop();
   1.455 +        }
   1.456 +        if (_level->emptyLevel(level)) {
   1.457 +          _level->liftToTop(level);
   1.458 +        }
   1.459 +      no_more_push:
   1.460 +        ;
   1.461 +      }
   1.462 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.463 +        if ((*_matching)[n] == INVALID) continue;
   1.464 +        Node u = _graph.target((*_matching)[n]);
   1.465 +        if ((*_indeg)[u] > 1) {
   1.466 +          _indeg->set(u, (*_indeg)[u] - 1);
   1.467 +          _matching->set(n, INVALID);
   1.468 +        }
   1.469 +      }
   1.470 +      if (postprocess) {
   1.471 +        postprocessing();
   1.472 +      }
   1.473 +    }
   1.474 +
   1.475 +    /// \brief Starts the algorithm and computes a perfect fractional
   1.476 +    /// matching
   1.477 +    ///
   1.478 +    /// The algorithm computes a perfect fractional matching. If it
   1.479 +    /// does not exists, then the algorithm returns false and the
   1.480 +    /// matching is undefined and the barrier.
   1.481 +    ///
   1.482 +    /// \param postprocess The algorithm computes first a matching
   1.483 +    /// which is a union of a matching with one value edges, cycles
   1.484 +    /// with half value edges and even length paths with half value
   1.485 +    /// edges. If the parameter is true, then after the push-relabel
   1.486 +    /// algorithm it postprocesses the matching to contain only
   1.487 +    /// matching edges and half value odd cycles.
   1.488 +    bool startPerfect(bool postprocess = true) {
   1.489 +      Node n;
   1.490 +      while ((n = _level->highestActive()) != INVALID) {
   1.491 +        int level = _level->highestActiveLevel();
   1.492 +        int new_level = _level->maxLevel();
   1.493 +        for (InArcIt a(_graph, n); a != INVALID; ++a) {
   1.494 +          Node u = _graph.source(a);
   1.495 +          if (n == u && !_allow_loops) continue;
   1.496 +          Node v = _graph.target((*_matching)[u]);
   1.497 +          if ((*_level)[v] < level) {
   1.498 +            _indeg->set(v, (*_indeg)[v] - 1);
   1.499 +            if ((*_indeg)[v] == 0) {
   1.500 +              _level->activate(v);
   1.501 +            }
   1.502 +            _matching->set(u, a);
   1.503 +            _indeg->set(n, (*_indeg)[n] + 1);
   1.504 +            _level->deactivate(n);
   1.505 +            goto no_more_push;
   1.506 +          } else if (new_level > (*_level)[v]) {
   1.507 +            new_level = (*_level)[v];
   1.508 +          }
   1.509 +        }
   1.510 +
   1.511 +        if (new_level + 1 < _level->maxLevel()) {
   1.512 +          _level->liftHighestActive(new_level + 1);
   1.513 +        } else {
   1.514 +          _level->liftHighestActiveToTop();
   1.515 +          _empty_level = _level->maxLevel() - 1;
   1.516 +          return false;
   1.517 +        }
   1.518 +        if (_level->emptyLevel(level)) {
   1.519 +          _level->liftToTop(level);
   1.520 +          _empty_level = level;
   1.521 +          return false;
   1.522 +        }
   1.523 +      no_more_push:
   1.524 +        ;
   1.525 +      }
   1.526 +      if (postprocess) {
   1.527 +        postprocessing();
   1.528 +      }
   1.529 +      return true;
   1.530 +    }
   1.531 +
   1.532 +    /// \brief Runs the algorithm
   1.533 +    ///
   1.534 +    /// Just a shortcut for the next code:
   1.535 +    ///\code
   1.536 +    /// init();
   1.537 +    /// start();
   1.538 +    ///\endcode
   1.539 +    void run(bool postprocess = true) {
   1.540 +      init();
   1.541 +      start(postprocess);
   1.542 +    }
   1.543 +
   1.544 +    /// \brief Runs the algorithm to find a perfect fractional matching
   1.545 +    ///
   1.546 +    /// Just a shortcut for the next code:
   1.547 +    ///\code
   1.548 +    /// init();
   1.549 +    /// startPerfect();
   1.550 +    ///\endcode
   1.551 +    bool runPerfect(bool postprocess = true) {
   1.552 +      init();
   1.553 +      return startPerfect(postprocess);
   1.554 +    }
   1.555 +
   1.556 +    ///@}
   1.557 +
   1.558 +    /// \name Query Functions
   1.559 +    /// The result of the %Matching algorithm can be obtained using these
   1.560 +    /// functions.\n
   1.561 +    /// Before the use of these functions,
   1.562 +    /// either run() or start() must be called.
   1.563 +    ///@{
   1.564 +
   1.565 +
   1.566 +    /// \brief Return the number of covered nodes in the matching.
   1.567 +    ///
   1.568 +    /// This function returns the number of covered nodes in the matching.
   1.569 +    ///
   1.570 +    /// \pre Either run() or start() must be called before using this function.
   1.571 +    int matchingSize() const {
   1.572 +      int num = 0;
   1.573 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.574 +        if ((*_matching)[n] != INVALID) {
   1.575 +          ++num;
   1.576 +        }
   1.577 +      }
   1.578 +      return num;
   1.579 +    }
   1.580 +
   1.581 +    /// \brief Returns a const reference to the matching map.
   1.582 +    ///
   1.583 +    /// Returns a const reference to the node map storing the found
   1.584 +    /// fractional matching. This method can be called after
   1.585 +    /// running the algorithm.
   1.586 +    ///
   1.587 +    /// \pre Either \ref run() or \ref init() must be called before
   1.588 +    /// using this function.
   1.589 +    const MatchingMap& matchingMap() const {
   1.590 +      return *_matching;
   1.591 +    }
   1.592 +
   1.593 +    /// \brief Return \c true if the given edge is in the matching.
   1.594 +    ///
   1.595 +    /// This function returns \c true if the given edge is in the
   1.596 +    /// found matching. The result is scaled by \ref primalScale
   1.597 +    /// "primal scale".
   1.598 +    ///
   1.599 +    /// \pre Either run() or start() must be called before using this function.
   1.600 +    int matching(const Edge& edge) const {
   1.601 +      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
   1.602 +        (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
   1.603 +    }
   1.604 +
   1.605 +    /// \brief Return the fractional matching arc (or edge) incident
   1.606 +    /// to the given node.
   1.607 +    ///
   1.608 +    /// This function returns one of the fractional matching arc (or
   1.609 +    /// edge) incident to the given node in the found matching or \c
   1.610 +    /// INVALID if the node is not covered by the matching or if the
   1.611 +    /// node is on an odd length cycle then it is the successor edge
   1.612 +    /// on the cycle.
   1.613 +    ///
   1.614 +    /// \pre Either run() or start() must be called before using this function.
   1.615 +    Arc matching(const Node& node) const {
   1.616 +      return (*_matching)[node];
   1.617 +    }
   1.618 +
   1.619 +    /// \brief Returns true if the node is in the barrier
   1.620 +    ///
   1.621 +    /// The barrier is a subset of the nodes. If the nodes in the
   1.622 +    /// barrier have less adjacent nodes than the size of the barrier,
   1.623 +    /// then at least as much nodes cannot be covered as the
   1.624 +    /// difference of the two subsets.
   1.625 +    bool barrier(const Node& node) const {
   1.626 +      return (*_level)[node] >= _empty_level;
   1.627 +    }
   1.628 +
   1.629 +    /// @}
   1.630 +
   1.631 +  };
   1.632 +
   1.633 +  /// \ingroup matching
   1.634 +  ///
   1.635 +  /// \brief Weighted fractional matching in general graphs
   1.636 +  ///
   1.637 +  /// This class provides an efficient implementation of fractional
   1.638 +  /// matching algorithm. The implementation is based on extensive use
   1.639 +  /// of priority queues and provides \f$O(nm\log n)\f$ time
   1.640 +  /// complexity.
   1.641 +  ///
   1.642 +  /// The maximum weighted fractional matching is a relaxation of the
   1.643 +  /// maximum weighted matching problem where the odd set constraints
   1.644 +  /// are omitted.
   1.645 +  /// It can be formulated with the following linear program.
   1.646 +  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   1.647 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
   1.648 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
   1.649 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   1.650 +  /// \f$X\f$. The result must be the union of a matching with one
   1.651 +  /// value edges and a set of odd length cycles with half value edges.
   1.652 +  ///
   1.653 +  /// The algorithm calculates an optimal fractional matching and a
   1.654 +  /// proof of the optimality. The solution of the dual problem can be
   1.655 +  /// used to check the result of the algorithm. The dual linear
   1.656 +  /// problem is the following.
   1.657 +  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
   1.658 +  /// \f[y_u \ge 0 \quad \forall u \in V\f]
   1.659 +  /// \f[\min \sum_{u \in V}y_u \f] ///
   1.660 +  ///
   1.661 +  /// The algorithm can be executed with the run() function.
   1.662 +  /// After it the matching (the primal solution) and the dual solution
   1.663 +  /// can be obtained using the query functions.
   1.664 +  ///
   1.665 +  /// If the value type is integer, then the primal and the dual
   1.666 +  /// solutions are multiplied by
   1.667 +  /// \ref MaxWeightedMatching::primalScale "2" and
   1.668 +  /// \ref MaxWeightedMatching::dualScale "4" respectively.
   1.669 +  ///
   1.670 +  /// \tparam GR The undirected graph type the algorithm runs on.
   1.671 +  /// \tparam WM The type edge weight map. The default type is
   1.672 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   1.673 +#ifdef DOXYGEN
   1.674 +  template <typename GR, typename WM>
   1.675 +#else
   1.676 +  template <typename GR,
   1.677 +            typename WM = typename GR::template EdgeMap<int> >
   1.678 +#endif
   1.679 +  class MaxWeightedFractionalMatching {
   1.680 +  public:
   1.681 +
   1.682 +    /// The graph type of the algorithm
   1.683 +    typedef GR Graph;
   1.684 +    /// The type of the edge weight map
   1.685 +    typedef WM WeightMap;
   1.686 +    /// The value type of the edge weights
   1.687 +    typedef typename WeightMap::Value Value;
   1.688 +
   1.689 +    /// The type of the matching map
   1.690 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
   1.691 +    MatchingMap;
   1.692 +
   1.693 +    /// \brief Scaling factor for primal solution
   1.694 +    ///
   1.695 +    /// Scaling factor for primal solution. It is equal to 2 or 1
   1.696 +    /// according to the value type.
   1.697 +    static const int primalScale =
   1.698 +      std::numeric_limits<Value>::is_integer ? 2 : 1;
   1.699 +
   1.700 +    /// \brief Scaling factor for dual solution
   1.701 +    ///
   1.702 +    /// Scaling factor for dual solution. It is equal to 4 or 1
   1.703 +    /// according to the value type.
   1.704 +    static const int dualScale =
   1.705 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
   1.706 +
   1.707 +  private:
   1.708 +
   1.709 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
   1.710 +
   1.711 +    typedef typename Graph::template NodeMap<Value> NodePotential;
   1.712 +
   1.713 +    const Graph& _graph;
   1.714 +    const WeightMap& _weight;
   1.715 +
   1.716 +    MatchingMap* _matching;
   1.717 +    NodePotential* _node_potential;
   1.718 +
   1.719 +    int _node_num;
   1.720 +    bool _allow_loops;
   1.721 +
   1.722 +    enum Status {
   1.723 +      EVEN = -1, MATCHED = 0, ODD = 1
   1.724 +    };
   1.725 +
   1.726 +    typedef typename Graph::template NodeMap<Status> StatusMap;
   1.727 +    StatusMap* _status;
   1.728 +
   1.729 +    typedef typename Graph::template NodeMap<Arc> PredMap;
   1.730 +    PredMap* _pred;
   1.731 +
   1.732 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
   1.733 +
   1.734 +    IntNodeMap *_tree_set_index;
   1.735 +    TreeSet *_tree_set;
   1.736 +
   1.737 +    IntNodeMap *_delta1_index;
   1.738 +    BinHeap<Value, IntNodeMap> *_delta1;
   1.739 +
   1.740 +    IntNodeMap *_delta2_index;
   1.741 +    BinHeap<Value, IntNodeMap> *_delta2;
   1.742 +
   1.743 +    IntEdgeMap *_delta3_index;
   1.744 +    BinHeap<Value, IntEdgeMap> *_delta3;
   1.745 +
   1.746 +    Value _delta_sum;
   1.747 +
   1.748 +    void createStructures() {
   1.749 +      _node_num = countNodes(_graph);
   1.750 +
   1.751 +      if (!_matching) {
   1.752 +        _matching = new MatchingMap(_graph);
   1.753 +      }
   1.754 +      if (!_node_potential) {
   1.755 +        _node_potential = new NodePotential(_graph);
   1.756 +      }
   1.757 +      if (!_status) {
   1.758 +        _status = new StatusMap(_graph);
   1.759 +      }
   1.760 +      if (!_pred) {
   1.761 +        _pred = new PredMap(_graph);
   1.762 +      }
   1.763 +      if (!_tree_set) {
   1.764 +        _tree_set_index = new IntNodeMap(_graph);
   1.765 +        _tree_set = new TreeSet(*_tree_set_index);
   1.766 +      }
   1.767 +      if (!_delta1) {
   1.768 +        _delta1_index = new IntNodeMap(_graph);
   1.769 +        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   1.770 +      }
   1.771 +      if (!_delta2) {
   1.772 +        _delta2_index = new IntNodeMap(_graph);
   1.773 +        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
   1.774 +      }
   1.775 +      if (!_delta3) {
   1.776 +        _delta3_index = new IntEdgeMap(_graph);
   1.777 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   1.778 +      }
   1.779 +    }
   1.780 +
   1.781 +    void destroyStructures() {
   1.782 +      if (_matching) {
   1.783 +        delete _matching;
   1.784 +      }
   1.785 +      if (_node_potential) {
   1.786 +        delete _node_potential;
   1.787 +      }
   1.788 +      if (_status) {
   1.789 +        delete _status;
   1.790 +      }
   1.791 +      if (_pred) {
   1.792 +        delete _pred;
   1.793 +      }
   1.794 +      if (_tree_set) {
   1.795 +        delete _tree_set_index;
   1.796 +        delete _tree_set;
   1.797 +      }
   1.798 +      if (_delta1) {
   1.799 +        delete _delta1_index;
   1.800 +        delete _delta1;
   1.801 +      }
   1.802 +      if (_delta2) {
   1.803 +        delete _delta2_index;
   1.804 +        delete _delta2;
   1.805 +      }
   1.806 +      if (_delta3) {
   1.807 +        delete _delta3_index;
   1.808 +        delete _delta3;
   1.809 +      }
   1.810 +    }
   1.811 +
   1.812 +    void matchedToEven(Node node, int tree) {
   1.813 +      _tree_set->insert(node, tree);
   1.814 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   1.815 +      _delta1->push(node, (*_node_potential)[node]);
   1.816 +
   1.817 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
   1.818 +        _delta2->erase(node);
   1.819 +      }
   1.820 +
   1.821 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
   1.822 +        Node v = _graph.source(a);
   1.823 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   1.824 +          dualScale * _weight[a];
   1.825 +        if (node == v) {
   1.826 +          if (_allow_loops && _graph.direction(a)) {
   1.827 +            _delta3->push(a, rw / 2);
   1.828 +          }
   1.829 +        } else if ((*_status)[v] == EVEN) {
   1.830 +          _delta3->push(a, rw / 2);
   1.831 +        } else if ((*_status)[v] == MATCHED) {
   1.832 +          if (_delta2->state(v) != _delta2->IN_HEAP) {
   1.833 +            _pred->set(v, a);
   1.834 +            _delta2->push(v, rw);
   1.835 +          } else if ((*_delta2)[v] > rw) {
   1.836 +            _pred->set(v, a);
   1.837 +            _delta2->decrease(v, rw);
   1.838 +          }
   1.839 +        }
   1.840 +      }
   1.841 +    }
   1.842 +
   1.843 +    void matchedToOdd(Node node, int tree) {
   1.844 +      _tree_set->insert(node, tree);
   1.845 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   1.846 +
   1.847 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
   1.848 +        _delta2->erase(node);
   1.849 +      }
   1.850 +    }
   1.851 +
   1.852 +    void evenToMatched(Node node, int tree) {
   1.853 +      _delta1->erase(node);
   1.854 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   1.855 +      Arc min = INVALID;
   1.856 +      Value minrw = std::numeric_limits<Value>::max();
   1.857 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
   1.858 +        Node v = _graph.source(a);
   1.859 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   1.860 +          dualScale * _weight[a];
   1.861 +
   1.862 +        if (node == v) {
   1.863 +          if (_allow_loops && _graph.direction(a)) {
   1.864 +            _delta3->erase(a);
   1.865 +          }
   1.866 +        } else if ((*_status)[v] == EVEN) {
   1.867 +          _delta3->erase(a);
   1.868 +          if (minrw > rw) {
   1.869 +            min = _graph.oppositeArc(a);
   1.870 +            minrw = rw;
   1.871 +          }
   1.872 +        } else if ((*_status)[v]  == MATCHED) {
   1.873 +          if ((*_pred)[v] == a) {
   1.874 +            Arc mina = INVALID;
   1.875 +            Value minrwa = std::numeric_limits<Value>::max();
   1.876 +            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
   1.877 +              Node va = _graph.target(aa);
   1.878 +              if ((*_status)[va] != EVEN ||
   1.879 +                  _tree_set->find(va) == tree) continue;
   1.880 +              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
   1.881 +                dualScale * _weight[aa];
   1.882 +              if (minrwa > rwa) {
   1.883 +                minrwa = rwa;
   1.884 +                mina = aa;
   1.885 +              }
   1.886 +            }
   1.887 +            if (mina != INVALID) {
   1.888 +              _pred->set(v, mina);
   1.889 +              _delta2->increase(v, minrwa);
   1.890 +            } else {
   1.891 +              _pred->set(v, INVALID);
   1.892 +              _delta2->erase(v);
   1.893 +            }
   1.894 +          }
   1.895 +        }
   1.896 +      }
   1.897 +      if (min != INVALID) {
   1.898 +        _pred->set(node, min);
   1.899 +        _delta2->push(node, minrw);
   1.900 +      } else {
   1.901 +        _pred->set(node, INVALID);
   1.902 +      }
   1.903 +    }
   1.904 +
   1.905 +    void oddToMatched(Node node) {
   1.906 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   1.907 +      Arc min = INVALID;
   1.908 +      Value minrw = std::numeric_limits<Value>::max();
   1.909 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
   1.910 +        Node v = _graph.source(a);
   1.911 +        if ((*_status)[v] != EVEN) continue;
   1.912 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   1.913 +          dualScale * _weight[a];
   1.914 +
   1.915 +        if (minrw > rw) {
   1.916 +          min = _graph.oppositeArc(a);
   1.917 +          minrw = rw;
   1.918 +        }
   1.919 +      }
   1.920 +      if (min != INVALID) {
   1.921 +        _pred->set(node, min);
   1.922 +        _delta2->push(node, minrw);
   1.923 +      } else {
   1.924 +        _pred->set(node, INVALID);
   1.925 +      }
   1.926 +    }
   1.927 +
   1.928 +    void alternatePath(Node even, int tree) {
   1.929 +      Node odd;
   1.930 +
   1.931 +      _status->set(even, MATCHED);
   1.932 +      evenToMatched(even, tree);
   1.933 +
   1.934 +      Arc prev = (*_matching)[even];
   1.935 +      while (prev != INVALID) {
   1.936 +        odd = _graph.target(prev);
   1.937 +        even = _graph.target((*_pred)[odd]);
   1.938 +        _matching->set(odd, (*_pred)[odd]);
   1.939 +        _status->set(odd, MATCHED);
   1.940 +        oddToMatched(odd);
   1.941 +
   1.942 +        prev = (*_matching)[even];
   1.943 +        _status->set(even, MATCHED);
   1.944 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
   1.945 +        evenToMatched(even, tree);
   1.946 +      }
   1.947 +    }
   1.948 +
   1.949 +    void destroyTree(int tree) {
   1.950 +      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
   1.951 +        if ((*_status)[n] == EVEN) {
   1.952 +          _status->set(n, MATCHED);
   1.953 +          evenToMatched(n, tree);
   1.954 +        } else if ((*_status)[n] == ODD) {
   1.955 +          _status->set(n, MATCHED);
   1.956 +          oddToMatched(n);
   1.957 +        }
   1.958 +      }
   1.959 +      _tree_set->eraseClass(tree);
   1.960 +    }
   1.961 +
   1.962 +
   1.963 +    void unmatchNode(const Node& node) {
   1.964 +      int tree = _tree_set->find(node);
   1.965 +
   1.966 +      alternatePath(node, tree);
   1.967 +      destroyTree(tree);
   1.968 +
   1.969 +      _matching->set(node, INVALID);
   1.970 +    }
   1.971 +
   1.972 +
   1.973 +    void augmentOnEdge(const Edge& edge) {
   1.974 +      Node left = _graph.u(edge);
   1.975 +      int left_tree = _tree_set->find(left);
   1.976 +
   1.977 +      alternatePath(left, left_tree);
   1.978 +      destroyTree(left_tree);
   1.979 +      _matching->set(left, _graph.direct(edge, true));
   1.980 +
   1.981 +      Node right = _graph.v(edge);
   1.982 +      int right_tree = _tree_set->find(right);
   1.983 +
   1.984 +      alternatePath(right, right_tree);
   1.985 +      destroyTree(right_tree);
   1.986 +      _matching->set(right, _graph.direct(edge, false));
   1.987 +    }
   1.988 +
   1.989 +    void augmentOnArc(const Arc& arc) {
   1.990 +      Node left = _graph.source(arc);
   1.991 +      _status->set(left, MATCHED);
   1.992 +      _matching->set(left, arc);
   1.993 +      _pred->set(left, arc);
   1.994 +
   1.995 +      Node right = _graph.target(arc);
   1.996 +      int right_tree = _tree_set->find(right);
   1.997 +
   1.998 +      alternatePath(right, right_tree);
   1.999 +      destroyTree(right_tree);
  1.1000 +      _matching->set(right, _graph.oppositeArc(arc));
  1.1001 +    }
  1.1002 +
  1.1003 +    void extendOnArc(const Arc& arc) {
  1.1004 +      Node base = _graph.target(arc);
  1.1005 +      int tree = _tree_set->find(base);
  1.1006 +
  1.1007 +      Node odd = _graph.source(arc);
  1.1008 +      _tree_set->insert(odd, tree);
  1.1009 +      _status->set(odd, ODD);
  1.1010 +      matchedToOdd(odd, tree);
  1.1011 +      _pred->set(odd, arc);
  1.1012 +
  1.1013 +      Node even = _graph.target((*_matching)[odd]);
  1.1014 +      _tree_set->insert(even, tree);
  1.1015 +      _status->set(even, EVEN);
  1.1016 +      matchedToEven(even, tree);
  1.1017 +    }
  1.1018 +
  1.1019 +    void cycleOnEdge(const Edge& edge, int tree) {
  1.1020 +      Node nca = INVALID;
  1.1021 +      std::vector<Node> left_path, right_path;
  1.1022 +
  1.1023 +      {
  1.1024 +        std::set<Node> left_set, right_set;
  1.1025 +        Node left = _graph.u(edge);
  1.1026 +        left_path.push_back(left);
  1.1027 +        left_set.insert(left);
  1.1028 +
  1.1029 +        Node right = _graph.v(edge);
  1.1030 +        right_path.push_back(right);
  1.1031 +        right_set.insert(right);
  1.1032 +
  1.1033 +        while (true) {
  1.1034 +
  1.1035 +          if (left_set.find(right) != left_set.end()) {
  1.1036 +            nca = right;
  1.1037 +            break;
  1.1038 +          }
  1.1039 +
  1.1040 +          if ((*_matching)[left] == INVALID) break;
  1.1041 +
  1.1042 +          left = _graph.target((*_matching)[left]);
  1.1043 +          left_path.push_back(left);
  1.1044 +          left = _graph.target((*_pred)[left]);
  1.1045 +          left_path.push_back(left);
  1.1046 +
  1.1047 +          left_set.insert(left);
  1.1048 +
  1.1049 +          if (right_set.find(left) != right_set.end()) {
  1.1050 +            nca = left;
  1.1051 +            break;
  1.1052 +          }
  1.1053 +
  1.1054 +          if ((*_matching)[right] == INVALID) break;
  1.1055 +
  1.1056 +          right = _graph.target((*_matching)[right]);
  1.1057 +          right_path.push_back(right);
  1.1058 +          right = _graph.target((*_pred)[right]);
  1.1059 +          right_path.push_back(right);
  1.1060 +
  1.1061 +          right_set.insert(right);
  1.1062 +
  1.1063 +        }
  1.1064 +
  1.1065 +        if (nca == INVALID) {
  1.1066 +          if ((*_matching)[left] == INVALID) {
  1.1067 +            nca = right;
  1.1068 +            while (left_set.find(nca) == left_set.end()) {
  1.1069 +              nca = _graph.target((*_matching)[nca]);
  1.1070 +              right_path.push_back(nca);
  1.1071 +              nca = _graph.target((*_pred)[nca]);
  1.1072 +              right_path.push_back(nca);
  1.1073 +            }
  1.1074 +          } else {
  1.1075 +            nca = left;
  1.1076 +            while (right_set.find(nca) == right_set.end()) {
  1.1077 +              nca = _graph.target((*_matching)[nca]);
  1.1078 +              left_path.push_back(nca);
  1.1079 +              nca = _graph.target((*_pred)[nca]);
  1.1080 +              left_path.push_back(nca);
  1.1081 +            }
  1.1082 +          }
  1.1083 +        }
  1.1084 +      }
  1.1085 +
  1.1086 +      alternatePath(nca, tree);
  1.1087 +      Arc prev;
  1.1088 +
  1.1089 +      prev = _graph.direct(edge, true);
  1.1090 +      for (int i = 0; left_path[i] != nca; i += 2) {
  1.1091 +        _matching->set(left_path[i], prev);
  1.1092 +        _status->set(left_path[i], MATCHED);
  1.1093 +        evenToMatched(left_path[i], tree);
  1.1094 +
  1.1095 +        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  1.1096 +        _status->set(left_path[i + 1], MATCHED);
  1.1097 +        oddToMatched(left_path[i + 1]);
  1.1098 +      }
  1.1099 +      _matching->set(nca, prev);
  1.1100 +
  1.1101 +      for (int i = 0; right_path[i] != nca; i += 2) {
  1.1102 +        _status->set(right_path[i], MATCHED);
  1.1103 +        evenToMatched(right_path[i], tree);
  1.1104 +
  1.1105 +        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  1.1106 +        _status->set(right_path[i + 1], MATCHED);
  1.1107 +        oddToMatched(right_path[i + 1]);
  1.1108 +      }
  1.1109 +
  1.1110 +      destroyTree(tree);
  1.1111 +    }
  1.1112 +
  1.1113 +    void extractCycle(const Arc &arc) {
  1.1114 +      Node left = _graph.source(arc);
  1.1115 +      Node odd = _graph.target((*_matching)[left]);
  1.1116 +      Arc prev;
  1.1117 +      while (odd != left) {
  1.1118 +        Node even = _graph.target((*_matching)[odd]);
  1.1119 +        prev = (*_matching)[odd];
  1.1120 +        odd = _graph.target((*_matching)[even]);
  1.1121 +        _matching->set(even, _graph.oppositeArc(prev));
  1.1122 +      }
  1.1123 +      _matching->set(left, arc);
  1.1124 +
  1.1125 +      Node right = _graph.target(arc);
  1.1126 +      int right_tree = _tree_set->find(right);
  1.1127 +      alternatePath(right, right_tree);
  1.1128 +      destroyTree(right_tree);
  1.1129 +      _matching->set(right, _graph.oppositeArc(arc));
  1.1130 +    }
  1.1131 +
  1.1132 +  public:
  1.1133 +
  1.1134 +    /// \brief Constructor
  1.1135 +    ///
  1.1136 +    /// Constructor.
  1.1137 +    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
  1.1138 +                                  bool allow_loops = true)
  1.1139 +      : _graph(graph), _weight(weight), _matching(0),
  1.1140 +      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  1.1141 +      _status(0),  _pred(0),
  1.1142 +      _tree_set_index(0), _tree_set(0),
  1.1143 +
  1.1144 +      _delta1_index(0), _delta1(0),
  1.1145 +      _delta2_index(0), _delta2(0),
  1.1146 +      _delta3_index(0), _delta3(0),
  1.1147 +
  1.1148 +      _delta_sum() {}
  1.1149 +
  1.1150 +    ~MaxWeightedFractionalMatching() {
  1.1151 +      destroyStructures();
  1.1152 +    }
  1.1153 +
  1.1154 +    /// \name Execution Control
  1.1155 +    /// The simplest way to execute the algorithm is to use the
  1.1156 +    /// \ref run() member function.
  1.1157 +
  1.1158 +    ///@{
  1.1159 +
  1.1160 +    /// \brief Initialize the algorithm
  1.1161 +    ///
  1.1162 +    /// This function initializes the algorithm.
  1.1163 +    void init() {
  1.1164 +      createStructures();
  1.1165 +
  1.1166 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1167 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
  1.1168 +        (*_delta2_index)[n] = _delta2->PRE_HEAP;
  1.1169 +      }
  1.1170 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1171 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1.1172 +      }
  1.1173 +
  1.1174 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1175 +        Value max = 0;
  1.1176 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1.1177 +          if (_graph.target(e) == n && !_allow_loops) continue;
  1.1178 +          if ((dualScale * _weight[e]) / 2 > max) {
  1.1179 +            max = (dualScale * _weight[e]) / 2;
  1.1180 +          }
  1.1181 +        }
  1.1182 +        _node_potential->set(n, max);
  1.1183 +        _delta1->push(n, max);
  1.1184 +
  1.1185 +        _tree_set->insert(n);
  1.1186 +
  1.1187 +        _matching->set(n, INVALID);
  1.1188 +        _status->set(n, EVEN);
  1.1189 +      }
  1.1190 +
  1.1191 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1192 +        Node left = _graph.u(e);
  1.1193 +        Node right = _graph.v(e);
  1.1194 +        if (left == right && !_allow_loops) continue;
  1.1195 +        _delta3->push(e, ((*_node_potential)[left] +
  1.1196 +                          (*_node_potential)[right] -
  1.1197 +                          dualScale * _weight[e]) / 2);
  1.1198 +      }
  1.1199 +    }
  1.1200 +
  1.1201 +    /// \brief Start the algorithm
  1.1202 +    ///
  1.1203 +    /// This function starts the algorithm.
  1.1204 +    ///
  1.1205 +    /// \pre \ref init() must be called before using this function.
  1.1206 +    void start() {
  1.1207 +      enum OpType {
  1.1208 +        D1, D2, D3
  1.1209 +      };
  1.1210 +
  1.1211 +      int unmatched = _node_num;
  1.1212 +      while (unmatched > 0) {
  1.1213 +        Value d1 = !_delta1->empty() ?
  1.1214 +          _delta1->prio() : std::numeric_limits<Value>::max();
  1.1215 +
  1.1216 +        Value d2 = !_delta2->empty() ?
  1.1217 +          _delta2->prio() : std::numeric_limits<Value>::max();
  1.1218 +
  1.1219 +        Value d3 = !_delta3->empty() ?
  1.1220 +          _delta3->prio() : std::numeric_limits<Value>::max();
  1.1221 +
  1.1222 +        _delta_sum = d3; OpType ot = D3;
  1.1223 +        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
  1.1224 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1.1225 +
  1.1226 +        switch (ot) {
  1.1227 +        case D1:
  1.1228 +          {
  1.1229 +            Node n = _delta1->top();
  1.1230 +            unmatchNode(n);
  1.1231 +            --unmatched;
  1.1232 +          }
  1.1233 +          break;
  1.1234 +        case D2:
  1.1235 +          {
  1.1236 +            Node n = _delta2->top();
  1.1237 +            Arc a = (*_pred)[n];
  1.1238 +            if ((*_matching)[n] == INVALID) {
  1.1239 +              augmentOnArc(a);
  1.1240 +              --unmatched;
  1.1241 +            } else {
  1.1242 +              Node v = _graph.target((*_matching)[n]);
  1.1243 +              if ((*_matching)[n] !=
  1.1244 +                  _graph.oppositeArc((*_matching)[v])) {
  1.1245 +                extractCycle(a);
  1.1246 +                --unmatched;
  1.1247 +              } else {
  1.1248 +                extendOnArc(a);
  1.1249 +              }
  1.1250 +            }
  1.1251 +          } break;
  1.1252 +        case D3:
  1.1253 +          {
  1.1254 +            Edge e = _delta3->top();
  1.1255 +
  1.1256 +            Node left = _graph.u(e);
  1.1257 +            Node right = _graph.v(e);
  1.1258 +
  1.1259 +            int left_tree = _tree_set->find(left);
  1.1260 +            int right_tree = _tree_set->find(right);
  1.1261 +
  1.1262 +            if (left_tree == right_tree) {
  1.1263 +              cycleOnEdge(e, left_tree);
  1.1264 +              --unmatched;
  1.1265 +            } else {
  1.1266 +              augmentOnEdge(e);
  1.1267 +              unmatched -= 2;
  1.1268 +            }
  1.1269 +          } break;
  1.1270 +        }
  1.1271 +      }
  1.1272 +    }
  1.1273 +
  1.1274 +    /// \brief Run the algorithm.
  1.1275 +    ///
  1.1276 +    /// This method runs the \c %MaxWeightedMatching algorithm.
  1.1277 +    ///
  1.1278 +    /// \note mwfm.run() is just a shortcut of the following code.
  1.1279 +    /// \code
  1.1280 +    ///   mwfm.init();
  1.1281 +    ///   mwfm.start();
  1.1282 +    /// \endcode
  1.1283 +    void run() {
  1.1284 +      init();
  1.1285 +      start();
  1.1286 +    }
  1.1287 +
  1.1288 +    /// @}
  1.1289 +
  1.1290 +    /// \name Primal Solution
  1.1291 +    /// Functions to get the primal solution, i.e. the maximum weighted
  1.1292 +    /// matching.\n
  1.1293 +    /// Either \ref run() or \ref start() function should be called before
  1.1294 +    /// using them.
  1.1295 +
  1.1296 +    /// @{
  1.1297 +
  1.1298 +    /// \brief Return the weight of the matching.
  1.1299 +    ///
  1.1300 +    /// This function returns the weight of the found matching. This
  1.1301 +    /// value is scaled by \ref primalScale "primal scale".
  1.1302 +    ///
  1.1303 +    /// \pre Either run() or start() must be called before using this function.
  1.1304 +    Value matchingWeight() const {
  1.1305 +      Value sum = 0;
  1.1306 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1307 +        if ((*_matching)[n] != INVALID) {
  1.1308 +          sum += _weight[(*_matching)[n]];
  1.1309 +        }
  1.1310 +      }
  1.1311 +      return sum * primalScale / 2;
  1.1312 +    }
  1.1313 +
  1.1314 +    /// \brief Return the number of covered nodes in the matching.
  1.1315 +    ///
  1.1316 +    /// This function returns the number of covered nodes in the matching.
  1.1317 +    ///
  1.1318 +    /// \pre Either run() or start() must be called before using this function.
  1.1319 +    int matchingSize() const {
  1.1320 +      int num = 0;
  1.1321 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1322 +        if ((*_matching)[n] != INVALID) {
  1.1323 +          ++num;
  1.1324 +        }
  1.1325 +      }
  1.1326 +      return num;
  1.1327 +    }
  1.1328 +
  1.1329 +    /// \brief Return \c true if the given edge is in the matching.
  1.1330 +    ///
  1.1331 +    /// This function returns \c true if the given edge is in the
  1.1332 +    /// found matching. The result is scaled by \ref primalScale
  1.1333 +    /// "primal scale".
  1.1334 +    ///
  1.1335 +    /// \pre Either run() or start() must be called before using this function.
  1.1336 +    Value matching(const Edge& edge) const {
  1.1337 +      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  1.1338 +        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
  1.1339 +        * primalScale / 2;
  1.1340 +    }
  1.1341 +
  1.1342 +    /// \brief Return the fractional matching arc (or edge) incident
  1.1343 +    /// to the given node.
  1.1344 +    ///
  1.1345 +    /// This function returns one of the fractional matching arc (or
  1.1346 +    /// edge) incident to the given node in the found matching or \c
  1.1347 +    /// INVALID if the node is not covered by the matching or if the
  1.1348 +    /// node is on an odd length cycle then it is the successor edge
  1.1349 +    /// on the cycle.
  1.1350 +    ///
  1.1351 +    /// \pre Either run() or start() must be called before using this function.
  1.1352 +    Arc matching(const Node& node) const {
  1.1353 +      return (*_matching)[node];
  1.1354 +    }
  1.1355 +
  1.1356 +    /// \brief Return a const reference to the matching map.
  1.1357 +    ///
  1.1358 +    /// This function returns a const reference to a node map that stores
  1.1359 +    /// the matching arc (or edge) incident to each node.
  1.1360 +    const MatchingMap& matchingMap() const {
  1.1361 +      return *_matching;
  1.1362 +    }
  1.1363 +
  1.1364 +    /// @}
  1.1365 +
  1.1366 +    /// \name Dual Solution
  1.1367 +    /// Functions to get the dual solution.\n
  1.1368 +    /// Either \ref run() or \ref start() function should be called before
  1.1369 +    /// using them.
  1.1370 +
  1.1371 +    /// @{
  1.1372 +
  1.1373 +    /// \brief Return the value of the dual solution.
  1.1374 +    ///
  1.1375 +    /// This function returns the value of the dual solution.
  1.1376 +    /// It should be equal to the primal value scaled by \ref dualScale
  1.1377 +    /// "dual scale".
  1.1378 +    ///
  1.1379 +    /// \pre Either run() or start() must be called before using this function.
  1.1380 +    Value dualValue() const {
  1.1381 +      Value sum = 0;
  1.1382 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1383 +        sum += nodeValue(n);
  1.1384 +      }
  1.1385 +      return sum;
  1.1386 +    }
  1.1387 +
  1.1388 +    /// \brief Return the dual value (potential) of the given node.
  1.1389 +    ///
  1.1390 +    /// This function returns the dual value (potential) of the given node.
  1.1391 +    ///
  1.1392 +    /// \pre Either run() or start() must be called before using this function.
  1.1393 +    Value nodeValue(const Node& n) const {
  1.1394 +      return (*_node_potential)[n];
  1.1395 +    }
  1.1396 +
  1.1397 +    /// @}
  1.1398 +
  1.1399 +  };
  1.1400 +
  1.1401 +  /// \ingroup matching
  1.1402 +  ///
  1.1403 +  /// \brief Weighted fractional perfect matching in general graphs
  1.1404 +  ///
  1.1405 +  /// This class provides an efficient implementation of fractional
  1.1406 +  /// matching algorithm. The implementation is based on extensive use
  1.1407 +  /// of priority queues and provides \f$O(nm\log n)\f$ time
  1.1408 +  /// complexity.
  1.1409 +  ///
  1.1410 +  /// The maximum weighted fractional perfect matching is a relaxation
  1.1411 +  /// of the maximum weighted perfect matching problem where the odd
  1.1412 +  /// set constraints are omitted.
  1.1413 +  /// It can be formulated with the following linear program.
  1.1414 +  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1.1415 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
  1.1416 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
  1.1417 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1.1418 +  /// \f$X\f$. The result must be the union of a matching with one
  1.1419 +  /// value edges and a set of odd length cycles with half value edges.
  1.1420 +  ///
  1.1421 +  /// The algorithm calculates an optimal fractional matching and a
  1.1422 +  /// proof of the optimality. The solution of the dual problem can be
  1.1423 +  /// used to check the result of the algorithm. The dual linear
  1.1424 +  /// problem is the following.
  1.1425 +  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
  1.1426 +  /// \f[\min \sum_{u \in V}y_u \f] ///
  1.1427 +  ///
  1.1428 +  /// The algorithm can be executed with the run() function.
  1.1429 +  /// After it the matching (the primal solution) and the dual solution
  1.1430 +  /// can be obtained using the query functions.
  1.1431 +
  1.1432 +  /// If the value type is integer, then the primal and the dual
  1.1433 +  /// solutions are multiplied by
  1.1434 +  /// \ref MaxWeightedMatching::primalScale "2" and
  1.1435 +  /// \ref MaxWeightedMatching::dualScale "4" respectively.
  1.1436 +  ///
  1.1437 +  /// \tparam GR The undirected graph type the algorithm runs on.
  1.1438 +  /// \tparam WM The type edge weight map. The default type is
  1.1439 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  1.1440 +#ifdef DOXYGEN
  1.1441 +  template <typename GR, typename WM>
  1.1442 +#else
  1.1443 +  template <typename GR,
  1.1444 +            typename WM = typename GR::template EdgeMap<int> >
  1.1445 +#endif
  1.1446 +  class MaxWeightedPerfectFractionalMatching {
  1.1447 +  public:
  1.1448 +
  1.1449 +    /// The graph type of the algorithm
  1.1450 +    typedef GR Graph;
  1.1451 +    /// The type of the edge weight map
  1.1452 +    typedef WM WeightMap;
  1.1453 +    /// The value type of the edge weights
  1.1454 +    typedef typename WeightMap::Value Value;
  1.1455 +
  1.1456 +    /// The type of the matching map
  1.1457 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
  1.1458 +    MatchingMap;
  1.1459 +
  1.1460 +    /// \brief Scaling factor for primal solution
  1.1461 +    ///
  1.1462 +    /// Scaling factor for primal solution. It is equal to 2 or 1
  1.1463 +    /// according to the value type.
  1.1464 +    static const int primalScale =
  1.1465 +      std::numeric_limits<Value>::is_integer ? 2 : 1;
  1.1466 +
  1.1467 +    /// \brief Scaling factor for dual solution
  1.1468 +    ///
  1.1469 +    /// Scaling factor for dual solution. It is equal to 4 or 1
  1.1470 +    /// according to the value type.
  1.1471 +    static const int dualScale =
  1.1472 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
  1.1473 +
  1.1474 +  private:
  1.1475 +
  1.1476 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
  1.1477 +
  1.1478 +    typedef typename Graph::template NodeMap<Value> NodePotential;
  1.1479 +
  1.1480 +    const Graph& _graph;
  1.1481 +    const WeightMap& _weight;
  1.1482 +
  1.1483 +    MatchingMap* _matching;
  1.1484 +    NodePotential* _node_potential;
  1.1485 +
  1.1486 +    int _node_num;
  1.1487 +    bool _allow_loops;
  1.1488 +
  1.1489 +    enum Status {
  1.1490 +      EVEN = -1, MATCHED = 0, ODD = 1
  1.1491 +    };
  1.1492 +
  1.1493 +    typedef typename Graph::template NodeMap<Status> StatusMap;
  1.1494 +    StatusMap* _status;
  1.1495 +
  1.1496 +    typedef typename Graph::template NodeMap<Arc> PredMap;
  1.1497 +    PredMap* _pred;
  1.1498 +
  1.1499 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
  1.1500 +
  1.1501 +    IntNodeMap *_tree_set_index;
  1.1502 +    TreeSet *_tree_set;
  1.1503 +
  1.1504 +    IntNodeMap *_delta2_index;
  1.1505 +    BinHeap<Value, IntNodeMap> *_delta2;
  1.1506 +
  1.1507 +    IntEdgeMap *_delta3_index;
  1.1508 +    BinHeap<Value, IntEdgeMap> *_delta3;
  1.1509 +
  1.1510 +    Value _delta_sum;
  1.1511 +
  1.1512 +    void createStructures() {
  1.1513 +      _node_num = countNodes(_graph);
  1.1514 +
  1.1515 +      if (!_matching) {
  1.1516 +        _matching = new MatchingMap(_graph);
  1.1517 +      }
  1.1518 +      if (!_node_potential) {
  1.1519 +        _node_potential = new NodePotential(_graph);
  1.1520 +      }
  1.1521 +      if (!_status) {
  1.1522 +        _status = new StatusMap(_graph);
  1.1523 +      }
  1.1524 +      if (!_pred) {
  1.1525 +        _pred = new PredMap(_graph);
  1.1526 +      }
  1.1527 +      if (!_tree_set) {
  1.1528 +        _tree_set_index = new IntNodeMap(_graph);
  1.1529 +        _tree_set = new TreeSet(*_tree_set_index);
  1.1530 +      }
  1.1531 +      if (!_delta2) {
  1.1532 +        _delta2_index = new IntNodeMap(_graph);
  1.1533 +        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
  1.1534 +      }
  1.1535 +      if (!_delta3) {
  1.1536 +        _delta3_index = new IntEdgeMap(_graph);
  1.1537 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  1.1538 +      }
  1.1539 +    }
  1.1540 +
  1.1541 +    void destroyStructures() {
  1.1542 +      if (_matching) {
  1.1543 +        delete _matching;
  1.1544 +      }
  1.1545 +      if (_node_potential) {
  1.1546 +        delete _node_potential;
  1.1547 +      }
  1.1548 +      if (_status) {
  1.1549 +        delete _status;
  1.1550 +      }
  1.1551 +      if (_pred) {
  1.1552 +        delete _pred;
  1.1553 +      }
  1.1554 +      if (_tree_set) {
  1.1555 +        delete _tree_set_index;
  1.1556 +        delete _tree_set;
  1.1557 +      }
  1.1558 +      if (_delta2) {
  1.1559 +        delete _delta2_index;
  1.1560 +        delete _delta2;
  1.1561 +      }
  1.1562 +      if (_delta3) {
  1.1563 +        delete _delta3_index;
  1.1564 +        delete _delta3;
  1.1565 +      }
  1.1566 +    }
  1.1567 +
  1.1568 +    void matchedToEven(Node node, int tree) {
  1.1569 +      _tree_set->insert(node, tree);
  1.1570 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  1.1571 +
  1.1572 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
  1.1573 +        _delta2->erase(node);
  1.1574 +      }
  1.1575 +
  1.1576 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1.1577 +        Node v = _graph.source(a);
  1.1578 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1.1579 +          dualScale * _weight[a];
  1.1580 +        if (node == v) {
  1.1581 +          if (_allow_loops && _graph.direction(a)) {
  1.1582 +            _delta3->push(a, rw / 2);
  1.1583 +          }
  1.1584 +        } else if ((*_status)[v] == EVEN) {
  1.1585 +          _delta3->push(a, rw / 2);
  1.1586 +        } else if ((*_status)[v] == MATCHED) {
  1.1587 +          if (_delta2->state(v) != _delta2->IN_HEAP) {
  1.1588 +            _pred->set(v, a);
  1.1589 +            _delta2->push(v, rw);
  1.1590 +          } else if ((*_delta2)[v] > rw) {
  1.1591 +            _pred->set(v, a);
  1.1592 +            _delta2->decrease(v, rw);
  1.1593 +          }
  1.1594 +        }
  1.1595 +      }
  1.1596 +    }
  1.1597 +
  1.1598 +    void matchedToOdd(Node node, int tree) {
  1.1599 +      _tree_set->insert(node, tree);
  1.1600 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  1.1601 +
  1.1602 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
  1.1603 +        _delta2->erase(node);
  1.1604 +      }
  1.1605 +    }
  1.1606 +
  1.1607 +    void evenToMatched(Node node, int tree) {
  1.1608 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  1.1609 +      Arc min = INVALID;
  1.1610 +      Value minrw = std::numeric_limits<Value>::max();
  1.1611 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1.1612 +        Node v = _graph.source(a);
  1.1613 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1.1614 +          dualScale * _weight[a];
  1.1615 +
  1.1616 +        if (node == v) {
  1.1617 +          if (_allow_loops && _graph.direction(a)) {
  1.1618 +            _delta3->erase(a);
  1.1619 +          }
  1.1620 +        } else if ((*_status)[v] == EVEN) {
  1.1621 +          _delta3->erase(a);
  1.1622 +          if (minrw > rw) {
  1.1623 +            min = _graph.oppositeArc(a);
  1.1624 +            minrw = rw;
  1.1625 +          }
  1.1626 +        } else if ((*_status)[v]  == MATCHED) {
  1.1627 +          if ((*_pred)[v] == a) {
  1.1628 +            Arc mina = INVALID;
  1.1629 +            Value minrwa = std::numeric_limits<Value>::max();
  1.1630 +            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
  1.1631 +              Node va = _graph.target(aa);
  1.1632 +              if ((*_status)[va] != EVEN ||
  1.1633 +                  _tree_set->find(va) == tree) continue;
  1.1634 +              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
  1.1635 +                dualScale * _weight[aa];
  1.1636 +              if (minrwa > rwa) {
  1.1637 +                minrwa = rwa;
  1.1638 +                mina = aa;
  1.1639 +              }
  1.1640 +            }
  1.1641 +            if (mina != INVALID) {
  1.1642 +              _pred->set(v, mina);
  1.1643 +              _delta2->increase(v, minrwa);
  1.1644 +            } else {
  1.1645 +              _pred->set(v, INVALID);
  1.1646 +              _delta2->erase(v);
  1.1647 +            }
  1.1648 +          }
  1.1649 +        }
  1.1650 +      }
  1.1651 +      if (min != INVALID) {
  1.1652 +        _pred->set(node, min);
  1.1653 +        _delta2->push(node, minrw);
  1.1654 +      } else {
  1.1655 +        _pred->set(node, INVALID);
  1.1656 +      }
  1.1657 +    }
  1.1658 +
  1.1659 +    void oddToMatched(Node node) {
  1.1660 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  1.1661 +      Arc min = INVALID;
  1.1662 +      Value minrw = std::numeric_limits<Value>::max();
  1.1663 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1.1664 +        Node v = _graph.source(a);
  1.1665 +        if ((*_status)[v] != EVEN) continue;
  1.1666 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1.1667 +          dualScale * _weight[a];
  1.1668 +
  1.1669 +        if (minrw > rw) {
  1.1670 +          min = _graph.oppositeArc(a);
  1.1671 +          minrw = rw;
  1.1672 +        }
  1.1673 +      }
  1.1674 +      if (min != INVALID) {
  1.1675 +        _pred->set(node, min);
  1.1676 +        _delta2->push(node, minrw);
  1.1677 +      } else {
  1.1678 +        _pred->set(node, INVALID);
  1.1679 +      }
  1.1680 +    }
  1.1681 +
  1.1682 +    void alternatePath(Node even, int tree) {
  1.1683 +      Node odd;
  1.1684 +
  1.1685 +      _status->set(even, MATCHED);
  1.1686 +      evenToMatched(even, tree);
  1.1687 +
  1.1688 +      Arc prev = (*_matching)[even];
  1.1689 +      while (prev != INVALID) {
  1.1690 +        odd = _graph.target(prev);
  1.1691 +        even = _graph.target((*_pred)[odd]);
  1.1692 +        _matching->set(odd, (*_pred)[odd]);
  1.1693 +        _status->set(odd, MATCHED);
  1.1694 +        oddToMatched(odd);
  1.1695 +
  1.1696 +        prev = (*_matching)[even];
  1.1697 +        _status->set(even, MATCHED);
  1.1698 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
  1.1699 +        evenToMatched(even, tree);
  1.1700 +      }
  1.1701 +    }
  1.1702 +
  1.1703 +    void destroyTree(int tree) {
  1.1704 +      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
  1.1705 +        if ((*_status)[n] == EVEN) {
  1.1706 +          _status->set(n, MATCHED);
  1.1707 +          evenToMatched(n, tree);
  1.1708 +        } else if ((*_status)[n] == ODD) {
  1.1709 +          _status->set(n, MATCHED);
  1.1710 +          oddToMatched(n);
  1.1711 +        }
  1.1712 +      }
  1.1713 +      _tree_set->eraseClass(tree);
  1.1714 +    }
  1.1715 +
  1.1716 +    void augmentOnEdge(const Edge& edge) {
  1.1717 +      Node left = _graph.u(edge);
  1.1718 +      int left_tree = _tree_set->find(left);
  1.1719 +
  1.1720 +      alternatePath(left, left_tree);
  1.1721 +      destroyTree(left_tree);
  1.1722 +      _matching->set(left, _graph.direct(edge, true));
  1.1723 +
  1.1724 +      Node right = _graph.v(edge);
  1.1725 +      int right_tree = _tree_set->find(right);
  1.1726 +
  1.1727 +      alternatePath(right, right_tree);
  1.1728 +      destroyTree(right_tree);
  1.1729 +      _matching->set(right, _graph.direct(edge, false));
  1.1730 +    }
  1.1731 +
  1.1732 +    void augmentOnArc(const Arc& arc) {
  1.1733 +      Node left = _graph.source(arc);
  1.1734 +      _status->set(left, MATCHED);
  1.1735 +      _matching->set(left, arc);
  1.1736 +      _pred->set(left, arc);
  1.1737 +
  1.1738 +      Node right = _graph.target(arc);
  1.1739 +      int right_tree = _tree_set->find(right);
  1.1740 +
  1.1741 +      alternatePath(right, right_tree);
  1.1742 +      destroyTree(right_tree);
  1.1743 +      _matching->set(right, _graph.oppositeArc(arc));
  1.1744 +    }
  1.1745 +
  1.1746 +    void extendOnArc(const Arc& arc) {
  1.1747 +      Node base = _graph.target(arc);
  1.1748 +      int tree = _tree_set->find(base);
  1.1749 +
  1.1750 +      Node odd = _graph.source(arc);
  1.1751 +      _tree_set->insert(odd, tree);
  1.1752 +      _status->set(odd, ODD);
  1.1753 +      matchedToOdd(odd, tree);
  1.1754 +      _pred->set(odd, arc);
  1.1755 +
  1.1756 +      Node even = _graph.target((*_matching)[odd]);
  1.1757 +      _tree_set->insert(even, tree);
  1.1758 +      _status->set(even, EVEN);
  1.1759 +      matchedToEven(even, tree);
  1.1760 +    }
  1.1761 +
  1.1762 +    void cycleOnEdge(const Edge& edge, int tree) {
  1.1763 +      Node nca = INVALID;
  1.1764 +      std::vector<Node> left_path, right_path;
  1.1765 +
  1.1766 +      {
  1.1767 +        std::set<Node> left_set, right_set;
  1.1768 +        Node left = _graph.u(edge);
  1.1769 +        left_path.push_back(left);
  1.1770 +        left_set.insert(left);
  1.1771 +
  1.1772 +        Node right = _graph.v(edge);
  1.1773 +        right_path.push_back(right);
  1.1774 +        right_set.insert(right);
  1.1775 +
  1.1776 +        while (true) {
  1.1777 +
  1.1778 +          if (left_set.find(right) != left_set.end()) {
  1.1779 +            nca = right;
  1.1780 +            break;
  1.1781 +          }
  1.1782 +
  1.1783 +          if ((*_matching)[left] == INVALID) break;
  1.1784 +
  1.1785 +          left = _graph.target((*_matching)[left]);
  1.1786 +          left_path.push_back(left);
  1.1787 +          left = _graph.target((*_pred)[left]);
  1.1788 +          left_path.push_back(left);
  1.1789 +
  1.1790 +          left_set.insert(left);
  1.1791 +
  1.1792 +          if (right_set.find(left) != right_set.end()) {
  1.1793 +            nca = left;
  1.1794 +            break;
  1.1795 +          }
  1.1796 +
  1.1797 +          if ((*_matching)[right] == INVALID) break;
  1.1798 +
  1.1799 +          right = _graph.target((*_matching)[right]);
  1.1800 +          right_path.push_back(right);
  1.1801 +          right = _graph.target((*_pred)[right]);
  1.1802 +          right_path.push_back(right);
  1.1803 +
  1.1804 +          right_set.insert(right);
  1.1805 +
  1.1806 +        }
  1.1807 +
  1.1808 +        if (nca == INVALID) {
  1.1809 +          if ((*_matching)[left] == INVALID) {
  1.1810 +            nca = right;
  1.1811 +            while (left_set.find(nca) == left_set.end()) {
  1.1812 +              nca = _graph.target((*_matching)[nca]);
  1.1813 +              right_path.push_back(nca);
  1.1814 +              nca = _graph.target((*_pred)[nca]);
  1.1815 +              right_path.push_back(nca);
  1.1816 +            }
  1.1817 +          } else {
  1.1818 +            nca = left;
  1.1819 +            while (right_set.find(nca) == right_set.end()) {
  1.1820 +              nca = _graph.target((*_matching)[nca]);
  1.1821 +              left_path.push_back(nca);
  1.1822 +              nca = _graph.target((*_pred)[nca]);
  1.1823 +              left_path.push_back(nca);
  1.1824 +            }
  1.1825 +          }
  1.1826 +        }
  1.1827 +      }
  1.1828 +
  1.1829 +      alternatePath(nca, tree);
  1.1830 +      Arc prev;
  1.1831 +
  1.1832 +      prev = _graph.direct(edge, true);
  1.1833 +      for (int i = 0; left_path[i] != nca; i += 2) {
  1.1834 +        _matching->set(left_path[i], prev);
  1.1835 +        _status->set(left_path[i], MATCHED);
  1.1836 +        evenToMatched(left_path[i], tree);
  1.1837 +
  1.1838 +        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  1.1839 +        _status->set(left_path[i + 1], MATCHED);
  1.1840 +        oddToMatched(left_path[i + 1]);
  1.1841 +      }
  1.1842 +      _matching->set(nca, prev);
  1.1843 +
  1.1844 +      for (int i = 0; right_path[i] != nca; i += 2) {
  1.1845 +        _status->set(right_path[i], MATCHED);
  1.1846 +        evenToMatched(right_path[i], tree);
  1.1847 +
  1.1848 +        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  1.1849 +        _status->set(right_path[i + 1], MATCHED);
  1.1850 +        oddToMatched(right_path[i + 1]);
  1.1851 +      }
  1.1852 +
  1.1853 +      destroyTree(tree);
  1.1854 +    }
  1.1855 +
  1.1856 +    void extractCycle(const Arc &arc) {
  1.1857 +      Node left = _graph.source(arc);
  1.1858 +      Node odd = _graph.target((*_matching)[left]);
  1.1859 +      Arc prev;
  1.1860 +      while (odd != left) {
  1.1861 +        Node even = _graph.target((*_matching)[odd]);
  1.1862 +        prev = (*_matching)[odd];
  1.1863 +        odd = _graph.target((*_matching)[even]);
  1.1864 +        _matching->set(even, _graph.oppositeArc(prev));
  1.1865 +      }
  1.1866 +      _matching->set(left, arc);
  1.1867 +
  1.1868 +      Node right = _graph.target(arc);
  1.1869 +      int right_tree = _tree_set->find(right);
  1.1870 +      alternatePath(right, right_tree);
  1.1871 +      destroyTree(right_tree);
  1.1872 +      _matching->set(right, _graph.oppositeArc(arc));
  1.1873 +    }
  1.1874 +
  1.1875 +  public:
  1.1876 +
  1.1877 +    /// \brief Constructor
  1.1878 +    ///
  1.1879 +    /// Constructor.
  1.1880 +    MaxWeightedPerfectFractionalMatching(const Graph& graph,
  1.1881 +                                         const WeightMap& weight,
  1.1882 +                                         bool allow_loops = true)
  1.1883 +      : _graph(graph), _weight(weight), _matching(0),
  1.1884 +      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  1.1885 +      _status(0),  _pred(0),
  1.1886 +      _tree_set_index(0), _tree_set(0),
  1.1887 +
  1.1888 +      _delta2_index(0), _delta2(0),
  1.1889 +      _delta3_index(0), _delta3(0),
  1.1890 +
  1.1891 +      _delta_sum() {}
  1.1892 +
  1.1893 +    ~MaxWeightedPerfectFractionalMatching() {
  1.1894 +      destroyStructures();
  1.1895 +    }
  1.1896 +
  1.1897 +    /// \name Execution Control
  1.1898 +    /// The simplest way to execute the algorithm is to use the
  1.1899 +    /// \ref run() member function.
  1.1900 +
  1.1901 +    ///@{
  1.1902 +
  1.1903 +    /// \brief Initialize the algorithm
  1.1904 +    ///
  1.1905 +    /// This function initializes the algorithm.
  1.1906 +    void init() {
  1.1907 +      createStructures();
  1.1908 +
  1.1909 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1910 +        (*_delta2_index)[n] = _delta2->PRE_HEAP;
  1.1911 +      }
  1.1912 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1913 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1.1914 +      }
  1.1915 +
  1.1916 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1917 +        Value max = - std::numeric_limits<Value>::max();
  1.1918 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1.1919 +          if (_graph.target(e) == n && !_allow_loops) continue;
  1.1920 +          if ((dualScale * _weight[e]) / 2 > max) {
  1.1921 +            max = (dualScale * _weight[e]) / 2;
  1.1922 +          }
  1.1923 +        }
  1.1924 +        _node_potential->set(n, max);
  1.1925 +
  1.1926 +        _tree_set->insert(n);
  1.1927 +
  1.1928 +        _matching->set(n, INVALID);
  1.1929 +        _status->set(n, EVEN);
  1.1930 +      }
  1.1931 +
  1.1932 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1933 +        Node left = _graph.u(e);
  1.1934 +        Node right = _graph.v(e);
  1.1935 +        if (left == right && !_allow_loops) continue;
  1.1936 +        _delta3->push(e, ((*_node_potential)[left] +
  1.1937 +                          (*_node_potential)[right] -
  1.1938 +                          dualScale * _weight[e]) / 2);
  1.1939 +      }
  1.1940 +    }
  1.1941 +
  1.1942 +    /// \brief Start the algorithm
  1.1943 +    ///
  1.1944 +    /// This function starts the algorithm.
  1.1945 +    ///
  1.1946 +    /// \pre \ref init() must be called before using this function.
  1.1947 +    bool start() {
  1.1948 +      enum OpType {
  1.1949 +        D2, D3
  1.1950 +      };
  1.1951 +
  1.1952 +      int unmatched = _node_num;
  1.1953 +      while (unmatched > 0) {
  1.1954 +        Value d2 = !_delta2->empty() ?
  1.1955 +          _delta2->prio() : std::numeric_limits<Value>::max();
  1.1956 +
  1.1957 +        Value d3 = !_delta3->empty() ?
  1.1958 +          _delta3->prio() : std::numeric_limits<Value>::max();
  1.1959 +
  1.1960 +        _delta_sum = d3; OpType ot = D3;
  1.1961 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1.1962 +
  1.1963 +        if (_delta_sum == std::numeric_limits<Value>::max()) {
  1.1964 +          return false;
  1.1965 +        }
  1.1966 +
  1.1967 +        switch (ot) {
  1.1968 +        case D2:
  1.1969 +          {
  1.1970 +            Node n = _delta2->top();
  1.1971 +            Arc a = (*_pred)[n];
  1.1972 +            if ((*_matching)[n] == INVALID) {
  1.1973 +              augmentOnArc(a);
  1.1974 +              --unmatched;
  1.1975 +            } else {
  1.1976 +              Node v = _graph.target((*_matching)[n]);
  1.1977 +              if ((*_matching)[n] !=
  1.1978 +                  _graph.oppositeArc((*_matching)[v])) {
  1.1979 +                extractCycle(a);
  1.1980 +                --unmatched;
  1.1981 +              } else {
  1.1982 +                extendOnArc(a);
  1.1983 +              }
  1.1984 +            }
  1.1985 +          } break;
  1.1986 +        case D3:
  1.1987 +          {
  1.1988 +            Edge e = _delta3->top();
  1.1989 +
  1.1990 +            Node left = _graph.u(e);
  1.1991 +            Node right = _graph.v(e);
  1.1992 +
  1.1993 +            int left_tree = _tree_set->find(left);
  1.1994 +            int right_tree = _tree_set->find(right);
  1.1995 +
  1.1996 +            if (left_tree == right_tree) {
  1.1997 +              cycleOnEdge(e, left_tree);
  1.1998 +              --unmatched;
  1.1999 +            } else {
  1.2000 +              augmentOnEdge(e);
  1.2001 +              unmatched -= 2;
  1.2002 +            }
  1.2003 +          } break;
  1.2004 +        }
  1.2005 +      }
  1.2006 +      return true;
  1.2007 +    }
  1.2008 +
  1.2009 +    /// \brief Run the algorithm.
  1.2010 +    ///
  1.2011 +    /// This method runs the \c %MaxWeightedMatching algorithm.
  1.2012 +    ///
  1.2013 +    /// \note mwfm.run() is just a shortcut of the following code.
  1.2014 +    /// \code
  1.2015 +    ///   mwpfm.init();
  1.2016 +    ///   mwpfm.start();
  1.2017 +    /// \endcode
  1.2018 +    bool run() {
  1.2019 +      init();
  1.2020 +      return start();
  1.2021 +    }
  1.2022 +
  1.2023 +    /// @}
  1.2024 +
  1.2025 +    /// \name Primal Solution
  1.2026 +    /// Functions to get the primal solution, i.e. the maximum weighted
  1.2027 +    /// matching.\n
  1.2028 +    /// Either \ref run() or \ref start() function should be called before
  1.2029 +    /// using them.
  1.2030 +
  1.2031 +    /// @{
  1.2032 +
  1.2033 +    /// \brief Return the weight of the matching.
  1.2034 +    ///
  1.2035 +    /// This function returns the weight of the found matching. This
  1.2036 +    /// value is scaled by \ref primalScale "primal scale".
  1.2037 +    ///
  1.2038 +    /// \pre Either run() or start() must be called before using this function.
  1.2039 +    Value matchingWeight() const {
  1.2040 +      Value sum = 0;
  1.2041 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.2042 +        if ((*_matching)[n] != INVALID) {
  1.2043 +          sum += _weight[(*_matching)[n]];
  1.2044 +        }
  1.2045 +      }
  1.2046 +      return sum * primalScale / 2;
  1.2047 +    }
  1.2048 +
  1.2049 +    /// \brief Return the number of covered nodes in the matching.
  1.2050 +    ///
  1.2051 +    /// This function returns the number of covered nodes in the matching.
  1.2052 +    ///
  1.2053 +    /// \pre Either run() or start() must be called before using this function.
  1.2054 +    int matchingSize() const {
  1.2055 +      int num = 0;
  1.2056 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.2057 +        if ((*_matching)[n] != INVALID) {
  1.2058 +          ++num;
  1.2059 +        }
  1.2060 +      }
  1.2061 +      return num;
  1.2062 +    }
  1.2063 +
  1.2064 +    /// \brief Return \c true if the given edge is in the matching.
  1.2065 +    ///
  1.2066 +    /// This function returns \c true if the given edge is in the
  1.2067 +    /// found matching. The result is scaled by \ref primalScale
  1.2068 +    /// "primal scale".
  1.2069 +    ///
  1.2070 +    /// \pre Either run() or start() must be called before using this function.
  1.2071 +    Value matching(const Edge& edge) const {
  1.2072 +      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  1.2073 +        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
  1.2074 +        * primalScale / 2;
  1.2075 +    }
  1.2076 +
  1.2077 +    /// \brief Return the fractional matching arc (or edge) incident
  1.2078 +    /// to the given node.
  1.2079 +    ///
  1.2080 +    /// This function returns one of the fractional matching arc (or
  1.2081 +    /// edge) incident to the given node in the found matching or \c
  1.2082 +    /// INVALID if the node is not covered by the matching or if the
  1.2083 +    /// node is on an odd length cycle then it is the successor edge
  1.2084 +    /// on the cycle.
  1.2085 +    ///
  1.2086 +    /// \pre Either run() or start() must be called before using this function.
  1.2087 +    Arc matching(const Node& node) const {
  1.2088 +      return (*_matching)[node];
  1.2089 +    }
  1.2090 +
  1.2091 +    /// \brief Return a const reference to the matching map.
  1.2092 +    ///
  1.2093 +    /// This function returns a const reference to a node map that stores
  1.2094 +    /// the matching arc (or edge) incident to each node.
  1.2095 +    const MatchingMap& matchingMap() const {
  1.2096 +      return *_matching;
  1.2097 +    }
  1.2098 +
  1.2099 +    /// @}
  1.2100 +
  1.2101 +    /// \name Dual Solution
  1.2102 +    /// Functions to get the dual solution.\n
  1.2103 +    /// Either \ref run() or \ref start() function should be called before
  1.2104 +    /// using them.
  1.2105 +
  1.2106 +    /// @{
  1.2107 +
  1.2108 +    /// \brief Return the value of the dual solution.
  1.2109 +    ///
  1.2110 +    /// This function returns the value of the dual solution.
  1.2111 +    /// It should be equal to the primal value scaled by \ref dualScale
  1.2112 +    /// "dual scale".
  1.2113 +    ///
  1.2114 +    /// \pre Either run() or start() must be called before using this function.
  1.2115 +    Value dualValue() const {
  1.2116 +      Value sum = 0;
  1.2117 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.2118 +        sum += nodeValue(n);
  1.2119 +      }
  1.2120 +      return sum;
  1.2121 +    }
  1.2122 +
  1.2123 +    /// \brief Return the dual value (potential) of the given node.
  1.2124 +    ///
  1.2125 +    /// This function returns the dual value (potential) of the given node.
  1.2126 +    ///
  1.2127 +    /// \pre Either run() or start() must be called before using this function.
  1.2128 +    Value nodeValue(const Node& n) const {
  1.2129 +      return (*_node_potential)[n];
  1.2130 +    }
  1.2131 +
  1.2132 +    /// @}
  1.2133 +
  1.2134 +  };
  1.2135 +
  1.2136 +} //END OF NAMESPACE LEMON
  1.2137 +
  1.2138 +#endif //LEMON_FRACTIONAL_MATCHING_H