Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

preflow.h

Go to the documentation of this file.
00001 /* -*- C++ -*-
00002  * src/lemon/preflow.h - Part of LEMON, a generic C++ optimization library
00003  *
00004  * Copyright (C) 2004 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
00005  * (Egervary Combinatorial Optimization Research Group, EGRES).
00006  *
00007  * Permission to use, modify and distribute this software is granted
00008  * provided that this copyright notice appears in all copies. For
00009  * precise terms see the accompanying LICENSE file.
00010  *
00011  * This software is provided "AS IS" with no warranty of any kind,
00012  * express or implied, and with no claim as to its suitability for any
00013  * purpose.
00014  *
00015  */
00016 
00017 #ifndef LEMON_PREFLOW_H
00018 #define LEMON_PREFLOW_H
00019 
00020 #include <vector>
00021 #include <queue>
00022 
00023 #include <lemon/invalid.h>
00024 #include <lemon/maps.h>
00025 #include <lemon/graph_utils.h>
00026 
00030 
00031 namespace lemon {
00032 
00035 
00037 
00062   template <typename Graph, typename Num,
00063             typename CapMap=typename Graph::template EdgeMap<Num>,
00064             typename FlowMap=typename Graph::template EdgeMap<Num> >
00065   class Preflow {
00066   protected:
00067     typedef typename Graph::Node Node;
00068     typedef typename Graph::NodeIt NodeIt;
00069     typedef typename Graph::EdgeIt EdgeIt;
00070     typedef typename Graph::OutEdgeIt OutEdgeIt;
00071     typedef typename Graph::InEdgeIt InEdgeIt;
00072 
00073     typedef typename Graph::template NodeMap<Node> NNMap;
00074     typedef typename std::vector<Node> VecNode;
00075 
00076     const Graph* g;
00077     Node s;
00078     Node t;
00079     const CapMap* capacity;
00080     FlowMap* flow;
00081     int n;      //the number of nodes of G
00082     
00083     typename Graph::template NodeMap<int> level;  
00084     typename Graph::template NodeMap<Num> excess;
00085 
00086     // constants used for heuristics
00087     static const int H0=20;
00088     static const int H1=1;
00089 
00090     public:
00091 
00093 
00105     enum FlowEnum{
00106       NO_FLOW,
00107       ZERO_FLOW,
00108       GEN_FLOW,
00109       PRE_FLOW
00110     };
00111 
00113 
00119     enum StatusEnum {
00120       AFTER_NOTHING,
00121       AFTER_PREFLOW_PHASE_1,      
00122       AFTER_PREFLOW_PHASE_2
00123     };
00124     
00125     protected: 
00126       FlowEnum flow_prop;
00127     StatusEnum status; // Do not needle this flag only if necessary.
00128     
00129   public: 
00131 
00141       Preflow(const Graph& _G, Node _s, Node _t, 
00142               const CapMap& _capacity, FlowMap& _flow) :
00143         g(&_G), s(_s), t(_t), capacity(&_capacity),
00144         flow(&_flow), n(countNodes(_G)), level(_G), excess(_G,0), 
00145         flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
00146 
00147 
00148                                                                               
00150 
00153     void run() {
00154       phase1(flow_prop);
00155       phase2();
00156     }
00157     
00159     
00168     void run(FlowEnum fp) {
00169       flow_prop=fp;
00170       run();
00171     }
00172       
00174 
00188     void phase1(FlowEnum fp)
00189     {
00190       flow_prop=fp;
00191       phase1();
00192     }
00193 
00194     
00196 
00205     void phase1()
00206     {
00207       int heur0=(int)(H0*n);  //time while running 'bound decrease'
00208       int heur1=(int)(H1*n);  //time while running 'highest label'
00209       int heur=heur1;         //starting time interval (#of relabels)
00210       int numrelabel=0;
00211 
00212       bool what_heur=1;
00213       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
00214 
00215       bool end=false;
00216       //Needed for 'bound decrease', true means no active 
00217       //nodes are above bound b.
00218 
00219       int k=n-2;  //bound on the highest level under n containing a node
00220       int b=k;    //bound on the highest level under n of an active node
00221 
00222       VecNode first(n, INVALID);
00223       NNMap next(*g, INVALID);
00224 
00225       NNMap left(*g, INVALID);
00226       NNMap right(*g, INVALID);
00227       VecNode level_list(n,INVALID);
00228       //List of the nodes in level i<n, set to n.
00229 
00230       preflowPreproc(first, next, level_list, left, right);
00231 
00232       //Push/relabel on the highest level active nodes.
00233       while ( true ) {
00234         if ( b == 0 ) {
00235           if ( !what_heur && !end && k > 0 ) {
00236             b=k;
00237             end=true;
00238           } else break;
00239         }
00240 
00241         if ( first[b]==INVALID ) --b;
00242         else {
00243           end=false;
00244           Node w=first[b];
00245           first[b]=next[w];
00246           int newlevel=push(w, next, first);
00247           if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list, 
00248                                        left, right, b, k, what_heur);
00249 
00250           ++numrelabel;
00251           if ( numrelabel >= heur ) {
00252             numrelabel=0;
00253             if ( what_heur ) {
00254               what_heur=0;
00255               heur=heur0;
00256               end=false;
00257             } else {
00258               what_heur=1;
00259               heur=heur1;
00260               b=k;
00261             }
00262           }
00263         }
00264       }
00265       flow_prop=PRE_FLOW;
00266       status=AFTER_PREFLOW_PHASE_1;
00267     }
00268     // Heuristics:
00269     //   2 phase
00270     //   gap
00271     //   list 'level_list' on the nodes on level i implemented by hand
00272     //   stack 'active' on the active nodes on level i      
00273     //   runs heuristic 'highest label' for H1*n relabels
00274     //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
00275     //   Parameters H0 and H1 are initialized to 20 and 1.
00276 
00277 
00279 
00287     void phase2()
00288     {
00289 
00290       int k=n-2;  //bound on the highest level under n containing a node
00291       int b=k;    //bound on the highest level under n of an active node
00292 
00293     
00294       VecNode first(n, INVALID);
00295       NNMap next(*g, INVALID); 
00296       level.set(s,0);
00297       std::queue<Node> bfs_queue;
00298       bfs_queue.push(s);
00299 
00300       while ( !bfs_queue.empty() ) {
00301 
00302         Node v=bfs_queue.front();
00303         bfs_queue.pop();
00304         int l=level[v]+1;
00305 
00306         for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
00307           if ( (*capacity)[e] <= (*flow)[e] ) continue;
00308           Node u=g->source(e);
00309           if ( level[u] >= n ) {
00310             bfs_queue.push(u);
00311             level.set(u, l);
00312             if ( excess[u] > 0 ) {
00313               next.set(u,first[l]);
00314               first[l]=u;
00315             }
00316           }
00317         }
00318 
00319         for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
00320           if ( 0 >= (*flow)[e] ) continue;
00321           Node u=g->target(e);
00322           if ( level[u] >= n ) {
00323             bfs_queue.push(u);
00324             level.set(u, l);
00325             if ( excess[u] > 0 ) {
00326               next.set(u,first[l]);
00327               first[l]=u;
00328             }
00329           }
00330         }
00331       }
00332       b=n-2;
00333 
00334       while ( true ) {
00335 
00336         if ( b == 0 ) break;
00337         if ( first[b]==INVALID ) --b;
00338         else {
00339           Node w=first[b];
00340           first[b]=next[w];
00341           int newlevel=push(w,next, first);
00342           
00343           //relabel
00344           if ( excess[w] > 0 ) {
00345             level.set(w,++newlevel);
00346             next.set(w,first[newlevel]);
00347             first[newlevel]=w;
00348             b=newlevel;
00349           }
00350         } 
00351       } // while(true)
00352       flow_prop=GEN_FLOW;
00353       status=AFTER_PREFLOW_PHASE_2;
00354     }
00355 
00357 
00361     Num flowValue() const {
00362       return excess[t];
00363     }
00364 
00365 
00367 
00374     template<typename _CutMap>
00375     void minCut(_CutMap& M) const {
00376       switch ( status ) {
00377         case AFTER_PREFLOW_PHASE_1:
00378         for(NodeIt v(*g); v!=INVALID; ++v) {
00379           if (level[v] < n) {
00380             M.set(v, false);
00381           } else {
00382             M.set(v, true);
00383           }
00384         }
00385         break;
00386         case AFTER_PREFLOW_PHASE_2:
00387         minMinCut(M);
00388         break;
00389         case AFTER_NOTHING:
00390         break;
00391       }
00392     }
00393 
00395 
00401     template<typename _CutMap>
00402     void minMinCut(_CutMap& M) const {
00403 
00404       std::queue<Node> queue;
00405       M.set(s,true);
00406       queue.push(s);
00407       
00408       while (!queue.empty()) {
00409         Node w=queue.front();
00410         queue.pop();
00411         
00412         for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00413           Node v=g->target(e);
00414           if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
00415             queue.push(v);
00416             M.set(v, true);
00417           }
00418         }
00419         
00420         for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00421           Node v=g->source(e);
00422           if (!M[v] && (*flow)[e] > 0 ) {
00423             queue.push(v);
00424             M.set(v, true);
00425           }
00426         }
00427       }
00428     }
00429     
00431 
00436     template<typename _CutMap>
00437     void maxMinCut(_CutMap& M) const {
00438 
00439       for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
00440 
00441       std::queue<Node> queue;
00442 
00443       M.set(t,false);
00444       queue.push(t);
00445 
00446       while (!queue.empty()) {
00447         Node w=queue.front();
00448         queue.pop();
00449 
00450         for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00451           Node v=g->source(e);
00452           if (M[v] && (*flow)[e] < (*capacity)[e] ) {
00453             queue.push(v);
00454             M.set(v, false);
00455           }
00456         }
00457 
00458         for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00459           Node v=g->target(e);
00460           if (M[v] && (*flow)[e] > 0 ) {
00461             queue.push(v);
00462             M.set(v, false);
00463           }
00464         }
00465       }
00466     }
00467 
00469 
00472     void setSource(Node _s) { 
00473       s=_s; 
00474       if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
00475       status=AFTER_NOTHING; 
00476     }
00477 
00479 
00482     void setTarget(Node _t) { 
00483       t=_t; 
00484       if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
00485       status=AFTER_NOTHING; 
00486     }
00487 
00489 
00492     void setCap(const CapMap& _cap) { 
00493       capacity=&_cap; 
00494       status=AFTER_NOTHING; 
00495     }
00496 
00498 
00501     void setFlow(FlowMap& _flow) { 
00502       flow=&_flow; 
00503       flow_prop=NO_FLOW;
00504       status=AFTER_NOTHING; 
00505     }
00506 
00507 
00508   private:
00509 
00510     int push(Node w, NNMap& next, VecNode& first) {
00511 
00512       int lev=level[w];
00513       Num exc=excess[w];
00514       int newlevel=n;       //bound on the next level of w
00515 
00516       for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00517         if ( (*flow)[e] >= (*capacity)[e] ) continue;
00518         Node v=g->target(e);
00519 
00520         if( lev > level[v] ) { //Push is allowed now
00521           
00522           if ( excess[v]<=0 && v!=t && v!=s ) {
00523             next.set(v,first[level[v]]);
00524             first[level[v]]=v;
00525           }
00526 
00527           Num cap=(*capacity)[e];
00528           Num flo=(*flow)[e];
00529           Num remcap=cap-flo;
00530           
00531           if ( remcap >= exc ) { //A nonsaturating push.
00532             
00533             flow->set(e, flo+exc);
00534             excess.set(v, excess[v]+exc);
00535             exc=0;
00536             break;
00537 
00538           } else { //A saturating push.
00539             flow->set(e, cap);
00540             excess.set(v, excess[v]+remcap);
00541             exc-=remcap;
00542           }
00543         } else if ( newlevel > level[v] ) newlevel = level[v];
00544       } //for out edges wv
00545 
00546       if ( exc > 0 ) {
00547         for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00548           
00549           if( (*flow)[e] <= 0 ) continue;
00550           Node v=g->source(e);
00551 
00552           if( lev > level[v] ) { //Push is allowed now
00553 
00554             if ( excess[v]<=0 && v!=t && v!=s ) {
00555               next.set(v,first[level[v]]);
00556               first[level[v]]=v;
00557             }
00558 
00559             Num flo=(*flow)[e];
00560 
00561             if ( flo >= exc ) { //A nonsaturating push.
00562 
00563               flow->set(e, flo-exc);
00564               excess.set(v, excess[v]+exc);
00565               exc=0;
00566               break;
00567             } else {  //A saturating push.
00568 
00569               excess.set(v, excess[v]+flo);
00570               exc-=flo;
00571               flow->set(e,0);
00572             }
00573           } else if ( newlevel > level[v] ) newlevel = level[v];
00574         } //for in edges vw
00575 
00576       } // if w still has excess after the out edge for cycle
00577 
00578       excess.set(w, exc);
00579       
00580       return newlevel;
00581     }
00582     
00583     
00584     
00585     void preflowPreproc(VecNode& first, NNMap& next, 
00586                         VecNode& level_list, NNMap& left, NNMap& right)
00587     {
00588       for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
00589       std::queue<Node> bfs_queue;
00590       
00591       if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
00592         //Reverse_bfs from t in the residual graph,
00593         //to find the starting level.
00594         level.set(t,0);
00595         bfs_queue.push(t);
00596         
00597         while ( !bfs_queue.empty() ) {
00598           
00599           Node v=bfs_queue.front();
00600           bfs_queue.pop();
00601           int l=level[v]+1;
00602           
00603           for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
00604             if ( (*capacity)[e] <= (*flow)[e] ) continue;
00605             Node w=g->source(e);
00606             if ( level[w] == n && w != s ) {
00607               bfs_queue.push(w);
00608               Node z=level_list[l];
00609               if ( z!=INVALID ) left.set(z,w);
00610               right.set(w,z);
00611               level_list[l]=w;
00612               level.set(w, l);
00613             }
00614           }
00615           
00616           for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
00617             if ( 0 >= (*flow)[e] ) continue;
00618             Node w=g->target(e);
00619             if ( level[w] == n && w != s ) {
00620               bfs_queue.push(w);
00621               Node z=level_list[l];
00622               if ( z!=INVALID ) left.set(z,w);
00623               right.set(w,z);
00624               level_list[l]=w;
00625               level.set(w, l);
00626             }
00627           }
00628         } //while
00629       } //if
00630 
00631 
00632       switch (flow_prop) {
00633         case NO_FLOW:  
00634         for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
00635         case ZERO_FLOW:
00636         for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
00637         
00638         //Reverse_bfs from t, to find the starting level.
00639         level.set(t,0);
00640         bfs_queue.push(t);
00641         
00642         while ( !bfs_queue.empty() ) {
00643           
00644           Node v=bfs_queue.front();
00645           bfs_queue.pop();
00646           int l=level[v]+1;
00647           
00648           for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
00649             Node w=g->source(e);
00650             if ( level[w] == n && w != s ) {
00651               bfs_queue.push(w);
00652               Node z=level_list[l];
00653               if ( z!=INVALID ) left.set(z,w);
00654               right.set(w,z);
00655               level_list[l]=w;
00656               level.set(w, l);
00657             }
00658           }
00659         }
00660         
00661         //the starting flow
00662         for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
00663           Num c=(*capacity)[e];
00664           if ( c <= 0 ) continue;
00665           Node w=g->target(e);
00666           if ( level[w] < n ) {
00667             if ( excess[w] <= 0 && w!=t ) { //putting into the stack
00668               next.set(w,first[level[w]]);
00669               first[level[w]]=w;
00670             }
00671             flow->set(e, c);
00672             excess.set(w, excess[w]+c);
00673           }
00674         }
00675         break;
00676 
00677         case GEN_FLOW:
00678         for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
00679         {
00680           Num exc=0;
00681           for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
00682           for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
00683           excess.set(t,exc);
00684         }
00685 
00686         //the starting flow
00687         for(OutEdgeIt e(*g,s); e!=INVALID; ++e) {
00688           Num rem=(*capacity)[e]-(*flow)[e];
00689           if ( rem <= 0 ) continue;
00690           Node w=g->target(e);
00691           if ( level[w] < n ) {
00692             if ( excess[w] <= 0 && w!=t ) { //putting into the stack
00693               next.set(w,first[level[w]]);
00694               first[level[w]]=w;
00695             }   
00696             flow->set(e, (*capacity)[e]);
00697             excess.set(w, excess[w]+rem);
00698           }
00699         }
00700         
00701         for(InEdgeIt e(*g,s); e!=INVALID; ++e) {
00702           if ( (*flow)[e] <= 0 ) continue;
00703           Node w=g->source(e);
00704           if ( level[w] < n ) {
00705             if ( excess[w] <= 0 && w!=t ) {
00706               next.set(w,first[level[w]]);
00707               first[level[w]]=w;
00708             }  
00709             excess.set(w, excess[w]+(*flow)[e]);
00710             flow->set(e, 0);
00711           }
00712         }
00713         break;
00714 
00715         case PRE_FLOW:  
00716         //the starting flow
00717         for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
00718           Num rem=(*capacity)[e]-(*flow)[e];
00719           if ( rem <= 0 ) continue;
00720           Node w=g->target(e);
00721           if ( level[w] < n ) flow->set(e, (*capacity)[e]);
00722         }
00723         
00724         for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
00725           if ( (*flow)[e] <= 0 ) continue;
00726           Node w=g->source(e);
00727           if ( level[w] < n ) flow->set(e, 0);
00728         }
00729         
00730         //computing the excess
00731         for(NodeIt w(*g); w!=INVALID; ++w) {
00732           Num exc=0;
00733           for(InEdgeIt e(*g,w); e!=INVALID; ++e) exc+=(*flow)[e];
00734           for(OutEdgeIt e(*g,w); e!=INVALID; ++e) exc-=(*flow)[e];
00735           excess.set(w,exc);
00736           
00737           //putting the active nodes into the stack
00738           int lev=level[w];
00739             if ( exc > 0 && lev < n && Node(w) != t ) {
00740               next.set(w,first[lev]);
00741               first[lev]=w;
00742             }
00743         }
00744         break;
00745       } //switch
00746     } //preflowPreproc
00747 
00748 
00749     void relabel(Node w, int newlevel, VecNode& first, NNMap& next, 
00750                  VecNode& level_list, NNMap& left,
00751                  NNMap& right, int& b, int& k, bool what_heur )
00752     {
00753 
00754       int lev=level[w];
00755 
00756       Node right_n=right[w];
00757       Node left_n=left[w];
00758 
00759       //unlacing starts
00760       if ( right_n!=INVALID ) {
00761         if ( left_n!=INVALID ) {
00762           right.set(left_n, right_n);
00763           left.set(right_n, left_n);
00764         } else {
00765           level_list[lev]=right_n;
00766           left.set(right_n, INVALID);
00767         }
00768       } else {
00769         if ( left_n!=INVALID ) {
00770           right.set(left_n, INVALID);
00771         } else {
00772           level_list[lev]=INVALID;
00773         }
00774       }
00775       //unlacing ends
00776 
00777       if ( level_list[lev]==INVALID ) {
00778 
00779         //gapping starts
00780         for (int i=lev; i!=k ; ) {
00781           Node v=level_list[++i];
00782           while ( v!=INVALID ) {
00783             level.set(v,n);
00784             v=right[v];
00785           }
00786           level_list[i]=INVALID;
00787           if ( !what_heur ) first[i]=INVALID;
00788         }
00789 
00790         level.set(w,n);
00791         b=lev-1;
00792         k=b;
00793         //gapping ends
00794 
00795       } else {
00796 
00797         if ( newlevel == n ) level.set(w,n);
00798         else {
00799           level.set(w,++newlevel);
00800           next.set(w,first[newlevel]);
00801           first[newlevel]=w;
00802           if ( what_heur ) b=newlevel;
00803           if ( k < newlevel ) ++k;      //now k=newlevel
00804           Node z=level_list[newlevel];
00805           if ( z!=INVALID ) left.set(z,w);
00806           right.set(w,z);
00807           left.set(w,INVALID);
00808           level_list[newlevel]=w;
00809         }
00810       }
00811     } //relabel
00812 
00813   }; 
00814 } //namespace lemon
00815 
00816 #endif //LEMON_PREFLOW_H
00817 
00818 
00819 
00820 

Generated on Mon Feb 21 15:02:22 2005 for LEMON by  doxygen 1.4.1