lemon/preflow.h
author deba
Tue, 24 Jan 2006 16:07:38 +0000
changeset 1901 723b2b81d900
parent 1875 98698b69a902
child 1953 d4f411003580
permissions -rw-r--r--
Lemon Graph Format uses label instead of id named map.
     1 /* -*- C++ -*-
     2  * lemon/preflow.h - Part of LEMON, a generic C++ optimization library
     3  *
     4  * Copyright (C) 2006 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     5  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     6  *
     7  * Permission to use, modify and distribute this software is granted
     8  * provided that this copyright notice appears in all copies. For
     9  * precise terms see the accompanying LICENSE file.
    10  *
    11  * This software is provided "AS IS" with no warranty of any kind,
    12  * express or implied, and with no claim as to its suitability for any
    13  * purpose.
    14  *
    15  */
    16 
    17 #ifndef LEMON_PREFLOW_H
    18 #define LEMON_PREFLOW_H
    19 
    20 #include <vector>
    21 #include <queue>
    22 
    23 #include <lemon/error.h>
    24 #include <lemon/invalid.h>
    25 #include <lemon/tolerance.h>
    26 #include <lemon/maps.h>
    27 #include <lemon/graph_utils.h>
    28 
    29 /// \file
    30 /// \ingroup flowalgs
    31 /// \brief Implementation of the preflow algorithm.
    32 
    33 namespace lemon {
    34 
    35   ///\ingroup flowalgs
    36   ///\brief %Preflow algorithms class.
    37   ///
    38   ///This class provides an implementation of the \e preflow \e
    39   ///algorithm producing a flow of maximum value in a directed
    40   ///graph. The preflow algorithms are the fastest known max flow algorithms
    41   ///up to now. The \e source node, the \e target node, the \e
    42   ///capacity of the edges and the \e starting \e flow value of the
    43   ///edges should be passed to the algorithm through the
    44   ///constructor. It is possible to change these quantities using the
    45   ///functions \ref source, \ref target, \ref capacityMap and \ref
    46   ///flowMap.
    47   ///
    48   ///After running \ref lemon::Preflow::phase1() "phase1()"
    49   ///or \ref lemon::Preflow::run() "run()", the maximal flow
    50   ///value can be obtained by calling \ref flowValue(). The minimum
    51   ///value cut can be written into a <tt>bool</tt> node map by
    52   ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
    53   ///the inclusionwise minimum and maximum of the minimum value cuts,
    54   ///resp.)
    55   ///
    56   ///\param Graph The directed graph type the algorithm runs on.
    57   ///\param Num The number type of the capacities and the flow values.
    58   ///\param CapacityMap The capacity map type.
    59   ///\param FlowMap The flow map type.
    60   ///
    61   ///\author Jacint Szabo 
    62   ///\todo Second template parameter is superfluous
    63   template <typename Graph, typename Num,
    64 	    typename CapacityMap=typename Graph::template EdgeMap<Num>,
    65             typename FlowMap=typename Graph::template EdgeMap<Num>,
    66 	    typename TOL=Tolerance<Num> >
    67   class Preflow {
    68   protected:
    69     typedef typename Graph::Node Node;
    70     typedef typename Graph::NodeIt NodeIt;
    71     typedef typename Graph::EdgeIt EdgeIt;
    72     typedef typename Graph::OutEdgeIt OutEdgeIt;
    73     typedef typename Graph::InEdgeIt InEdgeIt;
    74 
    75     typedef typename Graph::template NodeMap<Node> NNMap;
    76     typedef typename std::vector<Node> VecNode;
    77 
    78     const Graph* _g;
    79     Node _source;
    80     Node _target;
    81     const CapacityMap* _capacity;
    82     FlowMap* _flow;
    83 
    84     TOL surely;
    85     
    86     int _node_num;      //the number of nodes of G
    87     
    88     typename Graph::template NodeMap<int> level;  
    89     typename Graph::template NodeMap<Num> excess;
    90 
    91     // constants used for heuristics
    92     static const int H0=20;
    93     static const int H1=1;
    94 
    95   public:
    96 
    97     ///\ref Exception for the case when s=t.
    98 
    99     ///\ref Exception for the case when the source equals the target.
   100     class InvalidArgument : public lemon::LogicError {
   101     public:
   102       virtual const char* exceptionName() const {
   103 	return "lemon::Preflow::InvalidArgument";
   104       }
   105     };
   106     
   107     
   108     ///Indicates the property of the starting flow map.
   109     
   110     ///Indicates the property of the starting flow map.
   111     ///
   112     enum FlowEnum{
   113       ///indicates an unspecified edge map. \c flow will be 
   114       ///set to the constant zero flow in the beginning of
   115       ///the algorithm in this case.
   116       NO_FLOW,
   117       ///constant zero flow
   118       ZERO_FLOW,
   119       ///any flow, i.e. the sum of the in-flows equals to
   120       ///the sum of the out-flows in every node except the \c source and
   121       ///the \c target.
   122       GEN_FLOW,
   123       ///any preflow, i.e. the sum of the in-flows is at 
   124       ///least the sum of the out-flows in every node except the \c source.
   125       PRE_FLOW
   126     };
   127 
   128     ///Indicates the state of the preflow algorithm.
   129 
   130     ///Indicates the state of the preflow algorithm.
   131     ///
   132     enum StatusEnum {
   133       ///before running the algorithm or
   134       ///at an unspecified state.
   135       AFTER_NOTHING,
   136       ///right after running \ref phase1()
   137       AFTER_PREFLOW_PHASE_1,      
   138       ///after running \ref phase2()
   139       AFTER_PREFLOW_PHASE_2
   140     };
   141     
   142   protected: 
   143     FlowEnum flow_prop;
   144     StatusEnum status; // Do not needle this flag only if necessary.
   145     
   146   public: 
   147     ///The constructor of the class.
   148 
   149     ///The constructor of the class. 
   150     ///\param _gr The directed graph the algorithm runs on. 
   151     ///\param _s The source node.
   152     ///\param _t The target node.
   153     ///\param _cap The capacity of the edges. 
   154     ///\param _f The flow of the edges. 
   155     ///Except the graph, all of these parameters can be reset by
   156     ///calling \ref source, \ref target, \ref capacityMap and \ref
   157     ///flowMap, resp.
   158       Preflow(const Graph& _gr, Node _s, Node _t, 
   159 	      const CapacityMap& _cap, FlowMap& _f,
   160 	      const TOL &tol=TOL()) :
   161 	_g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
   162 	_flow(&_f), surely(tol),
   163 	_node_num(countNodes(_gr)), level(_gr), excess(_gr,0), 
   164 	flow_prop(NO_FLOW), status(AFTER_NOTHING) { 
   165 	if ( _source==_target )
   166 	  throw InvalidArgument();
   167       }
   168     
   169     ///Give a reference to the tolerance handler class
   170 
   171     ///Give a reference to the tolerance handler class
   172     ///\sa Tolerance
   173     TOL &tolerance() { return surely; }
   174 
   175     ///Runs the preflow algorithm.  
   176 
   177     ///Runs the preflow algorithm.
   178     ///
   179     void run() {
   180       phase1(flow_prop);
   181       phase2();
   182     }
   183     
   184     ///Runs the preflow algorithm.  
   185     
   186     ///Runs the preflow algorithm. 
   187     ///\pre The starting flow map must be
   188     /// - a constant zero flow if \c fp is \c ZERO_FLOW,
   189     /// - an arbitrary flow if \c fp is \c GEN_FLOW,
   190     /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
   191     /// - any map if \c fp is NO_FLOW.
   192     ///If the starting flow map is a flow or a preflow then 
   193     ///the algorithm terminates faster.
   194     void run(FlowEnum fp) {
   195       flow_prop=fp;
   196       run();
   197     }
   198       
   199     ///Runs the first phase of the preflow algorithm.
   200 
   201     ///The preflow algorithm consists of two phases, this method runs
   202     ///the first phase. After the first phase the maximum flow value
   203     ///and a minimum value cut can already be computed, although a
   204     ///maximum flow is not yet obtained. So after calling this method
   205     ///\ref flowValue returns the value of a maximum flow and \ref
   206     ///minCut returns a minimum cut.     
   207     ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
   208     ///value cuts unless calling \ref phase2.  
   209     ///\pre The starting flow must be 
   210     ///- a constant zero flow if \c fp is \c ZERO_FLOW, 
   211     ///- an arbitary flow if \c fp is \c GEN_FLOW, 
   212     ///- an arbitary preflow if \c fp is \c PRE_FLOW, 
   213     ///- any map if \c fp is NO_FLOW.
   214     void phase1(FlowEnum fp)
   215     {
   216       flow_prop=fp;
   217       phase1();
   218     }
   219 
   220     
   221     ///Runs the first phase of the preflow algorithm.
   222 
   223     ///The preflow algorithm consists of two phases, this method runs
   224     ///the first phase. After the first phase the maximum flow value
   225     ///and a minimum value cut can already be computed, although a
   226     ///maximum flow is not yet obtained. So after calling this method
   227     ///\ref flowValue returns the value of a maximum flow and \ref
   228     ///minCut returns a minimum cut.
   229     ///\warning \ref minMinCut() and \ref maxMinCut() do not
   230     ///give minimum value cuts unless calling \ref phase2().
   231     void phase1()
   232     {
   233       int heur0=(int)(H0*_node_num);  //time while running 'bound decrease'
   234       int heur1=(int)(H1*_node_num);  //time while running 'highest label'
   235       int heur=heur1;         //starting time interval (#of relabels)
   236       int numrelabel=0;
   237 
   238       bool what_heur=1;
   239       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   240 
   241       bool end=false;
   242       //Needed for 'bound decrease', true means no active 
   243       //nodes are above bound b.
   244 
   245       int k=_node_num-2;  //bound on the highest level under n containing a node
   246       int b=k;    //bound on the highest level under n of an active node
   247 
   248       VecNode first(_node_num, INVALID);
   249       NNMap next(*_g, INVALID);
   250 
   251       NNMap left(*_g, INVALID);
   252       NNMap right(*_g, INVALID);
   253       VecNode level_list(_node_num,INVALID);
   254       //List of the nodes in level i<n, set to n.
   255 
   256       preflowPreproc(first, next, level_list, left, right);
   257 
   258       //Push/relabel on the highest level active nodes.
   259       while ( true ) {
   260 	if ( b == 0 ) {
   261 	  if ( !what_heur && !end && k > 0 ) {
   262 	    b=k;
   263 	    end=true;
   264 	  } else break;
   265 	}
   266 
   267 	if ( first[b]==INVALID ) --b;
   268 	else {
   269 	  end=false;
   270 	  Node w=first[b];
   271 	  first[b]=next[w];
   272 	  int newlevel=push(w, next, first);
   273 	  if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list, 
   274 				       left, right, b, k, what_heur);
   275 
   276 	  ++numrelabel;
   277 	  if ( numrelabel >= heur ) {
   278 	    numrelabel=0;
   279 	    if ( what_heur ) {
   280 	      what_heur=0;
   281 	      heur=heur0;
   282 	      end=false;
   283 	    } else {
   284 	      what_heur=1;
   285 	      heur=heur1;
   286 	      b=k;
   287 	    }
   288 	  }
   289 	}
   290       }
   291       flow_prop=PRE_FLOW;
   292       status=AFTER_PREFLOW_PHASE_1;
   293     }
   294     // Heuristics:
   295     //   2 phase
   296     //   gap
   297     //   list 'level_list' on the nodes on level i implemented by hand
   298     //   stack 'active' on the active nodes on level i      
   299     //   runs heuristic 'highest label' for H1*n relabels
   300     //   runs heuristic 'bound decrease' for H0*n relabels,
   301     //        starts with 'highest label'
   302     //   Parameters H0 and H1 are initialized to 20 and 1.
   303 
   304 
   305     ///Runs the second phase of the preflow algorithm.
   306 
   307     ///The preflow algorithm consists of two phases, this method runs
   308     ///the second phase. After calling \ref phase1() and then
   309     ///\ref phase2(),
   310     /// \ref flowMap() return a maximum flow, \ref flowValue
   311     ///returns the value of a maximum flow, \ref minCut returns a
   312     ///minimum cut, while the methods \ref minMinCut and \ref
   313     ///maxMinCut return the inclusionwise minimum and maximum cuts of
   314     ///minimum value, resp.  \pre \ref phase1 must be called before.
   315     void phase2()
   316     {
   317 
   318       int k=_node_num-2;  //bound on the highest level under n containing a node
   319       int b=k;    //bound on the highest level under n of an active node
   320 
   321     
   322       VecNode first(_node_num, INVALID);
   323       NNMap next(*_g, INVALID); 
   324       level.set(_source,0);
   325       std::queue<Node> bfs_queue;
   326       bfs_queue.push(_source);
   327 
   328       while ( !bfs_queue.empty() ) {
   329 
   330 	Node v=bfs_queue.front();
   331 	bfs_queue.pop();
   332 	int l=level[v]+1;
   333 
   334 	for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
   335 	  if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
   336 	  Node u=_g->source(e);
   337 	  if ( level[u] >= _node_num ) {
   338 	    bfs_queue.push(u);
   339 	    level.set(u, l);
   340 	    if ( excess[u] > 0 ) {
   341 	      next.set(u,first[l]);
   342 	      first[l]=u;
   343 	    }
   344 	  }
   345 	}
   346 
   347 	for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
   348 	  if ( 0 >= (*_flow)[e] ) continue;
   349 	  Node u=_g->target(e);
   350 	  if ( level[u] >= _node_num ) {
   351 	    bfs_queue.push(u);
   352 	    level.set(u, l);
   353 	    if ( excess[u] > 0 ) {
   354 	      next.set(u,first[l]);
   355 	      first[l]=u;
   356 	    }
   357 	  }
   358 	}
   359       }
   360       b=_node_num-2;
   361 
   362       while ( true ) {
   363 
   364 	if ( b == 0 ) break;
   365 	if ( first[b]==INVALID ) --b;
   366 	else {
   367 	  Node w=first[b];
   368 	  first[b]=next[w];
   369 	  int newlevel=push(w,next, first);
   370 	  
   371 	  //relabel
   372 	  if ( excess[w] > 0 ) {
   373 	    level.set(w,++newlevel);
   374 	    next.set(w,first[newlevel]);
   375 	    first[newlevel]=w;
   376 	    b=newlevel;
   377 	  }
   378 	} 
   379       } // while(true)
   380       flow_prop=GEN_FLOW;
   381       status=AFTER_PREFLOW_PHASE_2;
   382     }
   383 
   384     /// Returns the value of the maximum flow.
   385 
   386     /// Returns the value of the maximum flow by returning the excess
   387     /// of the target node \c t. This value equals to the value of
   388     /// the maximum flow already after running \ref phase1.
   389     Num flowValue() const {
   390       return excess[_target];
   391     }
   392 
   393 
   394     ///Returns a minimum value cut.
   395 
   396     ///Sets \c M to the characteristic vector of a minimum value
   397     ///cut. This method can be called both after running \ref
   398     ///phase1 and \ref phase2. It is much faster after
   399     ///\ref phase1.  \pre M should be a bool-valued node-map. \pre
   400     ///If \ref minCut() is called after \ref phase2() then M should
   401     ///be initialized to false.
   402     template<typename _CutMap>
   403     void minCut(_CutMap& M) const {
   404       switch ( status ) {
   405 	case AFTER_PREFLOW_PHASE_1:
   406 	for(NodeIt v(*_g); v!=INVALID; ++v) {
   407 	  if (level[v] < _node_num) {
   408 	    M.set(v, false);
   409 	  } else {
   410 	    M.set(v, true);
   411 	  }
   412 	}
   413 	break;
   414 	case AFTER_PREFLOW_PHASE_2:
   415 	minMinCut(M);
   416 	break;
   417 	case AFTER_NOTHING:
   418 	break;
   419       }
   420     }
   421 
   422     ///Returns the inclusionwise minimum of the minimum value cuts.
   423 
   424     ///Sets \c M to the characteristic vector of the minimum value cut
   425     ///which is inclusionwise minimum. It is computed by processing a
   426     ///bfs from the source node \c s in the residual graph.  \pre M
   427     ///should be a node map of bools initialized to false.  \pre \ref
   428     ///phase2 should already be run.
   429     template<typename _CutMap>
   430     void minMinCut(_CutMap& M) const {
   431 
   432       std::queue<Node> queue;
   433       M.set(_source,true);
   434       queue.push(_source);
   435       
   436       while (!queue.empty()) {
   437 	Node w=queue.front();
   438 	queue.pop();
   439 	
   440 	for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   441 	  Node v=_g->target(e);
   442 	  if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
   443 	    queue.push(v);
   444 	    M.set(v, true);
   445 	  }
   446 	}
   447 	
   448 	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   449 	  Node v=_g->source(e);
   450 	  if (!M[v] && (*_flow)[e] > 0 ) {
   451 	    queue.push(v);
   452 	    M.set(v, true);
   453 	  }
   454 	}
   455       }
   456     }
   457     
   458     ///Returns the inclusionwise maximum of the minimum value cuts.
   459 
   460     ///Sets \c M to the characteristic vector of the minimum value cut
   461     ///which is inclusionwise maximum. It is computed by processing a
   462     ///backward bfs from the target node \c t in the residual graph.
   463     ///\pre \ref phase2() or run() should already be run.
   464     template<typename _CutMap>
   465     void maxMinCut(_CutMap& M) const {
   466 
   467       for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
   468 
   469       std::queue<Node> queue;
   470 
   471       M.set(_target,false);
   472       queue.push(_target);
   473 
   474       while (!queue.empty()) {
   475         Node w=queue.front();
   476 	queue.pop();
   477 
   478 	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   479 	  Node v=_g->source(e);
   480 	  if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
   481 	    queue.push(v);
   482 	    M.set(v, false);
   483 	  }
   484 	}
   485 
   486 	for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   487 	  Node v=_g->target(e);
   488 	  if (M[v] && (*_flow)[e] > 0 ) {
   489 	    queue.push(v);
   490 	    M.set(v, false);
   491 	  }
   492 	}
   493       }
   494     }
   495 
   496     ///Sets the source node to \c _s.
   497 
   498     ///Sets the source node to \c _s.
   499     /// 
   500     void source(Node _s) { 
   501       _source=_s; 
   502       if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
   503       status=AFTER_NOTHING; 
   504     }
   505 
   506     ///Returns the source node.
   507 
   508     ///Returns the source node.
   509     /// 
   510     Node source() const { 
   511       return _source;
   512     }
   513 
   514     ///Sets the target node to \c _t.
   515 
   516     ///Sets the target node to \c _t.
   517     ///
   518     void target(Node _t) { 
   519       _target=_t; 
   520       if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
   521       status=AFTER_NOTHING; 
   522     }
   523 
   524     ///Returns the target node.
   525 
   526     ///Returns the target node.
   527     /// 
   528     Node target() const { 
   529       return _target;
   530     }
   531 
   532     /// Sets the edge map of the capacities to _cap.
   533 
   534     /// Sets the edge map of the capacities to _cap.
   535     /// 
   536     void capacityMap(const CapacityMap& _cap) { 
   537       _capacity=&_cap; 
   538       status=AFTER_NOTHING; 
   539     }
   540     /// Returns a reference to capacity map.
   541 
   542     /// Returns a reference to capacity map.
   543     /// 
   544     const CapacityMap &capacityMap() const { 
   545       return *_capacity;
   546     }
   547 
   548     /// Sets the edge map of the flows to _flow.
   549 
   550     /// Sets the edge map of the flows to _flow.
   551     /// 
   552     void flowMap(FlowMap& _f) { 
   553       _flow=&_f; 
   554       flow_prop=NO_FLOW;
   555       status=AFTER_NOTHING; 
   556     }
   557      
   558     /// Returns a reference to flow map.
   559 
   560     /// Returns a reference to flow map.
   561     /// 
   562     const FlowMap &flowMap() const { 
   563       return *_flow;
   564     }
   565 
   566   private:
   567 
   568     int push(Node w, NNMap& next, VecNode& first) {
   569 
   570       int lev=level[w];
   571       Num exc=excess[w];
   572       int newlevel=_node_num;       //bound on the next level of w
   573 
   574       for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   575 	if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
   576 	Node v=_g->target(e);
   577 
   578 	if( lev > level[v] ) { //Push is allowed now
   579 	  
   580 	  if ( excess[v]<=0 && v!=_target && v!=_source ) {
   581 	    next.set(v,first[level[v]]);
   582 	    first[level[v]]=v;
   583 	  }
   584 
   585 	  Num cap=(*_capacity)[e];
   586 	  Num flo=(*_flow)[e];
   587 	  Num remcap=cap-flo;
   588 	  
   589 	  if ( remcap >= exc ) { //A nonsaturating push.
   590 	    
   591 	    _flow->set(e, flo+exc);
   592 	    excess.set(v, excess[v]+exc);
   593 	    exc=0;
   594 	    break;
   595 
   596 	  } else { //A saturating push.
   597 	    _flow->set(e, cap);
   598 	    excess.set(v, excess[v]+remcap);
   599 	    exc-=remcap;
   600 	  }
   601 	} else if ( newlevel > level[v] ) newlevel = level[v];
   602       } //for out edges wv
   603 
   604       if ( exc > 0 ) {
   605 	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   606 	  
   607 	  if( (*_flow)[e] <= 0 ) continue;
   608 	  Node v=_g->source(e);
   609 
   610 	  if( lev > level[v] ) { //Push is allowed now
   611 
   612 	    if ( excess[v]<=0 && v!=_target && v!=_source ) {
   613 	      next.set(v,first[level[v]]);
   614 	      first[level[v]]=v;
   615 	    }
   616 
   617 	    Num flo=(*_flow)[e];
   618 
   619 	    if ( flo >= exc ) { //A nonsaturating push.
   620 
   621 	      _flow->set(e, flo-exc);
   622 	      excess.set(v, excess[v]+exc);
   623 	      exc=0;
   624 	      break;
   625 	    } else {  //A saturating push.
   626 
   627 	      excess.set(v, excess[v]+flo);
   628 	      exc-=flo;
   629 	      _flow->set(e,0);
   630 	    }
   631 	  } else if ( newlevel > level[v] ) newlevel = level[v];
   632 	} //for in edges vw
   633 
   634       } // if w still has excess after the out edge for cycle
   635 
   636       excess.set(w, exc);
   637       
   638       return newlevel;
   639     }
   640     
   641     
   642     
   643     void preflowPreproc(VecNode& first, NNMap& next, 
   644 			VecNode& level_list, NNMap& left, NNMap& right)
   645     {
   646       for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
   647       std::queue<Node> bfs_queue;
   648       
   649       if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
   650 	//Reverse_bfs from t in the residual graph,
   651 	//to find the starting level.
   652 	level.set(_target,0);
   653 	bfs_queue.push(_target);
   654 	
   655 	while ( !bfs_queue.empty() ) {
   656 	  
   657 	  Node v=bfs_queue.front();
   658 	  bfs_queue.pop();
   659 	  int l=level[v]+1;
   660 	  
   661 	  for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
   662 	    if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
   663 	    Node w=_g->source(e);
   664 	    if ( level[w] == _node_num && w != _source ) {
   665 	      bfs_queue.push(w);
   666 	      Node z=level_list[l];
   667 	      if ( z!=INVALID ) left.set(z,w);
   668 	      right.set(w,z);
   669 	      level_list[l]=w;
   670 	      level.set(w, l);
   671 	    }
   672 	  }
   673 	  
   674 	  for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
   675 	    if ( 0 >= (*_flow)[e] ) continue;
   676 	    Node w=_g->target(e);
   677 	    if ( level[w] == _node_num && w != _source ) {
   678 	      bfs_queue.push(w);
   679 	      Node z=level_list[l];
   680 	      if ( z!=INVALID ) left.set(z,w);
   681 	      right.set(w,z);
   682 	      level_list[l]=w;
   683 	      level.set(w, l);
   684 	    }
   685 	  }
   686 	} //while
   687       } //if
   688 
   689 
   690       switch (flow_prop) {
   691 	case NO_FLOW:  
   692 	for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
   693 	case ZERO_FLOW:
   694 	for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
   695 	
   696 	//Reverse_bfs from t, to find the starting level.
   697 	level.set(_target,0);
   698 	bfs_queue.push(_target);
   699 	
   700 	while ( !bfs_queue.empty() ) {
   701 	  
   702 	  Node v=bfs_queue.front();
   703 	  bfs_queue.pop();
   704 	  int l=level[v]+1;
   705 	  
   706 	  for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
   707 	    Node w=_g->source(e);
   708 	    if ( level[w] == _node_num && w != _source ) {
   709 	      bfs_queue.push(w);
   710 	      Node z=level_list[l];
   711 	      if ( z!=INVALID ) left.set(z,w);
   712 	      right.set(w,z);
   713 	      level_list[l]=w;
   714 	      level.set(w, l);
   715 	    }
   716 	  }
   717 	}
   718 	
   719 	//the starting flow
   720 	for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
   721 	  Num c=(*_capacity)[e];
   722 	  if ( c <= 0 ) continue;
   723 	  Node w=_g->target(e);
   724 	  if ( level[w] < _node_num ) {
   725 	    if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
   726 	      next.set(w,first[level[w]]);
   727 	      first[level[w]]=w;
   728 	    }
   729 	    _flow->set(e, c);
   730 	    excess.set(w, excess[w]+c);
   731 	  }
   732 	}
   733 	break;
   734 
   735 	case GEN_FLOW:
   736 	for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
   737 	{
   738 	  Num exc=0;
   739 	  for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
   740 	  for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
   741 	  excess.set(_target,exc);
   742 	}
   743 
   744 	//the starting flow
   745 	for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)	{
   746 	  Num rem=(*_capacity)[e]-(*_flow)[e];
   747 	  if ( rem <= 0 ) continue;
   748 	  Node w=_g->target(e);
   749 	  if ( level[w] < _node_num ) {
   750 	    if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
   751 	      next.set(w,first[level[w]]);
   752 	      first[level[w]]=w;
   753 	    }   
   754 	    _flow->set(e, (*_capacity)[e]);
   755 	    excess.set(w, excess[w]+rem);
   756 	  }
   757 	}
   758 	
   759 	for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
   760 	  if ( (*_flow)[e] <= 0 ) continue;
   761 	  Node w=_g->source(e);
   762 	  if ( level[w] < _node_num ) {
   763 	    if ( excess[w] <= 0 && w!=_target ) {
   764 	      next.set(w,first[level[w]]);
   765 	      first[level[w]]=w;
   766 	    }  
   767 	    excess.set(w, excess[w]+(*_flow)[e]);
   768 	    _flow->set(e, 0);
   769 	  }
   770 	}
   771 	break;
   772 
   773 	case PRE_FLOW:	
   774 	//the starting flow
   775 	for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
   776 	  Num rem=(*_capacity)[e]-(*_flow)[e];
   777 	  if ( rem <= 0 ) continue;
   778 	  Node w=_g->target(e);
   779 	  if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
   780 	}
   781 	
   782 	for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
   783 	  if ( (*_flow)[e] <= 0 ) continue;
   784 	  Node w=_g->source(e);
   785 	  if ( level[w] < _node_num ) _flow->set(e, 0);
   786 	}
   787 	
   788 	//computing the excess
   789 	for(NodeIt w(*_g); w!=INVALID; ++w) {
   790 	  Num exc=0;
   791 	  for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
   792 	  for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
   793 	  excess.set(w,exc);
   794 	  
   795 	  //putting the active nodes into the stack
   796 	  int lev=level[w];
   797 	    if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
   798 	      next.set(w,first[lev]);
   799 	      first[lev]=w;
   800 	    }
   801 	}
   802 	break;
   803       } //switch
   804     } //preflowPreproc
   805 
   806 
   807     void relabel(Node w, int newlevel, VecNode& first, NNMap& next, 
   808 		 VecNode& level_list, NNMap& left,
   809 		 NNMap& right, int& b, int& k, bool what_heur )
   810     {
   811 
   812       int lev=level[w];
   813 
   814       Node right_n=right[w];
   815       Node left_n=left[w];
   816 
   817       //unlacing starts
   818       if ( right_n!=INVALID ) {
   819 	if ( left_n!=INVALID ) {
   820 	  right.set(left_n, right_n);
   821 	  left.set(right_n, left_n);
   822 	} else {
   823 	  level_list[lev]=right_n;
   824 	  left.set(right_n, INVALID);
   825 	}
   826       } else {
   827 	if ( left_n!=INVALID ) {
   828 	  right.set(left_n, INVALID);
   829 	} else {
   830 	  level_list[lev]=INVALID;
   831 	}
   832       }
   833       //unlacing ends
   834 
   835       if ( level_list[lev]==INVALID ) {
   836 
   837 	//gapping starts
   838 	for (int i=lev; i!=k ; ) {
   839 	  Node v=level_list[++i];
   840 	  while ( v!=INVALID ) {
   841 	    level.set(v,_node_num);
   842 	    v=right[v];
   843 	  }
   844 	  level_list[i]=INVALID;
   845 	  if ( !what_heur ) first[i]=INVALID;
   846 	}
   847 
   848 	level.set(w,_node_num);
   849 	b=lev-1;
   850 	k=b;
   851 	//gapping ends
   852 
   853       } else {
   854 
   855 	if ( newlevel == _node_num ) level.set(w,_node_num);
   856 	else {
   857 	  level.set(w,++newlevel);
   858 	  next.set(w,first[newlevel]);
   859 	  first[newlevel]=w;
   860 	  if ( what_heur ) b=newlevel;
   861 	  if ( k < newlevel ) ++k;      //now k=newlevel
   862 	  Node z=level_list[newlevel];
   863 	  if ( z!=INVALID ) left.set(z,w);
   864 	  right.set(w,z);
   865 	  left.set(w,INVALID);
   866 	  level_list[newlevel]=w;
   867 	}
   868       }
   869     } //relabel
   870 
   871   }; 
   872 
   873   ///\ingroup flowalgs
   874   ///\brief Function type interface for Preflow algorithm.
   875   ///
   876   ///Function type interface for Preflow algorithm.
   877   ///\sa Preflow
   878   template<class GR, class CM, class FM>
   879   Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
   880 			    typename GR::Node source,
   881 			    typename GR::Node target,
   882 			    const CM &cap,
   883 			    FM &flow
   884 			    )
   885   {
   886     return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
   887   }
   888 
   889 } //namespace lemon
   890 
   891 #endif //LEMON_PREFLOW_H