/* -*- 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 MaxWeightedMatching::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 = ↦

      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 is based on extensive use

  /// of 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.

  ///

  /// If the value type is integer, then the primal and the dual

  /// solutions are multiplied by

  /// \ref MaxWeightedMatching::primalScale "2" and

  /// \ref MaxWeightedMatching::dualScale "4" respectively.

  ///

  /// \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. It is equal to 2 or 1

    /// according to the value type.

    static const int primalScale =

      std::numeric_limits<Value>::is_integer ? 2 : 1;

    /// \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);