LEMON/LEMON-main Changeset - r874:d8ea85825e02

Repositories » LEMON » LEMON-main

Summary
Changelog
Files
Switch To
- loading...
Options
- Lightweight changelog
- Search
Pull Requests

Changeset - r874:d8ea85825e02

« 872:41d7ac
« 873:23da1b

875:07ec2b »

r874:d8ea85825e02

Tue, 16 Mar 2010 21:27:35

[Not Reviewed]

0 comments (0 inline)

alpar (Alpar Juttner)
alpar@cs.elte.hu

Merge #314

0 6 2

merge default

4 files changed with 3044 insertions and 176 deletions:

lemon/fractional_matching.h

2130

test/fractional_matching_test.cc

525

doc/groups.dox

lemon/Makefile.am

lemon/matching.h

350

172

test/CMakeLists.txt

test/Makefile.am

test/matching_test.cc

↑ Collapse diff ↑

lemon/fractional_matching.h

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#ifndef LEMON_FRACTIONAL_MATCHING_H
 		#define LEMON_FRACTIONAL_MATCHING_H
 		#include <vector>
 		#include <queue>
 		#include <set>
 		#include <limits>
 		#include <lemon/core.h>
 		#include <lemon/unionfind.h>
 		#include <lemon/bin_heap.h>
 		#include <lemon/maps.h>
 		#include <lemon/assert.h>
 		#include <lemon/elevator.h>
 		///\ingroup matching
 		///\file
 		///\brief Fractional matching algorithms in general graphs.
 		namespace lemon {
 		  /// \brief Default traits class of MaxFractionalMatching class.
 		  ///
 		  /// Default traits class of MaxFractionalMatching class.
 		  /// \tparam GR Graph type.
 		  template <typename GR>
 		  struct MaxFractionalMatchingDefaultTraits {
 		    /// \brief The type of the graph the algorithm runs on.
 		    typedef GR Graph;
 		    /// \brief The type of the map that stores the matching.
 		    ///
 		    /// The type of the map that stores the matching arcs.
 		    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
 		    typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
 		    /// \brief Instantiates a MatchingMap.
 		    ///
 		    /// This function instantiates a \ref MatchingMap.
 		    /// \param graph The graph for which we would like to define
 		    /// the matching map.
 		    static MatchingMap* createMatchingMap(const Graph& graph) {
 		      return new MatchingMap(graph);
+		    }
 		    /// \brief The elevator type used by MaxFractionalMatching algorithm.
 		    ///
 		    /// The elevator type used by MaxFractionalMatching algorithm.
 		    ///
 		    /// \sa Elevator
 		    /// \sa LinkedElevator
 		    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
 		    /// \brief Instantiates an Elevator.
 		    ///
 		    /// This function instantiates an \ref Elevator.
 		    /// \param graph The graph for which we would like to define
 		    /// the elevator.
 		    /// \param max_level The maximum level of the elevator.
 		    static Elevator* createElevator(const Graph& graph, int max_level) {
 		      return new Elevator(graph, max_level);
+		    }
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Max cardinality fractional matching
 		  ///
 		  /// This class provides an implementation of fractional matching
 		  /// algorithm based on push-relabel principle.
 		  ///
 		  /// The maximum cardinality fractional matching is a relaxation of the
 		  /// maximum cardinality matching problem where the odd set constraints
 		  /// are omitted.
 		  /// It can be formulated with the following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
 		  /// \f[x_e \ge 0\quad \forall e\in E\f]
 		  /// \f[\max \sum_{e\in E}x_e\f]
 		  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
 		  /// \f$X\f$. The result can be represented as the union of a
 		  /// matching with one value edges and a set of odd length cycles
 		  /// with half value edges.
 		  ///
 		  /// The algorithm calculates an optimal fractional matching and a
 		  /// barrier. The number of adjacents of any node set minus the size
 		  /// of node set is a lower bound on the uncovered nodes in the
 		  /// graph. For maximum matching a barrier is computed which
 		  /// maximizes this difference.
 		  ///
 		  /// The algorithm can be executed with the run() function.  After it
 		  /// the matching (the primal solution) and the barrier (the dual
 		  /// solution) can be obtained using the query functions.
 		  ///
 		  /// The primal solution is multiplied by
 		  /// \ref MaxFractionalMatching::primalScale "2".
 		  ///
 		  /// \tparam GR The undirected graph type the algorithm runs on.
 		#ifdef DOXYGEN
 		  template <typename GR, typename TR>
 		#else
 		  template <typename GR,
 		            typename TR = MaxFractionalMatchingDefaultTraits<GR> >
 		#endif
 		  class MaxFractionalMatching {
 		  public:
 		    /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
 		    /// class" of the algorithm.
 		    typedef TR Traits;
 		    /// The type of the graph the algorithm runs on.
 		    typedef typename TR::Graph Graph;
 		    /// The type of the matching map.
 		    typedef typename TR::MatchingMap MatchingMap;
 		    /// The type of the elevator.
 		    typedef typename TR::Elevator Elevator;
 		    /// \brief Scaling factor for primal solution
 		    ///
 		    /// Scaling factor for primal solution.
 		    static const int primalScale = 2;
 		  private:
 		    const Graph &_graph;
 		    int _node_num;
 		    bool _allow_loops;
 		    int _empty_level;
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    bool _local_matching;
 		    MatchingMap *_matching;
 		    bool _local_level;
 		    Elevator *_level;
 		    typedef typename Graph::template NodeMap<int> InDegMap;
 		    InDegMap *_indeg;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      if (!_matching) {
 		        _local_matching = true;
 		        _matching = Traits::createMatchingMap(_graph);
+		      }
 		      if (!_level) {
 		        _local_level = true;
 		        _level = Traits::createElevator(_graph, _node_num);
+		      }
 		      if (!_indeg) {
 		        _indeg = new InDegMap(_graph);
+		      }
+		    }
 		    void destroyStructures() {
 		      if (_local_matching) {
 		        delete _matching;
+		      }
 		      if (_local_level) {
 		        delete _level;
+		      }
 		      if (_indeg) {
 		        delete _indeg;
+		      }
+		    }
 		    void postprocessing() {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_indeg)[n] != 0) continue;
 		        _indeg->set(n, -1);
 		        Node u = n;
 		        while ((*_matching)[u] != INVALID) {
 		          Node v = _graph.target((*_matching)[u]);
 		          _indeg->set(v, -1);
 		          Arc a = _graph.oppositeArc((*_matching)[u]);
 		          u = _graph.target((*_matching)[v]);
 		          _indeg->set(u, -1);
 		          _matching->set(v, a);
+		        }
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_indeg)[n] != 1) continue;
 		        _indeg->set(n, -1);
 		        int num = 1;
 		        Node u = _graph.target((*_matching)[n]);
 		        while (u != n) {
 		          _indeg->set(u, -1);
 		          u = _graph.target((*_matching)[u]);
 		          ++num;
+		        }
 		        if (num % 2 == 0 && num > 2) {
 		          Arc prev = _graph.oppositeArc((*_matching)[n]);
 		          Node v = _graph.target((*_matching)[n]);
 		          u = _graph.target((*_matching)[v]);
 		          _matching->set(v, prev);
 		          while (u != n) {
 		            prev = _graph.oppositeArc((*_matching)[u]);
 		            v = _graph.target((*_matching)[u]);
 		            u = _graph.target((*_matching)[v]);
 		            _matching->set(v, prev);
+		          }
+		        }
+		      }
+		    }
 		  public:
 		    typedef MaxFractionalMatching Create;
 		    ///\name Named Template Parameters
 		    ///@{
 		    template <typename T>
 		    struct SetMatchingMapTraits : public Traits {
 		      typedef T MatchingMap;
 		      static MatchingMap *createMatchingMap(const Graph&) {
 		        LEMON_ASSERT(false, "MatchingMap is not initialized");
 		        return 0; // ignore warnings
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// MatchingMap type
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting MatchingMap
 		    /// type.
 		    template <typename T>
 		    struct SetMatchingMap
 		      : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
 		      typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
 		    };
 		    template <typename T>
 		    struct SetElevatorTraits : public Traits {
 		      typedef T Elevator;
 		      static Elevator *createElevator(const Graph&, int) {
 		        LEMON_ASSERT(false, "Elevator is not initialized");
 		        return 0; // ignore warnings
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// Elevator type
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting Elevator
 		    /// type. If this named parameter is used, then an external
 		    /// elevator object must be passed to the algorithm using the
 		    /// \ref elevator(Elevator&) "elevator()" function before calling
 		    /// \ref run() or \ref init().
 		    /// \sa SetStandardElevator
 		    template <typename T>
 		    struct SetElevator
 		      : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
 		      typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
 		    };
 		    template <typename T>
 		    struct SetStandardElevatorTraits : public Traits {
 		      typedef T Elevator;
 		      static Elevator *createElevator(const Graph& graph, int max_level) {
 		        return new Elevator(graph, max_level);
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// Elevator type with automatic allocation
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting Elevator
 		    /// type with automatic allocation.
 		    /// The Elevator should have standard constructor interface to be
 		    /// able to automatically created by the algorithm (i.e. the
 		    /// graph and the maximum level should be passed to it).
 		    /// However an external elevator object could also be passed to the
 		    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
 		    /// before calling \ref run() or \ref init().
 		    /// \sa SetElevator
 		    template <typename T>
 		    struct SetStandardElevator
 		      : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
 		      typedef MaxFractionalMatching<Graph,
 		                                    SetStandardElevatorTraits<T> > Create;
 		    };
 		    /// @}
 		  protected:
 		    MaxFractionalMatching() {}
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    ///
 		    MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
 		      : _graph(graph), _allow_loops(allow_loops),
 		        _local_matching(false), _matching(0),
 		        _local_level(false), _level(0),  _indeg(0)
 		    {}
 		    ~MaxFractionalMatching() {
 		      destroyStructures();
+		    }
 		    /// \brief Sets the matching map.
 		    ///
 		    /// Sets the matching map.
 		    /// If you don't use this function before calling \ref run() or
 		    /// \ref init(), an instance will be allocated automatically.
 		    /// The destructor deallocates this automatically allocated map,
 		    /// of course.
 		    /// \return <tt>(*this)</tt>
 		    MaxFractionalMatching& matchingMap(MatchingMap& map) {
 		      if (_local_matching) {
 		        delete _matching;
 		        _local_matching = false;
+		      }
 		      _matching = &map;
 		      return *this;
+		    }
 		    /// \brief Sets the elevator used by algorithm.
 		    ///
 		    /// Sets the elevator used by algorithm.
 		    /// If you don't use this function before calling \ref run() or
 		    /// \ref init(), an instance will be allocated automatically.
 		    /// The destructor deallocates this automatically allocated elevator,
 		    /// of course.
 		    /// \return <tt>(*this)</tt>
 		    MaxFractionalMatching& elevator(Elevator& elevator) {
 		      if (_local_level) {
 		        delete _level;
 		        _local_level = false;
+		      }
 		      _level = &elevator;
 		      return *this;
+		    }
 		    /// \brief Returns a const reference to the elevator.
 		    ///
 		    /// Returns a const reference to the elevator.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    const Elevator& elevator() const {
 		      return *_level;
+		    }
 		    /// \name Execution control
 		    /// The simplest way to execute the algorithm is to use one of the
 		    /// member functions called \c run(). \n
 		    /// If you need more control on the execution, first
 		    /// you must call \ref init() and then one variant of the start()
 		    /// member.
 		    /// @{
 		    /// \brief Initializes the internal data structures.
 		    ///
 		    /// Initializes the internal data structures and sets the initial
 		    /// matching.
 		    void init() {
 		      createStructures();
 		      _level->initStart();
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        _indeg->set(n, 0);
 		        _matching->set(n, INVALID);
 		        _level->initAddItem(n);
+		      }
 		      _level->initFinish();
 		      _empty_level = _node_num;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        for (OutArcIt a(_graph, n); a != INVALID; ++a) {
 		          if (_graph.target(a) == n && !_allow_loops) continue;
 		          _matching->set(n, a);
 		          Node v = _graph.target((*_matching)[n]);
 		          _indeg->set(v, (*_indeg)[v] + 1);
 		          break;
+		        }
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_indeg)[n] == 0) {
 		          _level->activate(n);
+		        }
+		      }
+		    }
 		    /// \brief Starts the algorithm and computes a fractional matching
 		    ///
 		    /// The algorithm computes a maximum fractional matching.
 		    ///
 		    /// \param postprocess The algorithm computes first a matching
 		    /// which is a union of a matching with one value edges, cycles
 		    /// with half value edges and even length paths with half value
 		    /// edges. If the parameter is true, then after the push-relabel
 		    /// algorithm it postprocesses the matching to contain only
 		    /// matching edges and half value odd cycles.
 		    void start(bool postprocess = true) {
 		      Node n;
 		      while ((n = _level->highestActive()) != INVALID) {
 		        int level = _level->highestActiveLevel();
 		        int new_level = _level->maxLevel();
 		        for (InArcIt a(_graph, n); a != INVALID; ++a) {
 		          Node u = _graph.source(a);
 		          if (n == u && !_allow_loops) continue;
 		          Node v = _graph.target((*_matching)[u]);
 		          if ((*_level)[v] < level) {
 		            _indeg->set(v, (*_indeg)[v] - 1);
 		            if ((*_indeg)[v] == 0) {
 		              _level->activate(v);
+		            }
 		            _matching->set(u, a);
 		            _indeg->set(n, (*_indeg)[n] + 1);
 		            _level->deactivate(n);
 		            goto no_more_push;
 		          } else if (new_level > (*_level)[v]) {
 		            new_level = (*_level)[v];
+		          }
+		        }
 		        if (new_level + 1 < _level->maxLevel()) {
 		          _level->liftHighestActive(new_level + 1);
 		        } else {
 		          _level->liftHighestActiveToTop();
+		        }
 		        if (_level->emptyLevel(level)) {
 		          _level->liftToTop(level);
+		        }
 		      no_more_push:
+		        ;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] == INVALID) continue;
 		        Node u = _graph.target((*_matching)[n]);
 		        if ((*_indeg)[u] > 1) {
 		          _indeg->set(u, (*_indeg)[u] - 1);
 		          _matching->set(n, INVALID);
+		        }
+		      }
 		      if (postprocess) {
 		        postprocessing();
+		      }
+		    }
 		    /// \brief Starts the algorithm and computes a perfect fractional
 		    /// matching
 		    ///
 		    /// The algorithm computes a perfect fractional matching. If it
 		    /// does not exists, then the algorithm returns false and the
 		    /// matching is undefined and the barrier.
 		    ///
 		    /// \param postprocess The algorithm computes first a matching
 		    /// which is a union of a matching with one value edges, cycles
 		    /// with half value edges and even length paths with half value
 		    /// edges. If the parameter is true, then after the push-relabel
 		    /// algorithm it postprocesses the matching to contain only
 		    /// matching edges and half value odd cycles.
 		    bool startPerfect(bool postprocess = true) {
 		      Node n;
 		      while ((n = _level->highestActive()) != INVALID) {
 		        int level = _level->highestActiveLevel();
 		        int new_level = _level->maxLevel();
 		        for (InArcIt a(_graph, n); a != INVALID; ++a) {
 		          Node u = _graph.source(a);
 		          if (n == u && !_allow_loops) continue;
 		          Node v = _graph.target((*_matching)[u]);
 		          if ((*_level)[v] < level) {
 		            _indeg->set(v, (*_indeg)[v] - 1);
 		            if ((*_indeg)[v] == 0) {
 		              _level->activate(v);
+		            }
 		            _matching->set(u, a);
 		            _indeg->set(n, (*_indeg)[n] + 1);
 		            _level->deactivate(n);
 		            goto no_more_push;
 		          } else if (new_level > (*_level)[v]) {
 		            new_level = (*_level)[v];
+		          }
+		        }
 		        if (new_level + 1 < _level->maxLevel()) {
 		          _level->liftHighestActive(new_level + 1);
 		        } else {
 		          _level->liftHighestActiveToTop();
 		          _empty_level = _level->maxLevel() - 1;
 		          return false;
+		        }
 		        if (_level->emptyLevel(level)) {
 		          _level->liftToTop(level);
 		          _empty_level = level;
 		          return false;
+		        }
 		      no_more_push:
+		        ;
+		      }
 		      if (postprocess) {
 		        postprocessing();
+		      }
 		      return true;
+		    }
 		    /// \brief Runs the algorithm
 		    ///
 		    /// Just a shortcut for the next code:
 		    ///\code
 		    /// init();
 		    /// start();
 		    ///\endcode
 		    void run(bool postprocess = true) {
 		      init();
 		      start(postprocess);
+		    }
 		    /// \brief Runs the algorithm to find a perfect fractional matching
 		    ///
 		    /// Just a shortcut for the next code:
 		    ///\code
 		    /// init();
 		    /// startPerfect();
 		    ///\endcode
 		    bool runPerfect(bool postprocess = true) {
 		      init();
 		      return startPerfect(postprocess);
+		    }
 		    ///@}
 		    /// \name Query Functions
 		    /// The result of the %Matching algorithm can be obtained using these
 		    /// functions.\n
 		    /// Before the use of these functions,
 		    /// either run() or start() must be called.
 		    ///@{
 		    /// \brief Return the number of covered nodes in the matching.
 		    ///
 		    /// This function returns the number of covered nodes in the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matchingSize() const {
 		      int num = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++num;
+		        }
+		      }
 		      return num;
+		    }
 		    /// \brief Returns a const reference to the matching map.
 		    ///
 		    /// Returns a const reference to the node map storing the found
 		    /// fractional matching. This method can be called after
 		    /// running the algorithm.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    const MatchingMap& matchingMap() const {
 		      return *_matching;
+		    }
 		    /// \brief Return \c true if the given edge is in the matching.
 		    ///
 		    /// This function returns \c true if the given edge is in the
 		    /// found matching. The result is scaled by \ref primalScale
 		    /// "primal scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matching(const Edge& edge) const {
 		      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
 		        (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
+		    }
 		    /// \brief Return the fractional matching arc (or edge) incident
 		    /// to the given node.
 		    ///
 		    /// This function returns one of the fractional matching arc (or
 		    /// edge) incident to the given node in the found matching or \c
 		    /// INVALID if the node is not covered by the matching or if the
 		    /// node is on an odd length cycle then it is the successor edge
 		    /// on the cycle.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Returns true if the node is in the barrier
 		    ///
 		    /// The barrier is a subset of the nodes. If the nodes in the
 		    /// barrier have less adjacent nodes than the size of the barrier,
 		    /// then at least as much nodes cannot be covered as the
 		    /// difference of the two subsets.
 		    bool barrier(const Node& node) const {
 		      return (*_level)[node] >= _empty_level;
+		    }
 		    /// @}
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Weighted fractional matching in general graphs
 		  ///
 		  /// This class provides an efficient implementation of fractional
 		  /// matching algorithm. The implementation uses priority queues and
 		  /// provides \f$O(nm\log n)\f$ time complexity.
 		  ///
 		  /// The maximum weighted fractional matching is a relaxation of the
 		  /// maximum weighted matching problem where the odd set constraints
 		  /// are omitted.
 		  /// It can be formulated with the following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
 		  /// \f[x_e \ge 0\quad \forall e\in E\f]
 		  /// \f[\max \sum_{e\in E}x_ew_e\f]
 		  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
 		  /// \f$X\f$. The result must be the union of a matching with one
 		  /// value edges and a set of odd length cycles with half value edges.
 		  ///
 		  /// The algorithm calculates an optimal fractional matching and a
 		  /// proof of the optimality. The solution of the dual problem can be
 		  /// used to check the result of the algorithm. The dual linear
 		  /// problem is the following.
 		  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
 		  /// \f[y_u \ge 0 \quad \forall u \in V\f]
 		  /// \f[\min \sum_{u \in V}y_u \f]
 		  ///
 		  /// The algorithm can be executed with the run() function.
 		  /// After it the matching (the primal solution) and the dual solution
 		  /// can be obtained using the query functions.
 		  ///
 		  /// The primal solution is multiplied by
 		  /// \ref MaxWeightedFractionalMatching::primalScale "2".
 		  /// If the value type is integer, then the dual
 		  /// solution is scaled by
 		  /// \ref MaxWeightedFractionalMatching::dualScale "4".
 		  ///
 		  /// \tparam GR The undirected graph type the algorithm runs on.
 		  /// \tparam WM The type edge weight map. The default type is
 		  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 		#ifdef DOXYGEN
 		  template <typename GR, typename WM>
 		#else
 		  template <typename GR,
 		            typename WM = typename GR::template EdgeMap<int> >
 		#endif
 		  class MaxWeightedFractionalMatching {
 		  public:
 		    /// The graph type of the algorithm
 		    typedef GR Graph;
 		    /// The type of the edge weight map
 		    typedef WM WeightMap;
 		    /// The value type of the edge weights
 		    typedef typename WeightMap::Value Value;
 		    /// The type of the matching map
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		    /// \brief Scaling factor for primal solution
 		    ///
 		    /// Scaling factor for primal solution.
 		    static const int primalScale = 2;
 		    /// \brief Scaling factor for dual solution
 		    ///
 		    /// Scaling factor for dual solution. It is equal to 4 or 1
 		    /// according to the value type.
 		    static const int dualScale =
 		      std::numeric_limits<Value>::is_integer ? 4 : 1;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef typename Graph::template NodeMap<Value> NodePotential;
 		    const Graph& _graph;
 		    const WeightMap& _weight;
 		    MatchingMap* _matching;
 		    NodePotential* _node_potential;
 		    int _node_num;
 		    bool _allow_loops;
 		    enum Status {
 		      EVEN = -1, MATCHED = 0, ODD = 1
 		    };
 		    typedef typename Graph::template NodeMap<Status> StatusMap;
 		    StatusMap* _status;
 		    typedef typename Graph::template NodeMap<Arc> PredMap;
 		    PredMap* _pred;
 		    typedef ExtendFindEnum<IntNodeMap> TreeSet;
 		    IntNodeMap *_tree_set_index;
 		    TreeSet *_tree_set;
 		    IntNodeMap *_delta1_index;
 		    BinHeap<Value, IntNodeMap> *_delta1;
 		    IntNodeMap *_delta2_index;
 		    BinHeap<Value, IntNodeMap> *_delta2;
 		    IntEdgeMap *_delta3_index;
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    Value _delta_sum;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_node_potential) {
 		        _node_potential = new NodePotential(_graph);
+		      }
 		      if (!_status) {
 		        _status = new StatusMap(_graph);
+		      }
 		      if (!_pred) {
 		        _pred = new PredMap(_graph);
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntNodeMap(_graph);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      if (!_delta1) {
 		        _delta1_index = new IntNodeMap(_graph);
 		        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
+		      }
 		      if (!_delta2) {
 		        _delta2_index = new IntNodeMap(_graph);
 		        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
+		      }
 		      if (!_delta3) {
 		        _delta3_index = new IntEdgeMap(_graph);
 		        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+		      }
+		    }
 		    void destroyStructures() {
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_node_potential) {
 		        delete _node_potential;
+		      }
 		      if (_status) {
 		        delete _status;
+		      }
 		      if (_pred) {
 		        delete _pred;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
 		      if (_delta1) {
 		        delete _delta1_index;
 		        delete _delta1;
+		      }
 		      if (_delta2) {
 		        delete _delta2_index;
 		        delete _delta2;
+		      }
 		      if (_delta3) {
 		        delete _delta3_index;
 		        delete _delta3;
+		      }
+		    }
 		    void matchedToEven(Node node, int tree) {
 		      _tree_set->insert(node, tree);
 		      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
 		      _delta1->push(node, (*_node_potential)[node]);
 		      if (_delta2->state(node) == _delta2->IN_HEAP) {
 		        _delta2->erase(node);
+		      }
 		      for (InArcIt a(_graph, node); a != INVALID; ++a) {
 		        Node v = _graph.source(a);
 		        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
 		          dualScale * _weight[a];
 		        if (node == v) {
 		          if (_allow_loops && _graph.direction(a)) {
 		            _delta3->push(a, rw / 2);
+		          }
 		        } else if ((*_status)[v] == EVEN) {
 		          _delta3->push(a, rw / 2);
 		        } else if ((*_status)[v] == MATCHED) {
 		          if (_delta2->state(v) != _delta2->IN_HEAP) {
 		            _pred->set(v, a);
 		            _delta2->push(v, rw);
 		          } else if ((*_delta2)[v] > rw) {
 		            _pred->set(v, a);
 		            _delta2->decrease(v, rw);
+		          }
+		        }
+		      }
+		    }
 		    void matchedToOdd(Node node, int tree) {
 		      _tree_set->insert(node, tree);
 		      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
 		      if (_delta2->state(node) == _delta2->IN_HEAP) {
 		        _delta2->erase(node);
+		      }
+		    }
 		    void evenToMatched(Node node, int tree) {
 		      _delta1->erase(node);
 		      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
 		      Arc min = INVALID;
 		      Value minrw = std::numeric_limits<Value>::max();
 		      for (InArcIt a(_graph, node); a != INVALID; ++a) {
 		        Node v = _graph.source(a);
 		        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
 		          dualScale * _weight[a];
 		        if (node == v) {
 		          if (_allow_loops && _graph.direction(a)) {
 		            _delta3->erase(a);
+		          }
 		        } else if ((*_status)[v] == EVEN) {
 		          _delta3->erase(a);
 		          if (minrw > rw) {
 		            min = _graph.oppositeArc(a);
 		            minrw = rw;
+		          }
 		        } else if ((*_status)[v]  == MATCHED) {
 		          if ((*_pred)[v] == a) {
 		            Arc mina = INVALID;
 		            Value minrwa = std::numeric_limits<Value>::max();
 		            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
 		              Node va = _graph.target(aa);
 		              if ((*_status)[va] != EVEN ||
 		                  _tree_set->find(va) == tree) continue;
 		              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
 		                dualScale * _weight[aa];
 		              if (minrwa > rwa) {
 		                minrwa = rwa;
 		                mina = aa;
+		              }
+		            }
 		            if (mina != INVALID) {
 		              _pred->set(v, mina);
 		              _delta2->increase(v, minrwa);
 		            } else {
 		              _pred->set(v, INVALID);
 		              _delta2->erase(v);
+		            }
+		          }
+		        }
+		      }
 		      if (min != INVALID) {
 		        _pred->set(node, min);
 		        _delta2->push(node, minrw);
 		      } else {
 		        _pred->set(node, INVALID);
+		      }
+		    }
 		    void oddToMatched(Node node) {
 		      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
 		      Arc min = INVALID;
 		      Value minrw = std::numeric_limits<Value>::max();
 		      for (InArcIt a(_graph, node); a != INVALID; ++a) {
 		        Node v = _graph.source(a);
 		        if ((*_status)[v] != EVEN) continue;
 		        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
 		          dualScale * _weight[a];
 		        if (minrw > rw) {
 		          min = _graph.oppositeArc(a);
 		          minrw = rw;
+		        }
+		      }
 		      if (min != INVALID) {
 		        _pred->set(node, min);
 		        _delta2->push(node, minrw);
 		      } else {
 		        _pred->set(node, INVALID);
+		      }
+		    }
 		    void alternatePath(Node even, int tree) {
 		      Node odd;
 		      _status->set(even, MATCHED);
 		      evenToMatched(even, tree);
 		      Arc prev = (*_matching)[even];
 		      while (prev != INVALID) {
 		        odd = _graph.target(prev);
 		        even = _graph.target((*_pred)[odd]);
 		        _matching->set(odd, (*_pred)[odd]);
 		        _status->set(odd, MATCHED);
 		        oddToMatched(odd);
 		        prev = (*_matching)[even];
 		        _status->set(even, MATCHED);
 		        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
 		        evenToMatched(even, tree);
+		      }
+		    }
 		    void destroyTree(int tree) {
 		      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
 		        if ((*_status)[n] == EVEN) {
 		          _status->set(n, MATCHED);
 		          evenToMatched(n, tree);
 		        } else if ((*_status)[n] == ODD) {
 		          _status->set(n, MATCHED);
 		          oddToMatched(n);
+		        }
+		      }
 		      _tree_set->eraseClass(tree);
+		    }
 		    void unmatchNode(const Node& node) {
 		      int tree = _tree_set->find(node);
 		      alternatePath(node, tree);
 		      destroyTree(tree);
 		      _matching->set(node, INVALID);
+		    }
 		    void augmentOnEdge(const Edge& edge) {
 		      Node left = _graph.u(edge);
 		      int left_tree = _tree_set->find(left);
 		      alternatePath(left, left_tree);
 		      destroyTree(left_tree);
 		      _matching->set(left, _graph.direct(edge, true));
 		      Node right = _graph.v(edge);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      _matching->set(right, _graph.direct(edge, false));
+		    }
 		    void augmentOnArc(const Arc& arc) {
 		      Node left = _graph.source(arc);
 		      _status->set(left, MATCHED);
 		      _matching->set(left, arc);
 		      _pred->set(left, arc);
 		      Node right = _graph.target(arc);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      _matching->set(right, _graph.oppositeArc(arc));
+		    }
 		    void extendOnArc(const Arc& arc) {
 		      Node base = _graph.target(arc);
 		      int tree = _tree_set->find(base);
 		      Node odd = _graph.source(arc);
 		      _tree_set->insert(odd, tree);
 		      _status->set(odd, ODD);
 		      matchedToOdd(odd, tree);
 		      _pred->set(odd, arc);
 		      Node even = _graph.target((*_matching)[odd]);
 		      _tree_set->insert(even, tree);
 		      _status->set(even, EVEN);
 		      matchedToEven(even, tree);
+		    }
 		    void cycleOnEdge(const Edge& edge, int tree) {
 		      Node nca = INVALID;
 		      std::vector<Node> left_path, right_path;
+		      {
 		        std::set<Node> left_set, right_set;
 		        Node left = _graph.u(edge);
 		        left_path.push_back(left);
 		        left_set.insert(left);
 		        Node right = _graph.v(edge);
 		        right_path.push_back(right);
 		        right_set.insert(right);
 		        while (true) {
 		          if (left_set.find(right) != left_set.end()) {
 		            nca = right;
 		            break;
+		          }
 		          if ((*_matching)[left] == INVALID) break;
 		          left = _graph.target((*_matching)[left]);
 		          left_path.push_back(left);
 		          left = _graph.target((*_pred)[left]);
 		          left_path.push_back(left);
 		          left_set.insert(left);
 		          if (right_set.find(left) != right_set.end()) {
 		            nca = left;
 		            break;
+		          }
 		          if ((*_matching)[right] == INVALID) break;
 		          right = _graph.target((*_matching)[right]);
 		          right_path.push_back(right);
 		          right = _graph.target((*_pred)[right]);
 		          right_path.push_back(right);
 		          right_set.insert(right);
+		        }
 		        if (nca == INVALID) {
 		          if ((*_matching)[left] == INVALID) {
 		            nca = right;
 		            while (left_set.find(nca) == left_set.end()) {
 		              nca = _graph.target((*_matching)[nca]);
 		              right_path.push_back(nca);
 		              nca = _graph.target((*_pred)[nca]);
 		              right_path.push_back(nca);
+		            }
 		          } else {
 		            nca = left;
 		            while (right_set.find(nca) == right_set.end()) {
 		              nca = _graph.target((*_matching)[nca]);
 		              left_path.push_back(nca);
 		              nca = _graph.target((*_pred)[nca]);
 		              left_path.push_back(nca);
+		            }
+		          }
+		        }
+		      }
 		      alternatePath(nca, tree);
 		      Arc prev;
 		      prev = _graph.direct(edge, true);
 		      for (int i = 0; left_path[i] != nca; i += 2) {
 		        _matching->set(left_path[i], prev);
 		        _status->set(left_path[i], MATCHED);
 		        evenToMatched(left_path[i], tree);
 		        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
 		        _status->set(left_path[i + 1], MATCHED);
 		        oddToMatched(left_path[i + 1]);
+		      }
 		      _matching->set(nca, prev);
 		      for (int i = 0; right_path[i] != nca; i += 2) {
 		        _status->set(right_path[i], MATCHED);
 		        evenToMatched(right_path[i], tree);
 		        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
 		        _status->set(right_path[i + 1], MATCHED);
 		        oddToMatched(right_path[i + 1]);
+		      }
 		      destroyTree(tree);
+		    }
 		    void extractCycle(const Arc &arc) {
 		      Node left = _graph.source(arc);
 		      Node odd = _graph.target((*_matching)[left]);
 		      Arc prev;
 		      while (odd != left) {
 		        Node even = _graph.target((*_matching)[odd]);
 		        prev = (*_matching)[odd];
 		        odd = _graph.target((*_matching)[even]);
 		        _matching->set(even, _graph.oppositeArc(prev));
+		      }
 		      _matching->set(left, arc);
 		      Node right = _graph.target(arc);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      _matching->set(right, _graph.oppositeArc(arc));
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
 		                                  bool allow_loops = true)
 		      : _graph(graph), _weight(weight), _matching(0),
 		      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
 		      _status(0),  _pred(0),
 		      _tree_set_index(0), _tree_set(0),
 		      _delta1_index(0), _delta1(0),
 		      _delta2_index(0), _delta2(0),
 		      _delta3_index(0), _delta3(0),
 		      _delta_sum() {}
 		    ~MaxWeightedFractionalMatching() {
 		      destroyStructures();
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \ref run() member function.
 		    ///@{
 		    /// \brief Initialize the algorithm
 		    ///
 		    /// This function initializes the algorithm.
 		    void init() {
 		      createStructures();
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_delta1_index)[n] = _delta1->PRE_HEAP;
 		        (*_delta2_index)[n] = _delta2->PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = 0;
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n && !_allow_loops) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
 		            max = (dualScale * _weight[e]) / 2;
+		          }
+		        }
 		        _node_potential->set(n, max);
 		        _delta1->push(n, max);
 		        _tree_set->insert(n);
 		        _matching->set(n, INVALID);
 		        _status->set(n, EVEN);
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        Node left = _graph.u(e);
 		        Node right = _graph.v(e);
 		        if (left == right && !_allow_loops) continue;
 		        _delta3->push(e, ((*_node_potential)[left] +
 		                          (*_node_potential)[right] -
 		                          dualScale * _weight[e]) / 2);
+		      }
+		    }
 		    /// \brief Start the algorithm
 		    ///
 		    /// This function starts the algorithm.
 		    ///
 		    /// \pre \ref init() must be called before using this function.
 		    void start() {
 		      enum OpType {
 		        D1, D2, D3
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		        Value d1 = !_delta1->empty() ?
 		          _delta1->prio() : std::numeric_limits<Value>::max();
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
 		        _delta_sum = d3; OpType ot = D3;
 		        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
 		        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
 		        switch (ot) {
 		        case D1:
+		          {
 		            Node n = _delta1->top();
 		            unmatchNode(n);
 		            --unmatched;
+		          }
 		          break;
 		        case D2:
+		          {
 		            Node n = _delta2->top();
 		            Arc a = (*_pred)[n];
 		            if ((*_matching)[n] == INVALID) {
 		              augmentOnArc(a);
 		              --unmatched;
 		            } else {
 		              Node v = _graph.target((*_matching)[n]);
 		              if ((*_matching)[n] !=
 		                  _graph.oppositeArc((*_matching)[v])) {
 		                extractCycle(a);
 		                --unmatched;
 		              } else {
 		                extendOnArc(a);
+		              }
+		            }
 		          } break;
 		        case D3:
+		          {
 		            Edge e = _delta3->top();
 		            Node left = _graph.u(e);
 		            Node right = _graph.v(e);
 		            int left_tree = _tree_set->find(left);
 		            int right_tree = _tree_set->find(right);
 		            if (left_tree == right_tree) {
 		              cycleOnEdge(e, left_tree);
 		              --unmatched;
 		            } else {
 		              augmentOnEdge(e);
 		              unmatched -= 2;
+		            }
 		          } break;
+		        }
+		      }
+		    }
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
 		    ///
 		    /// \note mwfm.run() is just a shortcut of the following code.
 		    /// \code
 		    ///   mwfm.init();
 		    ///   mwfm.start();
 		    /// \endcode
 		    void run() {
 		      init();
 		      start();
+		    }
 		    /// @}
 		    /// \name Primal Solution
 		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// matching.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the weight of the matching.
 		    ///
 		    /// This function returns the weight of the found matching. This
 		    /// value is scaled by \ref primalScale "primal scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value matchingWeight() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          sum += _weight[(*_matching)[n]];
+		        }
+		      }
 		      return sum * primalScale / 2;
+		    }
 		    /// \brief Return the number of covered nodes in the matching.
 		    ///
 		    /// This function returns the number of covered nodes in the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matchingSize() const {
 		      int num = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++num;
+		        }
+		      }
 		      return num;
+		    }
 		    /// \brief Return \c true if the given edge is in the matching.
 		    ///
 		    /// This function returns \c true if the given edge is in the
 		    /// found matching. The result is scaled by \ref primalScale
 		    /// "primal scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matching(const Edge& edge) const {
 		      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
 		        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
+		    }
 		    /// \brief Return the fractional matching arc (or edge) incident
 		    /// to the given node.
 		    ///
 		    /// This function returns one of the fractional matching arc (or
 		    /// edge) incident to the given node in the found matching or \c
 		    /// INVALID if the node is not covered by the matching or if the
 		    /// node is on an odd length cycle then it is the successor edge
 		    /// on the cycle.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Return a const reference to the matching map.
 		    ///
 		    /// This function returns a const reference to a node map that stores
 		    /// the matching arc (or edge) incident to each node.
 		    const MatchingMap& matchingMap() const {
 		      return *_matching;
+		    }
 		    /// @}
 		    /// \name Dual Solution
 		    /// Functions to get the dual solution.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the value of the dual solution.
 		    ///
 		    /// This function returns the value of the dual solution.
 		    /// It should be equal to the primal value scaled by \ref dualScale
 		    /// "dual scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value dualValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        sum += nodeValue(n);
+		      }
 		      return sum;
+		    }
 		    /// \brief Return the dual value (potential) of the given node.
 		    ///
 		    /// This function returns the dual value (potential) of the given node.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value nodeValue(const Node& n) const {
 		      return (*_node_potential)[n];
+		    }
 		    /// @}
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Weighted fractional perfect matching in general graphs
 		  ///
 		  /// This class provides an efficient implementation of fractional
 		  /// matching algorithm. The implementation uses priority queues and
 		  /// provides \f$O(nm\log n)\f$ time complexity.
 		  ///
 		  /// The maximum weighted fractional perfect matching is a relaxation
 		  /// of the maximum weighted perfect matching problem where the odd
 		  /// set constraints are omitted.
 		  /// It can be formulated with the following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
 		  /// \f[x_e \ge 0\quad \forall e\in E\f]
 		  /// \f[\max \sum_{e\in E}x_ew_e\f]
 		  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
 		  /// \f$X\f$. The result must be the union of a matching with one
 		  /// value edges and a set of odd length cycles with half value edges.
 		  ///
 		  /// The algorithm calculates an optimal fractional matching and a
 		  /// proof of the optimality. The solution of the dual problem can be
 		  /// used to check the result of the algorithm. The dual linear
 		  /// problem is the following.
 		  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
 		  /// \f[\min \sum_{u \in V}y_u \f]
 		  ///
 		  /// The algorithm can be executed with the run() function.
 		  /// After it the matching (the primal solution) and the dual solution
 		  /// can be obtained using the query functions.
 		  ///
 		  /// The primal solution is multiplied by
 		  /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
 		  /// If the value type is integer, then the dual
 		  /// solution is scaled by
 		  /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
 		  ///
 		  /// \tparam GR The undirected graph type the algorithm runs on.
 		  /// \tparam WM The type edge weight map. The default type is
 		  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 		#ifdef DOXYGEN
 		  template <typename GR, typename WM>
 		#else
 		  template <typename GR,
 		            typename WM = typename GR::template EdgeMap<int> >
 		#endif
 		  class MaxWeightedPerfectFractionalMatching {
 		  public:
 		    /// The graph type of the algorithm
 		    typedef GR Graph;
 		    /// The type of the edge weight map
 		    typedef WM WeightMap;
 		    /// The value type of the edge weights
 		    typedef typename WeightMap::Value Value;
 		    /// The type of the matching map
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		    /// \brief Scaling factor for primal solution
 		    ///
 		    /// Scaling factor for primal solution.
 		    static const int primalScale = 2;
 		    /// \brief Scaling factor for dual solution
 		    ///
 		    /// Scaling factor for dual solution. It is equal to 4 or 1
 		    /// according to the value type.
 		    static const int dualScale =
 		      std::numeric_limits<Value>::is_integer ? 4 : 1;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef typename Graph::template NodeMap<Value> NodePotential;
 		    const Graph& _graph;
 		    const WeightMap& _weight;
 		    MatchingMap* _matching;
 		    NodePotential* _node_potential;
 		    int _node_num;
 		    bool _allow_loops;
 		    enum Status {
 		      EVEN = -1, MATCHED = 0, ODD = 1
 		    };
 		    typedef typename Graph::template NodeMap<Status> StatusMap;
 		    StatusMap* _status;
 		    typedef typename Graph::template NodeMap<Arc> PredMap;
 		    PredMap* _pred;
 		    typedef ExtendFindEnum<IntNodeMap> TreeSet;
 		    IntNodeMap *_tree_set_index;
 		    TreeSet *_tree_set;
 		    IntNodeMap *_delta2_index;
 		    BinHeap<Value, IntNodeMap> *_delta2;
 		    IntEdgeMap *_delta3_index;
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    Value _delta_sum;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_node_potential) {
 		        _node_potential = new NodePotential(_graph);
+		      }
 		      if (!_status) {
 		        _status = new StatusMap(_graph);
+		      }
 		      if (!_pred) {
 		        _pred = new PredMap(_graph);
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntNodeMap(_graph);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      if (!_delta2) {
 		        _delta2_index = new IntNodeMap(_graph);
 		        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
+		      }
 		      if (!_delta3) {
 		        _delta3_index = new IntEdgeMap(_graph);
 		        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+		      }
+		    }
 		    void destroyStructures() {
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_node_potential) {
 		        delete _node_potential;
+		      }
 		      if (_status) {
 		        delete _status;
+		      }
 		      if (_pred) {
 		        delete _pred;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
 		      if (_delta2) {
 		        delete _delta2_index;
 		        delete _delta2;
+		      }
 		      if (_delta3) {
 		        delete _delta3_index;
 		        delete _delta3;
+		      }
+		    }
 		    void matchedToEven(Node node, int tree) {
 		      _tree_set->insert(node, tree);
 		      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
 		      if (_delta2->state(node) == _delta2->IN_HEAP) {
 		        _delta2->erase(node);
+		      }
 		      for (InArcIt a(_graph, node); a != INVALID; ++a) {
 		        Node v = _graph.source(a);
 		        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
 		          dualScale * _weight[a];
 		        if (node == v) {
 		          if (_allow_loops && _graph.direction(a)) {
 		            _delta3->push(a, rw / 2);
+		          }
 		        } else if ((*_status)[v] == EVEN) {
 		          _delta3->push(a, rw / 2);
 		        } else if ((*_status)[v] == MATCHED) {
 		          if (_delta2->state(v) != _delta2->IN_HEAP) {
 		            _pred->set(v, a);
 		            _delta2->push(v, rw);
 		          } else if ((*_delta2)[v] > rw) {
 		            _pred->set(v, a);
 		            _delta2->decrease(v, rw);
+		          }
+		        }
+		      }
+		    }
 		    void matchedToOdd(Node node, int tree) {
 		      _tree_set->insert(node, tree);
 		      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
 		      if (_delta2->state(node) == _delta2->IN_HEAP) {
 		        _delta2->erase(node);
+		      }
+		    }
 		    void evenToMatched(Node node, int tree) {
 		      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
 		      Arc min = INVALID;
 		      Value minrw = std::numeric_limits<Value>::max();
 		      for (InArcIt a(_graph, node); a != INVALID; ++a) {
 		        Node v = _graph.source(a);
 		        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
 		          dualScale * _weight[a];
 		        if (node == v) {
 		          if (_allow_loops && _graph.direction(a)) {
 		            _delta3->erase(a);
+		          }
 		        } else if ((*_status)[v] == EVEN) {
 		          _delta3->erase(a);
 		          if (minrw > rw) {
 		            min = _graph.oppositeArc(a);
 		            minrw = rw;
+		          }
 		        } else if ((*_status)[v]  == MATCHED) {
 		          if ((*_pred)[v] == a) {
 		            Arc mina = INVALID;
 		            Value minrwa = std::numeric_limits<Value>::max();
 		            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
 		              Node va = _graph.target(aa);
 		              if ((*_status)[va] != EVEN ||
 		                  _tree_set->find(va) == tree) continue;
 		              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
 		                dualScale * _weight[aa];
 		              if (minrwa > rwa) {
 		                minrwa = rwa;
 		                mina = aa;
+		              }
+		            }
 		            if (mina != INVALID) {
 		              _pred->set(v, mina);
 		              _delta2->increase(v, minrwa);
 		            } else {
 		              _pred->set(v, INVALID);
 		              _delta2->erase(v);
+		            }
+		          }
+		        }
+		      }
 		      if (min != INVALID) {
 		        _pred->set(node, min);
 		        _delta2->push(node, minrw);
 		      } else {
 		        _pred->set(node, INVALID);
+		      }
+		    }
 		    void oddToMatched(Node node) {
 		      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
 		      Arc min = INVALID;
 		      Value minrw = std::numeric_limits<Value>::max();
 		      for (InArcIt a(_graph, node); a != INVALID; ++a) {
 		        Node v = _graph.source(a);
 		        if ((*_status)[v] != EVEN) continue;
 		        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
 		          dualScale * _weight[a];
 		        if (minrw > rw) {
 		          min = _graph.oppositeArc(a);
 		          minrw = rw;
+		        }
+		      }
 		      if (min != INVALID) {
 		        _pred->set(node, min);
 		        _delta2->push(node, minrw);
 		      } else {
 		        _pred->set(node, INVALID);
+		      }
+		    }
 		    void alternatePath(Node even, int tree) {
 		      Node odd;
 		      _status->set(even, MATCHED);
 		      evenToMatched(even, tree);
 		      Arc prev = (*_matching)[even];
 		      while (prev != INVALID) {
 		        odd = _graph.target(prev);
 		        even = _graph.target((*_pred)[odd]);
 		        _matching->set(odd, (*_pred)[odd]);
 		        _status->set(odd, MATCHED);
 		        oddToMatched(odd);
 		        prev = (*_matching)[even];
 		        _status->set(even, MATCHED);
 		        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
 		        evenToMatched(even, tree);
+		      }
+		    }
 		    void destroyTree(int tree) {
 		      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
 		        if ((*_status)[n] == EVEN) {
 		          _status->set(n, MATCHED);
 		          evenToMatched(n, tree);
 		        } else if ((*_status)[n] == ODD) {
 		          _status->set(n, MATCHED);
 		          oddToMatched(n);
+		        }
+		      }
 		      _tree_set->eraseClass(tree);
+		    }
 		    void augmentOnEdge(const Edge& edge) {
 		      Node left = _graph.u(edge);
 		      int left_tree = _tree_set->find(left);
 		      alternatePath(left, left_tree);
 		      destroyTree(left_tree);
 		      _matching->set(left, _graph.direct(edge, true));
 		      Node right = _graph.v(edge);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      _matching->set(right, _graph.direct(edge, false));
+		    }
 		    void augmentOnArc(const Arc& arc) {
 		      Node left = _graph.source(arc);
 		      _status->set(left, MATCHED);
 		      _matching->set(left, arc);
 		      _pred->set(left, arc);
 		      Node right = _graph.target(arc);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      _matching->set(right, _graph.oppositeArc(arc));
+		    }
 		    void extendOnArc(const Arc& arc) {
 		      Node base = _graph.target(arc);
 		      int tree = _tree_set->find(base);
 		      Node odd = _graph.source(arc);
 		      _tree_set->insert(odd, tree);
 		      _status->set(odd, ODD);
 		      matchedToOdd(odd, tree);
 		      _pred->set(odd, arc);
 		      Node even = _graph.target((*_matching)[odd]);
 		      _tree_set->insert(even, tree);
 		      _status->set(even, EVEN);
 		      matchedToEven(even, tree);
+		    }
 		    void cycleOnEdge(const Edge& edge, int tree) {
 		      Node nca = INVALID;
 		      std::vector<Node> left_path, right_path;
+		      {
 		        std::set<Node> left_set, right_set;
 		        Node left = _graph.u(edge);
 		        left_path.push_back(left);
 		        left_set.insert(left);
 		        Node right = _graph.v(edge);
 		        right_path.push_back(right);
 		        right_set.insert(right);
 		        while (true) {
 		          if (left_set.find(right) != left_set.end()) {
 		            nca = right;
 		            break;
+		          }
 		          if ((*_matching)[left] == INVALID) break;
 		          left = _graph.target((*_matching)[left]);
 		          left_path.push_back(left);
 		          left = _graph.target((*_pred)[left]);
 		          left_path.push_back(left);
 		          left_set.insert(left);
 		          if (right_set.find(left) != right_set.end()) {
 		            nca = left;
 		            break;
+		          }
 		          if ((*_matching)[right] == INVALID) break;
 		          right = _graph.target((*_matching)[right]);
 		          right_path.push_back(right);
 		          right = _graph.target((*_pred)[right]);
 		          right_path.push_back(right);
 		          right_set.insert(right);
+		        }
 		        if (nca == INVALID) {
 		          if ((*_matching)[left] == INVALID) {
 		            nca = right;
 		            while (left_set.find(nca) == left_set.end()) {
 		              nca = _graph.target((*_matching)[nca]);
 		              right_path.push_back(nca);
 		              nca = _graph.target((*_pred)[nca]);
 		              right_path.push_back(nca);
+		            }
 		          } else {
 		            nca = left;
 		            while (right_set.find(nca) == right_set.end()) {
 		              nca = _graph.target((*_matching)[nca]);
 		              left_path.push_back(nca);
 		              nca = _graph.target((*_pred)[nca]);
 		              left_path.push_back(nca);
+		            }
+		          }
+		        }
+		      }
 		      alternatePath(nca, tree);
 		      Arc prev;
 		      prev = _graph.direct(edge, true);
 		      for (int i = 0; left_path[i] != nca; i += 2) {
 		        _matching->set(left_path[i], prev);
 		        _status->set(left_path[i], MATCHED);
 		        evenToMatched(left_path[i], tree);
 		        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
 		        _status->set(left_path[i + 1], MATCHED);
 		        oddToMatched(left_path[i + 1]);
+		      }
 		      _matching->set(nca, prev);
 		      for (int i = 0; right_path[i] != nca; i += 2) {
 		        _status->set(right_path[i], MATCHED);
 		        evenToMatched(right_path[i], tree);
 		        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
 		        _status->set(right_path[i + 1], MATCHED);
 		        oddToMatched(right_path[i + 1]);
+		      }
 		      destroyTree(tree);
+		    }
 		    void extractCycle(const Arc &arc) {
 		      Node left = _graph.source(arc);
 		      Node odd = _graph.target((*_matching)[left]);
 		      Arc prev;
 		      while (odd != left) {
 		        Node even = _graph.target((*_matching)[odd]);
 		        prev = (*_matching)[odd];
 		        odd = _graph.target((*_matching)[even]);
 		        _matching->set(even, _graph.oppositeArc(prev));
+		      }
 		      _matching->set(left, arc);
 		      Node right = _graph.target(arc);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      _matching->set(right, _graph.oppositeArc(arc));
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxWeightedPerfectFractionalMatching(const Graph& graph,
 		                                         const WeightMap& weight,
 		                                         bool allow_loops = true)
 		      : _graph(graph), _weight(weight), _matching(0),
 		      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
 		      _status(0),  _pred(0),
 		      _tree_set_index(0), _tree_set(0),
 		      _delta2_index(0), _delta2(0),
 		      _delta3_index(0), _delta3(0),
 		      _delta_sum() {}
 		    ~MaxWeightedPerfectFractionalMatching() {
 		      destroyStructures();
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \ref run() member function.
 		    ///@{
 		    /// \brief Initialize the algorithm
 		    ///
 		    /// This function initializes the algorithm.
 		    void init() {
 		      createStructures();
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_delta2_index)[n] = _delta2->PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = - std::numeric_limits<Value>::max();
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n && !_allow_loops) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
 		            max = (dualScale * _weight[e]) / 2;
+		          }
+		        }
 		        _node_potential->set(n, max);
 		        _tree_set->insert(n);
 		        _matching->set(n, INVALID);
 		        _status->set(n, EVEN);
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        Node left = _graph.u(e);
 		        Node right = _graph.v(e);
 		        if (left == right && !_allow_loops) continue;
 		        _delta3->push(e, ((*_node_potential)[left] +
 		                          (*_node_potential)[right] -
 		                          dualScale * _weight[e]) / 2);
+		      }
+		    }
 		    /// \brief Start the algorithm
 		    ///
 		    /// This function starts the algorithm.
 		    ///
 		    /// \pre \ref init() must be called before using this function.
 		    bool start() {
 		      enum OpType {
 		        D2, D3
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
 		        _delta_sum = d3; OpType ot = D3;
 		        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
 		        if (_delta_sum == std::numeric_limits<Value>::max()) {
 		          return false;
+		        }
 		        switch (ot) {
 		        case D2:
+		          {
 		            Node n = _delta2->top();
 		            Arc a = (*_pred)[n];
 		            if ((*_matching)[n] == INVALID) {
 		              augmentOnArc(a);
 		              --unmatched;
 		            } else {
 		              Node v = _graph.target((*_matching)[n]);
 		              if ((*_matching)[n] !=
 		                  _graph.oppositeArc((*_matching)[v])) {
 		                extractCycle(a);
 		                --unmatched;
 		              } else {
 		                extendOnArc(a);
+		              }
+		            }
 		          } break;
 		        case D3:
+		          {
 		            Edge e = _delta3->top();
 		            Node left = _graph.u(e);
 		            Node right = _graph.v(e);
 		            int left_tree = _tree_set->find(left);
 		            int right_tree = _tree_set->find(right);
 		            if (left_tree == right_tree) {
 		              cycleOnEdge(e, left_tree);
 		              --unmatched;
 		            } else {
 		              augmentOnEdge(e);
 		              unmatched -= 2;
+		            }
 		          } break;
+		        }
+		      }
 		      return true;
+		    }
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This method runs the \c %MaxWeightedPerfectFractionalMatching
 		    /// algorithm.
 		    ///
 		    /// \note mwfm.run() is just a shortcut of the following code.
 		    /// \code
 		    ///   mwpfm.init();
 		    ///   mwpfm.start();
 		    /// \endcode
 		    bool run() {
 		      init();
 		      return start();
+		    }
 		    /// @}
 		    /// \name Primal Solution
 		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// matching.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the weight of the matching.
 		    ///
 		    /// This function returns the weight of the found matching. This
 		    /// value is scaled by \ref primalScale "primal scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value matchingWeight() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          sum += _weight[(*_matching)[n]];
+		        }
+		      }
 		      return sum * primalScale / 2;
+		    }
 		    /// \brief Return the number of covered nodes in the matching.
 		    ///
 		    /// This function returns the number of covered nodes in the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matchingSize() const {
 		      int num = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++num;
+		        }
+		      }
 		      return num;
+		    }
 		    /// \brief Return \c true if the given edge is in the matching.
 		    ///
 		    /// This function returns \c true if the given edge is in the
 		    /// found matching. The result is scaled by \ref primalScale
 		    /// "primal scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matching(const Edge& edge) const {
 		      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
 		        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
+		    }
 		    /// \brief Return the fractional matching arc (or edge) incident
 		    /// to the given node.
 		    ///
 		    /// This function returns one of the fractional matching arc (or
 		    /// edge) incident to the given node in the found matching or \c
 		    /// INVALID if the node is not covered by the matching or if the
 		    /// node is on an odd length cycle then it is the successor edge
 		    /// on the cycle.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Return a const reference to the matching map.
 		    ///
 		    /// This function returns a const reference to a node map that stores
 		    /// the matching arc (or edge) incident to each node.
 		    const MatchingMap& matchingMap() const {
 		      return *_matching;
+		    }
 		    /// @}
 		    /// \name Dual Solution
 		    /// Functions to get the dual solution.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the value of the dual solution.
 		    ///
 		    /// This function returns the value of the dual solution.
 		    /// It should be equal to the primal value scaled by \ref dualScale
 		    /// "dual scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value dualValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        sum += nodeValue(n);
+		      }
 		      return sum;
+		    }
 		    /// \brief Return the dual value (potential) of the given node.
 		    ///
 		    /// This function returns the dual value (potential) of the given node.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value nodeValue(const Node& n) const {
 		      return (*_node_potential)[n];
+		    }
 		    /// @}
 		  };
 		} //END OF NAMESPACE LEMON
 		#endif //LEMON_FRACTIONAL_MATCHING_H

test/fractional_matching_test.cc

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#include <iostream>
 		#include <sstream>
 		#include <vector>
 		#include <queue>
 		#include <cstdlib>
 		#include <lemon/fractional_matching.h>
 		#include <lemon/smart_graph.h>
 		#include <lemon/concepts/graph.h>
 		#include <lemon/concepts/maps.h>
 		#include <lemon/lgf_reader.h>
 		#include <lemon/math.h>
 		#include "test_tools.h"
 		using namespace std;
 		using namespace lemon;
 		GRAPH_TYPEDEFS(SmartGraph);
 		const int lgfn = 4;
 		const std::string lgf[lgfn] = {
 		  "@nodes\n"
 		  "label\n"
 		  "0\n"
 		  "1\n"
 		  "2\n"
 		  "3\n"
 		  "4\n"
 		  "5\n"
 		  "6\n"
 		  "7\n"
 		  "@edges\n"
 		  "     label  weight\n"
 		  "7 4  0      984\n"
 		  "0 7  1      73\n"
 		  "7 1  2      204\n"
 		  "2 3  3      583\n"
 		  "2 7  4      565\n"
 		  "2 1  5      582\n"
 		  "0 4  6      551\n"
 		  "2 5  7      385\n"
 		  "1 5  8      561\n"
 		  "5 3  9      484\n"
 		  "7 5  10     904\n"
 		  "3 6  11     47\n"
 		  "7 6  12     888\n"
 		  "3 0  13     747\n"
 		  "6 1  14     310\n",
 		  "@nodes\n"
 		  "label\n"
 		  "0\n"
 		  "1\n"
 		  "2\n"
 		  "3\n"
 		  "4\n"
 		  "5\n"
 		  "6\n"
 		  "7\n"
 		  "@edges\n"
 		  "     label  weight\n"
 		  "2 5  0      710\n"
 		  "0 5  1      241\n"
 		  "2 4  2      856\n"
 		  "2 6  3      762\n"
 		  "4 1  4      747\n"
 		  "6 1  5      962\n"
 		  "4 7  6      723\n"
 		  "1 7  7      661\n"
 		  "2 3  8      376\n"
 		  "1 0  9      416\n"
 		  "6 7  10     391\n",
 		  "@nodes\n"
 		  "label\n"
 		  "0\n"
 		  "1\n"
 		  "2\n"
 		  "3\n"
 		  "4\n"
 		  "5\n"
 		  "6\n"
 		  "7\n"
 		  "@edges\n"
 		  "     label  weight\n"
 		  "6 2  0      553\n"
 		  "0 7  1      653\n"
 		  "6 3  2      22\n"
 		  "4 7  3      846\n"
 		  "7 2  4      981\n"
 		  "7 6  5      250\n"
 		  "5 2  6      539\n",
 		  "@nodes\n"
 		  "label\n"
 		  "0\n"
 		  "@edges\n"
 		  "     label  weight\n"
 		  "0 0  0      100\n"
 		};
 		void checkMaxFractionalMatchingCompile()
+		{
 		  typedef concepts::Graph Graph;
 		  typedef Graph::Node Node;
 		  typedef Graph::Edge Edge;
 		  Graph g;
 		  Node n;
 		  Edge e;
 		  MaxFractionalMatching<Graph> mat_test(g);
 		  const MaxFractionalMatching<Graph>&
 		    const_mat_test = mat_test;
 		  mat_test.init();
 		  mat_test.start();
 		  mat_test.start(true);
 		  mat_test.startPerfect();
 		  mat_test.startPerfect(true);
 		  mat_test.run();
 		  mat_test.run(true);
 		  mat_test.runPerfect();
 		  mat_test.runPerfect(true);
 		  const_mat_test.matchingSize();
 		  const_mat_test.matching(e);
 		  const_mat_test.matching(n);
 		  const MaxFractionalMatching<Graph>::MatchingMap& mmap =
 		    const_mat_test.matchingMap();
 		  e = mmap[n];
 		  const_mat_test.barrier(n);
+		}
 		void checkMaxWeightedFractionalMatchingCompile()
+		{
 		  typedef concepts::Graph Graph;
 		  typedef Graph::Node Node;
 		  typedef Graph::Edge Edge;
 		  typedef Graph::EdgeMap<int> WeightMap;
 		  Graph g;
 		  Node n;
 		  Edge e;
 		  WeightMap w(g);
 		  MaxWeightedFractionalMatching<Graph> mat_test(g, w);
 		  const MaxWeightedFractionalMatching<Graph>&
 		    const_mat_test = mat_test;
 		  mat_test.init();
 		  mat_test.start();
 		  mat_test.run();
 		  const_mat_test.matchingWeight();
 		  const_mat_test.matchingSize();
 		  const_mat_test.matching(e);
 		  const_mat_test.matching(n);
 		  const MaxWeightedFractionalMatching<Graph>::MatchingMap& mmap =
 		    const_mat_test.matchingMap();
 		  e = mmap[n];
 		  const_mat_test.dualValue();
 		  const_mat_test.nodeValue(n);
+		}
 		void checkMaxWeightedPerfectFractionalMatchingCompile()
+		{
 		  typedef concepts::Graph Graph;
 		  typedef Graph::Node Node;
 		  typedef Graph::Edge Edge;
 		  typedef Graph::EdgeMap<int> WeightMap;
 		  Graph g;
 		  Node n;
 		  Edge e;
 		  WeightMap w(g);
 		  MaxWeightedPerfectFractionalMatching<Graph> mat_test(g, w);
 		  const MaxWeightedPerfectFractionalMatching<Graph>&
 		    const_mat_test = mat_test;
 		  mat_test.init();
 		  mat_test.start();
 		  mat_test.run();
 		  const_mat_test.matchingWeight();
 		  const_mat_test.matching(e);
 		  const_mat_test.matching(n);
 		  const MaxWeightedPerfectFractionalMatching<Graph>::MatchingMap& mmap =
 		    const_mat_test.matchingMap();
 		  e = mmap[n];
 		  const_mat_test.dualValue();
 		  const_mat_test.nodeValue(n);
+		}
 		void checkFractionalMatching(const SmartGraph& graph,
 		                             const MaxFractionalMatching<SmartGraph>& mfm,
 		                             bool allow_loops = true) {
 		  int pv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    int indeg = 0;
 		    for (InArcIt a(graph, n); a != INVALID; ++a) {
 		      if (mfm.matching(graph.source(a)) == a) {
 		        ++indeg;
+		      }
+		    }
 		    if (mfm.matching(n) != INVALID) {
 		      check(indeg == 1, "Invalid matching");
 		      ++pv;
 		    } else {
 		      check(indeg == 0, "Invalid matching");
+		    }
+		  }
 		  check(pv == mfm.matchingSize(), "Wrong matching size");
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    check((e == mfm.matching(graph.u(e)) ? 1 : 0) +
 		          (e == mfm.matching(graph.v(e)) ? 1 : 0) ==
 		          mfm.matching(e), "Invalid matching");
+		  }
 		  SmartGraph::NodeMap<bool> processed(graph, false);
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    if (processed[n]) continue;
 		    processed[n] = true;
 		    if (mfm.matching(n) == INVALID) continue;
 		    int num = 1;
 		    Node v = graph.target(mfm.matching(n));
 		    while (v != n) {
 		      processed[v] = true;
 		      ++num;
 		      v = graph.target(mfm.matching(v));
+		    }
 		    check(num == 2 || num % 2 == 1, "Wrong cycle size");
 		    check(allow_loops || num != 1, "Wrong cycle size");
+		  }
 		  int anum = 0, bnum = 0;
 		  SmartGraph::NodeMap<bool> neighbours(graph, false);
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    if (!mfm.barrier(n)) continue;
 		    ++anum;
 		    for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
 		      Node u = graph.source(a);
 		      if (!allow_loops && u == n) continue;
 		      if (!neighbours[u]) {
 		        neighbours[u] = true;
 		        ++bnum;
+		      }
+		    }
+		  }
 		  check(anum - bnum + mfm.matchingSize() == countNodes(graph),
 		        "Wrong barrier");
+		}
 		void checkPerfectFractionalMatching(const SmartGraph& graph,
 		                             const MaxFractionalMatching<SmartGraph>& mfm,
 		                             bool perfect, bool allow_loops = true) {
 		  if (perfect) {
 		    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		      int indeg = 0;
 		      for (InArcIt a(graph, n); a != INVALID; ++a) {
 		        if (mfm.matching(graph.source(a)) == a) {
 		          ++indeg;
+		        }
+		      }
 		      check(mfm.matching(n) != INVALID, "Invalid matching");
 		      check(indeg == 1, "Invalid matching");
+		    }
 		    for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		      check((e == mfm.matching(graph.u(e)) ? 1 : 0) +
 		            (e == mfm.matching(graph.v(e)) ? 1 : 0) ==
 		            mfm.matching(e), "Invalid matching");
+		    }
 		  } else {
 		    int anum = 0, bnum = 0;
 		    SmartGraph::NodeMap<bool> neighbours(graph, false);
 		    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		      if (!mfm.barrier(n)) continue;
 		      ++anum;
 		      for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
 		        Node u = graph.source(a);
 		        if (!allow_loops && u == n) continue;
 		        if (!neighbours[u]) {
 		          neighbours[u] = true;
 		          ++bnum;
+		        }
+		      }
+		    }
 		    check(anum - bnum > 0, "Wrong barrier");
+		  }
+		}
 		void checkWeightedFractionalMatching(const SmartGraph& graph,
 		                   const SmartGraph::EdgeMap<int>& weight,
 		                   const MaxWeightedFractionalMatching<SmartGraph>& mwfm,
 		                   bool allow_loops = true) {
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
 		    int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e))
 		      - weight[e] * mwfm.dualScale;
 		    check(rw >= 0, "Negative reduced weight");
 		    check(rw == 0 || !mwfm.matching(e),
 		          "Non-zero reduced weight on matching edge");
+		  }
 		  int pv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    int indeg = 0;
 		    for (InArcIt a(graph, n); a != INVALID; ++a) {
 		      if (mwfm.matching(graph.source(a)) == a) {
 		        ++indeg;
+		      }
+		    }
 		    check(indeg <= 1, "Invalid matching");
 		    if (mwfm.matching(n) != INVALID) {
 		      check(mwfm.nodeValue(n) >= 0, "Invalid node value");
 		      check(indeg == 1, "Invalid matching");
 		      pv += weight[mwfm.matching(n)];
 		      SmartGraph::Node o = graph.target(mwfm.matching(n));
 		    } else {
 		      check(mwfm.nodeValue(n) == 0, "Invalid matching");
 		      check(indeg == 0, "Invalid matching");
+		    }
+		  }
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    check((e == mwfm.matching(graph.u(e)) ? 1 : 0) +
 		          (e == mwfm.matching(graph.v(e)) ? 1 : 0) ==
 		          mwfm.matching(e), "Invalid matching");
+		  }
 		  int dv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    dv += mwfm.nodeValue(n);
+		  }
 		  check(pv * mwfm.dualScale == dv * 2, "Wrong duality");
 		  SmartGraph::NodeMap<bool> processed(graph, false);
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    if (processed[n]) continue;
 		    processed[n] = true;
 		    if (mwfm.matching(n) == INVALID) continue;
 		    int num = 1;
 		    Node v = graph.target(mwfm.matching(n));
 		    while (v != n) {
 		      processed[v] = true;
 		      ++num;
 		      v = graph.target(mwfm.matching(v));
+		    }
 		    check(num == 2 || num % 2 == 1, "Wrong cycle size");
 		    check(allow_loops || num != 1, "Wrong cycle size");
+		  }
 		  return;
+		}
 		void checkWeightedPerfectFractionalMatching(const SmartGraph& graph,
 		                const SmartGraph::EdgeMap<int>& weight,
 		                const MaxWeightedPerfectFractionalMatching<SmartGraph>& mwpfm,
 		                bool allow_loops = true) {
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
 		    int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e))
 		      - weight[e] * mwpfm.dualScale;
 		    check(rw >= 0, "Negative reduced weight");
 		    check(rw == 0 || !mwpfm.matching(e),
 		          "Non-zero reduced weight on matching edge");
+		  }
 		  int pv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    int indeg = 0;
 		    for (InArcIt a(graph, n); a != INVALID; ++a) {
 		      if (mwpfm.matching(graph.source(a)) == a) {
 		        ++indeg;
+		      }
+		    }
 		    check(mwpfm.matching(n) != INVALID, "Invalid perfect matching");
 		    check(indeg == 1, "Invalid perfect matching");
 		    pv += weight[mwpfm.matching(n)];
 		    SmartGraph::Node o = graph.target(mwpfm.matching(n));
+		  }
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    check((e == mwpfm.matching(graph.u(e)) ? 1 : 0) +
 		          (e == mwpfm.matching(graph.v(e)) ? 1 : 0) ==
 		          mwpfm.matching(e), "Invalid matching");
+		  }
 		  int dv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    dv += mwpfm.nodeValue(n);
+		  }
 		  check(pv * mwpfm.dualScale == dv * 2, "Wrong duality");
 		  SmartGraph::NodeMap<bool> processed(graph, false);
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    if (processed[n]) continue;
 		    processed[n] = true;
 		    if (mwpfm.matching(n) == INVALID) continue;
 		    int num = 1;
 		    Node v = graph.target(mwpfm.matching(n));
 		    while (v != n) {
 		      processed[v] = true;
 		      ++num;
 		      v = graph.target(mwpfm.matching(v));
+		    }
 		    check(num == 2 || num % 2 == 1, "Wrong cycle size");
 		    check(allow_loops || num != 1, "Wrong cycle size");
+		  }
 		  return;
+		}
 		int main() {
 		  for (int i = 0; i < lgfn; ++i) {
 		    SmartGraph graph;
 		    SmartGraph::EdgeMap<int> weight(graph);
 		    istringstream lgfs(lgf[i]);
 		    graphReader(graph, lgfs).
 		      edgeMap("weight", weight).run();
 		    bool perfect_with_loops;
+		    {
 		      MaxFractionalMatching<SmartGraph> mfm(graph, true);
 		      mfm.run();
 		      checkFractionalMatching(graph, mfm, true);
 		      perfect_with_loops = mfm.matchingSize() == countNodes(graph);
+		    }
 		    bool perfect_without_loops;
+		    {
 		      MaxFractionalMatching<SmartGraph> mfm(graph, false);
 		      mfm.run();
 		      checkFractionalMatching(graph, mfm, false);
 		      perfect_without_loops = mfm.matchingSize() == countNodes(graph);
+		    }
+		    {
 		      MaxFractionalMatching<SmartGraph> mfm(graph, true);
 		      bool result = mfm.runPerfect();
 		      checkPerfectFractionalMatching(graph, mfm, result, true);
 		      check(result == perfect_with_loops, "Wrong perfect matching");
+		    }
+		    {
 		      MaxFractionalMatching<SmartGraph> mfm(graph, false);
 		      bool result = mfm.runPerfect();
 		      checkPerfectFractionalMatching(graph, mfm, result, false);
 		      check(result == perfect_without_loops, "Wrong perfect matching");
+		    }
+		    {
 		      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, true);
 		      mwfm.run();
 		      checkWeightedFractionalMatching(graph, weight, mwfm, true);
+		    }
+		    {
 		      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, false);
 		      mwfm.run();
 		      checkWeightedFractionalMatching(graph, weight, mwfm, false);
+		    }
+		    {
 		      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
 		                                                             true);
 		      bool perfect = mwpfm.run();
 		      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
 		            "Perfect matching found");
 		      check(perfect == perfect_with_loops, "Wrong perfect matching");
 		      if (perfect) {
 		        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true);
+		      }
+		    }
+		    {
 		      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
 		                                                             false);
 		      bool perfect = mwpfm.run();
 		      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
 		            "Perfect matching found");
 		      check(perfect == perfect_without_loops, "Wrong perfect matching");
 		      if (perfect) {
 		        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false);
+		      }
+		    }
+		  }
 		  return 0;
+		}

doc/groups.dox

Show inline comments

@@ -333,384 +333,391 @@
 		 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
 		   from a source node when arc lenghts can be either positive or negative,
 		   but the digraph should not contain directed cycles with negative total
 		   length.
 		 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
 		   for solving the \e all-pairs \e shortest \e paths \e problem when arc
 		   lenghts can be either positive or negative, but the digraph should
 		   not contain directed cycles with negative total length.
 		 - \ref Suurballe A successive shortest path algorithm for finding
 		   arc-disjoint paths between two nodes having minimum total length.
 		*/
 		/**
 		@defgroup spantree Minimum Spanning Tree Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum cost spanning trees and arborescences.
 		This group contains the algorithms for finding minimum cost spanning
 		trees and arborescences \ref clrs01algorithms.
 		*/
 		/**
 		@defgroup max_flow Maximum Flow Algorithms
 		@ingroup algs
 		\brief Algorithms for finding maximum flows.
 		This group contains the algorithms for finding maximum flows and
 		feasible circulations \ref clrs01algorithms, \ref amo93networkflows.
 		The \e maximum \e flow \e problem is to find a flow of maximum value between
 		a single source and a single target. Formally, there is a \f$G=(V,A)\f$
 		digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
 		\f$s, t \in V\f$ source and target nodes.
 		A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
 		following optimization problem.
 		\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
 		\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
 		    \quad \forall u\in V\setminus\{s,t\} \f]
 		\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
 		LEMON contains several algorithms for solving maximum flow problems:
 		- \ref EdmondsKarp Edmonds-Karp algorithm
 		  \ref edmondskarp72theoretical.
 		- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm
 		  \ref goldberg88newapproach.
 		- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees
 		  \ref dinic70algorithm, \ref sleator83dynamic.
 		- \ref GoldbergTarjan !Preflow push-relabel algorithm with dynamic trees
 		  \ref goldberg88newapproach, \ref sleator83dynamic.
 		In most cases the \ref Preflow algorithm provides the
 		fastest method for computing a maximum flow. All implementations
 		also provide functions to query the minimum cut, which is the dual
 		problem of maximum flow.
 		\ref Circulation is a preflow push-relabel algorithm implemented directly
 		for finding feasible circulations, which is a somewhat different problem,
 		but it is strongly related to maximum flow.
 		For more information, see \ref Circulation.
 		*/
 		/**
 		@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum cost flows and circulations.
 		This group contains the algorithms for finding minimum cost flows and
 		circulations \ref amo93networkflows. For more information about this
 		problem and its dual solution, see \ref min_cost_flow
 		"Minimum Cost Flow Problem".
 		LEMON contains several algorithms for this problem.
 		 - \ref NetworkSimplex Primal Network Simplex algorithm with various
 		   pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex.
 		 - \ref CostScaling Cost Scaling algorithm based on push/augment and
 		   relabel operations \ref goldberg90approximation, \ref goldberg97efficient,
 		   \ref bunnagel98efficient.
 		 - \ref CapacityScaling Capacity Scaling algorithm based on the successive
 		   shortest path method \ref edmondskarp72theoretical.
 		 - \ref CycleCanceling Cycle-Canceling algorithms, two of which are
 		   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
 		In general NetworkSimplex is the most efficient implementation,
 		but in special cases other algorithms could be faster.
 		For example, if the total supply and/or capacities are rather small,
 		CapacityScaling is usually the fastest algorithm (without effective scaling).
 		*/
 		/**
 		@defgroup min_cut Minimum Cut Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum cut in graphs.
 		This group contains the algorithms for finding minimum cut in graphs.
 		The \e minimum \e cut \e problem is to find a non-empty and non-complete
 		\f$X\f$ subset of the nodes with minimum overall capacity on
 		outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a
 		\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
 		cut is the \f$X\f$ solution of the next optimization problem:
 		\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
 		    \sum_{uv\in A: u\in X, v\not\in X}cap(uv) \f]
 		LEMON contains several algorithms related to minimum cut problems:
 		- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
 		  in directed graphs.
 		- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
 		  calculating minimum cut in undirected graphs.
 		- \ref GomoryHu "Gomory-Hu tree computation" for calculating
 		  all-pairs minimum cut in undirected graphs.
 		If you want to find minimum cut just between two distinict nodes,
 		see the \ref max_flow "maximum flow problem".
 		*/
 		/**
 		@defgroup min_mean_cycle Minimum Mean Cycle Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum mean cycles.
 		This group contains the algorithms for finding minimum mean cycles
 		\ref clrs01algorithms, \ref amo93networkflows.
 		The \e minimum \e mean \e cycle \e problem is to find a directed cycle
 		of minimum mean length (cost) in a digraph.
 		The mean length of a cycle is the average length of its arcs, i.e. the
 		ratio between the total length of the cycle and the number of arcs on it.
 		This problem has an important connection to \e conservative \e length
 		\e functions, too. A length function on the arcs of a digraph is called
 		conservative if and only if there is no directed cycle of negative total
 		length. For an arbitrary length function, the negative of the minimum
 		cycle mean is the smallest \f$\epsilon\f$ value so that increasing the
 		arc lengths uniformly by \f$\epsilon\f$ results in a conservative length
 		function.
 		LEMON contains three algorithms for solving the minimum mean cycle problem:
 		- \ref Karp "Karp"'s original algorithm \ref amo93networkflows,
 		  \ref dasdan98minmeancycle.
 		- \ref HartmannOrlin "Hartmann-Orlin"'s algorithm, which is an improved
 		  version of Karp's algorithm \ref dasdan98minmeancycle.
 		- \ref Howard "Howard"'s policy iteration algorithm
 		  \ref dasdan98minmeancycle.
 		In practice, the Howard algorithm proved to be by far the most efficient
 		one, though the best known theoretical bound on its running time is
 		exponential.
 		Both Karp and HartmannOrlin algorithms run in time O(ne) and use space
 		O(n<sup>2</sup>+e), but the latter one is typically faster due to the
 		applied early termination scheme.
 		*/
 		/**
 		@defgroup matching Matching Algorithms
 		@ingroup algs
 		\brief Algorithms for finding matchings in graphs and bipartite graphs.
 		This group contains the algorithms for calculating
 		matchings in graphs and bipartite graphs. The general matching problem is
 		finding a subset of the edges for which each node has at most one incident
 		edge.
 		There are several different algorithms for calculate matchings in
 		graphs.  The matching problems in bipartite graphs are generally
 		easier than in general graphs. The goal of the matching optimization
 		can be finding maximum cardinality, maximum weight or minimum cost
 		matching. The search can be constrained to find perfect or
 		maximum cardinality matching.
 		The matching algorithms implemented in LEMON:
 		- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
 		  for calculating maximum cardinality matching in bipartite graphs.
 		- \ref PrBipartiteMatching Push-relabel algorithm
 		  for calculating maximum cardinality matching in bipartite graphs.
 		- \ref MaxWeightedBipartiteMatching
 		  Successive shortest path algorithm for calculating maximum weighted
 		  matching and maximum weighted bipartite matching in bipartite graphs.
 		- \ref MinCostMaxBipartiteMatching
 		  Successive shortest path algorithm for calculating minimum cost maximum
 		  matching in bipartite graphs.
 		- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
 		  maximum cardinality matching in general graphs.
 		- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
 		  maximum weighted matching in general graphs.
 		- \ref MaxWeightedPerfectMatching
 		  Edmond's blossom shrinking algorithm for calculating maximum weighted
 		  perfect matching in general graphs.
 		- \ref MaxFractionalMatching Push-relabel algorithm for calculating
 		  maximum cardinality fractional matching in general graphs.
 		- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
 		  maximum weighted fractional matching in general graphs.
 		- \ref MaxWeightedPerfectFractionalMatching
 		  Augmenting path algorithm for calculating maximum weighted
 		  perfect fractional matching in general graphs.
 		\image html matching.png
 		\image latex matching.eps "Min Cost Perfect Matching" width=\textwidth
 		*/
 		/**
 		@defgroup graph_properties Connectivity and Other Graph Properties
 		@ingroup algs
 		\brief Algorithms for discovering the graph properties
 		This group contains the algorithms for discovering the graph properties
 		like connectivity, bipartiteness, euler property, simplicity etc.
 		\image html connected_components.png
 		\image latex connected_components.eps "Connected components" width=\textwidth
 		*/
 		/**
 		@defgroup planar Planarity Embedding and Drawing
 		@ingroup algs
 		\brief Algorithms for planarity checking, embedding and drawing
 		This group contains the algorithms for planarity checking,
 		embedding and drawing.
 		\image html planar.png
 		\image latex planar.eps "Plane graph" width=\textwidth
 		*/
 		/**
 		@defgroup approx Approximation Algorithms
 		@ingroup algs
 		\brief Approximation algorithms.
 		This group contains the approximation and heuristic algorithms
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup auxalg Auxiliary Algorithms
 		@ingroup algs
 		\brief Auxiliary algorithms implemented in LEMON.
 		This group contains some algorithms implemented in LEMON
 		in order to make it easier to implement complex algorithms.
 		*/
 		/**
 		@defgroup gen_opt_group General Optimization Tools
 		\brief This group contains some general optimization frameworks
 		implemented in LEMON.
 		This group contains some general optimization frameworks
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup lp_group LP and MIP Solvers
 		@ingroup gen_opt_group
 		\brief LP and MIP solver interfaces for LEMON.
 		This group contains LP and MIP solver interfaces for LEMON.
 		Various LP solvers could be used in the same manner with this
 		high-level interface.
 		The currently supported solvers are \ref glpk, \ref clp, \ref cbc,
 		\ref cplex, \ref soplex.
 		*/
 		/**
 		@defgroup lp_utils Tools for Lp and Mip Solvers
 		@ingroup lp_group
 		\brief Helper tools to the Lp and Mip solvers.
 		This group adds some helper tools to general optimization framework
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup metah Metaheuristics
 		@ingroup gen_opt_group
 		\brief Metaheuristics for LEMON library.
 		This group contains some metaheuristic optimization tools.
 		*/
 		/**
 		@defgroup utils Tools and Utilities
 		\brief Tools and utilities for programming in LEMON
 		Tools and utilities for programming in LEMON.
 		*/
 		/**
 		@defgroup gutils Basic Graph Utilities
 		@ingroup utils
 		\brief Simple basic graph utilities.
 		This group contains some simple basic graph utilities.
 		*/
 		/**
 		@defgroup misc Miscellaneous Tools
 		@ingroup utils
 		\brief Tools for development, debugging and testing.
 		This group contains several useful tools for development,
 		debugging and testing.
 		*/
 		/**
 		@defgroup timecount Time Measuring and Counting
 		@ingroup misc
 		\brief Simple tools for measuring the performance of algorithms.
 		This group contains simple tools for measuring the performance
 		of algorithms.
 		*/
 		/**
 		@defgroup exceptions Exceptions
 		@ingroup utils
 		\brief Exceptions defined in LEMON.
 		This group contains the exceptions defined in LEMON.
 		*/
 		/**
 		@defgroup io_group Input-Output
 		\brief Graph Input-Output methods
 		This group contains the tools for importing and exporting graphs
 		and graph related data. Now it supports the \ref lgf-format
 		"LEMON Graph Format", the \c DIMACS format and the encapsulated
 		postscript (EPS) format.
 		*/
 		/**
 		@defgroup lemon_io LEMON Graph Format
 		@ingroup io_group
 		\brief Reading and writing LEMON Graph Format.
 		This group contains methods for reading and writing
 		\ref lgf-format "LEMON Graph Format".
 		*/
 		/**
 		@defgroup eps_io Postscript Exporting
 		@ingroup io_group
 		\brief General \c EPS drawer and graph exporter
 		This group contains general \c EPS drawing methods and special
 		graph exporting tools.
 		*/
 		/**
 		@defgroup dimacs_group DIMACS Format
 		@ingroup io_group
 		\brief Read and write files in DIMACS format
 		Tools to read a digraph from or write it to a file in DIMACS format data.
 		*/
 		/**
 		@defgroup nauty_group NAUTY Format
 		@ingroup io_group
 		\brief Read \e Nauty format
 		Tool to read graphs from \e Nauty format data.
 		*/
 		/**
 		@defgroup concept Concepts
 		\brief Skeleton classes and concept checking classes
 		This group contains the data/algorithm skeletons and concept checking
 		classes implemented in LEMON.
 		The purpose of the classes in this group is fourfold.
 		- These classes contain the documentations of the %concepts. In order
 		  to avoid document multiplications, an implementation of a concept
 		  simply refers to the corresponding concept class.
 		- These classes declare every functions, <tt>typedef</tt>s etc. an
 		  implementation of the %concepts should provide, however completely
 		  without implementations and real data structures behind the
 		  interface. On the other hand they should provide nothing else. All
 		  the algorithms working on a data structure meeting a certain concept
 		  should compile with these classes. (Though it will not run properly,
 		  of course.) In this way it is easily to check if an algorithm
 		  doesn't use any extra feature of a certain implementation.

lemon/Makefile.am

Show inline comments

 		EXTRA_DIST += \
 			lemon/lemon.pc.in \
 			lemon/CMakeLists.txt \
 			lemon/config.h.cmake
 		pkgconfig_DATA += lemon/lemon.pc
 		lib_LTLIBRARIES += lemon/libemon.la
 		lemon_libemon_la_SOURCES = \
 			lemon/arg_parser.cc \
 			lemon/base.cc \
 			lemon/color.cc \
 			lemon/lp_base.cc \
 			lemon/lp_skeleton.cc \
 			lemon/random.cc \
 			lemon/bits/windows.cc
 		nodist_lemon_HEADERS = lemon/config.h
 		lemon_libemon_la_CXXFLAGS = \
 			$(AM_CXXFLAGS) \
 			$(GLPK_CFLAGS) \
 			$(CPLEX_CFLAGS) \
 			$(SOPLEX_CXXFLAGS) \
 			$(CLP_CXXFLAGS) \
 			$(CBC_CXXFLAGS)
 		lemon_libemon_la_LDFLAGS = \
 			$(GLPK_LIBS) \
 			$(CPLEX_LIBS) \
 			$(SOPLEX_LIBS) \
 			$(CLP_LIBS) \
 			$(CBC_LIBS)
 		if HAVE_GLPK
 		lemon_libemon_la_SOURCES += lemon/glpk.cc
 		endif
 		if HAVE_CPLEX
 		lemon_libemon_la_SOURCES += lemon/cplex.cc
 		endif
 		if HAVE_SOPLEX
 		lemon_libemon_la_SOURCES += lemon/soplex.cc
 		endif
 		if HAVE_CLP
 		lemon_libemon_la_SOURCES += lemon/clp.cc
 		endif
 		if HAVE_CBC
 		lemon_libemon_la_SOURCES += lemon/cbc.cc
 		endif
 		lemon_HEADERS += \
 			lemon/adaptors.h \
 			lemon/arg_parser.h \
 			lemon/assert.h \
 			lemon/bellman_ford.h \
 			lemon/bfs.h \
 			lemon/bin_heap.h \
 			lemon/binomial_heap.h \
 			lemon/bucket_heap.h \
 			lemon/capacity_scaling.h \
 			lemon/cbc.h \
 			lemon/circulation.h \
 			lemon/clp.h \
 			lemon/color.h \
 			lemon/concept_check.h \
 			lemon/connectivity.h \
 			lemon/core.h \
 			lemon/cost_scaling.h \
 			lemon/counter.h \
 			lemon/cplex.h \
 			lemon/cycle_canceling.h \
 			lemon/dfs.h \
 			lemon/dheap.h \
 			lemon/dijkstra.h \
 			lemon/dim2.h \
 			lemon/dimacs.h \
 			lemon/edge_set.h \
 			lemon/elevator.h \
 			lemon/error.h \
 			lemon/euler.h \
 			lemon/fib_heap.h \
 			lemon/fractional_matching.h \
 			lemon/full_graph.h \
 			lemon/glpk.h \
 			lemon/gomory_hu.h \
 			lemon/graph_to_eps.h \
 			lemon/grid_graph.h \
 			lemon/hartmann_orlin_mmc.h \
 			lemon/howard_mmc.h \
 			lemon/hypercube_graph.h \
 			lemon/karp_mmc.h \
 			lemon/kruskal.h \
 			lemon/hao_orlin.h \
 			lemon/lgf_reader.h \
 			lemon/lgf_writer.h \
 			lemon/list_graph.h \
 			lemon/lp.h \
 			lemon/lp_base.h \
 			lemon/lp_skeleton.h \
 			lemon/maps.h \
 			lemon/matching.h \
 			lemon/math.h \
 			lemon/min_cost_arborescence.h \
 			lemon/nauty_reader.h \
 			lemon/network_simplex.h \
 			lemon/pairing_heap.h \
 			lemon/path.h \
 			lemon/planarity.h \
 			lemon/preflow.h \
 			lemon/quad_heap.h \
 			lemon/radix_heap.h \
 			lemon/radix_sort.h \
 			lemon/random.h \
 			lemon/smart_graph.h \
 			lemon/soplex.h \
 			lemon/static_graph.h \
 			lemon/suurballe.h \
 			lemon/time_measure.h \
 			lemon/tolerance.h \
 			lemon/unionfind.h \
 			lemon/bits/windows.h
 		bits_HEADERS += \
 			lemon/bits/alteration_notifier.h \
 			lemon/bits/array_map.h \
 			lemon/bits/bezier.h \
 			lemon/bits/default_map.h \
 			lemon/bits/edge_set_extender.h \
 			lemon/bits/enable_if.h \
 			lemon/bits/graph_adaptor_extender.h \
 			lemon/bits/graph_extender.h \
 			lemon/bits/map_extender.h \
 			lemon/bits/path_dump.h \
 			lemon/bits/solver_bits.h \
 			lemon/bits/traits.h \
 			lemon/bits/variant.h \
 			lemon/bits/vector_map.h
 		concept_HEADERS += \
 			lemon/concepts/digraph.h \
 			lemon/concepts/graph.h \
 			lemon/concepts/graph_components.h \
 			lemon/concepts/heap.h \
 			lemon/concepts/maps.h \
 			lemon/concepts/path.h

lemon/matching.h

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#ifndef LEMON_MAX_MATCHING_H
 		#define LEMON_MAX_MATCHING_H
 		#ifndef LEMON_MATCHING_H
 		#define LEMON_MATCHING_H
 		#include <vector>
 		#include <queue>
 		#include <set>
 		#include <limits>
 		#include <lemon/core.h>
 		#include <lemon/unionfind.h>
 		#include <lemon/bin_heap.h>
 		#include <lemon/maps.h>
 		#include <lemon/fractional_matching.h>
 		///\ingroup matching
 		///\file
 		///\brief Maximum matching algorithms in general graphs.
 		namespace lemon {
 		  /// \ingroup matching
 		  ///
 		  /// \brief Maximum cardinality matching in general graphs
 		  ///
 		  /// This class implements Edmonds' alternating forest matching algorithm
 		  /// for finding a maximum cardinality matching in a general undirected graph.
 		  /// It can be started from an arbitrary initial matching
 		  /// (the default is the empty one).
 		  ///
 		  /// The dual solution of the problem is a map of the nodes to
 		  /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
 		  /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
 		  /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
 		  /// with factor-critical components, the nodes in \c ODD/A form the
 		  /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
 		  /// a perfect matching. The number of the factor-critical components
 		  /// minus the number of barrier nodes is a lower bound on the
 		  /// unmatched nodes, and the matching is optimal if and only if this bound is
 		  /// tight. This decomposition can be obtained using \ref status() or
 		  /// \ref statusMap() after running the algorithm.
 		  ///
 		  /// \tparam GR The undirected graph type the algorithm runs on.
 		  template <typename GR>
 		  class MaxMatching {
 		  public:
 		    /// The graph type of the algorithm
 		    typedef GR Graph;
 		    /// The type of the matching map
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		    ///\brief Status constants for Gallai-Edmonds decomposition.
 		    ///
 		    ///These constants are used for indicating the Gallai-Edmonds
 		    ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
 		    ///induce a subgraph with factor-critical components, the nodes with
 		    ///status \c ODD (or \c A) form the canonical barrier, and the nodes
 		    ///with status \c MATCHED (or \c C) induce a subgraph having a
 		    ///perfect matching.
 		    enum Status {
 		      EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
 		      D = 1,
 		      MATCHED = 0,    ///< = 0. (\c C is an alias for \c MATCHED.)
 		      C = 0,
 		      ODD = -1,       ///< = -1. (\c A is an alias for \c ODD.)
 		      A = -1,
 		      UNMATCHED = -2  ///< = -2.
 		    };
 		    /// The type of the status map
 		    typedef typename Graph::template NodeMap<Status> StatusMap;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef UnionFindEnum<IntNodeMap> BlossomSet;
 		    typedef ExtendFindEnum<IntNodeMap> TreeSet;
 		    typedef RangeMap<Node> NodeIntMap;
 		    typedef MatchingMap EarMap;
 		    typedef std::vector<Node> NodeQueue;
 		    const Graph& _graph;
 		    MatchingMap* _matching;
 		    StatusMap* _status;
 		    EarMap* _ear;
 		    IntNodeMap* _blossom_set_index;
 		    BlossomSet* _blossom_set;
 		    NodeIntMap* _blossom_rep;
 		    IntNodeMap* _tree_set_index;
 		    TreeSet* _tree_set;
 		    NodeQueue _node_queue;
 		    int _process, _postpone, _last;
 		    int _node_num;
 		  private:
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_status) {
 		        _status = new StatusMap(_graph);
+		      }
 		      if (!_ear) {
 		        _ear = new EarMap(_graph);
+		      }
 		      if (!_blossom_set) {
 		        _blossom_set_index = new IntNodeMap(_graph);
 		        _blossom_set = new BlossomSet(*_blossom_set_index);
+		      }
 		      if (!_blossom_rep) {
 		        _blossom_rep = new NodeIntMap(_node_num);
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntNodeMap(_graph);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      _node_queue.resize(_node_num);
+		    }
 		    void destroyStructures() {
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_status) {
 		        delete _status;
+		      }
 		      if (_ear) {
 		        delete _ear;
+		      }
 		      if (_blossom_set) {
 		        delete _blossom_set;
 		        delete _blossom_set_index;
+		      }
 		      if (_blossom_rep) {
 		        delete _blossom_rep;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
+		    }
 		    void processDense(const Node& n) {
 		      _process = _postpone = _last = 0;
 		      _node_queue[_last++] = n;
 		      while (_process != _last) {
 		        Node u = _node_queue[_process++];
 		        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
 		          Node v = _graph.target(a);
 		          if ((*_status)[v] == MATCHED) {
 		            extendOnArc(a);
 		          } else if ((*_status)[v] == UNMATCHED) {
 		            augmentOnArc(a);
 		            return;
+		          }
+		        }
+		      }
 		      while (_postpone != _last) {
 		        Node u = _node_queue[_postpone++];
 		        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
 		          Node v = _graph.target(a);
 		          if ((*_status)[v] == EVEN) {
 		            if (_blossom_set->find(u) != _blossom_set->find(v)) {
 		              shrinkOnEdge(a);
+		            }
+		          }
 		          while (_process != _last) {
 		            Node w = _node_queue[_process++];
 		            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
 		              Node x = _graph.target(b);
 		              if ((*_status)[x] == MATCHED) {
 		                extendOnArc(b);
 		              } else if ((*_status)[x] == UNMATCHED) {
 		                augmentOnArc(b);
 		                return;
+		              }
+		            }
+		          }
+		        }
+		      }
+		    }
 		    void processSparse(const Node& n) {
 		      _process = _last = 0;
 		      _node_queue[_last++] = n;
 		      while (_process != _last) {
 		        Node u = _node_queue[_process++];
 		        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
 		          Node v = _graph.target(a);
 		          if ((*_status)[v] == EVEN) {
@@ -556,1875 +557,1910 @@
 		    /// \brief Return the size (cardinality) of the matching.
 		    ///
 		    /// This function returns the size (cardinality) of the current matching.
 		    /// After run() it returns the size of the maximum matching in the graph.
 		    int matchingSize() const {
 		      int size = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++size;
+		        }
+		      }
 		      return size / 2;
+		    }
 		    /// \brief Return \c true if the given edge is in the matching.
 		    ///
 		    /// This function returns \c true if the given edge is in the current
 		    /// matching.
 		    bool matching(const Edge& edge) const {
 		      return edge == (*_matching)[_graph.u(edge)];
+		    }
 		    /// \brief Return the matching arc (or edge) incident to the given node.
 		    ///
 		    /// This function returns the matching arc (or edge) incident to the
 		    /// given node in the current matching or \c INVALID if the node is
 		    /// not covered by the matching.
 		    Arc matching(const Node& n) const {
 		      return (*_matching)[n];
+		    }
 		    /// \brief Return a const reference to the matching map.
 		    ///
 		    /// This function returns a const reference to a node map that stores
 		    /// the matching arc (or edge) incident to each node.
 		    const MatchingMap& matchingMap() const {
 		      return *_matching;
+		    }
 		    /// \brief Return the mate of the given node.
 		    ///
 		    /// This function returns the mate of the given node in the current
 		    /// matching or \c INVALID if the node is not covered by the matching.
 		    Node mate(const Node& n) const {
 		      return (*_matching)[n] != INVALID ?
 		        _graph.target((*_matching)[n]) : INVALID;
+		    }
 		    /// @}
 		    /// \name Dual Solution
 		    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
 		    /// decomposition.
 		    /// @{
 		    /// \brief Return the status of the given node in the Edmonds-Gallai
 		    /// decomposition.
 		    ///
 		    /// This function returns the \ref Status "status" of the given node
 		    /// in the Edmonds-Gallai decomposition.
 		    Status status(const Node& n) const {
 		      return (*_status)[n];
+		    }
 		    /// \brief Return a const reference to the status map, which stores
 		    /// the Edmonds-Gallai decomposition.
 		    ///
 		    /// This function returns a const reference to a node map that stores the
 		    /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
 		    const StatusMap& statusMap() const {
 		      return *_status;
+		    }
 		    /// \brief Return \c true if the given node is in the barrier.
 		    ///
 		    /// This function returns \c true if the given node is in the barrier.
 		    bool barrier(const Node& n) const {
 		      return (*_status)[n] == ODD;
+		    }
 		    /// @}
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Weighted matching in general graphs
 		  ///
 		  /// This class provides an efficient implementation of Edmond's
 		  /// maximum weighted matching algorithm. The implementation is based
 		  /// on extensive use of priority queues and provides
 		  /// \f$O(nm\log n)\f$ time complexity.
 		  ///
 		  /// The maximum weighted matching problem is to find a subset of the
 		  /// edges in an undirected graph with maximum overall weight for which
 		  /// each node has at most one incident edge.
 		  /// It can be formulated with the following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
 		  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
 		      \quad \forall B\in\mathcal{O}\f] */
 		  /// \f[x_e \ge 0\quad \forall e\in E\f]
 		  /// \f[\max \sum_{e\in E}x_ew_e\f]
 		  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
 		  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
 		  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
 		  /// subsets of the nodes.
 		  ///
 		  /// The algorithm calculates an optimal matching and a proof of the
 		  /// optimality. The solution of the dual problem can be used to check
 		  /// the result of the algorithm. The dual linear problem is the
 		  /// following.
 		  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
 		      z_B \ge w_{uv} \quad \forall uv\in E\f] */
 		  /// \f[y_u \ge 0 \quad \forall u \in V\f]
 		  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
 		  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
 		      \frac{\vert B \vert - 1}{2}z_B\f] */
 		  ///
 		  /// The algorithm can be executed with the run() function.
 		  /// After it the matching (the primal solution) and the dual solution
 		  /// can be obtained using the query functions and the
 		  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
 		  /// which is able to iterate on the nodes of a blossom.
 		  /// If the value type is integer, then the dual solution is multiplied
 		  /// by \ref MaxWeightedMatching::dualScale "4".
 		  ///
 		  /// \tparam GR The undirected graph type the algorithm runs on.
 		  /// \tparam WM The type edge weight map. The default type is
 		  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 		#ifdef DOXYGEN
 		  template <typename GR, typename WM>
 		#else
 		  template <typename GR,
 		            typename WM = typename GR::template EdgeMap<int> >
 		#endif
 		  class MaxWeightedMatching {
 		  public:
 		    /// The graph type of the algorithm
 		    typedef GR Graph;
 		    /// The type of the edge weight map
 		    typedef WM WeightMap;
 		    /// The value type of the edge weights
 		    typedef typename WeightMap::Value Value;
 		    /// The type of the matching map
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		    /// \brief Scaling factor for dual solution
 		    ///
 		    /// Scaling factor for dual solution. It is equal to 4 or 1
 		    /// according to the value type.
 		    static const int dualScale =
 		      std::numeric_limits<Value>::is_integer ? 4 : 1;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef typename Graph::template NodeMap<Value> NodePotential;
 		    typedef std::vector<Node> BlossomNodeList;
 		    struct BlossomVariable {
 		      int begin, end;
 		      Value value;
 		      BlossomVariable(int _begin, int _end, Value _value)
 		        : begin(_begin), end(_end), value(_value) {}
 		    };
 		    typedef std::vector<BlossomVariable> BlossomPotential;
 		    const Graph& _graph;
 		    const WeightMap& _weight;
 		    MatchingMap* _matching;
 		    NodePotential* _node_potential;
 		    BlossomPotential _blossom_potential;
 		    BlossomNodeList _blossom_node_list;
 		    int _node_num;
 		    int _blossom_num;
 		    typedef RangeMap<int> IntIntMap;
 		    enum Status {
-		      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
 		      EVEN = -1, MATCHED = 0, ODD = 1
 		    };
 		    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
 		    struct BlossomData {
 		      int tree;
 		      Status status;
 		      Arc pred, next;
 		      Value pot, offset;
 		      Node base;
 		    };
 		    IntNodeMap *_blossom_index;
 		    BlossomSet *_blossom_set;
 		    RangeMap<BlossomData>* _blossom_data;
 		    IntNodeMap *_node_index;
 		    IntArcMap *_node_heap_index;
 		    struct NodeData {
 		      NodeData(IntArcMap& node_heap_index)
 		        : heap(node_heap_index) {}
 		      int blossom;
 		      Value pot;
 		      BinHeap<Value, IntArcMap> heap;
 		      std::map<int, Arc> heap_index;
 		      int tree;
 		    };
 		    RangeMap<NodeData>* _node_data;
 		    typedef ExtendFindEnum<IntIntMap> TreeSet;
 		    IntIntMap *_tree_set_index;
 		    TreeSet *_tree_set;
 		    IntNodeMap *_delta1_index;
 		    BinHeap<Value, IntNodeMap> *_delta1;
 		    IntIntMap *_delta2_index;
 		    BinHeap<Value, IntIntMap> *_delta2;
 		    IntEdgeMap *_delta3_index;
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    IntIntMap *_delta4_index;
 		    BinHeap<Value, IntIntMap> *_delta4;
 		    Value _delta_sum;
 		    int _unmatched;
 		    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
 		    FractionalMatching *_fractional;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_node_potential) {
 		        _node_potential = new NodePotential(_graph);
+		      }
 		      if (!_blossom_set) {
 		        _blossom_index = new IntNodeMap(_graph);
 		        _blossom_set = new BlossomSet(*_blossom_index);
 		        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
+		      }
 		      if (!_node_index) {
 		        _node_index = new IntNodeMap(_graph);
 		        _node_heap_index = new IntArcMap(_graph);
 		        _node_data = new RangeMap<NodeData>(_node_num,
 		                                              NodeData(*_node_heap_index));
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntIntMap(_blossom_num);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      if (!_delta1) {
 		        _delta1_index = new IntNodeMap(_graph);
 		        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
+		      }
 		      if (!_delta2) {
 		        _delta2_index = new IntIntMap(_blossom_num);
 		        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
+		      }
 		      if (!_delta3) {
 		        _delta3_index = new IntEdgeMap(_graph);
 		        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+		      }
 		      if (!_delta4) {
 		        _delta4_index = new IntIntMap(_blossom_num);
 		        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+		      }
+		    }
 		    void destroyStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_node_potential) {
 		        delete _node_potential;
+		      }
 		      if (_blossom_set) {
 		        delete _blossom_index;
 		        delete _blossom_set;
 		        delete _blossom_data;
+		      }
 		      if (_node_index) {
 		        delete _node_index;
 		        delete _node_heap_index;
 		        delete _node_data;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
 		      if (_delta1) {
 		        delete _delta1_index;
 		        delete _delta1;
+		      }
 		      if (_delta2) {
 		        delete _delta2_index;
 		        delete _delta2;
+		      }
 		      if (_delta3) {
 		        delete _delta3_index;
 		        delete _delta3;
+		      }
 		      if (_delta4) {
 		        delete _delta4_index;
 		        delete _delta4;
+		      }
+		    }
 		    void matchedToEven(int blossom, int tree) {
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      if (!_blossom_set->trivial(blossom)) {
 		        (*_blossom_data)[blossom].pot -=
 * (_delta_sum - (*_blossom_data)[blossom].offset);
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        _blossom_set->increase(n, std::numeric_limits<Value>::max());
 		        int ni = (*_node_index)[n];
 		        (*_node_data)[ni].heap.clear();
 		        (*_node_data)[ni].heap_index.clear();
 		        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
 		        _delta1->push(n, (*_node_data)[ni].pot);
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
 		              _delta3->push(e, rw / 2);
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP) {
 		              _delta3->push(e, rw);
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              if ((*_node_data)[vi].heap[it->second] > rw) {
 		                (*_node_data)[vi].heap.replace(it->second, e);
 		                (*_node_data)[vi].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[vi].heap.push(e, rw);
 		              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+		            }
 		            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
 		              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
 		              if ((*_blossom_data)[vb].status == MATCHED) {
 		                if (_delta2->state(vb) != _delta2->IN_HEAP) {
 		                  _delta2->push(vb, _blossom_set->classPrio(vb) -
 		                               (*_blossom_data)[vb].offset);
 		                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
 		                           (*_blossom_data)[vb].offset){
 		                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
 		                                   (*_blossom_data)[vb].offset);
+		                }
+		              }
+		            }
+		          }
+		        }
+		      }
 		      (*_blossom_data)[blossom].offset = 0;
+		    }
 		    void matchedToOdd(int blossom) {
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      (*_blossom_data)[blossom].offset += _delta_sum;
 		      if (!_blossom_set->trivial(blossom)) {
 		        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
 		                     (*_blossom_data)[blossom].offset);
+		      }
+		    }
 		    void evenToMatched(int blossom, int tree) {
 		      if (!_blossom_set->trivial(blossom)) {
 		        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        int ni = (*_node_index)[n];
 		        (*_node_data)[ni].pot -= _delta_sum;
 		        _delta1->erase(n);
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if (vb == blossom) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		          } else if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		            int vt = _tree_set->find(vb);
 		            if (vt != tree) {
 		              Arc r = _graph.oppositeArc(e);
 		              typename std::map<int, Arc>::iterator it =
 		                (*_node_data)[ni].heap_index.find(vt);
 		              if (it != (*_node_data)[ni].heap_index.end()) {
 		                if ((*_node_data)[ni].heap[it->second] > rw) {
 		                  (*_node_data)[ni].heap.replace(it->second, r);
 		                  (*_node_data)[ni].heap.decrease(r, rw);
 		                  it->second = r;
+		                }
 		              } else {
 		                (*_node_data)[ni].heap.push(r, rw);
 		                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+		              }
 		              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
 		                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
 		                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
 		                               (*_blossom_data)[blossom].offset);
 		                } else if ((*_delta2)[blossom] >
 		                           _blossom_set->classPrio(blossom) -
 		                           (*_blossom_data)[blossom].offset){
 		                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
 		                                   (*_blossom_data)[blossom].offset);
+		                }
+		              }
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              (*_node_data)[vi].heap.erase(it->second);
 		              (*_node_data)[vi].heap_index.erase(it);
 		              if ((*_node_data)[vi].heap.empty()) {
 		                _blossom_set->increase(v, std::numeric_limits<Value>::max());
 		              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
 		                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
+		              }
 		              if ((*_blossom_data)[vb].status == MATCHED) {
 		                if (_blossom_set->classPrio(vb) ==
 		                    std::numeric_limits<Value>::max()) {
 		                  _delta2->erase(vb);
 		                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
 		                           (*_blossom_data)[vb].offset) {
 		                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
 		                                   (*_blossom_data)[vb].offset);
+		                }
+		              }
+		            }
+		          }
+		        }
+		      }
+		    }
 		    void oddToMatched(int blossom) {
 		      (*_blossom_data)[blossom].offset -= _delta_sum;
 		      if (_blossom_set->classPrio(blossom) !=
 		          std::numeric_limits<Value>::max()) {
 		        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
 		                       (*_blossom_data)[blossom].offset);
+		      }
 		      if (!_blossom_set->trivial(blossom)) {
 		        _delta4->erase(blossom);
+		      }
+		    }
 		    void oddToEven(int blossom, int tree) {
 		      if (!_blossom_set->trivial(blossom)) {
 		        _delta4->erase(blossom);
 		        (*_blossom_data)[blossom].pot -=
 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        int ni = (*_node_index)[n];
 		        _blossom_set->increase(n, std::numeric_limits<Value>::max());
 		        (*_node_data)[ni].heap.clear();
 		        (*_node_data)[ni].heap_index.clear();
 		        (*_node_data)[ni].pot +=
 * _delta_sum - (*_blossom_data)[blossom].offset;
 		        _delta1->push(n, (*_node_data)[ni].pot);
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
 		              _delta3->push(e, rw / 2);
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP) {
 		              _delta3->push(e, rw);
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              if ((*_node_data)[vi].heap[it->second] > rw) {
 		                (*_node_data)[vi].heap.replace(it->second, e);
 		                (*_node_data)[vi].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[vi].heap.push(e, rw);
 		              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+		            }
 		            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
 		              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
 		              if ((*_blossom_data)[vb].status == MATCHED) {
 		                if (_delta2->state(vb) != _delta2->IN_HEAP) {
 		                  _delta2->push(vb, _blossom_set->classPrio(vb) -
 		                               (*_blossom_data)[vb].offset);
 		                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
 		                           (*_blossom_data)[vb].offset) {
 		                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
 		                                   (*_blossom_data)[vb].offset);
+		                }
+		              }
+		            }
+		          }
+		        }
+		      }
 		      (*_blossom_data)[blossom].offset = 0;
+		    }
 		    void matchedToUnmatched(int blossom) {
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        int ni = (*_node_index)[n];
 		        _blossom_set->increase(n, std::numeric_limits<Value>::max());
 		        (*_node_data)[ni].heap.clear();
 		        (*_node_data)[ni].heap_index.clear();
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.target(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP) {
 		              _delta3->push(e, rw);
+		            }
+		          }
+		        }
+		      }
+		    }
 		    void unmatchedToMatched(int blossom) {
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        int ni = (*_node_index)[n];
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if (vb == blossom) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		          } else if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		            int vt = _tree_set->find(vb);
 		            Arc r = _graph.oppositeArc(e);
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[ni].heap_index.find(vt);
 		            if (it != (*_node_data)[ni].heap_index.end()) {
 		              if ((*_node_data)[ni].heap[it->second] > rw) {
 		                (*_node_data)[ni].heap.replace(it->second, r);
 		                (*_node_data)[ni].heap.decrease(r, rw);
 		                it->second = r;
+		              }
 		            } else {
 		              (*_node_data)[ni].heap.push(r, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+		            }
 		            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
 		              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
 		                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
 		                             (*_blossom_data)[blossom].offset);
 		              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
 		                         (*_blossom_data)[blossom].offset){
 		                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
 		                                 (*_blossom_data)[blossom].offset);
+		              }
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
+		          }
+		        }
+		      }
+		    }
 		    void alternatePath(int even, int tree) {
 		      int odd;
 		      evenToMatched(even, tree);
 		      (*_blossom_data)[even].status = MATCHED;
 		      while ((*_blossom_data)[even].pred != INVALID) {
 		        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
 		        (*_blossom_data)[odd].status = MATCHED;
 		        oddToMatched(odd);
 		        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
 		        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
 		        (*_blossom_data)[even].status = MATCHED;
 		        evenToMatched(even, tree);
 		        (*_blossom_data)[even].next =
 		          _graph.oppositeArc((*_blossom_data)[odd].pred);
+		      }
+		    }
 		    void destroyTree(int tree) {
 		      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
 		        if ((*_blossom_data)[b].status == EVEN) {
 		          (*_blossom_data)[b].status = MATCHED;
 		          evenToMatched(b, tree);
 		        } else if ((*_blossom_data)[b].status == ODD) {
 		          (*_blossom_data)[b].status = MATCHED;
 		          oddToMatched(b);
+		        }
+		      }
 		      _tree_set->eraseClass(tree);
+		    }
 		    void unmatchNode(const Node& node) {
 		      int blossom = _blossom_set->find(node);
 		      int tree = _tree_set->find(blossom);
 		      alternatePath(blossom, tree);
 		      destroyTree(tree);
 		      (*_blossom_data)[blossom].status = UNMATCHED;
 		      (*_blossom_data)[blossom].base = node;
 		      matchedToUnmatched(blossom);
+		    }
+		      (*_blossom_data)[blossom].next = INVALID;
+		    }
 		    void augmentOnEdge(const Edge& edge) {
 		      int left = _blossom_set->find(_graph.u(edge));
 		      int right = _blossom_set->find(_graph.v(edge));
 		      if ((*_blossom_data)[left].status == EVEN) {
 		        int left_tree = _tree_set->find(left);
 		        alternatePath(left, left_tree);
 		        destroyTree(left_tree);
 		      } else {
 		        (*_blossom_data)[left].status = MATCHED;
 		        unmatchedToMatched(left);
+		      }
 		      if ((*_blossom_data)[right].status == EVEN) {
 		        int right_tree = _tree_set->find(right);
 		        alternatePath(right, right_tree);
 		        destroyTree(right_tree);
 		      } else {
 		        (*_blossom_data)[right].status = MATCHED;
 		        unmatchedToMatched(right);
+		      }
 		      (*_blossom_data)[left].next = _graph.direct(edge, true);
 		      (*_blossom_data)[right].next = _graph.direct(edge, false);
+		    }
 		    void augmentOnArc(const Arc& arc) {
 		      int left = _blossom_set->find(_graph.source(arc));
 		      int right = _blossom_set->find(_graph.target(arc));
 		      (*_blossom_data)[left].status = MATCHED;
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      (*_blossom_data)[left].next = arc;
 		      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
+		    }
 		    void extendOnArc(const Arc& arc) {
 		      int base = _blossom_set->find(_graph.target(arc));
 		      int tree = _tree_set->find(base);
 		      int odd = _blossom_set->find(_graph.source(arc));
 		      _tree_set->insert(odd, tree);
 		      (*_blossom_data)[odd].status = ODD;
 		      matchedToOdd(odd);
 		      (*_blossom_data)[odd].pred = arc;
 		      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
 		      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
 		      _tree_set->insert(even, tree);
 		      (*_blossom_data)[even].status = EVEN;
 		      matchedToEven(even, tree);
+		    }
 		    void shrinkOnEdge(const Edge& edge, int tree) {
 		      int nca = -1;
 		      std::vector<int> left_path, right_path;
+		      {
 		        std::set<int> left_set, right_set;
 		        int left = _blossom_set->find(_graph.u(edge));
 		        left_path.push_back(left);
 		        left_set.insert(left);
 		        int right = _blossom_set->find(_graph.v(edge));
 		        right_path.push_back(right);
 		        right_set.insert(right);
 		        while (true) {
 		          if ((*_blossom_data)[left].pred == INVALID) break;
 		          left =
 		            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
 		          left_path.push_back(left);
 		          left =
 		            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
 		          left_path.push_back(left);
 		          left_set.insert(left);
 		          if (right_set.find(left) != right_set.end()) {
 		            nca = left;
 		            break;
+		          }
 		          if ((*_blossom_data)[right].pred == INVALID) break;
 		          right =
 		            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
 		          right_path.push_back(right);
 		          right =
 		            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
 		          right_path.push_back(right);
 		          right_set.insert(right);
 		          if (left_set.find(right) != left_set.end()) {
 		            nca = right;
 		            break;
+		          }
+		        }
 		        if (nca == -1) {
 		          if ((*_blossom_data)[left].pred == INVALID) {
 		            nca = right;
 		            while (left_set.find(nca) == left_set.end()) {
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              right_path.push_back(nca);
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              right_path.push_back(nca);
+		            }
 		          } else {
 		            nca = left;
 		            while (right_set.find(nca) == right_set.end()) {
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              left_path.push_back(nca);
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              left_path.push_back(nca);
+		            }
+		          }
+		        }
+		      }
 		      std::vector<int> subblossoms;
 		      Arc prev;
 		      prev = _graph.direct(edge, true);
 		      for (int i = 0; left_path[i] != nca; i += 2) {
 		        subblossoms.push_back(left_path[i]);
 		        (*_blossom_data)[left_path[i]].next = prev;
 		        _tree_set->erase(left_path[i]);
 		        subblossoms.push_back(left_path[i + 1]);
 		        (*_blossom_data)[left_path[i + 1]].status = EVEN;
 		        oddToEven(left_path[i + 1], tree);
 		        _tree_set->erase(left_path[i + 1]);
 		        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
+		      }
 		      int k = 0;
 		      while (right_path[k] != nca) ++k;
 		      subblossoms.push_back(nca);
 		      (*_blossom_data)[nca].next = prev;
 		      for (int i = k - 2; i >= 0; i -= 2) {
 		        subblossoms.push_back(right_path[i + 1]);
 		        (*_blossom_data)[right_path[i + 1]].status = EVEN;
 		        oddToEven(right_path[i + 1], tree);
 		        _tree_set->erase(right_path[i + 1]);
 		        (*_blossom_data)[right_path[i + 1]].next =
 		          (*_blossom_data)[right_path[i + 1]].pred;
 		        subblossoms.push_back(right_path[i]);
 		        _tree_set->erase(right_path[i]);
+		      }
 		      int surface =
 		        _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+		        }
 		        (*_blossom_data)[subblossoms[i]].status = MATCHED;
+		      }
 		      (*_blossom_data)[surface].pot = -2 * _delta_sum;
 		      (*_blossom_data)[surface].offset = 0;
 		      (*_blossom_data)[surface].status = EVEN;
 		      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
 		      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
 		      _tree_set->insert(surface, tree);
 		      _tree_set->erase(nca);
+		    }
 		    void splitBlossom(int blossom) {
 		      Arc next = (*_blossom_data)[blossom].next;
 		      Arc pred = (*_blossom_data)[blossom].pred;
 		      int tree = _tree_set->find(blossom);
 		      (*_blossom_data)[blossom].status = MATCHED;
 		      oddToMatched(blossom);
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      std::vector<int> subblossoms;
 		      _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		      Value offset = (*_blossom_data)[blossom].offset;
 		      int b = _blossom_set->find(_graph.source(pred));
 		      int d = _blossom_set->find(_graph.source(next));
 		      int ib = -1, id = -1;
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (subblossoms[i] == b) ib = i;
 		        if (subblossoms[i] == d) id = i;
 		        (*_blossom_data)[subblossoms[i]].offset = offset;
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+		        }
 		        if (_blossom_set->classPrio(subblossoms[i]) !=
 		            std::numeric_limits<Value>::max()) {
 		          _delta2->push(subblossoms[i],
 		                        _blossom_set->classPrio(subblossoms[i]) -
 		                        (*_blossom_data)[subblossoms[i]].offset);
+		        }
+		      }
 		      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
 		        for (int i = (id + 1) % subblossoms.size();
 		             i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].pred = pred;
 		          (*_blossom_data)[sb].next =
 		                           _graph.oppositeArc((*_blossom_data)[tb].next);
 		          pred = (*_blossom_data)[ub].next;
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+		        }
 		        (*_blossom_data)[subblossoms[id]].status = ODD;
 		        matchedToOdd(subblossoms[id]);
 		        _tree_set->insert(subblossoms[id], tree);
 		        (*_blossom_data)[subblossoms[id]].next = next;
 		        (*_blossom_data)[subblossoms[id]].pred = pred;
 		      } else {
 		        for (int i = (ib + 1) % subblossoms.size();
 		             i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].next = next;
 		          (*_blossom_data)[sb].pred =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred =
 		            (*_blossom_data)[tb].next =
 		            _graph.oppositeArc((*_blossom_data)[ub].next);
 		          next = (*_blossom_data)[ub].next;
+		        }
 		        (*_blossom_data)[subblossoms[ib]].status = ODD;
 		        matchedToOdd(subblossoms[ib]);
 		        _tree_set->insert(subblossoms[ib], tree);
 		        (*_blossom_data)[subblossoms[ib]].next = next;
 		        (*_blossom_data)[subblossoms[ib]].pred = pred;
+		      }
 		      _tree_set->erase(blossom);
+		    }
 		    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
 		      if (_blossom_set->trivial(blossom)) {
 		        int bi = (*_node_index)[base];
 		        Value pot = (*_node_data)[bi].pot;
 		        (*_matching)[base] = matching;
 		        _blossom_node_list.push_back(base);
 		        (*_node_potential)[base] = pot;
 		      } else {
 		        Value pot = (*_blossom_data)[blossom].pot;
 		        int bn = _blossom_node_list.size();
 		        std::vector<int> subblossoms;
 		        _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		        int b = _blossom_set->find(base);
 		        int ib = -1;
 		        for (int i = 0; i < int(subblossoms.size()); ++i) {
 		          if (subblossoms[i] == b) { ib = i; break; }
+		        }
 		        for (int i = 1; i < int(subblossoms.size()); i += 2) {
 		          int sb = subblossoms[(ib + i) % subblossoms.size()];
 		          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
 		          Arc m = (*_blossom_data)[tb].next;
 		          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
 		          extractBlossom(tb, _graph.source(m), m);
+		        }
 		        extractBlossom(subblossoms[ib], base, matching);
 		        int en = _blossom_node_list.size();
 		        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
+		      }
+		    }
 		    void extractMatching() {
 		      std::vector<int> blossoms;
 		      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
 		        blossoms.push_back(c);
+		      }
 		      for (int i = 0; i < int(blossoms.size()); ++i) {
-		        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
+		        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
 		          Value offset = (*_blossom_data)[blossoms[i]].offset;
 		          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
 		          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
 		               n != INVALID; ++n) {
 		            (*_node_data)[(*_node_index)[n]].pot -= offset;
+		          }
 		          Arc matching = (*_blossom_data)[blossoms[i]].next;
 		          Node base = _graph.source(matching);
 		          extractBlossom(blossoms[i], base, matching);
 		        } else {
 		          Node base = (*_blossom_data)[blossoms[i]].base;
 		          extractBlossom(blossoms[i], base, INVALID);
+		        }
+		      }
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
 		      : _graph(graph), _weight(weight), _matching(0),
 		        _node_potential(0), _blossom_potential(), _blossom_node_list(),
 		        _node_num(0), _blossom_num(0),
 		        _blossom_index(0), _blossom_set(0), _blossom_data(0),
 		        _node_index(0), _node_heap_index(0), _node_data(0),
 		        _tree_set_index(0), _tree_set(0),
 		        _delta1_index(0), _delta1(0),
 		        _delta2_index(0), _delta2(0),
 		        _delta3_index(0), _delta3(0),
 		        _delta4_index(0), _delta4(0),
-		        _delta_sum() {}
+		        _delta_sum(), _unmatched(0),
 		        _fractional(0)
 		    {}
 		    ~MaxWeightedMatching() {
 		      destroyStructures();
 		      if (_fractional) {
 		        delete _fractional;
+		      }
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \ref run() member function.
 		    ///@{
 		    /// \brief Initialize the algorithm
 		    ///
 		    /// This function initializes the algorithm.
 		    void init() {
 		      createStructures();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = _node_num;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = 0;
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
 		            max = (dualScale * _weight[e]) / 2;
+		          }
+		        }
 		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = max;
 		        _delta1->push(n, max);
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        _tree_set->insert(blossom);
 		        (*_blossom_data)[blossom].status = EVEN;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = INVALID;
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int ti = (*_node_index)[_graph.v(e)];
 		        if (_graph.u(e) != _graph.v(e)) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
+		    }
 		    /// \brief Initialize the algorithm with fractional matching
 		    ///
 		    /// This function initializes the algorithm with a fractional
 		    /// matching. This initialization is also called jumpstart heuristic.
 		    void fractionalInit() {
 		      createStructures();
 		      if (_fractional == 0) {
 		        _fractional = new FractionalMatching(_graph, _weight, false);
+		      }
 		      _fractional->run();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = 0;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value pot = _fractional->nodeValue(n);
 		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = pot;
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        (*_blossom_data)[blossom].status = MATCHED;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = _fractional->matching(n);
 		        if (_fractional->matching(n) == INVALID) {
 		          (*_blossom_data)[blossom].base = n;
+		        }
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      typename Graph::template NodeMap<bool> processed(_graph, false);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if (processed[n]) continue;
 		        processed[n] = true;
 		        if (_fractional->matching(n) == INVALID) continue;
 		        int num = 1;
 		        Node v = _graph.target(_fractional->matching(n));
 		        while (n != v) {
 		          processed[v] = true;
 		          v = _graph.target(_fractional->matching(v));
 		          ++num;
+		        }
 		        if (num % 2 == 1) {
 		          std::vector<int> subblossoms(num);
 		          subblossoms[--num] = _blossom_set->find(n);
 		          _delta1->push(n, _fractional->nodeValue(n));
 		          v = _graph.target(_fractional->matching(n));
 		          while (n != v) {
 		            subblossoms[--num] = _blossom_set->find(v);
 		            _delta1->push(v, _fractional->nodeValue(v));
 		            v = _graph.target(_fractional->matching(v));
+		          }
 		          int surface =
 		            _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		          (*_blossom_data)[surface].status = EVEN;
 		          (*_blossom_data)[surface].pred = INVALID;
 		          (*_blossom_data)[surface].next = INVALID;
 		          (*_blossom_data)[surface].pot = 0;
 		          (*_blossom_data)[surface].offset = 0;
 		          _tree_set->insert(surface);
 		          ++_unmatched;
+		        }
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int sb = _blossom_set->find(_graph.u(e));
 		        int ti = (*_node_index)[_graph.v(e)];
 		        int tb = _blossom_set->find(_graph.v(e));
 		        if ((*_blossom_data)[sb].status == EVEN &&
 		            (*_blossom_data)[tb].status == EVEN && sb != tb) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        int nb = _blossom_set->find(n);
 		        if ((*_blossom_data)[nb].status != MATCHED) continue;
 		        int ni = (*_node_index)[n];
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.target(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            int vt = _tree_set->find(vb);
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[ni].heap_index.find(vt);
 		            if (it != (*_node_data)[ni].heap_index.end()) {
 		              if ((*_node_data)[ni].heap[it->second] > rw) {
 		                (*_node_data)[ni].heap.replace(it->second, e);
 		                (*_node_data)[ni].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[ni].heap.push(e, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+		            }
+		          }
+		        }
 		        if (!(*_node_data)[ni].heap.empty()) {
 		          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		          _delta2->push(nb, _blossom_set->classPrio(nb));
+		        }
+		      }
+		    }
 		    /// \brief Start the algorithm
 		    ///
 		    /// This function starts the algorithm.
 		    ///
-		    /// \pre \ref init() must be called before using this function.
+		    /// \pre \ref init() or \ref fractionalInit() must be called
 		    /// before using this function.
 		    void start() {
 		      enum OpType {
 		        D1, D2, D3, D4
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		      while (_unmatched > 0) {
 		        Value d1 = !_delta1->empty() ?
 		          _delta1->prio() : std::numeric_limits<Value>::max();
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
 		        Value d4 = !_delta4->empty() ?
 		          _delta4->prio() : std::numeric_limits<Value>::max();
-		        _delta_sum = d1; OpType ot = D1;
+		        _delta_sum = d3; OpType ot = D3;
 		        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
 		        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
 		        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
 		        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
 		        switch (ot) {
 		        case D1:
+		          {
 		            Node n = _delta1->top();
 		            unmatchNode(n);
-		            --unmatched;
+		            --_unmatched;
+		          }
 		          break;
 		        case D2:
+		          {
 		            int blossom = _delta2->top();
 		            Node n = _blossom_set->classTop(blossom);
 		            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
 		            extendOnArc(e);
 		            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
 		            if ((*_blossom_data)[blossom].next == INVALID) {
 		              augmentOnArc(a);
 		              --_unmatched;
 		            } else {
 		              extendOnArc(a);
+		            }
+		          }
 		          break;
 		        case D3:
+		          {
 		            Edge e = _delta3->top();
 		            int left_blossom = _blossom_set->find(_graph.u(e));
 		            int right_blossom = _blossom_set->find(_graph.v(e));
 		            if (left_blossom == right_blossom) {
 		              _delta3->pop();
 		            } else {
 		              int left_tree;
 		              if ((*_blossom_data)[left_blossom].status == EVEN) {
 		                left_tree = _tree_set->find(left_blossom);
 		              } else {
 		                left_tree = -1;
 		                ++unmatched;
+		              }
 		              int right_tree;
 		              if ((*_blossom_data)[right_blossom].status == EVEN) {
 		                right_tree = _tree_set->find(right_blossom);
 		              } else {
 		                right_tree = -1;
 		                ++unmatched;
+		              }
 		              int left_tree = _tree_set->find(left_blossom);
 		              int right_tree = _tree_set->find(right_blossom);
 		              if (left_tree == right_tree) {
 		                shrinkOnEdge(e, left_tree);
 		              } else {
 		                augmentOnEdge(e);
-		                unmatched -= 2;
+		                _unmatched -= 2;
+		              }
+		            }
 		          } break;
 		        case D4:
 		          splitBlossom(_delta4->top());
 		          break;
+		        }
+		      }
 		      extractMatching();
+		    }
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This method runs the \c %MaxWeightedMatching algorithm.
 		    ///
 		    /// \note mwm.run() is just a shortcut of the following code.
 		    /// \code
-		    ///   mwm.init();
+		    ///   mwm.fractionalInit();
 		    ///   mwm.start();
 		    /// \endcode
 		    void run() {
-		      init();
+		      fractionalInit();
 		      start();
+		    }
 		    /// @}
 		    /// \name Primal Solution
 		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// matching.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the weight of the matching.
 		    ///
 		    /// This function returns the weight of the found matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value matchingWeight() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          sum += _weight[(*_matching)[n]];
+		        }
+		      }
-		      return sum /= 2;
 		      return sum / 2;
+		    }
 		    /// \brief Return the size (cardinality) of the matching.
 		    ///
 		    /// This function returns the size (cardinality) of the found matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    int matchingSize() const {
 		      int num = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++num;
+		        }
+		      }
 		      return num /= 2;
+		    }
 		    /// \brief Return \c true if the given edge is in the matching.
 		    ///
 		    /// This function returns \c true if the given edge is in the found
 		    /// matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    bool matching(const Edge& edge) const {
 		      return edge == (*_matching)[_graph.u(edge)];
+		    }
 		    /// \brief Return the matching arc (or edge) incident to the given node.
 		    ///
 		    /// This function returns the matching arc (or edge) incident to the
 		    /// given node in the found matching or \c INVALID if the node is
 		    /// not covered by the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Return a const reference to the matching map.
 		    ///
 		    /// This function returns a const reference to a node map that stores
 		    /// the matching arc (or edge) incident to each node.
 		    const MatchingMap& matchingMap() const {
 		      return *_matching;
+		    }
 		    /// \brief Return the mate of the given node.
 		    ///
 		    /// This function returns the mate of the given node in the found
 		    /// matching or \c INVALID if the node is not covered by the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Node mate(const Node& node) const {
 		      return (*_matching)[node] != INVALID ?
 		        _graph.target((*_matching)[node]) : INVALID;
+		    }
 		    /// @}
 		    /// \name Dual Solution
 		    /// Functions to get the dual solution.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the value of the dual solution.
 		    ///
 		    /// This function returns the value of the dual solution.
 		    /// It should be equal to the primal value scaled by \ref dualScale
 		    /// "dual scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value dualValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        sum += nodeValue(n);
+		      }
 		      for (int i = 0; i < blossomNum(); ++i) {
 		        sum += blossomValue(i) * (blossomSize(i) / 2);
+		      }
 		      return sum;
+		    }
 		    /// \brief Return the dual value (potential) of the given node.
 		    ///
 		    /// This function returns the dual value (potential) of the given node.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value nodeValue(const Node& n) const {
 		      return (*_node_potential)[n];
+		    }
 		    /// \brief Return the number of the blossoms in the basis.
 		    ///
 		    /// This function returns the number of the blossoms in the basis.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    /// \see BlossomIt
 		    int blossomNum() const {
 		      return _blossom_potential.size();
+		    }
 		    /// \brief Return the number of the nodes in the given blossom.
 		    ///
 		    /// This function returns the number of the nodes in the given blossom.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    /// \see BlossomIt
 		    int blossomSize(int k) const {
 		      return _blossom_potential[k].end - _blossom_potential[k].begin;
+		    }
 		    /// \brief Return the dual value (ptential) of the given blossom.
 		    ///
 		    /// This function returns the dual value (ptential) of the given blossom.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value blossomValue(int k) const {
 		      return _blossom_potential[k].value;
+		    }
 		    /// \brief Iterator for obtaining the nodes of a blossom.
 		    ///
 		    /// This class provides an iterator for obtaining the nodes of the
 		    /// given blossom. It lists a subset of the nodes.
 		    /// Before using this iterator, you must allocate a
 		    /// MaxWeightedMatching class and execute it.
 		    class BlossomIt {
 		    public:
 		      /// \brief Constructor.
 		      ///
 		      /// Constructor to get the nodes of the given variable.
 		      ///
 		      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
 		      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
 		      /// called before initializing this iterator.
 		      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
 		        : _algorithm(&algorithm)
+		      {
 		        _index = _algorithm->_blossom_potential[variable].begin;
 		        _last = _algorithm->_blossom_potential[variable].end;
+		      }
 		      /// \brief Conversion to \c Node.
 		      ///
 		      /// Conversion to \c Node.
 		      operator Node() const {
 		        return _algorithm->_blossom_node_list[_index];
+		      }
 		      /// \brief Increment operator.
 		      ///
 		      /// Increment operator.
 		      BlossomIt& operator++() {
 		        ++_index;
 		        return *this;
+		      }
 		      /// \brief Validity checking
 		      ///
 		      /// Checks whether the iterator is invalid.
 		      bool operator==(Invalid) const { return _index == _last; }
 		      /// \brief Validity checking
 		      ///
 		      /// Checks whether the iterator is valid.
 		      bool operator!=(Invalid) const { return _index != _last; }
 		    private:
 		      const MaxWeightedMatching* _algorithm;
 		      int _last;
 		      int _index;
 		    };
 		    /// @}
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Weighted perfect matching in general graphs
 		  ///
 		  /// This class provides an efficient implementation of Edmond's
 		  /// maximum weighted perfect matching algorithm. The implementation
 		  /// is based on extensive use of priority queues and provides
 		  /// \f$O(nm\log n)\f$ time complexity.
 		  ///
 		  /// The maximum weighted perfect matching problem is to find a subset of
 		  /// the edges in an undirected graph with maximum overall weight for which
 		  /// each node has exactly one incident edge.
 		  /// It can be formulated with the following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
 		  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
 		      \quad \forall B\in\mathcal{O}\f] */
 		  /// \f[x_e \ge 0\quad \forall e\in E\f]
 		  /// \f[\max \sum_{e\in E}x_ew_e\f]
 		  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
 		  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
 		  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
 		  /// subsets of the nodes.
 		  ///
 		  /// The algorithm calculates an optimal matching and a proof of the
 		  /// optimality. The solution of the dual problem can be used to check
 		  /// the result of the algorithm. The dual linear problem is the
 		  /// following.
 		  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
 		      w_{uv} \quad \forall uv\in E\f] */
 		  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
 		  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
 		      \frac{\vert B \vert - 1}{2}z_B\f] */
 		  ///
 		  /// The algorithm can be executed with the run() function.
 		  /// After it the matching (the primal solution) and the dual solution
 		  /// can be obtained using the query functions and the
 		  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
 		  /// which is able to iterate on the nodes of a blossom.
 		  /// If the value type is integer, then the dual solution is multiplied
 		  /// by \ref MaxWeightedMatching::dualScale "4".
 		  ///
 		  /// \tparam GR The undirected graph type the algorithm runs on.
 		  /// \tparam WM The type edge weight map. The default type is
 		  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 		#ifdef DOXYGEN
 		  template <typename GR, typename WM>
 		#else
 		  template <typename GR,
 		            typename WM = typename GR::template EdgeMap<int> >
 		#endif
 		  class MaxWeightedPerfectMatching {
 		  public:
 		    /// The graph type of the algorithm
 		    typedef GR Graph;
 		    /// The type of the edge weight map
 		    typedef WM WeightMap;
 		    /// The value type of the edge weights
 		    typedef typename WeightMap::Value Value;
 		    /// \brief Scaling factor for dual solution
 		    ///
 		    /// Scaling factor for dual solution, it is equal to 4 or 1
 		    /// according to the value type.
 		    static const int dualScale =
 		      std::numeric_limits<Value>::is_integer ? 4 : 1;
 		    /// The type of the matching map
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef typename Graph::template NodeMap<Value> NodePotential;
 		    typedef std::vector<Node> BlossomNodeList;
 		    struct BlossomVariable {
 		      int begin, end;
 		      Value value;
 		      BlossomVariable(int _begin, int _end, Value _value)
 		        : begin(_begin), end(_end), value(_value) {}
 		    };
 		    typedef std::vector<BlossomVariable> BlossomPotential;
 		    const Graph& _graph;
 		    const WeightMap& _weight;
 		    MatchingMap* _matching;
 		    NodePotential* _node_potential;
 		    BlossomPotential _blossom_potential;
 		    BlossomNodeList _blossom_node_list;
 		    int _node_num;
 		    int _blossom_num;
 		    typedef RangeMap<int> IntIntMap;
 		    enum Status {
 		      EVEN = -1, MATCHED = 0, ODD = 1
 		    };
 		    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
 		    struct BlossomData {
 		      int tree;
 		      Status status;
 		      Arc pred, next;
 		      Value pot, offset;
 		    };
 		    IntNodeMap *_blossom_index;
 		    BlossomSet *_blossom_set;
 		    RangeMap<BlossomData>* _blossom_data;
 		    IntNodeMap *_node_index;
 		    IntArcMap *_node_heap_index;
 		    struct NodeData {
 		      NodeData(IntArcMap& node_heap_index)
 		        : heap(node_heap_index) {}
 		      int blossom;
 		      Value pot;
 		      BinHeap<Value, IntArcMap> heap;
 		      std::map<int, Arc> heap_index;
 		      int tree;
 		    };
 		    RangeMap<NodeData>* _node_data;
 		    typedef ExtendFindEnum<IntIntMap> TreeSet;
 		    IntIntMap *_tree_set_index;
 		    TreeSet *_tree_set;
 		    IntIntMap *_delta2_index;
 		    BinHeap<Value, IntIntMap> *_delta2;
 		    IntEdgeMap *_delta3_index;
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    IntIntMap *_delta4_index;
 		    BinHeap<Value, IntIntMap> *_delta4;
 		    Value _delta_sum;
 		    int _unmatched;
 		    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
 		    FractionalMatching;
 		    FractionalMatching *_fractional;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_node_potential) {
 		        _node_potential = new NodePotential(_graph);
+		      }
 		      if (!_blossom_set) {
 		        _blossom_index = new IntNodeMap(_graph);
 		        _blossom_set = new BlossomSet(*_blossom_index);
 		        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
+		      }
 		      if (!_node_index) {
 		        _node_index = new IntNodeMap(_graph);
 		        _node_heap_index = new IntArcMap(_graph);
 		        _node_data = new RangeMap<NodeData>(_node_num,
 		                                            NodeData(*_node_heap_index));
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntIntMap(_blossom_num);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      if (!_delta2) {
 		        _delta2_index = new IntIntMap(_blossom_num);
 		        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
+		      }
 		      if (!_delta3) {
 		        _delta3_index = new IntEdgeMap(_graph);
 		        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+		      }
 		      if (!_delta4) {
 		        _delta4_index = new IntIntMap(_blossom_num);
 		        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+		      }
+		    }
 		    void destroyStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_node_potential) {
 		        delete _node_potential;
+		      }
 		      if (_blossom_set) {
 		        delete _blossom_index;
 		        delete _blossom_set;
 		        delete _blossom_data;
+		      }
 		      if (_node_index) {
 		        delete _node_index;
 		        delete _node_heap_index;
 		        delete _node_data;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
 		      if (_delta2) {
 		        delete _delta2_index;
 		        delete _delta2;
+		      }
 		      if (_delta3) {
 		        delete _delta3_index;
 		        delete _delta3;
+		      }
 		      if (_delta4) {
 		        delete _delta4_index;
 		        delete _delta4;
+		      }
+		    }
 		    void matchedToEven(int blossom, int tree) {
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      if (!_blossom_set->trivial(blossom)) {
 		        (*_blossom_data)[blossom].pot -=
 * (_delta_sum - (*_blossom_data)[blossom].offset);
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        _blossom_set->increase(n, std::numeric_limits<Value>::max());
 		        int ni = (*_node_index)[n];
 		        (*_node_data)[ni].heap.clear();
 		        (*_node_data)[ni].heap_index.clear();
 		        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
 		              _delta3->push(e, rw / 2);
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              if ((*_node_data)[vi].heap[it->second] > rw) {
 		                (*_node_data)[vi].heap.replace(it->second, e);
 		                (*_node_data)[vi].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[vi].heap.push(e, rw);
 		              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+		            }
 		            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
 		              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
 		              if ((*_blossom_data)[vb].status == MATCHED) {
 		                if (_delta2->state(vb) != _delta2->IN_HEAP) {
 		                  _delta2->push(vb, _blossom_set->classPrio(vb) -
 		                               (*_blossom_data)[vb].offset);
 		                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
 		                           (*_blossom_data)[vb].offset){
 		                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
 		                                   (*_blossom_data)[vb].offset);
+		                }
+		              }
+		            }
+		          }
+		        }
+		      }
 		      (*_blossom_data)[blossom].offset = 0;
+		    }
 		    void matchedToOdd(int blossom) {
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      (*_blossom_data)[blossom].offset += _delta_sum;
 		      if (!_blossom_set->trivial(blossom)) {
 		        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
 		                     (*_blossom_data)[blossom].offset);
+		      }
+		    }
 		    void evenToMatched(int blossom, int tree) {
 		      if (!_blossom_set->trivial(blossom)) {
 		        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        int ni = (*_node_index)[n];
 		        (*_node_data)[ni].pot -= _delta_sum;
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if (vb == blossom) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		          } else if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		            int vt = _tree_set->find(vb);
 		            if (vt != tree) {
 		              Arc r = _graph.oppositeArc(e);
 		              typename std::map<int, Arc>::iterator it =
 		                (*_node_data)[ni].heap_index.find(vt);
 		              if (it != (*_node_data)[ni].heap_index.end()) {
 		                if ((*_node_data)[ni].heap[it->second] > rw) {
 		                  (*_node_data)[ni].heap.replace(it->second, r);
 		                  (*_node_data)[ni].heap.decrease(r, rw);
 		                  it->second = r;
+		                }
 		              } else {
 		                (*_node_data)[ni].heap.push(r, rw);
 		                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+		              }
 		              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
 		                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
 		                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
 		                               (*_blossom_data)[blossom].offset);
 		                } else if ((*_delta2)[blossom] >
 		                           _blossom_set->classPrio(blossom) -
 		                           (*_blossom_data)[blossom].offset){
 		                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
 		                                   (*_blossom_data)[blossom].offset);
+		                }
+		              }
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              (*_node_data)[vi].heap.erase(it->second);
 		              (*_node_data)[vi].heap_index.erase(it);
 		              if ((*_node_data)[vi].heap.empty()) {
 		                _blossom_set->increase(v, std::numeric_limits<Value>::max());
 		              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
 		                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
+		              }
 		              if ((*_blossom_data)[vb].status == MATCHED) {
@@ -2719,526 +2755,668 @@
 		      _tree_set->insert(surface, tree);
 		      _tree_set->erase(nca);
+		    }
 		    void splitBlossom(int blossom) {
 		      Arc next = (*_blossom_data)[blossom].next;
 		      Arc pred = (*_blossom_data)[blossom].pred;
 		      int tree = _tree_set->find(blossom);
 		      (*_blossom_data)[blossom].status = MATCHED;
 		      oddToMatched(blossom);
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      std::vector<int> subblossoms;
 		      _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		      Value offset = (*_blossom_data)[blossom].offset;
 		      int b = _blossom_set->find(_graph.source(pred));
 		      int d = _blossom_set->find(_graph.source(next));
 		      int ib = -1, id = -1;
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (subblossoms[i] == b) ib = i;
 		        if (subblossoms[i] == d) id = i;
 		        (*_blossom_data)[subblossoms[i]].offset = offset;
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+		        }
 		        if (_blossom_set->classPrio(subblossoms[i]) !=
 		            std::numeric_limits<Value>::max()) {
 		          _delta2->push(subblossoms[i],
 		                        _blossom_set->classPrio(subblossoms[i]) -
 		                        (*_blossom_data)[subblossoms[i]].offset);
+		        }
+		      }
 		      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
 		        for (int i = (id + 1) % subblossoms.size();
 		             i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].pred = pred;
 		          (*_blossom_data)[sb].next =
 		                           _graph.oppositeArc((*_blossom_data)[tb].next);
 		          pred = (*_blossom_data)[ub].next;
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+		        }
 		        (*_blossom_data)[subblossoms[id]].status = ODD;
 		        matchedToOdd(subblossoms[id]);
 		        _tree_set->insert(subblossoms[id], tree);
 		        (*_blossom_data)[subblossoms[id]].next = next;
 		        (*_blossom_data)[subblossoms[id]].pred = pred;
 		      } else {
 		        for (int i = (ib + 1) % subblossoms.size();
 		             i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].next = next;
 		          (*_blossom_data)[sb].pred =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred =
 		            (*_blossom_data)[tb].next =
 		            _graph.oppositeArc((*_blossom_data)[ub].next);
 		          next = (*_blossom_data)[ub].next;
+		        }
 		        (*_blossom_data)[subblossoms[ib]].status = ODD;
 		        matchedToOdd(subblossoms[ib]);
 		        _tree_set->insert(subblossoms[ib], tree);
 		        (*_blossom_data)[subblossoms[ib]].next = next;
 		        (*_blossom_data)[subblossoms[ib]].pred = pred;
+		      }
 		      _tree_set->erase(blossom);
+		    }
 		    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
 		      if (_blossom_set->trivial(blossom)) {
 		        int bi = (*_node_index)[base];
 		        Value pot = (*_node_data)[bi].pot;
 		        (*_matching)[base] = matching;
 		        _blossom_node_list.push_back(base);
 		        (*_node_potential)[base] = pot;
 		      } else {
 		        Value pot = (*_blossom_data)[blossom].pot;
 		        int bn = _blossom_node_list.size();
 		        std::vector<int> subblossoms;
 		        _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		        int b = _blossom_set->find(base);
 		        int ib = -1;
 		        for (int i = 0; i < int(subblossoms.size()); ++i) {
 		          if (subblossoms[i] == b) { ib = i; break; }
+		        }
 		        for (int i = 1; i < int(subblossoms.size()); i += 2) {
 		          int sb = subblossoms[(ib + i) % subblossoms.size()];
 		          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
 		          Arc m = (*_blossom_data)[tb].next;
 		          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
 		          extractBlossom(tb, _graph.source(m), m);
+		        }
 		        extractBlossom(subblossoms[ib], base, matching);
 		        int en = _blossom_node_list.size();
 		        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
+		      }
+		    }
 		    void extractMatching() {
 		      std::vector<int> blossoms;
 		      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
 		        blossoms.push_back(c);
+		      }
 		      for (int i = 0; i < int(blossoms.size()); ++i) {
 		        Value offset = (*_blossom_data)[blossoms[i]].offset;
 		        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
 		        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
 		             n != INVALID; ++n) {
 		          (*_node_data)[(*_node_index)[n]].pot -= offset;
+		        }
 		        Arc matching = (*_blossom_data)[blossoms[i]].next;
 		        Node base = _graph.source(matching);
 		        extractBlossom(blossoms[i], base, matching);
+		      }
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
 		      : _graph(graph), _weight(weight), _matching(0),
 		        _node_potential(0), _blossom_potential(), _blossom_node_list(),
 		        _node_num(0), _blossom_num(0),
 		        _blossom_index(0), _blossom_set(0), _blossom_data(0),
 		        _node_index(0), _node_heap_index(0), _node_data(0),
 		        _tree_set_index(0), _tree_set(0),
 		        _delta2_index(0), _delta2(0),
 		        _delta3_index(0), _delta3(0),
 		        _delta4_index(0), _delta4(0),
-		        _delta_sum() {}
+		        _delta_sum(), _unmatched(0),
 		        _fractional(0)
 		    {}
 		    ~MaxWeightedPerfectMatching() {
 		      destroyStructures();
 		      if (_fractional) {
 		        delete _fractional;
+		      }
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \ref run() member function.
 		    ///@{
 		    /// \brief Initialize the algorithm
 		    ///
 		    /// This function initializes the algorithm.
 		    void init() {
 		      createStructures();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = _node_num;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = - std::numeric_limits<Value>::max();
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
 		            max = (dualScale * _weight[e]) / 2;
+		          }
+		        }
 		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = max;
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        _tree_set->insert(blossom);
 		        (*_blossom_data)[blossom].status = EVEN;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = INVALID;
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int ti = (*_node_index)[_graph.v(e)];
 		        if (_graph.u(e) != _graph.v(e)) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
+		    }
 		    /// \brief Initialize the algorithm with fractional matching
 		    ///
 		    /// This function initializes the algorithm with a fractional
 		    /// matching. This initialization is also called jumpstart heuristic.
 		    void fractionalInit() {
 		      createStructures();
 		      if (_fractional == 0) {
 		        _fractional = new FractionalMatching(_graph, _weight, false);
+		      }
 		      if (!_fractional->run()) {
 		        _unmatched = -1;
 		        return;
+		      }
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = 0;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value pot = _fractional->nodeValue(n);
 		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = pot;
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        (*_blossom_data)[blossom].status = MATCHED;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = _fractional->matching(n);
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      typename Graph::template NodeMap<bool> processed(_graph, false);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if (processed[n]) continue;
 		        processed[n] = true;
 		        if (_fractional->matching(n) == INVALID) continue;
 		        int num = 1;
 		        Node v = _graph.target(_fractional->matching(n));
 		        while (n != v) {
 		          processed[v] = true;
 		          v = _graph.target(_fractional->matching(v));
 		          ++num;
+		        }
 		        if (num % 2 == 1) {
 		          std::vector<int> subblossoms(num);
 		          subblossoms[--num] = _blossom_set->find(n);
 		          v = _graph.target(_fractional->matching(n));
 		          while (n != v) {
 		            subblossoms[--num] = _blossom_set->find(v);
 		            v = _graph.target(_fractional->matching(v));
+		          }
 		          int surface =
 		            _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		          (*_blossom_data)[surface].status = EVEN;
 		          (*_blossom_data)[surface].pred = INVALID;
 		          (*_blossom_data)[surface].next = INVALID;
 		          (*_blossom_data)[surface].pot = 0;
 		          (*_blossom_data)[surface].offset = 0;
 		          _tree_set->insert(surface);
 		          ++_unmatched;
+		        }
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int sb = _blossom_set->find(_graph.u(e));
 		        int ti = (*_node_index)[_graph.v(e)];
 		        int tb = _blossom_set->find(_graph.v(e));
 		        if ((*_blossom_data)[sb].status == EVEN &&
 		            (*_blossom_data)[tb].status == EVEN && sb != tb) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        int nb = _blossom_set->find(n);
 		        if ((*_blossom_data)[nb].status != MATCHED) continue;
 		        int ni = (*_node_index)[n];
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.target(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            int vt = _tree_set->find(vb);
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[ni].heap_index.find(vt);
 		            if (it != (*_node_data)[ni].heap_index.end()) {
 		              if ((*_node_data)[ni].heap[it->second] > rw) {
 		                (*_node_data)[ni].heap.replace(it->second, e);
 		                (*_node_data)[ni].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[ni].heap.push(e, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+		            }
+		          }
+		        }
 		        if (!(*_node_data)[ni].heap.empty()) {
 		          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		          _delta2->push(nb, _blossom_set->classPrio(nb));
+		        }
+		      }
+		    }
 		    /// \brief Start the algorithm
 		    ///
 		    /// This function starts the algorithm.
 		    ///
-		    /// \pre \ref init() must be called before using this function.
+		    /// \pre \ref init() or \ref fractionalInit() must be called before
 		    /// using this function.
 		    bool start() {
 		      enum OpType {
 		        D2, D3, D4
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		      if (_unmatched == -1) return false;
 		      while (_unmatched > 0) {
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
 		        Value d4 = !_delta4->empty() ?
 		          _delta4->prio() : std::numeric_limits<Value>::max();
 		        _delta_sum = d2; OpType ot = D2;
 		        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
 		        _delta_sum = d3; OpType ot = D3;
 		        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
 		        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
 		        if (_delta_sum == std::numeric_limits<Value>::max()) {
 		          return false;
+		        }
 		        switch (ot) {
 		        case D2:
+		          {
 		            int blossom = _delta2->top();
 		            Node n = _blossom_set->classTop(blossom);
 		            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
 		            extendOnArc(e);
+		          }
 		          break;
 		        case D3:
+		          {
 		            Edge e = _delta3->top();
 		            int left_blossom = _blossom_set->find(_graph.u(e));
 		            int right_blossom = _blossom_set->find(_graph.v(e));
 		            if (left_blossom == right_blossom) {
 		              _delta3->pop();
 		            } else {
 		              int left_tree = _tree_set->find(left_blossom);
 		              int right_tree = _tree_set->find(right_blossom);
 		              if (left_tree == right_tree) {
 		                shrinkOnEdge(e, left_tree);
 		              } else {
 		                augmentOnEdge(e);
-		                unmatched -= 2;
+		                _unmatched -= 2;
+		              }
+		            }
 		          } break;
 		        case D4:
 		          splitBlossom(_delta4->top());
 		          break;
+		        }
+		      }
 		      extractMatching();
 		      return true;
+		    }
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
 		    ///
 		    /// \note mwpm.run() is just a shortcut of the following code.
 		    /// \code
-		    ///   mwpm.init();
+		    ///   mwpm.fractionalInit();
 		    ///   mwpm.start();
 		    /// \endcode
 		    bool run() {
-		      init();
+		      fractionalInit();
 		      return start();
+		    }
 		    /// @}
 		    /// \name Primal Solution
 		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// perfect matching.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the weight of the matching.
 		    ///
 		    /// This function returns the weight of the found matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value matchingWeight() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          sum += _weight[(*_matching)[n]];
+		        }
+		      }
-		      return sum /= 2;
 		      return sum / 2;
+		    }
 		    /// \brief Return \c true if the given edge is in the matching.
 		    ///
 		    /// This function returns \c true if the given edge is in the found
 		    /// matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    bool matching(const Edge& edge) const {
 		      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
+		    }
 		    /// \brief Return the matching arc (or edge) incident to the given node.
 		    ///
 		    /// This function returns the matching arc (or edge) incident to the
 		    /// given node in the found matching or \c INVALID if the node is
 		    /// not covered by the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Return a const reference to the matching map.
 		    ///
 		    /// This function returns a const reference to a node map that stores
 		    /// the matching arc (or edge) incident to each node.
 		    const MatchingMap& matchingMap() const {
 		      return *_matching;
+		    }
 		    /// \brief Return the mate of the given node.
 		    ///
 		    /// This function returns the mate of the given node in the found
 		    /// matching or \c INVALID if the node is not covered by the matching.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Node mate(const Node& node) const {
 		      return _graph.target((*_matching)[node]);
+		    }
 		    /// @}
 		    /// \name Dual Solution
 		    /// Functions to get the dual solution.\n
 		    /// Either \ref run() or \ref start() function should be called before
 		    /// using them.
 		    /// @{
 		    /// \brief Return the value of the dual solution.
 		    ///
 		    /// This function returns the value of the dual solution.
 		    /// It should be equal to the primal value scaled by \ref dualScale
 		    /// "dual scale".
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value dualValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        sum += nodeValue(n);
+		      }
 		      for (int i = 0; i < blossomNum(); ++i) {
 		        sum += blossomValue(i) * (blossomSize(i) / 2);
+		      }
 		      return sum;
+		    }
 		    /// \brief Return the dual value (potential) of the given node.
 		    ///
 		    /// This function returns the dual value (potential) of the given node.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value nodeValue(const Node& n) const {
 		      return (*_node_potential)[n];
+		    }
 		    /// \brief Return the number of the blossoms in the basis.
 		    ///
 		    /// This function returns the number of the blossoms in the basis.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    /// \see BlossomIt
 		    int blossomNum() const {
 		      return _blossom_potential.size();
+		    }
 		    /// \brief Return the number of the nodes in the given blossom.
 		    ///
 		    /// This function returns the number of the nodes in the given blossom.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    /// \see BlossomIt
 		    int blossomSize(int k) const {
 		      return _blossom_potential[k].end - _blossom_potential[k].begin;
+		    }
 		    /// \brief Return the dual value (ptential) of the given blossom.
 		    ///
 		    /// This function returns the dual value (ptential) of the given blossom.
 		    ///
 		    /// \pre Either run() or start() must be called before using this function.
 		    Value blossomValue(int k) const {
 		      return _blossom_potential[k].value;
+		    }
 		    /// \brief Iterator for obtaining the nodes of a blossom.
 		    ///
 		    /// This class provides an iterator for obtaining the nodes of the
 		    /// given blossom. It lists a subset of the nodes.
 		    /// Before using this iterator, you must allocate a
 		    /// MaxWeightedPerfectMatching class and execute it.
 		    class BlossomIt {
 		    public:
 		      /// \brief Constructor.
 		      ///
 		      /// Constructor to get the nodes of the given variable.
 		      ///
 		      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
 		      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
 		      /// must be called before initializing this iterator.
 		      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
 		        : _algorithm(&algorithm)
+		      {
 		        _index = _algorithm->_blossom_potential[variable].begin;
 		        _last = _algorithm->_blossom_potential[variable].end;
+		      }
 		      /// \brief Conversion to \c Node.
 		      ///
 		      /// Conversion to \c Node.
 		      operator Node() const {
 		        return _algorithm->_blossom_node_list[_index];
+		      }
 		      /// \brief Increment operator.
 		      ///
 		      /// Increment operator.
 		      BlossomIt& operator++() {
 		        ++_index;
 		        return *this;
+		      }
 		      /// \brief Validity checking
 		      ///
 		      /// This function checks whether the iterator is invalid.
 		      bool operator==(Invalid) const { return _index == _last; }
 		      /// \brief Validity checking
 		      ///
 		      /// This function checks whether the iterator is valid.
 		      bool operator!=(Invalid) const { return _index != _last; }
 		    private:
 		      const MaxWeightedPerfectMatching* _algorithm;
 		      int _last;
 		      int _index;
 		    };
 		    /// @}
 		  };
 		} //END OF NAMESPACE LEMON
-		#endif //LEMON_MAX_MATCHING_H
+		#endif //LEMON_MATCHING_H

test/CMakeLists.txt

Show inline comments

 		INCLUDE_DIRECTORIES(
 		  ${PROJECT_SOURCE_DIR}
 		  ${PROJECT_BINARY_DIR}
+		)
 		LINK_DIRECTORIES(
 		  ${PROJECT_BINARY_DIR}/lemon
+		)
 		SET(TESTS
 		  adaptors_test
 		  bellman_ford_test
 		  bfs_test
 		  circulation_test
 		  connectivity_test
 		  counter_test
 		  dfs_test
 		  digraph_test
 		  dijkstra_test
 		  dim_test
 		  edge_set_test
 		  error_test
 		  euler_test
 		  fractional_matching_test
 		  gomory_hu_test
 		  graph_copy_test
 		  graph_test
 		  graph_utils_test
 		  hao_orlin_test
 		  heap_test
 		  kruskal_test
 		  maps_test
 		  matching_test
 		  min_cost_arborescence_test
 		  min_cost_flow_test
 		  min_mean_cycle_test
 		  path_test
 		  planarity_test
 		  preflow_test
 		  radix_sort_test
 		  random_test
 		  suurballe_test
 		  time_measure_test
 		  unionfind_test
+		)
 		IF(LEMON_HAVE_LP)
 		  ADD_EXECUTABLE(lp_test lp_test.cc)
 		  SET(LP_TEST_LIBS lemon)
 		  IF(LEMON_HAVE_GLPK)
 		    SET(LP_TEST_LIBS ${LP_TEST_LIBS} ${GLPK_LIBRARIES})
 		  ENDIF()
 		  IF(LEMON_HAVE_CPLEX)
 		    SET(LP_TEST_LIBS ${LP_TEST_LIBS} ${CPLEX_LIBRARIES})
 		  ENDIF()
 		  IF(LEMON_HAVE_CLP)
 		    SET(LP_TEST_LIBS ${LP_TEST_LIBS} ${COIN_CLP_LIBRARIES})
 		  ENDIF()
 		  TARGET_LINK_LIBRARIES(lp_test ${LP_TEST_LIBS})
 		  ADD_TEST(lp_test lp_test)
 		  IF(WIN32 AND LEMON_HAVE_GLPK)
 		    GET_TARGET_PROPERTY(TARGET_LOC lp_test LOCATION)
 		    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
 		    ADD_CUSTOM_COMMAND(TARGET lp_test POST_BUILD
 		      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/glpk.dll ${TARGET_PATH}
 		      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/libltdl3.dll ${TARGET_PATH}
 		      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/zlib1.dll ${TARGET_PATH}
+		    )
 		  ENDIF()
 		  IF(WIN32 AND LEMON_HAVE_CPLEX)
 		    GET_TARGET_PROPERTY(TARGET_LOC lp_test LOCATION)
 		    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
 		    ADD_CUSTOM_COMMAND(TARGET lp_test POST_BUILD
 		      COMMAND ${CMAKE_COMMAND} -E copy ${CPLEX_BIN_DIR}/cplex91.dll ${TARGET_PATH}
+		    )
 		  ENDIF()
 		ENDIF()
 		IF(LEMON_HAVE_MIP)
 		  ADD_EXECUTABLE(mip_test mip_test.cc)
 		  SET(MIP_TEST_LIBS lemon)
 		  IF(LEMON_HAVE_GLPK)
 		    SET(MIP_TEST_LIBS ${MIP_TEST_LIBS} ${GLPK_LIBRARIES})
 		  ENDIF()
 		  IF(LEMON_HAVE_CPLEX)
 		    SET(MIP_TEST_LIBS ${MIP_TEST_LIBS} ${CPLEX_LIBRARIES})
 		  ENDIF()
 		  IF(LEMON_HAVE_CBC)
 		    SET(MIP_TEST_LIBS ${MIP_TEST_LIBS} ${COIN_CBC_LIBRARIES})
 		  ENDIF()
 		  TARGET_LINK_LIBRARIES(mip_test ${MIP_TEST_LIBS})
 		  ADD_TEST(mip_test mip_test)
 		  IF(WIN32 AND LEMON_HAVE_GLPK)
 		    GET_TARGET_PROPERTY(TARGET_LOC mip_test LOCATION)
 		    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
 		    ADD_CUSTOM_COMMAND(TARGET mip_test POST_BUILD
 		      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/glpk.dll ${TARGET_PATH}
 		      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/libltdl3.dll ${TARGET_PATH}
 		      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/zlib1.dll ${TARGET_PATH}
+		    )
 		  ENDIF()
 		  IF(WIN32 AND LEMON_HAVE_CPLEX)
 		    GET_TARGET_PROPERTY(TARGET_LOC mip_test LOCATION)
 		    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
 		    ADD_CUSTOM_COMMAND(TARGET mip_test POST_BUILD
 		      COMMAND ${CMAKE_COMMAND} -E copy ${CPLEX_BIN_DIR}/cplex91.dll ${TARGET_PATH}
+		    )
 		  ENDIF()
 		ENDIF()
 		FOREACH(TEST_NAME ${TESTS})
 		  ADD_EXECUTABLE(${TEST_NAME} ${TEST_NAME}.cc)
 		  TARGET_LINK_LIBRARIES(${TEST_NAME} lemon)
 		  ADD_TEST(${TEST_NAME} ${TEST_NAME})
 		ENDFOREACH()

test/Makefile.am

Show inline comments

 		if USE_VALGRIND
 		TESTS_ENVIRONMENT=$(top_srcdir)/scripts/valgrind-wrapper.sh
 		endif
 		EXTRA_DIST += \
 			test/CMakeLists.txt
 		noinst_HEADERS += \
 			test/graph_test.h \
 			test/test_tools.h
 		check_PROGRAMS += \
 			test/adaptors_test \
 			test/bellman_ford_test \
 			test/bfs_test \
 			test/circulation_test \
 			test/connectivity_test \
 			test/counter_test \
 			test/dfs_test \
 			test/digraph_test \
 			test/dijkstra_test \
 			test/dim_test \
 			test/edge_set_test \
 			test/error_test \
 			test/euler_test \
 			test/fractional_matching_test \
 			test/gomory_hu_test \
 			test/graph_copy_test \
 			test/graph_test \
 			test/graph_utils_test \
 			test/hao_orlin_test \
 			test/heap_test \
 			test/kruskal_test \
 			test/maps_test \
 			test/matching_test \
 			test/min_cost_arborescence_test \
 			test/min_cost_flow_test \
 			test/min_mean_cycle_test \
 			test/path_test \
 			test/planarity_test \
 			test/preflow_test \
 			test/radix_sort_test \
 			test/random_test \
 			test/suurballe_test \
 			test/test_tools_fail \
 			test/test_tools_pass \
 			test/time_measure_test \
 			test/unionfind_test
 		test_test_tools_pass_DEPENDENCIES = demo
 		if HAVE_LP
 		check_PROGRAMS += test/lp_test
 		endif HAVE_LP
 		if HAVE_MIP
 		check_PROGRAMS += test/mip_test
 		endif HAVE_MIP
 		TESTS += $(check_PROGRAMS)
 		XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
 		test_adaptors_test_SOURCES = test/adaptors_test.cc
 		test_bellman_ford_test_SOURCES = test/bellman_ford_test.cc
 		test_bfs_test_SOURCES = test/bfs_test.cc
 		test_circulation_test_SOURCES = test/circulation_test.cc
 		test_counter_test_SOURCES = test/counter_test.cc
 		test_connectivity_test_SOURCES = test/connectivity_test.cc
 		test_dfs_test_SOURCES = test/dfs_test.cc
 		test_digraph_test_SOURCES = test/digraph_test.cc
 		test_dijkstra_test_SOURCES = test/dijkstra_test.cc
 		test_dim_test_SOURCES = test/dim_test.cc
 		test_edge_set_test_SOURCES = test/edge_set_test.cc
 		test_error_test_SOURCES = test/error_test.cc
 		test_euler_test_SOURCES = test/euler_test.cc
 		test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
 		test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
 		test_graph_copy_test_SOURCES = test/graph_copy_test.cc
 		test_graph_test_SOURCES = test/graph_test.cc
 		test_graph_utils_test_SOURCES = test/graph_utils_test.cc
 		test_heap_test_SOURCES = test/heap_test.cc
 		test_kruskal_test_SOURCES = test/kruskal_test.cc
 		test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
 		test_lp_test_SOURCES = test/lp_test.cc
 		test_maps_test_SOURCES = test/maps_test.cc
 		test_mip_test_SOURCES = test/mip_test.cc
 		test_matching_test_SOURCES = test/matching_test.cc
 		test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
 		test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
 		test_min_mean_cycle_test_SOURCES = test/min_mean_cycle_test.cc
 		test_path_test_SOURCES = test/path_test.cc
 		test_planarity_test_SOURCES = test/planarity_test.cc
 		test_preflow_test_SOURCES = test/preflow_test.cc
 		test_radix_sort_test_SOURCES = test/radix_sort_test.cc
 		test_suurballe_test_SOURCES = test/suurballe_test.cc
 		test_random_test_SOURCES = test/random_test.cc
 		test_test_tools_fail_SOURCES = test/test_tools_fail.cc
 		test_test_tools_pass_SOURCES = test/test_tools_pass.cc
 		test_time_measure_test_SOURCES = test/time_measure_test.cc
 		test_unionfind_test_SOURCES = test/unionfind_test.cc

test/matching_test.cc

Show inline comments

@@ -212,213 +212,237 @@
 		  const_mat_test.matching(e);
 		  const_mat_test.matching(n);
 		  const MaxWeightedPerfectMatching<Graph>::MatchingMap& mmap =
 		    const_mat_test.matchingMap();
 		  e = mmap[n];
 		  const_mat_test.mate(n);
 		  int k = 0;
 		  const_mat_test.dualValue();
 		  const_mat_test.nodeValue(n);
 		  const_mat_test.blossomNum();
 		  const_mat_test.blossomSize(k);
 		  const_mat_test.blossomValue(k);
+		}
 		void checkMatching(const SmartGraph& graph,
 		                   const MaxMatching<SmartGraph>& mm) {
 		  int num = 0;
 		  IntNodeMap comp_index(graph);
 		  UnionFind<IntNodeMap> comp(comp_index);
 		  int barrier_num = 0;
 		  for (NodeIt n(graph); n != INVALID; ++n) {
 		    check(mm.status(n) == MaxMatching<SmartGraph>::EVEN ||
 		          mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
 		    if (mm.status(n) == MaxMatching<SmartGraph>::ODD) {
 		      ++barrier_num;
 		    } else {
 		      comp.insert(n);
+		    }
+		  }
 		  for (EdgeIt e(graph); e != INVALID; ++e) {
 		    if (mm.matching(e)) {
 		      check(e == mm.matching(graph.u(e)), "Wrong matching");
 		      check(e == mm.matching(graph.v(e)), "Wrong matching");
 		      ++num;
+		    }
 		    check(mm.status(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
 		          mm.status(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
 		          "Wrong Gallai-Edmonds decomposition");
 		    check(mm.status(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
 		          mm.status(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
 		          "Wrong Gallai-Edmonds decomposition");
 		    if (mm.status(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
 		        mm.status(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
 		      comp.join(graph.u(e), graph.v(e));
+		    }
+		  }
 		  std::set<int> comp_root;
 		  int odd_comp_num = 0;
 		  for (NodeIt n(graph); n != INVALID; ++n) {
 		    if (mm.status(n) != MaxMatching<SmartGraph>::ODD) {
 		      int root = comp.find(n);
 		      if (comp_root.find(root) == comp_root.end()) {
 		        comp_root.insert(root);
 		        if (comp.size(n) % 2 == 1) {
 		          ++odd_comp_num;
+		        }
+		      }
+		    }
+		  }
 		  check(mm.matchingSize() == num, "Wrong matching");
 		  check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
 		         "Wrong matching");
 		  return;
+		}
 		void checkWeightedMatching(const SmartGraph& graph,
 		                   const SmartGraph::EdgeMap<int>& weight,
 		                   const MaxWeightedMatching<SmartGraph>& mwm) {
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    if (graph.u(e) == graph.v(e)) continue;
 		    int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
 		    for (int i = 0; i < mwm.blossomNum(); ++i) {
 		      bool s = false, t = false;
 		      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
 		           n != INVALID; ++n) {
 		        if (graph.u(e) == n) s = true;
 		        if (graph.v(e) == n) t = true;
+		      }
 		      if (s == true && t == true) {
 		        rw += mwm.blossomValue(i);
+		      }
+		    }
 		    rw -= weight[e] * mwm.dualScale;
 		    check(rw >= 0, "Negative reduced weight");
 		    check(rw == 0 || !mwm.matching(e),
 		          "Non-zero reduced weight on matching edge");
+		  }
 		  int pv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    if (mwm.matching(n) != INVALID) {
 		      check(mwm.nodeValue(n) >= 0, "Invalid node value");
 		      pv += weight[mwm.matching(n)];
 		      SmartGraph::Node o = graph.target(mwm.matching(n));
 		      check(mwm.mate(n) == o, "Invalid matching");
 		      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
 		            "Invalid matching");
 		    } else {
 		      check(mwm.mate(n) == INVALID, "Invalid matching");
 		      check(mwm.nodeValue(n) == 0, "Invalid matching");
+		    }
+		  }
 		  int dv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    dv += mwm.nodeValue(n);
+		  }
 		  for (int i = 0; i < mwm.blossomNum(); ++i) {
 		    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
 		    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
 		    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
+		  }
 		  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
 		  return;
+		}
 		void checkWeightedPerfectMatching(const SmartGraph& graph,
 		                          const SmartGraph::EdgeMap<int>& weight,
 		                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
 		  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
 		    if (graph.u(e) == graph.v(e)) continue;
 		    int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
 		    for (int i = 0; i < mwpm.blossomNum(); ++i) {
 		      bool s = false, t = false;
 		      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
 		           n != INVALID; ++n) {
 		        if (graph.u(e) == n) s = true;
 		        if (graph.v(e) == n) t = true;
+		      }
 		      if (s == true && t == true) {
 		        rw += mwpm.blossomValue(i);
+		      }
+		    }
 		    rw -= weight[e] * mwpm.dualScale;
 		    check(rw >= 0, "Negative reduced weight");
 		    check(rw == 0 || !mwpm.matching(e),
 		          "Non-zero reduced weight on matching edge");
+		  }
 		  int pv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    check(mwpm.matching(n) != INVALID, "Non perfect");
 		    pv += weight[mwpm.matching(n)];
 		    SmartGraph::Node o = graph.target(mwpm.matching(n));
 		    check(mwpm.mate(n) == o, "Invalid matching");
 		    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
 		          "Invalid matching");
+		  }
 		  int dv = 0;
 		  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
 		    dv += mwpm.nodeValue(n);
+		  }
 		  for (int i = 0; i < mwpm.blossomNum(); ++i) {
 		    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
 		    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
 		    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
+		  }
 		  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
 		  return;
+		}
 		int main() {
 		  for (int i = 0; i < lgfn; ++i) {
 		    SmartGraph graph;
 		    SmartGraph::EdgeMap<int> weight(graph);
 		    istringstream lgfs(lgf[i]);
 		    graphReader(graph, lgfs).
 		      edgeMap("weight", weight).run();
 		    bool perfect;
+		    {
 		    MaxMatching<SmartGraph> mm(graph);
 		    mm.run();
 		    checkMatching(graph, mm);
 		      perfect = 2 * mm.matchingSize() == countNodes(graph);
+		    }
+		    {
 		    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
 		    mwm.run();
 		    checkWeightedMatching(graph, weight, mwm);
+		    }
+		    {
 		      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
 		      mwm.init();
 		      mwm.start();
 		      checkWeightedMatching(graph, weight, mwm);
+		    }
+		    {
 		    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
-		    bool perfect = mwpm.run();
+		      bool result = mwpm.run();
 		    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
 		          "Perfect matching found");
+		      check(result == perfect, "Perfect matching found");
 		    if (perfect) {
 		      checkWeightedPerfectMatching(graph, weight, mwpm);
+		    }
+		  }
+		    {
 		      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
 		      mwpm.init();
 		      bool result = mwpm.start();
 		      check(result == perfect, "Perfect matching found");
 		      if (perfect) {
 		        checkWeightedPerfectMatching(graph, weight, mwpm);
+		      }
+		    }
+		  }
 		  return 0;
+		}

0 comments (0 inline)

RhodeCode

Login to your account