preflow.h

Go to the documentation of this file.
00001 /* -*- C++ -*-
00002  *
00003  * This file is a part of LEMON, a generic C++ optimization library
00004  *
00005  * Copyright (C) 2003-2006
00006  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
00007  * (Egervary Research Group on Combinatorial Optimization, EGRES).
00008  *
00009  * Permission to use, modify and distribute this software is granted
00010  * provided that this copyright notice appears in all copies. For
00011  * precise terms see the accompanying LICENSE file.
00012  *
00013  * This software is provided "AS IS" with no warranty of any kind,
00014  * express or implied, and with no claim as to its suitability for any
00015  * purpose.
00016  *
00017  */
00018 
00019 #ifndef LEMON_PREFLOW_H
00020 #define LEMON_PREFLOW_H
00021 
00022 #include <vector>
00023 #include <queue>
00024 
00025 #include <lemon/error.h>
00026 #include <lemon/invalid.h>
00027 #include <lemon/tolerance.h>
00028 #include <lemon/maps.h>
00029 #include <lemon/graph_utils.h>
00030 
00034 
00035 namespace lemon {
00036 
00065   template <typename Graph, typename Num,
00066             typename CapacityMap=typename Graph::template EdgeMap<Num>,
00067             typename FlowMap=typename Graph::template EdgeMap<Num>,
00068             typename TOL=Tolerance<Num> >
00069   class Preflow {
00070   protected:
00071     typedef typename Graph::Node Node;
00072     typedef typename Graph::NodeIt NodeIt;
00073     typedef typename Graph::EdgeIt EdgeIt;
00074     typedef typename Graph::OutEdgeIt OutEdgeIt;
00075     typedef typename Graph::InEdgeIt InEdgeIt;
00076 
00077     typedef typename Graph::template NodeMap<Node> NNMap;
00078     typedef typename std::vector<Node> VecNode;
00079 
00080     const Graph* _g;
00081     Node _source;
00082     Node _target;
00083     const CapacityMap* _capacity;
00084     FlowMap* _flow;
00085 
00086     TOL surely;
00087     
00088     int _node_num;      //the number of nodes of G
00089     
00090     typename Graph::template NodeMap<int> level;  
00091     typename Graph::template NodeMap<Num> excess;
00092 
00093     // constants used for heuristics
00094     static const int H0=20;
00095     static const int H1=1;
00096 
00097   public:
00098 
00100 
00102     class InvalidArgument : public lemon::LogicError {
00103     public:
00104       virtual const char* exceptionName() const {
00105         return "lemon::Preflow::InvalidArgument";
00106       }
00107     };
00108     
00109     
00111     
00114     enum FlowEnum{
00118       NO_FLOW,
00120       ZERO_FLOW,
00124       GEN_FLOW,
00127       PRE_FLOW
00128     };
00129 
00131 
00134     enum StatusEnum {
00137       AFTER_NOTHING,
00139       AFTER_PREFLOW_PHASE_1,      
00141       AFTER_PREFLOW_PHASE_2
00142     };
00143     
00144   protected: 
00145     FlowEnum flow_prop;
00146     StatusEnum status; // Do not needle this flag only if necessary.
00147     
00148   public: 
00150 
00161       Preflow(const Graph& _gr, Node _s, Node _t, 
00162               const CapacityMap& _cap, FlowMap& _f,
00163               const TOL &tol=TOL()) :
00164         _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
00165         _flow(&_f), surely(tol),
00166         _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), 
00167         flow_prop(NO_FLOW), status(AFTER_NOTHING) { 
00168         if ( _source==_target )
00169           throw InvalidArgument();
00170       }
00171     
00173 
00176     TOL &tolerance() { return surely; }
00177 
00179 
00182     void run() {
00183       phase1(flow_prop);
00184       phase2();
00185     }
00186     
00188     
00197     void run(FlowEnum fp) {
00198       flow_prop=fp;
00199       run();
00200     }
00201       
00203 
00217     void phase1(FlowEnum fp)
00218     {
00219       flow_prop=fp;
00220       phase1();
00221     }
00222 
00223     
00225 
00234     void phase1()
00235     {
00236       int heur0=(int)(H0*_node_num);  //time while running 'bound decrease'
00237       int heur1=(int)(H1*_node_num);  //time while running 'highest label'
00238       int heur=heur1;         //starting time interval (#of relabels)
00239       int numrelabel=0;
00240 
00241       bool what_heur=1;
00242       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
00243 
00244       bool end=false;
00245       //Needed for 'bound decrease', true means no active 
00246       //nodes are above bound b.
00247 
00248       int k=_node_num-2;  //bound on the highest level under n containing a node
00249       int b=k;    //bound on the highest level under n of an active node
00250 
00251       VecNode first(_node_num, INVALID);
00252       NNMap next(*_g, INVALID);
00253 
00254       NNMap left(*_g, INVALID);
00255       NNMap right(*_g, INVALID);
00256       VecNode level_list(_node_num,INVALID);
00257       //List of the nodes in level i<n, set to n.
00258 
00259       preflowPreproc(first, next, level_list, left, right);
00260 
00261       //Push/relabel on the highest level active nodes.
00262       while ( true ) {
00263         if ( b == 0 ) {
00264           if ( !what_heur && !end && k > 0 ) {
00265             b=k;
00266             end=true;
00267           } else break;
00268         }
00269 
00270         if ( first[b]==INVALID ) --b;
00271         else {
00272           end=false;
00273           Node w=first[b];
00274           first[b]=next[w];
00275           int newlevel=push(w, next, first);
00276           if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list, 
00277                                        left, right, b, k, what_heur);
00278 
00279           ++numrelabel;
00280           if ( numrelabel >= heur ) {
00281             numrelabel=0;
00282             if ( what_heur ) {
00283               what_heur=0;
00284               heur=heur0;
00285               end=false;
00286             } else {
00287               what_heur=1;
00288               heur=heur1;
00289               b=k;
00290             }
00291           }
00292         }
00293       }
00294       flow_prop=PRE_FLOW;
00295       status=AFTER_PREFLOW_PHASE_1;
00296     }
00297     // Heuristics:
00298     //   2 phase
00299     //   gap
00300     //   list 'level_list' on the nodes on level i implemented by hand
00301     //   stack 'active' on the active nodes on level i      
00302     //   runs heuristic 'highest label' for H1*n relabels
00303     //   runs heuristic 'bound decrease' for H0*n relabels,
00304     //        starts with 'highest label'
00305     //   Parameters H0 and H1 are initialized to 20 and 1.
00306 
00307 
00309 
00318     void phase2()
00319     {
00320 
00321       int k=_node_num-2;  //bound on the highest level under n containing a node
00322       int b=k;    //bound on the highest level under n of an active node
00323 
00324     
00325       VecNode first(_node_num, INVALID);
00326       NNMap next(*_g, INVALID); 
00327       level.set(_source,0);
00328       std::queue<Node> bfs_queue;
00329       bfs_queue.push(_source);
00330 
00331       while ( !bfs_queue.empty() ) {
00332 
00333         Node v=bfs_queue.front();
00334         bfs_queue.pop();
00335         int l=level[v]+1;
00336 
00337         for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
00338           if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
00339           Node u=_g->source(e);
00340           if ( level[u] >= _node_num ) {
00341             bfs_queue.push(u);
00342             level.set(u, l);
00343             if ( excess[u] > 0 ) {
00344               next.set(u,first[l]);
00345               first[l]=u;
00346             }
00347           }
00348         }
00349 
00350         for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
00351           if ( 0 >= (*_flow)[e] ) continue;
00352           Node u=_g->target(e);
00353           if ( level[u] >= _node_num ) {
00354             bfs_queue.push(u);
00355             level.set(u, l);
00356             if ( excess[u] > 0 ) {
00357               next.set(u,first[l]);
00358               first[l]=u;
00359             }
00360           }
00361         }
00362       }
00363       b=_node_num-2;
00364 
00365       while ( true ) {
00366 
00367         if ( b == 0 ) break;
00368         if ( first[b]==INVALID ) --b;
00369         else {
00370           Node w=first[b];
00371           first[b]=next[w];
00372           int newlevel=push(w,next, first);
00373           
00374           //relabel
00375           if ( excess[w] > 0 ) {
00376             level.set(w,++newlevel);
00377             next.set(w,first[newlevel]);
00378             first[newlevel]=w;
00379             b=newlevel;
00380           }
00381         } 
00382       } // while(true)
00383       flow_prop=GEN_FLOW;
00384       status=AFTER_PREFLOW_PHASE_2;
00385     }
00386 
00388 
00392     Num flowValue() const {
00393       return excess[_target];
00394     }
00395 
00396 
00398 
00405     template<typename _CutMap>
00406     void minCut(_CutMap& M) const {
00407       switch ( status ) {
00408         case AFTER_PREFLOW_PHASE_1:
00409         for(NodeIt v(*_g); v!=INVALID; ++v) {
00410           if (level[v] < _node_num) {
00411             M.set(v, false);
00412           } else {
00413             M.set(v, true);
00414           }
00415         }
00416         break;
00417         case AFTER_PREFLOW_PHASE_2:
00418         minMinCut(M);
00419         break;
00420         case AFTER_NOTHING:
00421         break;
00422       }
00423     }
00424 
00426 
00432     template<typename _CutMap>
00433     void minMinCut(_CutMap& M) const {
00434 
00435       std::queue<Node> queue;
00436       M.set(_source,true);
00437       queue.push(_source);
00438       
00439       while (!queue.empty()) {
00440         Node w=queue.front();
00441         queue.pop();
00442         
00443         for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00444           Node v=_g->target(e);
00445           if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
00446             queue.push(v);
00447             M.set(v, true);
00448           }
00449         }
00450         
00451         for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00452           Node v=_g->source(e);
00453           if (!M[v] && (*_flow)[e] > 0 ) {
00454             queue.push(v);
00455             M.set(v, true);
00456           }
00457         }
00458       }
00459     }
00460     
00462 
00467     template<typename _CutMap>
00468     void maxMinCut(_CutMap& M) const {
00469 
00470       for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
00471 
00472       std::queue<Node> queue;
00473 
00474       M.set(_target,false);
00475       queue.push(_target);
00476 
00477       while (!queue.empty()) {
00478         Node w=queue.front();
00479         queue.pop();
00480 
00481         for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00482           Node v=_g->source(e);
00483           if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
00484             queue.push(v);
00485             M.set(v, false);
00486           }
00487         }
00488 
00489         for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00490           Node v=_g->target(e);
00491           if (M[v] && (*_flow)[e] > 0 ) {
00492             queue.push(v);
00493             M.set(v, false);
00494           }
00495         }
00496       }
00497     }
00498 
00500 
00503     void source(Node _s) { 
00504       _source=_s; 
00505       if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
00506       status=AFTER_NOTHING; 
00507     }
00508 
00510 
00513     Node source() const { 
00514       return _source;
00515     }
00516 
00518 
00521     void target(Node _t) { 
00522       _target=_t; 
00523       if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
00524       status=AFTER_NOTHING; 
00525     }
00526 
00528 
00531     Node target() const { 
00532       return _target;
00533     }
00534 
00536 
00539     void capacityMap(const CapacityMap& _cap) { 
00540       _capacity=&_cap; 
00541       status=AFTER_NOTHING; 
00542     }
00544 
00547     const CapacityMap &capacityMap() const { 
00548       return *_capacity;
00549     }
00550 
00552 
00555     void flowMap(FlowMap& _f) { 
00556       _flow=&_f; 
00557       flow_prop=NO_FLOW;
00558       status=AFTER_NOTHING; 
00559     }
00560      
00562 
00565     const FlowMap &flowMap() const { 
00566       return *_flow;
00567     }
00568 
00569   private:
00570 
00571     int push(Node w, NNMap& next, VecNode& first) {
00572 
00573       int lev=level[w];
00574       Num exc=excess[w];
00575       int newlevel=_node_num;       //bound on the next level of w
00576 
00577       for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00578         if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
00579         Node v=_g->target(e);
00580 
00581         if( lev > level[v] ) { //Push is allowed now
00582           
00583           if ( excess[v]<=0 && v!=_target && v!=_source ) {
00584             next.set(v,first[level[v]]);
00585             first[level[v]]=v;
00586           }
00587 
00588           Num cap=(*_capacity)[e];
00589           Num flo=(*_flow)[e];
00590           Num remcap=cap-flo;
00591           
00592           if ( remcap >= exc ) { //A nonsaturating push.
00593             
00594             _flow->set(e, flo+exc);
00595             excess.set(v, excess[v]+exc);
00596             exc=0;
00597             break;
00598 
00599           } else { //A saturating push.
00600             _flow->set(e, cap);
00601             excess.set(v, excess[v]+remcap);
00602             exc-=remcap;
00603           }
00604         } else if ( newlevel > level[v] ) newlevel = level[v];
00605       } //for out edges wv
00606 
00607       if ( exc > 0 ) {
00608         for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00609           
00610           if( (*_flow)[e] <= 0 ) continue;
00611           Node v=_g->source(e);
00612 
00613           if( lev > level[v] ) { //Push is allowed now
00614 
00615             if ( excess[v]<=0 && v!=_target && v!=_source ) {
00616               next.set(v,first[level[v]]);
00617               first[level[v]]=v;
00618             }
00619 
00620             Num flo=(*_flow)[e];
00621 
00622             if ( flo >= exc ) { //A nonsaturating push.
00623 
00624               _flow->set(e, flo-exc);
00625               excess.set(v, excess[v]+exc);
00626               exc=0;
00627               break;
00628             } else {  //A saturating push.
00629 
00630               excess.set(v, excess[v]+flo);
00631               exc-=flo;
00632               _flow->set(e,0);
00633             }
00634           } else if ( newlevel > level[v] ) newlevel = level[v];
00635         } //for in edges vw
00636 
00637       } // if w still has excess after the out edge for cycle
00638 
00639       excess.set(w, exc);
00640       
00641       return newlevel;
00642     }
00643     
00644     
00645     
00646     void preflowPreproc(VecNode& first, NNMap& next, 
00647                         VecNode& level_list, NNMap& left, NNMap& right)
00648     {
00649       for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
00650       std::queue<Node> bfs_queue;
00651       
00652       if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
00653         //Reverse_bfs from t in the residual graph,
00654         //to find the starting level.
00655         level.set(_target,0);
00656         bfs_queue.push(_target);
00657         
00658         while ( !bfs_queue.empty() ) {
00659           
00660           Node v=bfs_queue.front();
00661           bfs_queue.pop();
00662           int l=level[v]+1;
00663           
00664           for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
00665             if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
00666             Node w=_g->source(e);
00667             if ( level[w] == _node_num && w != _source ) {
00668               bfs_queue.push(w);
00669               Node z=level_list[l];
00670               if ( z!=INVALID ) left.set(z,w);
00671               right.set(w,z);
00672               level_list[l]=w;
00673               level.set(w, l);
00674             }
00675           }
00676           
00677           for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
00678             if ( 0 >= (*_flow)[e] ) continue;
00679             Node w=_g->target(e);
00680             if ( level[w] == _node_num && w != _source ) {
00681               bfs_queue.push(w);
00682               Node z=level_list[l];
00683               if ( z!=INVALID ) left.set(z,w);
00684               right.set(w,z);
00685               level_list[l]=w;
00686               level.set(w, l);
00687             }
00688           }
00689         } //while
00690       } //if
00691 
00692 
00693       switch (flow_prop) {
00694         case NO_FLOW:  
00695         for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
00696         case ZERO_FLOW:
00697         for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
00698         
00699         //Reverse_bfs from t, to find the starting level.
00700         level.set(_target,0);
00701         bfs_queue.push(_target);
00702         
00703         while ( !bfs_queue.empty() ) {
00704           
00705           Node v=bfs_queue.front();
00706           bfs_queue.pop();
00707           int l=level[v]+1;
00708           
00709           for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
00710             Node w=_g->source(e);
00711             if ( level[w] == _node_num && w != _source ) {
00712               bfs_queue.push(w);
00713               Node z=level_list[l];
00714               if ( z!=INVALID ) left.set(z,w);
00715               right.set(w,z);
00716               level_list[l]=w;
00717               level.set(w, l);
00718             }
00719           }
00720         }
00721         
00722         //the starting flow
00723         for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
00724           Num c=(*_capacity)[e];
00725           if ( c <= 0 ) continue;
00726           Node w=_g->target(e);
00727           if ( level[w] < _node_num ) {
00728             if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
00729               next.set(w,first[level[w]]);
00730               first[level[w]]=w;
00731             }
00732             _flow->set(e, c);
00733             excess.set(w, excess[w]+c);
00734           }
00735         }
00736         break;
00737 
00738         case GEN_FLOW:
00739         for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
00740         {
00741           Num exc=0;
00742           for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
00743           for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
00744           excess.set(_target,exc);
00745         }
00746 
00747         //the starting flow
00748         for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)  {
00749           Num rem=(*_capacity)[e]-(*_flow)[e];
00750           if ( rem <= 0 ) continue;
00751           Node w=_g->target(e);
00752           if ( level[w] < _node_num ) {
00753             if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
00754               next.set(w,first[level[w]]);
00755               first[level[w]]=w;
00756             }   
00757             _flow->set(e, (*_capacity)[e]);
00758             excess.set(w, excess[w]+rem);
00759           }
00760         }
00761         
00762         for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
00763           if ( (*_flow)[e] <= 0 ) continue;
00764           Node w=_g->source(e);
00765           if ( level[w] < _node_num ) {
00766             if ( excess[w] <= 0 && w!=_target ) {
00767               next.set(w,first[level[w]]);
00768               first[level[w]]=w;
00769             }  
00770             excess.set(w, excess[w]+(*_flow)[e]);
00771             _flow->set(e, 0);
00772           }
00773         }
00774         break;
00775 
00776         case PRE_FLOW:  
00777         //the starting flow
00778         for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
00779           Num rem=(*_capacity)[e]-(*_flow)[e];
00780           if ( rem <= 0 ) continue;
00781           Node w=_g->target(e);
00782           if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
00783         }
00784         
00785         for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
00786           if ( (*_flow)[e] <= 0 ) continue;
00787           Node w=_g->source(e);
00788           if ( level[w] < _node_num ) _flow->set(e, 0);
00789         }
00790         
00791         //computing the excess
00792         for(NodeIt w(*_g); w!=INVALID; ++w) {
00793           Num exc=0;
00794           for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
00795           for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
00796           excess.set(w,exc);
00797           
00798           //putting the active nodes into the stack
00799           int lev=level[w];
00800             if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
00801               next.set(w,first[lev]);
00802               first[lev]=w;
00803             }
00804         }
00805         break;
00806       } //switch
00807     } //preflowPreproc
00808 
00809 
00810     void relabel(Node w, int newlevel, VecNode& first, NNMap& next, 
00811                  VecNode& level_list, NNMap& left,
00812                  NNMap& right, int& b, int& k, bool what_heur )
00813     {
00814 
00815       int lev=level[w];
00816 
00817       Node right_n=right[w];
00818       Node left_n=left[w];
00819 
00820       //unlacing starts
00821       if ( right_n!=INVALID ) {
00822         if ( left_n!=INVALID ) {
00823           right.set(left_n, right_n);
00824           left.set(right_n, left_n);
00825         } else {
00826           level_list[lev]=right_n;
00827           left.set(right_n, INVALID);
00828         }
00829       } else {
00830         if ( left_n!=INVALID ) {
00831           right.set(left_n, INVALID);
00832         } else {
00833           level_list[lev]=INVALID;
00834         }
00835       }
00836       //unlacing ends
00837 
00838       if ( level_list[lev]==INVALID ) {
00839 
00840         //gapping starts
00841         for (int i=lev; i!=k ; ) {
00842           Node v=level_list[++i];
00843           while ( v!=INVALID ) {
00844             level.set(v,_node_num);
00845             v=right[v];
00846           }
00847           level_list[i]=INVALID;
00848           if ( !what_heur ) first[i]=INVALID;
00849         }
00850 
00851         level.set(w,_node_num);
00852         b=lev-1;
00853         k=b;
00854         //gapping ends
00855 
00856       } else {
00857 
00858         if ( newlevel == _node_num ) level.set(w,_node_num);
00859         else {
00860           level.set(w,++newlevel);
00861           next.set(w,first[newlevel]);
00862           first[newlevel]=w;
00863           if ( what_heur ) b=newlevel;
00864           if ( k < newlevel ) ++k;      //now k=newlevel
00865           Node z=level_list[newlevel];
00866           if ( z!=INVALID ) left.set(z,w);
00867           right.set(w,z);
00868           left.set(w,INVALID);
00869           level_list[newlevel]=w;
00870         }
00871       }
00872     } //relabel
00873 
00874   }; 
00875 
00881   template<class GR, class CM, class FM>
00882   Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
00883                             typename GR::Node source,
00884                             typename GR::Node target,
00885                             const CM &cap,
00886                             FM &flow
00887                             )
00888   {
00889     return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
00890   }
00891 
00892 } //namespace lemon
00893 
00894 #endif //LEMON_PREFLOW_H

Generated on Fri Feb 3 18:39:40 2006 for LEMON by  doxygen 1.4.6