author Balazs Dezso Fri, 25 Sep 2009 21:51:36 +0200 changeset 948 636dadefe1e6 parent 947 0513ccfea967 child 949 61120524af27
 doc/groups.dox file | annotate | diff | comparison | revisions lemon/Makefile.am file | annotate | diff | comparison | revisions lemon/fractional_matching.h file | annotate | diff | comparison | revisions test/CMakeLists.txt file | annotate | diff | comparison | revisions test/Makefile.am file | annotate | diff | comparison | revisions test/fractional_matching_test.cc file | annotate | diff | comparison | revisions
     1.1 --- a/doc/groups.dox	Sun Sep 20 21:38:24 2009 +0200
1.2 +++ b/doc/groups.dox	Fri Sep 25 21:51:36 2009 +0200
1.3 @@ -349,7 +349,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.12 @@ -470,6 +470,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 bipartite_matching.png
1.25  \image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth

     2.1 --- a/lemon/Makefile.am	Sun Sep 20 21:38:24 2009 +0200
2.2 +++ b/lemon/Makefile.am	Fri Sep 25 21:51:36 2009 +0200
2.3 @@ -81,6 +81,7 @@
2.4  	lemon/euler.h \
2.5  	lemon/fib_heap.h \
2.6  	lemon/fourary_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	Fri Sep 25 21:51:36 2009 +0200
3.3 @@ -0,0 +1,2135 @@
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.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.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 MaxWeightedMatching::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.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 is based on extensive use
3.639 +  /// of priority queues and provides \f$O(nm\log n)\f$ time
3.640 +  /// complexity.
3.641 +  ///
3.642 +  /// The maximum weighted fractional matching is a relaxation of the
3.643 +  /// maximum weighted matching problem where the odd set constraints
3.644 +  /// are omitted.
3.645 +  /// It can be formulated with the following linear program.
3.646 +  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
3.647 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
3.648 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
3.649 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.650 +  /// \f$X\f$. The result must be the union of a matching with one
3.651 +  /// value edges and a set of odd length cycles with half value edges.
3.652 +  ///
3.653 +  /// The algorithm calculates an optimal fractional matching and a
3.654 +  /// proof of the optimality. The solution of the dual problem can be
3.655 +  /// used to check the result of the algorithm. The dual linear
3.656 +  /// problem is the following.
3.657 +  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
3.658 +  /// \f[y_u \ge 0 \quad \forall u \in V\f]
3.659 +  /// \f[\min \sum_{u \in V}y_u \f] ///
3.660 +  ///
3.661 +  /// The algorithm can be executed with the run() function.
3.662 +  /// After it the matching (the primal solution) and the dual solution
3.663 +  /// can be obtained using the query functions.
3.664 +  ///
3.665 +  /// If the value type is integer, then the primal and the dual
3.666 +  /// solutions are multiplied by
3.667 +  /// \ref MaxWeightedMatching::primalScale "2" and
3.668 +  /// \ref MaxWeightedMatching::dualScale "4" respectively.
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. It is equal to 2 or 1
3.696 +    /// according to the value type.
3.697 +    static const int primalScale =
3.698 +      std::numeric_limits<Value>::is_integer ? 2 : 1;
3.699 +
3.700 +    /// \brief Scaling factor for dual solution
3.701 +    ///
3.702 +    /// Scaling factor for dual solution. It is equal to 4 or 1
3.703 +    /// according to the value type.
3.704 +    static const int dualScale =
3.705 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
3.706 +
3.707 +  private:
3.708 +
3.709 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
3.710 +
3.711 +    typedef typename Graph::template NodeMap<Value> NodePotential;
3.712 +
3.713 +    const Graph& _graph;
3.714 +    const WeightMap& _weight;
3.715 +
3.716 +    MatchingMap* _matching;
3.717 +    NodePotential* _node_potential;
3.718 +
3.719 +    int _node_num;
3.720 +    bool _allow_loops;
3.721 +
3.722 +    enum Status {
3.723 +      EVEN = -1, MATCHED = 0, ODD = 1
3.724 +    };
3.725 +
3.726 +    typedef typename Graph::template NodeMap<Status> StatusMap;
3.727 +    StatusMap* _status;
3.728 +
3.729 +    typedef typename Graph::template NodeMap<Arc> PredMap;
3.730 +    PredMap* _pred;
3.731 +
3.732 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
3.733 +
3.734 +    IntNodeMap *_tree_set_index;
3.735 +    TreeSet *_tree_set;
3.736 +
3.737 +    IntNodeMap *_delta1_index;
3.738 +    BinHeap<Value, IntNodeMap> *_delta1;
3.739 +
3.740 +    IntNodeMap *_delta2_index;
3.741 +    BinHeap<Value, IntNodeMap> *_delta2;
3.742 +
3.743 +    IntEdgeMap *_delta3_index;
3.744 +    BinHeap<Value, IntEdgeMap> *_delta3;
3.745 +
3.746 +    Value _delta_sum;
3.747 +
3.748 +    void createStructures() {
3.749 +      _node_num = countNodes(_graph);
3.750 +
3.751 +      if (!_matching) {
3.752 +        _matching = new MatchingMap(_graph);
3.753 +      }
3.754 +      if (!_node_potential) {
3.755 +        _node_potential = new NodePotential(_graph);
3.756 +      }
3.757 +      if (!_status) {
3.758 +        _status = new StatusMap(_graph);
3.759 +      }
3.760 +      if (!_pred) {
3.761 +        _pred = new PredMap(_graph);
3.762 +      }
3.763 +      if (!_tree_set) {
3.764 +        _tree_set_index = new IntNodeMap(_graph);
3.765 +        _tree_set = new TreeSet(*_tree_set_index);
3.766 +      }
3.767 +      if (!_delta1) {
3.768 +        _delta1_index = new IntNodeMap(_graph);
3.769 +        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
3.770 +      }
3.771 +      if (!_delta2) {
3.772 +        _delta2_index = new IntNodeMap(_graph);
3.773 +        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
3.774 +      }
3.775 +      if (!_delta3) {
3.776 +        _delta3_index = new IntEdgeMap(_graph);
3.777 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
3.778 +      }
3.779 +    }
3.780 +
3.781 +    void destroyStructures() {
3.782 +      if (_matching) {
3.783 +        delete _matching;
3.784 +      }
3.785 +      if (_node_potential) {
3.786 +        delete _node_potential;
3.787 +      }
3.788 +      if (_status) {
3.789 +        delete _status;
3.790 +      }
3.791 +      if (_pred) {
3.792 +        delete _pred;
3.793 +      }
3.794 +      if (_tree_set) {
3.795 +        delete _tree_set_index;
3.796 +        delete _tree_set;
3.797 +      }
3.798 +      if (_delta1) {
3.799 +        delete _delta1_index;
3.800 +        delete _delta1;
3.801 +      }
3.802 +      if (_delta2) {
3.803 +        delete _delta2_index;
3.804 +        delete _delta2;
3.805 +      }
3.806 +      if (_delta3) {
3.807 +        delete _delta3_index;
3.808 +        delete _delta3;
3.809 +      }
3.810 +    }
3.811 +
3.812 +    void matchedToEven(Node node, int tree) {
3.813 +      _tree_set->insert(node, tree);
3.814 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.815 +      _delta1->push(node, (*_node_potential)[node]);
3.816 +
3.817 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
3.818 +        _delta2->erase(node);
3.819 +      }
3.820 +
3.821 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.822 +        Node v = _graph.source(a);
3.823 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.824 +          dualScale * _weight[a];
3.825 +        if (node == v) {
3.826 +          if (_allow_loops && _graph.direction(a)) {
3.827 +            _delta3->push(a, rw / 2);
3.828 +          }
3.829 +        } else if ((*_status)[v] == EVEN) {
3.830 +          _delta3->push(a, rw / 2);
3.831 +        } else if ((*_status)[v] == MATCHED) {
3.832 +          if (_delta2->state(v) != _delta2->IN_HEAP) {
3.833 +            _pred->set(v, a);
3.834 +            _delta2->push(v, rw);
3.835 +          } else if ((*_delta2)[v] > rw) {
3.836 +            _pred->set(v, a);
3.837 +            _delta2->decrease(v, rw);
3.838 +          }
3.839 +        }
3.840 +      }
3.841 +    }
3.842 +
3.843 +    void matchedToOdd(Node node, int tree) {
3.844 +      _tree_set->insert(node, tree);
3.845 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.846 +
3.847 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
3.848 +        _delta2->erase(node);
3.849 +      }
3.850 +    }
3.851 +
3.852 +    void evenToMatched(Node node, int tree) {
3.853 +      _delta1->erase(node);
3.854 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.855 +      Arc min = INVALID;
3.856 +      Value minrw = std::numeric_limits<Value>::max();
3.857 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.858 +        Node v = _graph.source(a);
3.859 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.860 +          dualScale * _weight[a];
3.861 +
3.862 +        if (node == v) {
3.863 +          if (_allow_loops && _graph.direction(a)) {
3.864 +            _delta3->erase(a);
3.865 +          }
3.866 +        } else if ((*_status)[v] == EVEN) {
3.867 +          _delta3->erase(a);
3.868 +          if (minrw > rw) {
3.869 +            min = _graph.oppositeArc(a);
3.870 +            minrw = rw;
3.871 +          }
3.872 +        } else if ((*_status)[v]  == MATCHED) {
3.873 +          if ((*_pred)[v] == a) {
3.874 +            Arc mina = INVALID;
3.875 +            Value minrwa = std::numeric_limits<Value>::max();
3.876 +            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
3.877 +              Node va = _graph.target(aa);
3.878 +              if ((*_status)[va] != EVEN ||
3.879 +                  _tree_set->find(va) == tree) continue;
3.880 +              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
3.881 +                dualScale * _weight[aa];
3.882 +              if (minrwa > rwa) {
3.883 +                minrwa = rwa;
3.884 +                mina = aa;
3.885 +              }
3.886 +            }
3.887 +            if (mina != INVALID) {
3.888 +              _pred->set(v, mina);
3.889 +              _delta2->increase(v, minrwa);
3.890 +            } else {
3.891 +              _pred->set(v, INVALID);
3.892 +              _delta2->erase(v);
3.893 +            }
3.894 +          }
3.895 +        }
3.896 +      }
3.897 +      if (min != INVALID) {
3.898 +        _pred->set(node, min);
3.899 +        _delta2->push(node, minrw);
3.900 +      } else {
3.901 +        _pred->set(node, INVALID);
3.902 +      }
3.903 +    }
3.904 +
3.905 +    void oddToMatched(Node node) {
3.906 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.907 +      Arc min = INVALID;
3.908 +      Value minrw = std::numeric_limits<Value>::max();
3.909 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.910 +        Node v = _graph.source(a);
3.911 +        if ((*_status)[v] != EVEN) continue;
3.912 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.913 +          dualScale * _weight[a];
3.914 +
3.915 +        if (minrw > rw) {
3.916 +          min = _graph.oppositeArc(a);
3.917 +          minrw = rw;
3.918 +        }
3.919 +      }
3.920 +      if (min != INVALID) {
3.921 +        _pred->set(node, min);
3.922 +        _delta2->push(node, minrw);
3.923 +      } else {
3.924 +        _pred->set(node, INVALID);
3.925 +      }
3.926 +    }
3.927 +
3.928 +    void alternatePath(Node even, int tree) {
3.929 +      Node odd;
3.930 +
3.931 +      _status->set(even, MATCHED);
3.932 +      evenToMatched(even, tree);
3.933 +
3.934 +      Arc prev = (*_matching)[even];
3.935 +      while (prev != INVALID) {
3.936 +        odd = _graph.target(prev);
3.937 +        even = _graph.target((*_pred)[odd]);
3.938 +        _matching->set(odd, (*_pred)[odd]);
3.939 +        _status->set(odd, MATCHED);
3.940 +        oddToMatched(odd);
3.941 +
3.942 +        prev = (*_matching)[even];
3.943 +        _status->set(even, MATCHED);
3.944 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
3.945 +        evenToMatched(even, tree);
3.946 +      }
3.947 +    }
3.948 +
3.949 +    void destroyTree(int tree) {
3.950 +      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
3.951 +        if ((*_status)[n] == EVEN) {
3.952 +          _status->set(n, MATCHED);
3.953 +          evenToMatched(n, tree);
3.954 +        } else if ((*_status)[n] == ODD) {
3.955 +          _status->set(n, MATCHED);
3.956 +          oddToMatched(n);
3.957 +        }
3.958 +      }
3.959 +      _tree_set->eraseClass(tree);
3.960 +    }
3.961 +
3.962 +
3.963 +    void unmatchNode(const Node& node) {
3.964 +      int tree = _tree_set->find(node);
3.965 +
3.966 +      alternatePath(node, tree);
3.967 +      destroyTree(tree);
3.968 +
3.969 +      _matching->set(node, INVALID);
3.970 +    }
3.971 +
3.972 +
3.973 +    void augmentOnEdge(const Edge& edge) {
3.974 +      Node left = _graph.u(edge);
3.975 +      int left_tree = _tree_set->find(left);
3.976 +
3.977 +      alternatePath(left, left_tree);
3.978 +      destroyTree(left_tree);
3.979 +      _matching->set(left, _graph.direct(edge, true));
3.980 +
3.981 +      Node right = _graph.v(edge);
3.982 +      int right_tree = _tree_set->find(right);
3.983 +
3.984 +      alternatePath(right, right_tree);
3.985 +      destroyTree(right_tree);
3.986 +      _matching->set(right, _graph.direct(edge, false));
3.987 +    }
3.988 +
3.989 +    void augmentOnArc(const Arc& arc) {
3.990 +      Node left = _graph.source(arc);
3.991 +      _status->set(left, MATCHED);
3.992 +      _matching->set(left, arc);
3.993 +      _pred->set(left, arc);
3.994 +
3.995 +      Node right = _graph.target(arc);
3.996 +      int right_tree = _tree_set->find(right);
3.997 +
3.998 +      alternatePath(right, right_tree);
3.999 +      destroyTree(right_tree);
3.1000 +      _matching->set(right, _graph.oppositeArc(arc));
3.1001 +    }
3.1002 +
3.1003 +    void extendOnArc(const Arc& arc) {
3.1004 +      Node base = _graph.target(arc);
3.1005 +      int tree = _tree_set->find(base);
3.1006 +
3.1007 +      Node odd = _graph.source(arc);
3.1008 +      _tree_set->insert(odd, tree);
3.1009 +      _status->set(odd, ODD);
3.1010 +      matchedToOdd(odd, tree);
3.1011 +      _pred->set(odd, arc);
3.1012 +
3.1013 +      Node even = _graph.target((*_matching)[odd]);
3.1014 +      _tree_set->insert(even, tree);
3.1015 +      _status->set(even, EVEN);
3.1016 +      matchedToEven(even, tree);
3.1017 +    }
3.1018 +
3.1019 +    void cycleOnEdge(const Edge& edge, int tree) {
3.1020 +      Node nca = INVALID;
3.1021 +      std::vector<Node> left_path, right_path;
3.1022 +
3.1023 +      {
3.1024 +        std::set<Node> left_set, right_set;
3.1025 +        Node left = _graph.u(edge);
3.1026 +        left_path.push_back(left);
3.1027 +        left_set.insert(left);
3.1028 +
3.1029 +        Node right = _graph.v(edge);
3.1030 +        right_path.push_back(right);
3.1031 +        right_set.insert(right);
3.1032 +
3.1033 +        while (true) {
3.1034 +
3.1035 +          if (left_set.find(right) != left_set.end()) {
3.1036 +            nca = right;
3.1037 +            break;
3.1038 +          }
3.1039 +
3.1040 +          if ((*_matching)[left] == INVALID) break;
3.1041 +
3.1042 +          left = _graph.target((*_matching)[left]);
3.1043 +          left_path.push_back(left);
3.1044 +          left = _graph.target((*_pred)[left]);
3.1045 +          left_path.push_back(left);
3.1046 +
3.1047 +          left_set.insert(left);
3.1048 +
3.1049 +          if (right_set.find(left) != right_set.end()) {
3.1050 +            nca = left;
3.1051 +            break;
3.1052 +          }
3.1053 +
3.1054 +          if ((*_matching)[right] == INVALID) break;
3.1055 +
3.1056 +          right = _graph.target((*_matching)[right]);
3.1057 +          right_path.push_back(right);
3.1058 +          right = _graph.target((*_pred)[right]);
3.1059 +          right_path.push_back(right);
3.1060 +
3.1061 +          right_set.insert(right);
3.1062 +
3.1063 +        }
3.1064 +
3.1065 +        if (nca == INVALID) {
3.1066 +          if ((*_matching)[left] == INVALID) {
3.1067 +            nca = right;
3.1068 +            while (left_set.find(nca) == left_set.end()) {
3.1069 +              nca = _graph.target((*_matching)[nca]);
3.1070 +              right_path.push_back(nca);
3.1071 +              nca = _graph.target((*_pred)[nca]);
3.1072 +              right_path.push_back(nca);
3.1073 +            }
3.1074 +          } else {
3.1075 +            nca = left;
3.1076 +            while (right_set.find(nca) == right_set.end()) {
3.1077 +              nca = _graph.target((*_matching)[nca]);
3.1078 +              left_path.push_back(nca);
3.1079 +              nca = _graph.target((*_pred)[nca]);
3.1080 +              left_path.push_back(nca);
3.1081 +            }
3.1082 +          }
3.1083 +        }
3.1084 +      }
3.1085 +
3.1086 +      alternatePath(nca, tree);
3.1087 +      Arc prev;
3.1088 +
3.1089 +      prev = _graph.direct(edge, true);
3.1090 +      for (int i = 0; left_path[i] != nca; i += 2) {
3.1091 +        _matching->set(left_path[i], prev);
3.1092 +        _status->set(left_path[i], MATCHED);
3.1093 +        evenToMatched(left_path[i], tree);
3.1094 +
3.1095 +        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
3.1096 +        _status->set(left_path[i + 1], MATCHED);
3.1097 +        oddToMatched(left_path[i + 1]);
3.1098 +      }
3.1099 +      _matching->set(nca, prev);
3.1100 +
3.1101 +      for (int i = 0; right_path[i] != nca; i += 2) {
3.1102 +        _status->set(right_path[i], MATCHED);
3.1103 +        evenToMatched(right_path[i], tree);
3.1104 +
3.1105 +        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
3.1106 +        _status->set(right_path[i + 1], MATCHED);
3.1107 +        oddToMatched(right_path[i + 1]);
3.1108 +      }
3.1109 +
3.1110 +      destroyTree(tree);
3.1111 +    }
3.1112 +
3.1113 +    void extractCycle(const Arc &arc) {
3.1114 +      Node left = _graph.source(arc);
3.1115 +      Node odd = _graph.target((*_matching)[left]);
3.1116 +      Arc prev;
3.1117 +      while (odd != left) {
3.1118 +        Node even = _graph.target((*_matching)[odd]);
3.1119 +        prev = (*_matching)[odd];
3.1120 +        odd = _graph.target((*_matching)[even]);
3.1121 +        _matching->set(even, _graph.oppositeArc(prev));
3.1122 +      }
3.1123 +      _matching->set(left, arc);
3.1124 +
3.1125 +      Node right = _graph.target(arc);
3.1126 +      int right_tree = _tree_set->find(right);
3.1127 +      alternatePath(right, right_tree);
3.1128 +      destroyTree(right_tree);
3.1129 +      _matching->set(right, _graph.oppositeArc(arc));
3.1130 +    }
3.1131 +
3.1132 +  public:
3.1133 +
3.1134 +    /// \brief Constructor
3.1135 +    ///
3.1136 +    /// Constructor.
3.1137 +    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
3.1138 +                                  bool allow_loops = true)
3.1139 +      : _graph(graph), _weight(weight), _matching(0),
3.1140 +      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
3.1141 +      _status(0),  _pred(0),
3.1142 +      _tree_set_index(0), _tree_set(0),
3.1143 +
3.1144 +      _delta1_index(0), _delta1(0),
3.1145 +      _delta2_index(0), _delta2(0),
3.1146 +      _delta3_index(0), _delta3(0),
3.1147 +
3.1148 +      _delta_sum() {}
3.1149 +
3.1150 +    ~MaxWeightedFractionalMatching() {
3.1151 +      destroyStructures();
3.1152 +    }
3.1153 +
3.1154 +    /// \name Execution Control
3.1155 +    /// The simplest way to execute the algorithm is to use the
3.1156 +    /// \ref run() member function.
3.1157 +
3.1158 +    ///@{
3.1159 +
3.1160 +    /// \brief Initialize the algorithm
3.1161 +    ///
3.1162 +    /// This function initializes the algorithm.
3.1163 +    void init() {
3.1164 +      createStructures();
3.1165 +
3.1166 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1167 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
3.1168 +        (*_delta2_index)[n] = _delta2->PRE_HEAP;
3.1169 +      }
3.1170 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1171 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
3.1172 +      }
3.1173 +
3.1174 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1175 +        Value max = 0;
3.1176 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3.1177 +          if (_graph.target(e) == n && !_allow_loops) continue;
3.1178 +          if ((dualScale * _weight[e]) / 2 > max) {
3.1179 +            max = (dualScale * _weight[e]) / 2;
3.1180 +          }
3.1181 +        }
3.1182 +        _node_potential->set(n, max);
3.1183 +        _delta1->push(n, max);
3.1184 +
3.1185 +        _tree_set->insert(n);
3.1186 +
3.1187 +        _matching->set(n, INVALID);
3.1188 +        _status->set(n, EVEN);
3.1189 +      }
3.1190 +
3.1191 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1192 +        Node left = _graph.u(e);
3.1193 +        Node right = _graph.v(e);
3.1194 +        if (left == right && !_allow_loops) continue;
3.1195 +        _delta3->push(e, ((*_node_potential)[left] +
3.1196 +                          (*_node_potential)[right] -
3.1197 +                          dualScale * _weight[e]) / 2);
3.1198 +      }
3.1199 +    }
3.1200 +
3.1201 +    /// \brief Start the algorithm
3.1202 +    ///
3.1203 +    /// This function starts the algorithm.
3.1204 +    ///
3.1205 +    /// \pre \ref init() must be called before using this function.
3.1206 +    void start() {
3.1207 +      enum OpType {
3.1208 +        D1, D2, D3
3.1209 +      };
3.1210 +
3.1211 +      int unmatched = _node_num;
3.1212 +      while (unmatched > 0) {
3.1213 +        Value d1 = !_delta1->empty() ?
3.1214 +          _delta1->prio() : std::numeric_limits<Value>::max();
3.1215 +
3.1216 +        Value d2 = !_delta2->empty() ?
3.1217 +          _delta2->prio() : std::numeric_limits<Value>::max();
3.1218 +
3.1219 +        Value d3 = !_delta3->empty() ?
3.1220 +          _delta3->prio() : std::numeric_limits<Value>::max();
3.1221 +
3.1222 +        _delta_sum = d3; OpType ot = D3;
3.1223 +        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
3.1224 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3.1225 +
3.1226 +        switch (ot) {
3.1227 +        case D1:
3.1228 +          {
3.1229 +            Node n = _delta1->top();
3.1230 +            unmatchNode(n);
3.1231 +            --unmatched;
3.1232 +          }
3.1233 +          break;
3.1234 +        case D2:
3.1235 +          {
3.1236 +            Node n = _delta2->top();
3.1237 +            Arc a = (*_pred)[n];
3.1238 +            if ((*_matching)[n] == INVALID) {
3.1239 +              augmentOnArc(a);
3.1240 +              --unmatched;
3.1241 +            } else {
3.1242 +              Node v = _graph.target((*_matching)[n]);
3.1243 +              if ((*_matching)[n] !=
3.1244 +                  _graph.oppositeArc((*_matching)[v])) {
3.1245 +                extractCycle(a);
3.1246 +                --unmatched;
3.1247 +              } else {
3.1248 +                extendOnArc(a);
3.1249 +              }
3.1250 +            }
3.1251 +          } break;
3.1252 +        case D3:
3.1253 +          {
3.1254 +            Edge e = _delta3->top();
3.1255 +
3.1256 +            Node left = _graph.u(e);
3.1257 +            Node right = _graph.v(e);
3.1258 +
3.1259 +            int left_tree = _tree_set->find(left);
3.1260 +            int right_tree = _tree_set->find(right);
3.1261 +
3.1262 +            if (left_tree == right_tree) {
3.1263 +              cycleOnEdge(e, left_tree);
3.1264 +              --unmatched;
3.1265 +            } else {
3.1266 +              augmentOnEdge(e);
3.1267 +              unmatched -= 2;
3.1268 +            }
3.1269 +          } break;
3.1270 +        }
3.1271 +      }
3.1272 +    }
3.1273 +
3.1274 +    /// \brief Run the algorithm.
3.1275 +    ///
3.1276 +    /// This method runs the \c %MaxWeightedMatching algorithm.
3.1277 +    ///
3.1278 +    /// \note mwfm.run() is just a shortcut of the following code.
3.1279 +    /// \code
3.1280 +    ///   mwfm.init();
3.1281 +    ///   mwfm.start();
3.1282 +    /// \endcode
3.1283 +    void run() {
3.1284 +      init();
3.1285 +      start();
3.1286 +    }
3.1287 +
3.1288 +    /// @}
3.1289 +
3.1290 +    /// \name Primal Solution
3.1291 +    /// Functions to get the primal solution, i.e. the maximum weighted
3.1292 +    /// matching.\n
3.1293 +    /// Either \ref run() or \ref start() function should be called before
3.1294 +    /// using them.
3.1295 +
3.1296 +    /// @{
3.1297 +
3.1298 +    /// \brief Return the weight of the matching.
3.1299 +    ///
3.1300 +    /// This function returns the weight of the found matching. This
3.1301 +    /// value is scaled by \ref primalScale "primal scale".
3.1302 +    ///
3.1303 +    /// \pre Either run() or start() must be called before using this function.
3.1304 +    Value matchingWeight() const {
3.1305 +      Value sum = 0;
3.1306 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1307 +        if ((*_matching)[n] != INVALID) {
3.1308 +          sum += _weight[(*_matching)[n]];
3.1309 +        }
3.1310 +      }
3.1311 +      return sum * primalScale / 2;
3.1312 +    }
3.1313 +
3.1314 +    /// \brief Return the number of covered nodes in the matching.
3.1315 +    ///
3.1316 +    /// This function returns the number of covered nodes in the matching.
3.1317 +    ///
3.1318 +    /// \pre Either run() or start() must be called before using this function.
3.1319 +    int matchingSize() const {
3.1320 +      int num = 0;
3.1321 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1322 +        if ((*_matching)[n] != INVALID) {
3.1323 +          ++num;
3.1324 +        }
3.1325 +      }
3.1326 +      return num;
3.1327 +    }
3.1328 +
3.1329 +    /// \brief Return \c true if the given edge is in the matching.
3.1330 +    ///
3.1331 +    /// This function returns \c true if the given edge is in the
3.1332 +    /// found matching. The result is scaled by \ref primalScale
3.1333 +    /// "primal scale".
3.1334 +    ///
3.1335 +    /// \pre Either run() or start() must be called before using this function.
3.1336 +    Value matching(const Edge& edge) const {
3.1337 +      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
3.1338 +        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
3.1339 +        * primalScale / 2;
3.1340 +    }
3.1341 +
3.1342 +    /// \brief Return the fractional matching arc (or edge) incident
3.1343 +    /// to the given node.
3.1344 +    ///
3.1345 +    /// This function returns one of the fractional matching arc (or
3.1346 +    /// edge) incident to the given node in the found matching or \c
3.1347 +    /// INVALID if the node is not covered by the matching or if the
3.1348 +    /// node is on an odd length cycle then it is the successor edge
3.1349 +    /// on the cycle.
3.1350 +    ///
3.1351 +    /// \pre Either run() or start() must be called before using this function.
3.1352 +    Arc matching(const Node& node) const {
3.1353 +      return (*_matching)[node];
3.1354 +    }
3.1355 +
3.1356 +    /// \brief Return a const reference to the matching map.
3.1357 +    ///
3.1358 +    /// This function returns a const reference to a node map that stores
3.1359 +    /// the matching arc (or edge) incident to each node.
3.1360 +    const MatchingMap& matchingMap() const {
3.1361 +      return *_matching;
3.1362 +    }
3.1363 +
3.1364 +    /// @}
3.1365 +
3.1366 +    /// \name Dual Solution
3.1367 +    /// Functions to get the dual solution.\n
3.1368 +    /// Either \ref run() or \ref start() function should be called before
3.1369 +    /// using them.
3.1370 +
3.1371 +    /// @{
3.1372 +
3.1373 +    /// \brief Return the value of the dual solution.
3.1374 +    ///
3.1375 +    /// This function returns the value of the dual solution.
3.1376 +    /// It should be equal to the primal value scaled by \ref dualScale
3.1377 +    /// "dual scale".
3.1378 +    ///
3.1379 +    /// \pre Either run() or start() must be called before using this function.
3.1380 +    Value dualValue() const {
3.1381 +      Value sum = 0;
3.1382 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1383 +        sum += nodeValue(n);
3.1384 +      }
3.1385 +      return sum;
3.1386 +    }
3.1387 +
3.1388 +    /// \brief Return the dual value (potential) of the given node.
3.1389 +    ///
3.1390 +    /// This function returns the dual value (potential) of the given node.
3.1391 +    ///
3.1392 +    /// \pre Either run() or start() must be called before using this function.
3.1393 +    Value nodeValue(const Node& n) const {
3.1394 +      return (*_node_potential)[n];
3.1395 +    }
3.1396 +
3.1397 +    /// @}
3.1398 +
3.1399 +  };
3.1400 +
3.1401 +  /// \ingroup matching
3.1402 +  ///
3.1403 +  /// \brief Weighted fractional perfect matching in general graphs
3.1404 +  ///
3.1405 +  /// This class provides an efficient implementation of fractional
3.1406 +  /// matching algorithm. The implementation is based on extensive use
3.1407 +  /// of priority queues and provides \f$O(nm\log n)\f$ time
3.1408 +  /// complexity.
3.1409 +  ///
3.1410 +  /// The maximum weighted fractional perfect matching is a relaxation
3.1411 +  /// of the maximum weighted perfect matching problem where the odd
3.1412 +  /// set constraints are omitted.
3.1413 +  /// It can be formulated with the following linear program.
3.1414 +  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
3.1415 +  /// \f[x_e \ge 0\quad \forall e\in E\f]
3.1416 +  /// \f[\max \sum_{e\in E}x_ew_e\f]
3.1417 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.1418 +  /// \f$X\f$. The result must be the union of a matching with one
3.1419 +  /// value edges and a set of odd length cycles with half value edges.
3.1420 +  ///
3.1421 +  /// The algorithm calculates an optimal fractional matching and a
3.1422 +  /// proof of the optimality. The solution of the dual problem can be
3.1423 +  /// used to check the result of the algorithm. The dual linear
3.1424 +  /// problem is the following.
3.1425 +  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
3.1426 +  /// \f[\min \sum_{u \in V}y_u \f] ///
3.1427 +  ///
3.1428 +  /// The algorithm can be executed with the run() function.
3.1429 +  /// After it the matching (the primal solution) and the dual solution
3.1430 +  /// can be obtained using the query functions.
3.1431 +
3.1432 +  /// If the value type is integer, then the primal and the dual
3.1433 +  /// solutions are multiplied by
3.1434 +  /// \ref MaxWeightedMatching::primalScale "2" and
3.1435 +  /// \ref MaxWeightedMatching::dualScale "4" respectively.
3.1436 +  ///
3.1437 +  /// \tparam GR The undirected graph type the algorithm runs on.
3.1438 +  /// \tparam WM The type edge weight map. The default type is
3.1439 +  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
3.1440 +#ifdef DOXYGEN
3.1441 +  template <typename GR, typename WM>
3.1442 +#else
3.1443 +  template <typename GR,
3.1444 +            typename WM = typename GR::template EdgeMap<int> >
3.1445 +#endif
3.1446 +  class MaxWeightedPerfectFractionalMatching {
3.1447 +  public:
3.1448 +
3.1449 +    /// The graph type of the algorithm
3.1450 +    typedef GR Graph;
3.1451 +    /// The type of the edge weight map
3.1452 +    typedef WM WeightMap;
3.1453 +    /// The value type of the edge weights
3.1454 +    typedef typename WeightMap::Value Value;
3.1455 +
3.1456 +    /// The type of the matching map
3.1457 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
3.1458 +    MatchingMap;
3.1459 +
3.1460 +    /// \brief Scaling factor for primal solution
3.1461 +    ///
3.1462 +    /// Scaling factor for primal solution. It is equal to 2 or 1
3.1463 +    /// according to the value type.
3.1464 +    static const int primalScale =
3.1465 +      std::numeric_limits<Value>::is_integer ? 2 : 1;
3.1466 +
3.1467 +    /// \brief Scaling factor for dual solution
3.1468 +    ///
3.1469 +    /// Scaling factor for dual solution. It is equal to 4 or 1
3.1470 +    /// according to the value type.
3.1471 +    static const int dualScale =
3.1472 +      std::numeric_limits<Value>::is_integer ? 4 : 1;
3.1473 +
3.1474 +  private:
3.1475 +
3.1476 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
3.1477 +
3.1478 +    typedef typename Graph::template NodeMap<Value> NodePotential;
3.1479 +
3.1480 +    const Graph& _graph;
3.1481 +    const WeightMap& _weight;
3.1482 +
3.1483 +    MatchingMap* _matching;
3.1484 +    NodePotential* _node_potential;
3.1485 +
3.1486 +    int _node_num;
3.1487 +    bool _allow_loops;
3.1488 +
3.1489 +    enum Status {
3.1490 +      EVEN = -1, MATCHED = 0, ODD = 1
3.1491 +    };
3.1492 +
3.1493 +    typedef typename Graph::template NodeMap<Status> StatusMap;
3.1494 +    StatusMap* _status;
3.1495 +
3.1496 +    typedef typename Graph::template NodeMap<Arc> PredMap;
3.1497 +    PredMap* _pred;
3.1498 +
3.1499 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
3.1500 +
3.1501 +    IntNodeMap *_tree_set_index;
3.1502 +    TreeSet *_tree_set;
3.1503 +
3.1504 +    IntNodeMap *_delta2_index;
3.1505 +    BinHeap<Value, IntNodeMap> *_delta2;
3.1506 +
3.1507 +    IntEdgeMap *_delta3_index;
3.1508 +    BinHeap<Value, IntEdgeMap> *_delta3;
3.1509 +
3.1510 +    Value _delta_sum;
3.1511 +
3.1512 +    void createStructures() {
3.1513 +      _node_num = countNodes(_graph);
3.1514 +
3.1515 +      if (!_matching) {
3.1516 +        _matching = new MatchingMap(_graph);
3.1517 +      }
3.1518 +      if (!_node_potential) {
3.1519 +        _node_potential = new NodePotential(_graph);
3.1520 +      }
3.1521 +      if (!_status) {
3.1522 +        _status = new StatusMap(_graph);
3.1523 +      }
3.1524 +      if (!_pred) {
3.1525 +        _pred = new PredMap(_graph);
3.1526 +      }
3.1527 +      if (!_tree_set) {
3.1528 +        _tree_set_index = new IntNodeMap(_graph);
3.1529 +        _tree_set = new TreeSet(*_tree_set_index);
3.1530 +      }
3.1531 +      if (!_delta2) {
3.1532 +        _delta2_index = new IntNodeMap(_graph);
3.1533 +        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
3.1534 +      }
3.1535 +      if (!_delta3) {
3.1536 +        _delta3_index = new IntEdgeMap(_graph);
3.1537 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
3.1538 +      }
3.1539 +    }
3.1540 +
3.1541 +    void destroyStructures() {
3.1542 +      if (_matching) {
3.1543 +        delete _matching;
3.1544 +      }
3.1545 +      if (_node_potential) {
3.1546 +        delete _node_potential;
3.1547 +      }
3.1548 +      if (_status) {
3.1549 +        delete _status;
3.1550 +      }
3.1551 +      if (_pred) {
3.1552 +        delete _pred;
3.1553 +      }
3.1554 +      if (_tree_set) {
3.1555 +        delete _tree_set_index;
3.1556 +        delete _tree_set;
3.1557 +      }
3.1558 +      if (_delta2) {
3.1559 +        delete _delta2_index;
3.1560 +        delete _delta2;
3.1561 +      }
3.1562 +      if (_delta3) {
3.1563 +        delete _delta3_index;
3.1564 +        delete _delta3;
3.1565 +      }
3.1566 +    }
3.1567 +
3.1568 +    void matchedToEven(Node node, int tree) {
3.1569 +      _tree_set->insert(node, tree);
3.1570 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.1571 +
3.1572 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
3.1573 +        _delta2->erase(node);
3.1574 +      }
3.1575 +
3.1576 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.1577 +        Node v = _graph.source(a);
3.1578 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.1579 +          dualScale * _weight[a];
3.1580 +        if (node == v) {
3.1581 +          if (_allow_loops && _graph.direction(a)) {
3.1582 +            _delta3->push(a, rw / 2);
3.1583 +          }
3.1584 +        } else if ((*_status)[v] == EVEN) {
3.1585 +          _delta3->push(a, rw / 2);
3.1586 +        } else if ((*_status)[v] == MATCHED) {
3.1587 +          if (_delta2->state(v) != _delta2->IN_HEAP) {
3.1588 +            _pred->set(v, a);
3.1589 +            _delta2->push(v, rw);
3.1590 +          } else if ((*_delta2)[v] > rw) {
3.1591 +            _pred->set(v, a);
3.1592 +            _delta2->decrease(v, rw);
3.1593 +          }
3.1594 +        }
3.1595 +      }
3.1596 +    }
3.1597 +
3.1598 +    void matchedToOdd(Node node, int tree) {
3.1599 +      _tree_set->insert(node, tree);
3.1600 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.1601 +
3.1602 +      if (_delta2->state(node) == _delta2->IN_HEAP) {
3.1603 +        _delta2->erase(node);
3.1604 +      }
3.1605 +    }
3.1606 +
3.1607 +    void evenToMatched(Node node, int tree) {
3.1608 +      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.1609 +      Arc min = INVALID;
3.1610 +      Value minrw = std::numeric_limits<Value>::max();
3.1611 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.1612 +        Node v = _graph.source(a);
3.1613 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.1614 +          dualScale * _weight[a];
3.1615 +
3.1616 +        if (node == v) {
3.1617 +          if (_allow_loops && _graph.direction(a)) {
3.1618 +            _delta3->erase(a);
3.1619 +          }
3.1620 +        } else if ((*_status)[v] == EVEN) {
3.1621 +          _delta3->erase(a);
3.1622 +          if (minrw > rw) {
3.1623 +            min = _graph.oppositeArc(a);
3.1624 +            minrw = rw;
3.1625 +          }
3.1626 +        } else if ((*_status)[v]  == MATCHED) {
3.1627 +          if ((*_pred)[v] == a) {
3.1628 +            Arc mina = INVALID;
3.1629 +            Value minrwa = std::numeric_limits<Value>::max();
3.1630 +            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
3.1631 +              Node va = _graph.target(aa);
3.1632 +              if ((*_status)[va] != EVEN ||
3.1633 +                  _tree_set->find(va) == tree) continue;
3.1634 +              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
3.1635 +                dualScale * _weight[aa];
3.1636 +              if (minrwa > rwa) {
3.1637 +                minrwa = rwa;
3.1638 +                mina = aa;
3.1639 +              }
3.1640 +            }
3.1641 +            if (mina != INVALID) {
3.1642 +              _pred->set(v, mina);
3.1643 +              _delta2->increase(v, minrwa);
3.1644 +            } else {
3.1645 +              _pred->set(v, INVALID);
3.1646 +              _delta2->erase(v);
3.1647 +            }
3.1648 +          }
3.1649 +        }
3.1650 +      }
3.1651 +      if (min != INVALID) {
3.1652 +        _pred->set(node, min);
3.1653 +        _delta2->push(node, minrw);
3.1654 +      } else {
3.1655 +        _pred->set(node, INVALID);
3.1656 +      }
3.1657 +    }
3.1658 +
3.1659 +    void oddToMatched(Node node) {
3.1660 +      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.1661 +      Arc min = INVALID;
3.1662 +      Value minrw = std::numeric_limits<Value>::max();
3.1663 +      for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.1664 +        Node v = _graph.source(a);
3.1665 +        if ((*_status)[v] != EVEN) continue;
3.1666 +        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.1667 +          dualScale * _weight[a];
3.1668 +
3.1669 +        if (minrw > rw) {
3.1670 +          min = _graph.oppositeArc(a);
3.1671 +          minrw = rw;
3.1672 +        }
3.1673 +      }
3.1674 +      if (min != INVALID) {
3.1675 +        _pred->set(node, min);
3.1676 +        _delta2->push(node, minrw);
3.1677 +      } else {
3.1678 +        _pred->set(node, INVALID);
3.1679 +      }
3.1680 +    }
3.1681 +
3.1682 +    void alternatePath(Node even, int tree) {
3.1683 +      Node odd;
3.1684 +
3.1685 +      _status->set(even, MATCHED);
3.1686 +      evenToMatched(even, tree);
3.1687 +
3.1688 +      Arc prev = (*_matching)[even];
3.1689 +      while (prev != INVALID) {
3.1690 +        odd = _graph.target(prev);
3.1691 +        even = _graph.target((*_pred)[odd]);
3.1692 +        _matching->set(odd, (*_pred)[odd]);
3.1693 +        _status->set(odd, MATCHED);
3.1694 +        oddToMatched(odd);
3.1695 +
3.1696 +        prev = (*_matching)[even];
3.1697 +        _status->set(even, MATCHED);
3.1698 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
3.1699 +        evenToMatched(even, tree);
3.1700 +      }
3.1701 +    }
3.1702 +
3.1703 +    void destroyTree(int tree) {
3.1704 +      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
3.1705 +        if ((*_status)[n] == EVEN) {
3.1706 +          _status->set(n, MATCHED);
3.1707 +          evenToMatched(n, tree);
3.1708 +        } else if ((*_status)[n] == ODD) {
3.1709 +          _status->set(n, MATCHED);
3.1710 +          oddToMatched(n);
3.1711 +        }
3.1712 +      }
3.1713 +      _tree_set->eraseClass(tree);
3.1714 +    }
3.1715 +
3.1716 +    void augmentOnEdge(const Edge& edge) {
3.1717 +      Node left = _graph.u(edge);
3.1718 +      int left_tree = _tree_set->find(left);
3.1719 +
3.1720 +      alternatePath(left, left_tree);
3.1721 +      destroyTree(left_tree);
3.1722 +      _matching->set(left, _graph.direct(edge, true));
3.1723 +
3.1724 +      Node right = _graph.v(edge);
3.1725 +      int right_tree = _tree_set->find(right);
3.1726 +
3.1727 +      alternatePath(right, right_tree);
3.1728 +      destroyTree(right_tree);
3.1729 +      _matching->set(right, _graph.direct(edge, false));
3.1730 +    }
3.1731 +
3.1732 +    void augmentOnArc(const Arc& arc) {
3.1733 +      Node left = _graph.source(arc);
3.1734 +      _status->set(left, MATCHED);
3.1735 +      _matching->set(left, arc);
3.1736 +      _pred->set(left, arc);
3.1737 +
3.1738 +      Node right = _graph.target(arc);
3.1739 +      int right_tree = _tree_set->find(right);
3.1740 +
3.1741 +      alternatePath(right, right_tree);
3.1742 +      destroyTree(right_tree);
3.1743 +      _matching->set(right, _graph.oppositeArc(arc));
3.1744 +    }
3.1745 +
3.1746 +    void extendOnArc(const Arc& arc) {
3.1747 +      Node base = _graph.target(arc);
3.1748 +      int tree = _tree_set->find(base);
3.1749 +
3.1750 +      Node odd = _graph.source(arc);
3.1751 +      _tree_set->insert(odd, tree);
3.1752 +      _status->set(odd, ODD);
3.1753 +      matchedToOdd(odd, tree);
3.1754 +      _pred->set(odd, arc);
3.1755 +
3.1756 +      Node even = _graph.target((*_matching)[odd]);
3.1757 +      _tree_set->insert(even, tree);
3.1758 +      _status->set(even, EVEN);
3.1759 +      matchedToEven(even, tree);
3.1760 +    }
3.1761 +
3.1762 +    void cycleOnEdge(const Edge& edge, int tree) {
3.1763 +      Node nca = INVALID;
3.1764 +      std::vector<Node> left_path, right_path;
3.1765 +
3.1766 +      {
3.1767 +        std::set<Node> left_set, right_set;
3.1768 +        Node left = _graph.u(edge);
3.1769 +        left_path.push_back(left);
3.1770 +        left_set.insert(left);
3.1771 +
3.1772 +        Node right = _graph.v(edge);
3.1773 +        right_path.push_back(right);
3.1774 +        right_set.insert(right);
3.1775 +
3.1776 +        while (true) {
3.1777 +
3.1778 +          if (left_set.find(right) != left_set.end()) {
3.1779 +            nca = right;
3.1780 +            break;
3.1781 +          }
3.1782 +
3.1783 +          if ((*_matching)[left] == INVALID) break;
3.1784 +
3.1785 +          left = _graph.target((*_matching)[left]);
3.1786 +          left_path.push_back(left);
3.1787 +          left = _graph.target((*_pred)[left]);
3.1788 +          left_path.push_back(left);
3.1789 +
3.1790 +          left_set.insert(left);
3.1791 +
3.1792 +          if (right_set.find(left) != right_set.end()) {
3.1793 +            nca = left;
3.1794 +            break;
3.1795 +          }
3.1796 +
3.1797 +          if ((*_matching)[right] == INVALID) break;
3.1798 +
3.1799 +          right = _graph.target((*_matching)[right]);
3.1800 +          right_path.push_back(right);
3.1801 +          right = _graph.target((*_pred)[right]);
3.1802 +          right_path.push_back(right);
3.1803 +
3.1804 +          right_set.insert(right);
3.1805 +
3.1806 +        }
3.1807 +
3.1808 +        if (nca == INVALID) {
3.1809 +          if ((*_matching)[left] == INVALID) {
3.1810 +            nca = right;
3.1811 +            while (left_set.find(nca) == left_set.end()) {
3.1812 +              nca = _graph.target((*_matching)[nca]);
3.1813 +              right_path.push_back(nca);
3.1814 +              nca = _graph.target((*_pred)[nca]);
3.1815 +              right_path.push_back(nca);
3.1816 +            }
3.1817 +          } else {
3.1818 +            nca = left;
3.1819 +            while (right_set.find(nca) == right_set.end()) {
3.1820 +              nca = _graph.target((*_matching)[nca]);
3.1821 +              left_path.push_back(nca);
3.1822 +              nca = _graph.target((*_pred)[nca]);
3.1823 +              left_path.push_back(nca);
3.1824 +            }
3.1825 +          }
3.1826 +        }
3.1827 +      }
3.1828 +
3.1829 +      alternatePath(nca, tree);
3.1830 +      Arc prev;
3.1831 +
3.1832 +      prev = _graph.direct(edge, true);
3.1833 +      for (int i = 0; left_path[i] != nca; i += 2) {
3.1834 +        _matching->set(left_path[i], prev);
3.1835 +        _status->set(left_path[i], MATCHED);
3.1836 +        evenToMatched(left_path[i], tree);
3.1837 +
3.1838 +        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
3.1839 +        _status->set(left_path[i + 1], MATCHED);
3.1840 +        oddToMatched(left_path[i + 1]);
3.1841 +      }
3.1842 +      _matching->set(nca, prev);
3.1843 +
3.1844 +      for (int i = 0; right_path[i] != nca; i += 2) {
3.1845 +        _status->set(right_path[i], MATCHED);
3.1846 +        evenToMatched(right_path[i], tree);
3.1847 +
3.1848 +        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
3.1849 +        _status->set(right_path[i + 1], MATCHED);
3.1850 +        oddToMatched(right_path[i + 1]);
3.1851 +      }
3.1852 +
3.1853 +      destroyTree(tree);
3.1854 +    }
3.1855 +
3.1856 +    void extractCycle(const Arc &arc) {
3.1857 +      Node left = _graph.source(arc);
3.1858 +      Node odd = _graph.target((*_matching)[left]);
3.1859 +      Arc prev;
3.1860 +      while (odd != left) {
3.1861 +        Node even = _graph.target((*_matching)[odd]);
3.1862 +        prev = (*_matching)[odd];
3.1863 +        odd = _graph.target((*_matching)[even]);
3.1864 +        _matching->set(even, _graph.oppositeArc(prev));
3.1865 +      }
3.1866 +      _matching->set(left, arc);
3.1867 +
3.1868 +      Node right = _graph.target(arc);
3.1869 +      int right_tree = _tree_set->find(right);
3.1870 +      alternatePath(right, right_tree);
3.1871 +      destroyTree(right_tree);
3.1872 +      _matching->set(right, _graph.oppositeArc(arc));
3.1873 +    }
3.1874 +
3.1875 +  public:
3.1876 +
3.1877 +    /// \brief Constructor
3.1878 +    ///
3.1879 +    /// Constructor.
3.1880 +    MaxWeightedPerfectFractionalMatching(const Graph& graph,
3.1881 +                                         const WeightMap& weight,
3.1882 +                                         bool allow_loops = true)
3.1883 +      : _graph(graph), _weight(weight), _matching(0),
3.1884 +      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
3.1885 +      _status(0),  _pred(0),
3.1886 +      _tree_set_index(0), _tree_set(0),
3.1887 +
3.1888 +      _delta2_index(0), _delta2(0),
3.1889 +      _delta3_index(0), _delta3(0),
3.1890 +
3.1891 +      _delta_sum() {}
3.1892 +
3.1893 +    ~MaxWeightedPerfectFractionalMatching() {
3.1894 +      destroyStructures();
3.1895 +    }
3.1896 +
3.1897 +    /// \name Execution Control
3.1898 +    /// The simplest way to execute the algorithm is to use the
3.1899 +    /// \ref run() member function.
3.1900 +
3.1901 +    ///@{
3.1902 +
3.1903 +    /// \brief Initialize the algorithm
3.1904 +    ///
3.1905 +    /// This function initializes the algorithm.
3.1906 +    void init() {
3.1907 +      createStructures();
3.1908 +
3.1909 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1910 +        (*_delta2_index)[n] = _delta2->PRE_HEAP;
3.1911 +      }
3.1912 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1913 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
3.1914 +      }
3.1915 +
3.1916 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.1917 +        Value max = - std::numeric_limits<Value>::max();
3.1918 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3.1919 +          if (_graph.target(e) == n && !_allow_loops) continue;
3.1920 +          if ((dualScale * _weight[e]) / 2 > max) {
3.1921 +            max = (dualScale * _weight[e]) / 2;
3.1922 +          }
3.1923 +        }
3.1924 +        _node_potential->set(n, max);
3.1925 +
3.1926 +        _tree_set->insert(n);
3.1927 +
3.1928 +        _matching->set(n, INVALID);
3.1929 +        _status->set(n, EVEN);
3.1930 +      }
3.1931 +
3.1932 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1933 +        Node left = _graph.u(e);
3.1934 +        Node right = _graph.v(e);
3.1935 +        if (left == right && !_allow_loops) continue;
3.1936 +        _delta3->push(e, ((*_node_potential)[left] +
3.1937 +                          (*_node_potential)[right] -
3.1938 +                          dualScale * _weight[e]) / 2);
3.1939 +      }
3.1940 +    }
3.1941 +
3.1942 +    /// \brief Start the algorithm
3.1943 +    ///
3.1944 +    /// This function starts the algorithm.
3.1945 +    ///
3.1946 +    /// \pre \ref init() must be called before using this function.
3.1947 +    bool start() {
3.1948 +      enum OpType {
3.1949 +        D2, D3
3.1950 +      };
3.1951 +
3.1952 +      int unmatched = _node_num;
3.1953 +      while (unmatched > 0) {
3.1954 +        Value d2 = !_delta2->empty() ?
3.1955 +          _delta2->prio() : std::numeric_limits<Value>::max();
3.1956 +
3.1957 +        Value d3 = !_delta3->empty() ?
3.1958 +          _delta3->prio() : std::numeric_limits<Value>::max();
3.1959 +
3.1960 +        _delta_sum = d3; OpType ot = D3;
3.1961 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3.1962 +
3.1963 +        if (_delta_sum == std::numeric_limits<Value>::max()) {
3.1964 +          return false;
3.1965 +        }
3.1966 +
3.1967 +        switch (ot) {
3.1968 +        case D2:
3.1969 +          {
3.1970 +            Node n = _delta2->top();
3.1971 +            Arc a = (*_pred)[n];
3.1972 +            if ((*_matching)[n] == INVALID) {
3.1973 +              augmentOnArc(a);
3.1974 +              --unmatched;
3.1975 +            } else {
3.1976 +              Node v = _graph.target((*_matching)[n]);
3.1977 +              if ((*_matching)[n] !=
3.1978 +                  _graph.oppositeArc((*_matching)[v])) {
3.1979 +                extractCycle(a);
3.1980 +                --unmatched;
3.1981 +              } else {
3.1982 +                extendOnArc(a);
3.1983 +              }
3.1984 +            }
3.1985 +          } break;
3.1986 +        case D3:
3.1987 +          {
3.1988 +            Edge e = _delta3->top();
3.1989 +
3.1990 +            Node left = _graph.u(e);
3.1991 +            Node right = _graph.v(e);
3.1992 +
3.1993 +            int left_tree = _tree_set->find(left);
3.1994 +            int right_tree = _tree_set->find(right);
3.1995 +
3.1996 +            if (left_tree == right_tree) {
3.1997 +              cycleOnEdge(e, left_tree);
3.1998 +              --unmatched;
3.1999 +            } else {
3.2000 +              augmentOnEdge(e);
3.2001 +              unmatched -= 2;
3.2002 +            }
3.2003 +          } break;
3.2004 +        }
3.2005 +      }
3.2006 +      return true;
3.2007 +    }
3.2008 +
3.2009 +    /// \brief Run the algorithm.
3.2010 +    ///
3.2011 +    /// This method runs the \c %MaxWeightedMatching algorithm.
3.2012 +    ///
3.2013 +    /// \note mwfm.run() is just a shortcut of the following code.
3.2014 +    /// \code
3.2015 +    ///   mwpfm.init();
3.2016 +    ///   mwpfm.start();
3.2017 +    /// \endcode
3.2018 +    bool run() {
3.2019 +      init();
3.2020 +      return start();
3.2021 +    }
3.2022 +
3.2023 +    /// @}
3.2024 +
3.2025 +    /// \name Primal Solution
3.2026 +    /// Functions to get the primal solution, i.e. the maximum weighted
3.2027 +    /// matching.\n
3.2028 +    /// Either \ref run() or \ref start() function should be called before
3.2029 +    /// using them.
3.2030 +
3.2031 +    /// @{
3.2032 +
3.2033 +    /// \brief Return the weight of the matching.
3.2034 +    ///
3.2035 +    /// This function returns the weight of the found matching. This
3.2036 +    /// value is scaled by \ref primalScale "primal scale".
3.2037 +    ///
3.2038 +    /// \pre Either run() or start() must be called before using this function.
3.2039 +    Value matchingWeight() const {
3.2040 +      Value sum = 0;
3.2041 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.2042 +        if ((*_matching)[n] != INVALID) {
3.2043 +          sum += _weight[(*_matching)[n]];
3.2044 +        }
3.2045 +      }
3.2046 +      return sum * primalScale / 2;
3.2047 +    }
3.2048 +
3.2049 +    /// \brief Return the number of covered nodes in the matching.
3.2050 +    ///
3.2051 +    /// This function returns the number of covered nodes in the matching.
3.2052 +    ///
3.2053 +    /// \pre Either run() or start() must be called before using this function.
3.2054 +    int matchingSize() const {
3.2055 +      int num = 0;
3.2056 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.2057 +        if ((*_matching)[n] != INVALID) {
3.2058 +          ++num;
3.2059 +        }
3.2060 +      }
3.2061 +      return num;
3.2062 +    }
3.2063 +
3.2064 +    /// \brief Return \c true if the given edge is in the matching.
3.2065 +    ///
3.2066 +    /// This function returns \c true if the given edge is in the
3.2067 +    /// found matching. The result is scaled by \ref primalScale
3.2068 +    /// "primal scale".
3.2069 +    ///
3.2070 +    /// \pre Either run() or start() must be called before using this function.
3.2071 +    Value matching(const Edge& edge) const {
3.2072 +      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
3.2073 +        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
3.2074 +        * primalScale / 2;
3.2075 +    }
3.2076 +
3.2077 +    /// \brief Return the fractional matching arc (or edge) incident
3.2078 +    /// to the given node.
3.2079 +    ///
3.2080 +    /// This function returns one of the fractional matching arc (or
3.2081 +    /// edge) incident to the given node in the found matching or \c
3.2082 +    /// INVALID if the node is not covered by the matching or if the
3.2083 +    /// node is on an odd length cycle then it is the successor edge
3.2084 +    /// on the cycle.
3.2085 +    ///
3.2086 +    /// \pre Either run() or start() must be called before using this function.
3.2087 +    Arc matching(const Node& node) const {
3.2088 +      return (*_matching)[node];
3.2089 +    }
3.2090 +
3.2091 +    /// \brief Return a const reference to the matching map.
3.2092 +    ///
3.2093 +    /// This function returns a const reference to a node map that stores
3.2094 +    /// the matching arc (or edge) incident to each node.
3.2095 +    const MatchingMap& matchingMap() const {
3.2096 +      return *_matching;
3.2097 +    }
3.2098 +
3.2099 +    /// @}
3.2100 +
3.2101 +    /// \name Dual Solution
3.2102 +    /// Functions to get the dual solution.\n
3.2103 +    /// Either \ref run() or \ref start() function should be called before
3.2104 +    /// using them.
3.2105 +
3.2106 +    /// @{
3.2107 +
3.2108 +    /// \brief Return the value of the dual solution.
3.2109 +    ///
3.2110 +    /// This function returns the value of the dual solution.
3.2111 +    /// It should be equal to the primal value scaled by \ref dualScale
3.2112 +    /// "dual scale".
3.2113 +    ///
3.2114 +    /// \pre Either run() or start() must be called before using this function.
3.2115 +    Value dualValue() const {
3.2116 +      Value sum = 0;
3.2117 +      for (NodeIt n(_graph); n != INVALID; ++n) {
3.2118 +        sum += nodeValue(n);
3.2119 +      }
3.2120 +      return sum;
3.2121 +    }
3.2122 +
3.2123 +    /// \brief Return the dual value (potential) of the given node.
3.2124 +    ///
3.2125 +    /// This function returns the dual value (potential) of the given node.
3.2126 +    ///
3.2127 +    /// \pre Either run() or start() must be called before using this function.
3.2128 +    Value nodeValue(const Node& n) const {
3.2129 +      return (*_node_potential)[n];
3.2130 +    }
3.2131 +
3.2132 +    /// @}
3.2133 +
3.2134 +  };
3.2135 +
3.2136 +} //END OF NAMESPACE LEMON
3.2137 +
3.2138 +#endif //LEMON_FRACTIONAL_MATCHING_H

     4.1 --- a/test/CMakeLists.txt	Sun Sep 20 21:38:24 2009 +0200
4.2 +++ b/test/CMakeLists.txt	Fri Sep 25 21:51:36 2009 +0200
4.3 @@ -21,6 +21,7 @@
4.4    edge_set_test
4.5    error_test
4.6    euler_test
4.7 +  fractional_matching_test
4.8    gomory_hu_test
4.9    graph_copy_test
4.10    graph_test

     5.1 --- a/test/Makefile.am	Sun Sep 20 21:38:24 2009 +0200
5.2 +++ b/test/Makefile.am	Fri Sep 25 21:51:36 2009 +0200
5.3 @@ -19,6 +19,7 @@
5.4  	test/edge_set_test \
5.5  	test/error_test \
5.6  	test/euler_test \
5.7 +	test/fractional_matching_test \
5.8  	test/gomory_hu_test \
5.9  	test/graph_copy_test \
5.10  	test/graph_test \
5.11 @@ -65,6 +66,7 @@
5.12  test_edge_set_test_SOURCES = test/edge_set_test.cc
5.13  test_error_test_SOURCES = test/error_test.cc
5.14  test_euler_test_SOURCES = test/euler_test.cc
5.15 +test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
5.16  test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
5.17  test_graph_copy_test_SOURCES = test/graph_copy_test.cc
5.18  test_graph_test_SOURCES = test/graph_test.cc

     6.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/test/fractional_matching_test.cc	Fri Sep 25 21:51:36 2009 +0200
6.3 @@ -0,0 +1,502 @@
6.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
6.5 + *
6.6 + * This file is a part of LEMON, a generic C++ optimization library.
6.7 + *
6.8 + * Copyright (C) 2003-2009
6.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
6.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
6.11 + *
6.12 + * Permission to use, modify and distribute this software is granted
6.13 + * provided that this copyright notice appears in all copies. For
6.14 + * precise terms see the accompanying LICENSE file.
6.15 + *
6.16 + * This software is provided "AS IS" with no warranty of any kind,
6.17 + * express or implied, and with no claim as to its suitability for any
6.18 + * purpose.
6.19 + *
6.20 + */
6.21 +
6.22 +#include <iostream>
6.23 +#include <sstream>
6.24 +#include <vector>
6.25 +#include <queue>
6.26 +#include <cstdlib>
6.27 +
6.28 +#include <lemon/fractional_matching.h>
6.29 +#include <lemon/smart_graph.h>
6.30 +#include <lemon/concepts/graph.h>
6.31 +#include <lemon/concepts/maps.h>
6.33 +#include <lemon/math.h>
6.34 +
6.35 +#include "test_tools.h"
6.36 +
6.37 +using namespace std;
6.38 +using namespace lemon;
6.39 +
6.40 +GRAPH_TYPEDEFS(SmartGraph);
6.41 +
6.42 +
6.43 +const int lgfn = 4;
6.44 +const std::string lgf[lgfn] = {
6.45 +  "@nodes\n"
6.46 +  "label\n"
6.47 +  "0\n"
6.48 +  "1\n"
6.49 +  "2\n"
6.50 +  "3\n"
6.51 +  "4\n"
6.52 +  "5\n"
6.53 +  "6\n"
6.54 +  "7\n"
6.55 +  "@edges\n"
6.56 +  "     label  weight\n"
6.57 +  "7 4  0      984\n"
6.58 +  "0 7  1      73\n"
6.59 +  "7 1  2      204\n"
6.60 +  "2 3  3      583\n"
6.61 +  "2 7  4      565\n"
6.62 +  "2 1  5      582\n"
6.63 +  "0 4  6      551\n"
6.64 +  "2 5  7      385\n"
6.65 +  "1 5  8      561\n"
6.66 +  "5 3  9      484\n"
6.67 +  "7 5  10     904\n"
6.68 +  "3 6  11     47\n"
6.69 +  "7 6  12     888\n"
6.70 +  "3 0  13     747\n"
6.71 +  "6 1  14     310\n",
6.72 +
6.73 +  "@nodes\n"
6.74 +  "label\n"
6.75 +  "0\n"
6.76 +  "1\n"
6.77 +  "2\n"
6.78 +  "3\n"
6.79 +  "4\n"
6.80 +  "5\n"
6.81 +  "6\n"
6.82 +  "7\n"
6.83 +  "@edges\n"
6.84 +  "     label  weight\n"
6.85 +  "2 5  0      710\n"
6.86 +  "0 5  1      241\n"
6.87 +  "2 4  2      856\n"
6.88 +  "2 6  3      762\n"
6.89 +  "4 1  4      747\n"
6.90 +  "6 1  5      962\n"
6.91 +  "4 7  6      723\n"
6.92 +  "1 7  7      661\n"
6.93 +  "2 3  8      376\n"
6.94 +  "1 0  9      416\n"
6.95 +  "6 7  10     391\n",
6.96 +
6.97 +  "@nodes\n"
6.98 +  "label\n"
6.99 +  "0\n"
6.100 +  "1\n"
6.101 +  "2\n"
6.102 +  "3\n"
6.103 +  "4\n"
6.104 +  "5\n"
6.105 +  "6\n"
6.106 +  "7\n"
6.107 +  "@edges\n"
6.108 +  "     label  weight\n"
6.109 +  "6 2  0      553\n"
6.110 +  "0 7  1      653\n"
6.111 +  "6 3  2      22\n"
6.112 +  "4 7  3      846\n"
6.113 +  "7 2  4      981\n"
6.114 +  "7 6  5      250\n"
6.115 +  "5 2  6      539\n",
6.116 +
6.117 +  "@nodes\n"
6.118 +  "label\n"
6.119 +  "0\n"
6.120 +  "@edges\n"
6.121 +  "     label  weight\n"
6.122 +  "0 0  0      100\n"
6.123 +};
6.124 +
6.125 +void checkMaxFractionalMatchingCompile()
6.126 +{
6.127 +  typedef concepts::Graph Graph;
6.128 +  typedef Graph::Node Node;
6.129 +  typedef Graph::Edge Edge;
6.130 +
6.131 +  Graph g;
6.132 +  Node n;
6.133 +  Edge e;
6.134 +
6.135 +  MaxFractionalMatching<Graph> mat_test(g);
6.136 +  const MaxFractionalMatching<Graph>&
6.137 +    const_mat_test = mat_test;
6.138 +
6.139 +  mat_test.init();
6.140 +  mat_test.start();
6.141 +  mat_test.start(true);
6.142 +  mat_test.startPerfect();
6.143 +  mat_test.startPerfect(true);
6.144 +  mat_test.run();
6.145 +  mat_test.run(true);
6.146 +  mat_test.runPerfect();
6.147 +  mat_test.runPerfect(true);
6.148 +
6.149 +  const_mat_test.matchingSize();
6.150 +  const_mat_test.matching(e);
6.151 +  const_mat_test.matching(n);
6.152 +  const MaxFractionalMatching<Graph>::MatchingMap& mmap =
6.153 +    const_mat_test.matchingMap();
6.154 +  e = mmap[n];
6.155 +
6.156 +  const_mat_test.barrier(n);
6.157 +}
6.158 +
6.159 +void checkMaxWeightedFractionalMatchingCompile()
6.160 +{
6.161 +  typedef concepts::Graph Graph;
6.162 +  typedef Graph::Node Node;
6.163 +  typedef Graph::Edge Edge;
6.164 +  typedef Graph::EdgeMap<int> WeightMap;
6.165 +
6.166 +  Graph g;
6.167 +  Node n;
6.168 +  Edge e;
6.169 +  WeightMap w(g);
6.170 +
6.171 +  MaxWeightedFractionalMatching<Graph> mat_test(g, w);
6.172 +  const MaxWeightedFractionalMatching<Graph>&
6.173 +    const_mat_test = mat_test;
6.174 +
6.175 +  mat_test.init();
6.176 +  mat_test.start();
6.177 +  mat_test.run();
6.178 +
6.179 +  const_mat_test.matchingWeight();
6.180 +  const_mat_test.matchingSize();
6.181 +  const_mat_test.matching(e);
6.182 +  const_mat_test.matching(n);
6.183 +  const MaxWeightedFractionalMatching<Graph>::MatchingMap& mmap =
6.184 +    const_mat_test.matchingMap();
6.185 +  e = mmap[n];
6.186 +
6.187 +  const_mat_test.dualValue();
6.188 +  const_mat_test.nodeValue(n);
6.189 +}
6.190 +
6.191 +void checkMaxWeightedPerfectFractionalMatchingCompile()
6.192 +{
6.193 +  typedef concepts::Graph Graph;
6.194 +  typedef Graph::Node Node;
6.195 +  typedef Graph::Edge Edge;
6.196 +  typedef Graph::EdgeMap<int> WeightMap;
6.197 +
6.198 +  Graph g;
6.199 +  Node n;
6.200 +  Edge e;
6.201 +  WeightMap w(g);
6.202 +
6.203 +  MaxWeightedPerfectFractionalMatching<Graph> mat_test(g, w);
6.204 +  const MaxWeightedPerfectFractionalMatching<Graph>&
6.205 +    const_mat_test = mat_test;
6.206 +
6.207 +  mat_test.init();
6.208 +  mat_test.start();
6.209 +  mat_test.run();
6.210 +
6.211 +  const_mat_test.matchingWeight();
6.212 +  const_mat_test.matching(e);
6.213 +  const_mat_test.matching(n);
6.214 +  const MaxWeightedPerfectFractionalMatching<Graph>::MatchingMap& mmap =
6.215 +    const_mat_test.matchingMap();
6.216 +  e = mmap[n];
6.217 +
6.218 +  const_mat_test.dualValue();
6.219 +  const_mat_test.nodeValue(n);
6.220 +}
6.221 +
6.222 +void checkFractionalMatching(const SmartGraph& graph,
6.223 +                             const MaxFractionalMatching<SmartGraph>& mfm,
6.224 +                             bool allow_loops = true) {
6.225 +  int pv = 0;
6.226 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.227 +    int indeg = 0;
6.228 +    for (InArcIt a(graph, n); a != INVALID; ++a) {
6.229 +      if (mfm.matching(graph.source(a)) == a) {
6.230 +        ++indeg;
6.231 +      }
6.232 +    }
6.233 +    if (mfm.matching(n) != INVALID) {
6.234 +      check(indeg == 1, "Invalid matching");
6.235 +      ++pv;
6.236 +    } else {
6.237 +      check(indeg == 0, "Invalid matching");
6.238 +    }
6.239 +  }
6.240 +  check(pv == mfm.matchingSize(), "Wrong matching size");
6.241 +
6.242 +  SmartGraph::NodeMap<bool> processed(graph, false);
6.243 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.244 +    if (processed[n]) continue;
6.245 +    processed[n] = true;
6.246 +    if (mfm.matching(n) == INVALID) continue;
6.247 +    int num = 1;
6.248 +    Node v = graph.target(mfm.matching(n));
6.249 +    while (v != n) {
6.250 +      processed[v] = true;
6.251 +      ++num;
6.252 +      v = graph.target(mfm.matching(v));
6.253 +    }
6.254 +    check(num == 2 || num % 2 == 1, "Wrong cycle size");
6.255 +    check(allow_loops || num != 1, "Wrong cycle size");
6.256 +  }
6.257 +
6.258 +  int anum = 0, bnum = 0;
6.259 +  SmartGraph::NodeMap<bool> neighbours(graph, false);
6.260 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.261 +    if (!mfm.barrier(n)) continue;
6.262 +    ++anum;
6.263 +    for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
6.264 +      Node u = graph.source(a);
6.265 +      if (!allow_loops && u == n) continue;
6.266 +      if (!neighbours[u]) {
6.267 +        neighbours[u] = true;
6.268 +        ++bnum;
6.269 +      }
6.270 +    }
6.271 +  }
6.272 +  check(anum - bnum + mfm.matchingSize() == countNodes(graph),
6.273 +        "Wrong barrier");
6.274 +}
6.275 +
6.276 +void checkPerfectFractionalMatching(const SmartGraph& graph,
6.277 +                             const MaxFractionalMatching<SmartGraph>& mfm,
6.278 +                             bool perfect, bool allow_loops = true) {
6.279 +  if (perfect) {
6.280 +    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.281 +      int indeg = 0;
6.282 +      for (InArcIt a(graph, n); a != INVALID; ++a) {
6.283 +        if (mfm.matching(graph.source(a)) == a) {
6.284 +          ++indeg;
6.285 +        }
6.286 +      }
6.287 +      check(mfm.matching(n) != INVALID, "Invalid matching");
6.288 +      check(indeg == 1, "Invalid matching");
6.289 +    }
6.290 +  } else {
6.291 +    int anum = 0, bnum = 0;
6.292 +    SmartGraph::NodeMap<bool> neighbours(graph, false);
6.293 +    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.294 +      if (!mfm.barrier(n)) continue;
6.295 +      ++anum;
6.296 +      for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
6.297 +        Node u = graph.source(a);
6.298 +        if (!allow_loops && u == n) continue;
6.299 +        if (!neighbours[u]) {
6.300 +          neighbours[u] = true;
6.301 +          ++bnum;
6.302 +        }
6.303 +      }
6.304 +    }
6.305 +    check(anum - bnum > 0, "Wrong barrier");
6.306 +  }
6.307 +}
6.308 +
6.309 +void checkWeightedFractionalMatching(const SmartGraph& graph,
6.310 +                   const SmartGraph::EdgeMap<int>& weight,
6.311 +                   const MaxWeightedFractionalMatching<SmartGraph>& mwfm,
6.312 +                   bool allow_loops = true) {
6.313 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
6.314 +    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
6.315 +    int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e))
6.316 +      - weight[e] * mwfm.dualScale;
6.317 +
6.318 +    check(rw >= 0, "Negative reduced weight");
6.319 +    check(rw == 0 || !mwfm.matching(e),
6.320 +          "Non-zero reduced weight on matching edge");
6.321 +  }
6.322 +
6.323 +  int pv = 0;
6.324 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.325 +    int indeg = 0;
6.326 +    for (InArcIt a(graph, n); a != INVALID; ++a) {
6.327 +      if (mwfm.matching(graph.source(a)) == a) {
6.328 +        ++indeg;
6.329 +      }
6.330 +    }
6.331 +    check(indeg <= 1, "Invalid matching");
6.332 +    if (mwfm.matching(n) != INVALID) {
6.333 +      check(mwfm.nodeValue(n) >= 0, "Invalid node value");
6.334 +      check(indeg == 1, "Invalid matching");
6.335 +      pv += weight[mwfm.matching(n)];
6.336 +      SmartGraph::Node o = graph.target(mwfm.matching(n));
6.337 +    } else {
6.338 +      check(mwfm.nodeValue(n) == 0, "Invalid matching");
6.339 +      check(indeg == 0, "Invalid matching");
6.340 +    }
6.341 +  }
6.342 +
6.343 +  int dv = 0;
6.344 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.345 +    dv += mwfm.nodeValue(n);
6.346 +  }
6.347 +
6.348 +  check(pv * mwfm.dualScale == dv * 2, "Wrong duality");
6.349 +
6.350 +  SmartGraph::NodeMap<bool> processed(graph, false);
6.351 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.352 +    if (processed[n]) continue;
6.353 +    processed[n] = true;
6.354 +    if (mwfm.matching(n) == INVALID) continue;
6.355 +    int num = 1;
6.356 +    Node v = graph.target(mwfm.matching(n));
6.357 +    while (v != n) {
6.358 +      processed[v] = true;
6.359 +      ++num;
6.360 +      v = graph.target(mwfm.matching(v));
6.361 +    }
6.362 +    check(num == 2 || num % 2 == 1, "Wrong cycle size");
6.363 +    check(allow_loops || num != 1, "Wrong cycle size");
6.364 +  }
6.365 +
6.366 +  return;
6.367 +}
6.368 +
6.369 +void checkWeightedPerfectFractionalMatching(const SmartGraph& graph,
6.370 +                const SmartGraph::EdgeMap<int>& weight,
6.371 +                const MaxWeightedPerfectFractionalMatching<SmartGraph>& mwpfm,
6.372 +                bool allow_loops = true) {
6.373 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
6.374 +    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
6.375 +    int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e))
6.376 +      - weight[e] * mwpfm.dualScale;
6.377 +
6.378 +    check(rw >= 0, "Negative reduced weight");
6.379 +    check(rw == 0 || !mwpfm.matching(e),
6.380 +          "Non-zero reduced weight on matching edge");
6.381 +  }
6.382 +
6.383 +  int pv = 0;
6.384 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.385 +    int indeg = 0;
6.386 +    for (InArcIt a(graph, n); a != INVALID; ++a) {
6.387 +      if (mwpfm.matching(graph.source(a)) == a) {
6.388 +        ++indeg;
6.389 +      }
6.390 +    }
6.391 +    check(mwpfm.matching(n) != INVALID, "Invalid perfect matching");
6.392 +    check(indeg == 1, "Invalid perfect matching");
6.393 +    pv += weight[mwpfm.matching(n)];
6.394 +    SmartGraph::Node o = graph.target(mwpfm.matching(n));
6.395 +  }
6.396 +
6.397 +  int dv = 0;
6.398 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.399 +    dv += mwpfm.nodeValue(n);
6.400 +  }
6.401 +
6.402 +  check(pv * mwpfm.dualScale == dv * 2, "Wrong duality");
6.403 +
6.404 +  SmartGraph::NodeMap<bool> processed(graph, false);
6.405 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
6.406 +    if (processed[n]) continue;
6.407 +    processed[n] = true;
6.408 +    if (mwpfm.matching(n) == INVALID) continue;
6.409 +    int num = 1;
6.410 +    Node v = graph.target(mwpfm.matching(n));
6.411 +    while (v != n) {
6.412 +      processed[v] = true;
6.413 +      ++num;
6.414 +      v = graph.target(mwpfm.matching(v));
6.415 +    }
6.416 +    check(num == 2 || num % 2 == 1, "Wrong cycle size");
6.417 +    check(allow_loops || num != 1, "Wrong cycle size");
6.418 +  }
6.419 +
6.420 +  return;
6.421 +}
6.422 +
6.423 +
6.424 +int main() {
6.425 +
6.426 +  for (int i = 0; i < lgfn; ++i) {
6.427 +    SmartGraph graph;
6.428 +    SmartGraph::EdgeMap<int> weight(graph);
6.429 +
6.430 +    istringstream lgfs(lgf[i]);
6.432 +      edgeMap("weight", weight).run();
6.433 +
6.434 +    bool perfect_with_loops;
6.435 +    {
6.436 +      MaxFractionalMatching<SmartGraph> mfm(graph, true);
6.437 +      mfm.run();
6.438 +      checkFractionalMatching(graph, mfm, true);
6.439 +      perfect_with_loops = mfm.matchingSize() == countNodes(graph);
6.440 +    }
6.441 +
6.442 +    bool perfect_without_loops;
6.443 +    {
6.444 +      MaxFractionalMatching<SmartGraph> mfm(graph, false);
6.445 +      mfm.run();
6.446 +      checkFractionalMatching(graph, mfm, false);
6.447 +      perfect_without_loops = mfm.matchingSize() == countNodes(graph);
6.448 +    }
6.449 +
6.450 +    {
6.451 +      MaxFractionalMatching<SmartGraph> mfm(graph, true);
6.452 +      bool result = mfm.runPerfect();
6.453 +      checkPerfectFractionalMatching(graph, mfm, result, true);
6.454 +      check(result == perfect_with_loops, "Wrong perfect matching");
6.455 +    }
6.456 +
6.457 +    {
6.458 +      MaxFractionalMatching<SmartGraph> mfm(graph, false);
6.459 +      bool result = mfm.runPerfect();
6.460 +      checkPerfectFractionalMatching(graph, mfm, result, false);
6.461 +      check(result == perfect_without_loops, "Wrong perfect matching");
6.462 +    }
6.463 +
6.464 +    {
6.465 +      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, true);
6.466 +      mwfm.run();
6.467 +      checkWeightedFractionalMatching(graph, weight, mwfm, true);
6.468 +    }
6.469 +
6.470 +    {
6.471 +      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, false);
6.472 +      mwfm.run();
6.473 +      checkWeightedFractionalMatching(graph, weight, mwfm, false);
6.474 +    }
6.475 +
6.476 +    {
6.477 +      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
6.478 +                                                             true);
6.479 +      bool perfect = mwpfm.run();
6.480 +      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
6.481 +            "Perfect matching found");
6.482 +      check(perfect == perfect_with_loops, "Wrong perfect matching");
6.483 +
6.484 +      if (perfect) {
6.485 +        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true);
6.486 +      }
6.487 +    }
6.488 +
6.489 +    {
6.490 +      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
6.491 +                                                             false);
6.492 +      bool perfect = mwpfm.run();
6.493 +      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
6.494 +            "Perfect matching found");
6.495 +      check(perfect == perfect_without_loops, "Wrong perfect matching");
6.496 +
6.497 +      if (perfect) {
6.498 +        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false);
6.499 +      }
6.500 +    }
6.501 +
6.502 +  }
6.503 +
6.504 +  return 0;
6.505 +}