RhodeCode

              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;

              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 

        _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));

            }

          }

        }

        _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),

        _delta2_index(0), _delta2(0),

        _delta3_index(0), _delta3(0),

        _delta4_index(0), _delta4(0),

        _delta_sum(), _unmatched(0),

        _fractional(0)

    {}

    ~MaxWeightedPerfectMatching() {

      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 (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;

      _delta2->clear();

      _delta3->clear();

      _delta4->clear();

      _blossom_set->clear();

      _tree_set->clear();

      int index = 0;

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

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

      }

      if (!_fractional->run()) {

        _unmatched = -1;

        return;

      }

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

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

      }

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

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

      }

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

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

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

      }

      _unmatched = 0;

      int index = 0;

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

        Value pot = _fractional->nodeValue(n);

        (*_node_index)[n] = index;

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

        int blossom =

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

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

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

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

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

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

        ++index;

      }

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

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

        if (processed[n]) continue;

        processed[n] = true;

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

        int num = 1;

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

        while (n != v) {

          processed[v] = true;

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

          ++num;

        }

        if (num % 2 == 1) {

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

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

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

          while (n != v) {

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

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

          }

          int surface = 

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

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

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

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

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

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

          _tree_set->insert(surface);

          ++_unmatched;

        }

      }

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

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

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

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

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

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

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

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

                            dualScale * _weight[e]) / 2);

        }

      }

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

        int nb = _blossom_set->find(n);

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

        int ni = (*_node_index)[n];

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

          Node v = _graph.target(e);

          int vb = _blossom_set->find(v);

          int vi = (*_node_index)[v];

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

            dualScale * _weight[e];

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

            int vt = _tree_set->find(vb);

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

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

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

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

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

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

                it->second = e;

              }

            } else {

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

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

            }

          }

        }

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

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

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

        }