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

       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