Ticket #168: 3180fd18651b.patch

doc/groups.dox

@@ -523 +523 @@
   for calculating maximum cardinality matching in bipartite graphs.
 - \ref PrBipartiteMatching Push-relabel algorithm
   for calculating maximum cardinality matching in bipartite graphs.
-- \ref MaxWeightedBipartiteMatching
-  Successive shortest path algorithm for calculating maximum weighted
-  matching and maximum weighted bipartite matching in bipartite graphs.
+- \ref MaxWeightedBpMatching Multiple search tree augmenting path algorithm for
+  calculating maximum weighted matching in bipartite graphs.
 - \ref MinCostMaxBipartiteMatching
   Successive shortest path algorithm for calculating minimum cost maximum
   matching in bipartite graphs.

new file lemon/bpmatching.h

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2018
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_BPMATCHING_H
+#define LEMON_BPMATCHING_H
+
+#include <limits>
+
+#include <lemon/core.h>
+#include <lemon/unionfind.h>
+#include <lemon/bin_heap.h>
+#include <lemon/maps.h>
+
+///\ingroup matching
+///\file
+///\brief Maximum matching algorithms in bipartite graphs.
+
+namespace lemon {
+
+  /// \ingroup matching
+  ///
+  /// \brief Weighted matching in bipartite graphs
+  ///
+  /// This class provides an efficient implementation of multiple search tree
+  /// augmenting path matching algorithm. The implementation is based on
+  /// extensive use of priority queues and provides \f$O(nm\log n)\f$ time
+  /// complexity.
+  ///
+  /// The maximum weighted matching problem is to find a subset of the
+  /// edges in a bipartite graph with maximum overall weight for which
+  /// each node has at most one incident edge.
+  /// 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 algorithm calculates an optimal 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]
+  /// \tparam BPGR The bipartite graph type the algorithm runs on.
+  /// \tparam WM The type edge weight map. The default type is
+  /// \ref concepts::BpGraph::EdgeMap "BPGR::EdgeMap<int>".
+#ifdef DOXYGEN
+  template <typename BPGR, typename WM>
+#else
+  template <typename BPGR,
+            typename WM = typename BPGR::template EdgeMap<int> >
+#endif
+  class MaxWeightedBpMatching1 {
+  public:
+
+    /// The graph type of the algorithm
+    typedef BPGR BpGraph;
+    /// 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 BpGraph::template NodeMap<typename BpGraph::Arc>
+    MatchingMap;
+
+    /// \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_BPGRAPH_TYPEDEFS(BpGraph);
+
+    typedef typename BpGraph::template NodeMap<Value> NodePotential;
+
+    const BpGraph& _bpgraph;
+    const WeightMap& _weight;
+
+    MatchingMap* _matching;
+
+    NodePotential* _node_potential;
+
+    int _node_num;
+
+    enum Status {
+      EVEN = -1, MATCHED = 0, ODD = 1
+    };
+
+    typedef typename BpGraph::template NodeMap<Status> StatusMap;
+    StatusMap* _status;
+
+    typedef typename BpGraph::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;
+    int _unmatched;
+
+    void createStructures() {
+      _node_num = countNodes(_bpgraph);
+
+      if (!_matching) {
+        _matching = new MatchingMap(_bpgraph);
+      }
+
+      if (!_node_potential) {
+        _node_potential = new NodePotential(_bpgraph);
+      }
+
+      if (!_status) {
+        _status = new StatusMap(_bpgraph);
+      }
+
+      if (!_pred) {
+        _pred = new PredMap(_bpgraph);
+      }
+
+      if (!_tree_set) {
+        _tree_set_index = new IntNodeMap(_bpgraph);
+        _tree_set = new TreeSet(*_tree_set_index);
+      }
+
+      if (!_delta1) {
+        _delta1_index = new IntNodeMap(_bpgraph);
+        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
+      }
+
+      if (!_delta2) {
+        _delta2_index = new IntNodeMap(_bpgraph);
+        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
+      }
+
+      if (!_delta3) {
+        _delta3_index = new IntEdgeMap(_bpgraph);
+        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+      }
+    }
+
+    void destroyStructures() {
+      if (_matching) {
+        delete _matching;
+      }
+      if (_node_potential) {
+        delete _node_potential;
+      }
+      if (_status) {
+        delete _status;
+      }
+      if (_pred) {
+        delete _pred;
+      }
+      if (_tree_set) {
+        delete _tree_set_index;
+        delete _tree_set;
+      }
+      if (_delta2) {
+        delete _delta2_index;
+        delete _delta2;
+      }
+      if (_delta3) {
+        delete _delta3_index;
+        delete _delta3;
+      }
+    }
+
+    void matchedToEven(Node node, int tree) {
+      _tree_set->insert(node, tree);
+      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
+      _delta1->push(node, (*_node_potential)[node]);
+
+      if (_delta2->state(node) == _delta2->IN_HEAP) {
+        _delta2->erase(node);
+      }
+
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
+          dualScale * _weight[a];
+        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(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
+          dualScale * _weight[a];
+
+        if ((*_status)[v] == EVEN) {
+          _delta3->erase(a);
+          if (minrw > rw) {
+            min = _bpgraph.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(_bpgraph, v); aa != INVALID; ++aa) {
+              Node va = _bpgraph.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(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        if ((*_status)[v] != EVEN) continue;
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
+          dualScale * _weight[a];
+
+        if (minrw > rw) {
+          min = _bpgraph.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 = _bpgraph.target(prev);
+        even = _bpgraph.target((*_pred)[odd]);
+        _matching->set(odd, (*_pred)[odd]);
+        _status->set(odd, MATCHED);
+        oddToMatched(odd);
+
+        prev = (*_matching)[even];
+        _status->set(even, MATCHED);
+        _matching->set(even, _bpgraph.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 = _bpgraph.u(edge);
+      int left_tree = _tree_set->find(left);
+
+      alternatePath(left, left_tree);
+      destroyTree(left_tree);
+      _matching->set(left, _bpgraph.direct(edge, true));
+
+      Node right = _bpgraph.v(edge);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right, right_tree);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.direct(edge, false));
+    }
+
+    void augmentOnArc(const Arc& arc) {
+      Node left = _bpgraph.source(arc);
+      _status->set(left, MATCHED);
+      _matching->set(left, arc);
+      _pred->set(left, arc);
+
+      Node right = _bpgraph.target(arc);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right, right_tree);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.oppositeArc(arc));
+    }
+
+    void extendOnArc(const Arc& arc) {
+      Node base = _bpgraph.target(arc);
+      int tree = _tree_set->find(base);
+
+      Node odd = _bpgraph.source(arc);
+      _tree_set->insert(odd, tree);
+      _status->set(odd, ODD);
+      matchedToOdd(odd, tree);
+      _pred->set(odd, arc);
+
+      Node even = _bpgraph.target((*_matching)[odd]);
+      _tree_set->insert(even, tree);
+      _status->set(even, EVEN);
+      matchedToEven(even, tree);
+    }
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor.
+    MaxWeightedBpMatching1(const BpGraph& bpgraph, const WeightMap& weight)
+      : _bpgraph(bpgraph), _weight(weight), _matching(0),
+        _node_potential(0), _node_num(0),
+        _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(), _unmatched(0)
+    {}
+
+    ~MaxWeightedBpMatching1() {
+      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(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = _node_num;
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) / 2 > max) {
+            max = (dualScale * _weight[a]) / 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(_bpgraph); e != INVALID; ++e) {
+        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
+                          (*_node_potential)[_bpgraph.v(e)] -
+                          dualScale * _weight[e]) / 2);
+      }
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void redRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) > max) {
+            max = dualScale * _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+
+        _tree_set->insert(n);
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+
+        ++_unmatched;
+      }
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, MATCHED);
+
+        Arc min = INVALID;
+        Value minrw = std::numeric_limits<Value>::max();
+        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.source(a);
+          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
+                     dualScale * _weight[a];
+
+          if (minrw > rw) {
+            min = _bpgraph.oppositeArc(a);
+            minrw = rw;
+          }
+        }
+        if (min != INVALID) {
+          _pred->set(n, min);
+          _delta2->push(n, minrw);
+        } else {
+          _pred->set(n, INVALID);
+        }
+      }
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void blueRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) > max) {
+            max = dualScale * _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+
+        _tree_set->insert(n);
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+
+        ++_unmatched;
+      }
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, MATCHED);
+
+        Arc min = INVALID;
+        Value minrw = std::numeric_limits<Value>::max();
+        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.source(a);
+          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
+                     dualScale * _weight[a];
+
+          if (minrw > rw) {
+            min = _bpgraph.oppositeArc(a);
+            minrw = rw;
+          }
+        }
+        if (min != INVALID) {
+          _pred->set(n, min);
+          _delta2->push(n, minrw);
+        } else {
+          _pred->set(n, INVALID);
+        }
+      }
+    }
+
+    /// \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
+      };
+
+      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 {
+              extendOnArc(a);
+            }
+          }
+          break;
+        case D3:
+          {
+            Edge e = _delta3->top();
+            augmentOnEdge(e);
+            _unmatched -= 2;
+          }
+          break;
+        }
+      }
+    }
+
+    /// \brief Run the algorithm.
+    ///
+    /// This method runs the \c %MaxWeightedBpMatching1 algorithm.
+    ///
+    /// \note mwbpm.run() is just a shortcut of the following code.
+    /// \code
+    ///   mwbpm.init();
+    ///   mwbpm.start();
+    /// \endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Primal Solution
+    /// Functions to get the primal solution, i.e. the maximum weighted
+    /// bipartite 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value matchingWeight() const {
+      Value sum = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          sum += _weight[(*_matching)[n]];
+        }
+      }
+      return sum / 2;
+    }
+
+    /// \brief Return the size (cardinality) of the matching.
+    ///
+    /// This function returns the size (cardinality) of the found matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    int matchingSize() const {
+      int num = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++num;
+        }
+      }
+      return num / 2;
+    }
+
+    /// \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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpgraph.u(edge)];
+    }
+
+    /// \brief Return the matching arc (or edge) incident to the given node.
+    ///
+    /// This function returns the 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Arc matching(const Node& node) const {
+      return (*_matching)[node];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the found
+    /// matching or \c INVALID if the node is not covered by the matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Node mate(const Node& node) const {
+      return (*_matching)[node] != INVALID ?
+        _bpgraph.target((*_matching)[node]) : INVALID;
+    }
+
+    /// @}
+
+    /// \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(_bpgraph); n != INVALID; ++n) {
+        sum += nodeValue(n);
+      }
+      return sum;
+    }
+
+    /// \brief Return the dual value (potential) of the given node.
+    ///
+    /// This function returns the dual value (potential) of the given node.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// @}
+
+  };
+
+#ifdef DOXYGEN
+  template <typename BPGR, typename WM>
+#else
+  template <typename BPGR,
+            typename WM = typename BPGR::template EdgeMap<int> >
+#endif
+  class MaxWeightedBpMatching2 {
+  public:
+
+    /// The graph type of the algorithm
+    typedef BPGR BpGraph;
+    /// 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 BpGraph::template NodeMap<typename BpGraph::Arc>
+    MatchingMap;
+
+    /// \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 = 1;
+
+  private:
+
+    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
+
+    typedef typename BpGraph::template NodeMap<Value> NodePotential;
+
+    const BpGraph& _bpgraph;
+    const WeightMap& _weight;
+
+    MatchingMap* _matching;
+
+    NodePotential* _node_potential;
+
+    int _node_num;
+
+    typedef typename BpGraph::template NodeMap<Arc> PredMap;
+    PredMap* _pred;
+
+    typedef ExtendFindEnum<IntNodeMap> TreeSet;
+    IntNodeMap *_tree_set_index;
+    TreeSet *_tree_set;
+
+    IntNodeMap *_heap_index;
+    BinHeap<Value, IntNodeMap> *_heap;
+
+    Value _delta_sum;
+    int _unmatched;
+
+    void createStructures() {
+      _node_num = countNodes(_bpgraph);
+
+      if (!_matching) {
+        _matching = new MatchingMap(_bpgraph);
+      }
+
+      if (!_node_potential) {
+        _node_potential = new NodePotential(_bpgraph);
+      }
+
+      if (!_pred) {
+        _pred = new PredMap(_bpgraph);
+      }
+
+      if (!_tree_set) {
+        _tree_set_index = new IntNodeMap(_bpgraph);
+        _tree_set = new TreeSet(*_tree_set_index);
+      }
+
+      if (!_heap) {
+        _heap_index = new IntNodeMap(_bpgraph);
+        _heap = new BinHeap<Value, IntNodeMap>(*_heap_index);
+      }
+
+    }
+
+    void destroyStructures() {
+      if (_matching) {
+        delete _matching;
+      }
+      if (_node_potential) {
+        delete _node_potential;
+      }
+      if (_pred) {
+        delete _pred;
+      }
+      if (_tree_set) {
+        delete _tree_set_index;
+        delete _tree_set;
+      }
+      if (_heap) {
+        delete _heap_index;
+        delete _heap;
+      }
+    }
+
+    void matchedToEven(Node node, int tree) {
+      _tree_set->insert(node, tree);
+      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
+      _heap->push(node, (*_node_potential)[node]);
+
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] - _weight[a];
+        if (_heap->state(v) == _heap->PRE_HEAP) {
+          _pred->set(v, a);
+          _heap->push(v, rw);
+        } else if (_heap->state(v) == _heap->IN_HEAP && (*_heap)[v] > rw) {
+          _pred->set(v, a);
+          _heap->decrease(v, rw);
+        }
+      }
+    }
+
+    void matchedToOdd(Node node, int tree) {
+      _tree_set->insert(node, tree);
+      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
+
+      if (_heap->state(node) == _heap->IN_HEAP) {
+        _heap->erase(node);
+      } else {
+        _heap_index->set(node, _heap->POST_HEAP);
+      }
+    }
+
+    void evenToMatched(Node node, int tree) {
+      _heap->erase(node);
+      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        if (_heap->state(v) == _heap->IN_HEAP && (*_pred)[v] == a) {
+          Arc mina = INVALID;
+          Value minrwa = std::numeric_limits<Value>::max();
+          for (OutArcIt aa(_bpgraph, v); aa != INVALID; ++aa) {
+            Node va = _bpgraph.target(aa);
+            if (_heap->state(va) == _heap->POST_HEAP ||
+                _tree_set->find(va) == tree) continue;
+            Value rwa = (*_node_potential)[v] + (*_node_potential)[va] - _weight[aa];
+            if (minrwa > rwa) {
+              minrwa = rwa;
+              mina = aa;
+            }
+          }
+          if (mina != INVALID) {
+            _pred->set(v, mina);
+            _heap->increase(v, minrwa);
+          } else {
+            _pred->set(v, INVALID);
+            _heap->erase(v);
+            _heap_index->set(v, _heap->PRE_HEAP);
+          }
+        }
+      }
+    }
+
+    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(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        if (_heap->state(v) == _heap->POST_HEAP) continue;
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] - _weight[a];
+
+        if (minrw > rw) {
+          min = _bpgraph.oppositeArc(a);
+          minrw = rw;
+        }
+      }
+      if (min != INVALID) {
+        _pred->set(node, min);
+        _heap->push(node, minrw);
+      } else {
+        _pred->set(node, INVALID);
+        _heap_index->set(node, _heap->PRE_HEAP);
+      }
+    }
+
+    void alternatePath(Node even, int tree) {
+      Node odd;
+
+      evenToMatched(even, tree);
+
+      Arc prev = (*_matching)[even];
+      while (prev != INVALID) {
+        odd = _bpgraph.target(prev);
+        even = _bpgraph.target((*_pred)[odd]);
+        _matching->set(odd, (*_pred)[odd]);
+        oddToMatched(odd);
+
+        prev = (*_matching)[even];
+        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
+        evenToMatched(even, tree);
+      }
+    }
+
+    void destroyTree(int tree) {
+      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
+        if (_bpgraph.red(n)) {
+          if (_heap->state(n) != _heap->POST_HEAP) {
+            evenToMatched(n, tree);
+          }
+        } else {
+          if (_heap->state(n) == _heap->POST_HEAP) {
+            oddToMatched(n);
+          }
+        }
+      }
+      _tree_set->eraseClass(tree);
+    }
+
+    void unmatchNode(Node node) {
+      int tree = _tree_set->find(node);
+
+      alternatePath(node, tree);
+      destroyTree(tree);
+
+      _matching->set(node, INVALID);
+    }
+
+    void augmentOnArc(Arc arc) {
+      Node left = _bpgraph.source(arc);
+      _matching->set(left, arc);
+
+      Node right = _bpgraph.target(arc);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right, right_tree);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.oppositeArc(arc));
+    }
+
+    void extendOnArc(Arc arc) {
+      Node base = _bpgraph.target(arc);
+      int tree = _tree_set->find(base);
+
+      Node odd = _bpgraph.source(arc);
+      _tree_set->insert(odd, tree);
+      matchedToOdd(odd, tree);
+      _pred->set(odd, arc);
+
+      Node even = _bpgraph.target((*_matching)[odd]);
+      _tree_set->insert(even, tree);
+      matchedToEven(even, tree);
+    }
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor.
+    MaxWeightedBpMatching2(const BpGraph& bpgraph, const WeightMap& weight)
+      : _bpgraph(bpgraph), _weight(weight), _matching(0),
+        _node_potential(0), _node_num(0), _pred(0),
+        _tree_set_index(0), _tree_set(0),
+
+        _heap_index(0), _heap(0),
+
+        _delta_sum(), _unmatched(0)
+    {}
+
+    ~MaxWeightedBpMatching2() {
+      destroyStructures();
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void redRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        _heap_index->set(n, _heap->PRE_HEAP);
+      }
+
+      _heap->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if (_weight[a] > max) {
+            max = _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _heap->push(n, max);
+
+        _tree_set->insert(n);
+
+        ++_unmatched;
+      }
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _node_potential->set(n, 0);
+
+        Arc min = INVALID;
+        Value minrw = std::numeric_limits<Value>::max();
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.target(a);
+          Value rw = (*_node_potential)[n] + (*_node_potential)[v] - _weight[a];
+
+          if (minrw > rw) {
+            min = a;
+            minrw = rw;
+          }
+        }
+        if (min != INVALID) {
+          _pred->set(n, min);
+          _heap->push(n, minrw);
+        } else {
+          _pred->set(n, INVALID);
+        }
+      }
+    }
+
+    /// \brief Start the algorithm
+    ///
+    /// This function starts the algorithm.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (_unmatched > 0) {
+        _delta_sum = _heap->prio();
+        Node n = _heap->top();
+        if (_bpgraph.red(n)) {
+          unmatchNode(n);
+          --_unmatched;
+        } else if ((*_matching)[n] == INVALID) {
+          augmentOnArc((*_pred)[n]);
+          --_unmatched;
+        } else {
+          extendOnArc((*_pred)[n]);
+        }
+      }
+    }
+
+    /// \brief Run the algorithm.
+    ///
+    /// This method runs the \c %MaxWeightedBpMatching algorithm.
+    ///
+    /// \note mwbpm.run() is just a shortcut of the following code.
+    /// \code
+    ///   mwbpm.init();
+    ///   mwbpm.start();
+    /// \endcode
+    void run() {
+      redRootInit();
+      start();
+    }
+
+    /// @}
+
+    /// \name Primal Solution
+    /// Functions to get the primal solution, i.e. the maximum weighted
+    /// bipartite 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value matchingWeight() const {
+      Value sum = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          sum += _weight[(*_matching)[n]];
+        }
+      }
+      return sum / 2;
+    }
+
+    /// \brief Return the size (cardinality) of the matching.
+    ///
+    /// This function returns the size (cardinality) of the found matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    int matchingSize() const {
+      int num = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++num;
+        }
+      }
+      return num / 2;
+    }
+
+    /// \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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpgraph.u(edge)];
+    }
+
+    /// \brief Return the matching arc (or edge) incident to the given node.
+    ///
+    /// This function returns the 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Arc matching(const Node& node) const {
+      return (*_matching)[node];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the found
+    /// matching or \c INVALID if the node is not covered by the matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Node mate(const Node& node) const {
+      return (*_matching)[node] != INVALID ?
+        _bpgraph.target((*_matching)[node]) : INVALID;
+    }
+
+    /// @}
+
+    /// \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(_bpgraph); n != INVALID; ++n) {
+        sum += nodeValue(n);
+      }
+      return sum;
+    }
+
+    /// \brief Return the dual value (potential) of the given node.
+    ///
+    /// This function returns the dual value (potential) of the given node.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// @}
+
+  };
+
+#ifdef DOXYGEN
+  template <typename BPGR, typename WM>
+#else
+  template <typename BPGR,
+            typename WM = typename BPGR::template EdgeMap<int> >
+#endif
+  class MaxWeightedBpMatching3 {
+  public:
+
+    /// The graph type of the algorithm
+    typedef BPGR BpGraph;
+    /// 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 BpGraph::template NodeMap<typename BpGraph::Arc>
+    MatchingMap;
+
+    /// \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_BPGRAPH_TYPEDEFS(BpGraph);
+
+    typedef typename BpGraph::template NodeMap<Value> NodePotential;
+
+    const BpGraph& _bpgraph;
+    const WeightMap& _weight;
+
+    MatchingMap* _matching;
+
+    NodePotential* _node_potential;
+
+    int _node_num;
+
+    enum Status {
+      EVEN = -1, MATCHED = 0, ODD = 1
+    };
+
+    typedef typename BpGraph::template NodeMap<Status> StatusMap;
+    StatusMap* _status;
+
+    typedef typename BpGraph::template NodeMap<Arc> PredMap;
+    PredMap* _pred;
+
+    IntArcMap* _node_heap_index;
+    typedef typename BpGraph::template NodeMap<BinHeap<Value, IntArcMap>*> NodeHeapMap;
+    NodeHeapMap* _node_heap;
+
+    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;
+    int _unmatched;
+
+    void createStructures() {
+      _node_num = countNodes(_bpgraph);
+
+      if (!_matching) {
+        _matching = new MatchingMap(_bpgraph);
+      }
+
+      if (!_node_potential) {
+        _node_potential = new NodePotential(_bpgraph);
+      }
+
+      if (!_status) {
+        _status = new StatusMap(_bpgraph);
+      }
+
+      if (!_pred) {
+        _pred = new PredMap(_bpgraph);
+      }
+
+      if (!_node_heap) {
+        _node_heap_index = new IntArcMap(_bpgraph);
+        _node_heap = new NodeHeapMap(_bpgraph, 0);
+      } else {
+        for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+          if ((*_status)[n] == MATCHED) {
+            delete (*_node_heap)[n];
+          }
+        }
+      }
+
+      if (!_tree_set) {
+        _tree_set_index = new IntNodeMap(_bpgraph);
+        _tree_set = new TreeSet(*_tree_set_index);
+      }
+
+      if (!_delta1) {
+        _delta1_index = new IntNodeMap(_bpgraph);
+        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
+      }
+
+      if (!_delta2) {
+        _delta2_index = new IntNodeMap(_bpgraph);
+        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
+      }
+
+      if (!_delta3) {
+        _delta3_index = new IntEdgeMap(_bpgraph);
+        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+      }
+    }
+
+    void destroyStructures() {
+      if (_matching) {
+        delete _matching;
+      }
+      if (_node_potential) {
+        delete _node_potential;
+      }
+      if (_node_heap) {
+        for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+          if ((*_status)[n] == MATCHED) {
+            delete (*_node_heap)[n];
+          }
+        }
+        delete _node_heap_index;
+        delete _node_heap;
+      }
+      if (_status) {
+        delete _status;
+      }
+      if (_pred) {
+        delete _pred;
+      }
+      if (_tree_set) {
+        delete _tree_set_index;
+        delete _tree_set;
+      }
+      if (_delta2) {
+        delete _delta2_index;
+        delete _delta2;
+      }
+      if (_delta3) {
+        delete _delta3_index;
+        delete _delta3;
+      }
+    }
+
+    void matchedToEven(Node node, int tree) {
+      _tree_set->insert(node, tree);
+      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
+      _delta1->push(node, (*_node_potential)[node]);
+
+      if (_delta2->state(node) == _delta2->IN_HEAP) {
+        _delta2->erase(node);
+      }
+      delete (*_node_heap)[node];
+
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
+          dualScale * _weight[a];
+        if ((*_status)[v] == EVEN) {
+          _delta3->push(a, rw / 2);
+        } else if ((*_status)[v] == MATCHED) {
+          (*_node_heap)[v]->push(a, rw);
+          if (_delta2->state(v) != _delta2->IN_HEAP) {
+            _delta2->push(v, rw);
+          } else if ((*_delta2)[v] > rw) {
+            _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);
+        delete (*_node_heap)[node];
+      }
+    }
+
+    void evenToMatched(Node node) {
+      _delta1->erase(node);
+      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
+      (*_node_heap)[node] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
+          dualScale * _weight[a];
+
+        if ((*_status)[v] == EVEN) {
+          _delta3->erase(a);
+
+          (*_node_heap)[node]->push(_bpgraph.oppositeArc(a), rw);
+        } else if ((*_status)[v]  == MATCHED) {
+          (*_node_heap)[v]->erase(a);
+          if ((*_node_heap)[v]->empty()) {
+            _delta2->erase(v);
+          } else if ((*_node_heap)[v]->prio() > (*_delta2)[v]) {
+            _delta2->increase(v, (*_node_heap)[v]->prio());
+          }
+        }
+      }
+      if (!(*_node_heap)[node]->empty()) {
+        _delta2->push(node, (*_node_heap)[node]->prio());
+      }
+    }
+
+    void oddToMatched(Node node) {
+      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
+      (*_node_heap)[node] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        if ((*_status)[v] != EVEN) continue;
+        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
+          dualScale * _weight[a];
+        (*_node_heap)[node]->push(_bpgraph.oppositeArc(a), rw);
+      }
+      if (!(*_node_heap)[node]->empty()) {
+        _delta2->push(node, (*_node_heap)[node]->prio());
+      }
+    }
+
+    void alternatePath(Node even) {
+      Node odd;
+
+      _status->set(even, MATCHED);
+      evenToMatched(even);
+
+      Arc prev = (*_matching)[even];
+      while (prev != INVALID) {
+        odd = _bpgraph.target(prev);
+        even = _bpgraph.target((*_pred)[odd]);
+        _matching->set(odd, (*_pred)[odd]);
+        _status->set(odd, MATCHED);
+        oddToMatched(odd);
+
+        prev = (*_matching)[even];
+        _status->set(even, MATCHED);
+        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
+        evenToMatched(even);
+      }
+    }
+
+    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);
+        } 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);
+      destroyTree(tree);
+
+      _matching->set(node, INVALID);
+    }
+
+    void augmentOnEdge(const Edge& edge) {
+      Node left = _bpgraph.u(edge);
+      int left_tree = _tree_set->find(left);
+
+      alternatePath(left);
+      destroyTree(left_tree);
+      _matching->set(left, _bpgraph.direct(edge, true));
+
+      Node right = _bpgraph.v(edge);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.direct(edge, false));
+    }
+
+    void augmentOnArc(const Arc& arc) {
+      Node left = _bpgraph.source(arc);
+      _matching->set(left, arc);
+      _pred->set(left, arc);
+
+      Node right = _bpgraph.target(arc);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.oppositeArc(arc));
+    }
+
+    void extendOnArc(const Arc& arc) {
+      Node base = _bpgraph.target(arc);
+      int tree = _tree_set->find(base);
+
+      Node odd = _bpgraph.source(arc);
+      _tree_set->insert(odd, tree);
+      _status->set(odd, ODD);
+      matchedToOdd(odd, tree);
+      _pred->set(odd, arc);
+
+      Node even = _bpgraph.target((*_matching)[odd]);
+      _tree_set->insert(even, tree);
+      _status->set(even, EVEN);
+      matchedToEven(even, tree);
+    }
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor.
+    MaxWeightedBpMatching3(const BpGraph& bpgraph, const WeightMap& weight)
+      : _bpgraph(bpgraph), _weight(weight), _matching(0),
+        _node_potential(0), _node_num(0),
+        _status(0), _pred(0), _node_heap_index(0), _node_heap(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(), _unmatched(0)
+    {}
+
+    ~MaxWeightedBpMatching3() {
+      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(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = _node_num;
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) / 2 > max) {
+            max = (dualScale * _weight[a]) / 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(_bpgraph); e != INVALID; ++e) {
+        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
+                          (*_node_potential)[_bpgraph.v(e)] -
+                          dualScale * _weight[e]) / 2);
+      }
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void redRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) > max) {
+            max = dualScale * _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+
+        _tree_set->insert(n);
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+
+        ++_unmatched;
+      }
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, MATCHED);
+        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.target(a);
+          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
+                     dualScale * _weight[a];
+          (*_node_heap)[n]->push(a, rw);
+        }
+        if (!(*_node_heap)[n]->empty()) {
+          _delta2->push(n, (*_node_heap)[n]->prio());
+        }
+      }
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void blueRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) > max) {
+            max = dualScale * _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+
+        _tree_set->insert(n);
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+
+        ++_unmatched;
+      }
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, MATCHED);
+        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.target(a);
+          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
+                     dualScale * _weight[a];
+          (*_node_heap)[n]->push(a, rw);
+        }
+        if (!(*_node_heap)[n]->empty()) {
+          _delta2->push(n, (*_node_heap)[n]->prio());
+        }
+      }
+    }
+
+    /// \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
+      };
+
+      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 = (*_node_heap)[n]->top();
+            if ((*_matching)[n] == INVALID) {
+              augmentOnArc(a);
+              --_unmatched;
+            } else {
+              extendOnArc(a);
+            }
+          }
+          break;
+        case D3:
+          {
+            Edge e = _delta3->top();
+            augmentOnEdge(e);
+            _unmatched -= 2;
+          }
+          break;
+        }
+      }
+    }
+
+    /// \brief Run the algorithm.
+    ///
+    /// This method runs the \c %MaxWeightedBpMatching3 algorithm.
+    ///
+    /// \note mwbpm.run() is just a shortcut of the following code.
+    /// \code
+    ///   mwbpm.init();
+    ///   mwbpm.start();
+    /// \endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Primal Solution
+    /// Functions to get the primal solution, i.e. the maximum weighted
+    /// bipartite 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value matchingWeight() const {
+      Value sum = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          sum += _weight[(*_matching)[n]];
+        }
+      }
+      return sum / 2;
+    }
+
+    /// \brief Return the size (cardinality) of the matching.
+    ///
+    /// This function returns the size (cardinality) of the found matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    int matchingSize() const {
+      int num = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++num;
+        }
+      }
+      return num / 2;
+    }
+
+    /// \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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpgraph.u(edge)];
+    }
+
+    /// \brief Return the matching arc (or edge) incident to the given node.
+    ///
+    /// This function returns the 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Arc matching(const Node& node) const {
+      return (*_matching)[node];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the found
+    /// matching or \c INVALID if the node is not covered by the matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Node mate(const Node& node) const {
+      return (*_matching)[node] != INVALID ?
+        _bpgraph.target((*_matching)[node]) : INVALID;
+    }
+
+    /// @}
+
+    /// \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(_bpgraph); n != INVALID; ++n) {
+        sum += nodeValue(n);
+      }
+      return sum;
+    }
+
+    /// \brief Return the dual value (potential) of the given node.
+    ///
+    /// This function returns the dual value (potential) of the given node.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// @}
+
+  };
+
+#ifdef DOXYGEN
+  template <typename BPGR, typename WM>
+#else
+  template <typename BPGR,
+            typename WM = typename BPGR::template EdgeMap<int> >
+#endif
+  class MaxWeightedBpMatching4 {
+  public:
+
+    /// The graph type of the algorithm
+    typedef BPGR BpGraph;
+    /// 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 BpGraph::template NodeMap<typename BpGraph::Arc>
+    MatchingMap;
+
+    /// \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_BPGRAPH_TYPEDEFS(BpGraph);
+
+    typedef typename BpGraph::template NodeMap<Value> NodePotential;
+
+    const BpGraph& _bpgraph;
+    const WeightMap& _weight;
+
+    MatchingMap* _matching;
+
+    NodePotential* _node_potential;
+
+    int _node_num;
+
+    enum Status {
+      EVEN = -1, MATCHED = 0, ODD = 1
+    };
+
+    typedef typename BpGraph::template NodeMap<Status> StatusMap;
+    StatusMap* _status;
+
+    typedef typename BpGraph::template NodeMap<Arc> PredMap;
+    PredMap* _pred;
+
+    IntArcMap* _node_heap_index;
+    typedef typename BpGraph::template NodeMap<BinHeap<Value, IntArcMap>*> NodeHeapMap;
+    NodeHeapMap* _node_heap;
+
+    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;
+    int _unmatched;
+
+    void createStructures() {
+      _node_num = countNodes(_bpgraph);
+
+      if (!_matching) {
+        _matching = new MatchingMap(_bpgraph);
+      }
+
+      if (!_node_potential) {
+        _node_potential = new NodePotential(_bpgraph);
+      }
+
+      if (!_status) {
+        _status = new StatusMap(_bpgraph);
+      }
+
+      if (!_pred) {
+        _pred = new PredMap(_bpgraph);
+      }
+
+      if (!_node_heap) {
+        _node_heap_index = new IntArcMap(_bpgraph);
+        _node_heap = new NodeHeapMap(_bpgraph, 0);
+      }
+
+      if (!_tree_set) {
+        _tree_set_index = new IntNodeMap(_bpgraph);
+        _tree_set = new TreeSet(*_tree_set_index);
+      }
+
+      if (!_delta1) {
+        _delta1_index = new IntNodeMap(_bpgraph);
+        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
+      }
+
+      if (!_delta2) {
+        _delta2_index = new IntNodeMap(_bpgraph);
+        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
+      }
+
+      if (!_delta3) {
+        _delta3_index = new IntEdgeMap(_bpgraph);
+        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+      }
+    }
+
+    void destroyStructures() {
+      if (_matching) {
+        delete _matching;
+      }
+      if (_node_potential) {
+        delete _node_potential;
+      }
+      if (_node_heap) {
+        for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+          delete (*_node_heap)[n];
+        }
+        delete _node_heap_index;
+        delete _node_heap;
+      }
+      if (_status) {
+        delete _status;
+      }
+      if (_pred) {
+        delete _pred;
+      }
+      if (_tree_set) {
+        delete _tree_set_index;
+        delete _tree_set;
+      }
+      if (_delta2) {
+        delete _delta2_index;
+        delete _delta2;
+      }
+      if (_delta3) {
+        delete _delta3_index;
+        delete _delta3;
+      }
+    }
+
+    void matchedToEven(Node node, int tree) {
+      _tree_set->insert(node, tree);
+      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
+      _delta1->push(node, (*_node_potential)[node]);
+
+      if (_delta2->state(node) == _delta2->IN_HEAP) {
+        _delta2->erase(node);
+      }
+
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+        Value hrw = (*_node_potential)[node] - dualScale * _weight[a];
+        (*_node_heap)[v]->push(a, hrw);
+        Value rw = hrw + (*_node_potential)[v];
+        if ((*_status)[v] == EVEN) {
+          _delta3->push(a, rw / 2);
+        } else if ((*_status)[v] == MATCHED) {
+          if (_delta2->state(v) != _delta2->IN_HEAP) {
+            _delta2->push(v, rw);
+          } else if ((*_delta2)[v] > rw) {
+            _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) {
+      _delta1->erase(node);
+      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
+      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
+        Node v = _bpgraph.source(a);
+
+        (*_node_heap)[v]->erase(a);
+
+        if ((*_status)[v] == EVEN) {
+          _delta3->erase(a);
+
+        } else if ((*_status)[v]  == MATCHED) {
+          if ((*_node_heap)[v]->empty()) {
+            _delta2->erase(v);
+          } else if ((*_node_heap)[v]->prio() + (*_node_potential)[v] > (*_delta2)[v]) {
+            _delta2->increase(v, (*_node_heap)[v]->prio() + (*_node_potential)[v]);
+          }
+        }
+      }
+      if (!(*_node_heap)[node]->empty()) {
+        _delta2->push(
+            node, (*_node_heap)[node]->prio() + (*_node_potential)[node]);
+      }
+    }
+
+    void oddToMatched(Node node) {
+      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
+      if (!(*_node_heap)[node]->empty()) {
+        _delta2->push(node, (*_node_heap)[node]->prio() + (*_node_potential)[node]);
+      }
+    }
+
+    void alternatePath(Node even) {
+      Node odd;
+
+      _status->set(even, MATCHED);
+      evenToMatched(even);
+
+      Arc prev = (*_matching)[even];
+      while (prev != INVALID) {
+        odd = _bpgraph.target(prev);
+        even = _bpgraph.target((*_pred)[odd]);
+        _matching->set(odd, (*_pred)[odd]);
+        _status->set(odd, MATCHED);
+        oddToMatched(odd);
+
+        prev = (*_matching)[even];
+        _status->set(even, MATCHED);
+        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
+        evenToMatched(even);
+      }
+    }
+
+    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);
+        } 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);
+      destroyTree(tree);
+
+      _matching->set(node, INVALID);
+    }
+
+    void augmentOnEdge(const Edge& edge) {
+      Node left = _bpgraph.u(edge);
+      int left_tree = _tree_set->find(left);
+
+      alternatePath(left);
+      destroyTree(left_tree);
+      _matching->set(left, _bpgraph.direct(edge, true));
+
+      Node right = _bpgraph.v(edge);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.direct(edge, false));
+    }
+
+    void augmentOnArc(const Arc& arc) {
+      Node left = _bpgraph.source(arc);
+      _matching->set(left, arc);
+      _pred->set(left, arc);
+
+      Node right = _bpgraph.target(arc);
+      int right_tree = _tree_set->find(right);
+
+      alternatePath(right);
+      destroyTree(right_tree);
+      _matching->set(right, _bpgraph.oppositeArc(arc));
+    }
+
+    void extendOnArc(const Arc& arc) {
+      Node base = _bpgraph.target(arc);
+      int tree = _tree_set->find(base);
+
+      Node odd = _bpgraph.source(arc);
+      _tree_set->insert(odd, tree);
+      _status->set(odd, ODD);
+      matchedToOdd(odd, tree);
+      _pred->set(odd, arc);
+
+      Node even = _bpgraph.target((*_matching)[odd]);
+      _tree_set->insert(even, tree);
+      _status->set(even, EVEN);
+      matchedToEven(even, tree);
+    }
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor.
+    MaxWeightedBpMatching4(const BpGraph& bpgraph, const WeightMap& weight)
+      : _bpgraph(bpgraph), _weight(weight), _matching(0),
+        _node_potential(0), _node_num(0),
+        _status(0), _pred(0), _node_heap_index(0), _node_heap(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(), _unmatched(0)
+    {}
+
+    ~MaxWeightedBpMatching4() {
+      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(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = _node_num;
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) / 2 > max) {
+            max = (dualScale * _weight[a]) / 2;
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+        _tree_set->insert(n);
+
+        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          (*_node_heap)[_bpgraph.source(a)]->push(
+              a, (*_node_potential)[n] - dualScale * _weight[a]);
+        }
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
+                          (*_node_potential)[_bpgraph.v(e)] -
+                          dualScale * _weight[e]) / 2);
+      }
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void redRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) > max) {
+            max = dualScale * _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+
+        _tree_set->insert(n);
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+
+        ++_unmatched;
+      }
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, MATCHED);
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.target(a);
+          Value rw = (*_node_potential)[v] - dualScale * _weight[a];
+          (*_node_heap)[n]->push(a, rw);
+        }
+        if (!(*_node_heap)[n]->empty()) {
+          _delta2->push(n, (*_node_heap)[n]->prio() + (*_node_potential)[n]);
+        }
+      }
+    }
+
+    /// \brief Initialize the algorithm
+    ///
+    /// This function initializes the algorithm.
+    void blueRootInit() {
+      createStructures();
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+        (*_delta2_index)[n] = _delta1->PRE_HEAP;
+        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = 0;
+      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = 0;
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) > max) {
+            max = dualScale * _weight[a];
+          }
+        }
+        _node_potential->set(n, max);
+        _delta1->push(n, max);
+
+        _tree_set->insert(n);
+
+        _matching->set(n, INVALID);
+        _status->set(n, EVEN);
+
+        ++_unmatched;
+      }
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, MATCHED);
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          Node v = _bpgraph.target(a);
+          Value rw = (*_node_potential)[v] - dualScale * _weight[a];
+          (*_node_heap)[n]->push(a, rw);
+        }
+        if (!(*_node_heap)[n]->empty()) {
+          _delta2->push(n, (*_node_heap)[n]->prio() + (*_node_potential)[n]);
+        }
+      }
+    }
+
+    /// \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
+      };
+
+      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 = (*_node_heap)[n]->top();
+            if ((*_matching)[n] == INVALID) {
+              augmentOnArc(a);
+              --_unmatched;
+            } else {
+              extendOnArc(a);
+            }
+          }
+          break;
+        case D3:
+          {
+            Edge e = _delta3->top();
+            augmentOnEdge(e);
+            _unmatched -= 2;
+          }
+          break;
+        }
+      }
+    }
+
+    /// \brief Run the algorithm.
+    ///
+    /// This method runs the \c %MaxWeightedBpMatching4 algorithm.
+    ///
+    /// \note mwbpm.run() is just a shortcut of the following code.
+    /// \code
+    ///   mwbpm.init();
+    ///   mwbpm.start();
+    /// \endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Primal Solution
+    /// Functions to get the primal solution, i.e. the maximum weighted
+    /// bipartite 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value matchingWeight() const {
+      Value sum = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          sum += _weight[(*_matching)[n]];
+        }
+      }
+      return sum / 2;
+    }
+
+    /// \brief Return the size (cardinality) of the matching.
+    ///
+    /// This function returns the size (cardinality) of the found matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    int matchingSize() const {
+      int num = 0;
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++num;
+        }
+      }
+      return num / 2;
+    }
+
+    /// \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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpgraph.u(edge)];
+    }
+
+    /// \brief Return the matching arc (or edge) incident to the given node.
+    ///
+    /// This function returns the 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.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Arc matching(const Node& node) const {
+      return (*_matching)[node];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the found
+    /// matching or \c INVALID if the node is not covered by the matching.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Node mate(const Node& node) const {
+      return (*_matching)[node] != INVALID ?
+        _bpgraph.target((*_matching)[node]) : INVALID;
+    }
+
+    /// @}
+
+    /// \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(_bpgraph); n != INVALID; ++n) {
+        sum += nodeValue(n);
+      }
+      return sum;
+    }
+
+    /// \brief Return the dual value (potential) of the given node.
+    ///
+    /// This function returns the dual value (potential) of the given node.
+    ///
+    /// \pre Either run() or start() must be called before using this function.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// @}
+
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_BPMATCHING_H

lemon/matching.h

@@ -1183 +1183 @@
         (*_blossom_data)[even].next =
           _graph.oppositeArc((*_blossom_data)[odd].pred);
       }
-
     }
 
     void destroyTree(int tree) {
@@ -1586 +1585 @@
         _delta_sum(), _unmatched(0),
 
         _fractional(0)
+
     {}
 
     ~MaxWeightedMatching() {

test/CMakeLists.txt

@@ -17 +17 @@
   bellman_ford_test
   bfs_test
   bpgraph_test
+  bpmatching_test
   circulation_test
   connectivity_test
   counter_test

new file test/bpmatching_test.cc

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2012
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <iostream>
+#include <sstream>
+#include <list>
+
+#include <lemon/random.h>
+#include <lemon/bpmatching.h>
+#include <lemon/matching.h>
+#include <lemon/smart_graph.h>
+#include <lemon/concepts/bpgraph.h>
+#include <lemon/concepts/maps.h>
+#include <lemon/lgf_reader.h>
+#include <lemon/time_measure.h>
+#include <lemon/network_simplex.h>
+#include <lemon/cost_scaling.h>
+#include <lemon/capacity_scaling.h>
+
+#include "test_tools.h"
+
+using namespace std;
+using namespace lemon;
+
+const int lgfn = 6;
+const string lgf[lgfn] = {
+  // Empty graph
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label weight\n"
+  "\n",
+
+  // No red nodes
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "1\n"
+  "@edges\n"
+  "  label weight\n"
+  "\n",
+
+  // No blue nodes
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label weight\n"
+  "\n",
+
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "@blue_nodes\n"
+  "label\n"
+  "2\n"
+  "\n"
+  "@edges\n"
+  "  label weight\n"
+  "1 2 1 1\n"
+  "\n",
+
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "3\n"
+  "4\n"
+  "\n"
+  "@edges\n"
+  "  label weight\n"
+  "1 3 1 5\n"
+  "1 4 2 3\n"
+  "2 3 3 1\n"
+  "2 4 4 2\n"
+  "\n",
+
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "3\n"
+  "4\n"
+  "\n"
+  "@edges\n"
+  "  label weight\n"
+  "1 3 1 3\n"
+  "1 4 2 2\n"
+  "2 3 3 2\n"
+  "\n",
+};
+
+void checkMaxWeightedBpMatchingCompiles() {
+  typedef concepts::BpGraph BpGraph;
+  typedef BpGraph::Node Node;
+  typedef BpGraph::Edge Edge;
+  typedef BpGraph::EdgeMap<int> IntEdgeMap;
+  typedef MaxWeightedBpMatching1<BpGraph, IntEdgeMap> MWBPM;
+
+  int v;
+  Node n;
+  Edge e;
+  bool b;
+
+  BpGraph graph;
+  IntEdgeMap weight(graph);
+  MWBPM mwbpm(graph, weight);
+  const MWBPM& cmwbpm = mwbpm;
+
+  mwbpm.init();
+  mwbpm.start();
+  mwbpm.run();
+
+  v = cmwbpm.matchingWeight();
+  v = cmwbpm.matchingSize();
+  e = cmwbpm.matching(n);
+  b = cmwbpm.matching(e);
+  n = cmwbpm.mate(n);
+  v = cmwbpm.dualValue();
+  v = cmwbpm.nodeValue(n);
+
+  ::lemon::ignore_unused_variable_warning(v, b);
+}
+
+
+template <typename BpMatching>
+void checkMaxWeightedBpMatching(const SmartBpGraph& bpgraph,
+                                const SmartBpGraph::EdgeMap<int>& weight,
+                                const BpMatching& mwbpm) {
+  int pv = 0;
+  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
+    int rw = mwbpm.nodeValue(bpgraph.redNode(e)) +
+             mwbpm.nodeValue(bpgraph.blueNode(e));
+    rw -= weight[e] * mwbpm.dualScale;
+    check(rw >= 0, "Negative reduced weight");
+    check(rw == 0 || !mwbpm.matching(e),
+          "Non-zero reduced weight on matching edge");
+    if (mwbpm.matching(e)) {
+      pv += weight[e];
+    }
+  }
+  int dv = 0;
+  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
+    if (mwbpm.matching(n) != INVALID) {
+      check(mwbpm.nodeValue(n) >= 0, "Invalid node value");
+      SmartBpGraph::Node o = bpgraph.target(mwbpm.matching(n));
+      check(mwbpm.mate(n) == o, "Invalid matching");
+      check(mwbpm.matching(n) == bpgraph.oppositeArc(mwbpm.matching(o)),
+            "Invalid matching");
+    } else {
+      check(mwbpm.mate(n) == INVALID, "Invalid matching");
+      check(mwbpm.nodeValue(n) == 0, "Invalid matching");
+    }
+  }
+
+  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
+    dv += mwbpm.nodeValue(n);
+  }
+
+  check(pv * mwbpm.dualScale == dv, "Wrong duality");
+}
+
+void generateRandomBpGraph(int redNum, int blueNum, int edgeNum, int min, int max,
+                           SmartBpGraph& bpgraph, SmartBpGraph::EdgeMap<int>& weight) {
+  std::vector<SmartBpGraph::RedNode> redNodes;
+  for (int i = 0; i < redNum; ++i) {
+    redNodes.push_back(bpgraph.addRedNode());
+  }
+  std::vector<SmartBpGraph::BlueNode> blueNodes;
+  for (int i = 0; i < blueNum; ++i) {
+    blueNodes.push_back(bpgraph.addBlueNode());
+  }
+  for (int i = 0; i < edgeNum; ++i) {
+    SmartBpGraph::Edge e = bpgraph.addEdge(
+        redNodes[lemon::rnd[redNum]], blueNodes[lemon::rnd[blueNum]]);
+    weight[e] = lemon::rnd[max - min] + min;
+  }
+}
+
+void generateFullBpGraph(int redNum, int blueNum, int min, int max,
+                         SmartBpGraph& bpgraph, SmartBpGraph::EdgeMap<int>& weight) {
+  std::vector<SmartBpGraph::RedNode> redNodes;
+  for (int i = 0; i < redNum; ++i) {
+    redNodes.push_back(bpgraph.addRedNode());
+  }
+  std::vector<SmartBpGraph::BlueNode> blueNodes;
+  for (int i = 0; i < blueNum; ++i) {
+    blueNodes.push_back(bpgraph.addBlueNode());
+  }
+  for (int i = 0; i < redNum; ++i) {
+    for (int j = 0; j < blueNum; ++j) {
+      SmartBpGraph::Edge e = bpgraph.addEdge(
+          redNodes[i], blueNodes[j]);
+      weight[e] = lemon::rnd[max - min] + min;
+    }
+  }
+}
+
+void generateSkewBpGraph(int partNodeNum, int deg,
+                         SmartBpGraph& bpgraph,
+                         SmartBpGraph::EdgeMap<int>& weight) {
+  std::vector<SmartBpGraph::RedNode> redNodes;
+  for (int i = 0; i < partNodeNum; ++i) {
+    redNodes.push_back(bpgraph.addRedNode());
+  }
+  for (int i = 0; i < partNodeNum; ++i) {
+    std::swap(redNodes[i], redNodes[lemon::rnd[partNodeNum - i] + i]);
+  }
+  std::vector<SmartBpGraph::BlueNode> blueNodes;
+  for (int i = 0; i < partNodeNum; ++i) {
+    blueNodes.push_back(bpgraph.addBlueNode());
+  }
+  for (int i = 0; i < partNodeNum; ++i) {
+    std::swap(blueNodes[i], blueNodes[lemon::rnd[partNodeNum - i] + i]);
+  }
+  for (int d = 0; d < deg; ++d) {
+    for (int i = d; i < partNodeNum; ++i) {
+      SmartBpGraph::Edge e = bpgraph.addEdge(
+          redNodes[i], blueNodes[i - d]);
+      weight[e] = partNodeNum + d;
+    }
+  }
+}
+
+void convertToMinCostFlow(const SmartBpGraph& bpgraph,
+                          const SmartBpGraph::EdgeMap<int>& weight,
+                          SmartDigraph& digraph,
+                          SmartDigraph::NodeMap<int>& supply,
+                          SmartDigraph::ArcMap<int>& cost) {
+  SmartDigraph::Node outer = digraph.addNode();
+  supply.set(outer, countBlueNodes(bpgraph) - countRedNodes(bpgraph));
+
+  SmartBpGraph::NodeMap<SmartDigraph::Node> nodeRef(bpgraph);
+  for (SmartBpGraph::RedNodeIt n(bpgraph); n != INVALID; ++n) {
+    nodeRef[n] = digraph.addNode();
+    supply.set(nodeRef[n], 1);
+    cost.set(digraph.addArc(nodeRef[n], outer), 0);
+  }
+  for (SmartBpGraph::BlueNodeIt n(bpgraph); n != INVALID; ++n) {
+    nodeRef[n] = digraph.addNode();
+    supply.set(nodeRef[n], -1);
+    cost.set(digraph.addArc(outer, nodeRef[n]), 0);
+  }
+
+  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
+    cost.set(
+        digraph.addArc(nodeRef[bpgraph.redNode(e)], nodeRef[bpgraph.blueNode(e)]),
+        - weight[e]);
+  }
+}
+
+
+void runAll(const SmartBpGraph& bpgraph,
+            const SmartBpGraph::EdgeMap<int>& weight,
+            bool includeMcf=true) {
+  SmartDigraph mcfDigraph;
+  SmartDigraph::NodeMap<int> mcfSupply(mcfDigraph);
+  SmartDigraph::ArcMap<int> mcfCost(mcfDigraph);
+  convertToMinCostFlow(bpgraph, weight, mcfDigraph, mcfSupply, mcfCost);
+
+  {
+    Timer t;
+    MaxWeightedMatching<SmartBpGraph> mwm(bpgraph, weight);
+    mwm.init();
+    mwm.start();
+    std::cerr << "MaxWeightedMatching/init: " << t << std::endl;
+  }
+  {
+    Timer t;
+    MaxWeightedMatching<SmartBpGraph> mwm(bpgraph, weight);
+    mwm.fractionalInit();
+    mwm.start();
+    std::cerr << "MaxWeightedMatching/fractionalInit: " << t << std::endl;
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching1<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.init();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching1/init: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching1<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.redRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching1/redRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching1<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.blueRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching1/blueRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching2<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.redRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching2/redRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching3<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.init();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching3/init: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching3<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.redRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching3/redRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching3<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.blueRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching3/blueRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching4<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.init();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching4/init: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching4<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.redRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching4/redRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  {
+    Timer t;
+    MaxWeightedBpMatching4<SmartBpGraph> mwbpm(bpgraph, weight);
+    mwbpm.blueRootInit();
+    mwbpm.start();
+    std::cerr << "MaxWeightedBpMatching4/blueRootInit: " << t << std::endl;
+    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+  }
+  if (includeMcf) {
+    Timer t;
+    NetworkSimplex<SmartDigraph> ns(mcfDigraph);
+    ns.costMap(mcfCost);
+    ns.supplyMap(mcfSupply);
+    ns.upperMap(constMap<SmartDigraph::Arc, int>(1));
+    ns.run();
+    std::cerr << "NetworkSimplex: " << t << std::endl;
+  }
+  if (includeMcf) {
+    Timer t;
+    NetworkSimplex<SmartDigraph> ns(mcfDigraph);
+    ns.costMap(mcfCost);
+    ns.supplyMap(mcfSupply);
+    ns.upperMap(constMap<SmartDigraph::Arc, int>(1));
+    ns.run(NetworkSimplex<SmartDigraph>::ALTERING_LIST);
+    std::cerr << "NetworkSimplex/ALTERING_LIST: " << t << std::endl;
+  }
+  if (includeMcf) {
+    Timer t;
+    CapacityScaling<SmartDigraph> cs(mcfDigraph);
+    cs.costMap(mcfCost);
+    cs.supplyMap(mcfSupply);
+    cs.upperMap(constMap<SmartDigraph::Arc, int>(1));
+    cs.run();
+    std::cerr << "CapacityScaling: " << t << std::endl;
+  }
+  if (includeMcf) {
+    Timer t;
+    CostScaling<SmartDigraph> cs(mcfDigraph);
+    cs.costMap(mcfCost);
+    cs.supplyMap(mcfSupply);
+    cs.upperMap(constMap<SmartDigraph::Arc, int>(1));
+    cs.run();
+    std::cerr << "CostScaling: " << t << std::endl;
+  }
+}
+
+
+int main() {
+  // // Checking hard wired extremal graphs
+  // for (int i=0; i<lgfn; ++i) {
+  //   SmartBpGraph bpgraph;
+  //   SmartBpGraph::EdgeMap<int> weight(bpgraph);
+  //   istringstream lgfs(lgf[i]);
+  //   bpGraphReader(bpgraph, lgfs)
+  //       .edgeMap("weight", weight).run();
+  //   runAll(bpgraph, weight);
+  // }
+
+  std::cerr << "Random red=10000 blue=10000 edge=10000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(10000, 10000, 10000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=20000 blue=20000 edge=20000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(20000, 20000, 20000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=40000 blue=40000 edge=40000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(40000, 40000, 40000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=100000 blue=100000 edge=100000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(100000, 100000, 100000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=1000000 blue=1000000 edge=1000000 range=0-1000" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(1000000, 1000000, 1000000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=10000 blue=10000 edge=10000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(10000, 10000, 10000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=20000 blue=20000 edge=20000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(20000, 20000, 20000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=40000 blue=40000 edge=40000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(40000, 40000, 40000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=100000 blue=100000 edge=100000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(100000, 100000, 100000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=1000000 blue=1000000 edge=1000000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(1000000, 1000000, 1000000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(1000000, 1000000, 2000000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=10000 blue=10000 edge=20000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(10000, 10000, 20000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=20000 blue=20000 edge=40000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(20000, 20000, 40000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=100000 blue=100000 edge=200000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(100000, 100000, 200000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(1000000, 1000000, 2000000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=10000 blue=10000 edge=20000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(10000, 10000, 20000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=20000 blue=20000 edge=40000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(20000, 20000, 40000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=100000 blue=100000 edge=200000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(100000, 100000, 200000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(1000000, 1000000, 2000000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=1000-1100" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(1000000, 1000000, 2000000, 1000, 1100, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=10000 blue=10000 edge=200000 range=0-10000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(10000, 10000, 200000, 0, 10000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=20000 blue=20000 edge=400000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(20000, 20000, 400000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=40000 blue=40000 edge=800000 range=0-1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(40000, 40000, 800000, 0, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Random red=100000 blue=100000 edge=2000000 range=0-10000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateRandomBpGraph(100000, 100000, 2000000, 0, 10000, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=10000 blue=10000 deg=2" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(10000, 2, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=20000 blue=20000 deg=2" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(20000, 2, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=100000 blue=100000 deg=2" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(100000, 2, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=1000000 blue=1000000 deg=2" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(1000000, 2, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=10000 blue=10000 deg=3" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(10000, 3, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=20000 blue=20000 deg=3" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(20000, 3, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=100000 blue=100000 deg=3" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(100000, 3, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=1000000 blue=1000000 deg=3" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(1000000, 3, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=10000 blue=10000 deg=20" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(10000, 20, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=20000 blue=20000 deg=20" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(20000, 20, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=100000 blue=100000 deg=20" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(100000, 20, bpgraph, weight);
+    runAll(bpgraph, weight, false);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Skew red=1000 blue=1000 deg=1000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateSkewBpGraph(1000, 1000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Full red=4000 blue=1000 range=10000-11000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateFullBpGraph(4000, 1000, 10000, 11000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  std::cerr << "Full red=1000 blue=4000 range=10000-11000" << std::endl;
+  std::cerr << "=================================================================" << std::endl;
+  for (int i = 0; i < 3; ++i) {
+    std::cerr << "===== case " << i << ":" << std::endl;
+    SmartBpGraph bpgraph;
+    SmartBpGraph::EdgeMap<int> weight(bpgraph);
+    generateFullBpGraph(1000, 4000, 10000, 11000, bpgraph, weight);
+    runAll(bpgraph, weight);
+  }
+  std::cerr << std::endl;
+
+  return 0;
+}