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