COIN-OR::LEMON - Graph Library

About

Downloads

Documentation

Development

Issue Tracker

Browse Tickets

Custom Query

Roadmap

Browse Source

Local Info

Mailing Lists

Search

Changeset 485:7f461ab4af1a in lemon-0.x


Ignore:
Timestamp:
04/29/04 19:34:42 (20 years ago)
Author:
marci
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@642
Message:

Some docu in MaxFlow? class, jacint/max_flow.h

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/work/jacint/max_flow.h

    r480 r485  
    22
    33/*
    4 Heuristics:
    5  2 phase
    6  gap
    7  list 'level_list' on the nodes on level i implemented by hand
    8  stack 'active' on the active nodes on level i
    9  runs heuristic 'highest label' for H1*n relabels
    10  runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
     4  Heuristics:
     5  2 phase
     6  gap
     7  list 'level_list' on the nodes on level i implemented by hand
     8  stack 'active' on the active nodes on level i
     9  runs heuristic 'highest label' for H1*n relabels
     10  runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
    1111 
    12 Parameters H0 and H1 are initialized to 20 and 1.
    13 
    14 Constructors:
    15 
    16 Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
    17      FlowMap is not constant zero, and should be true if it is
    18 
    19 Members:
    20 
    21 void run()
    22 
    23 Num flowValue() : returns the value of a maximum flow
    24 
    25 void minMinCut(CutMap& M) : sets M to the characteristic vector of the
    26      minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
    27 
    28 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
    29      maximum min cut. M should be a map of bools initialized to false.
    30 
    31 void minCut(CutMap& M) : sets M to the characteristic vector of
    32      a min cut. M should be a map of bools initialized to false.
     12  Parameters H0 and H1 are initialized to 20 and 1.
     13
     14  Constructors:
     15
     16  Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
     17  FlowMap is not constant zero, and should be true if it is
     18
     19  Members:
     20
     21  void run()
     22
     23  Num flowValue() : returns the value of a maximum flow
     24
     25  void minMinCut(CutMap& M) : sets M to the characteristic vector of the
     26  minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
     27
     28  void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
     29  maximum min cut. M should be a map of bools initialized to false.
     30
     31  void minCut(CutMap& M) : sets M to the characteristic vector of
     32  a min cut. M should be a map of bools initialized to false.
    3333
    3434*/
    … …  
    5353namespace hugo {
    5454
     55  ///\author Marton Makai, Jacint Szabo
    5556  template <typename Graph, typename Num,
    5657            typename CapMap=typename Graph::template EdgeMap<Num>,
    … …  
    9899      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
    99100
     101    /// A max flow algorithm is run.
     102    ///\pre the flow have to be 0 at the beginning.
    100103    void run() {
    101104      preflow( ZERO_FLOW );
    102105    }
    103106   
     107    /// A preflow algorithm is run. The initial edge-set have to be a flow,
     108    /// or from a preflow, according to \c fe.
    104109    void preflow( flowEnum fe ) {
    105110      preflowPhase0(fe);
    … …  
    107112    }
    108113
     114    /// Run the first phase of preflow, starting from a 0 flow, from a flow,
     115    /// or from a preflow, according to \c fe.
    109116    void preflowPhase0( flowEnum fe );
    110117
     118    /// Second phase of preflow.
    111119    void preflowPhase1();
    112120
     121    /// Starting from a flow, this method searches for an augmenting path
     122    /// according to the Edmonds-Karp algorithm
     123    /// and augments the flow on if any.
    113124    bool augmentOnShortestPath();
    114125
     126    /// Starting from a flow, this method searches for an augmenting blockin
     127    /// flow according to Dinits' algorithm and augments the flow on if any.
     128    /// The blocking flow is computed in a physically constructed
     129    /// residual graph of type \c Mutablegraph.
    115130    template<typename MutableGraph> bool augmentOnBlockingFlow();
    116131
     132    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
     133    /// residual graph is not constructed physically.
    117134    bool augmentOnBlockingFlow2();
    118135
    … …  
    127144    }
    128145
    129     //should be used only between preflowPhase0 and preflowPhase1
     146    /// Should be used between preflowPhase0 and preflowPhase1.
     147    ///\todo We have to make some status variable which shows the actual state
     148    /// of the class. This enables us to determine which methods are valid
     149    /// for MinCut computation
    130150    template<typename _CutMap>
    131151    void actMinCut(_CutMap& M) {
    132152      NodeIt v;
    133       for(g->first(v); g->valid(v); g->next(v))
    134       if ( level[v] < n ) {
    135         M.set(v,false);
    136       } else {
    137         M.set(v,true);
    138       }
    139     }
    140 
    141 
    142 
    143     /*
    144       Returns the minimum min cut, by a bfs from s in the residual graph.
    145     */
     153      for(g->first(v); g->valid(v); g->next(v)) {
     154        if ( level[v] < n ) {
     155          M.set(v,false);
     156        } else {
     157          M.set(v,true);
     158        }
     159      }
     160    }
     161
     162
     163    /// The unique inclusionwise minimum cut is computed by
     164    /// processing a bfs from s in the residual graph.
     165    ///\pre flow have to be a max flow otherwise it will the whole node-set.
    146166    template<typename _CutMap>
    147167    void minMinCut(_CutMap& M) {
    … …  
    177197
    178198
    179  
    180     /*
    181       Returns the maximum min cut, by a reverse bfs
    182       from t in the residual graph.
    183     */
    184    
     199    /// The unique inclusionwise maximum cut is computed by
     200    /// processing a reverse bfs from t in the residual graph.
     201    ///\pre flow have to be a max flow otherwise it will be empty.
    185202    template<typename _CutMap>
    186203    void maxMinCut(_CutMap& M) {
    … …  
    222239
    223240
     241    /// A minimum cut is computed.
    224242    template<typename CutMap>
    225     void minCut(CutMap& M) {
    226       minMinCut(M);
    227     }
    228 
     243    void minCut(CutMap& M) { minMinCut(M); }
     244
     245    ///
     246    void resetSource(Node _s) {s=_s;}
     247    ///
    229248    void resetTarget(Node _t) {t=_t;}
    230     void resetSource(Node _s) {s=_s;}
    231249   
    232     void resetCap(const CapMap& _cap) {
    233       capacity=&_cap;
    234     }
     250    /// capacity-map is changed.
     251    void resetCap(const CapMap& _cap) { capacity=&_cap; }
    235252   
    236     void resetFlow(FlowMap& _flow) {
    237       flow=&_flow;
    238     }
     253    /// flow-map is changed.
     254    void resetFlow(FlowMap& _flow) { flow=&_flow; }
    239255
    240256
    … …  
    315331     
    316332      return newlevel;
    317      }
     333    }
    318334
    319335
    … …  
    321337                          VecNode& level_list, NNMap& left, NNMap& right ) {
    322338
    323       std::queue<Node> bfs_queue;
    324      
    325       switch ( fe ) {
    326       case ZERO_FLOW:
    327         {
    328           //Reverse_bfs from t, to find the starting level.
    329           level.set(t,0);
    330           bfs_queue.push(t);
    331        
    332           while (!bfs_queue.empty()) {
    333            
    334             Node v=bfs_queue.front();   
    335             bfs_queue.pop();
    336             int l=level[v]+1;
    337            
    338             InEdgeIt e;
    339             for(g->first(e,v); g->valid(e); g->next(e)) {
    340               Node w=g->tail(e);
    341               if ( level[w] == n && w != s ) {
    342                 bfs_queue.push(w);
    343                 Node first=level_list[l];
    344                 if ( g->valid(first) ) left.set(first,w);
    345                 right.set(w,first);
    346                 level_list[l]=w;
    347                 level.set(w, l);
    348               }
    349             }
    350           }
    351          
    352           //the starting flow
    353           OutEdgeIt e;
    354           for(g->first(e,s); g->valid(e); g->next(e))
    355             {
    356               Num c=(*capacity)[e];
    357               if ( c <= 0 ) continue;
    358               Node w=g->head(e);
    359               if ( level[w] < n ) {       
    360                 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
    361                 flow->set(e, c);
    362                 excess.set(w, excess[w]+c);
    363               }
    364             }
    365           break;
    366         }
    367        
    368       case GEN_FLOW:
    369       case PREFLOW:
    370         {
    371           //Reverse_bfs from t in the residual graph,
    372           //to find the starting level.
    373           level.set(t,0);
    374           bfs_queue.push(t);
    375          
    376           while (!bfs_queue.empty()) {
    377            
    378             Node v=bfs_queue.front();   
    379             bfs_queue.pop();
    380             int l=level[v]+1;
    381            
    382             InEdgeIt e;
    383             for(g->first(e,v); g->valid(e); g->next(e)) {
    384               if ( (*capacity)[e] <= (*flow)[e] ) continue;
    385               Node w=g->tail(e);
    386               if ( level[w] == n && w != s ) {
    387                 bfs_queue.push(w);
    388                 Node first=level_list[l];
    389                 if ( g->valid(first) ) left.set(first,w);
    390                 right.set(w,first);
    391                 level_list[l]=w;
    392                 level.set(w, l);
    393               }
    394             }
    395            
    396             OutEdgeIt f;
    397             for(g->first(f,v); g->valid(f); g->next(f)) {
    398               if ( 0 >= (*flow)[f] ) continue;
    399               Node w=g->head(f);
    400               if ( level[w] == n && w != s ) {
    401                 bfs_queue.push(w);
    402                 Node first=level_list[l];
    403                 if ( g->valid(first) ) left.set(first,w);
    404                 right.set(w,first);
    405                 level_list[l]=w;
    406                 level.set(w, l);
    407               }
    408             }
    409           }
    410          
    411          
    412           //the starting flow
    413           OutEdgeIt e;
    414           for(g->first(e,s); g->valid(e); g->next(e))
    415             {
    416               Num rem=(*capacity)[e]-(*flow)[e];
    417               if ( rem <= 0 ) continue;
    418               Node w=g->head(e);
    419               if ( level[w] < n ) {       
    420                 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
    421                 flow->set(e, (*capacity)[e]);
    422                 excess.set(w, excess[w]+rem);
    423               }
    424             }
    425          
    426           InEdgeIt f;
    427           for(g->first(f,s); g->valid(f); g->next(f))
    428             {
    429               if ( (*flow)[f] <= 0 ) continue;
    430               Node w=g->tail(f);
    431               if ( level[w] < n ) {       
    432                 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
    433                 excess.set(w, excess[w]+(*flow)[f]);
    434                 flow->set(f, 0);
    435               }
    436             } 
    437           break;
    438         } //case PREFLOW
    439       }
    440     } //preflowPreproc
     339                            std::queue<Node> bfs_queue;
     340     
     341                            switch ( fe ) {
     342                            case ZERO_FLOW:
     343                              {
     344                                //Reverse_bfs from t, to find the starting level.
     345                                level.set(t,0);
     346                                bfs_queue.push(t);
     347       
     348                                while (!bfs_queue.empty()) {
     349           
     350                                  Node v=bfs_queue.front();     
     351                                  bfs_queue.pop();
     352                                  int l=level[v]+1;
     353           
     354                                  InEdgeIt e;
     355                                  for(g->first(e,v); g->valid(e); g->next(e)) {
     356                                    Node w=g->tail(e);
     357                                    if ( level[w] == n && w != s ) {
     358                                      bfs_queue.push(w);
     359                                      Node first=level_list[l];
     360                                      if ( g->valid(first) ) left.set(first,w);
     361                                      right.set(w,first);
     362                                      level_list[l]=w;
     363                                      level.set(w, l);
     364                                    }
     365                                  }
     366                                }
     367         
     368                                //the starting flow
     369                                OutEdgeIt e;
     370                                for(g->first(e,s); g->valid(e); g->next(e))
     371                                  {
     372                                    Num c=(*capacity)[e];
     373                                    if ( c <= 0 ) continue;
     374                                    Node w=g->head(e);
     375                                    if ( level[w] < n ) {         
     376                                      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
     377                                      flow->set(e, c);
     378                                      excess.set(w, excess[w]+c);
     379                                    }
     380                                  }
     381                                break;
     382                              }
     383       
     384                            case GEN_FLOW:
     385                            case PREFLOW:
     386                              {
     387                                //Reverse_bfs from t in the residual graph,
     388                                //to find the starting level.
     389                                level.set(t,0);
     390                                bfs_queue.push(t);
     391         
     392                                while (!bfs_queue.empty()) {
     393           
     394                                  Node v=bfs_queue.front();     
     395                                  bfs_queue.pop();
     396                                  int l=level[v]+1;
     397           
     398                                  InEdgeIt e;
     399                                  for(g->first(e,v); g->valid(e); g->next(e)) {
     400                                    if ( (*capacity)[e] <= (*flow)[e] ) continue;
     401                                    Node w=g->tail(e);
     402                                    if ( level[w] == n && w != s ) {
     403                                      bfs_queue.push(w);
     404                                      Node first=level_list[l];
     405                                      if ( g->valid(first) ) left.set(first,w);
     406                                      right.set(w,first);
     407                                      level_list[l]=w;
     408                                      level.set(w, l);
     409                                    }
     410                                  }
     411           
     412                                  OutEdgeIt f;
     413                                  for(g->first(f,v); g->valid(f); g->next(f)) {
     414                                    if ( 0 >= (*flow)[f] ) continue;
     415                                    Node w=g->head(f);
     416                                    if ( level[w] == n && w != s ) {
     417                                      bfs_queue.push(w);
     418                                      Node first=level_list[l];
     419                                      if ( g->valid(first) ) left.set(first,w);
     420                                      right.set(w,first);
     421                                      level_list[l]=w;
     422                                      level.set(w, l);
     423                                    }
     424                                  }
     425                                }
     426         
     427         
     428                                //the starting flow
     429                                OutEdgeIt e;
     430                                for(g->first(e,s); g->valid(e); g->next(e))
     431                                  {
     432                                    Num rem=(*capacity)[e]-(*flow)[e];
     433                                    if ( rem <= 0 ) continue;
     434                                    Node w=g->head(e);
     435                                    if ( level[w] < n ) {         
     436                                      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
     437                                      flow->set(e, (*capacity)[e]);
     438                                      excess.set(w, excess[w]+rem);
     439                                    }
     440                                  }
     441         
     442                                InEdgeIt f;
     443                                for(g->first(f,s); g->valid(f); g->next(f))
     444                                  {
     445                                    if ( (*flow)[f] <= 0 ) continue;
     446                                    Node w=g->tail(f);
     447                                    if ( level[w] < n ) {         
     448                                      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
     449                                      excess.set(w, excess[w]+(*flow)[f]);
     450                                      flow->set(f, 0);
     451                                    }
     452                                  } 
     453                                break;
     454                              } //case PREFLOW
     455                            }
     456                          } //preflowPreproc
    441457
    442458
    … …  
    522538      }
    523539      int operator[](const typename MapGraphWrapper::Node& n)
    524         { return dist[n]; }
    525 //       int get(const typename MapGraphWrapper::Node& n) const {
    526 //      return dist[n]; }
    527 //       bool get(const typename MapGraphWrapper::Edge& e) const {
    528 //      return (dist.get(g->tail(e))<dist.get(g->head(e))); }
     540      { return dist[n]; }
     541      //       int get(const typename MapGraphWrapper::Node& n) const {
     542      //        return dist[n]; }
     543      //       bool get(const typename MapGraphWrapper::Edge& e) const {
     544      //        return (dist.get(g->tail(e))<dist.get(g->head(e))); }
    529545      bool operator[](const typename MapGraphWrapper::Edge& e) const {
    530546        return (dist[g->tail(e)]<dist[g->head(e)]);
    … …  
    539555  {
    540556     
    541       int heur0=(int)(H0*n);  //time while running 'bound decrease'
    542       int heur1=(int)(H1*n);  //time while running 'highest label'
    543       int heur=heur1;         //starting time interval (#of relabels)
    544       int numrelabel=0;
     557    int heur0=(int)(H0*n);  //time while running 'bound decrease'
     558    int heur1=(int)(H1*n);  //time while running 'highest label'
     559    int heur=heur1;         //starting time interval (#of relabels)
     560    int numrelabel=0;
    545561     
    546       bool what_heur=1;       
    547       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
    548 
    549       bool end=false;     
    550       //Needed for 'bound decrease', true means no active nodes are above bound b.
    551 
    552       int k=n-2;  //bound on the highest level under n containing a node
    553       int b=k;    //bound on the highest level under n of an active node
    554      
    555       VecStack active(n);
    556      
    557       NNMap left(*g, INVALID);
    558       NNMap right(*g, INVALID);
    559       VecNode level_list(n,INVALID);
    560       //List of the nodes in level i<n, set to n.
    561 
    562       NodeIt v;
    563       for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
    564       //setting each node to level n
    565      
    566       switch ( fe ) {
    567       case PREFLOW:
    568         {
    569           //counting the excess
    570           NodeIt v;
    571           for(g->first(v); g->valid(v); g->next(v)) {
    572             Num exc=0;
    573          
    574             InEdgeIt e;
    575             for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
    576             OutEdgeIt f;
    577             for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
    578            
    579             excess.set(v,exc);   
    580            
    581             //putting the active nodes into the stack
    582             int lev=level[v];
    583             if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
     562    bool what_heur=1;       
     563    //It is 0 in case 'bound decrease' and 1 in case 'highest label'
     564
     565    bool end=false;     
     566    //Needed for 'bound decrease', true means no active nodes are above bound b.
     567
     568    int k=n-2;  //bound on the highest level under n containing a node
     569    int b=k;    //bound on the highest level under n of an active node
     570     
     571    VecStack active(n);
     572     
     573    NNMap left(*g, INVALID);
     574    NNMap right(*g, INVALID);
     575    VecNode level_list(n,INVALID);
     576    //List of the nodes in level i<n, set to n.
     577
     578    NodeIt v;
     579    for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
     580    //setting each node to level n
     581     
     582    switch ( fe ) {
     583    case PREFLOW:
     584      {
     585        //counting the excess
     586        NodeIt v;
     587        for(g->first(v); g->valid(v); g->next(v)) {
     588          Num exc=0;
     589         
     590          InEdgeIt e;
     591          for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
     592          OutEdgeIt f;
     593          for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
     594           
     595          excess.set(v,exc);     
     596           
     597          //putting the active nodes into the stack
     598          int lev=level[v];
     599          if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
     600        }
     601        break;
     602      }
     603    case GEN_FLOW:
     604      {
     605        //Counting the excess of t
     606        Num exc=0;
     607         
     608        InEdgeIt e;
     609        for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
     610        OutEdgeIt f;
     611        for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
     612         
     613        excess.set(t,exc);     
     614         
     615        break;
     616      }
     617    default:
     618      break;
     619    }
     620     
     621    preflowPreproc( fe, active, level_list, left, right );
     622    //End of preprocessing
     623     
     624     
     625    //Push/relabel on the highest level active nodes.
     626    while ( true ) {
     627      if ( b == 0 ) {
     628        if ( !what_heur && !end && k > 0 ) {
     629          b=k;
     630          end=true;
     631        } else break;
     632      }
     633       
     634      if ( active[b].empty() ) --b;
     635      else {
     636        end=false; 
     637        Node w=active[b].top();
     638        active[b].pop();
     639        int newlevel=push(w,active);
     640        if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
     641                                     left, right, b, k, what_heur);
     642         
     643        ++numrelabel;
     644        if ( numrelabel >= heur ) {
     645          numrelabel=0;
     646          if ( what_heur ) {
     647            what_heur=0;
     648            heur=heur0;
     649            end=false;
     650          } else {
     651            what_heur=1;
     652            heur=heur1;
     653            b=k;
    584654          }
    585           break;
    586         }
    587       case GEN_FLOW:
    588         {
    589           //Counting the excess of t
    590           Num exc=0;
    591          
    592           InEdgeIt e;
    593           for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
    594           OutEdgeIt f;
    595           for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
    596          
    597           excess.set(t,exc);   
    598          
    599           break;
    600         }
    601       default:
    602         break;
    603       }
    604      
    605       preflowPreproc( fe, active, level_list, left, right );
    606       //End of preprocessing
    607      
    608      
    609       //Push/relabel on the highest level active nodes.
    610       while ( true ) {
    611         if ( b == 0 ) {
    612           if ( !what_heur && !end && k > 0 ) {
    613             b=k;
    614             end=true;
    615           } else break;
    616         }
    617        
    618         if ( active[b].empty() ) --b;
    619         else {
    620           end=false; 
    621           Node w=active[b].top();
    622           active[b].pop();
    623           int newlevel=push(w,active);
    624           if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
    625                                        left, right, b, k, what_heur);
    626          
    627           ++numrelabel;
    628           if ( numrelabel >= heur ) {
    629             numrelabel=0;
    630             if ( what_heur ) {
    631               what_heur=0;
    632               heur=heur0;
    633               end=false;
    634             } else {
    635               what_heur=1;
    636               heur=heur1;
    637               b=k;
    638             }
    639           }
    640         }
     655        }
    641656      }
    642     }
     657    }
     658  }
    643659
    644660
    … …  
    648664  {
    649665     
    650       int k=n-2;  //bound on the highest level under n containing a node
    651       int b=k;    //bound on the highest level under n of an active node
    652      
    653       VecStack active(n);
    654       level.set(s,0);
    655       std::queue<Node> bfs_queue;
    656       bfs_queue.push(s);
    657            
    658       while (!bfs_queue.empty()) {
    659        
    660         Node v=bfs_queue.front();       
    661         bfs_queue.pop();
    662         int l=level[v]+1;
     666    int k=n-2;  //bound on the highest level under n containing a node
     667    int b=k;    //bound on the highest level under n of an active node
     668     
     669    VecStack active(n);
     670    level.set(s,0);
     671    std::queue<Node> bfs_queue;
     672    bfs_queue.push(s);
     673           
     674    while (!bfs_queue.empty()) {
     675       
     676      Node v=bfs_queue.front();
     677      bfs_queue.pop();
     678      int l=level[v]+1;
    663679             
    664         InEdgeIt e;
    665         for(g->first(e,v); g->valid(e); g->next(e)) {
    666           if ( (*capacity)[e] <= (*flow)[e] ) continue;
    667           Node u=g->tail(e);
    668           if ( level[u] >= n ) {
    669             bfs_queue.push(u);
    670             level.set(u, l);
    671             if ( excess[u] > 0 ) active[l].push(u);
    672           }
    673         }
    674        
    675         OutEdgeIt f;
    676         for(g->first(f,v); g->valid(f); g->next(f)) {
    677           if ( 0 >= (*flow)[f] ) continue;
    678           Node u=g->head(f);
    679           if ( level[u] >= n ) {
    680             bfs_queue.push(u);
    681             level.set(u, l);
    682             if ( excess[u] > 0 ) active[l].push(u);
    683           }
    684         }
    685       }
    686       b=n-2;
    687 
    688       while ( true ) {
    689        
    690         if ( b == 0 ) break;
    691 
    692         if ( active[b].empty() ) --b;
    693         else {
    694           Node w=active[b].top();
    695           active[b].pop();
    696           int newlevel=push(w,active);   
    697 
    698           //relabel
    699           if ( excess[w] > 0 ) {
    700             level.set(w,++newlevel);
    701             active[newlevel].push(w);
    702             b=newlevel;
    703           }
    704         }  // if stack[b] is nonempty
    705       } // while(true)
    706     }
     680      InEdgeIt e;
     681      for(g->first(e,v); g->valid(e); g->next(e)) {
     682        if ( (*capacity)[e] <= (*flow)[e] ) continue;
     683        Node u=g->tail(e);
     684        if ( level[u] >= n ) {
     685          bfs_queue.push(u);
     686          level.set(u, l);
     687          if ( excess[u] > 0 ) active[l].push(u);
     688        }
     689      }
     690       
     691      OutEdgeIt f;
     692      for(g->first(f,v); g->valid(f); g->next(f)) {
     693        if ( 0 >= (*flow)[f] ) continue;
     694        Node u=g->head(f);
     695        if ( level[u] >= n ) {
     696          bfs_queue.push(u);
     697          level.set(u, l);
     698          if ( excess[u] > 0 ) active[l].push(u);
     699        }
     700      }
     701    }
     702    b=n-2;
     703
     704    while ( true ) {
     705       
     706      if ( b == 0 ) break;
     707
     708      if ( active[b].empty() ) --b;
     709      else {
     710        Node w=active[b].top();
     711        active[b].pop();
     712        int newlevel=push(w,active);     
     713
     714        //relabel
     715        if ( excess[w] > 0 ) {
     716          level.set(w,++newlevel);
     717          active[newlevel].push(w);
     718          b=newlevel;
     719        }
     720      }  // if stack[b] is nonempty
     721    } // while(true)
     722  }
    707723
    708724
    … …  
    711727  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
    712728  {
    713       ResGW res_graph(*g, *capacity, *flow);
    714       bool _augment=false;
    715      
    716       //ReachedMap level(res_graph);
    717       FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
    718       BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    719       bfs.pushAndSetReached(s);
    720        
    721       typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
    722       pred.set(s, INVALID);
    723      
    724       typename ResGW::template NodeMap<Num> free(res_graph);
    725        
    726       //searching for augmenting path
    727       while ( !bfs.finished() ) {
    728         ResGWOutEdgeIt e=bfs;
    729         if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
    730           Node v=res_graph.tail(e);
    731           Node w=res_graph.head(e);
    732           pred.set(w, e);
    733           if (res_graph.valid(pred[v])) {
    734             free.set(w, std::min(free[v], res_graph.resCap(e)));
    735           } else {
    736             free.set(w, res_graph.resCap(e));
    737           }
    738           if (res_graph.head(e)==t) { _augment=true; break; }
    739         }
    740        
    741         ++bfs;
    742       } //end of searching augmenting path
    743 
    744       if (_augment) {
    745         Node n=t;
    746         Num augment_value=free[t];
    747         while (res_graph.valid(pred[n])) {
    748           ResGWEdge e=pred[n];
    749           res_graph.augment(e, augment_value);
    750           n=res_graph.tail(e);
    751         }
    752       }
    753 
    754       return _augment;
    755     }
     729    ResGW res_graph(*g, *capacity, *flow);
     730    bool _augment=false;
     731     
     732    //ReachedMap level(res_graph);
     733    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
     734    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
     735    bfs.pushAndSetReached(s);
     736       
     737    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
     738    pred.set(s, INVALID);
     739     
     740    typename ResGW::template NodeMap<Num> free(res_graph);
     741       
     742    //searching for augmenting path
     743    while ( !bfs.finished() ) {
     744      ResGWOutEdgeIt e=bfs;
     745      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
     746        Node v=res_graph.tail(e);
     747        Node w=res_graph.head(e);
     748        pred.set(w, e);
     749        if (res_graph.valid(pred[v])) {
     750          free.set(w, std::min(free[v], res_graph.resCap(e)));
     751        } else {
     752          free.set(w, res_graph.resCap(e));
     753        }
     754        if (res_graph.head(e)==t) { _augment=true; break; }
     755      }
     756       
     757      ++bfs;
     758    } //end of searching augmenting path
     759
     760    if (_augment) {
     761      Node n=t;
     762      Num augment_value=free[t];
     763      while (res_graph.valid(pred[n])) {
     764        ResGWEdge e=pred[n];
     765        res_graph.augment(e, augment_value);
     766        n=res_graph.tail(e);
     767      }
     768    }
     769
     770    return _augment;
     771  }
    756772
    757773
    … …  
    767783  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
    768784  {     
    769       typedef MutableGraph MG;
    770       bool _augment=false;
    771 
    772       ResGW res_graph(*g, *capacity, *flow);
    773 
    774       //bfs for distances on the residual graph
    775       //ReachedMap level(res_graph);
    776       FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
    777       BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    778       bfs.pushAndSetReached(s);
    779       typename ResGW::template NodeMap<int>
    780         dist(res_graph); //filled up with 0's
    781 
    782       //F will contain the physical copy of the residual graph
    783       //with the set of edges which are on shortest paths
    784       MG F;
    785       typename ResGW::template NodeMap<typename MG::Node>
    786         res_graph_to_F(res_graph);
    787       {
    788         typename ResGW::NodeIt n;
    789         for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
    790           res_graph_to_F.set(n, F.addNode());
    791         }
    792       }
    793 
    794       typename MG::Node sF=res_graph_to_F[s];
    795       typename MG::Node tF=res_graph_to_F[t];
    796       typename MG::template EdgeMap<ResGWEdge> original_edge(F);
    797       typename MG::template EdgeMap<Num> residual_capacity(F);
    798 
    799       while ( !bfs.finished() ) {
    800         ResGWOutEdgeIt e=bfs;
    801         if (res_graph.valid(e)) {
    802           if (bfs.isBNodeNewlyReached()) {
    803             dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
     785    typedef MutableGraph MG;
     786    bool _augment=false;
     787
     788    ResGW res_graph(*g, *capacity, *flow);
     789
     790    //bfs for distances on the residual graph
     791    //ReachedMap level(res_graph);
     792    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
     793    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
     794    bfs.pushAndSetReached(s);
     795    typename ResGW::template NodeMap<int>
     796      dist(res_graph); //filled up with 0's
     797
     798    //F will contain the physical copy of the residual graph
     799    //with the set of edges which are on shortest paths
     800    MG F;
     801    typename ResGW::template NodeMap<typename MG::Node>
     802      res_graph_to_F(res_graph);
     803    {
     804      typename ResGW::NodeIt n;
     805      for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
     806        res_graph_to_F.set(n, F.addNode());
     807      }
     808    }
     809
     810    typename MG::Node sF=res_graph_to_F[s];
     811    typename MG::Node tF=res_graph_to_F[t];
     812    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
     813    typename MG::template EdgeMap<Num> residual_capacity(F);
     814
     815    while ( !bfs.finished() ) {
     816      ResGWOutEdgeIt e=bfs;
     817      if (res_graph.valid(e)) {
     818        if (bfs.isBNodeNewlyReached()) {
     819          dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
     820          typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
     821          original_edge.update();
     822          original_edge.set(f, e);
     823          residual_capacity.update();
     824          residual_capacity.set(f, res_graph.resCap(e));
     825        } else {
     826          if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
    804827            typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
    805828            original_edge.update();
    … …  
    807830            residual_capacity.update();
    808831            residual_capacity.set(f, res_graph.resCap(e));
     832          }
     833        }
     834      }
     835      ++bfs;
     836    } //computing distances from s in the residual graph
     837
     838    bool __augment=true;
     839
     840    while (__augment) {
     841      __augment=false;
     842      //computing blocking flow with dfs
     843      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
     844      typename MG::template NodeMap<typename MG::Edge> pred(F);
     845      pred.set(sF, INVALID);
     846      //invalid iterators for sources
     847
     848      typename MG::template NodeMap<Num> free(F);
     849
     850      dfs.pushAndSetReached(sF);     
     851      while (!dfs.finished()) {
     852        ++dfs;
     853        if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
     854          if (dfs.isBNodeNewlyReached()) {
     855            typename MG::Node v=F.aNode(dfs);
     856            typename MG::Node w=F.bNode(dfs);
     857            pred.set(w, dfs);
     858            if (F.valid(pred[v])) {
     859              free.set(w, std::min(free[v], residual_capacity[dfs]));
     860            } else {
     861              free.set(w, residual_capacity[dfs]);
     862            }
     863            if (w==tF) {
     864              __augment=true;
     865              _augment=true;
     866              break;
     867            }
     868             
    809869          } else {
    810             if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
    811               typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
    812               original_edge.update();
    813               original_edge.set(f, e);
    814               residual_capacity.update();
    815               residual_capacity.set(f, res_graph.resCap(e));
    816             }
     870            F.erase(/*typename MG::OutEdgeIt*/(dfs));
    817871          }
    818         }
    819         ++bfs;
    820       } //computing distances from s in the residual graph
    821 
    822       bool __augment=true;
    823 
    824       while (__augment) {
    825         __augment=false;
    826         //computing blocking flow with dfs
    827         DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
    828         typename MG::template NodeMap<typename MG::Edge> pred(F);
    829         pred.set(sF, INVALID);
    830         //invalid iterators for sources
    831 
    832         typename MG::template NodeMap<Num> free(F);
    833 
    834         dfs.pushAndSetReached(sF);     
    835         while (!dfs.finished()) {
    836           ++dfs;
    837           if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
    838             if (dfs.isBNodeNewlyReached()) {
    839               typename MG::Node v=F.aNode(dfs);
    840               typename MG::Node w=F.bNode(dfs);
    841               pred.set(w, dfs);
    842               if (F.valid(pred[v])) {
    843                 free.set(w, std::min(free[v], residual_capacity[dfs]));
    844               } else {
    845                 free.set(w, residual_capacity[dfs]);
    846               }
    847               if (w==tF) {
    848                 __augment=true;
    849                 _augment=true;
    850                 break;
    851               }
    852              
    853             } else {
    854               F.erase(/*typename MG::OutEdgeIt*/(dfs));
    855             }
    856           }
    857         }
    858 
    859         if (__augment) {
    860           typename MG::Node n=tF;
    861           Num augment_value=free[tF];
    862           while (F.valid(pred[n])) {
    863             typename MG::Edge e=pred[n];
    864             res_graph.augment(original_edge[e], augment_value);
    865             n=F.tail(e);
    866             if (residual_capacity[e]==augment_value)
    867               F.erase(e);
    868             else
    869               residual_capacity.set(e, residual_capacity[e]-augment_value);
    870           }
    871         }
    872        
    873       }
     872        }
     873      }
     874
     875      if (__augment) {
     876        typename MG::Node n=tF;
     877        Num augment_value=free[tF];
     878        while (F.valid(pred[n])) {
     879          typename MG::Edge e=pred[n];
     880          res_graph.augment(original_edge[e], augment_value);
     881          n=F.tail(e);
     882          if (residual_capacity[e]==augment_value)
     883            F.erase(e);
     884          else
     885            residual_capacity.set(e, residual_capacity[e]-augment_value);
     886        }
     887      }
     888       
     889    }
    874890           
    875       return _augment;
    876     }
     891    return _augment;
     892  }
    877893
    878894
    … …  
    884900  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
    885901  {
    886       bool _augment=false;
    887 
    888       ResGW res_graph(*g, *capacity, *flow);
    889      
    890       //ReachedMap level(res_graph);
    891       FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
    892       BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    893 
    894       bfs.pushAndSetReached(s);
    895       DistanceMap<ResGW> dist(res_graph);
    896       while ( !bfs.finished() ) {
    897         ResGWOutEdgeIt e=bfs;
    898         if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
    899           dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
    900         }
    901         ++bfs;
    902       } //computing distances from s in the residual graph
     902    bool _augment=false;
     903
     904    ResGW res_graph(*g, *capacity, *flow);
     905     
     906    //ReachedMap level(res_graph);
     907    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
     908    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
     909
     910    bfs.pushAndSetReached(s);
     911    DistanceMap<ResGW> dist(res_graph);
     912    while ( !bfs.finished() ) {
     913      ResGWOutEdgeIt e=bfs;
     914      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
     915        dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
     916      }
     917      ++bfs;
     918    } //computing distances from s in the residual graph
    903919
    904920      //Subgraph containing the edges on some shortest paths
    905       ConstMap<typename ResGW::Node, bool> true_map(true);
    906       typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
    907         DistanceMap<ResGW> > FilterResGW;
    908       FilterResGW filter_res_graph(res_graph, true_map, dist);
    909 
    910       //Subgraph, which is able to delete edges which are already
    911       //met by the dfs
    912       typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
    913         first_out_edges(filter_res_graph);
    914       typename FilterResGW::NodeIt v;
    915       for(filter_res_graph.first(v); filter_res_graph.valid(v);
    916           filter_res_graph.next(v))
     921    ConstMap<typename ResGW::Node, bool> true_map(true);
     922    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
     923      DistanceMap<ResGW> > FilterResGW;
     924    FilterResGW filter_res_graph(res_graph, true_map, dist);
     925
     926    //Subgraph, which is able to delete edges which are already
     927    //met by the dfs
     928    typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
     929      first_out_edges(filter_res_graph);
     930    typename FilterResGW::NodeIt v;
     931    for(filter_res_graph.first(v); filter_res_graph.valid(v);
     932        filter_res_graph.next(v))
    917933      {
    918934        typename FilterResGW::OutEdgeIt e;
    … …  
    920936        first_out_edges.set(v, e);
    921937      }
    922       typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
    923         template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
    924       ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
    925 
    926       bool __augment=true;
    927 
    928       while (__augment) {
    929 
    930         __augment=false;
    931         //computing blocking flow with dfs
    932         DfsIterator< ErasingResGW,
    933           typename ErasingResGW::template NodeMap<bool> >
    934           dfs(erasing_res_graph);
    935         typename ErasingResGW::
    936           template NodeMap<typename ErasingResGW::OutEdgeIt>
    937           pred(erasing_res_graph);
    938         pred.set(s, INVALID);
    939         //invalid iterators for sources
    940 
    941         typename ErasingResGW::template NodeMap<Num>
    942           free1(erasing_res_graph);
    943 
    944         dfs.pushAndSetReached(
    945           typename ErasingResGW::Node(
    946             typename FilterResGW::Node(
    947               typename ResGW::Node(s)
    948               )
    949             )
    950           );
    951         while (!dfs.finished()) {
    952           ++dfs;
    953           if (erasing_res_graph.valid(
    954                 typename ErasingResGW::OutEdgeIt(dfs)))
     938    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
     939      template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
     940    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
     941
     942    bool __augment=true;
     943
     944    while (__augment) {
     945
     946      __augment=false;
     947      //computing blocking flow with dfs
     948      DfsIterator< ErasingResGW,
     949        typename ErasingResGW::template NodeMap<bool> >
     950        dfs(erasing_res_graph);
     951      typename ErasingResGW::
     952        template NodeMap<typename ErasingResGW::OutEdgeIt>
     953        pred(erasing_res_graph);
     954      pred.set(s, INVALID);
     955      //invalid iterators for sources
     956
     957      typename ErasingResGW::template NodeMap<Num>
     958        free1(erasing_res_graph);
     959
     960      dfs.pushAndSetReached(
     961                            typename ErasingResGW::Node(
     962                                                        typename FilterResGW::Node(
     963                                                                                   typename ResGW::Node(s)
     964                                                                                   )
     965                                                        )
     966                            );
     967      while (!dfs.finished()) {
     968        ++dfs;
     969        if (erasing_res_graph.valid(
     970                                    typename ErasingResGW::OutEdgeIt(dfs)))
    955971          {
    956972            if (dfs.isBNodeNewlyReached()) {
    … …  
    962978              if (erasing_res_graph.valid(pred[v])) {
    963979                free1.set(w, std::min(free1[v], res_graph.resCap(
    964                                       typename ErasingResGW::OutEdgeIt(dfs))));
     980                                                                typename ErasingResGW::OutEdgeIt(dfs))));
    965981              } else {
    966982                free1.set(w, res_graph.resCap(
    967                            typename ErasingResGW::OutEdgeIt(dfs)));
     983                                              typename ErasingResGW::OutEdgeIt(dfs)));
    968984              }
    969985             
    … …  
    977993            }
    978994          }
    979         }       
    980 
    981         if (__augment) {
    982           typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
    983 //        typename ResGW::NodeMap<Num> a(res_graph);
    984 //        typename ResGW::Node b;
    985 //        Num j=a[b];
    986 //        typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
    987 //        typename FilterResGW::Node b1;
    988 //        Num j1=a1[b1];
    989 //        typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
    990 //        typename ErasingResGW::Node b2;
    991 //        Num j2=a2[b2];
    992           Num augment_value=free1[n];
    993           while (erasing_res_graph.valid(pred[n])) {
    994             typename ErasingResGW::OutEdgeIt e=pred[n];
    995             res_graph.augment(e, augment_value);
    996             n=erasing_res_graph.tail(e);
    997             if (res_graph.resCap(e)==0)
    998               erasing_res_graph.erase(e);
    999         }
    1000       }
    1001      
    1002       } //while (__augment)
     995      }
     996
     997      if (__augment) {
     998        typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
     999        //        typename ResGW::NodeMap<Num> a(res_graph);
     1000        //        typename ResGW::Node b;
     1001        //        Num j=a[b];
     1002        //        typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
     1003        //        typename FilterResGW::Node b1;
     1004        //        Num j1=a1[b1];
     1005        //        typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
     1006        //        typename ErasingResGW::Node b2;
     1007        //        Num j2=a2[b2];
     1008        Num augment_value=free1[n];
     1009        while (erasing_res_graph.valid(pred[n])) {
     1010          typename ErasingResGW::OutEdgeIt e=pred[n];
     1011          res_graph.augment(e, augment_value);
     1012          n=erasing_res_graph.tail(e);
     1013          if (res_graph.resCap(e)==0)
     1014            erasing_res_graph.erase(e);
     1015        }
     1016      }
     1017     
     1018    } //while (__augment)
    10031019           
    1004       return _augment;
    1005     }
     1020    return _augment;
     1021  }
    10061022
    10071023
Note: See TracChangeset for help on using the changeset viewer.