Ticket #168: bpm_pack.patch

lemon/bp_matching.h

@@ -45 +45 @@
     /// The graph type of the algorithm
     typedef BGR BpGraph;
     /// The type of the matching map on the red nodes
-    typedef typename BpGraph::template RedMap<typename BpGraph::Edge>
+    typedef typename BpGraph::template RedNodeMap<typename BpGraph::Edge>
     RedMatchingMap;
     /// The type of the matching map on the blue nodes
-    typedef typename BpGraph::template BlueMap<typename BpGraph::Edge>
+    typedef typename BpGraph::template BlueNodeMap<typename BpGraph::Edge>
     BlueMatchingMap;
 
   private:
@@ -118 +118 @@
       RM = new RedMatchingMap(*G, INVALID);
       BM = new BlueMatchingMap(*G, INVALID);
       
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         for( IncEdgeIt eit(*G, n); eit != INVALID; ++eit ) {
           if( (*BM)[ G->blueNode(eit) ] == INVALID ){
             (*BM)[G->blueNode(eit) ] = eit;
@@ -164 +164 @@
     /// called before using this function.
     bool step()
     {
-      queue<Node> sources;
+      queue<RedNode> sources;
       BoolNodeMap _processed(*G, false);
       EdgeNodeMap _pred(*G, INVALID);
       BoolNodeMap _reached(*G, false);
-      vector<Node> H;
+      vector<RedNode> H;
       bool isAlternatePath = false;
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         if( (*RM)[n] == INVALID ){
           sources.push( n );
         }
       }
         
       while ( !sources.empty() && !isAlternatePath ){
-        Node s = sources.front();
+        RedNode s = sources.front();
         if(_processed[s]) {sources.pop(); continue;}
         _processed[s] = true;
         _reached[s] = true;
         H.push_back( s );
-        Node m;
+        BlueNode m;
         for(IncEdgeIt e(*G,s);e!=INVALID;++e) {
           if(!_reached[m=G->blueNode(e)]) {
             _reached[m]=true;
             _pred[m]=e;
             if( (*BM)[m] != INVALID ) {
-              Node n;
+              RedNode n;
               if (!_reached[n = G->redNode((*BM)[m])] ) {
                 sources.push( n );
                 _pred[n] = (*BM)[m];
@@ -257 +257 @@
     int matchingSize() const
     {
       int size = 0;
-      for ( BlueIt n(*G); n != INVALID; ++n) {
+      for ( BlueNodeIt n(*G); n != INVALID; ++n) {
         if ( (*BM)[n] != INVALID) {
           ++size;
         }
@@ -281 +281 @@
     /// not covered by the matching.
     Edge matching(const Node& n) const
     {
-      return G->blue(n) ? (*BM)[n] : (*RM)[n];
+      return G->blue(n) ? (*BM)[G->asBlueNodeUnsafe(n)] : (*RM)[G->asRedNodeUnsafe(n)];
     }
 
     /// \brief Return a const reference to the blue matching map.
@@ -310 +310 @@
     {
       Edge e;
       if( G->blue(node) ){
-        if( ( e = (*BM)[node]) != INVALID ) return G->redNode(e);
+        if( ( e = (*BM)[G->asBlueNodeUnsafe(node)]) != INVALID ) return G->redNode(e);
         else return INVALID;
       }
       else {
-        if( ( e = (*RM)[node]) != INVALID ) return G->blueNode(e);
+        if( ( e = (*RM)[G->asRedNodeUnsafe(node)]) != INVALID ) return G->blueNode(e);
         else return INVALID;
       }
     }
@@ -351 +351 @@
     typedef typename WeightMap::Value Value;
     
     /// The type of the red matching map
-    typedef typename BpGraph::template RedMap<typename BpGraph::Edge>
+    typedef typename BpGraph::template RedNodeMap<typename BpGraph::Edge>
     RedMatchingMap;
     /// The type of the blue matching map
-    typedef typename BpGraph::template BlueMap<typename BpGraph::Edge>
+    typedef typename BpGraph::template BlueNodeMap<typename BpGraph::Edge>
     BlueMatchingMap;
 
   private:
@@ -420 +420 @@
       BM = new BlueMatchingMap(*G, INVALID);
       
       int nodeCount = 0;
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         Value max_C = 0;
         bool start = true;
         vector<Edge> maxEdges;
@@ -446 +446 @@
         }
         nodeCount++;
       }
-      for( BlueIt n(*G); n != INVALID; ++n ){
+      for( BlueNodeIt n(*G); n != INVALID; ++n ){
         Value minDif = 0;
         bool start = true;
         vector<Edge> minEdges;
@@ -482 +482 @@
     HungAlgStatus step()
     {
       bool perfect = true;
-      for( BlueIt n(*G); n != INVALID; ++n ){
+      for( BlueNodeIt n(*G); n != INVALID; ++n ){
         if( (*BM)[n] == INVALID ) {
           perfect = false;
           break;
         }
       }
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         if( (*RM)[n] == INVALID ) {
           perfect = false;
           break;
         }
       }
       if( perfect ) return HUNGALG_SUCCESS;
-      queue<Node> sources;
+      queue<RedNode> sources;
       BoolNodeMap _processed(*G, false);
       EdgeNodeMap _pred(*G, INVALID);
       BoolNodeMap _reached(*G, false);
-      vector<Node> H;
+      vector<RedNode> H;
       bool isAlternatePath = false;
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         if( (*RM)[n] == INVALID ){
           sources.push( n );
         }
       }
       
       while ( !sources.empty() && !isAlternatePath ){
-        Node s = sources.front();
+        RedNode s = sources.front();
         if(_processed[s]) {sources.pop(); continue;}
         _processed[s] = true;
         _reached[s] = true;
         H.push_back( s );
-        Node m;
+        BlueNode m;
         for(IncEdgeIt e(*G,s);e!=INVALID;++e) {
           if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
           if(!_reached[m=G->blueNode(e)]) {
             _reached[m]=true;
             _pred[m]=e;
             if( (*BM)[m] != INVALID ) {
-              Node n;
+              RedNode n;
               if (!_reached[n = G->redNode((*BM)[m])] ) {
                 sources.push( n );
                 _pred[n] = (*BM)[m];
@@ -537 +537 @@
       }
         
       if( !isAlternatePath ) {
-        vector<Node> GammaH;
-        BoolBlueMap GammaH_map( *G, false );
-        for (typename vector<Node>::iterator nit = H.begin();
+        vector<BlueNode> GammaH;
+        BoolBlueNodeMap GammaH_map( *G, false );
+        for (typename vector<RedNode>::iterator nit = H.begin();
              nit!=H.end(); ++nit){
           for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
             if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
-            Node n = G->blueNode( e );
+            BlueNode n = G->blueNode( e );
             if( !GammaH_map[n] ){
               GammaH.push_back(n);
               GammaH_map[n] = true;
@@ -552 +552 @@
         }
         Value gamma = 0;
         bool start = true;
-        for (typename vector<Node>::iterator nit = H.begin();
+        for (typename vector<RedNode>::iterator nit = H.begin();
              nit!=H.end(); ++nit){
           for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
-            Node n = G->oppositeNode( *nit, e );
+            BlueNode n = G->oppositeNode( *nit, e );
             if( !GammaH_map[ n ] ){
               Value dif = PI[ G->u(e) ] + PI[ G->v(e) ] - (*c)[e];
               if( start ){ gamma=dif; start=false; continue; }
@@ -567 +567 @@
             
         if(!gamma) return HUNGALG_ERROR;
             
-        for (typename vector<Node>::iterator nit = H.begin();
+        for (typename vector<RedNode>::iterator nit = H.begin();
              nit!=H.end(); ++nit){
           PI[*nit] -= gamma;
         }
-        for (typename vector<Node>::iterator nit = GammaH.begin();
+        for (typename vector<BlueNode>::iterator nit = GammaH.begin();
              nit!=GammaH.end(); ++nit){
           PI[*nit] += gamma;
         }
@@ -631 +631 @@
     {
       Value retVal = 0;
       Edge e;
-      for( BlueIt n(*G); n != INVALID; ++n )
+      for( BlueNodeIt n(*G); n != INVALID; ++n )
         if( ( e = (*BM)[n] ) != INVALID ) retVal += (*c)[e];
       return retVal;
     }
@@ -656 +656 @@
     /// \pre Either run() or start() must be called before using this function.
     Edge matching(const Node& node) const
     {
-      return G->blue(node) ? (*BM)[node] : (*RM)[node];
+      return G->blue(node) ? (*BM)[G->asBlueNodeUnsafe(node)] : (*RM)[G->asRedNodeUnsafe(node)];
     }
 
     /// \brief Return a const reference to the blue matching map.
@@ -687 +687 @@
     {
       Edge e;
       if( G->blue(node) ){
-        if( ( e = (*BM)[node]) != INVALID ) return G->redNode(e);
+        if( ( e = (*BM)[G->asBlueNodeUnsafe(node)]) != INVALID ) return G->redNode(e);
         else return INVALID;
       }
       else {
-        if( ( e = (*RM)[node]) != INVALID ) return G->blueNode(e);
+        if( ( e = (*RM)[G->asRedNodeUnsafe(node)]) != INVALID ) return G->blueNode(e);
         else return INVALID;
       }
     }
@@ -738 +738 @@
     /// The value type of the edge weights
     typedef typename WeightMap::Value Value;
     /// The type of the red matching map
-    typedef typename BpGraph::template RedMap<typename BpGraph::Edge>
+    typedef typename BpGraph::template RedNodeMap<typename BpGraph::Edge>
     RedMatchingMap;
     /// The type of the blue matching map
-    typedef typename BpGraph::template BlueMap<typename BpGraph::Edge>
+    typedef typename BpGraph::template BlueNodeMap<typename BpGraph::Edge>
     BlueMatchingMap;
     
   private:
@@ -750 +750 @@
     
     typedef typename BpGraph::template NodeMap<Value> NodeMap;
     
-    typedef typename BpGraph::template RedMap<BlueNode> BlueRedMap;
+    typedef typename BpGraph::template RedNodeMap<BlueNode> BlueRedMap;
     
-    typedef typename BpGraph::template BlueMap<RedNode> RedBlueMap;
+    typedef typename BpGraph::template BlueNodeMap<RedNode> RedBlueMap;
 
     typedef typename BpGraph::template NodeMap<Edge> EdgeNodeMap;
     
@@ -770 +770 @@
     RedBlueMap *RBM;
     BlueRedMap *BRM;
     
-    void _processMins( bool &processedMins, queue<Node> &sources, NodeNodeMap &_pred, BoolNodeMap &_reached, BlueNode &first )
+    void _processMins( bool &processedMins, queue<RedNode> &sources, NodeNodeMap &_pred, BoolNodeMap &_reached, BlueNode &first )
     {
       if( processedMins ) return;
       processedMins = true;
-      for( BlueIt m(*G); m!=INVALID; ++m ) {
+      for( BlueNodeIt m(*G); m!=INVALID; ++m ) {
         if( PI[m] == -MinPI ){
           _reached[m]=true;
           if( (*RBM)[m] != INVALID ) {
-            Node n;
+            RedNode n;
             if (!_reached[n = (*RBM)[m]] ) {
               sources.push( n );
               _pred[n] = m;
@@ -788 +788 @@
         }
       }
         
-      for( RedIt m(*G); m!=INVALID; ++m ) {
+      for( RedNodeIt m(*G); m!=INVALID; ++m ) {
         if( PI[m] == MinPI && !_reached[ m ] && (*BRM)[m] == INVALID ) {
           sources.push( m );
           _pred[m] = INVALID;
@@ -798 +798 @@
       
     }
     
-    void _setMatching( NodeNodeMap &pred, Node n, BlueNode &first ) {
-      Node node = n;
-      Node prevNode = pred[n];
+    void _setMatching( NodeNodeMap &pred, BlueNode n, BlueNode &first ) {
+      BlueNode node = n;
+      RedNode prevNode = G->asRedNodeUnsafe(pred[n]);
       while( 1 ) {
         if( prevNode == INVALID ){
           if( PI[node] == -MinPI ){
@@ -811 +811 @@
         else if( (*RBM)[node] != prevNode ) {
           (*BRM)[prevNode] = node;
           (*RBM)[node] = prevNode;
-          node = pred[prevNode];
+          node = G->asBlueNodeUnsafe(pred[prevNode]);
         }
         if( node == INVALID ) {
           if( ( PI[prevNode] == MinPI && (*BRM)[prevNode] == INVALID ) || first == INVALID ) return;
           else node = first;
         }
-        prevNode = pred[node];
+        prevNode = G->asRedNodeUnsafe(pred[node]);
       }
     }
     
@@ -858 +858 @@
       BRM = new BlueRedMap(*G, INVALID);
       RBM = new RedBlueMap(*G, INVALID);
       
-      for( BlueIt n(*G); n != INVALID; ++n )
+      for( BlueNodeIt n(*G); n != INVALID; ++n )
       {
         Value max_C = 0;
-        vector<Node> maxNodes;
+        vector<RedNode> maxNodes;
         for( IncEdgeIt e( *G, n ); e!=INVALID; ++e ) {
           if( (*c)[e] > max_C ) {
             max_C = (*c)[e];
@@ -873 +873 @@
           }
         }
         PI[n] = max_C;
-        for (typename vector<Node>::iterator eit = maxNodes.begin();
+        for (typename vector<RedNode>::iterator eit = maxNodes.begin();
           eit!=maxNodes.end(); ++eit)
         {
           if( (*BRM)[*eit] == INVALID )
@@ -884 +884 @@
           }
         }
       }
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         PI[n] = 0;
       }
     }
@@ -901 +901 @@
       for( NodeIt n(*G); n != INVALID; ++n ){
         
         if( G->blue(n) ) {
-          if( (*RBM)[n] == INVALID && PI[n] != -MinPI ){
+          if( (*RBM)[G->asBlueNodeUnsafe(n)] == INVALID && PI[n] != -MinPI ){
             Continue = true;
           }
         }
         else {
-          if( (*BRM)[n] == INVALID && PI[n] != MinPI ){
+          if( (*BRM)[G->asRedNodeUnsafe(n)] == INVALID && PI[n] != MinPI ){
             Continue = true;
           }
         }
       }
       if( !Continue ){
-        for( RedIt n(*G); n != INVALID; ++n ) {
+        for( RedNodeIt n(*G); n != INVALID; ++n ) {
           PI[n] -= MinPI;
           
           for( IncEdgeIt e(*G,n); e != INVALID; ++e ) {
             if( (*BRM)[n] == G->blueNode(e) ) (*RM)[n] = e;
           }
         }
-        for( BlueIt n(*G); n != INVALID; ++n ) {
+        for( BlueNodeIt n(*G); n != INVALID; ++n ) {
           PI[n] += MinPI;
           
           for( IncEdgeIt e(*G,n); e != INVALID; ++e ) {
@@ -929 +929 @@
         
         return true;
       }
-      queue<Node> sources;
+      queue<RedNode> sources;
       BoolNodeMap _processed(*G, false);
       NodeNodeMap _pred(*G, INVALID);
       BoolNodeMap _reached(*G, false);
-      vector<Node> H;
+      vector<RedNode> H;
       bool isAlternatePath = false;
       bool notMinUncovered = false;
       BlueNode first = INVALID;
       
-      for( BlueIt n(*G); n != INVALID; ++n ){
+      for( BlueNodeIt n(*G); n != INVALID; ++n ){
         if( (*RBM)[n] == INVALID && PI[n] != -MinPI ){
           notMinUncovered = true;
           break;
         }
       }
       
-      for( RedIt n(*G); n != INVALID; ++n ){
+      for( RedNodeIt n(*G); n != INVALID; ++n ){
         if( (*BRM)[n] == INVALID ){
           if( notMinUncovered || PI[n] != MinPI ) {
             sources.push( n );
@@ -961 +961 @@
       }
       
       while ( !sources.empty() && !isAlternatePath ){
-        Node s = sources.front();
+        RedNode s = sources.front();
         if(_processed[s]) { sources.pop(); continue; }
         if(PI[s] == MinPI && (*BRM)[s] == INVALID ) _processMins( processedMins, sources, _pred, _reached, first );
         _processed[s] = true;
         _reached[s] = true;
         H.push_back( s );
-        Node m;
+        BlueNode m;
         for( IncEdgeIt e(*G,s); e!=INVALID; ++e ) {
           if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
           if( _reached[m=G->blueNode(e)] ) continue;
           _reached[m]=true;
           _pred[m]=s;
           if( (*RBM)[m] != INVALID ) {
-            Node n;
+            RedNode n;
             if (!_reached[n = (*RBM)[m]] ) {
               sources.push( n );
               _pred[n] = m;
@@ -994 +994 @@
       }
       
       if( !isAlternatePath ) {
-        vector<Node> GammaH;
-        BoolBlueMap GammaH_map( *G, false );
-        for (typename vector<Node>::iterator nit = H.begin();
+        vector<BlueNode> GammaH;
+        BoolBlueNodeMap GammaH_map( *G, false );
+        for (typename vector<RedNode>::iterator nit = H.begin();
              nit!=H.end(); ++nit){
           for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
             if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
-            Node n = G->blueNode( e );
+            BlueNode n = G->blueNode( e );
             if( !GammaH_map[n] ){
               GammaH.push_back(n);
               GammaH_map[n] = true;
@@ -1008 +1008 @@
           }
         }
         if( processedMins ) {
-          for( BlueIt n(*G); n != INVALID; ++n ) {
+          for( BlueNodeIt n(*G); n != INVALID; ++n ) {
             if( PI[n] == -MinPI && !GammaH_map[n] ) {
               GammaH.push_back(n);
               GammaH_map[n] = true;
@@ -1018 +1018 @@
         
         Value gamma = 0;
         bool start = true;
-        for (typename vector<Node>::iterator nit = H.begin();
+        for (typename vector<RedNode>::iterator nit = H.begin();
              nit!=H.end(); ++nit){
           for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
             if((*c)[e]<0) continue;
-            Node n = G->oppositeNode( *nit, e );
+            BlueNode n = G->asBlueNodeUnsafe(G->oppositeNode( *nit, e ));
             if( !GammaH_map[ n ] ){
               Value dif = PI[ G->u(e) ] + PI[ G->v(e) ] - (*c)[e];
               if( start ){ gamma=dif; start=false; continue; }
@@ -1033 +1033 @@
         }
         
         if( processedMins ) {
-          for( BlueIt n(*G); n != INVALID; ++n ) {
+          for( BlueNodeIt n(*G); n != INVALID; ++n ) {
             if( !GammaH_map[n] && ( PI[n] + MinPI < gamma || !gamma ) ) gamma = PI[n] + MinPI;
           }
         }
             
-        for (typename vector<Node>::iterator nit = H.begin();
+        for (typename vector<RedNode>::iterator nit = H.begin();
              nit!=H.end(); ++nit){
           PI[*nit] -= gamma;
         }
         if( processedMins ) MinPI -= gamma;
-        for (typename vector<Node>::iterator nit = GammaH.begin();
+        for (typename vector<BlueNode>::iterator nit = GammaH.begin();
              nit!=GammaH.end(); ++nit){
           PI[*nit] += gamma;
         }
@@ -1097 +1097 @@
     {
       Value retVal = 0;
       Edge e;
-      for( BlueIt n(*G); n != INVALID; ++n ) {
+      for( BlueNodeIt n(*G); n != INVALID; ++n ) {
         if( ( e = (*BM)[n] ) != INVALID ) retVal += (*c)[e];
       }
       return retVal;
@@ -1111 +1111 @@
     int matchingSize() const
     {
       int size = 0;
-      for ( BlueIt n(*G); n != INVALID; ++n) {
+      for ( BlueNodeIt n(*G); n != INVALID; ++n) {
         if ( (*BM)[n] != INVALID) {
           ++size;
         }
@@ -1139 +1139 @@
     /// \pre Either run() or start() must be called before using this function.
     Edge matching( const Node& node ) const
     {
-      return G->blue(node) ? (*BM)[node] : (*RM)[node];
+      return G->blue(node) ? (*BM)[G->asBlueNodeUnsafe(node)] : (*RM)[G->asRedNodeUnsafe(node)];
     }
 
     /// \brief Return a const reference to the blue matching map.
@@ -1170 +1170 @@
     {
       Edge e;
       if( G->blue(node) ){
-        if( ( e = (*BM)[node]) != INVALID ) return G->redNode(e);
+        if( ( e = (*BM)[G->asBlueNodeUnsafe(node)]) != INVALID ) return G->redNode(e);
         else return INVALID;
       }
       else {
-        if( ( e = (*RM)[node]) != INVALID ) return G->blueNode(e);
+        if( ( e = (*RM)[G->asRedNodeUnsafe(node)]) != INVALID ) return G->blueNode(e);
         else return INVALID;
       }
     }

lemon/lgf_reader.h

@@ -2953 +2953 @@
             throw FormatError(msg.str());
           }
 
-          e = _graph.addEdge(source, target);
+          // It is checked that source is red and
+          // target is blue, so this should be safe:
+          e = _graph.addEdge(_graph.asRedNodeUnsafe(source),
+                             _graph.asBlueNodeUnsafe(target));
           if (label_index != -1)
             _edge_index.insert(std::make_pair(tokens[label_index], e));
         } else {

lemon/elevator.h

@@ -2 +2 @@
  *
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2003-2009
+ * Copyright (C) 2003-2012
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
  *
@@ -108 +108 @@
     Elevator(const GR &graph,int max_level) :
       _g(graph),
       _max_level(max_level),
-      _item_num(_max_level),
+      _item_num(countItems<GR, Item>(graph)),
       _where(graph),
       _level(graph,0),
-      _items(_max_level),
+      _items(_item_num),
       _first(_max_level+2),
       _last_active(_max_level+2),
       _highest_active(-1) {}
@@ -526 +526 @@
     ///\param max_level The maximum allowed level.
     ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
     LinkedElevator(const GR& graph, int max_level)
-      : _graph(graph), _max_level(max_level), _item_num(_max_level),
+      : _graph(graph), _max_level(max_level),
+        _item_num(countItems<GR, Item>(graph)),
         _first(_max_level + 1), _last(_max_level + 1),
         _prev(graph), _next(graph),
         _highest_active(-1), _level(graph), _active(graph) {}

new file lemon/hopcroft_karp.h

@@ +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.
+ *
+ */
+
+#ifndef LEMON_HOPCROFT_KARP_H
+#define LEMON_HOPCROFT_KARP_H
+
+#include <lemon/core.h>
+#include <vector>
+#include <queue>
+#include <list>
+
+/// \ingroup matching
+/// \file
+/// \brief Hopcroft-Karp algorithm.
+namespace lemon {
+
+  /// \brief The Hopcroft-Karp bipartite matching algorithm
+  ///
+  /// Finds maximal matching in a given bipartite
+  /// graph using the Hopcroft-Karp algorithm,
+  /// having \f$O(e\sqrt{n})\f$ complexity.
+  template <typename BPG>
+  class HopcroftKarp {
+  public:
+    /// Type of the bipartite graph
+    typedef BPG BpGraph;
+    /// Type of the matching map
+    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
+  private:
+    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
+
+    const BpGraph& _bpg;
+    MatchingMap* _matching;
+    bool _local_matching;
+
+  protected:
+    void createStructures() {
+      if (!_matching) {
+        _matching = new MatchingMap(_bpg, INVALID);
+        _local_matching = true;
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_matching) {
+        delete _matching;
+      }
+    }
+
+    HopcroftKarp() {}
+
+  public:
+    /// \brief Constructor
+    ///
+    /// Constructs the class on the given bipartite graph.
+    HopcroftKarp(const BpGraph& bpg) : _bpg(bpg),
+                                       _matching(0),
+                                       _local_matching(false)
+    {}
+
+
+    /// \brief Destructor
+    ///
+    /// Destructor.
+    ~HopcroftKarp() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the matching map
+    ///
+    /// Sets the matching map. It should contain a valid matching,
+    /// for example the empty matching is fine.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
+    ///
+    /// \return <tt>(*this)</tt>
+    HopcroftKarp& matchingMap(MatchingMap& map) {
+      _matching = &map;
+      _local_matching = false;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the matching map
+    ///
+    /// Returns a const reference to the matching map, which contains
+    /// the matching edge for every node (and \c INVALID for the
+    /// unmatched nodes).
+    const MatchingMap& matchingMap() const {
+      return *_matching;
+    }
+
+
+    /// \brief Initializes the algorithm
+    ///
+    /// Allocates the matching map if necessary.
+    ///
+    /// \pre \ref init() is not called.
+    void init() {
+      createStructures();
+    }
+
+    /// \brief Smarter initialization of the matching.
+    ///
+    /// Allocates the matching map if necessary, and initializes
+    /// the algorithm with a trivially found matching: iterating
+    /// on every edge, take it in if both endnodes are unmatched.
+    ///
+    /// \return Size of the found matching.
+    ///
+    /// \note heuristicInit() replaces init() regarding the preconditions.
+    int heuristicInit() {
+      createStructures();
+      int size = 0;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        bool matched = false;
+        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
+          if ((*_matching)[_bpg.target(oi)] == INVALID) {
+            _matching->set(_bpg.source(oi), oi);
+            _matching->set(_bpg.target(oi), oi);
+            matched = true;
+            ++size;
+          }
+        }
+      }
+
+      return size;
+    }
+
+    /// \brief Initialize the matching from a map.
+    ///
+    /// Allocates the matching map if necessary, and
+    /// initializes the algorithm with a matching given as a bool edgemap,
+    /// in which an edge is true if it is in the matching.
+    ///
+    /// If the given matching is invalid, some edges will be left out.
+    ///
+    /// \return \c false if the given matching is invalid.
+    ///
+    /// \note matchingInit() replaces init() regarding the preconditions.
+    bool matchingInit(const BoolEdgeMap& matching) {
+      createStructures();
+      bool valid = true;
+      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
+        if (matching[it]) {
+          Node red = _bpg.redNode(it);
+          Node blue = _bpg.blueNode(it);
+          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
+            valid = false;
+          } else {
+            _matching->set(red, it);
+            _matching->set(blue, it);
+          }
+        }
+      }
+      return valid;
+    }
+
+    /// \brief Executes an augmenting phase
+    ///
+    /// Searches a maximal set of vertex disjoint shortest alternating paths.
+    /// Meaning:
+    ///  - alternating: connents an unmatched red- and an unmatched blue node,
+    ///    and exactly every second edge is in the current matching;
+    ///  - shortest: contain a minimal number of edges;
+    ///  - vertex disjoint: every vertex belong to at most one path;
+    ///  - maximal set: a set of path is maximal when it is not expandable
+    ///    further (and not when there does not exist a set with
+    ///    more vertex disjoint shortest alternating paths).
+    ///
+    /// After a maximal set is found, it applies the augmenting paths,
+    /// so edges of the matching are taken out, the others are put in
+    /// the matching.
+    ///
+    /// \return The length of the found augmenting paths, or -1 when none
+    /// found, which occurs if and only if the current matching is maximal.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    int augment() {
+      IntBlueNodeMap level(_bpg, -2);
+      std::vector<RedNode> act_red_nodes;
+
+      for (RedNodeIt r_it(_bpg); r_it != INVALID; ++r_it) {
+        if ((*_matching)[r_it] == INVALID) {
+          act_red_nodes.push_back(r_it);
+        }
+      }
+      // This vector will contain the end nodes of the possible
+      // augmenting paths.
+      std::vector<BlueNode> path_ends;
+
+      // Layer counter
+      int cur_level = -1;
+
+      // Starting from the unmatched red nodes search for unmatched
+      // blue nodes, using Bfs (but only on alternating paths).
+      // Blue nodes with level==-2 are unreached.
+      while (path_ends.empty()) {
+        cur_level += 2;
+        std::vector<RedNode> next_red_nodes;
+        for (typename std::vector<RedNode>::iterator r=act_red_nodes.begin();
+             r != act_red_nodes.end(); ++r) {
+          for (OutArcIt it(_bpg, *r); it!=INVALID; ++it) {
+            BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+            if (level[blue] != -2) continue;
+            level[blue] = cur_level;
+            if ((*_matching)[blue] == INVALID) {
+              path_ends.push_back(blue);
+            } else {
+              next_red_nodes
+                .push_back(_bpg.redNode((*_matching)[blue]));
+            }
+          }
+        }
+
+        if (path_ends.empty() && next_red_nodes.empty())
+          return -1;
+
+        act_red_nodes.swap(next_red_nodes);
+        next_red_nodes.clear();
+      }
+      // length of found augmenting paths
+      int max_level = cur_level;
+
+      // Stack for the Dfs
+      std::vector<OutArcIt> stack;
+      // Raise this flag when a shortest augmenting path is found.
+      bool path_found;
+      // Searching backward, starting from the previously found
+      // blue nodes, we build vertex disjoint alternating paths.
+      // Blue nodes with level!=-2 are unreached.
+      typename std::vector<BlueNode>::iterator end = path_ends.begin();
+
+      while (end != path_ends.end()) {
+        stack.clear();
+        stack.push_back(OutArcIt(_bpg, *end));
+        cur_level = max_level;
+        path_found = false;
+
+        while (!stack.empty() && !path_found) {
+          OutArcIt &o = stack.back();
+          // Search an arc that goes one level higher.
+          // If none found, retreat.
+          while (o != INVALID &&
+                 (*_matching)[_bpg.target(o)] != INVALID &&
+                 level[_bpg.blueNode((*_matching)[_bpg.target(o)])] !=
+                 cur_level - 2) {
+            ++o;
+          }
+          if (o==INVALID) {
+            stack.pop_back();
+            cur_level += 2;
+          } else if ((*_matching)[_bpg.target(o)] == INVALID) {
+            path_found = true;
+          } else {
+            BlueNode b = _bpg.blueNode((*_matching)[_bpg.target(o)]);
+            stack.push_back(OutArcIt(_bpg, b));
+            level[b] = -2;        // Do not visit this node again
+            cur_level -= 2;
+          }
+        }
+        if (path_found) {
+          OutArcIt o;
+          while (!stack.empty()) {
+            o = stack.back();
+            _matching->set(_bpg.source(o), o);
+            _matching->set(_bpg.target(o), o);
+            stack.pop_back();
+          }
+        }
+        ++end;
+      }
+      return max_level;
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// It runs augmenting phases until the optimal solution is reached.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (augment() > 0) {}
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// hk.run() is just a shorthand for:
+    ///
+    ///\code
+    /// hk.heuristicInit();
+    /// hk.start();
+    ///\endcode
+    void run() {
+      heuristicInit();
+      start();
+    }
+
+    /// \brief Size of the matching
+    ///
+    /// Returns the size of the current matching.
+    /// Function runs in \f$O(n)\f$ time.
+    int matchingSize() const {
+      int size = 0;
+      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
+        if ((*_matching)[it] != INVALID) ++size;
+      }
+      return size;
+    }
+
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the current
+    /// matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpg.redNode(edge)];
+    }
+
+    /// \brief Return the matching edge incident to the given node.
+    ///
+    /// This function returns the matching edge incident to the
+    /// given node in the current matching or \c INVALID if the node is
+    /// not covered by the matching.
+    Edge matching(const Node& n) const {
+      return (*_matching)[n];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the current
+    /// matching or \c INVALID if the node is not covered by the matching.
+    Node mate(const Node& n) const {
+      return (*_matching)[n] != INVALID ?
+        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
+    }
+
+    /// \brief Query a vertex cover.
+    ///
+    /// This function finds a vertex cover based on the current matching,
+    /// and in \c cover sets node in the cover \c true and \c false
+    /// otherwise. In case of maximal matching, this will be a minimal
+    /// vertex cover.
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return The number of nodes in the cover.
+    int cover(BoolNodeMap &cover) const {
+      int count = 0;
+      std::queue<RedNode> q;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
+        cover[b] = false;
+
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          cover[r] = false;
+          q.push(r);
+        } else {
+          cover[r] = true;
+          ++count;
+        }
+      }
+
+      while (!q.empty()) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (cover[blue]) continue;
+
+          cover[blue] = true;
+          ++count;
+
+          if ((*_matching)[blue] != INVALID) {
+            RedNode nr = _bpg.redNode((*_matching)[blue]);
+            q.push(nr);
+            cover[nr] = false;
+            --count;
+          }
+        }
+        q.pop();
+      }
+
+      return count;
+    }
+
+    /// \brief Query a barrier of red nodes.
+    ///
+    /// This function tries to create a barrier of red nodes, which should:
+    ///  - contain every unmatched red nodes,
+    ///  - contain exactly that many matched red nodes as much blue nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool redBarrier(BoolRedNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<RedNode> q;
+
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          barrier[r] = true;
+          q.push(r);
+        } else {
+          barrier[r] = false;
+        }
+      }
+
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+
+          if (exist &= ((*_matching)[blue] != INVALID)) {
+            RedNode r = _bpg.redNode((*_matching)[blue]);
+            if (!barrier[r]) {
+              q.push(r);
+              barrier[r] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+
+      return exist;
+    }
+
+    /// \brief Query a barrier of blue nodes.
+    ///
+    /// This function tries to create a barrier of blue nodes, which should:
+    ///  - contain every unmatched blue nodes,
+    ///  - contain exactly that many matched blue nodes as much red nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool blueBarrier(BoolBlueNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<BlueNode> q;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        if ((*_matching)[b] == INVALID) {
+          barrier[b] = true;
+          q.push(b);
+        } else {
+          barrier[b] = false;
+        }
+      }
+
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
+
+          if (exist &= ((*_matching)[r] != INVALID)) {
+            BlueNode b = _bpg.blueNode((*_matching)[r]);
+            if (!barrier[b]) {
+              q.push(b);
+              barrier[b] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+
+      return exist;
+    }
+  };
+
+}
+#endif

test/CMakeLists.txt

@@ -34 +34 @@
   graph_utils_test
   hao_orlin_test
   heap_test
+  hopcroft_karp_test
   kruskal_test
   lgf_test
   maps_test

test/Makefile.am

@@ -32 +32 @@
         test/graph_utils_test \
         test/hao_orlin_test \
         test/heap_test \
+        test/hopcroft_karp_test \
         test/kruskal_test \
         test/lgf_test \
         test/maps_test \
@@ -87 +88 @@
 test_graph_utils_test_SOURCES = test/graph_utils_test.cc
 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
 test_heap_test_SOURCES = test/heap_test.cc
+test_hopcroft_karp_test_SOURCES = test/hopcroft_karp_test.cc
 test_kruskal_test_SOURCES = test/kruskal_test.cc
 test_lgf_test_SOURCES = test/lgf_test.cc
 test_lp_test_SOURCES = test/lp_test.cc

new file test/hopcroft_karp_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/hopcroft_karp.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 = 3;
+const string lgf[lgfn] = {
+  // Empty graph
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+
+  // No red nodes
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "1\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+
+  // No blue nodes
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+};
+
+const string u_lgf =
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "@blue_nodes\n"
+  "label\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "label\n"
+  "1 4 1\n"
+  "1 5 2\n"
+  "2 5 3\n"
+  "2 6 4\n"
+  "2 7 5\n"
+  "3 7 6\n"
+;
+
+
+void checkHopcroftKarpCompile() {
+  typedef concepts::BpGraph BpGraph;
+  BPGRAPH_TYPEDEFS(BpGraph);
+  typedef HopcroftKarp<BpGraph> HK;
+
+  BpGraph graph;
+  RedNode red;
+  BlueNode blue;
+  Node node;
+  Edge edge;
+  BoolEdgeMap init_matching(graph);
+  HK hk_test(graph);
+  const HK& const_hk_test = hk_test;
+  HK::MatchingMap matching_map(graph, INVALID);
+  BoolNodeMap cov(graph);
+  int cover_size;
+  BoolRedNodeMap rb(graph);
+  BoolBlueNodeMap bb(graph);
+
+  hk_test.matchingMap(matching_map);
+  hk_test.init();
+  hk_test.heuristicInit();
+  hk_test.matchingInit(init_matching);
+
+  hk_test.start();
+  hk_test.augment();
+  hk_test.run();
+
+  const_hk_test.matchingSize();
+  const_hk_test.matching(node);
+  const_hk_test.matching(edge);
+  const_hk_test.mate(red);
+  const_hk_test.mate(blue);
+  const_hk_test.mate(node);
+  const HK::MatchingMap& max_matching = const_hk_test.matchingMap();
+  edge = max_matching[node];
+  cover_size = const_hk_test.cover(cov);
+  const_hk_test.redBarrier(rb);
+  const_hk_test.blueBarrier(bb);
+}
+
+typedef SmartBpGraph BpGraph;
+typedef HopcroftKarp<BpGraph> HK;
+BPGRAPH_TYPEDEFS(BpGraph);
+
+void checkMatchingValidity(const BpGraph& bpg, const HK& hk) {
+  BoolBlueNodeMap matched(bpg, false);
+  for (RedNodeIt it(bpg); it != INVALID; ++it) {
+    if (hk.matching(it) != INVALID) {
+      check(hk.mate(hk.mate(it)) == it, "Wrong matching");
+      matched.set(bpg.asBlueNode(hk.mate(it)), true);
+    }
+  }
+  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
+    // != is used as logical xor here
+    check(matched[it] != (hk.matching(it) == INVALID), "Wrong matching");
+  }
+
+}
+
+void checkCover(const BpGraph& bpg, const HK& hk) {
+  int cov_size;
+  BoolNodeMap cov(bpg);
+  cov_size = hk.cover(cov);
+  BoolEdgeMap covered(bpg, false);
+  for (NodeIt n(bpg); n!=INVALID; ++n) {
+    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
+      covered[e] = true;
+    }
+  }
+  bool all = true;
+  for (EdgeIt e(bpg); e!=INVALID; ++e) {
+    all&=covered[e];
+  }
+
+  check(all, "Invalid cover");
+  check(hk.matchingSize() == cov_size, "Suboptimal matching.");
+}
+
+void checkBarriers(const BpGraph& bpg, const HK& hk) {
+  BoolRedNodeMap rb(bpg);
+  BoolBlueNodeMap bb(bpg);
+
+  check(hk.redBarrier(rb), "Suboptimal matching.");
+  check(hk.blueBarrier(bb), "Suboptimal matching.");
+
+  BoolNodeMap reached(bpg, false);
+
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (!rb[r]) continue;
+    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (!bb[b]) continue;
+    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+
+  int ur = 0, ub = 0, rc = 0, bc = 0;
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (rb[r]) ++rc;
+    if (hk.mate(r) == INVALID) ++ur;
+    if (reached[r]) --bc;
+  }
+
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (bb[b]) ++bc;
+    if (hk.mate(b) == INVALID) ++ub;
+    if (reached[b]) --rc;
+  }
+
+  check(rc == ur, "Red barrier is wrong.");
+  check(bc == ub, "Blue barrier is wrong.");
+}
+
+void checkMatching(const BpGraph& bpg, const HK& hk) {
+  checkMatchingValidity(bpg, hk);
+  checkBarriers(bpg, hk);
+  checkCover(bpg, hk);
+}
+
+
+int main() {
+  // Checking hard wired extremal graphs
+  for (int i=0; i<lgfn; ++i) {
+    BpGraph bpg;
+    istringstream lgfs(lgf[i]);
+    bpGraphReader(bpg, lgfs).run();
+    HK hk(bpg);
+    hk.run();
+    checkMatching(bpg, hk);
+  }
+
+  // Checking some use cases
+  BpGraph bpg;
+  istringstream u_lgfs(u_lgf);
+  bpGraphReader(bpg, u_lgfs).run();
+  {
+    HK hk(bpg);
+    hk.init();
+    hk.augment();
+    hk.start();
+    checkMatching(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    hk.heuristicInit();
+    hk.augment();
+    hk.start();
+    checkMatching(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    HK::MatchingMap m(bpg, INVALID);
+    hk.matchingMap(m);
+    hk.run();
+    checkMatching(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    BpGraph::EdgeMap<bool> init_matching(bpg, false);
+    BpGraph::Edge e;
+    bpg.first(e);
+    init_matching[e] = true;
+    check(hk.matchingInit(init_matching), "Wrong result from matchingInit()");
+    hk.start();
+    checkMatching(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    HK::MatchingMap m(bpg, INVALID);
+    hk.matchingMap(m);
+    BpGraph::EdgeMap<bool> init_matching(bpg, true);
+    check(!hk.matchingInit(init_matching), "Wrong result from matchingInit()");
+    hk.start();
+    checkMatching(bpg, hk);
+  }
+  // Checking some random generated graphs
+  const int random_test_n = 10;
+  const int max_red_n = 1000;
+  const int min_red_n = 10;
+  const int max_blue_n = 1000;
+  const int min_blue_n = 10;
+  for (int i=0; i<random_test_n; ++i) {
+    int red_n = rnd.integer(min_red_n, max_red_n);
+    int blue_n = rnd.integer(min_blue_n, max_blue_n);
+    double prob = rnd();
+    SmartBpGraph bpg;
+    for (int i=0; i<red_n; ++i) {
+      bpg.addRedNode();
+    }
+    for (int i=0; i<blue_n; ++i) {
+      bpg.addBlueNode();
+    }
+    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (rnd() < prob/10.) {
+          bpg.addEdge(r,b);
+        }
+      }
+    }
+
+    HK hk(bpg);
+    hk.run();
+    checkMatching(bpg, hk);
+  }
+
+
+  return 0;
+}
+
+

new file lemon/preflow_bp_matching.h

@@ +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.
+ *
+ */
+
+#ifndef LEMON_PREFLOW_BP_MATCHING_H
+#define LEMON_PREFLOW_BP_MATCHING_H
+
+#include <lemon/core.h>
+#include <lemon/elevator.h>
+#include <vector>
+#include <list>
+#include <queue>
+
+/// \ingroup matching
+/// \file
+/// \brief A preflow based bipartite matching algorithm.
+namespace lemon {
+
+  /// \brief A preflow based bipartite matching algorithm
+  ///
+  /// Finds maximum cardinality matching in a given bipartite
+  /// graph using specialized preflow algorithm.
+  template <typename BPG>
+  class PreflowBpMatching {
+  public:
+    /// Type of the bipartite graph
+    typedef BPG BpGraph;
+    /// Type of the matching map
+    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
+    /// Type of the elevator.
+    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
+
+  private:
+    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
+
+  protected:
+    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
+    typedef std::list<typename BPG::Arc> CurrentArcSet;
+    typedef typename
+    std::list<typename BPG::Arc>::iterator CurrentArcIterator;
+
+    const BpGraph& _bpg;
+    MatchingMap* _matching;
+    bool _local_matching;
+    Elevator_t* _elevator;
+    bool _local_elevator;
+    typename BPG::template BlueNodeMap<CurrentArcSet> _cas;
+    typename BPG::template BlueNodeMap<CurrentArcIterator> _ca;
+    Preflow _flow;
+
+    void createStructures() {
+      if (!_matching) {
+        _matching = new MatchingMap(_bpg, INVALID);
+        _local_matching = true;
+      }
+      if (!_elevator) {
+        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
+                                                  countBlueNodes(_bpg)));
+        _local_elevator = true;
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_matching) {
+        delete _matching;
+      }
+      if (_local_elevator) {
+        delete _elevator;
+      }
+    }
+
+    /// \brief Initialize the elevator
+    ///
+    /// When this function is called, the elevator is initialized with
+    /// the exact distance labels corresponding the current matching.
+    void initElevator() {
+
+      // Set the initial flow:
+      // choose arbitrary outgoing arc for every unmatched red node.
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        Arc a = (*_matching)[r] != INVALID ?
+          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
+        if (a != INVALID) {
+          _flow[r] = a;
+          _cas[_bpg.asBlueNodeUnsafe(_bpg.target(a))].push_front(a);
+        }
+      }
+      // Now we start a backward Bfs from the unmatched blue nodes in the
+      // residual graph, to supply the exact distance labels to the elevator.
+      std::vector<BlueNode> act_blue_nodes, next_blue_nodes;
+      BoolBlueNodeMap reached(_bpg, false);
+
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (_cas[b].empty()) {
+          act_blue_nodes.push_back(b);
+          reached.set(b, true);
+        }
+      }
+
+      _elevator->initStart();
+      while (!act_blue_nodes.empty()) {
+        for (typename std::vector<BlueNode>::iterator b_it =
+               act_blue_nodes.begin();
+             b_it != act_blue_nodes.end(); ++b_it) {
+
+          _elevator->initAddItem(*b_it);
+          for (OutArcIt oi(_bpg, *b_it); oi!=INVALID; ++oi) {
+            RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(oi));
+            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
+            if (b != *b_it && !reached[b]) {
+              next_blue_nodes.push_back(b);
+              reached.set(b, true);
+            }
+          }
+        }
+        _elevator->initNewLevel();
+        next_blue_nodes.swap(act_blue_nodes);
+        next_blue_nodes.clear();
+      }
+      _elevator->initFinish();
+
+      // Activate blue nodes with excess, i. e. those
+      // which get flow from more than one red node.
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (_cas[b].size() > 1 && (*_elevator)[b] < _elevator->maxLevel()) {
+          _elevator->activate(b);
+        }
+
+        _ca[b] = _cas[b].begin();
+      }
+    }
+
+    PreflowBpMatching() {}
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructs the class on the given bipartite graph.
+    PreflowBpMatching(const BpGraph& bpg) : _bpg(bpg),
+                                            _matching(0),
+                                            _local_matching(false),
+                                            _elevator(0),
+                                            _local_elevator(false),
+                                            _cas(bpg),
+                                            _ca(bpg),
+                                            _flow(bpg, INVALID)
+    {}
+
+    /// \brief Destructor
+    ///
+    /// Destructor.
+    ~PreflowBpMatching() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the matching map
+    ///
+    /// Sets the matching map. It should contain a valid matching,
+    /// for example the empty matching is fine.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
+    ///
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching& matchingMap(MatchingMap& map) {
+      _matching = &map;
+      _local_matching = false;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the matching map.
+    ///
+    /// Returns a const reference to the matching map, which contains
+    /// the matching edge for each node (or INVALID for unmatched nodes).
+    const MatchingMap& matchingMap() const {
+      return *_matching;
+    }
+
+    /// \brief Sets the elevator used by algorithm.
+    ///
+    /// Sets the elevator.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated elevator,
+    /// of course.
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching& elevator(Elevator_t& elevator) {
+      _elevator = &elevator;
+      _local_elevator = false;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the elevator.
+    ///
+    /// Returns a const reference to the elevator, which keeps track
+    /// of the distance labels of nodes.
+    const Elevator_t& elevator() const {
+      return *_elevator;
+    }
+
+    /// \brief Initializes the algorithm
+    ///
+    /// Allocates the elevator and the matching map if necessary,
+    /// and initializes the elevator.
+    ///
+    /// \pre \ref init() is not called.
+    void init() {
+      createStructures();
+      initElevator();
+    }
+
+    /// \brief Smarter initialization of the matching.
+    ///
+    /// Allocates the matching map if necessary, and initializes
+    /// the algorithm with a trivially found matching: iterating
+    /// on every edge, take it in if both endnodes are unmatched.
+    ///
+    /// \return Size of the so found matching.
+    ///
+    /// \note heuristicInit() replaces init() regarding the preconditions.
+    int heuristicInit() {
+      createStructures();
+      int size = 0;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        bool matched = false;
+        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
+          if ((*_matching)[_bpg.target(oi)] == INVALID) {
+            _matching->set(_bpg.source(oi), oi);
+            _matching->set(_bpg.target(oi), oi);
+            matched = true;
+            ++size;
+          }
+        }
+      }
+
+      initElevator();
+      return size;
+    }
+
+    /// \brief Initialize the matching from a map.
+    ///
+    /// Allocates the matching map if necessary, and
+    /// initializes the algorithm with a matching given as a bool edgemap,
+    /// in which an edge is true if it is in the matching.
+    /// Also initializes the elevator.
+    ///
+    /// \return \c true if the matching is valid.
+    ///
+    /// \note matchingInit() replaces init() regarding the preconditions.
+    bool matchingInit(const BoolEdgeMap& matching) {
+      createStructures();
+
+      bool valid = true;
+      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
+        if (matching[it]) {
+          Node red = _bpg.redNode(it);
+          Node blue = _bpg.blueNode(it);
+          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
+            valid = false;
+          } else {
+            _matching->set(red, it);
+            _matching->set(blue, it);
+          }
+        }
+      }
+
+      initElevator();
+      return valid;
+    }
+
+
+    /// \brief Perform a push/relabel operation on a node
+    ///
+    /// Perform a push/relabel operation on a node.
+    ///
+    /// \pre \ref init() is called, and the node is active.
+    void push_relabel(const BlueNode& chosen) {
+      Arc new_flow = INVALID;
+      while (_ca[chosen] != _cas[chosen].end() && new_flow == INVALID) {
+        OutArcIt r_oi(_bpg, _bpg.source(*_ca[chosen]));
+        while (r_oi != INVALID && new_flow == INVALID) {
+          if ((*_elevator)[_bpg.asBlueNodeUnsafe(_bpg.target(r_oi))] ==
+              (*_elevator)[chosen] - 1) {
+            new_flow = r_oi;
+          }
+          ++r_oi;
+        }
+        if (new_flow == INVALID) ++_ca[chosen];
+      }
+
+      // If an incident blue node is found on one level lower,
+      // then perform a push, i. e. flip the flow from the
+      // current arc to 'new_flow'.
+      if (new_flow != INVALID) {
+        _ca[chosen] = _cas[chosen].erase(_ca[chosen]);
+        if (_cas[chosen].size() == 1) {
+          _elevator->deactivate(chosen);
+        }
+        _flow[_bpg.asRedNodeUnsafe(_bpg.source(new_flow))] = new_flow;
+        BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(new_flow));
+        _cas[b].push_front(new_flow);
+        if (_cas[b].size() == 2) {
+          _elevator->activate(b);
+        }
+      }
+
+      // If none found, then relabel.
+      else {
+        int min = _elevator->maxLevel() + 1;
+        for (CurrentArcIterator it1 = _cas[chosen].begin();
+             it1 != _cas[chosen].end();
+             ++it1) {
+          for (OutArcIt it2(_bpg, _bpg.asRedNodeUnsafe(_bpg.source(*it1)));
+               it2 != INVALID;
+               ++it2) {
+            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(it2));
+            if (b != chosen && min > (*_elevator)[b]) {
+              min = (*_elevator)[b];
+            }
+          }
+        }
+
+        int chosen_level = (*_elevator)[chosen];
+        if (min >= _elevator->maxLevel() - 1) {
+          _elevator->lift(chosen, _elevator->maxLevel());
+          _elevator->deactivate(chosen);
+        } else {
+          _elevator->lift(chosen, min+1);
+        }
+        if (_elevator->emptyLevel(chosen_level)) {
+          _elevator->liftToTop(chosen_level);
+        }
+
+        // Restore the current arc iterator.
+        _ca[chosen] = _cas[chosen].begin();
+      }
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// Executes push/relabel operation on the highest active node
+    /// until there is. After start(), writeMatching() should be called to
+    /// remove unnecessary flow and choose a valid maximum matching.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (_elevator->highestActive() != INVALID)
+        push_relabel(_elevator->highestActive());
+    }
+
+    /// \brief Deduces and writes a matching from the preflow
+    ///
+    /// The underlying preflow algorithm does not use the matching map,
+    /// so before using the matching query functions ( \ref mate(), \ref
+    /// matchingSize(), etc.), you should explicitly call this function
+    /// to write a matching based on the current preflow to the matching
+    /// map. When the preflow results an ambiguous matching (a node has
+    /// incoming flow on more than one arc), it chooses arbitrary one.
+    void writeMatching() {
+      for (NodeIt it(_bpg); it != INVALID; ++it) {
+        _matching->set(it, INVALID);
+      }
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (!_cas[b].empty()) {
+          Arc a = _cas[b].front();
+          _matching->set(_bpg.source(a), a);
+          _matching->set(_bpg.target(a), a);
+        }
+      }
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// bp1.run() is just a shorthand for:
+    ///
+    ///\code
+    /// bp1.heuristicInit();
+    /// bp1.start();
+    /// bp1.writeMatching();
+    ///\endcode
+    void run() {
+      heuristicInit();
+      start();
+      writeMatching();
+    }
+
+    /// \brief Size of the matching
+    ///
+    /// Returns the size of the current matching.
+    /// Function runs in \f$O(n)\f$ time.
+    int matchingSize() const {
+      int size = 0;
+      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
+        if ((*_matching)[it] != INVALID) ++size;
+      }
+      return size;
+    }
+
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the current
+    /// matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpg.redNode(edge)];
+    }
+
+    /// \brief Return the matching edge incident to the given node.
+    ///
+    /// This function returns the matching edge incident to the
+    /// given node in the current matching or \c INVALID if the node is
+    /// not covered by the matching.
+    Edge matching(const Node& n) const {
+      return (*_matching)[n];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the current
+    /// matching or \c INVALID if the node is not covered by the matching.
+    Node mate(const Node& n) const {
+      return (*_matching)[n] != INVALID ?
+        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
+    }
+
+    /// \brief Query a vertex cover.
+    ///
+    /// This function finds a vertex cover based on the current matching,
+    /// and in \c cover sets node in the cover \c true and \c false
+    /// otherwise. In case of maximal matching, this will be a minimal
+    /// vertex cover.
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return The number of nodes in the cover.
+    int cover(BoolNodeMap &cover) const {
+      int count = 0;
+      std::queue<RedNode> q;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
+        cover[b] = false;
+
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          cover[r] = false;
+          q.push(r);
+        } else {
+          cover[r] = true;
+          ++count;
+        }
+      }
+
+      while (!q.empty()) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (cover[blue]) continue;
+
+          cover[blue] = true;
+          ++count;
+
+          if ((*_matching)[blue] != INVALID) {
+            RedNode nr = _bpg.redNode((*_matching)[blue]);
+            q.push(nr);
+            cover[nr] = false;
+            --count;
+          }
+        }
+        q.pop();
+      }
+
+      return count;
+    }
+
+    /// \brief Query a barrier of red nodes.
+    ///
+    /// This function tries to create a barrier of red nodes, which should:
+    ///  - contain every unmatched red nodes,
+    ///  - contain exactly that many matched red nodes as much blue nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool redBarrier(BoolRedNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<RedNode> q;
+
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          barrier[r] = true;
+          q.push(r);
+        } else {
+          barrier[r] = false;
+        }
+      }
+
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+
+          if (exist &= ((*_matching)[blue] != INVALID)) {
+            RedNode r = _bpg.redNode((*_matching)[blue]);
+            if (!barrier[r]) {
+              q.push(r);
+              barrier[r] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+
+      return exist;
+    }
+
+    /// \brief Query a barrier of blue nodes.
+    ///
+    /// This function tries to create a barrier of blue nodes, which should:
+    ///  - contain every unmatched blue nodes,
+    ///  - contain exactly that many matched blue nodes as much red nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool blueBarrier(BoolBlueNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<BlueNode> q;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        if ((*_matching)[b] == INVALID) {
+          barrier[b] = true;
+          q.push(b);
+        } else {
+          barrier[b] = false;
+        }
+      }
+
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
+
+          if (exist &= ((*_matching)[r] != INVALID)) {
+            BlueNode b = _bpg.blueNode((*_matching)[r]);
+            if (!barrier[b]) {
+              q.push(b);
+              barrier[b] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+
+      return exist;
+    }
+  };
+
+}
+#endif

test/CMakeLists.txt

@@ -48 +48 @@
   path_test
   planarity_test
   preflow_test
+  preflow_bp_matching_test
   radix_sort_test
   random_test
   suurballe_test

test/Makefile.am

@@ -46 +46 @@
         test/path_test \
         test/planarity_test \
         test/preflow_test \
+        test/preflow_bp_matching_test \
         test/radix_sort_test \
         test/random_test \
         test/suurballe_test \
@@ -104 +105 @@
 test_path_test_SOURCES = test/path_test.cc
 test_planarity_test_SOURCES = test/planarity_test.cc
 test_preflow_test_SOURCES = test/preflow_test.cc
+test_preflow_bp_matching_test_SOURCES = test/preflow_bp_matching_test.cc
 test_radix_sort_test_SOURCES = test/radix_sort_test.cc
 test_suurballe_test_SOURCES = test/suurballe_test.cc
 test_random_test_SOURCES = test/random_test.cc

new file test/preflow_bp_matching_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/preflow_bp_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 = 3;
+const string lgf[lgfn] = {
+  // Empty graph
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+
+  // No red nodes
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "1\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+
+  // No blue nodes
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+};
+
+const string u_lgf =
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "@blue_nodes\n"
+  "label\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "label\n"
+  "1 4 1\n"
+  "1 5 2\n"
+  "2 5 3\n"
+  "2 6 4\n"
+  "2 7 5\n"
+  "3 7 6\n"
+;
+
+
+
+void checkPreflowBpMatchingCompile() {
+  typedef concepts::BpGraph BpGraph;
+  BPGRAPH_TYPEDEFS(BpGraph);
+  typedef PreflowBpMatching<BpGraph> PBM;
+
+  BpGraph graph;
+  RedNode red;
+  BlueNode blue;
+  Node node;
+  Edge edge;
+  BoolEdgeMap init_matching(graph);
+  PBM pbm_test(graph);
+  const PBM& const_pbm_test = pbm_test;
+  PBM::MatchingMap matching_map(graph, INVALID);
+  PBM::Elevator_t elevator(graph);
+  BpGraph::NodeMap<bool> cov(graph);
+  int cover_size;
+  BoolRedNodeMap rb(graph);
+  BoolBlueNodeMap bb(graph);
+
+  pbm_test.matchingMap(matching_map);
+  pbm_test.elevator(elevator);
+  pbm_test.init();
+  pbm_test.heuristicInit();
+  pbm_test.matchingInit(init_matching);
+
+  pbm_test.start();
+  pbm_test.push_relabel(blue);
+  pbm_test.run();
+  pbm_test.writeMatching();
+
+  const_pbm_test.matchingSize();
+  const_pbm_test.matching(node);
+  const_pbm_test.matching(edge);
+  const_pbm_test.mate(red);
+  const_pbm_test.mate(blue);
+  const_pbm_test.mate(node);
+  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
+  edge = max_matching[node];
+  const PBM::Elevator_t& elev = const_pbm_test.elevator();
+  blue = elev.highestActive();
+  cover_size = const_pbm_test.cover(cov);
+  const_pbm_test.redBarrier(rb);
+  const_pbm_test.blueBarrier(bb);
+}
+
+typedef SmartBpGraph BpGraph;
+typedef PreflowBpMatching<BpGraph> PBM;
+BPGRAPH_TYPEDEFS(BpGraph);
+
+void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
+  BoolBlueNodeMap matched(bpg, false);
+  for (RedNodeIt it(bpg); it != INVALID; ++it) {
+    if (pbm.matching(it) != INVALID) {
+      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
+      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
+    }
+  }
+  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
+    // != is used as logical xor here
+    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
+  }
+}
+
+void checkCover(const BpGraph& bpg, const PBM& pbm) {
+  int cov_size;
+  BoolNodeMap cov(bpg);
+  cov_size = pbm.cover(cov);
+  BoolEdgeMap covered(bpg, false);
+  for (NodeIt n(bpg); n!=INVALID; ++n) {
+    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
+      covered[e] = true;
+    }
+  }
+  bool all = true;
+  for (EdgeIt e(bpg); e!=INVALID; ++e) {
+    all&=covered[e];
+  }
+
+  check(all, "Invalid cover");
+  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
+}
+
+void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
+  BoolRedNodeMap rb(bpg);
+  BoolBlueNodeMap bb(bpg);
+
+  check(pbm.redBarrier(rb), "Suboptimal matching.");
+  check(pbm.blueBarrier(bb), "Suboptimal matching.");
+
+  BoolNodeMap reached(bpg, false);
+
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (!rb[r]) continue;
+    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (!bb[b]) continue;
+    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+
+  int ur = 0, ub = 0, rc = 0, bc = 0;
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (rb[r]) ++rc;
+    if (pbm.mate(r) == INVALID) ++ur;
+    if (reached[r]) --bc;
+  }
+
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (bb[b]) ++bc;
+    if (pbm.mate(b) == INVALID) ++ub;
+    if (reached[b]) --rc;
+  }
+
+  check(rc == ur, "Red barrier is wrong.");
+  check(bc == ub, "Blue barrier is wrong.");
+}
+
+void checkMatching(const BpGraph& bpg, const PBM& pbm) {
+  checkMatchingValidity(bpg, pbm);
+  checkBarriers(bpg, pbm);
+  checkCover(bpg, pbm);
+}
+
+
+int main() {
+  // Checking hard wired extremal graphs
+  for (int i=0; i<lgfn; ++i) {
+    BpGraph bpg;
+    istringstream lgfs(lgf[i]);
+    bpGraphReader(bpg, lgfs).run();
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+
+  // Checking some use cases
+  BpGraph bpg;
+  istringstream u_lgfs(u_lgf);
+  bpGraphReader(bpg, u_lgfs).run();
+
+  {
+    PBM pbm(bpg);
+    pbm.init();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.push_relabel(b);
+    }
+    pbm.start();
+    pbm.writeMatching();
+    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
+    checkMatching(bpg, pbm);
+  }
+
+  {
+    PBM pbm(bpg);
+    pbm.heuristicInit();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.push_relabel(b);
+    }
+    pbm.start();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    BpGraph::EdgeMap<bool> init_matching(bpg, false);
+    BpGraph::Edge e;
+    bpg.first(e);
+    init_matching[e] = true;
+    check(pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    BpGraph::EdgeMap<bool> init_matching(bpg, true);
+    check(!pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+
+  // Checking some random generated graphs
+  const int random_test_n = 10;
+  const int max_red_n = 1000;
+  const int min_red_n = 10;
+  const int max_blue_n = 1000;
+  const int min_blue_n = 10;
+  for (int i=0; i<random_test_n; ++i) {
+    int red_n = rnd.integer(min_red_n, max_red_n);
+    int blue_n = rnd.integer(min_blue_n, max_blue_n);
+    double prob = rnd();
+
+    SmartBpGraph bpg;
+    for (int j=0; j<red_n; ++j) {
+      bpg.addRedNode();
+    }
+    for (int j=0; j<blue_n; ++j) {
+      bpg.addBlueNode();
+    }
+    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (rnd() < prob/10.) {
+          bpg.addEdge(r,b);
+        }
+      }
+    }
+
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+
+    // do not always use the highest active node
+    PBM pbm2(bpg);
+    PBM::Elevator_t elev(bpg);
+    pbm2.elevator(elev);
+    pbm2.init();
+    while (elev.highestActive() != INVALID) {
+      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (elev.active(b)) pbm2.push_relabel(b);
+      }
+    }
+    pbm2.writeMatching();
+    checkMatching(bpg, pbm2);
+  }
+
+
+  return 0;
+}
+
+

new file lemon/preflow_bp_matching_2.h

@@ +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.
+ *
+ */
+
+#ifndef LEMON_PREFLOW_BP_MATCHING_2_H
+#define LEMON_PREFLOW_BP_MATCHING_2_H
+
+#include <lemon/core.h>
+#include <lemon/elevator.h>
+#include <vector>
+#include <queue>
+
+/// \ingroup matching
+/// \file
+/// \brief A preflow based bipartite matching algorithm.
+namespace lemon {
+
+  /// \brief A preflow based bipartite matching algorithm
+  ///
+  /// Finds maximum cardinality matching in a given bipartite
+  /// graph using specialized preflow algorithm.
+  template <typename BPG>
+  class PreflowBpMatching2 {
+  public:
+    /// Type of the bipartite graph
+    typedef BPG BpGraph;
+    /// Type of the matching map
+    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
+    /// Type of the elevator.
+    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
+
+  private:
+    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
+
+  protected:
+    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
+    typedef typename BPG::template BlueNodeMap<int> Excess;
+    typedef std::vector<typename BPG::BlueNode> BlueNodes;
+    typedef std::vector<typename BPG::Arc> BlueFlow;
+    typedef typename BPG::template BlueNodeMap<BlueFlow> IncomingFlows;
+
+    const BpGraph& _bpg;
+    MatchingMap* _matching;
+    bool _local_matching;
+    Elevator_t* _elevator;
+    bool _local_elevator;
+    Preflow _flow;
+    Excess _excess;
+
+    void createStructures() {
+      if (!_matching) {
+        _matching = new MatchingMap(_bpg, INVALID);
+        _local_matching = true;
+      }
+      if (!_elevator) {
+        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
+                                                  countBlueNodes(_bpg)));
+        _local_elevator = true;
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_matching) {
+        delete _matching;
+      }
+      if (_local_elevator) {
+        delete _elevator;
+      }
+    }
+
+    /// \brief Initialize the elevator
+    ///
+    /// When this function is called, the elevator is initialized with
+    /// the exact distance labels corresponding the current matching.
+    void initElevator() {
+      BlueNodes act_blue_nodes, next_blue_nodes;
+      BoolBlueNodeMap reached(_bpg, false);
+      IncomingFlows incoming(_bpg);
+
+      // Set the initial flow:
+      // choose arbitrary outgoing arc for every unmatched red node.
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        Arc a = (*_matching)[r] != INVALID ?
+          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
+        if (a != INVALID) {
+          _flow[r] = a;
+          BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(a));
+          incoming[b].push_back(a);
+          if (++_excess[b] == 2) {
+            act_blue_nodes.push_back(b);
+            reached.set(b, true);
+          }
+        }
+      }
+      // Now we start a Bfs from the blue nodes with excess
+      // (those which get flow from more than one red node) in the
+      // residual graph, to supply the exact distance labels to the elevator.
+      _elevator->initStart();
+      while (!act_blue_nodes.empty()) {
+        for (typename BlueNodes::iterator b_it = act_blue_nodes.begin();
+             b_it != act_blue_nodes.end(); ++b_it) {
+
+          _elevator->initAddItem(*b_it);
+          for (typename BlueFlow::iterator b_in = incoming[*b_it].begin();
+               b_in != incoming[*b_it].end(); ++b_in) {
+            RedNode r = _bpg.asRedNodeUnsafe(_bpg.source(*b_in));
+            for (OutArcIt r_oi(_bpg, r); r_oi!=INVALID; ++r_oi) {
+              BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(r_oi));
+              if (!reached[b]) {
+                next_blue_nodes.push_back(b);
+                reached.set(b, true);
+              }
+            }
+          }
+        }
+        _elevator->initNewLevel();
+        next_blue_nodes.swap(act_blue_nodes);
+        next_blue_nodes.clear();
+      }
+      _elevator->initFinish();
+
+      // Activate blue nodes with 0 excess, i. e. those
+      // which does not get any flow.
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (_excess[b] == 0  && (*_elevator)[b] < _elevator->maxLevel()) {
+          _elevator->activate(b);
+        }
+      }
+    }
+
+    PreflowBpMatching2() {}
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructs the class on the given bipartite graph.
+    PreflowBpMatching2(const BpGraph& bpg) : _bpg(bpg),
+                                             _matching(0),
+                                             _local_matching(false),
+                                             _elevator(0),
+                                             _local_elevator(false),
+                                             _flow(bpg, INVALID),
+                                             _excess(_bpg, 0)
+
+    {}
+
+    /// \brief Destructor
+    ///
+    /// Destructor.
+    ~PreflowBpMatching2() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the matching map
+    ///
+    /// Sets the matching map. It should contain a valid matching,
+    /// for example the empty matching is fine.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
+    ///
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching2& matchingMap(MatchingMap& map) {
+      _matching = &map;
+      _local_matching = false;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the matching map.
+    ///
+    /// Returns a const reference to the matching map, which contains
+    /// the matching edge for each node (or INVALID for unmatched nodes).
+    const MatchingMap& matchingMap() const {
+      return *_matching;
+    }
+
+    /// \brief Sets the elevator used by algorithm.
+    ///
+    /// Sets the elevator.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated elevator,
+    /// of course.
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching2& elevator(Elevator_t& elevator) {
+      _elevator = &elevator;
+      _local_elevator = false;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the elevator.
+    ///
+    /// Returns a const reference to the elevator, which keeps track
+    /// of the distance labels of nodes.
+    const Elevator_t& elevator() const {
+      return *_elevator;
+    }
+
+    /// \brief Initializes the algorithm
+    ///
+    /// Allocates the elevator and the matching map if necessary,
+    /// and initializes the elevator.
+    ///
+    /// \pre \ref init() is not called.
+    void init() {
+      createStructures();
+      initElevator();
+    }
+
+    /// \brief Smarter initialization of the matching.
+    ///
+    /// Allocates the matching map if necessary, and initializes
+    /// the algorithm with a trivially found matching: iterating
+    /// on every edge, take it in if both endnodes are unmatched.
+    ///
+    /// \return Size of the so found matching.
+    ///
+    /// \note heuristicInit() replaces init() regarding the preconditions.
+    int heuristicInit() {
+      createStructures();
+      int size = 0;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        bool matched = false;
+        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
+          if ((*_matching)[_bpg.target(oi)] == INVALID) {
+            _matching->set(_bpg.source(oi), oi);
+            _matching->set(_bpg.target(oi), oi);
+            matched = true;
+            ++size;
+          }
+        }
+      }
+
+      initElevator();
+      return size;
+    }
+
+    /// \brief Initialize the matching from a map.
+    ///
+    /// Allocates the matching map if necessary, and
+    /// initializes the algorithm with a matching given as a bool edgemap,
+    /// in which an edge is true if it is in the matching.
+    /// Also initializes the elevator.
+    ///
+    /// \return \c true if the matching is valid.
+    ///
+    /// \note matchingInit() replaces init() regarding the preconditions.
+    bool matchingInit(const BoolEdgeMap& matching) {
+      createStructures();
+
+      bool valid = true;
+      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
+        if (matching[it]) {
+          Node red = _bpg.redNode(it);
+          Node blue = _bpg.blueNode(it);
+          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
+            valid = false;
+          } else {
+            _matching->set(red, it);
+            _matching->set(blue, it);
+          }
+        }
+      }
+
+      initElevator();
+      return valid;
+    }
+
+
+    /// \brief Perform a pull/relabel operation on a node
+    ///
+    /// Perform a pull/relabel operation on a node:
+    /// search an incident blue node on one level lower,
+    /// and flip its flow to the 'chosen', or relabel
+    /// if none found.
+    ///
+    /// \pre \ref init() is called, and the node is active.
+    void pull_relabel(const BlueNode& chosen) {
+
+      int min_inc_level = _elevator->maxLevel() + 1;
+      for (OutArcIt b_oi(_bpg, chosen); b_oi!=INVALID; ++b_oi) {
+        RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(b_oi));
+        BlueNode lower = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
+        if ((*_elevator)[(lower)] < (*_elevator)[chosen]) {
+          if (--_excess[lower] == 0) {
+            _elevator->activate(lower);
+          }
+          _flow[r] = _bpg.direct(b_oi, true);
+          ++_excess[chosen];
+          _elevator->deactivate(chosen);
+          break;
+        } else {
+          if ((*_elevator)[lower] < min_inc_level &&
+              lower != chosen) {
+            min_inc_level = (*_elevator)[lower];
+          }
+        }
+      }
+      // If none found, then relabel.
+      if (_excess[chosen] == 0) {
+        int chosen_level = (*_elevator)[chosen];
+        if (min_inc_level >= _elevator->maxLevel() - 1) {
+          _elevator->lift(chosen, _elevator->maxLevel());
+          _elevator->deactivate(chosen);
+        } else {
+          _elevator->lift(chosen, min_inc_level+1);
+        }
+        if (_elevator->emptyLevel(chosen_level)) {
+          _elevator->liftToTop(chosen_level);
+        }
+      }
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// Executes pull/relabel operation on the highest active node
+    /// until there is. After start(), writeMatching() should be called to
+    /// remove unnecessary flow and choose a valid maximum matching.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (_elevator->highestActive() != INVALID)
+        pull_relabel(_elevator->highestActive());
+    }
+
+    /// \brief Deduces and writes a matching from the preflow
+    ///
+    /// The underlying preflow algorithm does not use the matching map,
+    /// so before using the matching query functions ( \ref mate(), \ref
+    /// matchingSize(), etc.), you should explicitly call this function
+    /// to write a matching based on the current preflow to the matching
+    /// map. When the preflow results an ambiguous matching (a node has
+    /// incoming flow on more than one arc), it chooses arbitrary one.
+    void writeMatching() {
+      for (NodeIt it(_bpg); it != INVALID; ++it) {
+        _matching->set(it, INVALID);
+      }
+      for (RedNodeIt r(_bpg); r != INVALID; ++r) {
+        if (_flow[r] != INVALID) {
+          BlueNode mate = _bpg.blueNode(_flow[r]);
+          if (_excess[mate] > 1) {
+            --_excess[mate];
+          } else {
+            _matching->set(mate, _flow[r]);
+            _matching->set(r, _flow[r]);
+          }
+        }
+      }
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// bp1.run() is just a shorthand for:
+    ///
+    ///\code
+    /// bp1.heuristicInit();
+    /// bp1.start();
+    /// bp1.writeMatching();
+    ///\endcode
+    void run() {
+      heuristicInit();
+      start();
+      writeMatching();
+    }
+
+    /// \brief Size of the matching
+    ///
+    /// Returns the size of the current matching.
+    /// Function runs in \f$O(n)\f$ time.
+    int matchingSize() const {
+      int size = 0;
+      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
+        if ((*_matching)[it] != INVALID) ++size;
+      }
+      return size;
+    }
+
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the current
+    /// matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpg.redNode(edge)];
+    }
+
+    /// \brief Return the matching edge incident to the given node.
+    ///
+    /// This function returns the matching edge incident to the
+    /// given node in the current matching or \c INVALID if the node is
+    /// not covered by the matching.
+    Edge matching(const Node& n) const {
+      return (*_matching)[n];
+    }
+
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the current
+    /// matching or \c INVALID if the node is not covered by the matching.
+    Node mate(const Node& n) const {
+      return (*_matching)[n] != INVALID ?
+        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
+    }
+
+    /// \brief Query a vertex cover.
+    ///
+    /// This function finds a vertex cover based on the current matching,
+    /// and in \c cover sets node in the cover \c true and \c false
+    /// otherwise. In case of maximal matching, this will be a minimal
+    /// vertex cover.
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return The number of nodes in the cover.
+    int cover(BoolNodeMap &cover) const {
+      int count = 0;
+      std::queue<RedNode> q;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
+        cover[b] = false;
+
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          cover[r] = false;
+          q.push(r);
+        } else {
+          cover[r] = true;
+          ++count;
+        }
+      }
+
+      while (!q.empty()) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (cover[blue]) continue;
+
+          cover[blue] = true;
+          ++count;
+
+          if ((*_matching)[blue] != INVALID) {
+            RedNode nr = _bpg.redNode((*_matching)[blue]);
+            q.push(nr);
+            cover[nr] = false;
+            --count;
+          }
+        }
+        q.pop();
+      }
+
+      return count;
+    }
+
+    /// \brief Query a barrier of red nodes.
+    ///
+    /// This function tries to create a barrier of red nodes, which should:
+    ///  - contain every unmatched red nodes,
+    ///  - contain exactly that many matched red nodes as much blue nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool redBarrier(BoolRedNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<RedNode> q;
+
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          barrier[r] = true;
+          q.push(r);
+        } else {
+          barrier[r] = false;
+        }
+      }
+
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+
+          if (exist &= ((*_matching)[blue] != INVALID)) {
+            RedNode r = _bpg.redNode((*_matching)[blue]);
+            if (!barrier[r]) {
+              q.push(r);
+              barrier[r] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+
+      return exist;
+    }
+
+    /// \brief Query a barrier of blue nodes.
+    ///
+    /// This function tries to create a barrier of blue nodes, which should:
+    ///  - contain every unmatched blue nodes,
+    ///  - contain exactly that many matched blue nodes as much red nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool blueBarrier(BoolBlueNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<BlueNode> q;
+
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        if ((*_matching)[b] == INVALID) {
+          barrier[b] = true;
+          q.push(b);
+        } else {
+          barrier[b] = false;
+        }
+      }
+
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
+
+          if (exist &= ((*_matching)[r] != INVALID)) {
+            BlueNode b = _bpg.blueNode((*_matching)[r]);
+            if (!barrier[b]) {
+              q.push(b);
+              barrier[b] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+
+      return exist;
+    }
+  };
+
+}
+#endif

test/CMakeLists.txt

@@ -48 +48 @@
   path_test
   planarity_test
   preflow_test
+  preflow_bp_matching_2_test
   preflow_bp_matching_test
   radix_sort_test
   random_test

test/Makefile.am

@@ -46 +46 @@
         test/path_test \
         test/planarity_test \
         test/preflow_test \
+        test/preflow_bp_matching_2_test \
         test/preflow_bp_matching_test \
         test/radix_sort_test \
         test/random_test \
@@ -105 +106 @@
 test_path_test_SOURCES = test/path_test.cc
 test_planarity_test_SOURCES = test/planarity_test.cc
 test_preflow_test_SOURCES = test/preflow_test.cc
+test_preflow_bp_matching_2_test_SOURCES = test/preflow_bp_matching_2_test.cc
 test_preflow_bp_matching_test_SOURCES = test/preflow_bp_matching_test.cc
 test_radix_sort_test_SOURCES = test/radix_sort_test.cc
 test_suurballe_test_SOURCES = test/suurballe_test.cc

new file test/preflow_bp_matching_2_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/preflow_bp_matching_2.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 = 3;
+const string lgf[lgfn] = {
+  // Empty graph
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+
+  // No red nodes
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "1\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+
+  // No blue nodes
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+};
+
+const string u_lgf =
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "@blue_nodes\n"
+  "label\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "label\n"
+  "1 4 1\n"
+  "1 5 2\n"
+  "2 5 3\n"
+  "2 6 4\n"
+  "2 7 5\n"
+  "3 7 6\n"
+;
+
+
+
+void checkPreflowBpMatching2Compile() {
+  typedef concepts::BpGraph BpGraph;
+  BPGRAPH_TYPEDEFS(BpGraph);
+  typedef PreflowBpMatching2<BpGraph> PBM;
+
+  BpGraph graph;
+  RedNode red;
+  BlueNode blue;
+  Node node;
+  Edge edge;
+  BoolEdgeMap init_matching(graph);
+  PBM pbm_test(graph);
+  const PBM& const_pbm_test = pbm_test;
+  PBM::MatchingMap matching_map(graph, INVALID);
+  PBM::Elevator_t elevator(graph);
+  BpGraph::NodeMap<bool> cov(graph);
+  int cover_size;
+  BoolRedNodeMap rb(graph);
+  BoolBlueNodeMap bb(graph);
+
+  pbm_test.matchingMap(matching_map);
+  pbm_test.elevator(elevator);
+  pbm_test.init();
+  pbm_test.heuristicInit();
+  pbm_test.matchingInit(init_matching);
+
+  pbm_test.start();
+  pbm_test.pull_relabel(blue);
+  pbm_test.run();
+  pbm_test.writeMatching();
+
+  const_pbm_test.matchingSize();
+  const_pbm_test.matching(node);
+  const_pbm_test.matching(edge);
+  const_pbm_test.mate(red);
+  const_pbm_test.mate(blue);
+  const_pbm_test.mate(node);
+  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
+  edge = max_matching[node];
+  const PBM::Elevator_t& elev = const_pbm_test.elevator();
+  blue = elev.highestActive();
+  cover_size = const_pbm_test.cover(cov);
+  const_pbm_test.redBarrier(rb);
+  const_pbm_test.blueBarrier(bb);
+}
+
+typedef SmartBpGraph BpGraph;
+typedef PreflowBpMatching2<BpGraph> PBM;
+BPGRAPH_TYPEDEFS(BpGraph);
+
+void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
+  BoolBlueNodeMap matched(bpg, false);
+  for (RedNodeIt it(bpg); it != INVALID; ++it) {
+    if (pbm.matching(it) != INVALID) {
+      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
+      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
+    }
+  }
+  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
+    // != is used as logical xor here
+    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
+  }
+}
+
+void checkCover(const BpGraph& bpg, const PBM& pbm) {
+  int cov_size;
+  BoolNodeMap cov(bpg);
+  cov_size = pbm.cover(cov);
+  BoolEdgeMap covered(bpg, false);
+  for (NodeIt n(bpg); n!=INVALID; ++n) {
+    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
+      covered[e] = true;
+    }
+  }
+  bool all = true;
+  for (EdgeIt e(bpg); e!=INVALID; ++e) {
+    all&=covered[e];
+  }
+
+  check(all, "Invalid cover");
+  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
+}
+void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
+  BoolRedNodeMap rb(bpg);
+  BoolBlueNodeMap bb(bpg);
+
+  check(pbm.redBarrier(rb), "Suboptimal matching.");
+  check(pbm.blueBarrier(bb), "Suboptimal matching.");
+
+  BoolNodeMap reached(bpg, false);
+
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (!rb[r]) continue;
+    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (!bb[b]) continue;
+    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+
+  int ur = 0, ub = 0, rc = 0, bc = 0;
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (rb[r]) ++rc;
+    if (pbm.mate(r) == INVALID) ++ur;
+    if (reached[r]) --bc;
+  }
+
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (bb[b]) ++bc;
+    if (pbm.mate(b) == INVALID) ++ub;
+    if (reached[b]) --rc;
+  }
+
+  check(rc == ur, "Red barrier is wrong.");
+  check(bc == ub, "Blue barrier is wrong.");
+}
+
+void checkMatching(const BpGraph& bpg, const PBM& pbm) {
+  checkMatchingValidity(bpg, pbm);
+  checkBarriers(bpg, pbm);
+  checkCover(bpg, pbm);
+}
+
+
+int main() {
+  // Checking hard wired extremal graphs
+  for (int i=0; i<lgfn; ++i) {
+    BpGraph bpg;
+    istringstream lgfs(lgf[i]);
+    bpGraphReader(bpg, lgfs).run();
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+
+  // Checking some use cases
+  BpGraph bpg;
+  istringstream u_lgfs(u_lgf);
+  bpGraphReader(bpg, u_lgfs).run();
+
+  {
+    PBM pbm(bpg);
+    pbm.init();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.pull_relabel(b);
+    }
+    pbm.start();
+    pbm.writeMatching();
+    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
+    checkMatching(bpg, pbm);
+  }
+
+  {
+    PBM pbm(bpg);
+    pbm.heuristicInit();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.pull_relabel(b);
+    }
+    pbm.start();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    pbm.run();
+    checkMatching(bpg, pbm);
+
+  }
+  {
+    PBM pbm(bpg);
+    BpGraph::EdgeMap<bool> init_matching(bpg, false);
+    BpGraph::Edge e;
+    bpg.first(e);
+    init_matching[e] = true;
+    check(pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    BpGraph::EdgeMap<bool> init_matching(bpg, true);
+    check(!pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+
+  // Checking some random generated graphs
+  const int random_test_n = 10;
+  const int max_red_n = 1000;
+  const int min_red_n = 10;
+  const int max_blue_n = 1000;
+  const int min_blue_n = 10;
+  for (int i=0; i<random_test_n; ++i) {
+    int red_n = rnd.integer(min_red_n, max_red_n);
+    int blue_n = rnd.integer(min_blue_n, max_blue_n);
+    double prob = rnd();
+
+    SmartBpGraph bpg;
+    for (int j=0; j<red_n; ++j) {
+      bpg.addRedNode();
+    }
+    for (int j=0; j<blue_n; ++j) {
+      bpg.addBlueNode();
+    }
+    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (rnd() < prob/10.) {
+          bpg.addEdge(r,b);
+        }
+      }
+    }
+
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+
+    // do not always use the highest active node
+    PBM pbm2(bpg);
+    PBM::Elevator_t elev(bpg);
+    pbm2.elevator(elev);
+    pbm2.init();
+    while (elev.highestActive() != INVALID) {
+      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (elev.active(b)) pbm2.pull_relabel(b);
+      }
+    }
+    pbm2.writeMatching();
+    checkMatching(bpg, pbm2);
+  }
+
+
+  return 0;
+}
+
+

lemon/dimacs.h

@@ -2 +2 @@
  *
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2003-2010
+ * Copyright (C) 2003-2012
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
  *
@@ -45 +45 @@
       MIN,   ///< DIMACS file type for minimum cost flow problems.
       MAX,   ///< DIMACS file type for maximum flow problems.
       SP,    ///< DIMACS file type for shostest path problems.
-      MAT    ///< DIMACS file type for plain graphs and matching problems.
+      MAT,   ///< DIMACS file type for plain graphs and matching problems.
+      ASN,   ///< DIMACS file type for assignment problems.
+      UBM,   ///< DIMACS file type for unweighted bipartite matching problems.
+      WBM    ///< DIMACS file type for weighted bipartite matching problems.
     };
     ///The file type
     Type type;
@@ -82 +85 @@
               else if(problem=="max") r.type=DimacsDescriptor::MAX;
               else if(problem=="sp") r.type=DimacsDescriptor::SP;
               else if(problem=="mat") r.type=DimacsDescriptor::MAT;
+              else if(problem=="asn") r.type=DimacsDescriptor::ASN;
+              else if(problem=="ubm") r.type=DimacsDescriptor::UBM;
+              else if(problem=="wbm") r.type=DimacsDescriptor::WBM;
               else throw FormatError("Unknown problem type");
               return r;
             }
@@ -101 +107 @@
   }
 
 
+
+  /// \brief DIMACS unweighted bipartite matching problem reader function.
+  ///
+  /// This function reads an unweighted bipartite matching problem instance
+  /// from DIMACS format, i.e. from a DIMACS file having a line starting with
+  /// \code
+  ///   p ubm
+  /// \endcode
+  /// and creates the corresponding bipartite graph.
+  ///
+  /// At the beginning, \c g is cleared by \c g.clear().
+  ///
+  /// If the file type was previously evaluated by dimacsType(), then
+  /// the descriptor struct should be given by the \c dest parameter.
+  template <typename BpGraph>
+  void readDimacsUbm(std::istream& is,
+                     BpGraph &g,
+                     DimacsDescriptor desc=DimacsDescriptor())
+  {
+    g.clear();
+
+    if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+    if(desc.type!=DimacsDescriptor::UBM)
+      throw FormatError("Problem type mismatch");
+
+    std::vector<typename BpGraph::Node> nodes(desc.nodeNum + 1, INVALID);
+    typename BpGraph::Edge e;
+    std::string problem, str;
+    char c;
+    int i, j;
+
+    // eat up comment lines and problem line
+    while (is >> c && (c=='c' || c=='p')) {
+      getline(is, str);
+    }
+    is.unget();
+
+    // reading red nodes
+    while (is >> c && c == 'n') {
+      is >> i;
+      getline(is, str);
+      nodes[i] = g.addRedNode();
+    }
+    is.unget();
+
+    // remaining nodes are blue:
+    for (int k = 1; k <= desc.nodeNum; ++k) {
+      if (nodes[k] == INVALID) nodes[k] = g.addBlueNode();
+    }
+
+    // reading edges
+    while (is >> c && c == 'a') {
+      is >> i >> j;
+      getline(is, str);
+      e = g.addEdge(g.asRedNode(nodes[i]), g.asBlueNode(nodes[j]));
+    }
+  }
+
+
+  /// \brief DIMACS weighted bipartite matching problem reader function.
+  ///
+  /// This function reads a weighted bipartite matching problem instance from
+  /// DIMACS format, i.e. from a DIMACS file having a line starting with
+  /// \code
+  ///   p wbm
+  /// \endcode
+  /// and creates the corresponding bipartite graph and weight map.
+  ///
+  /// At the beginning, \c g is cleared by \c g.clear().
+  ///
+  /// If the file type was previously evaluated by dimacsType(), then
+  /// the descriptor struct should be given by the \c dest parameter.
+  template <typename BpGraph, typename WeightMap>
+  void readDimacsWbm(std::istream& is,
+                     BpGraph &g,
+                     WeightMap &m,
+                     DimacsDescriptor desc=DimacsDescriptor())
+  {
+    g.clear();
+
+    if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+    if(desc.type!=DimacsDescriptor::WBM)
+      throw FormatError("Problem type mismatch");
+
+    std::vector<typename BpGraph::Node> nodes(desc.nodeNum + 1, INVALID);
+    typename BpGraph::Edge e;
+    std::string problem, str;
+    char c;
+    int i, j;
+    typename WeightMap::Value w;
+
+    // eat up comment lines and problem line
+    while (is >> c && (c=='c' || c=='p')) {
+      getline(is, str);
+    }
+    is.unget();
+
+    // reading red nodes
+    while (is >> c && c == 'n') {
+      is >> i;
+      getline(is, str);
+      nodes[i] = g.addRedNode();
+    }
+    is.unget();
+
+    // remaining nodes are blue:
+    for (int k = 1; k <= desc.nodeNum; ++k) {
+      if (nodes[k] == INVALID) nodes[k] = g.addBlueNode();
+    }
+
+    // reading edges and weights
+    while (is >> c && c == 'a') {
+      is >> i >> j >> w;
+      getline(is, str);
+      e = g.addEdge(g.asRedNode(nodes[i]), g.asBlueNode(nodes[j]));
+      m.set(e, w);
+    }
+  }
+
+
+  /// \brief DIMACS assignment problem reader function.
+  ///
+  /// This function reads an assignment problem instance from DIMACS format,
+  /// i.e. from a DIMACS file having a line starting with
+  /// \code
+  ///   p asn
+  /// \endcode
+  /// and creates the corresponding digraph, arc cost- and supply map.
+  ///
+  /// At the beginning, \c g is cleared by \c g.clear().
+  /// Cost of the arcs are written to the cost map,
+  /// and the supply map contains the supply of the nodes, which is
+  /// 1 for sources and -1 for sink nodes.
+  ///
+  /// If the file type was previously evaluated by dimacsType(), then
+  /// the descriptor struct should be given by the \c dest parameter.
+  template <typename Digraph, typename CostMap>
+  void readDimacsAsn(std::istream& is,
+                     Digraph &g,
+                     CostMap &cost,
+                     typename Digraph::template NodeMap<int> &supply,
+                     DimacsDescriptor desc=DimacsDescriptor())
+  {
+    g.clear();
+
+    if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+    if(desc.type!=DimacsDescriptor::ASN)
+      throw FormatError("Problem type mismatch");
+
+    std::vector<typename Digraph::Node> nodes(desc.nodeNum + 1, INVALID);
+    typename Digraph::Arc e;
+    std::string problem, str;
+    char c;
+    int i, j;
+    typename CostMap::Value co;
+
+    for (int k = 1; k <= desc.nodeNum; ++k) {
+      nodes[k] = g.addNode();
+      supply.set(nodes[k], -1);
+    }
+
+    // eat up comment lines and problem line
+    while (is >> c && (c=='c' || c=='p')) {
+      getline(is, str);
+    }
+    is.unget();
+
+    // reading source nodes
+    while (is >> c && c == 'n') {
+      is >> i;
+      getline(is, str);
+      supply.set(nodes[i], 1);
+    }
+    is.unget();
+
+    // reading arcs
+    while (is >> c && c == 'a') {
+      is >> i >> j >> co;
+      getline(is, str);
+      e = g.addArc(nodes[i], nodes[j]);
+      cost.set(e, co);
+    }
+  }
+
   /// \brief DIMACS minimum cost flow reader function.
   ///
   /// This function reads a minimum cost flow instance from DIMACS format,

tools/CMakeLists.txt

@@ -16 +16 @@
 ADD_EXECUTABLE(dimacs-solver dimacs-solver.cc)
 TARGET_LINK_LIBRARIES(dimacs-solver lemon)
 
+ADD_EXECUTABLE(weighted-bp-gen weighted-bp-gen.cc)
+TARGET_LINK_LIBRARIES(weighted-bp-gen lemon)
+
+ADD_EXECUTABLE(unweighted-bp-gen unweighted-bp-gen.cc)
+TARGET_LINK_LIBRARIES(unweighted-bp-gen lemon)
+
+ADD_EXECUTABLE(bp-matching-benchmark bp-matching-benchmark.cc)
+TARGET_LINK_LIBRARIES(bp-matching-benchmark lemon)
+
 INSTALL(
   TARGETS lgf-gen dimacs-to-lgf dimacs-solver
   RUNTIME DESTINATION bin

tools/Makefile.am

@@ -4 +4 @@
 if WANT_TOOLS
 
 bin_PROGRAMS += \
+        tools/bp-matching-benchmark \
         tools/dimacs-solver \
         tools/dimacs-to-lgf \
-        tools/lgf-gen
+        tools/lgf-gen \
+        tools/unweighted-bp-gen \
+        tools/weighted-bp-gen
 
 dist_bin_SCRIPTS += tools/lemon-0.x-to-1.x.sh
 
 endif WANT_TOOLS
 
+tools_bp_matching_benchmark_SOURCES = tools/bp-matching-benchmark.cc
 tools_dimacs_solver_SOURCES = tools/dimacs-solver.cc
 tools_dimacs_to_lgf_SOURCES = tools/dimacs-to-lgf.cc
 tools_lgf_gen_SOURCES = tools/lgf-gen.cc
+tools_unweighted_bp_gen_SOURCES = tools/unweighted-bp-gen.cc
+tools_weighted_bp_gen_SOURCES = tools/weighted-bp-gen.cc

new file tools/bp-matching-benchmark.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 <lemon/dimacs.h>
+#include <lemon/maps.h>
+#include <lemon/smart_graph.h>
+#include <lemon/time_measure.h>
+#include <lemon/hopcroft_karp.h>
+#include <lemon/matching.h>
+#include <lemon/preflow.h>
+#include <lemon/preflow_bp_matching.h>
+#include <lemon/preflow_bp_matching_2.h>
+#include <lemon/bp_matching.h>
+#include <lemon/cost_scaling.h>
+#include <lemon/capacity_scaling.h>
+#include <lemon/network_simplex.h>
+#include <iostream>
+
+using namespace std;
+using namespace lemon;
+
+void check(bool c, const string &msg) {
+  if (!c) {
+    std::cerr << "Error: " << msg << "\n";
+  }
+}
+
+template <typename BpGraph, typename Digraph>
+void create_network(const BpGraph &bpg,
+                    Digraph &network,
+                    typename Digraph::Node &source,
+                    typename Digraph::Node &sink)
+{
+  network.clear();
+  source = network.addNode();
+  sink = network.addNode();
+
+  typename BpGraph::template NodeMap<typename Digraph::Node> xref(bpg);
+  for (typename BpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+    xref.set(r, network.addNode());
+    network.addArc(source, xref[r]);
+  }
+  for (typename BpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    xref.set(b, network.addNode());
+    network.addArc(xref[b], sink);
+  }
+
+  for (typename BpGraph::EdgeIt e(bpg); e!=INVALID; ++e) {
+    network.addArc(xref[bpg.redNode(e)], xref[bpg.blueNode(e)]);
+  }
+}
+
+template <typename BpGraph, typename BpCostMap,
+          typename Digraph, typename DgCostMap, typename SupplyMap>
+void create_weighted_net(const BpGraph &bpg,
+                         const BpCostMap &bpcost,
+                         Digraph &dg,
+                         DgCostMap &dgcost,
+                         SupplyMap &supply)
+{
+  dg.clear();
+  typename Digraph::Node source = dg.addNode();
+  supply.set(source, countBlueNodes(bpg) - countRedNodes(bpg));
+
+  typename BpGraph::template NodeMap<typename Digraph::Node> xref(bpg);
+
+  for (typename BpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+    xref.set(r, dg.addNode());
+    dgcost.set(dg.addArc(xref[r], source), 0);
+    supply.set(xref[r], 1);
+  }
+  for (typename BpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    xref.set(b, dg.addNode());
+    dgcost.set(dg.addArc(source, xref[b]), 0);
+    supply.set(xref[b], -1);
+  }
+
+  for (typename BpGraph::EdgeIt e(bpg); e!=INVALID; ++e) {
+    typename Digraph::Arc a =
+      dg.addArc(xref[bpg.redNode(e)], xref[bpg.blueNode(e)]);
+    dgcost.set(a, -bpcost[e]);
+  }
+}
+
+typedef SmartDigraph Digraph;
+typedef SmartBpGraph BpGraph;
+
+///\ingroup tools
+///\file
+///\brief Benchmark program for weighted and unweighted bipartite matching.
+///
+/// The program input is a Dimacs file with a problem line starting
+/// \code
+///   p ubm
+/// \endcode
+/// or
+/// \code
+///   p wbm
+/// \endcode
+///
+/// Depending on the input, it runs the \ref lemon::HopcroftKarp
+/// "Hopcroft-Karp", the \ref lemon::MaxMatching "general matching"
+/// and three algorithms based on preflow (\ref lemon::Preflow "preflow
+/// on extended graph", \ref lemon::PreflowBpMatching,
+/// \ref lemon::PreflowBpMatching2") for finding
+/// a maximum cardinality bipartite matching, or the
+/// \ref lemon::MaxWeightedBpMatching "max. weighted matching
+/// for sparse bipartite graphs", the \ref lemon::MaxWeightedMatching
+/// "max. weighted matching for general graph", \ref
+/// lemon::CostScaling "cost scaling", \ref lemon::CapacityScaling
+/// "capacity scaling" and \ref lemon::NetworkSimplex "network
+/// simplex" algorithms to find the maximum weighted matching.
+int main() {
+  DimacsDescriptor desc = dimacsType(cin);
+  if (desc.type == DimacsDescriptor::UBM) {
+    // unweighted bipartite matching
+    Timer hk_t(false),
+      mm_t(false),
+      pf_t(false),
+      pbm_t(false),
+      pbm2_t(false);
+
+    BpGraph *bpg = new BpGraph();
+    readDimacsUbm(cin, *bpg, desc);
+
+    {
+      hk_t.start();
+      HopcroftKarp<BpGraph> hk(*bpg);
+      hk.run();
+      hk_t.stop();
+      BpGraph::NodeMap<bool> m(*bpg);
+      check(hk.matchingSize() == hk.cover(m),
+            "Suboptimal matching in Hopcroft-Karp.");
+    }
+    {
+      mm_t.start();
+      MaxMatching<BpGraph> mm(*bpg);
+      mm.run();
+      mm_t.stop();
+    }
+    {
+      pbm_t.start();
+      PreflowBpMatching<BpGraph> pbm(*bpg);
+      pbm.run();
+      pbm_t.stop();
+      BpGraph::NodeMap<bool> m(*bpg);
+      check(pbm.matchingSize() == pbm.cover(m),
+            "Suboptimal matching in push-relabel.");
+    }
+    {
+      pbm2_t.start();
+      PreflowBpMatching2<BpGraph> pbm2(*bpg);
+      pbm2.run();
+      pbm2_t.stop();
+      BpGraph::NodeMap<bool> m(*bpg);
+      check(pbm2.matchingSize() == pbm2.cover(m),
+            "Suboptimal matching in pull-relabel.");
+    }
+
+    Digraph *network = new Digraph();
+    Digraph::Node source, target;
+    create_network(*bpg, *network, source, target);
+    int rednum = countRedNodes(*bpg),
+      bluenum = countBlueNodes(*bpg),
+      edgenum = desc.edgeNum;
+
+    delete bpg;
+
+    {
+      pf_t.start();
+      Preflow<Digraph, ConstMap<Digraph::Arc, long> >
+        pf(*network, constMap<Digraph::Arc, long>(1), source, target);
+      pf.run();
+      pf_t.stop();
+    }
+    delete network;
+
+    cout << "Benchmarking in unweighted case, using a bipartite graph\n"
+         << "with " << rednum << " red, " << bluenum << " blue nodes, and "
+         << edgenum << " edges.\n"
+         << "--------------------------------------------------------\n"
+         << "Algorithm used               Time\n"
+         << "Hopcroft-Karp            " << hk_t.realTime() << "\n"
+         << "General matching         " << mm_t.realTime() << "\n"
+         << "Push-relabel matching    " << pbm_t.realTime() << "\n"
+         << "Pull-relabel matching    " << pbm2_t.realTime() << "\n"
+         << "Preflow                  " << pf_t.realTime() << endl;
+
+  } else if (desc.type == DimacsDescriptor::WBM) {
+    // weighted bipartite matching
+    Timer bm_t(false),
+      gm_t(false),
+      cas_t(false),
+      cos_t(false),
+      ns_t(false);
+
+    BpGraph *bpg = new BpGraph();
+    BpGraph::EdgeMap<long> *weight = new BpGraph::EdgeMap<long>(*bpg);
+    readDimacsWbm(cin, *bpg, *weight, desc);
+
+    {
+      bm_t.start();
+      MaxWeightedBpMatching<BpGraph, BpGraph::EdgeMap<long> > bm(*bpg, *weight);
+      bm.run();
+      bm_t.stop();
+    }
+    {
+      gm_t.start();
+      MaxWeightedMatching<BpGraph, BpGraph::EdgeMap<long> > gm(*bpg, *weight);
+      gm.run();
+      gm_t.stop();
+    }
+
+    Digraph *g = new Digraph();
+    Digraph::ArcMap<long> *cost = new Digraph::ArcMap<long>(*g);
+    Digraph::NodeMap<long> *supply = new Digraph::NodeMap<long>(*g);
+    create_weighted_net(*bpg, *weight, *g, *cost, *supply);
+
+    int rednum = countRedNodes(*bpg),
+      bluenum = countBlueNodes(*bpg),
+      edgenum = desc.edgeNum;
+
+    delete bpg;
+    delete weight;
+    {
+      cas_t.start();
+      CapacityScaling<Digraph> cas(*g);
+      cas.costMap(*cost)
+        .supplyMap(*supply)
+        .upperMap(constMap<Digraph::Arc, long>(1));
+      cas.run();
+      cas_t.stop();
+    }
+    {
+      cos_t.start();
+      CostScaling<Digraph> cos(*g);
+      cos.costMap(*cost)
+        .supplyMap(*supply)
+        .upperMap(constMap<Digraph::Arc, long>(1));
+      cos.run();
+      cos_t.stop();
+    }
+    {
+      ns_t.start();
+      NetworkSimplex<Digraph> ns(*g);
+      ns.costMap(*cost)
+        .supplyMap(*supply)
+        .upperMap(constMap<Digraph::Arc, long>(1));
+      ns.run();
+      ns_t.stop();
+    }
+    delete g;
+    delete cost;
+    delete supply;
+
+    cout << "Benchmarking on weighted bipartite matching problem on a graph\n"
+         << "with " << rednum << " red, " << bluenum << " blue nodes and "
+         << edgenum << " edges\n"
+         << "--------------------------------------------------------\n"
+         << "Algorithm used             Time\n"
+         << "Bipartite matching       " << bm_t.realTime() << "\n"
+         << "General matching         " << gm_t.realTime() << "\n"
+         << "Capacity scaling         " << cas_t.realTime() << "\n"
+         << "Cost scaling             " << cos_t.realTime() << "\n"
+         << "Network simplex          " << ns_t.realTime() << endl;
+
+  } else {
+    cerr << "Wrong problem type." << endl;
+    return 1;
+  }
+
+  return 0;
+}

new file tools/unweighted-bp-gen.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 <lemon/arg_parser.h>
+#include <lemon/random.h>
+
+using namespace lemon;
+using namespace std;
+
+///\ingroup tools
+///\file
+///\brief Random bipartite graph generator.
+///
+/// This program generates an unweighted bipartite matching problem.
+/// The output is in Dimacs format, with problem type of 'ubm'. The
+/// enumerated nodes form one set of the bipartite graph, and
+/// the edges do not have any label.
+int main(int argc, char** argv) {
+  int seed = 1,
+    red = 5,
+    blue = 5,
+    density = 10;
+
+  ArgParser arg_p(argc, argv);
+
+  arg_p.refOption("seed",
+                  "Seed of the random number generator (default: 1)",
+                  seed)
+    .refOption("red",
+               "Number of red nodes (default: 5)",
+               red)
+    .refOption("blue",
+               "Number of blue nodes (default: 5)",
+               blue)
+    .refOption("density",
+               "Number of edges (default: 10)",
+               density);
+
+  arg_p.synonym("s", "seed")
+    .synonym("r", "red")
+    .synonym("b", "blue")
+    .synonym("d", "density");
+
+  arg_p.parse();
+
+  rnd.seed(seed);
+
+
+  long max_edge = static_cast<long>(red)*blue;
+  long edge_left = density;
+
+  if (max_edge < density) {
+    cerr << "Too many edges!\n"
+         << "At most " << max_edge << " edges can be placed in the graph.\n";
+    return 1;
+  }
+
+  cout << "c Random bipartite graph for bipartite matching\n"
+       << "c----------------------------------------------\n"
+       << "c Input parameters:\n"
+       << "c   Random seed: " << seed << "\n"
+       << "c   Number of nodes: " << red + blue << "\n"
+       << "c     Number of red nodes:  " << red << "\n"
+       << "c     Number of blue nodes: " << blue << "\n"
+       << "c   Number of edges: " << density << "\n"
+       << "c----------------------------------------------\n";
+  cout << "p ubm " << red + blue << " " << density << "\n";
+
+  // red nodes:
+  for (int r=1; r <= red; ++r) {
+    cout << "n " << r << "\n";
+  }
+
+  // edges:
+  for (int r = 1; r <= red; ++r) {
+    for (int b = red+1; b <= red + blue; ++b) {
+      if (rnd.boolean(static_cast<long double>(edge_left) / max_edge)) {
+        --edge_left;
+        cout << "a " << r << " " << b << "\n";
+      }
+      --max_edge;
+    }
+  }
+
+  return 0;
+}
+

new file tools/weighted-bp-gen.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 <lemon/arg_parser.h>
+#include <lemon/random.h>
+
+using namespace lemon;
+using namespace std;
+
+///\ingroup tools
+///\file
+///\brief Random weighted bipartite graph generator.
+///
+/// This program generates a weighted bipartite matching problem.
+/// The output is in Dimacs format, with problem type of 'wbm'. The
+/// enumerated nodes form one set of the bipartite graph, and
+/// the edge labels denote weights.
+int main(int argc, char** argv) {
+  int seed = 1,
+    red = 5,
+    blue = 5,
+    density = 10,
+    minweight = 10,
+    maxweight = 100;
+
+  ArgParser arg_p(argc, argv);
+
+  arg_p.refOption("seed",
+                  "Seed of the random number generator (default: 1)",
+                  seed)
+    .refOption("red",
+               "Number of red nodes (default: 5)",
+               red)
+    .refOption("blue",
+               "Number of blue nodes (default: 5)",
+               blue)
+    .refOption("density",
+               "Number of edges (default: 10)",
+               density)
+    .refOption("minweight",
+               "Minimal weight of edges (default: 10)",
+               minweight)
+    .refOption("maxweight",
+               "Maximal weight of edges (default: 100",
+               maxweight);
+
+  arg_p.synonym("s", "seed")
+    .synonym("r", "red")
+    .synonym("b", "blue")
+    .synonym("d", "density")
+    .synonym("min", "minweight")
+    .synonym("max", "maxweight");
+
+  arg_p.parse();
+
+  rnd.seed(seed);
+
+
+  long max_edge = static_cast<long>(red)*blue;
+  long edge_left = density;
+
+  if (max_edge < density) {
+    cerr << "Too many edges!\n"
+         << "At most " << max_edge << " edges can be placed in the graph.\n";
+    return 1;
+  }
+
+  cout << "c Random bipartite graph for weighted bipartite matching\n"
+       << "c----------------------------------------------\n"
+       << "c Input parameters:\n"
+       << "c   Random seed: " << seed << "\n"
+       << "c   Number of nodes: " << red + blue << "\n"
+       << "c     Number of red nodes:  " << red << "\n"
+       << "c     Number of blue nodes: " << blue << "\n"
+       << "c   Number of edges: " << density << "\n"
+       << "c   Minimal weight of edges: " << minweight << "\n"
+       << "c   Maximal weight of edges: " << maxweight << "\n"
+       << "c----------------------------------------------\n";
+  cout << "p wbm " << red + blue << " " << density << "\n";
+
+  // red nodes:
+  for (int r=1; r <= red; ++r) {
+    cout << "n " << r << "\n";
+  }
+
+  // edges:
+  for (int r = 1; r <= red; ++r) {
+    for (int b = red+1; b <= red + blue; ++b) {
+      if (rnd.boolean(static_cast<long double>(edge_left) / max_edge)) {
+        --edge_left;
+        cout << "a " << r << " " << b << " "
+             << rnd.integer(minweight, maxweight+1) << "\n";
+      }
+      --max_edge;
+    }
+  }
+
+  return 0;
+}
+