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 3246 insertions and 378 deletions:

lemon/fractional_matching.h

2130

test/fractional_matching_test.cc

525

doc/groups.dox

lemon/Makefile.am

lemon/matching.h

543

365

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

@@ -387,5 +387,5 @@
 		problem of maximum flow.
-		\ref Circulation is a preflow push-relabel algorithm implemented directly
 		\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.
@@ -523,4 +523,11 @@
 		  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

lemon/Makefile.am

Show inline comments

@@ -85,4 +85,5 @@
 			lemon/euler.h \
 			lemon/fib_heap.h \
 			lemon/fractional_matching.h \
 			lemon/full_graph.h \
 			lemon/glpk.h \

lemon/matching.h

Show inline comments

@@ -17,6 +17,6 @@
 		 */
 		#ifndef LEMON_MAX_MATCHING_H
 		#define LEMON_MAX_MATCHING_H
 		#ifndef LEMON_MATCHING_H
 		#define LEMON_MATCHING_H
 		#include <vector>
@@ -29,4 +29,5 @@
 		#include <lemon/bin_heap.h>
 		#include <lemon/maps.h>
 		#include <lemon/fractional_matching.h>
 		///\ingroup matching
@@ -42,5 +43,5 @@
 		  /// 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
 		  /// It can be started from an arbitrary initial matching
 		  /// (the default is the empty one).
 		  ///
@@ -70,9 +71,9 @@
 		    ///\brief Status constants for Gallai-Edmonds decomposition.
 		    ///
-		    ///These constants are used for indicating the Gallai-Edmonds
 		    ///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
 		    ///with status \c MATCHED (or \c C) induce a subgraph having a
 		    ///perfect matching.
 		    enum Status {
@@ -513,5 +514,5 @@
+		    }
-		    /// \brief Start Edmonds' algorithm with a heuristic improvement
 		    /// \brief Start Edmonds' algorithm with a heuristic improvement
 		    /// for dense graphs
 		    ///
@@ -535,6 +536,6 @@
 		    /// \brief Run Edmonds' algorithm
 		    ///
 		    /// This function runs Edmonds' algorithm. An additional heuristic of
 		    /// postponing shrinks is used for relatively dense graphs
 		    /// This function runs Edmonds' algorithm. An additional heuristic of
 		    /// postponing shrinks is used for relatively dense graphs
 		    /// (for which <tt>m>=2*n</tt> holds).
 		    void run() {
@@ -557,5 +558,5 @@
 		    /// \brief Return the size (cardinality) of the matching.
 		    ///
-		    /// This function returns the size (cardinality) of the current 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 {
@@ -571,5 +572,5 @@
 		    /// \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
 		    /// This function returns \c true if the given edge is in the current
 		    /// matching.
 		    bool matching(const Edge& edge) const {
@@ -580,5 +581,5 @@
 		    ///
 		    /// This function returns the matching arc (or edge) incident to the
-		    /// given node in the current matching or \c INVALID if the node is
 		    /// given node in the current matching or \c INVALID if the node is
 		    /// not covered by the matching.
 		    Arc matching(const Node& n) const {
@@ -596,5 +597,5 @@
 		    /// \brief Return the mate of the given node.
 		    ///
-		    /// This function returns the mate of the given node in the current
 		    /// 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 {
@@ -606,5 +607,5 @@
 		    /// \name Dual Solution
-		    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
 		    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
 		    /// decomposition.
@@ -649,6 +650,6 @@
 		  /// \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
 		  /// 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.
@@ -674,14 +675,14 @@
 		      \frac{\vert B \vert - 1}{2}z_B\f] */
 		  ///
-		  /// The algorithm can be executed with the run() function.
 		  /// 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.
+		  /// 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
 		  /// \tparam WM The type edge weight map. The default type is
 		  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 		#ifdef DOXYGEN
@@ -746,5 +747,5 @@
 		    enum Status {
-		      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
 		      EVEN = -1, MATCHED = 0, ODD = 1
 		    };
@@ -798,4 +799,8 @@
 		    Value _delta_sum;
 		    int _unmatched;
 		    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
 		    FractionalMatching *_fractional;
 		    void createStructures() {
@@ -845,7 +850,4 @@
 		    void destroyStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (_matching) {
 		        delete _matching;
@@ -923,8 +925,4 @@
 		              _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 =
@@ -950,200 +948,4 @@
 		                               (*_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) -
@@ -1158,21 +960,29 @@
+		    }
 		    void matchedToUnmatched(int blossom) {
 		    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];
 		        _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);
 		        (*_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];
@@ -1181,7 +991,72 @@
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP) {
-		              _delta3->push(e, rw);
+		          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) {
 		                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);
+		                }
+		              }
+		            }
+		          }
@@ -1190,9 +1065,38 @@
+		    }
-		    void unmatchedToMatched(int blossom) {
+		    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);
@@ -1203,52 +1107,42 @@
 		            dualScale * _weight[e];
 		          if (vb == blossom) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
-		              _delta3->erase(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 == EVEN) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		            int vt = _tree_set->find(vb);
-		            Arc r = _graph.oppositeArc(e);
+		          } else {
 		            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;
+		              (*_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)[ni].heap.push(r, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
 		              (*_node_data)[vi].heap.push(e, rw);
 		              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+		            }
 		            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);
 		            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);
+		                }
+		              }
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
+		          }
+		        }
+		      }
 		      (*_blossom_data)[blossom].offset = 0;
+		    }
@@ -1295,10 +1189,8 @@
 		      destroyTree(tree);
 		      (*_blossom_data)[blossom].status = UNMATCHED;
 		      (*_blossom_data)[blossom].base = node;
-		      matchedToUnmatched(blossom);
+		      (*_blossom_data)[blossom].next = INVALID;
+		    }
 		    void augmentOnEdge(const Edge& edge) {
@@ -1306,21 +1198,11 @@
 		      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);
+		      }
+		      int left_tree = _tree_set->find(left);
 		      alternatePath(left, left_tree);
 		      destroyTree(left_tree);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      (*_blossom_data)[left].next = _graph.direct(edge, true);
@@ -1328,4 +1210,19 @@
+		    }
 		    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));
@@ -1530,5 +1427,5 @@
 		          (*_blossom_data)[sb].pred = pred;
 		          (*_blossom_data)[sb].next =
-		                           _graph.oppositeArc((*_blossom_data)[tb].next);
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
 		          pred = (*_blossom_data)[ub].next;
@@ -1630,5 +1527,5 @@
 		      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;
@@ -1668,8 +1565,14 @@
 		        _delta4_index(0), _delta4(0),
-		        _delta_sum() {}
+		        _delta_sum(), _unmatched(0),
 		        _fractional(0)
 		    {}
 		    ~MaxWeightedMatching() {
 		      destroyStructures();
 		      if (_fractional) {
 		        delete _fractional;
+		      }
+		    }
@@ -1700,4 +1603,6 @@
+		      }
 		      _unmatched = _node_num;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
@@ -1734,9 +1639,147 @@
+		    }
 		    /// \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 {
@@ -1744,6 +1787,5 @@
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		      while (_unmatched > 0) {
 		        Value d1 = !_delta1->empty() ?
 		          _delta1->prio() : std::numeric_limits<Value>::max();
@@ -1758,10 +1800,9 @@
 		          _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:
@@ -1769,5 +1810,5 @@
 		            Node n = _delta1->top();
 		            unmatchNode(n);
-		            --unmatched;
+		            --_unmatched;
+		          }
 		          break;
@@ -1776,6 +1817,11 @@
 		            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;
@@ -1790,18 +1836,6 @@
 		              _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) {
@@ -1809,5 +1843,5 @@
 		              } else {
 		                augmentOnEdge(e);
-		                unmatched -= 2;
+		                _unmatched -= 2;
+		              }
+		            }
@@ -1827,9 +1861,9 @@
 		    /// \note mwm.run() is just a shortcut of the following code.
 		    /// \code
-		    ///   mwm.init();
+		    ///   mwm.fractionalInit();
 		    ///   mwm.start();
 		    /// \endcode
 		    void run() {
-		      init();
+		      fractionalInit();
 		      start();
+		    }
@@ -1838,5 +1872,5 @@
 		    /// \name Primal Solution
-		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// matching.\n
 		    /// Either \ref run() or \ref start() function should be called before
@@ -1857,5 +1891,5 @@
+		        }
+		      }
-		      return sum /= 2;
 		      return sum / 2;
+		    }
@@ -1877,5 +1911,5 @@
 		    /// \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
 		    /// This function returns \c true if the given edge is in the found
 		    /// matching.
 		    ///
@@ -1888,5 +1922,5 @@
 		    ///
 		    /// This function returns the matching arc (or edge) incident to the
-		    /// given node in the found matching or \c INVALID if the node is
 		    /// given node in the found matching or \c INVALID if the node is
 		    /// not covered by the matching.
 		    ///
@@ -1906,5 +1940,5 @@
 		    /// \brief Return the mate of the given node.
 		    ///
-		    /// This function returns the mate of the given node in the found
 		    /// 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.
 		    ///
@@ -1926,6 +1960,6 @@
 		    /// \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
 		    /// This function returns the value of the dual solution.
 		    /// It should be equal to the primal value scaled by \ref dualScale
 		    /// "dual scale".
 		    ///
@@ -1982,7 +2016,7 @@
 		    /// \brief Iterator for obtaining the nodes of a blossom.
 		    ///
-		    /// This class provides an iterator for obtaining the nodes of the
 		    /// 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
 		    /// Before using this iterator, you must allocate a
 		    /// MaxWeightedMatching class and execute it.
 		    class BlossomIt {
@@ -1993,6 +2027,6 @@
 		      /// Constructor to get the nodes of the given variable.
 		      ///
 		      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
 		      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
 		      /// \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)
@@ -2047,6 +2081,6 @@
 		  /// \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
 		  /// 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.
@@ -2071,14 +2105,14 @@
 		      \frac{\vert B \vert - 1}{2}z_B\f] */
 		  ///
-		  /// The algorithm can be executed with the run() function.
 		  /// 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.
+		  /// 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
 		  /// \tparam WM The type edge weight map. The default type is
 		  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 		#ifdef DOXYGEN
@@ -2191,4 +2225,9 @@
 		    Value _delta_sum;
 		    int _unmatched;
 		    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
 		    FractionalMatching;
 		    FractionalMatching *_fractional;
 		    void createStructures() {
@@ -2234,7 +2273,4 @@
 		    void destroyStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (_matching) {
 		        delete _matching;
@@ -2909,8 +2945,14 @@
 		        _delta4_index(0), _delta4(0),
-		        _delta_sum() {}
+		        _delta_sum(), _unmatched(0),
 		        _fractional(0)
 		    {}
 		    ~MaxWeightedPerfectMatching() {
 		      destroyStructures();
 		      if (_fractional) {
 		        delete _fractional;
+		      }
+		    }
@@ -2938,4 +2980,6 @@
+		      }
 		      _unmatched = _node_num;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
@@ -2971,9 +3015,142 @@
+		    }
 		    /// \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 {
@@ -2981,6 +3158,7 @@
 		      };
 		      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();
@@ -2992,6 +3170,6 @@
 		          _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; }
@@ -3026,5 +3204,5 @@
 		              } else {
 		                augmentOnEdge(e);
-		                unmatched -= 2;
+		                _unmatched -= 2;
+		              }
+		            }
@@ -3045,9 +3223,9 @@
 		    /// \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();
+		    }
@@ -3056,5 +3234,5 @@
 		    /// \name Primal Solution
-		    /// Functions to get the primal solution, i.e. the maximum weighted
 		    /// 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
@@ -3075,10 +3253,10 @@
+		        }
+		      }
-		      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
 		    /// This function returns \c true if the given edge is in the found
 		    /// matching.
 		    ///
@@ -3091,5 +3269,5 @@
 		    ///
 		    /// This function returns the matching arc (or edge) incident to the
-		    /// given node in the found matching or \c INVALID if the node is
 		    /// given node in the found matching or \c INVALID if the node is
 		    /// not covered by the matching.
 		    ///
@@ -3109,5 +3287,5 @@
 		    /// \brief Return the mate of the given node.
 		    ///
-		    /// This function returns the mate of the given node in the found
 		    /// 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.
 		    ///
@@ -3128,6 +3306,6 @@
 		    /// \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
 		    /// This function returns the value of the dual solution.
 		    /// It should be equal to the primal value scaled by \ref dualScale
 		    /// "dual scale".
 		    ///
@@ -3184,7 +3362,7 @@
 		    /// \brief Iterator for obtaining the nodes of a blossom.
 		    ///
-		    /// This class provides an iterator for obtaining the nodes of the
 		    /// 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
 		    /// Before using this iterator, you must allocate a
 		    /// MaxWeightedPerfectMatching class and execute it.
 		    class BlossomIt {
@@ -3195,6 +3373,6 @@
 		      /// Constructor to get the nodes of the given variable.
 		      ///
 		      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
 		      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
 		      /// \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)
@@ -3242,3 +3420,3 @@
 		} //END OF NAMESPACE LEMON
-		#endif //LEMON_MAX_MATCHING_H
+		#endif //LEMON_MATCHING_H

test/CMakeLists.txt

Show inline comments

@@ -22,4 +22,5 @@
 		  error_test
 		  euler_test
 		  fractional_matching_test
 		  gomory_hu_test
 		  graph_copy_test

test/Makefile.am

Show inline comments

@@ -24,4 +24,5 @@
 			test/error_test \
 			test/euler_test \
 			test/fractional_matching_test \
 			test/gomory_hu_test \
 			test/graph_copy_test \
@@ -72,4 +73,5 @@
 		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/matching_test.cc

Show inline comments

@@ -402,20 +402,44 @@
 		      edgeMap("weight", weight).run();
 		    MaxMatching<SmartGraph> mm(graph);
 		    mm.run();
-		    checkMatching(graph, mm);
+		    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.run();
 		      checkWeightedMatching(graph, weight, mwm);
+		    }
 		    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
 		    bool perfect = mwpm.run();
+		    {
 		      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
 		      mwm.init();
 		      mwm.start();
 		      checkWeightedMatching(graph, weight, mwm);
+		    }
 		    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
 		          "Perfect matching found");
+		    {
 		      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
 		      bool result = mwpm.run();
 		      check(result == perfect, "Perfect matching found");
 		      if (perfect) {
 		        checkWeightedPerfectMatching(graph, weight, mwpm);
+		      }
+		    }
 		    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);
+		      }
+		    }
+		  }

0 comments (0 inline)

RhodeCode

Login to your account