RhodeCode

      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.

  ///

  ///

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

    ///

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

    }

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

    }