RhodeCode

      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

      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.

    /// @{

          (*_blossom_data)[tb].status = EVEN;

          matchedToEven(tb, tree);

          _tree_set->insert(tb, tree);

          (*_blossom_data)[tb].pred =

            (*_blossom_data)[tb].next =

            _graph.oppositeArc((*_blossom_data)[ub].next);

          next = (*_blossom_data)[ub].next;

        }

        (*_blossom_data)[subblossoms[ib]].status = ODD;

        matchedToOdd(subblossoms[ib]);

        _tree_set->insert(subblossoms[ib], tree);

        (*_blossom_data)[subblossoms[ib]].next = next;

        (*_blossom_data)[subblossoms[ib]].pred = pred;

      }

      _tree_set->erase(blossom);

    }

    void extractBlossom(int blossom, const Node& base, const Arc& matching) {

      if (_blossom_set->trivial(blossom)) {

        int bi = (*_node_index)[base];

        Value pot = (*_node_data)[bi].pot;

        (*_matching)[base] = matching;

        _blossom_node_list.push_back(base);

        (*_node_potential)[base] = pot;

      } else {

        Value pot = (*_blossom_data)[blossom].pot;

        int bn = _blossom_node_list.size();

        std::vector<int> subblossoms;

        _blossom_set->split(blossom, std::back_inserter(subblossoms));

        int b = _blossom_set->find(base);

        int ib = -1;

        for (int i = 0; i < int(subblossoms.size()); ++i) {

          if (subblossoms[i] == b) { ib = i; break; }

        }

        for (int i = 1; i < int(subblossoms.size()); i += 2) {

          int sb = subblossoms[(ib + i) % subblossoms.size()];

          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];

          Arc m = (*_blossom_data)[tb].next;

          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));

          extractBlossom(tb, _graph.source(m), m);

        }

        extractBlossom(subblossoms[ib], base, matching);

        int en = _blossom_node_list.size();

        _blossom_potential.push_back(BlossomVariable(bn, en, pot));

      }

    }

    void extractMatching() {

      std::vector<int> blossoms;

      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {

        blossoms.push_back(c);

      }

      for (int i = 0; i < int(blossoms.size()); ++i) {

        if ((*_blossom_data)[blossoms[i]].next != INVALID) {

          Value offset = (*_blossom_data)[blossoms[i]].offset;

          (*_blossom_data)[blossoms[i]].pot += 2 * offset;

          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);

               n != INVALID; ++n) {

            (*_node_data)[(*_node_index)[n]].pot -= offset;

          }

          Arc matching = (*_blossom_data)[blossoms[i]].next;

          Node base = _graph.source(matching);

          extractBlossom(blossoms[i], base, matching);

        } else {

          Node base = (*_blossom_data)[blossoms[i]].base;

          extractBlossom(blossoms[i], base, INVALID);

        }

      }

    }

  public:

    /// \brief Constructor

    ///

    /// Constructor.

    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)

      : _graph(graph), _weight(weight), _matching(0),

        _node_potential(0), _blossom_potential(), _blossom_node_list(),

        _node_num(0), _blossom_num(0),

        _blossom_index(0), _blossom_set(0), _blossom_data(0),

        _node_index(0), _node_heap_index(0), _node_data(0),

        _tree_set_index(0), _tree_set(0),

        _delta1_index(0), _delta1(0),

        _delta2_index(0), _delta2(0),

        _delta3_index(0), _delta3(0),

        _delta4_index(0), _delta4(0),

        _delta_sum(), _unmatched(0),

        _fractional(0)

    {}

    ~MaxWeightedMatching() {

      destroyStructures();

      if (_fractional) {

        delete _fractional;

      }

    }

    /// \name Execution Control

    /// The simplest way to execute the algorithm is to use the

    /// \ref run() member function.

    ///@{

    /// \brief Initialize the algorithm

    ///

    /// This function initializes the algorithm.

    void init() {

      createStructures();

      _blossom_node_list.clear();

      _blossom_potential.clear();

      for (ArcIt e(_graph); e != INVALID; ++e) {

        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;

      }

      for (NodeIt n(_graph); n != INVALID; ++n) {

        (*_delta1_index)[n] = _delta1->PRE_HEAP;

      }

      for (EdgeIt e(_graph); e != INVALID; ++e) {

        (*_delta3_index)[e] = _delta3->PRE_HEAP;

      }

      for (int i = 0; i < _blossom_num; ++i) {

        (*_delta2_index)[i] = _delta2->PRE_HEAP;

        (*_delta4_index)[i] = _delta4->PRE_HEAP;

      }

      _unmatched = _node_num;

      _delta1->clear();

      _delta2->clear();

      _delta3->clear();

      _delta4->clear();

      _blossom_set->clear();

      _tree_set->clear();

      int index = 0;

      for (NodeIt n(_graph); n != INVALID; ++n) {

        Value max = 0;

        for (OutArcIt e(_graph, n); e != INVALID; ++e) {

          if (_graph.target(e) == n) continue;

          if ((dualScale * _weight[e]) / 2 > max) {

            max = (dualScale * _weight[e]) / 2;

          }

        }

        (*_node_index)[n] = index;

        (*_node_data)[index].heap_index.clear();

        (*_node_data)[index].heap.clear();

        (*_node_data)[index].pot = max;

        _delta1->push(n, max);

        int blossom =

          _blossom_set->insert(n, std::numeric_limits<Value>::max());

        _tree_set->insert(blossom);

        (*_blossom_data)[blossom].status = EVEN;

        (*_blossom_data)[blossom].pred = INVALID;

        (*_blossom_data)[blossom].next = INVALID;

        (*_blossom_data)[blossom].pot = 0;

        (*_blossom_data)[blossom].offset = 0;

        ++index;

      }

      for (EdgeIt e(_graph); e != INVALID; ++e) {

        int si = (*_node_index)[_graph.u(e)];

        int ti = (*_node_index)[_graph.v(e)];

        if (_graph.u(e) != _graph.v(e)) {

          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -

                            dualScale * _weight[e]) / 2);

        }

      }

    }

    /// \brief Initialize the algorithm with fractional matching

    ///

    /// This function initializes the algorithm with a fractional

    /// matching. This initialization is also called jumpstart heuristic.

    void fractionalInit() {

      createStructures();

      if (_fractional == 0) {

        _fractional = new FractionalMatching(_graph, _weight, false);

      }

      _fractional->run();

      for (ArcIt e(_graph); e != INVALID; ++e) {

        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;

      }

      for (NodeIt n(_graph); n != INVALID; ++n) {

        (*_delta1_index)[n] = _delta1->PRE_HEAP;

      }

      for (EdgeIt e(_graph); e != INVALID; ++e) {

        (*_delta3_index)[e] = _delta3->PRE_HEAP;

      }

      for (int i = 0; i < _blossom_num; ++i) {

        (*_delta2_index)[i] = _delta2->PRE_HEAP;

        (*_delta4_index)[i] = _delta4->PRE_HEAP;

      }

      _unmatched = 0;

      int index = 0;

      for (NodeIt n(_graph); n != INVALID; ++n) {

        Value pot = _fractional->nodeValue(n);

        (*_node_index)[n] = index;

        (*_node_data)[index].pot = pot;

        int blossom =

          _blossom_set->insert(n, std::numeric_limits<Value>::max());

        (*_blossom_data)[blossom].status = MATCHED;

        (*_blossom_data)[blossom].pred = INVALID;

        (*_blossom_data)[blossom].next = _fractional->matching(n);

        if (_fractional->matching(n) == INVALID) {

          (*_blossom_data)[blossom].base = n;

        }

        (*_blossom_data)[blossom].pot = 0;

        (*_blossom_data)[blossom].offset = 0;

        ++index;

      }

      typename Graph::template NodeMap<bool> processed(_graph, false);

      for (NodeIt n(_graph); n != INVALID; ++n) {

        if (processed[n]) continue;

        processed[n] = true;

        if (_fractional->matching(n) == INVALID) continue;

        int num = 1;

        Node v = _graph.target(_fractional->matching(n));

        while (n != v) {

          processed[v] = true;

          v = _graph.target(_fractional->matching(v));

          ++num;

        }

        if (num % 2 == 1) {

          std::vector<int> subblossoms(num);

          subblossoms[--num] = _blossom_set->find(n);

          _delta1->push(n, _fractional->nodeValue(n));

          v = _graph.target(_fractional->matching(n));

          while (n != v) {

            subblossoms[--num] = _blossom_set->find(v);

            _delta1->push(v, _fractional->nodeValue(v));

            v = _graph.target(_fractional->matching(v));            

          }

          int surface = 

            _blossom_set->join(subblossoms.begin(), subblossoms.end());

          (*_blossom_data)[surface].status = EVEN;

          (*_blossom_data)[surface].pred = INVALID;

          (*_blossom_data)[surface].next = INVALID;

          (*_blossom_data)[surface].pot = 0;

          (*_blossom_data)[surface].offset = 0;

          _tree_set->insert(surface);

          ++_unmatched;

        }

      }

      for (EdgeIt e(_graph); e != INVALID; ++e) {

        int si = (*_node_index)[_graph.u(e)];

        int sb = _blossom_set->find(_graph.u(e));

        int ti = (*_node_index)[_graph.v(e)];

        int tb = _blossom_set->find(_graph.v(e));

        if ((*_blossom_data)[sb].status == EVEN &&

            (*_blossom_data)[tb].status == EVEN && sb != tb) {

          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -

                            dualScale * _weight[e]) / 2);

        }

      }

      for (NodeIt n(_graph); n != INVALID; ++n) {

        int nb = _blossom_set->find(n);

        if ((*_blossom_data)[nb].status != MATCHED) continue;

        int ni = (*_node_index)[n];

        for (OutArcIt e(_graph, n); e != INVALID; ++e) {

          Node v = _graph.target(e);

          int vb = _blossom_set->find(v);

          int vi = (*_node_index)[v];

          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

            dualScale * _weight[e];

          if ((*_blossom_data)[vb].status == EVEN) {

            int vt = _tree_set->find(vb);

            typename std::map<int, Arc>::iterator it =

              (*_node_data)[ni].heap_index.find(vt);

            if (it != (*_node_data)[ni].heap_index.end()) {

              if ((*_node_data)[ni].heap[it->second] > rw) {

                (*_node_data)[ni].heap.replace(it->second, e);

                (*_node_data)[ni].heap.decrease(e, rw);

                it->second = e;

              }

            } else {

              (*_node_data)[ni].heap.push(e, rw);

              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));

            }

          }

        }

        if (!(*_node_data)[ni].heap.empty()) {

          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());

          _delta2->push(nb, _blossom_set->classPrio(nb));

        }

      }

    }