src/work/jacint/max_flow.h
author marci
Fri, 30 Apr 2004 17:10:01 +0000
changeset 499 767f3da8ce0e
parent 487 11ad69691d18
child 510 72143568cadc
permissions -rw-r--r--
A bipartite graph template can be used as BipartiteGraph<ListGraph>.
     1 // -*- C++ -*-
     2 
     3 /*
     4   Heuristics: 
     5   2 phase
     6   gap
     7   list 'level_list' on the nodes on level i implemented by hand
     8   stack 'active' on the active nodes on level i
     9   runs heuristic 'highest label' for H1*n relabels
    10   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
    11  
    12   Parameters H0 and H1 are initialized to 20 and 1.
    13 
    14   Constructors:
    15 
    16   Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if 
    17   FlowMap is not constant zero, and should be true if it is
    18 
    19   Members:
    20 
    21   void run()
    22 
    23   Num flowValue() : returns the value of a maximum flow
    24 
    25   void minMinCut(CutMap& M) : sets M to the characteristic vector of the 
    26   minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
    27 
    28   void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 
    29   maximum min cut. M should be a map of bools initialized to false.
    30 
    31   void minCut(CutMap& M) : sets M to the characteristic vector of 
    32   a min cut. M should be a map of bools initialized to false.
    33 
    34 */
    35 
    36 #ifndef HUGO_MAX_FLOW_H
    37 #define HUGO_MAX_FLOW_H
    38 
    39 #define H0 20
    40 #define H1 1
    41 
    42 #include <vector>
    43 #include <queue>
    44 #include <stack>
    45 
    46 #include <graph_wrapper.h>
    47 #include <bfs_iterator.h>
    48 #include <invalid.h>
    49 #include <maps.h>
    50 #include <for_each_macros.h>
    51 
    52 /// \file
    53 /// \brief Dimacs file format reader.
    54 
    55 namespace hugo {
    56 
    57 
    58 //  ///\author Marton Makai, Jacint Szabo
    59   /// A class for computing max flows and related quantities.
    60   template <typename Graph, typename Num, 
    61 	    typename CapMap=typename Graph::template EdgeMap<Num>, 
    62             typename FlowMap=typename Graph::template EdgeMap<Num> >
    63   class MaxFlow {
    64     
    65     typedef typename Graph::Node Node;
    66     typedef typename Graph::NodeIt NodeIt;
    67     typedef typename Graph::OutEdgeIt OutEdgeIt;
    68     typedef typename Graph::InEdgeIt InEdgeIt;
    69 
    70     typedef typename std::vector<std::stack<Node> > VecStack;
    71     typedef typename Graph::template NodeMap<Node> NNMap;
    72     typedef typename std::vector<Node> VecNode;
    73 
    74     const Graph* g;
    75     Node s;
    76     Node t;
    77     const CapMap* capacity;  
    78     FlowMap* flow;
    79     int n;      //the number of nodes of G
    80     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    81     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    82     typedef typename ResGW::Edge ResGWEdge;
    83     //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    84     typedef typename Graph::template NodeMap<int> ReachedMap;
    85     ReachedMap level;
    86     //level works as a bool map in augmenting path algorithms 
    87     //and is used by bfs for storing reached information.
    88     //In preflow, it shows levels of nodes.
    89     //typename Graph::template NodeMap<int> level;    
    90     typename Graph::template NodeMap<Num> excess; 
    91 
    92   public:
    93  
    94     ///\todo Document this
    95     enum flowEnum{
    96       ZERO_FLOW=0,
    97       GEN_FLOW=1,
    98       PREFLOW=2
    99     };
   100 
   101     MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, 
   102 	    FlowMap& _flow) :
   103       g(&_G), s(_s), t(_t), capacity(&_capacity), 
   104       flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
   105 
   106     /// A max flow algorithm is run.
   107     ///\pre the flow have to be 0 at the beginning.
   108     void run() {
   109       preflow(ZERO_FLOW);
   110     }
   111     
   112     /// A preflow algorithm is run. 
   113     ///\pre The initial edge-map have to be a 
   114     /// zero flow if \c fe is \c ZERO_FLOW,
   115     /// a flow if \c fe is \c GEN_FLOW, 
   116     /// and a pre-flow it is \c PREFLOW.
   117     void preflow(flowEnum fe) {
   118       preflowPhase0(fe);
   119       preflowPhase1();
   120     }
   121 
   122     /// Run the first phase of preflow, starting from a 0 flow, from a flow, 
   123     /// or from a preflow, according to \c fe.
   124     void preflowPhase0( flowEnum fe );
   125 
   126     /// Second phase of preflow.
   127     void preflowPhase1();
   128 
   129     /// Starting from a flow, this method searches for an augmenting path 
   130     /// according to the Edmonds-Karp algorithm 
   131     /// and augments the flow on if any. 
   132     /// The return value shows if the augmentation was succesful.
   133     bool augmentOnShortestPath();
   134 
   135     /// Starting from a flow, this method searches for an augmenting blockin 
   136     /// flow according to Dinits' algorithm and augments the flow on if any. 
   137     /// The blocking flow is computed in a physically constructed 
   138     /// residual graph of type \c Mutablegraph.
   139     /// The return value show sif the augmentation was succesful.
   140     template<typename MutableGraph> bool augmentOnBlockingFlow();
   141 
   142     /// The same as \c augmentOnBlockingFlow<MutableGraph> but the 
   143     /// residual graph is not constructed physically.
   144     /// The return value shows if the augmentation was succesful.
   145     bool augmentOnBlockingFlow2();
   146 
   147     /// Returns the actual flow value.
   148     /// More precisely, it returns the negative excess of s, thus 
   149     /// this works also for preflows.
   150     Num flowValue() { 
   151       Num a=0;
   152       FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
   153       FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
   154       return a;
   155     }
   156 
   157     /// Should be used between preflowPhase0 and preflowPhase1.
   158     ///\todo We have to make some status variable which shows the actual state 
   159     /// of the class. This enables us to determine which methods are valid 
   160     /// for MinCut computation
   161     template<typename _CutMap>
   162     void actMinCut(_CutMap& M) {
   163       NodeIt v;
   164       for(g->first(v); g->valid(v); g->next(v)) {
   165 	if ( level[v] < n ) {
   166 	  M.set(v,false);
   167 	} else {
   168 	  M.set(v,true);
   169 	}
   170       }
   171     }
   172 
   173     /// The unique inclusionwise minimum cut is computed by 
   174     /// processing a bfs from s in the residual graph.
   175     ///\pre flow have to be a max flow otherwise it will the whole node-set.
   176     template<typename _CutMap>
   177     void minMinCut(_CutMap& M) {
   178     
   179       std::queue<Node> queue;
   180       
   181       M.set(s,true);      
   182       queue.push(s);
   183 
   184       while (!queue.empty()) {
   185         Node w=queue.front();
   186 	queue.pop();
   187 
   188 	OutEdgeIt e;
   189 	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   190 	  Node v=g->head(e);
   191 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   192 	    queue.push(v);
   193 	    M.set(v, true);
   194 	  }
   195 	} 
   196 
   197 	InEdgeIt f;
   198 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   199 	  Node v=g->tail(f);
   200 	  if (!M[v] && (*flow)[f] > 0 ) {
   201 	    queue.push(v);
   202 	    M.set(v, true);
   203 	  }
   204 	} 
   205       }
   206     }
   207 
   208 
   209     /// The unique inclusionwise maximum cut is computed by 
   210     /// processing a reverse bfs from t in the residual graph.
   211     ///\pre flow have to be a max flow otherwise it will be empty.
   212     template<typename _CutMap>
   213     void maxMinCut(_CutMap& M) {
   214 
   215       NodeIt v;
   216       for(g->first(v) ; g->valid(v); g->next(v)) {
   217 	M.set(v, true);
   218       }
   219 
   220       std::queue<Node> queue;
   221       
   222       M.set(t,false);        
   223       queue.push(t);
   224 
   225       while (!queue.empty()) {
   226         Node w=queue.front();
   227 	queue.pop();
   228 
   229 
   230 	InEdgeIt e;
   231 	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   232 	  Node v=g->tail(e);
   233 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   234 	    queue.push(v);
   235 	    M.set(v, false);
   236 	  }
   237 	}
   238 	
   239 	OutEdgeIt f;
   240 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   241 	  Node v=g->head(f);
   242 	  if (M[v] && (*flow)[f] > 0 ) {
   243 	    queue.push(v);
   244 	    M.set(v, false);
   245 	  }
   246 	}
   247       }
   248     }
   249 
   250 
   251     /// A minimum cut is computed.
   252     template<typename CutMap>
   253     void minCut(CutMap& M) { minMinCut(M); }
   254 
   255     ///
   256     void resetSource(Node _s) { s=_s; }
   257     ///
   258     void resetTarget(Node _t) { t=_t; }
   259    
   260     /// capacity-map is changed.
   261     void resetCap(const CapMap& _cap) { capacity=&_cap; }
   262     
   263     /// flow-map is changed. 
   264     void resetFlow(FlowMap& _flow) { flow=&_flow; }
   265 
   266 
   267   private:
   268 
   269     int push(Node w, VecStack& active) {
   270       
   271       int lev=level[w];
   272       Num exc=excess[w];
   273       int newlevel=n;       //bound on the next level of w
   274 	  
   275       OutEdgeIt e;
   276       for(g->first(e,w); g->valid(e); g->next(e)) {
   277 	    
   278 	if ( (*flow)[e] >= (*capacity)[e] ) continue; 
   279 	Node v=g->head(e);            
   280 	    
   281 	if( lev > level[v] ) { //Push is allowed now
   282 	  
   283 	  if ( excess[v]<=0 && v!=t && v!=s ) {
   284 	    int lev_v=level[v];
   285 	    active[lev_v].push(v);
   286 	  }
   287 	  
   288 	  Num cap=(*capacity)[e];
   289 	  Num flo=(*flow)[e];
   290 	  Num remcap=cap-flo;
   291 	  
   292 	  if ( remcap >= exc ) { //A nonsaturating push.
   293 	    
   294 	    flow->set(e, flo+exc);
   295 	    excess.set(v, excess[v]+exc);
   296 	    exc=0;
   297 	    break; 
   298 	    
   299 	  } else { //A saturating push.
   300 	    flow->set(e, cap);
   301 	    excess.set(v, excess[v]+remcap);
   302 	    exc-=remcap;
   303 	  }
   304 	} else if ( newlevel > level[v] ) newlevel = level[v];
   305       } //for out edges wv 
   306       
   307       if ( exc > 0 ) {	
   308 	InEdgeIt e;
   309 	for(g->first(e,w); g->valid(e); g->next(e)) {
   310 	  
   311 	  if( (*flow)[e] <= 0 ) continue; 
   312 	  Node v=g->tail(e); 
   313 	  
   314 	  if( lev > level[v] ) { //Push is allowed now
   315 	    
   316 	    if ( excess[v]<=0 && v!=t && v!=s ) {
   317 	      int lev_v=level[v];
   318 	      active[lev_v].push(v);
   319 	    }
   320 	    
   321 	    Num flo=(*flow)[e];
   322 	    
   323 	    if ( flo >= exc ) { //A nonsaturating push.
   324 	      
   325 	      flow->set(e, flo-exc);
   326 	      excess.set(v, excess[v]+exc);
   327 	      exc=0;
   328 	      break; 
   329 	    } else {  //A saturating push.
   330 	      
   331 	      excess.set(v, excess[v]+flo);
   332 	      exc-=flo;
   333 	      flow->set(e,0);
   334 	    }  
   335 	  } else if ( newlevel > level[v] ) newlevel = level[v];
   336 	} //for in edges vw
   337 	
   338       } // if w still has excess after the out edge for cycle
   339       
   340       excess.set(w, exc);
   341       
   342       return newlevel;
   343     }
   344 
   345 
   346     void preflowPreproc ( flowEnum fe, VecStack& active, 
   347 			  VecNode& level_list, NNMap& left, NNMap& right ) {
   348 
   349 			    std::queue<Node> bfs_queue;
   350       
   351 			    switch ( fe ) {
   352 			    case ZERO_FLOW: 
   353 			      {
   354 				//Reverse_bfs from t, to find the starting level.
   355 				level.set(t,0);
   356 				bfs_queue.push(t);
   357 	
   358 				while (!bfs_queue.empty()) {
   359 	    
   360 				  Node v=bfs_queue.front();	
   361 				  bfs_queue.pop();
   362 				  int l=level[v]+1;
   363 	    
   364 				  InEdgeIt e;
   365 				  for(g->first(e,v); g->valid(e); g->next(e)) {
   366 				    Node w=g->tail(e);
   367 				    if ( level[w] == n && w != s ) {
   368 				      bfs_queue.push(w);
   369 				      Node first=level_list[l];
   370 				      if ( g->valid(first) ) left.set(first,w);
   371 				      right.set(w,first);
   372 				      level_list[l]=w;
   373 				      level.set(w, l);
   374 				    }
   375 				  }
   376 				}
   377 	  
   378 				//the starting flow
   379 				OutEdgeIt e;
   380 				for(g->first(e,s); g->valid(e); g->next(e)) 
   381 				  {
   382 				    Num c=(*capacity)[e];
   383 				    if ( c <= 0 ) continue;
   384 				    Node w=g->head(e);
   385 				    if ( level[w] < n ) {	  
   386 				      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   387 				      flow->set(e, c); 
   388 				      excess.set(w, excess[w]+c);
   389 				    }
   390 				  }
   391 				break;
   392 			      }
   393 	
   394 			    case GEN_FLOW:
   395 			    case PREFLOW:
   396 			      {
   397 				//Reverse_bfs from t in the residual graph, 
   398 				//to find the starting level.
   399 				level.set(t,0);
   400 				bfs_queue.push(t);
   401 	  
   402 				while (!bfs_queue.empty()) {
   403 	    
   404 				  Node v=bfs_queue.front();	
   405 				  bfs_queue.pop();
   406 				  int l=level[v]+1;
   407 	    
   408 				  InEdgeIt e;
   409 				  for(g->first(e,v); g->valid(e); g->next(e)) {
   410 				    if ( (*capacity)[e] <= (*flow)[e] ) continue;
   411 				    Node w=g->tail(e);
   412 				    if ( level[w] == n && w != s ) {
   413 				      bfs_queue.push(w);
   414 				      Node first=level_list[l];
   415 				      if ( g->valid(first) ) left.set(first,w);
   416 				      right.set(w,first);
   417 				      level_list[l]=w;
   418 				      level.set(w, l);
   419 				    }
   420 				  }
   421 	    
   422 				  OutEdgeIt f;
   423 				  for(g->first(f,v); g->valid(f); g->next(f)) {
   424 				    if ( 0 >= (*flow)[f] ) continue;
   425 				    Node w=g->head(f);
   426 				    if ( level[w] == n && w != s ) {
   427 				      bfs_queue.push(w);
   428 				      Node first=level_list[l];
   429 				      if ( g->valid(first) ) left.set(first,w);
   430 				      right.set(w,first);
   431 				      level_list[l]=w;
   432 				      level.set(w, l);
   433 				    }
   434 				  }
   435 				}
   436 	  
   437 	  
   438 				//the starting flow
   439 				OutEdgeIt e;
   440 				for(g->first(e,s); g->valid(e); g->next(e)) 
   441 				  {
   442 				    Num rem=(*capacity)[e]-(*flow)[e];
   443 				    if ( rem <= 0 ) continue;
   444 				    Node w=g->head(e);
   445 				    if ( level[w] < n ) {	  
   446 				      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   447 				      flow->set(e, (*capacity)[e]); 
   448 				      excess.set(w, excess[w]+rem);
   449 				    }
   450 				  }
   451 	  
   452 				InEdgeIt f;
   453 				for(g->first(f,s); g->valid(f); g->next(f)) 
   454 				  {
   455 				    if ( (*flow)[f] <= 0 ) continue;
   456 				    Node w=g->tail(f);
   457 				    if ( level[w] < n ) {	  
   458 				      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   459 				      excess.set(w, excess[w]+(*flow)[f]);
   460 				      flow->set(f, 0); 
   461 				    }
   462 				  }  
   463 				break;
   464 			      } //case PREFLOW
   465 			    }
   466 			  } //preflowPreproc
   467 
   468 
   469 
   470     void relabel(Node w, int newlevel, VecStack& active,  
   471 		 VecNode& level_list, NNMap& left, 
   472 		 NNMap& right, int& b, int& k, bool what_heur ) 
   473     {
   474 
   475       Num lev=level[w];	
   476       
   477       Node right_n=right[w];
   478       Node left_n=left[w];
   479       
   480       //unlacing starts
   481       if ( g->valid(right_n) ) {
   482 	if ( g->valid(left_n) ) {
   483 	  right.set(left_n, right_n);
   484 	  left.set(right_n, left_n);
   485 	} else {
   486 	  level_list[lev]=right_n;   
   487 	  left.set(right_n, INVALID);
   488 	} 
   489       } else {
   490 	if ( g->valid(left_n) ) {
   491 	  right.set(left_n, INVALID);
   492 	} else { 
   493 	  level_list[lev]=INVALID;   
   494 	} 
   495       } 
   496       //unlacing ends
   497 		
   498       if ( !g->valid(level_list[lev]) ) {
   499 	      
   500 	//gapping starts
   501 	for (int i=lev; i!=k ; ) {
   502 	  Node v=level_list[++i];
   503 	  while ( g->valid(v) ) {
   504 	    level.set(v,n);
   505 	    v=right[v];
   506 	  }
   507 	  level_list[i]=INVALID;
   508 	  if ( !what_heur ) {
   509 	    while ( !active[i].empty() ) {
   510 	      active[i].pop();    //FIXME: ezt szebben kene
   511 	    }
   512 	  }	     
   513 	}
   514 	
   515 	level.set(w,n);
   516 	b=lev-1;
   517 	k=b;
   518 	//gapping ends
   519 	
   520       } else {
   521 	
   522 	if ( newlevel == n ) level.set(w,n); 
   523 	else {
   524 	  level.set(w,++newlevel);
   525 	  active[newlevel].push(w);
   526 	  if ( what_heur ) b=newlevel;
   527 	  if ( k < newlevel ) ++k;      //now k=newlevel
   528 	  Node first=level_list[newlevel];
   529 	  if ( g->valid(first) ) left.set(first,w);
   530 	  right.set(w,first);
   531 	  left.set(w,INVALID);
   532 	  level_list[newlevel]=w;
   533 	}
   534       }
   535       
   536     } //relabel
   537 
   538 
   539     template<typename MapGraphWrapper> 
   540     class DistanceMap {
   541     protected:
   542       const MapGraphWrapper* g;
   543       typename MapGraphWrapper::template NodeMap<int> dist; 
   544     public:
   545       DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
   546       void set(const typename MapGraphWrapper::Node& n, int a) { 
   547 	dist.set(n, a); 
   548       }
   549       int operator[](const typename MapGraphWrapper::Node& n) 
   550       { return dist[n]; }
   551       //       int get(const typename MapGraphWrapper::Node& n) const { 
   552       // 	return dist[n]; }
   553       //       bool get(const typename MapGraphWrapper::Edge& e) const { 
   554       // 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
   555       bool operator[](const typename MapGraphWrapper::Edge& e) const { 
   556 	return (dist[g->tail(e)]<dist[g->head(e)]); 
   557       }
   558     };
   559     
   560   };
   561 
   562 
   563   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   564   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe ) 
   565   {
   566       
   567     int heur0=(int)(H0*n);  //time while running 'bound decrease' 
   568     int heur1=(int)(H1*n);  //time while running 'highest label'
   569     int heur=heur1;         //starting time interval (#of relabels)
   570     int numrelabel=0;
   571      
   572     bool what_heur=1;       
   573     //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   574 
   575     bool end=false;     
   576     //Needed for 'bound decrease', true means no active nodes are above bound b.
   577 
   578     int k=n-2;  //bound on the highest level under n containing a node
   579     int b=k;    //bound on the highest level under n of an active node
   580       
   581     VecStack active(n);
   582       
   583     NNMap left(*g, INVALID);
   584     NNMap right(*g, INVALID);
   585     VecNode level_list(n,INVALID);
   586     //List of the nodes in level i<n, set to n.
   587 
   588     NodeIt v;
   589     for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   590     //setting each node to level n
   591       
   592     switch ( fe ) {
   593     case PREFLOW:
   594       {
   595 	//counting the excess
   596 	NodeIt v;
   597 	for(g->first(v); g->valid(v); g->next(v)) {
   598 	  Num exc=0;
   599 	  
   600 	  InEdgeIt e;
   601 	  for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
   602 	  OutEdgeIt f;
   603 	  for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
   604 	    
   605 	  excess.set(v,exc);	  
   606 	    
   607 	  //putting the active nodes into the stack
   608 	  int lev=level[v];
   609 	  if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
   610 	}
   611 	break;
   612       }
   613     case GEN_FLOW:
   614       {
   615 	//Counting the excess of t
   616 	Num exc=0;
   617 	  
   618 	InEdgeIt e;
   619 	for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   620 	OutEdgeIt f;
   621 	for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
   622 	  
   623 	excess.set(t,exc);	
   624 	  
   625 	break;
   626       }
   627     default:
   628       break;
   629     }
   630       
   631     preflowPreproc( fe, active, level_list, left, right );
   632     //End of preprocessing 
   633       
   634       
   635     //Push/relabel on the highest level active nodes.
   636     while ( true ) {
   637       if ( b == 0 ) {
   638 	if ( !what_heur && !end && k > 0 ) {
   639 	  b=k;
   640 	  end=true;
   641 	} else break;
   642       }
   643 	
   644       if ( active[b].empty() ) --b; 
   645       else {
   646 	end=false;  
   647 	Node w=active[b].top();
   648 	active[b].pop();
   649 	int newlevel=push(w,active);
   650 	if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 
   651 				     left, right, b, k, what_heur);
   652 	  
   653 	++numrelabel; 
   654 	if ( numrelabel >= heur ) {
   655 	  numrelabel=0;
   656 	  if ( what_heur ) {
   657 	    what_heur=0;
   658 	    heur=heur0;
   659 	    end=false;
   660 	  } else {
   661 	    what_heur=1;
   662 	    heur=heur1;
   663 	    b=k; 
   664 	  }
   665 	}
   666       } 
   667     } 
   668   }
   669 
   670 
   671 
   672   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   673   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1() 
   674   {
   675       
   676     int k=n-2;  //bound on the highest level under n containing a node
   677     int b=k;    //bound on the highest level under n of an active node
   678       
   679     VecStack active(n);
   680     level.set(s,0);
   681     std::queue<Node> bfs_queue;
   682     bfs_queue.push(s);
   683 	    
   684     while (!bfs_queue.empty()) {
   685 	
   686       Node v=bfs_queue.front();	
   687       bfs_queue.pop();
   688       int l=level[v]+1;
   689 	      
   690       InEdgeIt e;
   691       for(g->first(e,v); g->valid(e); g->next(e)) {
   692 	if ( (*capacity)[e] <= (*flow)[e] ) continue;
   693 	Node u=g->tail(e);
   694 	if ( level[u] >= n ) { 
   695 	  bfs_queue.push(u);
   696 	  level.set(u, l);
   697 	  if ( excess[u] > 0 ) active[l].push(u);
   698 	}
   699       }
   700 	
   701       OutEdgeIt f;
   702       for(g->first(f,v); g->valid(f); g->next(f)) {
   703 	if ( 0 >= (*flow)[f] ) continue;
   704 	Node u=g->head(f);
   705 	if ( level[u] >= n ) { 
   706 	  bfs_queue.push(u);
   707 	  level.set(u, l);
   708 	  if ( excess[u] > 0 ) active[l].push(u);
   709 	}
   710       }
   711     }
   712     b=n-2;
   713 
   714     while ( true ) {
   715 	
   716       if ( b == 0 ) break;
   717 
   718       if ( active[b].empty() ) --b; 
   719       else {
   720 	Node w=active[b].top();
   721 	active[b].pop();
   722 	int newlevel=push(w,active);	  
   723 
   724 	//relabel
   725 	if ( excess[w] > 0 ) {
   726 	  level.set(w,++newlevel);
   727 	  active[newlevel].push(w);
   728 	  b=newlevel;
   729 	}
   730       }  // if stack[b] is nonempty
   731     } // while(true)
   732   }
   733 
   734 
   735 
   736   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   737   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() 
   738   {
   739     ResGW res_graph(*g, *capacity, *flow);
   740     bool _augment=false;
   741       
   742     //ReachedMap level(res_graph);
   743     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   744     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   745     bfs.pushAndSetReached(s);
   746 	
   747     typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); 
   748     pred.set(s, INVALID);
   749       
   750     typename ResGW::template NodeMap<Num> free(res_graph);
   751 	
   752     //searching for augmenting path
   753     while ( !bfs.finished() ) { 
   754       ResGWOutEdgeIt e=bfs;
   755       if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   756 	Node v=res_graph.tail(e);
   757 	Node w=res_graph.head(e);
   758 	pred.set(w, e);
   759 	if (res_graph.valid(pred[v])) {
   760 	  free.set(w, std::min(free[v], res_graph.resCap(e)));
   761 	} else {
   762 	  free.set(w, res_graph.resCap(e)); 
   763 	}
   764 	if (res_graph.head(e)==t) { _augment=true; break; }
   765       }
   766 	
   767       ++bfs;
   768     } //end of searching augmenting path
   769 
   770     if (_augment) {
   771       Node n=t;
   772       Num augment_value=free[t];
   773       while (res_graph.valid(pred[n])) { 
   774 	ResGWEdge e=pred[n];
   775 	res_graph.augment(e, augment_value); 
   776 	n=res_graph.tail(e);
   777       }
   778     }
   779 
   780     return _augment;
   781   }
   782 
   783 
   784 
   785 
   786 
   787 
   788 
   789 
   790 
   791   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   792   template<typename MutableGraph> 
   793   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() 
   794   {      
   795     typedef MutableGraph MG;
   796     bool _augment=false;
   797 
   798     ResGW res_graph(*g, *capacity, *flow);
   799 
   800     //bfs for distances on the residual graph
   801     //ReachedMap level(res_graph);
   802     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   803     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   804     bfs.pushAndSetReached(s);
   805     typename ResGW::template NodeMap<int> 
   806       dist(res_graph); //filled up with 0's
   807 
   808     //F will contain the physical copy of the residual graph
   809     //with the set of edges which are on shortest paths
   810     MG F;
   811     typename ResGW::template NodeMap<typename MG::Node> 
   812       res_graph_to_F(res_graph);
   813     {
   814       typename ResGW::NodeIt n;
   815       for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
   816 	res_graph_to_F.set(n, F.addNode());
   817       }
   818     }
   819 
   820     typename MG::Node sF=res_graph_to_F[s];
   821     typename MG::Node tF=res_graph_to_F[t];
   822     typename MG::template EdgeMap<ResGWEdge> original_edge(F);
   823     typename MG::template EdgeMap<Num> residual_capacity(F);
   824 
   825     while ( !bfs.finished() ) { 
   826       ResGWOutEdgeIt e=bfs;
   827       if (res_graph.valid(e)) {
   828 	if (bfs.isBNodeNewlyReached()) {
   829 	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   830 	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   831 	  original_edge.update();
   832 	  original_edge.set(f, e);
   833 	  residual_capacity.update();
   834 	  residual_capacity.set(f, res_graph.resCap(e));
   835 	} else {
   836 	  if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
   837 	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   838 	    original_edge.update();
   839 	    original_edge.set(f, e);
   840 	    residual_capacity.update();
   841 	    residual_capacity.set(f, res_graph.resCap(e));
   842 	  }
   843 	}
   844       }
   845       ++bfs;
   846     } //computing distances from s in the residual graph
   847 
   848     bool __augment=true;
   849 
   850     while (__augment) {
   851       __augment=false;
   852       //computing blocking flow with dfs
   853       DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
   854       typename MG::template NodeMap<typename MG::Edge> pred(F);
   855       pred.set(sF, INVALID);
   856       //invalid iterators for sources
   857 
   858       typename MG::template NodeMap<Num> free(F);
   859 
   860       dfs.pushAndSetReached(sF);      
   861       while (!dfs.finished()) {
   862 	++dfs;
   863 	if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
   864 	  if (dfs.isBNodeNewlyReached()) {
   865 	    typename MG::Node v=F.aNode(dfs);
   866 	    typename MG::Node w=F.bNode(dfs);
   867 	    pred.set(w, dfs);
   868 	    if (F.valid(pred[v])) {
   869 	      free.set(w, std::min(free[v], residual_capacity[dfs]));
   870 	    } else {
   871 	      free.set(w, residual_capacity[dfs]); 
   872 	    }
   873 	    if (w==tF) { 
   874 	      __augment=true; 
   875 	      _augment=true;
   876 	      break; 
   877 	    }
   878 	      
   879 	  } else {
   880 	    F.erase(/*typename MG::OutEdgeIt*/(dfs));
   881 	  }
   882 	} 
   883       }
   884 
   885       if (__augment) {
   886 	typename MG::Node n=tF;
   887 	Num augment_value=free[tF];
   888 	while (F.valid(pred[n])) { 
   889 	  typename MG::Edge e=pred[n];
   890 	  res_graph.augment(original_edge[e], augment_value); 
   891 	  n=F.tail(e);
   892 	  if (residual_capacity[e]==augment_value) 
   893 	    F.erase(e); 
   894 	  else 
   895 	    residual_capacity.set(e, residual_capacity[e]-augment_value);
   896 	}
   897       }
   898 	
   899     }
   900             
   901     return _augment;
   902   }
   903 
   904 
   905 
   906 
   907 
   908 
   909   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   910   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() 
   911   {
   912     bool _augment=false;
   913 
   914     ResGW res_graph(*g, *capacity, *flow);
   915       
   916     //ReachedMap level(res_graph);
   917     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   918     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   919 
   920     bfs.pushAndSetReached(s);
   921     DistanceMap<ResGW> dist(res_graph);
   922     while ( !bfs.finished() ) { 
   923       ResGWOutEdgeIt e=bfs;
   924       if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   925 	dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   926       }
   927       ++bfs;
   928     } //computing distances from s in the residual graph
   929 
   930       //Subgraph containing the edges on some shortest paths
   931     ConstMap<typename ResGW::Node, bool> true_map(true);
   932     typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, 
   933       DistanceMap<ResGW> > FilterResGW;
   934     FilterResGW filter_res_graph(res_graph, true_map, dist);
   935 
   936     //Subgraph, which is able to delete edges which are already 
   937     //met by the dfs
   938     typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> 
   939       first_out_edges(filter_res_graph);
   940     typename FilterResGW::NodeIt v;
   941     for(filter_res_graph.first(v); filter_res_graph.valid(v); 
   942 	filter_res_graph.next(v)) 
   943       {
   944  	typename FilterResGW::OutEdgeIt e;
   945  	filter_res_graph.first(e, v);
   946  	first_out_edges.set(v, e);
   947       }
   948     typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
   949       template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
   950     ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
   951 
   952     bool __augment=true;
   953 
   954     while (__augment) {
   955 
   956       __augment=false;
   957       //computing blocking flow with dfs
   958       DfsIterator< ErasingResGW, 
   959 	typename ErasingResGW::template NodeMap<bool> > 
   960 	dfs(erasing_res_graph);
   961       typename ErasingResGW::
   962 	template NodeMap<typename ErasingResGW::OutEdgeIt> 
   963 	pred(erasing_res_graph); 
   964       pred.set(s, INVALID);
   965       //invalid iterators for sources
   966 
   967       typename ErasingResGW::template NodeMap<Num> 
   968 	free1(erasing_res_graph);
   969 
   970       dfs.pushAndSetReached(
   971 			    typename ErasingResGW::Node(
   972 							typename FilterResGW::Node(
   973 										   typename ResGW::Node(s)
   974 										   )
   975 							)
   976 			    );
   977       while (!dfs.finished()) {
   978 	++dfs;
   979 	if (erasing_res_graph.valid(
   980 				    typename ErasingResGW::OutEdgeIt(dfs))) 
   981  	  { 
   982   	    if (dfs.isBNodeNewlyReached()) {
   983 	  
   984  	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
   985  	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
   986 
   987  	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
   988  	      if (erasing_res_graph.valid(pred[v])) {
   989  		free1.set(w, std::min(free1[v], res_graph.resCap(
   990 								 typename ErasingResGW::OutEdgeIt(dfs))));
   991  	      } else {
   992  		free1.set(w, res_graph.resCap(
   993 					      typename ErasingResGW::OutEdgeIt(dfs))); 
   994  	      }
   995 	      
   996  	      if (w==t) { 
   997  		__augment=true; 
   998  		_augment=true;
   999  		break; 
  1000  	      }
  1001  	    } else {
  1002  	      erasing_res_graph.erase(dfs);
  1003 	    }
  1004 	  }
  1005       }	
  1006 
  1007       if (__augment) {
  1008 	typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
  1009 	// 	  typename ResGW::NodeMap<Num> a(res_graph);
  1010 	// 	  typename ResGW::Node b;
  1011 	// 	  Num j=a[b];
  1012 	// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
  1013 	// 	  typename FilterResGW::Node b1;
  1014 	// 	  Num j1=a1[b1];
  1015 	// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
  1016 	// 	  typename ErasingResGW::Node b2;
  1017 	// 	  Num j2=a2[b2];
  1018 	Num augment_value=free1[n];
  1019 	while (erasing_res_graph.valid(pred[n])) { 
  1020 	  typename ErasingResGW::OutEdgeIt e=pred[n];
  1021 	  res_graph.augment(e, augment_value);
  1022 	  n=erasing_res_graph.tail(e);
  1023 	  if (res_graph.resCap(e)==0)
  1024 	    erasing_res_graph.erase(e);
  1025 	}
  1026       }
  1027       
  1028     } //while (__augment) 
  1029             
  1030     return _augment;
  1031   }
  1032 
  1033 
  1034 
  1035 
  1036 } //namespace hugo
  1037 
  1038 #endif //HUGO_MAX_FLOW_H
  1039 
  1040 
  1041 
  1042