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