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