Ticket #168: 5b4317f3ce36.patch

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-2019
+ * 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 MaxWeightedBpMatching {
+  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.
+    MaxWeightedBpMatching(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)
+    {}
+
+    ~MaxWeightedBpMatching() {
+      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] = _delta2->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] = _delta2->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] = _delta2->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 %MaxWeightedBpMatching 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 (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          sum += _weight[(*_matching)[n]];
+        }
+      }
+      return sum;
+    }
+
+    /// \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 (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++num;
+        }
+      }
+      return num;
+    }
+
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the found
+    /// matching.
+    ///
+    /// \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];
+    }
+
+    /// @}
+  };
+
+  /// \ingroup matching
+  ///
+  /// \brief Weighted perfect 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 exactly one incident edge.
+  /// It can be formulated with the following linear program.
+  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
+  /// \f[x_e \ge 0 \quad \forall e\in E\f]
+  /// \f[\max \sum_{e\in E}x_ew_e\f]
+  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+  /// \f$X\f$.
+  ///
+  /// The 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[\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 MaxWeightedPerfectBpMatching {
+  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 *_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 (!_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);
+
+      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) {
+      _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.
+    MaxWeightedPerfectBpMatching(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),
+
+        _delta2_index(0), _delta2(0),
+        _delta3_index(0), _delta3(0),
+
+        _delta_sum(), _unmatched(0)
+    {}
+
+    ~MaxWeightedPerfectBpMatching() {
+      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.
+    ///
+    /// \return If it is false, then the graph does not have a perfect matching.
+    bool init() {
+      createStructures();
+
+      if (countRedNodes(_bpgraph) != countBlueNodes(_bpgraph)) {
+        return false;
+      }
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        (*_delta2_index)[n] = _delta2->PRE_HEAP;
+      }
+      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+
+      _delta2->clear();
+      _delta3->clear();
+      _tree_set->clear();
+
+      _unmatched = _node_num;
+
+      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
+        Value max = - std::numeric_limits<Value>::max();
+        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
+          if ((dualScale * _weight[a]) / 2 > max) {
+            max = (dualScale * _weight[a]) / 2;
+          }
+        }
+        if (max == - std::numeric_limits<Value>::max()) {
+          return false;
+        }
+        _node_potential->set(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);
+      }
+      return true;
+    }
+
+    /// \brief Start the algorithm
+    ///
+    /// This function starts the algorithm.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    ///
+    /// \return True when a perfect matching is found.
+    bool start() {
+      enum OpType {
+        D2, D3
+      };
+
+      while (_unmatched > 0) {
+        Value d2 = !_delta2->empty() ?
+          _delta2->prio() : std::numeric_limits<Value>::max();
+
+        Value d3 = !_delta3->empty() ?
+          _delta3->prio() : std::numeric_limits<Value>::max();
+
+        _delta_sum = d3; OpType ot = D3;
+        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
+
+        if (_delta_sum == std::numeric_limits<Value>::max()) {
+          return false;
+        }
+
+        switch (ot) {
+        case D2:
+          {
+            Node n = _delta2->top();
+            Arc a = (*_pred)[n];
+            if ((*_matching)[n] == INVALID) {
+              augmentOnArc(a);
+              --_unmatched;
+            } else {
+              extendOnArc(a);
+            }
+          }
+          break;
+        case D3:
+          {
+            Edge e = _delta3->top();
+            augmentOnEdge(e);
+            _unmatched -= 2;
+          }
+          break;
+        }
+      }
+      return true;
+    }
+
+    /// \brief Run the algorithm.
+    ///
+    /// This method runs the \c %MaxWeightedPerfectBpMatching algorithm.
+    ///
+    /// \note mwpbpm.run() is just a shortcut of the following code.
+    /// \code
+    ///   return mwpbpm.init() && mwpbpm.start();
+    /// \endcode
+    ///
+    /// \return True when a perfect matching is found.
+    bool run() {
+      return 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 and a matching should be found.
+
+    /// @{
+
+    /// \brief Return the weight of the matching.
+    ///
+    /// This function returns the weight of the found matching.
+    ///
+    /// \pre A perfect matching has been found.
+    Value matchingWeight() const {
+      Value sum = 0;
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+        sum += _weight[(*_matching)[n]];
+      }
+      return sum;
+    }
+
+    /// \brief Return the size (cardinality) of the matching.
+    ///
+    /// This function returns the size (cardinality) of the found matching.
+    ///
+    /// \pre A perfect matching has been found.
+    int matchingSize() const {
+      int num = 0;
+      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
+          ++num;
+      }
+      return num;
+    }
+
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the found
+    /// matching.
+    ///
+    /// \pre A perfect matching has been found.
+    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 A perfect matching has been found.
+    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 and a matching should be found.
+
+    /// @{
+
+    /// \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 A perfect matching has been found.
+    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 A perfect matching has been found.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// @}
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_BPMATCHING_H

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 "test_tools.h"
+
+using namespace std;
+using namespace lemon;
+
+const int lgfn = 7;
+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",
+
+  "@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 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 MaxWeightedBpMatching<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.redRootInit();
+  mwbpm.blueRootInit();
+  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);
+}
+
+void checkMaxWeightedPerfectBpMatchingCompiles() {
+  typedef concepts::BpGraph BpGraph;
+  typedef BpGraph::Node Node;
+  typedef BpGraph::Edge Edge;
+  typedef BpGraph::EdgeMap<int> IntEdgeMap;
+  typedef MaxWeightedPerfectBpMatching<BpGraph, IntEdgeMap> MWPBPM;
+
+  int v;
+  Node n;
+  Edge e;
+  bool b;
+
+  BpGraph graph;
+  IntEdgeMap weight(graph);
+  MWPBPM mwpbpm(graph, weight);
+  const MWPBPM& cmwpbpm = mwpbpm;
+
+  mwpbpm.init();
+  mwpbpm.start();
+  mwpbpm.run();
+
+  v = cmwpbpm.matchingWeight();
+  v = cmwpbpm.matchingSize();
+  e = cmwpbpm.matching(n);
+  b = cmwpbpm.matching(e);
+  n = cmwpbpm.mate(n);
+  v = cmwpbpm.dualValue();
+  v = cmwpbpm.nodeValue(n);
+
+  ::lemon::ignore_unused_variable_warning(v, b);
+}
+
+
+void checkMaxWeightedBpMatching(
+    const SmartBpGraph& bpgraph,
+    const SmartBpGraph::EdgeMap<int>& weight,
+    const MaxWeightedBpMatching<SmartBpGraph>& 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 checkPerfectBpMatchingExists(
+    const SmartBpGraph& bpgraph,
+    bool exists) {
+  if (countRedNodes(bpgraph) == countBlueNodes(bpgraph)) {
+    // TODO: Use a bipartite maximum cardinality matching algorithm here.
+    MaxMatching<SmartBpGraph> mbpm(bpgraph);
+    mbpm.run();
+    check(exists == (countRedNodes(bpgraph) == mbpm.matchingSize()),
+          "Unexpected maximum cardinality matching size");
+  } else {
+    check(!exists, "Different size for red and blue node set");
+  }
+}
+
+void checkMaxWeightedPerfectBpMatching(
+    const SmartBpGraph& bpgraph,
+    const SmartBpGraph::EdgeMap<int>& weight,
+    const MaxWeightedPerfectBpMatching<SmartBpGraph>& mwpbpm) {
+  int pv = 0;
+  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
+    int rw = mwpbpm.nodeValue(bpgraph.redNode(e)) +
+             mwpbpm.nodeValue(bpgraph.blueNode(e));
+    rw -= weight[e] * mwpbpm.dualScale;
+    check(rw >= 0, "Negative reduced weight");
+    check(rw == 0 || !mwpbpm.matching(e),
+          "Non-zero reduced weight on matching edge");
+    if (mwpbpm.matching(e)) {
+      pv += weight[e];
+    }
+  }
+  int dv = 0;
+  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
+    check(mwpbpm.matching(n) != INVALID, "Matching is not perfect");
+    SmartBpGraph::Node o = bpgraph.target(mwpbpm.matching(n));
+    check(mwpbpm.mate(n) == o, "Invalid matching");
+    check(mwpbpm.matching(n) == bpgraph.oppositeArc(mwpbpm.matching(o)),
+          "Invalid matching");
+  }
+
+  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
+    dv += mwpbpm.nodeValue(n);
+  }
+
+  check(pv * mwpbpm.dualScale == dv, "Wrong duality");
+}
+
+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();
+
+    {
+      MaxWeightedBpMatching<SmartBpGraph> mwbpm(bpgraph, weight);
+      mwbpm.run();
+      checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+    }
+    {
+      MaxWeightedBpMatching<SmartBpGraph> mwbpm(bpgraph, weight);
+      mwbpm.redRootInit();
+      mwbpm.start();
+      checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+    }
+    {
+      MaxWeightedBpMatching<SmartBpGraph> mwbpm(bpgraph, weight);
+      mwbpm.blueRootInit();
+      mwbpm.start();
+      checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
+    }
+
+    {
+      MaxWeightedPerfectBpMatching<SmartBpGraph> mwpbpm(bpgraph, weight);
+      bool exists = mwpbpm.run();
+      checkPerfectBpMatchingExists(bpgraph, exists);
+      if (exists) {
+        checkMaxWeightedPerfectBpMatching(bpgraph, weight, mwpbpm);
+      }
+    }
+  }
+
+  return 0;
+}