COIN-OR::LEMON - Graph Library

Ticket #168: bpm_pack.patch

File bpm_pack.patch, 124.6 KB (added by Poroszkai Daniel, 9 years ago)
  • lemon/bp_matching.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1328397120 -3600
    # Node ID 2cf5b6b4bdcbba2495741f4aa7d6f54b02640dad
    # Parent  7b95667b0241007b7f725d0bf286f7bbdb82920e
    Update weighted bp matching to work with typesafe bipartite node sets
    
    diff --git a/lemon/bp_matching.h b/lemon/bp_matching.h
    a b  
    4545    /// The graph type of the algorithm
    4646    typedef BGR BpGraph;
    4747    /// The type of the matching map on the red nodes
    48     typedef typename BpGraph::template RedMap<typename BpGraph::Edge>
     48    typedef typename BpGraph::template RedNodeMap<typename BpGraph::Edge>
    4949    RedMatchingMap;
    5050    /// The type of the matching map on the blue nodes
    51     typedef typename BpGraph::template BlueMap<typename BpGraph::Edge>
     51    typedef typename BpGraph::template BlueNodeMap<typename BpGraph::Edge>
    5252    BlueMatchingMap;
    5353
    5454  private:
     
    118118      RM = new RedMatchingMap(*G, INVALID);
    119119      BM = new BlueMatchingMap(*G, INVALID);
    120120     
    121       for( RedIt n(*G); n != INVALID; ++n ){
     121      for( RedNodeIt n(*G); n != INVALID; ++n ){
    122122        for( IncEdgeIt eit(*G, n); eit != INVALID; ++eit ) {
    123123          if( (*BM)[ G->blueNode(eit) ] == INVALID ){
    124124            (*BM)[G->blueNode(eit) ] = eit;
     
    164164    /// called before using this function.
    165165    bool step()
    166166    {
    167       queue<Node> sources;
     167      queue<RedNode> sources;
    168168      BoolNodeMap _processed(*G, false);
    169169      EdgeNodeMap _pred(*G, INVALID);
    170170      BoolNodeMap _reached(*G, false);
    171       vector<Node> H;
     171      vector<RedNode> H;
    172172      bool isAlternatePath = false;
    173       for( RedIt n(*G); n != INVALID; ++n ){
     173      for( RedNodeIt n(*G); n != INVALID; ++n ){
    174174        if( (*RM)[n] == INVALID ){
    175175          sources.push( n );
    176176        }
    177177      }
    178178       
    179179      while ( !sources.empty() && !isAlternatePath ){
    180         Node s = sources.front();
     180        RedNode s = sources.front();
    181181        if(_processed[s]) {sources.pop(); continue;}
    182182        _processed[s] = true;
    183183        _reached[s] = true;
    184184        H.push_back( s );
    185         Node m;
     185        BlueNode m;
    186186        for(IncEdgeIt e(*G,s);e!=INVALID;++e) {
    187187          if(!_reached[m=G->blueNode(e)]) {
    188188            _reached[m]=true;
    189189            _pred[m]=e;
    190190            if( (*BM)[m] != INVALID ) {
    191               Node n;
     191              RedNode n;
    192192              if (!_reached[n = G->redNode((*BM)[m])] ) {
    193193                sources.push( n );
    194194                _pred[n] = (*BM)[m];
     
    257257    int matchingSize() const
    258258    {
    259259      int size = 0;
    260       for ( BlueIt n(*G); n != INVALID; ++n) {
     260      for ( BlueNodeIt n(*G); n != INVALID; ++n) {
    261261        if ( (*BM)[n] != INVALID) {
    262262          ++size;
    263263        }
     
    281281    /// not covered by the matching.
    282282    Edge matching(const Node& n) const
    283283    {
    284       return G->blue(n) ? (*BM)[n] : (*RM)[n];
     284      return G->blue(n) ? (*BM)[G->asBlueNodeUnsafe(n)] : (*RM)[G->asRedNodeUnsafe(n)];
    285285    }
    286286
    287287    /// \brief Return a const reference to the blue matching map.
     
    310310    {
    311311      Edge e;
    312312      if( G->blue(node) ){
    313         if( ( e = (*BM)[node]) != INVALID ) return G->redNode(e);
     313        if( ( e = (*BM)[G->asBlueNodeUnsafe(node)]) != INVALID ) return G->redNode(e);
    314314        else return INVALID;
    315315      }
    316316      else {
    317         if( ( e = (*RM)[node]) != INVALID ) return G->blueNode(e);
     317        if( ( e = (*RM)[G->asRedNodeUnsafe(node)]) != INVALID ) return G->blueNode(e);
    318318        else return INVALID;
    319319      }
    320320    }
     
    351351    typedef typename WeightMap::Value Value;
    352352   
    353353    /// The type of the red matching map
    354     typedef typename BpGraph::template RedMap<typename BpGraph::Edge>
     354    typedef typename BpGraph::template RedNodeMap<typename BpGraph::Edge>
    355355    RedMatchingMap;
    356356    /// The type of the blue matching map
    357     typedef typename BpGraph::template BlueMap<typename BpGraph::Edge>
     357    typedef typename BpGraph::template BlueNodeMap<typename BpGraph::Edge>
    358358    BlueMatchingMap;
    359359
    360360  private:
     
    420420      BM = new BlueMatchingMap(*G, INVALID);
    421421     
    422422      int nodeCount = 0;
    423       for( RedIt n(*G); n != INVALID; ++n ){
     423      for( RedNodeIt n(*G); n != INVALID; ++n ){
    424424        Value max_C = 0;
    425425        bool start = true;
    426426        vector<Edge> maxEdges;
     
    446446        }
    447447        nodeCount++;
    448448      }
    449       for( BlueIt n(*G); n != INVALID; ++n ){
     449      for( BlueNodeIt n(*G); n != INVALID; ++n ){
    450450        Value minDif = 0;
    451451        bool start = true;
    452452        vector<Edge> minEdges;
     
    482482    HungAlgStatus step()
    483483    {
    484484      bool perfect = true;
    485       for( BlueIt n(*G); n != INVALID; ++n ){
     485      for( BlueNodeIt n(*G); n != INVALID; ++n ){
    486486        if( (*BM)[n] == INVALID ) {
    487487          perfect = false;
    488488          break;
    489489        }
    490490      }
    491       for( RedIt n(*G); n != INVALID; ++n ){
     491      for( RedNodeIt n(*G); n != INVALID; ++n ){
    492492        if( (*RM)[n] == INVALID ) {
    493493          perfect = false;
    494494          break;
    495495        }
    496496      }
    497497      if( perfect ) return HUNGALG_SUCCESS;
    498       queue<Node> sources;
     498      queue<RedNode> sources;
    499499      BoolNodeMap _processed(*G, false);
    500500      EdgeNodeMap _pred(*G, INVALID);
    501501      BoolNodeMap _reached(*G, false);
    502       vector<Node> H;
     502      vector<RedNode> H;
    503503      bool isAlternatePath = false;
    504       for( RedIt n(*G); n != INVALID; ++n ){
     504      for( RedNodeIt n(*G); n != INVALID; ++n ){
    505505        if( (*RM)[n] == INVALID ){
    506506          sources.push( n );
    507507        }
    508508      }
    509509     
    510510      while ( !sources.empty() && !isAlternatePath ){
    511         Node s = sources.front();
     511        RedNode s = sources.front();
    512512        if(_processed[s]) {sources.pop(); continue;}
    513513        _processed[s] = true;
    514514        _reached[s] = true;
    515515        H.push_back( s );
    516         Node m;
     516        BlueNode m;
    517517        for(IncEdgeIt e(*G,s);e!=INVALID;++e) {
    518518          if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
    519519          if(!_reached[m=G->blueNode(e)]) {
    520520            _reached[m]=true;
    521521            _pred[m]=e;
    522522            if( (*BM)[m] != INVALID ) {
    523               Node n;
     523              RedNode n;
    524524              if (!_reached[n = G->redNode((*BM)[m])] ) {
    525525                sources.push( n );
    526526                _pred[n] = (*BM)[m];
     
    537537      }
    538538       
    539539      if( !isAlternatePath ) {
    540         vector<Node> GammaH;
    541         BoolBlueMap GammaH_map( *G, false );
    542         for (typename vector<Node>::iterator nit = H.begin();
     540        vector<BlueNode> GammaH;
     541        BoolBlueNodeMap GammaH_map( *G, false );
     542        for (typename vector<RedNode>::iterator nit = H.begin();
    543543             nit!=H.end(); ++nit){
    544544          for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
    545545            if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
    546             Node n = G->blueNode( e );
     546            BlueNode n = G->blueNode( e );
    547547            if( !GammaH_map[n] ){
    548548              GammaH.push_back(n);
    549549              GammaH_map[n] = true;
     
    552552        }
    553553        Value gamma = 0;
    554554        bool start = true;
    555         for (typename vector<Node>::iterator nit = H.begin();
     555        for (typename vector<RedNode>::iterator nit = H.begin();
    556556             nit!=H.end(); ++nit){
    557557          for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
    558             Node n = G->oppositeNode( *nit, e );
     558            BlueNode n = G->oppositeNode( *nit, e );
    559559            if( !GammaH_map[ n ] ){
    560560              Value dif = PI[ G->u(e) ] + PI[ G->v(e) ] - (*c)[e];
    561561              if( start ){ gamma=dif; start=false; continue; }
     
    567567           
    568568        if(!gamma) return HUNGALG_ERROR;
    569569           
    570         for (typename vector<Node>::iterator nit = H.begin();
     570        for (typename vector<RedNode>::iterator nit = H.begin();
    571571             nit!=H.end(); ++nit){
    572572          PI[*nit] -= gamma;
    573573        }
    574         for (typename vector<Node>::iterator nit = GammaH.begin();
     574        for (typename vector<BlueNode>::iterator nit = GammaH.begin();
    575575             nit!=GammaH.end(); ++nit){
    576576          PI[*nit] += gamma;
    577577        }
     
    631631    {
    632632      Value retVal = 0;
    633633      Edge e;
    634       for( BlueIt n(*G); n != INVALID; ++n )
     634      for( BlueNodeIt n(*G); n != INVALID; ++n )
    635635        if( ( e = (*BM)[n] ) != INVALID ) retVal += (*c)[e];
    636636      return retVal;
    637637    }
     
    656656    /// \pre Either run() or start() must be called before using this function.
    657657    Edge matching(const Node& node) const
    658658    {
    659       return G->blue(node) ? (*BM)[node] : (*RM)[node];
     659      return G->blue(node) ? (*BM)[G->asBlueNodeUnsafe(node)] : (*RM)[G->asRedNodeUnsafe(node)];
    660660    }
    661661
    662662    /// \brief Return a const reference to the blue matching map.
     
    687687    {
    688688      Edge e;
    689689      if( G->blue(node) ){
    690         if( ( e = (*BM)[node]) != INVALID ) return G->redNode(e);
     690        if( ( e = (*BM)[G->asBlueNodeUnsafe(node)]) != INVALID ) return G->redNode(e);
    691691        else return INVALID;
    692692      }
    693693      else {
    694         if( ( e = (*RM)[node]) != INVALID ) return G->blueNode(e);
     694        if( ( e = (*RM)[G->asRedNodeUnsafe(node)]) != INVALID ) return G->blueNode(e);
    695695        else return INVALID;
    696696      }
    697697    }
     
    738738    /// The value type of the edge weights
    739739    typedef typename WeightMap::Value Value;
    740740    /// The type of the red matching map
    741     typedef typename BpGraph::template RedMap<typename BpGraph::Edge>
     741    typedef typename BpGraph::template RedNodeMap<typename BpGraph::Edge>
    742742    RedMatchingMap;
    743743    /// The type of the blue matching map
    744     typedef typename BpGraph::template BlueMap<typename BpGraph::Edge>
     744    typedef typename BpGraph::template BlueNodeMap<typename BpGraph::Edge>
    745745    BlueMatchingMap;
    746746   
    747747  private:
     
    750750   
    751751    typedef typename BpGraph::template NodeMap<Value> NodeMap;
    752752   
    753     typedef typename BpGraph::template RedMap<BlueNode> BlueRedMap;
     753    typedef typename BpGraph::template RedNodeMap<BlueNode> BlueRedMap;
    754754   
    755     typedef typename BpGraph::template BlueMap<RedNode> RedBlueMap;
     755    typedef typename BpGraph::template BlueNodeMap<RedNode> RedBlueMap;
    756756
    757757    typedef typename BpGraph::template NodeMap<Edge> EdgeNodeMap;
    758758   
     
    770770    RedBlueMap *RBM;
    771771    BlueRedMap *BRM;
    772772   
    773     void _processMins( bool &processedMins, queue<Node> &sources, NodeNodeMap &_pred, BoolNodeMap &_reached, BlueNode &first )
     773    void _processMins( bool &processedMins, queue<RedNode> &sources, NodeNodeMap &_pred, BoolNodeMap &_reached, BlueNode &first )
    774774    {
    775775      if( processedMins ) return;
    776776      processedMins = true;
    777       for( BlueIt m(*G); m!=INVALID; ++m ) {
     777      for( BlueNodeIt m(*G); m!=INVALID; ++m ) {
    778778        if( PI[m] == -MinPI ){
    779779          _reached[m]=true;
    780780          if( (*RBM)[m] != INVALID ) {
    781             Node n;
     781            RedNode n;
    782782            if (!_reached[n = (*RBM)[m]] ) {
    783783              sources.push( n );
    784784              _pred[n] = m;
     
    788788        }
    789789      }
    790790       
    791       for( RedIt m(*G); m!=INVALID; ++m ) {
     791      for( RedNodeIt m(*G); m!=INVALID; ++m ) {
    792792        if( PI[m] == MinPI && !_reached[ m ] && (*BRM)[m] == INVALID ) {
    793793          sources.push( m );
    794794          _pred[m] = INVALID;
     
    798798     
    799799    }
    800800   
    801     void _setMatching( NodeNodeMap &pred, Node n, BlueNode &first ) {
    802       Node node = n;
    803       Node prevNode = pred[n];
     801    void _setMatching( NodeNodeMap &pred, BlueNode n, BlueNode &first ) {
     802      BlueNode node = n;
     803      RedNode prevNode = G->asRedNodeUnsafe(pred[n]);
    804804      while( 1 ) {
    805805        if( prevNode == INVALID ){
    806806          if( PI[node] == -MinPI ){
     
    811811        else if( (*RBM)[node] != prevNode ) {
    812812          (*BRM)[prevNode] = node;
    813813          (*RBM)[node] = prevNode;
    814           node = pred[prevNode];
     814          node = G->asBlueNodeUnsafe(pred[prevNode]);
    815815        }
    816816        if( node == INVALID ) {
    817817          if( ( PI[prevNode] == MinPI && (*BRM)[prevNode] == INVALID ) || first == INVALID ) return;
    818818          else node = first;
    819819        }
    820         prevNode = pred[node];
     820        prevNode = G->asRedNodeUnsafe(pred[node]);
    821821      }
    822822    }
    823823   
     
    858858      BRM = new BlueRedMap(*G, INVALID);
    859859      RBM = new RedBlueMap(*G, INVALID);
    860860     
    861       for( BlueIt n(*G); n != INVALID; ++n )
     861      for( BlueNodeIt n(*G); n != INVALID; ++n )
    862862      {
    863863        Value max_C = 0;
    864         vector<Node> maxNodes;
     864        vector<RedNode> maxNodes;
    865865        for( IncEdgeIt e( *G, n ); e!=INVALID; ++e ) {
    866866          if( (*c)[e] > max_C ) {
    867867            max_C = (*c)[e];
     
    873873          }
    874874        }
    875875        PI[n] = max_C;
    876         for (typename vector<Node>::iterator eit = maxNodes.begin();
     876        for (typename vector<RedNode>::iterator eit = maxNodes.begin();
    877877          eit!=maxNodes.end(); ++eit)
    878878        {
    879879          if( (*BRM)[*eit] == INVALID )
     
    884884          }
    885885        }
    886886      }
    887       for( RedIt n(*G); n != INVALID; ++n ){
     887      for( RedNodeIt n(*G); n != INVALID; ++n ){
    888888        PI[n] = 0;
    889889      }
    890890    }
     
    901901      for( NodeIt n(*G); n != INVALID; ++n ){
    902902       
    903903        if( G->blue(n) ) {
    904           if( (*RBM)[n] == INVALID && PI[n] != -MinPI ){
     904          if( (*RBM)[G->asBlueNodeUnsafe(n)] == INVALID && PI[n] != -MinPI ){
    905905            Continue = true;
    906906          }
    907907        }
    908908        else {
    909           if( (*BRM)[n] == INVALID && PI[n] != MinPI ){
     909          if( (*BRM)[G->asRedNodeUnsafe(n)] == INVALID && PI[n] != MinPI ){
    910910            Continue = true;
    911911          }
    912912        }
    913913      }
    914914      if( !Continue ){
    915         for( RedIt n(*G); n != INVALID; ++n ) {
     915        for( RedNodeIt n(*G); n != INVALID; ++n ) {
    916916          PI[n] -= MinPI;
    917917         
    918918          for( IncEdgeIt e(*G,n); e != INVALID; ++e ) {
    919919            if( (*BRM)[n] == G->blueNode(e) ) (*RM)[n] = e;
    920920          }
    921921        }
    922         for( BlueIt n(*G); n != INVALID; ++n ) {
     922        for( BlueNodeIt n(*G); n != INVALID; ++n ) {
    923923          PI[n] += MinPI;
    924924         
    925925          for( IncEdgeIt e(*G,n); e != INVALID; ++e ) {
     
    929929       
    930930        return true;
    931931      }
    932       queue<Node> sources;
     932      queue<RedNode> sources;
    933933      BoolNodeMap _processed(*G, false);
    934934      NodeNodeMap _pred(*G, INVALID);
    935935      BoolNodeMap _reached(*G, false);
    936       vector<Node> H;
     936      vector<RedNode> H;
    937937      bool isAlternatePath = false;
    938938      bool notMinUncovered = false;
    939939      BlueNode first = INVALID;
    940940     
    941       for( BlueIt n(*G); n != INVALID; ++n ){
     941      for( BlueNodeIt n(*G); n != INVALID; ++n ){
    942942        if( (*RBM)[n] == INVALID && PI[n] != -MinPI ){
    943943          notMinUncovered = true;
    944944          break;
    945945        }
    946946      }
    947947     
    948       for( RedIt n(*G); n != INVALID; ++n ){
     948      for( RedNodeIt n(*G); n != INVALID; ++n ){
    949949        if( (*BRM)[n] == INVALID ){
    950950          if( notMinUncovered || PI[n] != MinPI ) {
    951951            sources.push( n );
     
    961961      }
    962962     
    963963      while ( !sources.empty() && !isAlternatePath ){
    964         Node s = sources.front();
     964        RedNode s = sources.front();
    965965        if(_processed[s]) { sources.pop(); continue; }
    966966        if(PI[s] == MinPI && (*BRM)[s] == INVALID ) _processMins( processedMins, sources, _pred, _reached, first );
    967967        _processed[s] = true;
    968968        _reached[s] = true;
    969969        H.push_back( s );
    970         Node m;
     970        BlueNode m;
    971971        for( IncEdgeIt e(*G,s); e!=INVALID; ++e ) {
    972972          if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
    973973          if( _reached[m=G->blueNode(e)] ) continue;
    974974          _reached[m]=true;
    975975          _pred[m]=s;
    976976          if( (*RBM)[m] != INVALID ) {
    977             Node n;
     977            RedNode n;
    978978            if (!_reached[n = (*RBM)[m]] ) {
    979979              sources.push( n );
    980980              _pred[n] = m;
     
    994994      }
    995995     
    996996      if( !isAlternatePath ) {
    997         vector<Node> GammaH;
    998         BoolBlueMap GammaH_map( *G, false );
    999         for (typename vector<Node>::iterator nit = H.begin();
     997        vector<BlueNode> GammaH;
     998        BoolBlueNodeMap GammaH_map( *G, false );
     999        for (typename vector<RedNode>::iterator nit = H.begin();
    10001000             nit!=H.end(); ++nit){
    10011001          for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
    10021002            if( (*c)[e] != PI[G->blueNode(e)] + PI[G->redNode(e)] ) continue;
    1003             Node n = G->blueNode( e );
     1003            BlueNode n = G->blueNode( e );
    10041004            if( !GammaH_map[n] ){
    10051005              GammaH.push_back(n);
    10061006              GammaH_map[n] = true;
     
    10081008          }
    10091009        }
    10101010        if( processedMins ) {
    1011           for( BlueIt n(*G); n != INVALID; ++n ) {
     1011          for( BlueNodeIt n(*G); n != INVALID; ++n ) {
    10121012            if( PI[n] == -MinPI && !GammaH_map[n] ) {
    10131013              GammaH.push_back(n);
    10141014              GammaH_map[n] = true;
     
    10181018       
    10191019        Value gamma = 0;
    10201020        bool start = true;
    1021         for (typename vector<Node>::iterator nit = H.begin();
     1021        for (typename vector<RedNode>::iterator nit = H.begin();
    10221022             nit!=H.end(); ++nit){
    10231023          for( IncEdgeIt e( *G, *nit ); e!=INVALID; ++e ){
    10241024            if((*c)[e]<0) continue;
    1025             Node n = G->oppositeNode( *nit, e );
     1025            BlueNode n = G->asBlueNodeUnsafe(G->oppositeNode( *nit, e ));
    10261026            if( !GammaH_map[ n ] ){
    10271027              Value dif = PI[ G->u(e) ] + PI[ G->v(e) ] - (*c)[e];
    10281028              if( start ){ gamma=dif; start=false; continue; }
     
    10331033        }
    10341034       
    10351035        if( processedMins ) {
    1036           for( BlueIt n(*G); n != INVALID; ++n ) {
     1036          for( BlueNodeIt n(*G); n != INVALID; ++n ) {
    10371037            if( !GammaH_map[n] && ( PI[n] + MinPI < gamma || !gamma ) ) gamma = PI[n] + MinPI;
    10381038          }
    10391039        }
    10401040           
    1041         for (typename vector<Node>::iterator nit = H.begin();
     1041        for (typename vector<RedNode>::iterator nit = H.begin();
    10421042             nit!=H.end(); ++nit){
    10431043          PI[*nit] -= gamma;
    10441044        }
    10451045        if( processedMins ) MinPI -= gamma;
    1046         for (typename vector<Node>::iterator nit = GammaH.begin();
     1046        for (typename vector<BlueNode>::iterator nit = GammaH.begin();
    10471047             nit!=GammaH.end(); ++nit){
    10481048          PI[*nit] += gamma;
    10491049        }
     
    10971097    {
    10981098      Value retVal = 0;
    10991099      Edge e;
    1100       for( BlueIt n(*G); n != INVALID; ++n ) {
     1100      for( BlueNodeIt n(*G); n != INVALID; ++n ) {
    11011101        if( ( e = (*BM)[n] ) != INVALID ) retVal += (*c)[e];
    11021102      }
    11031103      return retVal;
     
    11111111    int matchingSize() const
    11121112    {
    11131113      int size = 0;
    1114       for ( BlueIt n(*G); n != INVALID; ++n) {
     1114      for ( BlueNodeIt n(*G); n != INVALID; ++n) {
    11151115        if ( (*BM)[n] != INVALID) {
    11161116          ++size;
    11171117        }
     
    11391139    /// \pre Either run() or start() must be called before using this function.
    11401140    Edge matching( const Node& node ) const
    11411141    {
    1142       return G->blue(node) ? (*BM)[node] : (*RM)[node];
     1142      return G->blue(node) ? (*BM)[G->asBlueNodeUnsafe(node)] : (*RM)[G->asRedNodeUnsafe(node)];
    11431143    }
    11441144
    11451145    /// \brief Return a const reference to the blue matching map.
     
    11701170    {
    11711171      Edge e;
    11721172      if( G->blue(node) ){
    1173         if( ( e = (*BM)[node]) != INVALID ) return G->redNode(e);
     1173        if( ( e = (*BM)[G->asBlueNodeUnsafe(node)]) != INVALID ) return G->redNode(e);
    11741174        else return INVALID;
    11751175      }
    11761176      else {
    1177         if( ( e = (*RM)[node]) != INVALID ) return G->blueNode(e);
     1177        if( ( e = (*RM)[G->asRedNodeUnsafe(node)]) != INVALID ) return G->blueNode(e);
    11781178        else return INVALID;
    11791179      }
    11801180    }
  • lemon/lgf_reader.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1328396684 -3600
    # Node ID ce9531988b727c5997b9ff84ad9f65f49bbabced
    # Parent  2cf5b6b4bdcbba2495741f4aa7d6f54b02640dad
    Update LGF reader to work with typesafe bipartite node sets
    
    diff --git a/lemon/lgf_reader.h b/lemon/lgf_reader.h
    a b  
    29532953            throw FormatError(msg.str());
    29542954          }
    29552955
    2956           e = _graph.addEdge(source, target);
     2956          // It is checked that source is red and
     2957          // target is blue, so this should be safe:
     2958          e = _graph.addEdge(_graph.asRedNodeUnsafe(source),
     2959                             _graph.asBlueNodeUnsafe(target));
    29572960          if (label_index != -1)
    29582961            _edge_index.insert(std::make_pair(tokens[label_index], e));
    29592962        } else {
  • lemon/elevator.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1329865937 -3600
    # Node ID 5cc53dbced80bc10b013763fcc17c840fa5e736c
    # Parent  ce9531988b727c5997b9ff84ad9f65f49bbabced
    Fix initialization bug in elevator
    
    diff --git a/lemon/elevator.h b/lemon/elevator.h
    a b  
    22 *
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2009
     5 * Copyright (C) 2003-2012
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
    88 *
     
    108108    Elevator(const GR &graph,int max_level) :
    109109      _g(graph),
    110110      _max_level(max_level),
    111       _item_num(_max_level),
     111      _item_num(countItems<GR, Item>(graph)),
    112112      _where(graph),
    113113      _level(graph,0),
    114       _items(_max_level),
     114      _items(_item_num),
    115115      _first(_max_level+2),
    116116      _last_active(_max_level+2),
    117117      _highest_active(-1) {}
     
    526526    ///\param max_level The maximum allowed level.
    527527    ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
    528528    LinkedElevator(const GR& graph, int max_level)
    529       : _graph(graph), _max_level(max_level), _item_num(_max_level),
     529      : _graph(graph), _max_level(max_level),
     530        _item_num(countItems<GR, Item>(graph)),
    530531        _first(_max_level + 1), _last(_max_level + 1),
    531532        _prev(graph), _next(graph),
    532533        _highest_active(-1), _level(graph), _active(graph) {}
  • new file lemon/hopcroft_karp.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1330302130 -3600
    # Node ID fd8741fc21438ea3c829a7d4f6af50d92b7ce426
    # Parent  5cc53dbced80bc10b013763fcc17c840fa5e736c
    Hopcroft-Karp bipartite matching
    
    diff --git a/lemon/hopcroft_karp.h b/lemon/hopcroft_karp.h
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_HOPCROFT_KARP_H
     20#define LEMON_HOPCROFT_KARP_H
     21
     22#include <lemon/core.h>
     23#include <vector>
     24#include <queue>
     25#include <list>
     26
     27/// \ingroup matching
     28/// \file
     29/// \brief Hopcroft-Karp algorithm.
     30namespace lemon {
     31
     32  /// \brief The Hopcroft-Karp bipartite matching algorithm
     33  ///
     34  /// Finds maximal matching in a given bipartite
     35  /// graph using the Hopcroft-Karp algorithm,
     36  /// having \f$O(e\sqrt{n})\f$ complexity.
     37  template <typename BPG>
     38  class HopcroftKarp {
     39  public:
     40    /// Type of the bipartite graph
     41    typedef BPG BpGraph;
     42    /// Type of the matching map
     43    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     44  private:
     45    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     46
     47    const BpGraph& _bpg;
     48    MatchingMap* _matching;
     49    bool _local_matching;
     50
     51  protected:
     52    void createStructures() {
     53      if (!_matching) {
     54        _matching = new MatchingMap(_bpg, INVALID);
     55        _local_matching = true;
     56      }
     57    }
     58
     59    void destroyStructures() {
     60      if (_local_matching) {
     61        delete _matching;
     62      }
     63    }
     64
     65    HopcroftKarp() {}
     66
     67  public:
     68    /// \brief Constructor
     69    ///
     70    /// Constructs the class on the given bipartite graph.
     71    HopcroftKarp(const BpGraph& bpg) : _bpg(bpg),
     72                                       _matching(0),
     73                                       _local_matching(false)
     74    {}
     75
     76
     77    /// \brief Destructor
     78    ///
     79    /// Destructor.
     80    ~HopcroftKarp() {
     81      destroyStructures();
     82    }
     83
     84    /// \brief Sets the matching map
     85    ///
     86    /// Sets the matching map. It should contain a valid matching,
     87    /// for example the empty matching is fine.
     88    /// If you don't use this function before calling \ref init(),
     89    /// an instance will be allocated automatically.
     90    /// The destructor deallocates this automatically allocated map,
     91    /// of course.
     92    ///
     93    /// \return <tt>(*this)</tt>
     94    HopcroftKarp& matchingMap(MatchingMap& map) {
     95      _matching = &map;
     96      _local_matching = false;
     97      return *this;
     98    }
     99
     100    /// \brief Returns a const reference to the matching map
     101    ///
     102    /// Returns a const reference to the matching map, which contains
     103    /// the matching edge for every node (and \c INVALID for the
     104    /// unmatched nodes).
     105    const MatchingMap& matchingMap() const {
     106      return *_matching;
     107    }
     108
     109
     110    /// \brief Initializes the algorithm
     111    ///
     112    /// Allocates the matching map if necessary.
     113    ///
     114    /// \pre \ref init() is not called.
     115    void init() {
     116      createStructures();
     117    }
     118
     119    /// \brief Smarter initialization of the matching.
     120    ///
     121    /// Allocates the matching map if necessary, and initializes
     122    /// the algorithm with a trivially found matching: iterating
     123    /// on every edge, take it in if both endnodes are unmatched.
     124    ///
     125    /// \return Size of the found matching.
     126    ///
     127    /// \note heuristicInit() replaces init() regarding the preconditions.
     128    int heuristicInit() {
     129      createStructures();
     130      int size = 0;
     131
     132      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     133        bool matched = false;
     134        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
     135          if ((*_matching)[_bpg.target(oi)] == INVALID) {
     136            _matching->set(_bpg.source(oi), oi);
     137            _matching->set(_bpg.target(oi), oi);
     138            matched = true;
     139            ++size;
     140          }
     141        }
     142      }
     143
     144      return size;
     145    }
     146
     147    /// \brief Initialize the matching from a map.
     148    ///
     149    /// Allocates the matching map if necessary, and
     150    /// initializes the algorithm with a matching given as a bool edgemap,
     151    /// in which an edge is true if it is in the matching.
     152    ///
     153    /// If the given matching is invalid, some edges will be left out.
     154    ///
     155    /// \return \c false if the given matching is invalid.
     156    ///
     157    /// \note matchingInit() replaces init() regarding the preconditions.
     158    bool matchingInit(const BoolEdgeMap& matching) {
     159      createStructures();
     160      bool valid = true;
     161      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
     162        if (matching[it]) {
     163          Node red = _bpg.redNode(it);
     164          Node blue = _bpg.blueNode(it);
     165          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
     166            valid = false;
     167          } else {
     168            _matching->set(red, it);
     169            _matching->set(blue, it);
     170          }
     171        }
     172      }
     173      return valid;
     174    }
     175
     176    /// \brief Executes an augmenting phase
     177    ///
     178    /// Searches a maximal set of vertex disjoint shortest alternating paths.
     179    /// Meaning:
     180    ///  - alternating: connents an unmatched red- and an unmatched blue node,
     181    ///    and exactly every second edge is in the current matching;
     182    ///  - shortest: contain a minimal number of edges;
     183    ///  - vertex disjoint: every vertex belong to at most one path;
     184    ///  - maximal set: a set of path is maximal when it is not expandable
     185    ///    further (and not when there does not exist a set with
     186    ///    more vertex disjoint shortest alternating paths).
     187    ///
     188    /// After a maximal set is found, it applies the augmenting paths,
     189    /// so edges of the matching are taken out, the others are put in
     190    /// the matching.
     191    ///
     192    /// \return The length of the found augmenting paths, or -1 when none
     193    /// found, which occurs if and only if the current matching is maximal.
     194    ///
     195    /// \pre \ref init() must be called before using this function.
     196    int augment() {
     197      IntBlueNodeMap level(_bpg, -2);
     198      std::vector<RedNode> act_red_nodes;
     199
     200      for (RedNodeIt r_it(_bpg); r_it != INVALID; ++r_it) {
     201        if ((*_matching)[r_it] == INVALID) {
     202          act_red_nodes.push_back(r_it);
     203        }
     204      }
     205      // This vector will contain the end nodes of the possible
     206      // augmenting paths.
     207      std::vector<BlueNode> path_ends;
     208
     209      // Layer counter
     210      int cur_level = -1;
     211
     212      // Starting from the unmatched red nodes search for unmatched
     213      // blue nodes, using Bfs (but only on alternating paths).
     214      // Blue nodes with level==-2 are unreached.
     215      while (path_ends.empty()) {
     216        cur_level += 2;
     217        std::vector<RedNode> next_red_nodes;
     218        for (typename std::vector<RedNode>::iterator r=act_red_nodes.begin();
     219             r != act_red_nodes.end(); ++r) {
     220          for (OutArcIt it(_bpg, *r); it!=INVALID; ++it) {
     221            BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     222            if (level[blue] != -2) continue;
     223            level[blue] = cur_level;
     224            if ((*_matching)[blue] == INVALID) {
     225              path_ends.push_back(blue);
     226            } else {
     227              next_red_nodes
     228                .push_back(_bpg.redNode((*_matching)[blue]));
     229            }
     230          }
     231        }
     232
     233        if (path_ends.empty() && next_red_nodes.empty())
     234          return -1;
     235
     236        act_red_nodes.swap(next_red_nodes);
     237        next_red_nodes.clear();
     238      }
     239      // length of found augmenting paths
     240      int max_level = cur_level;
     241
     242      // Stack for the Dfs
     243      std::vector<OutArcIt> stack;
     244      // Raise this flag when a shortest augmenting path is found.
     245      bool path_found;
     246      // Searching backward, starting from the previously found
     247      // blue nodes, we build vertex disjoint alternating paths.
     248      // Blue nodes with level!=-2 are unreached.
     249      typename std::vector<BlueNode>::iterator end = path_ends.begin();
     250
     251      while (end != path_ends.end()) {
     252        stack.clear();
     253        stack.push_back(OutArcIt(_bpg, *end));
     254        cur_level = max_level;
     255        path_found = false;
     256
     257        while (!stack.empty() && !path_found) {
     258          OutArcIt &o = stack.back();
     259          // Search an arc that goes one level higher.
     260          // If none found, retreat.
     261          while (o != INVALID &&
     262                 (*_matching)[_bpg.target(o)] != INVALID &&
     263                 level[_bpg.blueNode((*_matching)[_bpg.target(o)])] !=
     264                 cur_level - 2) {
     265            ++o;
     266          }
     267          if (o==INVALID) {
     268            stack.pop_back();
     269            cur_level += 2;
     270          } else if ((*_matching)[_bpg.target(o)] == INVALID) {
     271            path_found = true;
     272          } else {
     273            BlueNode b = _bpg.blueNode((*_matching)[_bpg.target(o)]);
     274            stack.push_back(OutArcIt(_bpg, b));
     275            level[b] = -2;        // Do not visit this node again
     276            cur_level -= 2;
     277          }
     278        }
     279        if (path_found) {
     280          OutArcIt o;
     281          while (!stack.empty()) {
     282            o = stack.back();
     283            _matching->set(_bpg.source(o), o);
     284            _matching->set(_bpg.target(o), o);
     285            stack.pop_back();
     286          }
     287        }
     288        ++end;
     289      }
     290      return max_level;
     291    }
     292
     293    /// \brief Executes the algorithm
     294    ///
     295    /// It runs augmenting phases until the optimal solution is reached.
     296    ///
     297    /// \pre \ref init() must be called before using this function.
     298    void start() {
     299      while (augment() > 0) {}
     300    }
     301
     302    /// \brief Runs the algorithm.
     303    ///
     304    /// hk.run() is just a shorthand for:
     305    ///
     306    ///\code
     307    /// hk.heuristicInit();
     308    /// hk.start();
     309    ///\endcode
     310    void run() {
     311      heuristicInit();
     312      start();
     313    }
     314
     315    /// \brief Size of the matching
     316    ///
     317    /// Returns the size of the current matching.
     318    /// Function runs in \f$O(n)\f$ time.
     319    int matchingSize() const {
     320      int size = 0;
     321      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
     322        if ((*_matching)[it] != INVALID) ++size;
     323      }
     324      return size;
     325    }
     326
     327    /// \brief Return \c true if the given edge is in the matching.
     328    ///
     329    /// This function returns \c true if the given edge is in the current
     330    /// matching.
     331    bool matching(const Edge& edge) const {
     332      return edge == (*_matching)[_bpg.redNode(edge)];
     333    }
     334
     335    /// \brief Return the matching edge incident to the given node.
     336    ///
     337    /// This function returns the matching edge incident to the
     338    /// given node in the current matching or \c INVALID if the node is
     339    /// not covered by the matching.
     340    Edge matching(const Node& n) const {
     341      return (*_matching)[n];
     342    }
     343
     344    /// \brief Return the mate of the given node.
     345    ///
     346    /// This function returns the mate of the given node in the current
     347    /// matching or \c INVALID if the node is not covered by the matching.
     348    Node mate(const Node& n) const {
     349      return (*_matching)[n] != INVALID ?
     350        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     351    }
     352
     353    /// \brief Query a vertex cover.
     354    ///
     355    /// This function finds a vertex cover based on the current matching,
     356    /// and in \c cover sets node in the cover \c true and \c false
     357    /// otherwise. In case of maximal matching, this will be a minimal
     358    /// vertex cover.
     359    /// Function runs in \f$O(n + e)\f$ time.
     360    ///
     361    /// \return The number of nodes in the cover.
     362    int cover(BoolNodeMap &cover) const {
     363      int count = 0;
     364      std::queue<RedNode> q;
     365
     366      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
     367        cover[b] = false;
     368
     369      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     370        if ((*_matching)[r] == INVALID) {
     371          cover[r] = false;
     372          q.push(r);
     373        } else {
     374          cover[r] = true;
     375          ++count;
     376        }
     377      }
     378
     379      while (!q.empty()) {
     380        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     381          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     382          if (cover[blue]) continue;
     383
     384          cover[blue] = true;
     385          ++count;
     386
     387          if ((*_matching)[blue] != INVALID) {
     388            RedNode nr = _bpg.redNode((*_matching)[blue]);
     389            q.push(nr);
     390            cover[nr] = false;
     391            --count;
     392          }
     393        }
     394        q.pop();
     395      }
     396
     397      return count;
     398    }
     399
     400    /// \brief Query a barrier of red nodes.
     401    ///
     402    /// This function tries to create a barrier of red nodes, which should:
     403    ///  - contain every unmatched red nodes,
     404    ///  - contain exactly that many matched red nodes as much blue nodes
     405    ///    there are in its neighbor set.
     406    ///
     407    /// Barrier exist if and only if the matching is maximal. If exist, the
     408    /// map is set true for nodes of the barrier, and false for the others.
     409    ///
     410    /// Function runs in \f$O(n + e)\f$ time.
     411    ///
     412    /// \return \c true if a barrier found.
     413    bool redBarrier(BoolRedNodeMap &barrier) const {
     414      bool exist = true;
     415      std::queue<RedNode> q;
     416
     417      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     418        if ((*_matching)[r] == INVALID) {
     419          barrier[r] = true;
     420          q.push(r);
     421        } else {
     422          barrier[r] = false;
     423        }
     424      }
     425
     426      while (!q.empty() && exist) {
     427        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     428          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     429
     430          if (exist &= ((*_matching)[blue] != INVALID)) {
     431            RedNode r = _bpg.redNode((*_matching)[blue]);
     432            if (!barrier[r]) {
     433              q.push(r);
     434              barrier[r] = true;
     435            }
     436          }
     437        }
     438        q.pop();
     439      }
     440
     441      return exist;
     442    }
     443
     444    /// \brief Query a barrier of blue nodes.
     445    ///
     446    /// This function tries to create a barrier of blue nodes, which should:
     447    ///  - contain every unmatched blue nodes,
     448    ///  - contain exactly that many matched blue nodes as much red nodes
     449    ///    there are in its neighbor set.
     450    ///
     451    /// Barrier exist if and only if the matching is maximal. If exist, the
     452    /// map is set true for nodes of the barrier, and false for the others.
     453    ///
     454    /// Function runs in \f$O(n + e)\f$ time.
     455    ///
     456    /// \return \c true if a barrier found.
     457    bool blueBarrier(BoolBlueNodeMap &barrier) const {
     458      bool exist = true;
     459      std::queue<BlueNode> q;
     460
     461      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     462        if ((*_matching)[b] == INVALID) {
     463          barrier[b] = true;
     464          q.push(b);
     465        } else {
     466          barrier[b] = false;
     467        }
     468      }
     469
     470      while (!q.empty() && exist) {
     471        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     472          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
     473
     474          if (exist &= ((*_matching)[r] != INVALID)) {
     475            BlueNode b = _bpg.blueNode((*_matching)[r]);
     476            if (!barrier[b]) {
     477              q.push(b);
     478              barrier[b] = true;
     479            }
     480          }
     481        }
     482        q.pop();
     483      }
     484
     485      return exist;
     486    }
     487  };
     488
     489}
     490#endif
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    3434  graph_utils_test
    3535  hao_orlin_test
    3636  heap_test
     37  hopcroft_karp_test
    3738  kruskal_test
    3839  lgf_test
    3940  maps_test
  • test/Makefile.am

    diff --git a/test/Makefile.am b/test/Makefile.am
    a b  
    3232        test/graph_utils_test \
    3333        test/hao_orlin_test \
    3434        test/heap_test \
     35        test/hopcroft_karp_test \
    3536        test/kruskal_test \
    3637        test/lgf_test \
    3738        test/maps_test \
     
    8788test_graph_utils_test_SOURCES = test/graph_utils_test.cc
    8889test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
    8990test_heap_test_SOURCES = test/heap_test.cc
     91test_hopcroft_karp_test_SOURCES = test/hopcroft_karp_test.cc
    9092test_kruskal_test_SOURCES = test/kruskal_test.cc
    9193test_lgf_test_SOURCES = test/lgf_test.cc
    9294test_lp_test_SOURCES = test/lp_test.cc
  • new file test/hopcroft_karp_test.cc

    diff --git a/test/hopcroft_karp_test.cc b/test/hopcroft_karp_test.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <list>
     22
     23#include <lemon/random.h>
     24#include <lemon/hopcroft_karp.h>
     25#include <lemon/smart_graph.h>
     26#include <lemon/concepts/bpgraph.h>
     27#include <lemon/concepts/maps.h>
     28#include <lemon/lgf_reader.h>
     29
     30#include "test_tools.h"
     31
     32using namespace std;
     33using namespace lemon;
     34
     35const int lgfn = 3;
     36const string lgf[lgfn] = {
     37  // Empty graph
     38  "@red_nodes\n"
     39  "label\n"
     40  "\n"
     41  "@blue_nodes\n"
     42  "label\n"
     43  "\n"
     44  "@edges\n"
     45  "  label\n"
     46  "\n",
     47
     48  // No red nodes
     49  "@red_nodes\n"
     50  "label\n"
     51  "\n"
     52  "@blue_nodes\n"
     53  "label\n"
     54  "1\n"
     55  "@edges\n"
     56  "  label\n"
     57  "\n",
     58
     59  // No blue nodes
     60  "@red_nodes\n"
     61  "label\n"
     62  "1\n"
     63  "@blue_nodes\n"
     64  "label\n"
     65  "\n"
     66  "@edges\n"
     67  "  label\n"
     68  "\n",
     69};
     70
     71const string u_lgf =
     72  "@red_nodes\n"
     73  "label\n"
     74  "1\n"
     75  "2\n"
     76  "3\n"
     77  "@blue_nodes\n"
     78  "label\n"
     79  "4\n"
     80  "5\n"
     81  "6\n"
     82  "7\n"
     83  "@edges\n"
     84  "label\n"
     85  "1 4 1\n"
     86  "1 5 2\n"
     87  "2 5 3\n"
     88  "2 6 4\n"
     89  "2 7 5\n"
     90  "3 7 6\n"
     91;
     92
     93
     94void checkHopcroftKarpCompile() {
     95  typedef concepts::BpGraph BpGraph;
     96  BPGRAPH_TYPEDEFS(BpGraph);
     97  typedef HopcroftKarp<BpGraph> HK;
     98
     99  BpGraph graph;
     100  RedNode red;
     101  BlueNode blue;
     102  Node node;
     103  Edge edge;
     104  BoolEdgeMap init_matching(graph);
     105  HK hk_test(graph);
     106  const HK& const_hk_test = hk_test;
     107  HK::MatchingMap matching_map(graph, INVALID);
     108  BoolNodeMap cov(graph);
     109  int cover_size;
     110  BoolRedNodeMap rb(graph);
     111  BoolBlueNodeMap bb(graph);
     112
     113  hk_test.matchingMap(matching_map);
     114  hk_test.init();
     115  hk_test.heuristicInit();
     116  hk_test.matchingInit(init_matching);
     117
     118  hk_test.start();
     119  hk_test.augment();
     120  hk_test.run();
     121
     122  const_hk_test.matchingSize();
     123  const_hk_test.matching(node);
     124  const_hk_test.matching(edge);
     125  const_hk_test.mate(red);
     126  const_hk_test.mate(blue);
     127  const_hk_test.mate(node);
     128  const HK::MatchingMap& max_matching = const_hk_test.matchingMap();
     129  edge = max_matching[node];
     130  cover_size = const_hk_test.cover(cov);
     131  const_hk_test.redBarrier(rb);
     132  const_hk_test.blueBarrier(bb);
     133}
     134
     135typedef SmartBpGraph BpGraph;
     136typedef HopcroftKarp<BpGraph> HK;
     137BPGRAPH_TYPEDEFS(BpGraph);
     138
     139void checkMatchingValidity(const BpGraph& bpg, const HK& hk) {
     140  BoolBlueNodeMap matched(bpg, false);
     141  for (RedNodeIt it(bpg); it != INVALID; ++it) {
     142    if (hk.matching(it) != INVALID) {
     143      check(hk.mate(hk.mate(it)) == it, "Wrong matching");
     144      matched.set(bpg.asBlueNode(hk.mate(it)), true);
     145    }
     146  }
     147  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
     148    // != is used as logical xor here
     149    check(matched[it] != (hk.matching(it) == INVALID), "Wrong matching");
     150  }
     151
     152}
     153
     154void checkCover(const BpGraph& bpg, const HK& hk) {
     155  int cov_size;
     156  BoolNodeMap cov(bpg);
     157  cov_size = hk.cover(cov);
     158  BoolEdgeMap covered(bpg, false);
     159  for (NodeIt n(bpg); n!=INVALID; ++n) {
     160    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
     161      covered[e] = true;
     162    }
     163  }
     164  bool all = true;
     165  for (EdgeIt e(bpg); e!=INVALID; ++e) {
     166    all&=covered[e];
     167  }
     168
     169  check(all, "Invalid cover");
     170  check(hk.matchingSize() == cov_size, "Suboptimal matching.");
     171}
     172
     173void checkBarriers(const BpGraph& bpg, const HK& hk) {
     174  BoolRedNodeMap rb(bpg);
     175  BoolBlueNodeMap bb(bpg);
     176
     177  check(hk.redBarrier(rb), "Suboptimal matching.");
     178  check(hk.blueBarrier(bb), "Suboptimal matching.");
     179
     180  BoolNodeMap reached(bpg, false);
     181
     182  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     183    if (!rb[r]) continue;
     184    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
     185      reached.set(bpg.target(a), true);
     186    }
     187  }
     188  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     189    if (!bb[b]) continue;
     190    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
     191      reached.set(bpg.target(a), true);
     192    }
     193  }
     194
     195  int ur = 0, ub = 0, rc = 0, bc = 0;
     196  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     197    if (rb[r]) ++rc;
     198    if (hk.mate(r) == INVALID) ++ur;
     199    if (reached[r]) --bc;
     200  }
     201
     202  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     203    if (bb[b]) ++bc;
     204    if (hk.mate(b) == INVALID) ++ub;
     205    if (reached[b]) --rc;
     206  }
     207
     208  check(rc == ur, "Red barrier is wrong.");
     209  check(bc == ub, "Blue barrier is wrong.");
     210}
     211
     212void checkMatching(const BpGraph& bpg, const HK& hk) {
     213  checkMatchingValidity(bpg, hk);
     214  checkBarriers(bpg, hk);
     215  checkCover(bpg, hk);
     216}
     217
     218
     219int main() {
     220  // Checking hard wired extremal graphs
     221  for (int i=0; i<lgfn; ++i) {
     222    BpGraph bpg;
     223    istringstream lgfs(lgf[i]);
     224    bpGraphReader(bpg, lgfs).run();
     225    HK hk(bpg);
     226    hk.run();
     227    checkMatching(bpg, hk);
     228  }
     229
     230  // Checking some use cases
     231  BpGraph bpg;
     232  istringstream u_lgfs(u_lgf);
     233  bpGraphReader(bpg, u_lgfs).run();
     234  {
     235    HK hk(bpg);
     236    hk.init();
     237    hk.augment();
     238    hk.start();
     239    checkMatching(bpg, hk);
     240  }
     241  {
     242    HK hk(bpg);
     243    hk.heuristicInit();
     244    hk.augment();
     245    hk.start();
     246    checkMatching(bpg, hk);
     247  }
     248  {
     249    HK hk(bpg);
     250    HK::MatchingMap m(bpg, INVALID);
     251    hk.matchingMap(m);
     252    hk.run();
     253    checkMatching(bpg, hk);
     254  }
     255  {
     256    HK hk(bpg);
     257    BpGraph::EdgeMap<bool> init_matching(bpg, false);
     258    BpGraph::Edge e;
     259    bpg.first(e);
     260    init_matching[e] = true;
     261    check(hk.matchingInit(init_matching), "Wrong result from matchingInit()");
     262    hk.start();
     263    checkMatching(bpg, hk);
     264  }
     265  {
     266    HK hk(bpg);
     267    HK::MatchingMap m(bpg, INVALID);
     268    hk.matchingMap(m);
     269    BpGraph::EdgeMap<bool> init_matching(bpg, true);
     270    check(!hk.matchingInit(init_matching), "Wrong result from matchingInit()");
     271    hk.start();
     272    checkMatching(bpg, hk);
     273  }
     274  // Checking some random generated graphs
     275  const int random_test_n = 10;
     276  const int max_red_n = 1000;
     277  const int min_red_n = 10;
     278  const int max_blue_n = 1000;
     279  const int min_blue_n = 10;
     280  for (int i=0; i<random_test_n; ++i) {
     281    int red_n = rnd.integer(min_red_n, max_red_n);
     282    int blue_n = rnd.integer(min_blue_n, max_blue_n);
     283    double prob = rnd();
     284    SmartBpGraph bpg;
     285    for (int i=0; i<red_n; ++i) {
     286      bpg.addRedNode();
     287    }
     288    for (int i=0; i<blue_n; ++i) {
     289      bpg.addBlueNode();
     290    }
     291    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     292      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     293        if (rnd() < prob/10.) {
     294          bpg.addEdge(r,b);
     295        }
     296      }
     297    }
     298
     299    HK hk(bpg);
     300    hk.run();
     301    checkMatching(bpg, hk);
     302  }
     303
     304
     305  return 0;
     306}
     307
     308
  • new file lemon/preflow_bp_matching.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1330302343 -3600
    # Node ID 07e1840cc91276dc3312ced8c57dff8eb924567e
    # Parent  fd8741fc21438ea3c829a7d4f6af50d92b7ce426
    Preflow based bipartite matching
    
    diff --git a/lemon/preflow_bp_matching.h b/lemon/preflow_bp_matching.h
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_PREFLOW_BP_MATCHING_H
     20#define LEMON_PREFLOW_BP_MATCHING_H
     21
     22#include <lemon/core.h>
     23#include <lemon/elevator.h>
     24#include <vector>
     25#include <list>
     26#include <queue>
     27
     28/// \ingroup matching
     29/// \file
     30/// \brief A preflow based bipartite matching algorithm.
     31namespace lemon {
     32
     33  /// \brief A preflow based bipartite matching algorithm
     34  ///
     35  /// Finds maximum cardinality matching in a given bipartite
     36  /// graph using specialized preflow algorithm.
     37  template <typename BPG>
     38  class PreflowBpMatching {
     39  public:
     40    /// Type of the bipartite graph
     41    typedef BPG BpGraph;
     42    /// Type of the matching map
     43    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     44    /// Type of the elevator.
     45    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
     46
     47  private:
     48    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
     49
     50  protected:
     51    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
     52    typedef std::list<typename BPG::Arc> CurrentArcSet;
     53    typedef typename
     54    std::list<typename BPG::Arc>::iterator CurrentArcIterator;
     55
     56    const BpGraph& _bpg;
     57    MatchingMap* _matching;
     58    bool _local_matching;
     59    Elevator_t* _elevator;
     60    bool _local_elevator;
     61    typename BPG::template BlueNodeMap<CurrentArcSet> _cas;
     62    typename BPG::template BlueNodeMap<CurrentArcIterator> _ca;
     63    Preflow _flow;
     64
     65    void createStructures() {
     66      if (!_matching) {
     67        _matching = new MatchingMap(_bpg, INVALID);
     68        _local_matching = true;
     69      }
     70      if (!_elevator) {
     71        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
     72                                                  countBlueNodes(_bpg)));
     73        _local_elevator = true;
     74      }
     75    }
     76
     77    void destroyStructures() {
     78      if (_local_matching) {
     79        delete _matching;
     80      }
     81      if (_local_elevator) {
     82        delete _elevator;
     83      }
     84    }
     85
     86    /// \brief Initialize the elevator
     87    ///
     88    /// When this function is called, the elevator is initialized with
     89    /// the exact distance labels corresponding the current matching.
     90    void initElevator() {
     91
     92      // Set the initial flow:
     93      // choose arbitrary outgoing arc for every unmatched red node.
     94      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     95        Arc a = (*_matching)[r] != INVALID ?
     96          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
     97        if (a != INVALID) {
     98          _flow[r] = a;
     99          _cas[_bpg.asBlueNodeUnsafe(_bpg.target(a))].push_front(a);
     100        }
     101      }
     102      // Now we start a backward Bfs from the unmatched blue nodes in the
     103      // residual graph, to supply the exact distance labels to the elevator.
     104      std::vector<BlueNode> act_blue_nodes, next_blue_nodes;
     105      BoolBlueNodeMap reached(_bpg, false);
     106
     107      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     108        if (_cas[b].empty()) {
     109          act_blue_nodes.push_back(b);
     110          reached.set(b, true);
     111        }
     112      }
     113
     114      _elevator->initStart();
     115      while (!act_blue_nodes.empty()) {
     116        for (typename std::vector<BlueNode>::iterator b_it =
     117               act_blue_nodes.begin();
     118             b_it != act_blue_nodes.end(); ++b_it) {
     119
     120          _elevator->initAddItem(*b_it);
     121          for (OutArcIt oi(_bpg, *b_it); oi!=INVALID; ++oi) {
     122            RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(oi));
     123            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
     124            if (b != *b_it && !reached[b]) {
     125              next_blue_nodes.push_back(b);
     126              reached.set(b, true);
     127            }
     128          }
     129        }
     130        _elevator->initNewLevel();
     131        next_blue_nodes.swap(act_blue_nodes);
     132        next_blue_nodes.clear();
     133      }
     134      _elevator->initFinish();
     135
     136      // Activate blue nodes with excess, i. e. those
     137      // which get flow from more than one red node.
     138      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     139        if (_cas[b].size() > 1 && (*_elevator)[b] < _elevator->maxLevel()) {
     140          _elevator->activate(b);
     141        }
     142
     143        _ca[b] = _cas[b].begin();
     144      }
     145    }
     146
     147    PreflowBpMatching() {}
     148
     149  public:
     150
     151    /// \brief Constructor
     152    ///
     153    /// Constructs the class on the given bipartite graph.
     154    PreflowBpMatching(const BpGraph& bpg) : _bpg(bpg),
     155                                            _matching(0),
     156                                            _local_matching(false),
     157                                            _elevator(0),
     158                                            _local_elevator(false),
     159                                            _cas(bpg),
     160                                            _ca(bpg),
     161                                            _flow(bpg, INVALID)
     162    {}
     163
     164    /// \brief Destructor
     165    ///
     166    /// Destructor.
     167    ~PreflowBpMatching() {
     168      destroyStructures();
     169    }
     170
     171    /// \brief Sets the matching map
     172    ///
     173    /// Sets the matching map. It should contain a valid matching,
     174    /// for example the empty matching is fine.
     175    /// If you don't use this function before calling \ref init(),
     176    /// an instance will be allocated automatically.
     177    /// The destructor deallocates this automatically allocated map,
     178    /// of course.
     179    ///
     180    /// \return <tt>(*this)</tt>
     181    PreflowBpMatching& matchingMap(MatchingMap& map) {
     182      _matching = &map;
     183      _local_matching = false;
     184      return *this;
     185    }
     186
     187    /// \brief Returns a const reference to the matching map.
     188    ///
     189    /// Returns a const reference to the matching map, which contains
     190    /// the matching edge for each node (or INVALID for unmatched nodes).
     191    const MatchingMap& matchingMap() const {
     192      return *_matching;
     193    }
     194
     195    /// \brief Sets the elevator used by algorithm.
     196    ///
     197    /// Sets the elevator.
     198    /// If you don't use this function before calling \ref init(),
     199    /// an instance will be allocated automatically.
     200    /// The destructor deallocates this automatically allocated elevator,
     201    /// of course.
     202    /// \return <tt>(*this)</tt>
     203    PreflowBpMatching& elevator(Elevator_t& elevator) {
     204      _elevator = &elevator;
     205      _local_elevator = false;
     206      return *this;
     207    }
     208
     209    /// \brief Returns a const reference to the elevator.
     210    ///
     211    /// Returns a const reference to the elevator, which keeps track
     212    /// of the distance labels of nodes.
     213    const Elevator_t& elevator() const {
     214      return *_elevator;
     215    }
     216
     217    /// \brief Initializes the algorithm
     218    ///
     219    /// Allocates the elevator and the matching map if necessary,
     220    /// and initializes the elevator.
     221    ///
     222    /// \pre \ref init() is not called.
     223    void init() {
     224      createStructures();
     225      initElevator();
     226    }
     227
     228    /// \brief Smarter initialization of the matching.
     229    ///
     230    /// Allocates the matching map if necessary, and initializes
     231    /// the algorithm with a trivially found matching: iterating
     232    /// on every edge, take it in if both endnodes are unmatched.
     233    ///
     234    /// \return Size of the so found matching.
     235    ///
     236    /// \note heuristicInit() replaces init() regarding the preconditions.
     237    int heuristicInit() {
     238      createStructures();
     239      int size = 0;
     240
     241      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     242        bool matched = false;
     243        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
     244          if ((*_matching)[_bpg.target(oi)] == INVALID) {
     245            _matching->set(_bpg.source(oi), oi);
     246            _matching->set(_bpg.target(oi), oi);
     247            matched = true;
     248            ++size;
     249          }
     250        }
     251      }
     252
     253      initElevator();
     254      return size;
     255    }
     256
     257    /// \brief Initialize the matching from a map.
     258    ///
     259    /// Allocates the matching map if necessary, and
     260    /// initializes the algorithm with a matching given as a bool edgemap,
     261    /// in which an edge is true if it is in the matching.
     262    /// Also initializes the elevator.
     263    ///
     264    /// \return \c true if the matching is valid.
     265    ///
     266    /// \note matchingInit() replaces init() regarding the preconditions.
     267    bool matchingInit(const BoolEdgeMap& matching) {
     268      createStructures();
     269
     270      bool valid = true;
     271      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
     272        if (matching[it]) {
     273          Node red = _bpg.redNode(it);
     274          Node blue = _bpg.blueNode(it);
     275          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
     276            valid = false;
     277          } else {
     278            _matching->set(red, it);
     279            _matching->set(blue, it);
     280          }
     281        }
     282      }
     283
     284      initElevator();
     285      return valid;
     286    }
     287
     288
     289    /// \brief Perform a push/relabel operation on a node
     290    ///
     291    /// Perform a push/relabel operation on a node.
     292    ///
     293    /// \pre \ref init() is called, and the node is active.
     294    void push_relabel(const BlueNode& chosen) {
     295      Arc new_flow = INVALID;
     296      while (_ca[chosen] != _cas[chosen].end() && new_flow == INVALID) {
     297        OutArcIt r_oi(_bpg, _bpg.source(*_ca[chosen]));
     298        while (r_oi != INVALID && new_flow == INVALID) {
     299          if ((*_elevator)[_bpg.asBlueNodeUnsafe(_bpg.target(r_oi))] ==
     300              (*_elevator)[chosen] - 1) {
     301            new_flow = r_oi;
     302          }
     303          ++r_oi;
     304        }
     305        if (new_flow == INVALID) ++_ca[chosen];
     306      }
     307
     308      // If an incident blue node is found on one level lower,
     309      // then perform a push, i. e. flip the flow from the
     310      // current arc to 'new_flow'.
     311      if (new_flow != INVALID) {
     312        _ca[chosen] = _cas[chosen].erase(_ca[chosen]);
     313        if (_cas[chosen].size() == 1) {
     314          _elevator->deactivate(chosen);
     315        }
     316        _flow[_bpg.asRedNodeUnsafe(_bpg.source(new_flow))] = new_flow;
     317        BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(new_flow));
     318        _cas[b].push_front(new_flow);
     319        if (_cas[b].size() == 2) {
     320          _elevator->activate(b);
     321        }
     322      }
     323
     324      // If none found, then relabel.
     325      else {
     326        int min = _elevator->maxLevel() + 1;
     327        for (CurrentArcIterator it1 = _cas[chosen].begin();
     328             it1 != _cas[chosen].end();
     329             ++it1) {
     330          for (OutArcIt it2(_bpg, _bpg.asRedNodeUnsafe(_bpg.source(*it1)));
     331               it2 != INVALID;
     332               ++it2) {
     333            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(it2));
     334            if (b != chosen && min > (*_elevator)[b]) {
     335              min = (*_elevator)[b];
     336            }
     337          }
     338        }
     339
     340        int chosen_level = (*_elevator)[chosen];
     341        if (min >= _elevator->maxLevel() - 1) {
     342          _elevator->lift(chosen, _elevator->maxLevel());
     343          _elevator->deactivate(chosen);
     344        } else {
     345          _elevator->lift(chosen, min+1);
     346        }
     347        if (_elevator->emptyLevel(chosen_level)) {
     348          _elevator->liftToTop(chosen_level);
     349        }
     350
     351        // Restore the current arc iterator.
     352        _ca[chosen] = _cas[chosen].begin();
     353      }
     354    }
     355
     356    /// \brief Executes the algorithm
     357    ///
     358    /// Executes push/relabel operation on the highest active node
     359    /// until there is. After start(), writeMatching() should be called to
     360    /// remove unnecessary flow and choose a valid maximum matching.
     361    ///
     362    /// \pre \ref init() must be called before using this function.
     363    void start() {
     364      while (_elevator->highestActive() != INVALID)
     365        push_relabel(_elevator->highestActive());
     366    }
     367
     368    /// \brief Deduces and writes a matching from the preflow
     369    ///
     370    /// The underlying preflow algorithm does not use the matching map,
     371    /// so before using the matching query functions ( \ref mate(), \ref
     372    /// matchingSize(), etc.), you should explicitly call this function
     373    /// to write a matching based on the current preflow to the matching
     374    /// map. When the preflow results an ambiguous matching (a node has
     375    /// incoming flow on more than one arc), it chooses arbitrary one.
     376    void writeMatching() {
     377      for (NodeIt it(_bpg); it != INVALID; ++it) {
     378        _matching->set(it, INVALID);
     379      }
     380      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     381        if (!_cas[b].empty()) {
     382          Arc a = _cas[b].front();
     383          _matching->set(_bpg.source(a), a);
     384          _matching->set(_bpg.target(a), a);
     385        }
     386      }
     387    }
     388
     389    /// \brief Runs the algorithm.
     390    ///
     391    /// bp1.run() is just a shorthand for:
     392    ///
     393    ///\code
     394    /// bp1.heuristicInit();
     395    /// bp1.start();
     396    /// bp1.writeMatching();
     397    ///\endcode
     398    void run() {
     399      heuristicInit();
     400      start();
     401      writeMatching();
     402    }
     403
     404    /// \brief Size of the matching
     405    ///
     406    /// Returns the size of the current matching.
     407    /// Function runs in \f$O(n)\f$ time.
     408    int matchingSize() const {
     409      int size = 0;
     410      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
     411        if ((*_matching)[it] != INVALID) ++size;
     412      }
     413      return size;
     414    }
     415
     416    /// \brief Return \c true if the given edge is in the matching.
     417    ///
     418    /// This function returns \c true if the given edge is in the current
     419    /// matching.
     420    bool matching(const Edge& edge) const {
     421      return edge == (*_matching)[_bpg.redNode(edge)];
     422    }
     423
     424    /// \brief Return the matching edge incident to the given node.
     425    ///
     426    /// This function returns the matching edge incident to the
     427    /// given node in the current matching or \c INVALID if the node is
     428    /// not covered by the matching.
     429    Edge matching(const Node& n) const {
     430      return (*_matching)[n];
     431    }
     432
     433    /// \brief Return the mate of the given node.
     434    ///
     435    /// This function returns the mate of the given node in the current
     436    /// matching or \c INVALID if the node is not covered by the matching.
     437    Node mate(const Node& n) const {
     438      return (*_matching)[n] != INVALID ?
     439        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     440    }
     441
     442    /// \brief Query a vertex cover.
     443    ///
     444    /// This function finds a vertex cover based on the current matching,
     445    /// and in \c cover sets node in the cover \c true and \c false
     446    /// otherwise. In case of maximal matching, this will be a minimal
     447    /// vertex cover.
     448    /// Function runs in \f$O(n + e)\f$ time.
     449    ///
     450    /// \return The number of nodes in the cover.
     451    int cover(BoolNodeMap &cover) const {
     452      int count = 0;
     453      std::queue<RedNode> q;
     454
     455      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
     456        cover[b] = false;
     457
     458      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     459        if ((*_matching)[r] == INVALID) {
     460          cover[r] = false;
     461          q.push(r);
     462        } else {
     463          cover[r] = true;
     464          ++count;
     465        }
     466      }
     467
     468      while (!q.empty()) {
     469        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     470          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     471          if (cover[blue]) continue;
     472
     473          cover[blue] = true;
     474          ++count;
     475
     476          if ((*_matching)[blue] != INVALID) {
     477            RedNode nr = _bpg.redNode((*_matching)[blue]);
     478            q.push(nr);
     479            cover[nr] = false;
     480            --count;
     481          }
     482        }
     483        q.pop();
     484      }
     485
     486      return count;
     487    }
     488
     489    /// \brief Query a barrier of red nodes.
     490    ///
     491    /// This function tries to create a barrier of red nodes, which should:
     492    ///  - contain every unmatched red nodes,
     493    ///  - contain exactly that many matched red nodes as much blue nodes
     494    ///    there are in its neighbor set.
     495    ///
     496    /// Barrier exist if and only if the matching is maximal. If exist, the
     497    /// map is set true for nodes of the barrier, and false for the others.
     498    ///
     499    /// Function runs in \f$O(n + e)\f$ time.
     500    ///
     501    /// \return \c true if a barrier found.
     502    bool redBarrier(BoolRedNodeMap &barrier) const {
     503      bool exist = true;
     504      std::queue<RedNode> q;
     505
     506      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     507        if ((*_matching)[r] == INVALID) {
     508          barrier[r] = true;
     509          q.push(r);
     510        } else {
     511          barrier[r] = false;
     512        }
     513      }
     514
     515      while (!q.empty() && exist) {
     516        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     517          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     518
     519          if (exist &= ((*_matching)[blue] != INVALID)) {
     520            RedNode r = _bpg.redNode((*_matching)[blue]);
     521            if (!barrier[r]) {
     522              q.push(r);
     523              barrier[r] = true;
     524            }
     525          }
     526        }
     527        q.pop();
     528      }
     529
     530      return exist;
     531    }
     532
     533    /// \brief Query a barrier of blue nodes.
     534    ///
     535    /// This function tries to create a barrier of blue nodes, which should:
     536    ///  - contain every unmatched blue nodes,
     537    ///  - contain exactly that many matched blue nodes as much red nodes
     538    ///    there are in its neighbor set.
     539    ///
     540    /// Barrier exist if and only if the matching is maximal. If exist, the
     541    /// map is set true for nodes of the barrier, and false for the others.
     542    ///
     543    /// Function runs in \f$O(n + e)\f$ time.
     544    ///
     545    /// \return \c true if a barrier found.
     546    bool blueBarrier(BoolBlueNodeMap &barrier) const {
     547      bool exist = true;
     548      std::queue<BlueNode> q;
     549
     550      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     551        if ((*_matching)[b] == INVALID) {
     552          barrier[b] = true;
     553          q.push(b);
     554        } else {
     555          barrier[b] = false;
     556        }
     557      }
     558
     559      while (!q.empty() && exist) {
     560        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     561          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
     562
     563          if (exist &= ((*_matching)[r] != INVALID)) {
     564            BlueNode b = _bpg.blueNode((*_matching)[r]);
     565            if (!barrier[b]) {
     566              q.push(b);
     567              barrier[b] = true;
     568            }
     569          }
     570        }
     571        q.pop();
     572      }
     573
     574      return exist;
     575    }
     576  };
     577
     578}
     579#endif
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    4848  path_test
    4949  planarity_test
    5050  preflow_test
     51  preflow_bp_matching_test
    5152  radix_sort_test
    5253  random_test
    5354  suurballe_test
  • test/Makefile.am

    diff --git a/test/Makefile.am b/test/Makefile.am
    a b  
    4646        test/path_test \
    4747        test/planarity_test \
    4848        test/preflow_test \
     49        test/preflow_bp_matching_test \
    4950        test/radix_sort_test \
    5051        test/random_test \
    5152        test/suurballe_test \
     
    104105test_path_test_SOURCES = test/path_test.cc
    105106test_planarity_test_SOURCES = test/planarity_test.cc
    106107test_preflow_test_SOURCES = test/preflow_test.cc
     108test_preflow_bp_matching_test_SOURCES = test/preflow_bp_matching_test.cc
    107109test_radix_sort_test_SOURCES = test/radix_sort_test.cc
    108110test_suurballe_test_SOURCES = test/suurballe_test.cc
    109111test_random_test_SOURCES = test/random_test.cc
  • new file test/preflow_bp_matching_test.cc

    diff --git a/test/preflow_bp_matching_test.cc b/test/preflow_bp_matching_test.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <list>
     22
     23#include <lemon/random.h>
     24#include <lemon/preflow_bp_matching.h>
     25#include <lemon/smart_graph.h>
     26#include <lemon/concepts/bpgraph.h>
     27#include <lemon/concepts/maps.h>
     28#include <lemon/lgf_reader.h>
     29
     30#include "test_tools.h"
     31
     32using namespace std;
     33using namespace lemon;
     34
     35const int lgfn = 3;
     36const string lgf[lgfn] = {
     37  // Empty graph
     38  "@red_nodes\n"
     39  "label\n"
     40  "\n"
     41  "@blue_nodes\n"
     42  "label\n"
     43  "\n"
     44  "@edges\n"
     45  "  label\n"
     46  "\n",
     47
     48  // No red nodes
     49  "@red_nodes\n"
     50  "label\n"
     51  "\n"
     52  "@blue_nodes\n"
     53  "label\n"
     54  "1\n"
     55  "@edges\n"
     56  "  label\n"
     57  "\n",
     58
     59  // No blue nodes
     60  "@red_nodes\n"
     61  "label\n"
     62  "1\n"
     63  "@blue_nodes\n"
     64  "label\n"
     65  "\n"
     66  "@edges\n"
     67  "  label\n"
     68  "\n",
     69};
     70
     71const string u_lgf =
     72  "@red_nodes\n"
     73  "label\n"
     74  "1\n"
     75  "2\n"
     76  "3\n"
     77  "@blue_nodes\n"
     78  "label\n"
     79  "4\n"
     80  "5\n"
     81  "6\n"
     82  "7\n"
     83  "@edges\n"
     84  "label\n"
     85  "1 4 1\n"
     86  "1 5 2\n"
     87  "2 5 3\n"
     88  "2 6 4\n"
     89  "2 7 5\n"
     90  "3 7 6\n"
     91;
     92
     93
     94
     95void checkPreflowBpMatchingCompile() {
     96  typedef concepts::BpGraph BpGraph;
     97  BPGRAPH_TYPEDEFS(BpGraph);
     98  typedef PreflowBpMatching<BpGraph> PBM;
     99
     100  BpGraph graph;
     101  RedNode red;
     102  BlueNode blue;
     103  Node node;
     104  Edge edge;
     105  BoolEdgeMap init_matching(graph);
     106  PBM pbm_test(graph);
     107  const PBM& const_pbm_test = pbm_test;
     108  PBM::MatchingMap matching_map(graph, INVALID);
     109  PBM::Elevator_t elevator(graph);
     110  BpGraph::NodeMap<bool> cov(graph);
     111  int cover_size;
     112  BoolRedNodeMap rb(graph);
     113  BoolBlueNodeMap bb(graph);
     114
     115  pbm_test.matchingMap(matching_map);
     116  pbm_test.elevator(elevator);
     117  pbm_test.init();
     118  pbm_test.heuristicInit();
     119  pbm_test.matchingInit(init_matching);
     120
     121  pbm_test.start();
     122  pbm_test.push_relabel(blue);
     123  pbm_test.run();
     124  pbm_test.writeMatching();
     125
     126  const_pbm_test.matchingSize();
     127  const_pbm_test.matching(node);
     128  const_pbm_test.matching(edge);
     129  const_pbm_test.mate(red);
     130  const_pbm_test.mate(blue);
     131  const_pbm_test.mate(node);
     132  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
     133  edge = max_matching[node];
     134  const PBM::Elevator_t& elev = const_pbm_test.elevator();
     135  blue = elev.highestActive();
     136  cover_size = const_pbm_test.cover(cov);
     137  const_pbm_test.redBarrier(rb);
     138  const_pbm_test.blueBarrier(bb);
     139}
     140
     141typedef SmartBpGraph BpGraph;
     142typedef PreflowBpMatching<BpGraph> PBM;
     143BPGRAPH_TYPEDEFS(BpGraph);
     144
     145void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
     146  BoolBlueNodeMap matched(bpg, false);
     147  for (RedNodeIt it(bpg); it != INVALID; ++it) {
     148    if (pbm.matching(it) != INVALID) {
     149      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
     150      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
     151    }
     152  }
     153  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
     154    // != is used as logical xor here
     155    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
     156  }
     157}
     158
     159void checkCover(const BpGraph& bpg, const PBM& pbm) {
     160  int cov_size;
     161  BoolNodeMap cov(bpg);
     162  cov_size = pbm.cover(cov);
     163  BoolEdgeMap covered(bpg, false);
     164  for (NodeIt n(bpg); n!=INVALID; ++n) {
     165    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
     166      covered[e] = true;
     167    }
     168  }
     169  bool all = true;
     170  for (EdgeIt e(bpg); e!=INVALID; ++e) {
     171    all&=covered[e];
     172  }
     173
     174  check(all, "Invalid cover");
     175  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
     176}
     177
     178void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
     179  BoolRedNodeMap rb(bpg);
     180  BoolBlueNodeMap bb(bpg);
     181
     182  check(pbm.redBarrier(rb), "Suboptimal matching.");
     183  check(pbm.blueBarrier(bb), "Suboptimal matching.");
     184
     185  BoolNodeMap reached(bpg, false);
     186
     187  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     188    if (!rb[r]) continue;
     189    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
     190      reached.set(bpg.target(a), true);
     191    }
     192  }
     193  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     194    if (!bb[b]) continue;
     195    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
     196      reached.set(bpg.target(a), true);
     197    }
     198  }
     199
     200  int ur = 0, ub = 0, rc = 0, bc = 0;
     201  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     202    if (rb[r]) ++rc;
     203    if (pbm.mate(r) == INVALID) ++ur;
     204    if (reached[r]) --bc;
     205  }
     206
     207  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     208    if (bb[b]) ++bc;
     209    if (pbm.mate(b) == INVALID) ++ub;
     210    if (reached[b]) --rc;
     211  }
     212
     213  check(rc == ur, "Red barrier is wrong.");
     214  check(bc == ub, "Blue barrier is wrong.");
     215}
     216
     217void checkMatching(const BpGraph& bpg, const PBM& pbm) {
     218  checkMatchingValidity(bpg, pbm);
     219  checkBarriers(bpg, pbm);
     220  checkCover(bpg, pbm);
     221}
     222
     223
     224int main() {
     225  // Checking hard wired extremal graphs
     226  for (int i=0; i<lgfn; ++i) {
     227    BpGraph bpg;
     228    istringstream lgfs(lgf[i]);
     229    bpGraphReader(bpg, lgfs).run();
     230    PBM pbm(bpg);
     231    pbm.run();
     232    checkMatching(bpg, pbm);
     233  }
     234
     235  // Checking some use cases
     236  BpGraph bpg;
     237  istringstream u_lgfs(u_lgf);
     238  bpGraphReader(bpg, u_lgfs).run();
     239
     240  {
     241    PBM pbm(bpg);
     242    pbm.init();
     243    BpGraph::BlueNode b = pbm.elevator().highestActive();
     244    if (b != INVALID) {
     245      pbm.push_relabel(b);
     246    }
     247    pbm.start();
     248    pbm.writeMatching();
     249    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
     250    checkMatching(bpg, pbm);
     251  }
     252
     253  {
     254    PBM pbm(bpg);
     255    pbm.heuristicInit();
     256    BpGraph::BlueNode b = pbm.elevator().highestActive();
     257    if (b != INVALID) {
     258      pbm.push_relabel(b);
     259    }
     260    pbm.start();
     261    checkMatching(bpg, pbm);
     262  }
     263  {
     264    PBM pbm(bpg);
     265    PBM::MatchingMap m(bpg, INVALID);
     266    pbm.matchingMap(m);
     267    pbm.run();
     268    checkMatching(bpg, pbm);
     269  }
     270  {
     271    PBM pbm(bpg);
     272    BpGraph::EdgeMap<bool> init_matching(bpg, false);
     273    BpGraph::Edge e;
     274    bpg.first(e);
     275    init_matching[e] = true;
     276    check(pbm.matchingInit(init_matching),
     277          "Wrong result from matchingInit()");
     278    pbm.start();
     279    pbm.writeMatching();
     280    checkMatching(bpg, pbm);
     281  }
     282  {
     283    PBM pbm(bpg);
     284    PBM::MatchingMap m(bpg, INVALID);
     285    pbm.matchingMap(m);
     286    BpGraph::EdgeMap<bool> init_matching(bpg, true);
     287    check(!pbm.matchingInit(init_matching),
     288          "Wrong result from matchingInit()");
     289    pbm.start();
     290    pbm.writeMatching();
     291    checkMatching(bpg, pbm);
     292  }
     293
     294  // Checking some random generated graphs
     295  const int random_test_n = 10;
     296  const int max_red_n = 1000;
     297  const int min_red_n = 10;
     298  const int max_blue_n = 1000;
     299  const int min_blue_n = 10;
     300  for (int i=0; i<random_test_n; ++i) {
     301    int red_n = rnd.integer(min_red_n, max_red_n);
     302    int blue_n = rnd.integer(min_blue_n, max_blue_n);
     303    double prob = rnd();
     304
     305    SmartBpGraph bpg;
     306    for (int j=0; j<red_n; ++j) {
     307      bpg.addRedNode();
     308    }
     309    for (int j=0; j<blue_n; ++j) {
     310      bpg.addBlueNode();
     311    }
     312    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     313      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     314        if (rnd() < prob/10.) {
     315          bpg.addEdge(r,b);
     316        }
     317      }
     318    }
     319
     320    PBM pbm(bpg);
     321    pbm.run();
     322    checkMatching(bpg, pbm);
     323
     324    // do not always use the highest active node
     325    PBM pbm2(bpg);
     326    PBM::Elevator_t elev(bpg);
     327    pbm2.elevator(elev);
     328    pbm2.init();
     329    while (elev.highestActive() != INVALID) {
     330      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     331        if (elev.active(b)) pbm2.push_relabel(b);
     332      }
     333    }
     334    pbm2.writeMatching();
     335    checkMatching(bpg, pbm2);
     336  }
     337
     338
     339  return 0;
     340}
     341
     342
  • new file lemon/preflow_bp_matching_2.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1330302480 -3600
    # Node ID 7ae506e088394579bddd0d4c28e76a19459b04ee
    # Parent  07e1840cc91276dc3312ced8c57dff8eb924567e
    Other variant of preflow based bipartite matching
    
    diff --git a/lemon/preflow_bp_matching_2.h b/lemon/preflow_bp_matching_2.h
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_PREFLOW_BP_MATCHING_2_H
     20#define LEMON_PREFLOW_BP_MATCHING_2_H
     21
     22#include <lemon/core.h>
     23#include <lemon/elevator.h>
     24#include <vector>
     25#include <queue>
     26
     27/// \ingroup matching
     28/// \file
     29/// \brief A preflow based bipartite matching algorithm.
     30namespace lemon {
     31
     32  /// \brief A preflow based bipartite matching algorithm
     33  ///
     34  /// Finds maximum cardinality matching in a given bipartite
     35  /// graph using specialized preflow algorithm.
     36  template <typename BPG>
     37  class PreflowBpMatching2 {
     38  public:
     39    /// Type of the bipartite graph
     40    typedef BPG BpGraph;
     41    /// Type of the matching map
     42    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     43    /// Type of the elevator.
     44    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
     45
     46  private:
     47    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
     48
     49  protected:
     50    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
     51    typedef typename BPG::template BlueNodeMap<int> Excess;
     52    typedef std::vector<typename BPG::BlueNode> BlueNodes;
     53    typedef std::vector<typename BPG::Arc> BlueFlow;
     54    typedef typename BPG::template BlueNodeMap<BlueFlow> IncomingFlows;
     55
     56    const BpGraph& _bpg;
     57    MatchingMap* _matching;
     58    bool _local_matching;
     59    Elevator_t* _elevator;
     60    bool _local_elevator;
     61    Preflow _flow;
     62    Excess _excess;
     63
     64    void createStructures() {
     65      if (!_matching) {
     66        _matching = new MatchingMap(_bpg, INVALID);
     67        _local_matching = true;
     68      }
     69      if (!_elevator) {
     70        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
     71                                                  countBlueNodes(_bpg)));
     72        _local_elevator = true;
     73      }
     74    }
     75
     76    void destroyStructures() {
     77      if (_local_matching) {
     78        delete _matching;
     79      }
     80      if (_local_elevator) {
     81        delete _elevator;
     82      }
     83    }
     84
     85    /// \brief Initialize the elevator
     86    ///
     87    /// When this function is called, the elevator is initialized with
     88    /// the exact distance labels corresponding the current matching.
     89    void initElevator() {
     90      BlueNodes act_blue_nodes, next_blue_nodes;
     91      BoolBlueNodeMap reached(_bpg, false);
     92      IncomingFlows incoming(_bpg);
     93
     94      // Set the initial flow:
     95      // choose arbitrary outgoing arc for every unmatched red node.
     96      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     97        Arc a = (*_matching)[r] != INVALID ?
     98          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
     99        if (a != INVALID) {
     100          _flow[r] = a;
     101          BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(a));
     102          incoming[b].push_back(a);
     103          if (++_excess[b] == 2) {
     104            act_blue_nodes.push_back(b);
     105            reached.set(b, true);
     106          }
     107        }
     108      }
     109      // Now we start a Bfs from the blue nodes with excess
     110      // (those which get flow from more than one red node) in the
     111      // residual graph, to supply the exact distance labels to the elevator.
     112      _elevator->initStart();
     113      while (!act_blue_nodes.empty()) {
     114        for (typename BlueNodes::iterator b_it = act_blue_nodes.begin();
     115             b_it != act_blue_nodes.end(); ++b_it) {
     116
     117          _elevator->initAddItem(*b_it);
     118          for (typename BlueFlow::iterator b_in = incoming[*b_it].begin();
     119               b_in != incoming[*b_it].end(); ++b_in) {
     120            RedNode r = _bpg.asRedNodeUnsafe(_bpg.source(*b_in));
     121            for (OutArcIt r_oi(_bpg, r); r_oi!=INVALID; ++r_oi) {
     122              BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(r_oi));
     123              if (!reached[b]) {
     124                next_blue_nodes.push_back(b);
     125                reached.set(b, true);
     126              }
     127            }
     128          }
     129        }
     130        _elevator->initNewLevel();
     131        next_blue_nodes.swap(act_blue_nodes);
     132        next_blue_nodes.clear();
     133      }
     134      _elevator->initFinish();
     135
     136      // Activate blue nodes with 0 excess, i. e. those
     137      // which does not get any flow.
     138      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     139        if (_excess[b] == 0  && (*_elevator)[b] < _elevator->maxLevel()) {
     140          _elevator->activate(b);
     141        }
     142      }
     143    }
     144
     145    PreflowBpMatching2() {}
     146
     147  public:
     148
     149    /// \brief Constructor
     150    ///
     151    /// Constructs the class on the given bipartite graph.
     152    PreflowBpMatching2(const BpGraph& bpg) : _bpg(bpg),
     153                                             _matching(0),
     154                                             _local_matching(false),
     155                                             _elevator(0),
     156                                             _local_elevator(false),
     157                                             _flow(bpg, INVALID),
     158                                             _excess(_bpg, 0)
     159
     160    {}
     161
     162    /// \brief Destructor
     163    ///
     164    /// Destructor.
     165    ~PreflowBpMatching2() {
     166      destroyStructures();
     167    }
     168
     169    /// \brief Sets the matching map
     170    ///
     171    /// Sets the matching map. It should contain a valid matching,
     172    /// for example the empty matching is fine.
     173    /// If you don't use this function before calling \ref init(),
     174    /// an instance will be allocated automatically.
     175    /// The destructor deallocates this automatically allocated map,
     176    /// of course.
     177    ///
     178    /// \return <tt>(*this)</tt>
     179    PreflowBpMatching2& matchingMap(MatchingMap& map) {
     180      _matching = &map;
     181      _local_matching = false;
     182      return *this;
     183    }
     184
     185    /// \brief Returns a const reference to the matching map.
     186    ///
     187    /// Returns a const reference to the matching map, which contains
     188    /// the matching edge for each node (or INVALID for unmatched nodes).
     189    const MatchingMap& matchingMap() const {
     190      return *_matching;
     191    }
     192
     193    /// \brief Sets the elevator used by algorithm.
     194    ///
     195    /// Sets the elevator.
     196    /// If you don't use this function before calling \ref init(),
     197    /// an instance will be allocated automatically.
     198    /// The destructor deallocates this automatically allocated elevator,
     199    /// of course.
     200    /// \return <tt>(*this)</tt>
     201    PreflowBpMatching2& elevator(Elevator_t& elevator) {
     202      _elevator = &elevator;
     203      _local_elevator = false;
     204      return *this;
     205    }
     206
     207    /// \brief Returns a const reference to the elevator.
     208    ///
     209    /// Returns a const reference to the elevator, which keeps track
     210    /// of the distance labels of nodes.
     211    const Elevator_t& elevator() const {
     212      return *_elevator;
     213    }
     214
     215    /// \brief Initializes the algorithm
     216    ///
     217    /// Allocates the elevator and the matching map if necessary,
     218    /// and initializes the elevator.
     219    ///
     220    /// \pre \ref init() is not called.
     221    void init() {
     222      createStructures();
     223      initElevator();
     224    }
     225
     226    /// \brief Smarter initialization of the matching.
     227    ///
     228    /// Allocates the matching map if necessary, and initializes
     229    /// the algorithm with a trivially found matching: iterating
     230    /// on every edge, take it in if both endnodes are unmatched.
     231    ///
     232    /// \return Size of the so found matching.
     233    ///
     234    /// \note heuristicInit() replaces init() regarding the preconditions.
     235    int heuristicInit() {
     236      createStructures();
     237      int size = 0;
     238
     239      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     240        bool matched = false;
     241        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
     242          if ((*_matching)[_bpg.target(oi)] == INVALID) {
     243            _matching->set(_bpg.source(oi), oi);
     244            _matching->set(_bpg.target(oi), oi);
     245            matched = true;
     246            ++size;
     247          }
     248        }
     249      }
     250
     251      initElevator();
     252      return size;
     253    }
     254
     255    /// \brief Initialize the matching from a map.
     256    ///
     257    /// Allocates the matching map if necessary, and
     258    /// initializes the algorithm with a matching given as a bool edgemap,
     259    /// in which an edge is true if it is in the matching.
     260    /// Also initializes the elevator.
     261    ///
     262    /// \return \c true if the matching is valid.
     263    ///
     264    /// \note matchingInit() replaces init() regarding the preconditions.
     265    bool matchingInit(const BoolEdgeMap& matching) {
     266      createStructures();
     267
     268      bool valid = true;
     269      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
     270        if (matching[it]) {
     271          Node red = _bpg.redNode(it);
     272          Node blue = _bpg.blueNode(it);
     273          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
     274            valid = false;
     275          } else {
     276            _matching->set(red, it);
     277            _matching->set(blue, it);
     278          }
     279        }
     280      }
     281
     282      initElevator();
     283      return valid;
     284    }
     285
     286
     287    /// \brief Perform a pull/relabel operation on a node
     288    ///
     289    /// Perform a pull/relabel operation on a node:
     290    /// search an incident blue node on one level lower,
     291    /// and flip its flow to the 'chosen', or relabel
     292    /// if none found.
     293    ///
     294    /// \pre \ref init() is called, and the node is active.
     295    void pull_relabel(const BlueNode& chosen) {
     296
     297      int min_inc_level = _elevator->maxLevel() + 1;
     298      for (OutArcIt b_oi(_bpg, chosen); b_oi!=INVALID; ++b_oi) {
     299        RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(b_oi));
     300        BlueNode lower = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
     301        if ((*_elevator)[(lower)] < (*_elevator)[chosen]) {
     302          if (--_excess[lower] == 0) {
     303            _elevator->activate(lower);
     304          }
     305          _flow[r] = _bpg.direct(b_oi, true);
     306          ++_excess[chosen];
     307          _elevator->deactivate(chosen);
     308          break;
     309        } else {
     310          if ((*_elevator)[lower] < min_inc_level &&
     311              lower != chosen) {
     312            min_inc_level = (*_elevator)[lower];
     313          }
     314        }
     315      }
     316      // If none found, then relabel.
     317      if (_excess[chosen] == 0) {
     318        int chosen_level = (*_elevator)[chosen];
     319        if (min_inc_level >= _elevator->maxLevel() - 1) {
     320          _elevator->lift(chosen, _elevator->maxLevel());
     321          _elevator->deactivate(chosen);
     322        } else {
     323          _elevator->lift(chosen, min_inc_level+1);
     324        }
     325        if (_elevator->emptyLevel(chosen_level)) {
     326          _elevator->liftToTop(chosen_level);
     327        }
     328      }
     329    }
     330
     331    /// \brief Executes the algorithm
     332    ///
     333    /// Executes pull/relabel operation on the highest active node
     334    /// until there is. After start(), writeMatching() should be called to
     335    /// remove unnecessary flow and choose a valid maximum matching.
     336    ///
     337    /// \pre \ref init() must be called before using this function.
     338    void start() {
     339      while (_elevator->highestActive() != INVALID)
     340        pull_relabel(_elevator->highestActive());
     341    }
     342
     343    /// \brief Deduces and writes a matching from the preflow
     344    ///
     345    /// The underlying preflow algorithm does not use the matching map,
     346    /// so before using the matching query functions ( \ref mate(), \ref
     347    /// matchingSize(), etc.), you should explicitly call this function
     348    /// to write a matching based on the current preflow to the matching
     349    /// map. When the preflow results an ambiguous matching (a node has
     350    /// incoming flow on more than one arc), it chooses arbitrary one.
     351    void writeMatching() {
     352      for (NodeIt it(_bpg); it != INVALID; ++it) {
     353        _matching->set(it, INVALID);
     354      }
     355      for (RedNodeIt r(_bpg); r != INVALID; ++r) {
     356        if (_flow[r] != INVALID) {
     357          BlueNode mate = _bpg.blueNode(_flow[r]);
     358          if (_excess[mate] > 1) {
     359            --_excess[mate];
     360          } else {
     361            _matching->set(mate, _flow[r]);
     362            _matching->set(r, _flow[r]);
     363          }
     364        }
     365      }
     366    }
     367
     368    /// \brief Runs the algorithm.
     369    ///
     370    /// bp1.run() is just a shorthand for:
     371    ///
     372    ///\code
     373    /// bp1.heuristicInit();
     374    /// bp1.start();
     375    /// bp1.writeMatching();
     376    ///\endcode
     377    void run() {
     378      heuristicInit();
     379      start();
     380      writeMatching();
     381    }
     382
     383    /// \brief Size of the matching
     384    ///
     385    /// Returns the size of the current matching.
     386    /// Function runs in \f$O(n)\f$ time.
     387    int matchingSize() const {
     388      int size = 0;
     389      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
     390        if ((*_matching)[it] != INVALID) ++size;
     391      }
     392      return size;
     393    }
     394
     395    /// \brief Return \c true if the given edge is in the matching.
     396    ///
     397    /// This function returns \c true if the given edge is in the current
     398    /// matching.
     399    bool matching(const Edge& edge) const {
     400      return edge == (*_matching)[_bpg.redNode(edge)];
     401    }
     402
     403    /// \brief Return the matching edge incident to the given node.
     404    ///
     405    /// This function returns the matching edge incident to the
     406    /// given node in the current matching or \c INVALID if the node is
     407    /// not covered by the matching.
     408    Edge matching(const Node& n) const {
     409      return (*_matching)[n];
     410    }
     411
     412    /// \brief Return the mate of the given node.
     413    ///
     414    /// This function returns the mate of the given node in the current
     415    /// matching or \c INVALID if the node is not covered by the matching.
     416    Node mate(const Node& n) const {
     417      return (*_matching)[n] != INVALID ?
     418        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     419    }
     420
     421    /// \brief Query a vertex cover.
     422    ///
     423    /// This function finds a vertex cover based on the current matching,
     424    /// and in \c cover sets node in the cover \c true and \c false
     425    /// otherwise. In case of maximal matching, this will be a minimal
     426    /// vertex cover.
     427    /// Function runs in \f$O(n + e)\f$ time.
     428    ///
     429    /// \return The number of nodes in the cover.
     430    int cover(BoolNodeMap &cover) const {
     431      int count = 0;
     432      std::queue<RedNode> q;
     433
     434      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
     435        cover[b] = false;
     436
     437      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     438        if ((*_matching)[r] == INVALID) {
     439          cover[r] = false;
     440          q.push(r);
     441        } else {
     442          cover[r] = true;
     443          ++count;
     444        }
     445      }
     446
     447      while (!q.empty()) {
     448        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     449          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     450          if (cover[blue]) continue;
     451
     452          cover[blue] = true;
     453          ++count;
     454
     455          if ((*_matching)[blue] != INVALID) {
     456            RedNode nr = _bpg.redNode((*_matching)[blue]);
     457            q.push(nr);
     458            cover[nr] = false;
     459            --count;
     460          }
     461        }
     462        q.pop();
     463      }
     464
     465      return count;
     466    }
     467
     468    /// \brief Query a barrier of red nodes.
     469    ///
     470    /// This function tries to create a barrier of red nodes, which should:
     471    ///  - contain every unmatched red nodes,
     472    ///  - contain exactly that many matched red nodes as much blue nodes
     473    ///    there are in its neighbor set.
     474    ///
     475    /// Barrier exist if and only if the matching is maximal. If exist, the
     476    /// map is set true for nodes of the barrier, and false for the others.
     477    ///
     478    /// Function runs in \f$O(n + e)\f$ time.
     479    ///
     480    /// \return \c true if a barrier found.
     481    bool redBarrier(BoolRedNodeMap &barrier) const {
     482      bool exist = true;
     483      std::queue<RedNode> q;
     484
     485      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     486        if ((*_matching)[r] == INVALID) {
     487          barrier[r] = true;
     488          q.push(r);
     489        } else {
     490          barrier[r] = false;
     491        }
     492      }
     493
     494      while (!q.empty() && exist) {
     495        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     496          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     497
     498          if (exist &= ((*_matching)[blue] != INVALID)) {
     499            RedNode r = _bpg.redNode((*_matching)[blue]);
     500            if (!barrier[r]) {
     501              q.push(r);
     502              barrier[r] = true;
     503            }
     504          }
     505        }
     506        q.pop();
     507      }
     508
     509      return exist;
     510    }
     511
     512    /// \brief Query a barrier of blue nodes.
     513    ///
     514    /// This function tries to create a barrier of blue nodes, which should:
     515    ///  - contain every unmatched blue nodes,
     516    ///  - contain exactly that many matched blue nodes as much red nodes
     517    ///    there are in its neighbor set.
     518    ///
     519    /// Barrier exist if and only if the matching is maximal. If exist, the
     520    /// map is set true for nodes of the barrier, and false for the others.
     521    ///
     522    /// Function runs in \f$O(n + e)\f$ time.
     523    ///
     524    /// \return \c true if a barrier found.
     525    bool blueBarrier(BoolBlueNodeMap &barrier) const {
     526      bool exist = true;
     527      std::queue<BlueNode> q;
     528
     529      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     530        if ((*_matching)[b] == INVALID) {
     531          barrier[b] = true;
     532          q.push(b);
     533        } else {
     534          barrier[b] = false;
     535        }
     536      }
     537
     538      while (!q.empty() && exist) {
     539        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     540          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
     541
     542          if (exist &= ((*_matching)[r] != INVALID)) {
     543            BlueNode b = _bpg.blueNode((*_matching)[r]);
     544            if (!barrier[b]) {
     545              q.push(b);
     546              barrier[b] = true;
     547            }
     548          }
     549        }
     550        q.pop();
     551      }
     552
     553      return exist;
     554    }
     555  };
     556
     557}
     558#endif
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    4848  path_test
    4949  planarity_test
    5050  preflow_test
     51  preflow_bp_matching_2_test
    5152  preflow_bp_matching_test
    5253  radix_sort_test
    5354  random_test
  • test/Makefile.am

    diff --git a/test/Makefile.am b/test/Makefile.am
    a b  
    4646        test/path_test \
    4747        test/planarity_test \
    4848        test/preflow_test \
     49        test/preflow_bp_matching_2_test \
    4950        test/preflow_bp_matching_test \
    5051        test/radix_sort_test \
    5152        test/random_test \
     
    105106test_path_test_SOURCES = test/path_test.cc
    106107test_planarity_test_SOURCES = test/planarity_test.cc
    107108test_preflow_test_SOURCES = test/preflow_test.cc
     109test_preflow_bp_matching_2_test_SOURCES = test/preflow_bp_matching_2_test.cc
    108110test_preflow_bp_matching_test_SOURCES = test/preflow_bp_matching_test.cc
    109111test_radix_sort_test_SOURCES = test/radix_sort_test.cc
    110112test_suurballe_test_SOURCES = test/suurballe_test.cc
  • new file test/preflow_bp_matching_2_test.cc

    diff --git a/test/preflow_bp_matching_2_test.cc b/test/preflow_bp_matching_2_test.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <list>
     22
     23#include <lemon/random.h>
     24#include <lemon/preflow_bp_matching_2.h>
     25#include <lemon/smart_graph.h>
     26#include <lemon/concepts/bpgraph.h>
     27#include <lemon/concepts/maps.h>
     28#include <lemon/lgf_reader.h>
     29
     30#include "test_tools.h"
     31
     32using namespace std;
     33using namespace lemon;
     34
     35const int lgfn = 3;
     36const string lgf[lgfn] = {
     37  // Empty graph
     38  "@red_nodes\n"
     39  "label\n"
     40  "\n"
     41  "@blue_nodes\n"
     42  "label\n"
     43  "\n"
     44  "@edges\n"
     45  "  label\n"
     46  "\n",
     47
     48  // No red nodes
     49  "@red_nodes\n"
     50  "label\n"
     51  "\n"
     52  "@blue_nodes\n"
     53  "label\n"
     54  "1\n"
     55  "@edges\n"
     56  "  label\n"
     57  "\n",
     58
     59  // No blue nodes
     60  "@red_nodes\n"
     61  "label\n"
     62  "1\n"
     63  "@blue_nodes\n"
     64  "label\n"
     65  "\n"
     66  "@edges\n"
     67  "  label\n"
     68  "\n",
     69};
     70
     71const string u_lgf =
     72  "@red_nodes\n"
     73  "label\n"
     74  "1\n"
     75  "2\n"
     76  "3\n"
     77  "@blue_nodes\n"
     78  "label\n"
     79  "4\n"
     80  "5\n"
     81  "6\n"
     82  "7\n"
     83  "@edges\n"
     84  "label\n"
     85  "1 4 1\n"
     86  "1 5 2\n"
     87  "2 5 3\n"
     88  "2 6 4\n"
     89  "2 7 5\n"
     90  "3 7 6\n"
     91;
     92
     93
     94
     95void checkPreflowBpMatching2Compile() {
     96  typedef concepts::BpGraph BpGraph;
     97  BPGRAPH_TYPEDEFS(BpGraph);
     98  typedef PreflowBpMatching2<BpGraph> PBM;
     99
     100  BpGraph graph;
     101  RedNode red;
     102  BlueNode blue;
     103  Node node;
     104  Edge edge;
     105  BoolEdgeMap init_matching(graph);
     106  PBM pbm_test(graph);
     107  const PBM& const_pbm_test = pbm_test;
     108  PBM::MatchingMap matching_map(graph, INVALID);
     109  PBM::Elevator_t elevator(graph);
     110  BpGraph::NodeMap<bool> cov(graph);
     111  int cover_size;
     112  BoolRedNodeMap rb(graph);
     113  BoolBlueNodeMap bb(graph);
     114
     115  pbm_test.matchingMap(matching_map);
     116  pbm_test.elevator(elevator);
     117  pbm_test.init();
     118  pbm_test.heuristicInit();
     119  pbm_test.matchingInit(init_matching);
     120
     121  pbm_test.start();
     122  pbm_test.pull_relabel(blue);
     123  pbm_test.run();
     124  pbm_test.writeMatching();
     125
     126  const_pbm_test.matchingSize();
     127  const_pbm_test.matching(node);
     128  const_pbm_test.matching(edge);
     129  const_pbm_test.mate(red);
     130  const_pbm_test.mate(blue);
     131  const_pbm_test.mate(node);
     132  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
     133  edge = max_matching[node];
     134  const PBM::Elevator_t& elev = const_pbm_test.elevator();
     135  blue = elev.highestActive();
     136  cover_size = const_pbm_test.cover(cov);
     137  const_pbm_test.redBarrier(rb);
     138  const_pbm_test.blueBarrier(bb);
     139}
     140
     141typedef SmartBpGraph BpGraph;
     142typedef PreflowBpMatching2<BpGraph> PBM;
     143BPGRAPH_TYPEDEFS(BpGraph);
     144
     145void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
     146  BoolBlueNodeMap matched(bpg, false);
     147  for (RedNodeIt it(bpg); it != INVALID; ++it) {
     148    if (pbm.matching(it) != INVALID) {
     149      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
     150      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
     151    }
     152  }
     153  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
     154    // != is used as logical xor here
     155    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
     156  }
     157}
     158
     159void checkCover(const BpGraph& bpg, const PBM& pbm) {
     160  int cov_size;
     161  BoolNodeMap cov(bpg);
     162  cov_size = pbm.cover(cov);
     163  BoolEdgeMap covered(bpg, false);
     164  for (NodeIt n(bpg); n!=INVALID; ++n) {
     165    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
     166      covered[e] = true;
     167    }
     168  }
     169  bool all = true;
     170  for (EdgeIt e(bpg); e!=INVALID; ++e) {
     171    all&=covered[e];
     172  }
     173
     174  check(all, "Invalid cover");
     175  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
     176}
     177void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
     178  BoolRedNodeMap rb(bpg);
     179  BoolBlueNodeMap bb(bpg);
     180
     181  check(pbm.redBarrier(rb), "Suboptimal matching.");
     182  check(pbm.blueBarrier(bb), "Suboptimal matching.");
     183
     184  BoolNodeMap reached(bpg, false);
     185
     186  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     187    if (!rb[r]) continue;
     188    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
     189      reached.set(bpg.target(a), true);
     190    }
     191  }
     192  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     193    if (!bb[b]) continue;
     194    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
     195      reached.set(bpg.target(a), true);
     196    }
     197  }
     198
     199  int ur = 0, ub = 0, rc = 0, bc = 0;
     200  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     201    if (rb[r]) ++rc;
     202    if (pbm.mate(r) == INVALID) ++ur;
     203    if (reached[r]) --bc;
     204  }
     205
     206  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     207    if (bb[b]) ++bc;
     208    if (pbm.mate(b) == INVALID) ++ub;
     209    if (reached[b]) --rc;
     210  }
     211
     212  check(rc == ur, "Red barrier is wrong.");
     213  check(bc == ub, "Blue barrier is wrong.");
     214}
     215
     216void checkMatching(const BpGraph& bpg, const PBM& pbm) {
     217  checkMatchingValidity(bpg, pbm);
     218  checkBarriers(bpg, pbm);
     219  checkCover(bpg, pbm);
     220}
     221
     222
     223int main() {
     224  // Checking hard wired extremal graphs
     225  for (int i=0; i<lgfn; ++i) {
     226    BpGraph bpg;
     227    istringstream lgfs(lgf[i]);
     228    bpGraphReader(bpg, lgfs).run();
     229    PBM pbm(bpg);
     230    pbm.run();
     231    checkMatching(bpg, pbm);
     232  }
     233
     234  // Checking some use cases
     235  BpGraph bpg;
     236  istringstream u_lgfs(u_lgf);
     237  bpGraphReader(bpg, u_lgfs).run();
     238
     239  {
     240    PBM pbm(bpg);
     241    pbm.init();
     242    BpGraph::BlueNode b = pbm.elevator().highestActive();
     243    if (b != INVALID) {
     244      pbm.pull_relabel(b);
     245    }
     246    pbm.start();
     247    pbm.writeMatching();
     248    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
     249    checkMatching(bpg, pbm);
     250  }
     251
     252  {
     253    PBM pbm(bpg);
     254    pbm.heuristicInit();
     255    BpGraph::BlueNode b = pbm.elevator().highestActive();
     256    if (b != INVALID) {
     257      pbm.pull_relabel(b);
     258    }
     259    pbm.start();
     260    checkMatching(bpg, pbm);
     261  }
     262  {
     263    PBM pbm(bpg);
     264    PBM::MatchingMap m(bpg, INVALID);
     265    pbm.matchingMap(m);
     266    pbm.run();
     267    checkMatching(bpg, pbm);
     268
     269  }
     270  {
     271    PBM pbm(bpg);
     272    BpGraph::EdgeMap<bool> init_matching(bpg, false);
     273    BpGraph::Edge e;
     274    bpg.first(e);
     275    init_matching[e] = true;
     276    check(pbm.matchingInit(init_matching),
     277          "Wrong result from matchingInit()");
     278    pbm.start();
     279    pbm.writeMatching();
     280    checkMatching(bpg, pbm);
     281  }
     282  {
     283    PBM pbm(bpg);
     284    PBM::MatchingMap m(bpg, INVALID);
     285    pbm.matchingMap(m);
     286    BpGraph::EdgeMap<bool> init_matching(bpg, true);
     287    check(!pbm.matchingInit(init_matching),
     288          "Wrong result from matchingInit()");
     289    pbm.start();
     290    pbm.writeMatching();
     291    checkMatching(bpg, pbm);
     292  }
     293
     294  // Checking some random generated graphs
     295  const int random_test_n = 10;
     296  const int max_red_n = 1000;
     297  const int min_red_n = 10;
     298  const int max_blue_n = 1000;
     299  const int min_blue_n = 10;
     300  for (int i=0; i<random_test_n; ++i) {
     301    int red_n = rnd.integer(min_red_n, max_red_n);
     302    int blue_n = rnd.integer(min_blue_n, max_blue_n);
     303    double prob = rnd();
     304
     305    SmartBpGraph bpg;
     306    for (int j=0; j<red_n; ++j) {
     307      bpg.addRedNode();
     308    }
     309    for (int j=0; j<blue_n; ++j) {
     310      bpg.addBlueNode();
     311    }
     312    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     313      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     314        if (rnd() < prob/10.) {
     315          bpg.addEdge(r,b);
     316        }
     317      }
     318    }
     319
     320    PBM pbm(bpg);
     321    pbm.run();
     322    checkMatching(bpg, pbm);
     323
     324    // do not always use the highest active node
     325    PBM pbm2(bpg);
     326    PBM::Elevator_t elev(bpg);
     327    pbm2.elevator(elev);
     328    pbm2.init();
     329    while (elev.highestActive() != INVALID) {
     330      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     331        if (elev.active(b)) pbm2.pull_relabel(b);
     332      }
     333    }
     334    pbm2.writeMatching();
     335    checkMatching(bpg, pbm2);
     336  }
     337
     338
     339  return 0;
     340}
     341
     342
  • lemon/dimacs.h

    # HG changeset patch
    # User Daniel Poroszkai <poroszd@inf.elte.hu>
    # Date 1330303848 -3600
    # Node ID 06e679b37a9042e3f77b40b170904e3b59ef213d
    # Parent  7ae506e088394579bddd0d4c28e76a19459b04ee
    Tools for bipartite matching benchmarks
    
    diff --git a/lemon/dimacs.h b/lemon/dimacs.h
    a b  
    22 *
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2010
     5 * Copyright (C) 2003-2012
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
    88 *
     
    4545      MIN,   ///< DIMACS file type for minimum cost flow problems.
    4646      MAX,   ///< DIMACS file type for maximum flow problems.
    4747      SP,    ///< DIMACS file type for shostest path problems.
    48       MAT    ///< DIMACS file type for plain graphs and matching problems.
     48      MAT,   ///< DIMACS file type for plain graphs and matching problems.
     49      ASN,   ///< DIMACS file type for assignment problems.
     50      UBM,   ///< DIMACS file type for unweighted bipartite matching problems.
     51      WBM    ///< DIMACS file type for weighted bipartite matching problems.
    4952    };
    5053    ///The file type
    5154    Type type;
     
    8285              else if(problem=="max") r.type=DimacsDescriptor::MAX;
    8386              else if(problem=="sp") r.type=DimacsDescriptor::SP;
    8487              else if(problem=="mat") r.type=DimacsDescriptor::MAT;
     88              else if(problem=="asn") r.type=DimacsDescriptor::ASN;
     89              else if(problem=="ubm") r.type=DimacsDescriptor::UBM;
     90              else if(problem=="wbm") r.type=DimacsDescriptor::WBM;
    8591              else throw FormatError("Unknown problem type");
    8692              return r;
    8793            }
     
    101107  }
    102108
    103109
     110
     111  /// \brief DIMACS unweighted bipartite matching problem reader function.
     112  ///
     113  /// This function reads an unweighted bipartite matching problem instance
     114  /// from DIMACS format, i.e. from a DIMACS file having a line starting with
     115  /// \code
     116  ///   p ubm
     117  /// \endcode
     118  /// and creates the corresponding bipartite graph.
     119  ///
     120  /// At the beginning, \c g is cleared by \c g.clear().
     121  ///
     122  /// If the file type was previously evaluated by dimacsType(), then
     123  /// the descriptor struct should be given by the \c dest parameter.
     124  template <typename BpGraph>
     125  void readDimacsUbm(std::istream& is,
     126                     BpGraph &g,
     127                     DimacsDescriptor desc=DimacsDescriptor())
     128  {
     129    g.clear();
     130
     131    if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
     132    if(desc.type!=DimacsDescriptor::UBM)
     133      throw FormatError("Problem type mismatch");
     134
     135    std::vector<typename BpGraph::Node> nodes(desc.nodeNum + 1, INVALID);
     136    typename BpGraph::Edge e;
     137    std::string problem, str;
     138    char c;
     139    int i, j;
     140
     141    // eat up comment lines and problem line
     142    while (is >> c && (c=='c' || c=='p')) {
     143      getline(is, str);
     144    }
     145    is.unget();
     146
     147    // reading red nodes
     148    while (is >> c && c == 'n') {
     149      is >> i;
     150      getline(is, str);
     151      nodes[i] = g.addRedNode();
     152    }
     153    is.unget();
     154
     155    // remaining nodes are blue:
     156    for (int k = 1; k <= desc.nodeNum; ++k) {
     157      if (nodes[k] == INVALID) nodes[k] = g.addBlueNode();
     158    }
     159
     160    // reading edges
     161    while (is >> c && c == 'a') {
     162      is >> i >> j;
     163      getline(is, str);
     164      e = g.addEdge(g.asRedNode(nodes[i]), g.asBlueNode(nodes[j]));
     165    }
     166  }
     167
     168
     169  /// \brief DIMACS weighted bipartite matching problem reader function.
     170  ///
     171  /// This function reads a weighted bipartite matching problem instance from
     172  /// DIMACS format, i.e. from a DIMACS file having a line starting with
     173  /// \code
     174  ///   p wbm
     175  /// \endcode
     176  /// and creates the corresponding bipartite graph and weight map.
     177  ///
     178  /// At the beginning, \c g is cleared by \c g.clear().
     179  ///
     180  /// If the file type was previously evaluated by dimacsType(), then
     181  /// the descriptor struct should be given by the \c dest parameter.
     182  template <typename BpGraph, typename WeightMap>
     183  void readDimacsWbm(std::istream& is,
     184                     BpGraph &g,
     185                     WeightMap &m,
     186                     DimacsDescriptor desc=DimacsDescriptor())
     187  {
     188    g.clear();
     189
     190    if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
     191    if(desc.type!=DimacsDescriptor::WBM)
     192      throw FormatError("Problem type mismatch");
     193
     194    std::vector<typename BpGraph::Node> nodes(desc.nodeNum + 1, INVALID);
     195    typename BpGraph::Edge e;
     196    std::string problem, str;
     197    char c;
     198    int i, j;
     199    typename WeightMap::Value w;
     200
     201    // eat up comment lines and problem line
     202    while (is >> c && (c=='c' || c=='p')) {
     203      getline(is, str);
     204    }
     205    is.unget();
     206
     207    // reading red nodes
     208    while (is >> c && c == 'n') {
     209      is >> i;
     210      getline(is, str);
     211      nodes[i] = g.addRedNode();
     212    }
     213    is.unget();
     214
     215    // remaining nodes are blue:
     216    for (int k = 1; k <= desc.nodeNum; ++k) {
     217      if (nodes[k] == INVALID) nodes[k] = g.addBlueNode();
     218    }
     219
     220    // reading edges and weights
     221    while (is >> c && c == 'a') {
     222      is >> i >> j >> w;
     223      getline(is, str);
     224      e = g.addEdge(g.asRedNode(nodes[i]), g.asBlueNode(nodes[j]));
     225      m.set(e, w);
     226    }
     227  }
     228
     229
     230  /// \brief DIMACS assignment problem reader function.
     231  ///
     232  /// This function reads an assignment problem instance from DIMACS format,
     233  /// i.e. from a DIMACS file having a line starting with
     234  /// \code
     235  ///   p asn
     236  /// \endcode
     237  /// and creates the corresponding digraph, arc cost- and supply map.
     238  ///
     239  /// At the beginning, \c g is cleared by \c g.clear().
     240  /// Cost of the arcs are written to the cost map,
     241  /// and the supply map contains the supply of the nodes, which is
     242  /// 1 for sources and -1 for sink nodes.
     243  ///
     244  /// If the file type was previously evaluated by dimacsType(), then
     245  /// the descriptor struct should be given by the \c dest parameter.
     246  template <typename Digraph, typename CostMap>
     247  void readDimacsAsn(std::istream& is,
     248                     Digraph &g,
     249                     CostMap &cost,
     250                     typename Digraph::template NodeMap<int> &supply,
     251                     DimacsDescriptor desc=DimacsDescriptor())
     252  {
     253    g.clear();
     254
     255    if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
     256    if(desc.type!=DimacsDescriptor::ASN)
     257      throw FormatError("Problem type mismatch");
     258
     259    std::vector<typename Digraph::Node> nodes(desc.nodeNum + 1, INVALID);
     260    typename Digraph::Arc e;
     261    std::string problem, str;
     262    char c;
     263    int i, j;
     264    typename CostMap::Value co;
     265
     266    for (int k = 1; k <= desc.nodeNum; ++k) {
     267      nodes[k] = g.addNode();
     268      supply.set(nodes[k], -1);
     269    }
     270
     271    // eat up comment lines and problem line
     272    while (is >> c && (c=='c' || c=='p')) {
     273      getline(is, str);
     274    }
     275    is.unget();
     276
     277    // reading source nodes
     278    while (is >> c && c == 'n') {
     279      is >> i;
     280      getline(is, str);
     281      supply.set(nodes[i], 1);
     282    }
     283    is.unget();
     284
     285    // reading arcs
     286    while (is >> c && c == 'a') {
     287      is >> i >> j >> co;
     288      getline(is, str);
     289      e = g.addArc(nodes[i], nodes[j]);
     290      cost.set(e, co);
     291    }
     292  }
     293
    104294  /// \brief DIMACS minimum cost flow reader function.
    105295  ///
    106296  /// This function reads a minimum cost flow instance from DIMACS format,
  • tools/CMakeLists.txt

    diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt
    a b  
    1616ADD_EXECUTABLE(dimacs-solver dimacs-solver.cc)
    1717TARGET_LINK_LIBRARIES(dimacs-solver lemon)
    1818
     19ADD_EXECUTABLE(weighted-bp-gen weighted-bp-gen.cc)
     20TARGET_LINK_LIBRARIES(weighted-bp-gen lemon)
     21
     22ADD_EXECUTABLE(unweighted-bp-gen unweighted-bp-gen.cc)
     23TARGET_LINK_LIBRARIES(unweighted-bp-gen lemon)
     24
     25ADD_EXECUTABLE(bp-matching-benchmark bp-matching-benchmark.cc)
     26TARGET_LINK_LIBRARIES(bp-matching-benchmark lemon)
     27
    1928INSTALL(
    2029  TARGETS lgf-gen dimacs-to-lgf dimacs-solver
    2130  RUNTIME DESTINATION bin
  • tools/Makefile.am

    diff --git a/tools/Makefile.am b/tools/Makefile.am
    a b  
    44if WANT_TOOLS
    55
    66bin_PROGRAMS += \
     7        tools/bp-matching-benchmark \
    78        tools/dimacs-solver \
    89        tools/dimacs-to-lgf \
    9         tools/lgf-gen
     10        tools/lgf-gen \
     11        tools/unweighted-bp-gen \
     12        tools/weighted-bp-gen
    1013
    1114dist_bin_SCRIPTS += tools/lemon-0.x-to-1.x.sh
    1215
    1316endif WANT_TOOLS
    1417
     18tools_bp_matching_benchmark_SOURCES = tools/bp-matching-benchmark.cc
    1519tools_dimacs_solver_SOURCES = tools/dimacs-solver.cc
    1620tools_dimacs_to_lgf_SOURCES = tools/dimacs-to-lgf.cc
    1721tools_lgf_gen_SOURCES = tools/lgf-gen.cc
     22tools_unweighted_bp_gen_SOURCES = tools/unweighted-bp-gen.cc
     23tools_weighted_bp_gen_SOURCES = tools/weighted-bp-gen.cc
  • new file tools/bp-matching-benchmark.cc

    diff --git a/tools/bp-matching-benchmark.cc b/tools/bp-matching-benchmark.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <lemon/dimacs.h>
     20#include <lemon/maps.h>
     21#include <lemon/smart_graph.h>
     22#include <lemon/time_measure.h>
     23#include <lemon/hopcroft_karp.h>
     24#include <lemon/matching.h>
     25#include <lemon/preflow.h>
     26#include <lemon/preflow_bp_matching.h>
     27#include <lemon/preflow_bp_matching_2.h>
     28#include <lemon/bp_matching.h>
     29#include <lemon/cost_scaling.h>
     30#include <lemon/capacity_scaling.h>
     31#include <lemon/network_simplex.h>
     32#include <iostream>
     33
     34using namespace std;
     35using namespace lemon;
     36
     37void check(bool c, const string &msg) {
     38  if (!c) {
     39    std::cerr << "Error: " << msg << "\n";
     40  }
     41}
     42
     43template <typename BpGraph, typename Digraph>
     44void create_network(const BpGraph &bpg,
     45                    Digraph &network,
     46                    typename Digraph::Node &source,
     47                    typename Digraph::Node &sink)
     48{
     49  network.clear();
     50  source = network.addNode();
     51  sink = network.addNode();
     52
     53  typename BpGraph::template NodeMap<typename Digraph::Node> xref(bpg);
     54  for (typename BpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     55    xref.set(r, network.addNode());
     56    network.addArc(source, xref[r]);
     57  }
     58  for (typename BpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     59    xref.set(b, network.addNode());
     60    network.addArc(xref[b], sink);
     61  }
     62
     63  for (typename BpGraph::EdgeIt e(bpg); e!=INVALID; ++e) {
     64    network.addArc(xref[bpg.redNode(e)], xref[bpg.blueNode(e)]);
     65  }
     66}
     67
     68template <typename BpGraph, typename BpCostMap,
     69          typename Digraph, typename DgCostMap, typename SupplyMap>
     70void create_weighted_net(const BpGraph &bpg,
     71                         const BpCostMap &bpcost,
     72                         Digraph &dg,
     73                         DgCostMap &dgcost,
     74                         SupplyMap &supply)
     75{
     76  dg.clear();
     77  typename Digraph::Node source = dg.addNode();
     78  supply.set(source, countBlueNodes(bpg) - countRedNodes(bpg));
     79
     80  typename BpGraph::template NodeMap<typename Digraph::Node> xref(bpg);
     81
     82  for (typename BpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     83    xref.set(r, dg.addNode());
     84    dgcost.set(dg.addArc(xref[r], source), 0);
     85    supply.set(xref[r], 1);
     86  }
     87  for (typename BpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     88    xref.set(b, dg.addNode());
     89    dgcost.set(dg.addArc(source, xref[b]), 0);
     90    supply.set(xref[b], -1);
     91  }
     92
     93  for (typename BpGraph::EdgeIt e(bpg); e!=INVALID; ++e) {
     94    typename Digraph::Arc a =
     95      dg.addArc(xref[bpg.redNode(e)], xref[bpg.blueNode(e)]);
     96    dgcost.set(a, -bpcost[e]);
     97  }
     98}
     99
     100typedef SmartDigraph Digraph;
     101typedef SmartBpGraph BpGraph;
     102
     103///\ingroup tools
     104///\file
     105///\brief Benchmark program for weighted and unweighted bipartite matching.
     106///
     107/// The program input is a Dimacs file with a problem line starting
     108/// \code
     109///   p ubm
     110/// \endcode
     111/// or
     112/// \code
     113///   p wbm
     114/// \endcode
     115///
     116/// Depending on the input, it runs the \ref lemon::HopcroftKarp
     117/// "Hopcroft-Karp", the \ref lemon::MaxMatching "general matching"
     118/// and three algorithms based on preflow (\ref lemon::Preflow "preflow
     119/// on extended graph", \ref lemon::PreflowBpMatching,
     120/// \ref lemon::PreflowBpMatching2") for finding
     121/// a maximum cardinality bipartite matching, or the
     122/// \ref lemon::MaxWeightedBpMatching "max. weighted matching
     123/// for sparse bipartite graphs", the \ref lemon::MaxWeightedMatching
     124/// "max. weighted matching for general graph", \ref
     125/// lemon::CostScaling "cost scaling", \ref lemon::CapacityScaling
     126/// "capacity scaling" and \ref lemon::NetworkSimplex "network
     127/// simplex" algorithms to find the maximum weighted matching.
     128int main() {
     129  DimacsDescriptor desc = dimacsType(cin);
     130  if (desc.type == DimacsDescriptor::UBM) {
     131    // unweighted bipartite matching
     132    Timer hk_t(false),
     133      mm_t(false),
     134      pf_t(false),
     135      pbm_t(false),
     136      pbm2_t(false);
     137
     138    BpGraph *bpg = new BpGraph();
     139    readDimacsUbm(cin, *bpg, desc);
     140
     141    {
     142      hk_t.start();
     143      HopcroftKarp<BpGraph> hk(*bpg);
     144      hk.run();
     145      hk_t.stop();
     146      BpGraph::NodeMap<bool> m(*bpg);
     147      check(hk.matchingSize() == hk.cover(m),
     148            "Suboptimal matching in Hopcroft-Karp.");
     149    }
     150    {
     151      mm_t.start();
     152      MaxMatching<BpGraph> mm(*bpg);
     153      mm.run();
     154      mm_t.stop();
     155    }
     156    {
     157      pbm_t.start();
     158      PreflowBpMatching<BpGraph> pbm(*bpg);
     159      pbm.run();
     160      pbm_t.stop();
     161      BpGraph::NodeMap<bool> m(*bpg);
     162      check(pbm.matchingSize() == pbm.cover(m),
     163            "Suboptimal matching in push-relabel.");
     164    }
     165    {
     166      pbm2_t.start();
     167      PreflowBpMatching2<BpGraph> pbm2(*bpg);
     168      pbm2.run();
     169      pbm2_t.stop();
     170      BpGraph::NodeMap<bool> m(*bpg);
     171      check(pbm2.matchingSize() == pbm2.cover(m),
     172            "Suboptimal matching in pull-relabel.");
     173    }
     174
     175    Digraph *network = new Digraph();
     176    Digraph::Node source, target;
     177    create_network(*bpg, *network, source, target);
     178    int rednum = countRedNodes(*bpg),
     179      bluenum = countBlueNodes(*bpg),
     180      edgenum = desc.edgeNum;
     181
     182    delete bpg;
     183
     184    {
     185      pf_t.start();
     186      Preflow<Digraph, ConstMap<Digraph::Arc, long> >
     187        pf(*network, constMap<Digraph::Arc, long>(1), source, target);
     188      pf.run();
     189      pf_t.stop();
     190    }
     191    delete network;
     192
     193    cout << "Benchmarking in unweighted case, using a bipartite graph\n"
     194         << "with " << rednum << " red, " << bluenum << " blue nodes, and "
     195         << edgenum << " edges.\n"
     196         << "--------------------------------------------------------\n"
     197         << "Algorithm used               Time\n"
     198         << "Hopcroft-Karp            " << hk_t.realTime() << "\n"
     199         << "General matching         " << mm_t.realTime() << "\n"
     200         << "Push-relabel matching    " << pbm_t.realTime() << "\n"
     201         << "Pull-relabel matching    " << pbm2_t.realTime() << "\n"
     202         << "Preflow                  " << pf_t.realTime() << endl;
     203
     204  } else if (desc.type == DimacsDescriptor::WBM) {
     205    // weighted bipartite matching
     206    Timer bm_t(false),
     207      gm_t(false),
     208      cas_t(false),
     209      cos_t(false),
     210      ns_t(false);
     211
     212    BpGraph *bpg = new BpGraph();
     213    BpGraph::EdgeMap<long> *weight = new BpGraph::EdgeMap<long>(*bpg);
     214    readDimacsWbm(cin, *bpg, *weight, desc);
     215
     216    {
     217      bm_t.start();
     218      MaxWeightedBpMatching<BpGraph, BpGraph::EdgeMap<long> > bm(*bpg, *weight);
     219      bm.run();
     220      bm_t.stop();
     221    }
     222    {
     223      gm_t.start();
     224      MaxWeightedMatching<BpGraph, BpGraph::EdgeMap<long> > gm(*bpg, *weight);
     225      gm.run();
     226      gm_t.stop();
     227    }
     228
     229    Digraph *g = new Digraph();
     230    Digraph::ArcMap<long> *cost = new Digraph::ArcMap<long>(*g);
     231    Digraph::NodeMap<long> *supply = new Digraph::NodeMap<long>(*g);
     232    create_weighted_net(*bpg, *weight, *g, *cost, *supply);
     233
     234    int rednum = countRedNodes(*bpg),
     235      bluenum = countBlueNodes(*bpg),
     236      edgenum = desc.edgeNum;
     237
     238    delete bpg;
     239    delete weight;
     240    {
     241      cas_t.start();
     242      CapacityScaling<Digraph> cas(*g);
     243      cas.costMap(*cost)
     244        .supplyMap(*supply)
     245        .upperMap(constMap<Digraph::Arc, long>(1));
     246      cas.run();
     247      cas_t.stop();
     248    }
     249    {
     250      cos_t.start();
     251      CostScaling<Digraph> cos(*g);
     252      cos.costMap(*cost)
     253        .supplyMap(*supply)
     254        .upperMap(constMap<Digraph::Arc, long>(1));
     255      cos.run();
     256      cos_t.stop();
     257    }
     258    {
     259      ns_t.start();
     260      NetworkSimplex<Digraph> ns(*g);
     261      ns.costMap(*cost)
     262        .supplyMap(*supply)
     263        .upperMap(constMap<Digraph::Arc, long>(1));
     264      ns.run();
     265      ns_t.stop();
     266    }
     267    delete g;
     268    delete cost;
     269    delete supply;
     270
     271    cout << "Benchmarking on weighted bipartite matching problem on a graph\n"
     272         << "with " << rednum << " red, " << bluenum << " blue nodes and "
     273         << edgenum << " edges\n"
     274         << "--------------------------------------------------------\n"
     275         << "Algorithm used             Time\n"
     276         << "Bipartite matching       " << bm_t.realTime() << "\n"
     277         << "General matching         " << gm_t.realTime() << "\n"
     278         << "Capacity scaling         " << cas_t.realTime() << "\n"
     279         << "Cost scaling             " << cos_t.realTime() << "\n"
     280         << "Network simplex          " << ns_t.realTime() << endl;
     281
     282  } else {
     283    cerr << "Wrong problem type." << endl;
     284    return 1;
     285  }
     286
     287  return 0;
     288}
  • new file tools/unweighted-bp-gen.cc

    diff --git a/tools/unweighted-bp-gen.cc b/tools/unweighted-bp-gen.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <lemon/arg_parser.h>
     20#include <lemon/random.h>
     21
     22using namespace lemon;
     23using namespace std;
     24
     25///\ingroup tools
     26///\file
     27///\brief Random bipartite graph generator.
     28///
     29/// This program generates an unweighted bipartite matching problem.
     30/// The output is in Dimacs format, with problem type of 'ubm'. The
     31/// enumerated nodes form one set of the bipartite graph, and
     32/// the edges do not have any label.
     33int main(int argc, char** argv) {
     34  int seed = 1,
     35    red = 5,
     36    blue = 5,
     37    density = 10;
     38
     39  ArgParser arg_p(argc, argv);
     40
     41  arg_p.refOption("seed",
     42                  "Seed of the random number generator (default: 1)",
     43                  seed)
     44    .refOption("red",
     45               "Number of red nodes (default: 5)",
     46               red)
     47    .refOption("blue",
     48               "Number of blue nodes (default: 5)",
     49               blue)
     50    .refOption("density",
     51               "Number of edges (default: 10)",
     52               density);
     53
     54  arg_p.synonym("s", "seed")
     55    .synonym("r", "red")
     56    .synonym("b", "blue")
     57    .synonym("d", "density");
     58
     59  arg_p.parse();
     60
     61  rnd.seed(seed);
     62
     63
     64  long max_edge = static_cast<long>(red)*blue;
     65  long edge_left = density;
     66
     67  if (max_edge < density) {
     68    cerr << "Too many edges!\n"
     69         << "At most " << max_edge << " edges can be placed in the graph.\n";
     70    return 1;
     71  }
     72
     73  cout << "c Random bipartite graph for bipartite matching\n"
     74       << "c----------------------------------------------\n"
     75       << "c Input parameters:\n"
     76       << "c   Random seed: " << seed << "\n"
     77       << "c   Number of nodes: " << red + blue << "\n"
     78       << "c     Number of red nodes:  " << red << "\n"
     79       << "c     Number of blue nodes: " << blue << "\n"
     80       << "c   Number of edges: " << density << "\n"
     81       << "c----------------------------------------------\n";
     82  cout << "p ubm " << red + blue << " " << density << "\n";
     83
     84  // red nodes:
     85  for (int r=1; r <= red; ++r) {
     86    cout << "n " << r << "\n";
     87  }
     88
     89  // edges:
     90  for (int r = 1; r <= red; ++r) {
     91    for (int b = red+1; b <= red + blue; ++b) {
     92      if (rnd.boolean(static_cast<long double>(edge_left) / max_edge)) {
     93        --edge_left;
     94        cout << "a " << r << " " << b << "\n";
     95      }
     96      --max_edge;
     97    }
     98  }
     99
     100  return 0;
     101}
     102
  • new file tools/weighted-bp-gen.cc

    diff --git a/tools/weighted-bp-gen.cc b/tools/weighted-bp-gen.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <lemon/arg_parser.h>
     20#include <lemon/random.h>
     21
     22using namespace lemon;
     23using namespace std;
     24
     25///\ingroup tools
     26///\file
     27///\brief Random weighted bipartite graph generator.
     28///
     29/// This program generates a weighted bipartite matching problem.
     30/// The output is in Dimacs format, with problem type of 'wbm'. The
     31/// enumerated nodes form one set of the bipartite graph, and
     32/// the edge labels denote weights.
     33int main(int argc, char** argv) {
     34  int seed = 1,
     35    red = 5,
     36    blue = 5,
     37    density = 10,
     38    minweight = 10,
     39    maxweight = 100;
     40
     41  ArgParser arg_p(argc, argv);
     42
     43  arg_p.refOption("seed",
     44                  "Seed of the random number generator (default: 1)",
     45                  seed)
     46    .refOption("red",
     47               "Number of red nodes (default: 5)",
     48               red)
     49    .refOption("blue",
     50               "Number of blue nodes (default: 5)",
     51               blue)
     52    .refOption("density",
     53               "Number of edges (default: 10)",
     54               density)
     55    .refOption("minweight",
     56               "Minimal weight of edges (default: 10)",
     57               minweight)
     58    .refOption("maxweight",
     59               "Maximal weight of edges (default: 100",
     60               maxweight);
     61
     62  arg_p.synonym("s", "seed")
     63    .synonym("r", "red")
     64    .synonym("b", "blue")
     65    .synonym("d", "density")
     66    .synonym("min", "minweight")
     67    .synonym("max", "maxweight");
     68
     69  arg_p.parse();
     70
     71  rnd.seed(seed);
     72
     73
     74  long max_edge = static_cast<long>(red)*blue;
     75  long edge_left = density;
     76
     77  if (max_edge < density) {
     78    cerr << "Too many edges!\n"
     79         << "At most " << max_edge << " edges can be placed in the graph.\n";
     80    return 1;
     81  }
     82
     83  cout << "c Random bipartite graph for weighted bipartite matching\n"
     84       << "c----------------------------------------------\n"
     85       << "c Input parameters:\n"
     86       << "c   Random seed: " << seed << "\n"
     87       << "c   Number of nodes: " << red + blue << "\n"
     88       << "c     Number of red nodes:  " << red << "\n"
     89       << "c     Number of blue nodes: " << blue << "\n"
     90       << "c   Number of edges: " << density << "\n"
     91       << "c   Minimal weight of edges: " << minweight << "\n"
     92       << "c   Maximal weight of edges: " << maxweight << "\n"
     93       << "c----------------------------------------------\n";
     94  cout << "p wbm " << red + blue << " " << density << "\n";
     95
     96  // red nodes:
     97  for (int r=1; r <= red; ++r) {
     98    cout << "n " << r << "\n";
     99  }
     100
     101  // edges:
     102  for (int r = 1; r <= red; ++r) {
     103    for (int b = red+1; b <= red + blue; ++b) {
     104      if (rnd.boolean(static_cast<long double>(edge_left) / max_edge)) {
     105        --edge_left;
     106        cout << "a " << r << " " << b << " "
     107             << rnd.integer(minweight, maxweight+1) << "\n";
     108      }
     109      --max_edge;
     110    }
     111  }
     112
     113  return 0;
     114}
     115