(none)
authorjacint
Wed, 12 May 2004 10:51:53 +0000
changeset 6206e917be931af
parent 619 e09818232531
child 621 2db02d4a9e6e
(none)
src/work/jacint/max_save.h
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/work/jacint/max_save.h	Wed May 12 10:51:53 2004 +0000
     1.3 @@ -0,0 +1,1136 @@
     1.4 +// -*- C++ -*-
     1.5 +#ifndef HUGO_MAX_FLOW_H
     1.6 +#define HUGO_MAX_FLOW_H
     1.7 +
     1.8 +///\ingroup galgs
     1.9 +///\file
    1.10 +///\brief Maximum flow algorithm.
    1.11 +
    1.12 +#define H0 20
    1.13 +#define H1 1
    1.14 +
    1.15 +#include <vector>
    1.16 +#include <queue>
    1.17 +#include <stack>
    1.18 +
    1.19 +#include <graph_wrapper.h>
    1.20 +#include <bfs_iterator.h>
    1.21 +#include <invalid.h>
    1.22 +#include <maps.h>
    1.23 +#include <for_each_macros.h>
    1.24 +
    1.25 +/// \file
    1.26 +/// \brief Dimacs file format reader.
    1.27 +
    1.28 +namespace hugo {
    1.29 +
    1.30 +  /// \addtogroup galgs
    1.31 +  /// @{
    1.32 +
    1.33 +  ///Maximum flow algorithms class.
    1.34 +
    1.35 +  ///This class provides various algorithms for finding a flow of
    1.36 +  ///maximum value in a directed graph. The \e source node, the \e
    1.37 +  ///target node, the \e capacity of the edges and the \e starting \e
    1.38 +  ///flow value of the edges can be passed to the algorithm by the
    1.39 +  ///constructor. It is possible to change these quantities using the
    1.40 +  ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
    1.41 +  ///\ref resetFlow. Before any subsequent runs of any algorithm of
    1.42 +  ///the class \ref resetFlow should be called, otherwise it will
    1.43 +  ///start from a maximum flow.
    1.44 +
    1.45 +  ///After running an algorithm of the class, the maximum value of a
    1.46 +  ///value can be obtained by calling \ref flowValue(). The minimum
    1.47 +  ///value cut can be written into a \c node map of \c bools by
    1.48 +  ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
    1.49 +  ///the inclusionwise minimum and maximum of the minimum value
    1.50 +  ///cuts, resp.)
    1.51 +
    1.52 +  ///\param Graph The undirected graph type the algorithm runs on.
    1.53 +  ///\param Num The number type of the capacities and the flow values.
    1.54 +  ///\param The type of the capacity map.
    1.55 +  ///\param The type of the flow map.
    1.56 +
    1.57 +  ///\author Marton Makai, Jacint Szabo
    1.58 +  template <typename Graph, typename Num, 
    1.59 +	    typename CapMap=typename Graph::template EdgeMap<Num>, 
    1.60 +            typename FlowMap=typename Graph::template EdgeMap<Num> >
    1.61 +  class MaxFlow {
    1.62 +    
    1.63 +    typedef typename Graph::Node Node;
    1.64 +    typedef typename Graph::NodeIt NodeIt;
    1.65 +    typedef typename Graph::OutEdgeIt OutEdgeIt;
    1.66 +    typedef typename Graph::InEdgeIt InEdgeIt;
    1.67 +
    1.68 +    typedef typename std::vector<std::stack<Node> > VecStack;
    1.69 +    typedef typename Graph::template NodeMap<Node> NNMap;
    1.70 +    typedef typename std::vector<Node> VecNode;
    1.71 +    
    1.72 +    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    1.73 +    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    1.74 +    typedef typename ResGW::Edge ResGWEdge;
    1.75 +    //typedef typename ResGW::template NodeMap<bool> ReachedMap;  //fixme
    1.76 +    typedef typename Graph::template NodeMap<int> ReachedMap;
    1.77 +    
    1.78 +    const Graph* g;
    1.79 +    Node s;
    1.80 +    Node t;
    1.81 +    const CapMap* capacity;  
    1.82 +    FlowMap* flow;
    1.83 +    int n;          //the number of nodes of G
    1.84 +
    1.85 +    //level works as a bool map in augmenting path algorithms and is
    1.86 +    //used by bfs for storing reached information.  In preflow, it
    1.87 +    //shows the levels of nodes. 
    1.88 +    ReachedMap level;
    1.89 +    
    1.90 +    //excess is needed only in preflow
    1.91 +    typename Graph::template NodeMap<Num> excess; 
    1.92 +
    1.93 +
    1.94 +    //fixme
    1.95 +    //   protected:
    1.96 +    //     MaxFlow() { }
    1.97 +    //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, 
    1.98 +    // 	     FlowMap& _flow) 
    1.99 +    //       {
   1.100 +    // 	g=&_G; 
   1.101 +    // 	s=_s; 
   1.102 +    // 	t=_t; 
   1.103 +    // 	capacity=&_capacity;
   1.104 +    // 	flow=&_flow;
   1.105 +    // 	n=_G.nodeNum;
   1.106 +    // 	level.set (_G); //kellene vmi ilyesmi fv 
   1.107 +    // 	excess(_G,0); //itt is
   1.108 +    //       }
   1.109 +
   1.110 +  public:
   1.111 + 
   1.112 +    ///Indicates the property of the starting flow. 
   1.113 +
   1.114 +    ///Indicates the property of the starting flow. The meanings: 
   1.115 +    ///- \c ZERO_FLOW: constant zero flow
   1.116 +    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
   1.117 +    ///the sum of the out-flows in every node except the source and
   1.118 +    ///the target.
   1.119 +    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
   1.120 +    ///least the sum of the out-flows in every node except the source.
   1.121 +    enum flowEnum{
   1.122 +      ZERO_FLOW=0,
   1.123 +      GEN_FLOW=1,
   1.124 +      PRE_FLOW=2
   1.125 +    };
   1.126 +
   1.127 +    MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, 
   1.128 +	    FlowMap& _flow) :
   1.129 +      g(&_G), s(_s), t(_t), capacity(&_capacity), 
   1.130 +      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
   1.131 +
   1.132 +    ///Runs a maximum flow algorithm.
   1.133 +
   1.134 +    ///Runs a preflow algorithm, which is the fastest maximum flow
   1.135 +    ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
   1.136 +    ///\pre The starting flow must be a 
   1.137 +    /// - constant zero flow if \c fe is \c ZERO_FLOW,
   1.138 +    /// - an arbitary flow if \c fe is \c GEN_FLOW, 
   1.139 +    /// - an arbitary preflow if \c fe is \c PRE_FLOW.
   1.140 +    void run( flowEnum fe=ZERO_FLOW ) {
   1.141 +      preflow(fe);
   1.142 +    }
   1.143 +    
   1.144 +    ///Runs a preflow algorithm.
   1.145 +
   1.146 +    ///Runs a preflow algorithm. The preflow algorithms provide the
   1.147 +    ///fastest way to compute a maximum flow in a directed graph.
   1.148 +    ///\pre The starting flow must be a 
   1.149 +    /// - constant zero flow if \c fe is \c ZERO_FLOW,
   1.150 +    /// - an arbitary flow if \c fe is \c GEN_FLOW, 
   1.151 +    /// - an arbitary preflow if \c fe is \c PRE_FLOW.
   1.152 +    void preflow(flowEnum fe) {
   1.153 +      preflowPhase1(fe);
   1.154 +      preflowPhase2();
   1.155 +    }
   1.156 +    // Heuristics: 
   1.157 +    //   2 phase
   1.158 +    //   gap
   1.159 +    //   list 'level_list' on the nodes on level i implemented by hand
   1.160 +    //   stack 'active' on the active nodes on level i
   1.161 +    //   runs heuristic 'highest label' for H1*n relabels
   1.162 +    //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
   1.163 +    //   Parameters H0 and H1 are initialized to 20 and 1.
   1.164 +
   1.165 +    ///Runs the first phase of the preflow algorithm.
   1.166 +    
   1.167 +    ///The preflow algorithm consists of two phases, this method runs the
   1.168 +    ///first phase. After the first phase the maximum flow value and a
   1.169 +    ///minimum value cut can already be computed, though a maximum flow
   1.170 +    ///is net yet obtained. So after calling this method \ref flowValue
   1.171 +    ///and \ref actMinCut gives proper results. 
   1.172 +    ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
   1.173 +    ///give minimum value cuts unless calling \ref preflowPhase2.
   1.174 +    ///\pre The starting flow must be a 
   1.175 +    /// - constant zero flow if \c fe is \c ZERO_FLOW,
   1.176 +    /// - an arbitary flow if \c fe is \c GEN_FLOW, 
   1.177 +    /// - an arbitary preflow if \c fe is \c PRE_FLOW.
   1.178 +    void preflowPhase1( flowEnum fe );
   1.179 +
   1.180 +    ///Runs the second phase of the preflow algorithm.
   1.181 +    
   1.182 +    ///The preflow algorithm consists of two phases, this method runs
   1.183 +    ///the second phase. After calling \ref preflowPhase1 and then
   1.184 +    ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
   1.185 +    ///\ref minMinCut and \ref maxMinCut give proper results.
   1.186 +    ///\pre \ref preflowPhase1 must be called before.
   1.187 +    void preflowPhase2();
   1.188 +
   1.189 +    /// Starting from a flow, this method searches for an augmenting path 
   1.190 +    /// according to the Edmonds-Karp algorithm 
   1.191 +    /// and augments the flow on if any. 
   1.192 +    /// The return value shows if the augmentation was successful.
   1.193 +    bool augmentOnShortestPath();
   1.194 +
   1.195 +    /// Starting from a flow, this method searches for an augmenting blockin 
   1.196 +    /// flow according to Dinits' algorithm and augments the flow on if any. 
   1.197 +    /// The blocking flow is computed in a physically constructed 
   1.198 +    /// residual graph of type \c Mutablegraph.
   1.199 +    /// The return value show sif the augmentation was succesful.
   1.200 +    template<typename MutableGraph> bool augmentOnBlockingFlow();
   1.201 +
   1.202 +    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the 
   1.203 +    /// residual graph is not constructed physically.
   1.204 +    /// The return value shows if the augmentation was succesful.
   1.205 +    bool augmentOnBlockingFlow2();
   1.206 +
   1.207 +    /// Returns the actual flow value.
   1.208 +    /// More precisely, it returns the negative excess of s, thus 
   1.209 +    /// this works also for preflows.
   1.210 +    ///Can be called already after \ref preflowPhase1.
   1.211 +
   1.212 +    Num flowValue() { 
   1.213 +      Num a=0;
   1.214 +      FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
   1.215 +      FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
   1.216 +      return a;
   1.217 +      //marci figyu: excess[t] epp ezt adja preflow 0. fazisa utan
   1.218 +    }
   1.219 +
   1.220 +    ///Returns a minimum value cut after calling \ref preflowPhase1.
   1.221 +
   1.222 +    ///After the first phase of the preflow algorithm the maximum flow
   1.223 +    ///value and a minimum value cut can already be computed. This
   1.224 +    ///method can be called after running \ref preflowPhase1 for
   1.225 +    ///obtaining a minimum value cut.
   1.226 +    ///\warning: Gives proper result only right after calling \ref
   1.227 +    ///preflowPhase1.
   1.228 +    ///\todo We have to make some status variable which shows the actual state 
   1.229 +    /// of the class. This enables us to determine which methods are valid 
   1.230 +    /// for MinCut computation
   1.231 +    template<typename _CutMap>
   1.232 +    void actMinCut(_CutMap& M) {
   1.233 +      NodeIt v;
   1.234 +      for(g->first(v); g->valid(v); g->next(v)) {
   1.235 +	if ( level[v] < n ) {
   1.236 +	  M.set(v,false);
   1.237 +	} else {
   1.238 +	  M.set(v,true);
   1.239 +	}
   1.240 +      }
   1.241 +    }
   1.242 +    
   1.243 +    ///Returns the inclusionwise minimum of the minimum value cuts.
   1.244 +
   1.245 +    ///Sets \c M to the characteristic vector of the minimum value cut
   1.246 +    ///which is inclusionwise minimum. It is computed by processing
   1.247 +    ///a bfs from the source node \c s in the residual graph.
   1.248 +    ///\pre M should be a node map of bools initialized to false.
   1.249 +    ///\pre \c flow must be a maximum flow.
   1.250 +    template<typename _CutMap>
   1.251 +    void minMinCut(_CutMap& M) {
   1.252 +    
   1.253 +      std::queue<Node> queue;
   1.254 +      
   1.255 +      M.set(s,true);      
   1.256 +      queue.push(s);
   1.257 +
   1.258 +      while (!queue.empty()) {
   1.259 +        Node w=queue.front();
   1.260 +	queue.pop();
   1.261 +
   1.262 +	OutEdgeIt e;
   1.263 +	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   1.264 +	  Node v=g->head(e);
   1.265 +	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   1.266 +	    queue.push(v);
   1.267 +	    M.set(v, true);
   1.268 +	  }
   1.269 +	} 
   1.270 +
   1.271 +	InEdgeIt f;
   1.272 +	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   1.273 +	  Node v=g->tail(f);
   1.274 +	  if (!M[v] && (*flow)[f] > 0 ) {
   1.275 +	    queue.push(v);
   1.276 +	    M.set(v, true);
   1.277 +	  }
   1.278 +	} 
   1.279 +      }
   1.280 +    }
   1.281 +
   1.282 +
   1.283 +    ///Returns the inclusionwise maximum of the minimum value cuts.
   1.284 +
   1.285 +    ///Sets \c M to the characteristic vector of the minimum value cut
   1.286 +    ///which is inclusionwise maximum. It is computed by processing a
   1.287 +    ///backward bfs from the target node \c t in the residual graph.
   1.288 +    ///\pre M should be a node map of bools initialized to false.
   1.289 +    ///\pre \c flow must be a maximum flow.
   1.290 +    template<typename _CutMap>
   1.291 +    void maxMinCut(_CutMap& M) {
   1.292 +
   1.293 +      NodeIt v;
   1.294 +      for(g->first(v) ; g->valid(v); g->next(v)) {
   1.295 +	M.set(v, true);
   1.296 +      }
   1.297 +
   1.298 +      std::queue<Node> queue;
   1.299 +      
   1.300 +      M.set(t,false);        
   1.301 +      queue.push(t);
   1.302 +
   1.303 +      while (!queue.empty()) {
   1.304 +        Node w=queue.front();
   1.305 +	queue.pop();
   1.306 +
   1.307 +
   1.308 +	InEdgeIt e;
   1.309 +	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   1.310 +	  Node v=g->tail(e);
   1.311 +	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   1.312 +	    queue.push(v);
   1.313 +	    M.set(v, false);
   1.314 +	  }
   1.315 +	}
   1.316 +	
   1.317 +	OutEdgeIt f;
   1.318 +	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   1.319 +	  Node v=g->head(f);
   1.320 +	  if (M[v] && (*flow)[f] > 0 ) {
   1.321 +	    queue.push(v);
   1.322 +	    M.set(v, false);
   1.323 +	  }
   1.324 +	}
   1.325 +      }
   1.326 +    }
   1.327 +
   1.328 +
   1.329 +    ///Returns a minimum value cut.
   1.330 +
   1.331 +    ///Sets \c M to the characteristic vector of a minimum value cut.
   1.332 +    ///\pre M should be a node map of bools initialized to false.
   1.333 +    ///\pre \c flow must be a maximum flow.
   1.334 +    template<typename CutMap>
   1.335 +    void minCut(CutMap& M) { minMinCut(M); }
   1.336 +
   1.337 +    ///Resets the source node to \c _s.
   1.338 +
   1.339 +    ///Resets the source node to \c _s.
   1.340 +    ///
   1.341 +    void resetSource(Node _s) { s=_s; }
   1.342 +
   1.343 +
   1.344 +    ///Resets the target node to \c _t.
   1.345 +
   1.346 +    ///Resets the target node to \c _t.
   1.347 +    ///
   1.348 +    void resetTarget(Node _t) { t=_t; }
   1.349 +   
   1.350 +    /// Resets the edge map of the capacities to _cap.
   1.351 +
   1.352 +    /// Resets the edge map of the capacities to _cap.
   1.353 +    ///
   1.354 +    void resetCap(const CapMap& _cap) { capacity=&_cap; }
   1.355 +    
   1.356 +    /// Resets the edge map of the flows to _flow.
   1.357 +
   1.358 +    /// Resets the edge map of the flows to _flow.
   1.359 +    ///
   1.360 +    void resetFlow(FlowMap& _flow) { flow=&_flow; }
   1.361 +
   1.362 +
   1.363 +  private:
   1.364 +
   1.365 +    int push(Node w, VecStack& active) {
   1.366 +      
   1.367 +      int lev=level[w];
   1.368 +      Num exc=excess[w];
   1.369 +      int newlevel=n;       //bound on the next level of w
   1.370 +	  
   1.371 +      OutEdgeIt e;
   1.372 +      for(g->first(e,w); g->valid(e); g->next(e)) {
   1.373 +	    
   1.374 +	if ( (*flow)[e] >= (*capacity)[e] ) continue; 
   1.375 +	Node v=g->head(e);            
   1.376 +	    
   1.377 +	if( lev > level[v] ) { //Push is allowed now
   1.378 +	  
   1.379 +	  if ( excess[v]<=0 && v!=t && v!=s ) {
   1.380 +	    int lev_v=level[v];
   1.381 +	    active[lev_v].push(v);
   1.382 +	  }
   1.383 +	  
   1.384 +	  Num cap=(*capacity)[e];
   1.385 +	  Num flo=(*flow)[e];
   1.386 +	  Num remcap=cap-flo;
   1.387 +	  
   1.388 +	  if ( remcap >= exc ) { //A nonsaturating push.
   1.389 +	    
   1.390 +	    flow->set(e, flo+exc);
   1.391 +	    excess.set(v, excess[v]+exc);
   1.392 +	    exc=0;
   1.393 +	    break; 
   1.394 +	    
   1.395 +	  } else { //A saturating push.
   1.396 +	    flow->set(e, cap);
   1.397 +	    excess.set(v, excess[v]+remcap);
   1.398 +	    exc-=remcap;
   1.399 +	  }
   1.400 +	} else if ( newlevel > level[v] ) newlevel = level[v];
   1.401 +      } //for out edges wv 
   1.402 +      
   1.403 +      if ( exc > 0 ) {	
   1.404 +	InEdgeIt e;
   1.405 +	for(g->first(e,w); g->valid(e); g->next(e)) {
   1.406 +	  
   1.407 +	  if( (*flow)[e] <= 0 ) continue; 
   1.408 +	  Node v=g->tail(e); 
   1.409 +	  
   1.410 +	  if( lev > level[v] ) { //Push is allowed now
   1.411 +	    
   1.412 +	    if ( excess[v]<=0 && v!=t && v!=s ) {
   1.413 +	      int lev_v=level[v];
   1.414 +	      active[lev_v].push(v);
   1.415 +	    }
   1.416 +	    
   1.417 +	    Num flo=(*flow)[e];
   1.418 +	    
   1.419 +	    if ( flo >= exc ) { //A nonsaturating push.
   1.420 +	      
   1.421 +	      flow->set(e, flo-exc);
   1.422 +	      excess.set(v, excess[v]+exc);
   1.423 +	      exc=0;
   1.424 +	      break; 
   1.425 +	    } else {  //A saturating push.
   1.426 +	      
   1.427 +	      excess.set(v, excess[v]+flo);
   1.428 +	      exc-=flo;
   1.429 +	      flow->set(e,0);
   1.430 +	    }  
   1.431 +	  } else if ( newlevel > level[v] ) newlevel = level[v];
   1.432 +	} //for in edges vw
   1.433 +	
   1.434 +      } // if w still has excess after the out edge for cycle
   1.435 +      
   1.436 +      excess.set(w, exc);
   1.437 +      
   1.438 +      return newlevel;
   1.439 +    }
   1.440 +
   1.441 +
   1.442 +    void preflowPreproc ( flowEnum fe, VecStack& active, 
   1.443 +			  VecNode& level_list, NNMap& left, NNMap& right ) {
   1.444 +
   1.445 +			    std::queue<Node> bfs_queue;
   1.446 +      
   1.447 +			    switch ( fe ) {
   1.448 +			    case ZERO_FLOW: 
   1.449 +			      {
   1.450 +				//Reverse_bfs from t, to find the starting level.
   1.451 +				level.set(t,0);
   1.452 +				bfs_queue.push(t);
   1.453 +	
   1.454 +				while (!bfs_queue.empty()) {
   1.455 +	    
   1.456 +				  Node v=bfs_queue.front();	
   1.457 +				  bfs_queue.pop();
   1.458 +				  int l=level[v]+1;
   1.459 +	    
   1.460 +				  InEdgeIt e;
   1.461 +				  for(g->first(e,v); g->valid(e); g->next(e)) {
   1.462 +				    Node w=g->tail(e);
   1.463 +				    if ( level[w] == n && w != s ) {
   1.464 +				      bfs_queue.push(w);
   1.465 +				      Node first=level_list[l];
   1.466 +				      if ( g->valid(first) ) left.set(first,w);
   1.467 +				      right.set(w,first);
   1.468 +				      level_list[l]=w;
   1.469 +				      level.set(w, l);
   1.470 +				    }
   1.471 +				  }
   1.472 +				}
   1.473 +	  
   1.474 +				//the starting flow
   1.475 +				OutEdgeIt e;
   1.476 +				for(g->first(e,s); g->valid(e); g->next(e)) 
   1.477 +				  {
   1.478 +				    Num c=(*capacity)[e];
   1.479 +				    if ( c <= 0 ) continue;
   1.480 +				    Node w=g->head(e);
   1.481 +				    if ( level[w] < n ) {	  
   1.482 +				      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   1.483 +				      flow->set(e, c); 
   1.484 +				      excess.set(w, excess[w]+c);
   1.485 +				    }
   1.486 +				  }
   1.487 +				break;
   1.488 +			      }
   1.489 +	
   1.490 +			    case GEN_FLOW:
   1.491 +			    case PRE_FLOW:
   1.492 +			      {
   1.493 +				//Reverse_bfs from t in the residual graph, 
   1.494 +				//to find the starting level.
   1.495 +				level.set(t,0);
   1.496 +				bfs_queue.push(t);
   1.497 +	  
   1.498 +				while (!bfs_queue.empty()) {
   1.499 +	    
   1.500 +				  Node v=bfs_queue.front();	
   1.501 +				  bfs_queue.pop();
   1.502 +				  int l=level[v]+1;
   1.503 +	    
   1.504 +				  InEdgeIt e;
   1.505 +				  for(g->first(e,v); g->valid(e); g->next(e)) {
   1.506 +				    if ( (*capacity)[e] <= (*flow)[e] ) continue;
   1.507 +				    Node w=g->tail(e);
   1.508 +				    if ( level[w] == n && w != s ) {
   1.509 +				      bfs_queue.push(w);
   1.510 +				      Node first=level_list[l];
   1.511 +				      if ( g->valid(first) ) left.set(first,w);
   1.512 +				      right.set(w,first);
   1.513 +				      level_list[l]=w;
   1.514 +				      level.set(w, l);
   1.515 +				    }
   1.516 +				  }
   1.517 +	    
   1.518 +				  OutEdgeIt f;
   1.519 +				  for(g->first(f,v); g->valid(f); g->next(f)) {
   1.520 +				    if ( 0 >= (*flow)[f] ) continue;
   1.521 +				    Node w=g->head(f);
   1.522 +				    if ( level[w] == n && w != s ) {
   1.523 +				      bfs_queue.push(w);
   1.524 +				      Node first=level_list[l];
   1.525 +				      if ( g->valid(first) ) left.set(first,w);
   1.526 +				      right.set(w,first);
   1.527 +				      level_list[l]=w;
   1.528 +				      level.set(w, l);
   1.529 +				    }
   1.530 +				  }
   1.531 +				}
   1.532 +	  
   1.533 +	  
   1.534 +				//the starting flow
   1.535 +				OutEdgeIt e;
   1.536 +				for(g->first(e,s); g->valid(e); g->next(e)) 
   1.537 +				  {
   1.538 +				    Num rem=(*capacity)[e]-(*flow)[e];
   1.539 +				    if ( rem <= 0 ) continue;
   1.540 +				    Node w=g->head(e);
   1.541 +				    if ( level[w] < n ) {	  
   1.542 +				      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   1.543 +				      flow->set(e, (*capacity)[e]); 
   1.544 +				      excess.set(w, excess[w]+rem);
   1.545 +				    }
   1.546 +				  }
   1.547 +	  
   1.548 +				InEdgeIt f;
   1.549 +				for(g->first(f,s); g->valid(f); g->next(f)) 
   1.550 +				  {
   1.551 +				    if ( (*flow)[f] <= 0 ) continue;
   1.552 +				    Node w=g->tail(f);
   1.553 +				    if ( level[w] < n ) {	  
   1.554 +				      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   1.555 +				      excess.set(w, excess[w]+(*flow)[f]);
   1.556 +				      flow->set(f, 0); 
   1.557 +				    }
   1.558 +				  }  
   1.559 +				break;
   1.560 +			      } //case PRE_FLOW
   1.561 +			    }
   1.562 +			  } //preflowPreproc
   1.563 +
   1.564 +
   1.565 +
   1.566 +    void relabel(Node w, int newlevel, VecStack& active,  
   1.567 +		 VecNode& level_list, NNMap& left, 
   1.568 +		 NNMap& right, int& b, int& k, bool what_heur ) 
   1.569 +    {
   1.570 +
   1.571 +      Num lev=level[w];	
   1.572 +      
   1.573 +      Node right_n=right[w];
   1.574 +      Node left_n=left[w];
   1.575 +      
   1.576 +      //unlacing starts
   1.577 +      if ( g->valid(right_n) ) {
   1.578 +	if ( g->valid(left_n) ) {
   1.579 +	  right.set(left_n, right_n);
   1.580 +	  left.set(right_n, left_n);
   1.581 +	} else {
   1.582 +	  level_list[lev]=right_n;   
   1.583 +	  left.set(right_n, INVALID);
   1.584 +	} 
   1.585 +      } else {
   1.586 +	if ( g->valid(left_n) ) {
   1.587 +	  right.set(left_n, INVALID);
   1.588 +	} else { 
   1.589 +	  level_list[lev]=INVALID;   
   1.590 +	} 
   1.591 +      } 
   1.592 +      //unlacing ends
   1.593 +		
   1.594 +      if ( !g->valid(level_list[lev]) ) {
   1.595 +	      
   1.596 +	//gapping starts
   1.597 +	for (int i=lev; i!=k ; ) {
   1.598 +	  Node v=level_list[++i];
   1.599 +	  while ( g->valid(v) ) {
   1.600 +	    level.set(v,n);
   1.601 +	    v=right[v];
   1.602 +	  }
   1.603 +	  level_list[i]=INVALID;
   1.604 +	  if ( !what_heur ) {
   1.605 +	    while ( !active[i].empty() ) {
   1.606 +	      active[i].pop();    //FIXME: ezt szebben kene
   1.607 +	    }
   1.608 +	  }	     
   1.609 +	}
   1.610 +	
   1.611 +	level.set(w,n);
   1.612 +	b=lev-1;
   1.613 +	k=b;
   1.614 +	//gapping ends
   1.615 +	
   1.616 +      } else {
   1.617 +	
   1.618 +	if ( newlevel == n ) level.set(w,n); 
   1.619 +	else {
   1.620 +	  level.set(w,++newlevel);
   1.621 +	  active[newlevel].push(w);
   1.622 +	  if ( what_heur ) b=newlevel;
   1.623 +	  if ( k < newlevel ) ++k;      //now k=newlevel
   1.624 +	  Node first=level_list[newlevel];
   1.625 +	  if ( g->valid(first) ) left.set(first,w);
   1.626 +	  right.set(w,first);
   1.627 +	  left.set(w,INVALID);
   1.628 +	  level_list[newlevel]=w;
   1.629 +	}
   1.630 +      }
   1.631 +      
   1.632 +    } //relabel
   1.633 +
   1.634 +
   1.635 +    template<typename MapGraphWrapper> 
   1.636 +    class DistanceMap {
   1.637 +    protected:
   1.638 +      const MapGraphWrapper* g;
   1.639 +      typename MapGraphWrapper::template NodeMap<int> dist; 
   1.640 +    public:
   1.641 +      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
   1.642 +      void set(const typename MapGraphWrapper::Node& n, int a) { 
   1.643 +	dist.set(n, a); 
   1.644 +      }
   1.645 +      int operator[](const typename MapGraphWrapper::Node& n) 
   1.646 +      { return dist[n]; }
   1.647 +      //       int get(const typename MapGraphWrapper::Node& n) const { 
   1.648 +      // 	return dist[n]; }
   1.649 +      //       bool get(const typename MapGraphWrapper::Edge& e) const { 
   1.650 +      // 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
   1.651 +      bool operator[](const typename MapGraphWrapper::Edge& e) const { 
   1.652 +	return (dist[g->tail(e)]<dist[g->head(e)]); 
   1.653 +      }
   1.654 +    };
   1.655 +    
   1.656 +  };
   1.657 +
   1.658 +
   1.659 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   1.660 +  void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1( flowEnum fe ) 
   1.661 +  {
   1.662 +      
   1.663 +    int heur0=(int)(H0*n);  //time while running 'bound decrease' 
   1.664 +    int heur1=(int)(H1*n);  //time while running 'highest label'
   1.665 +    int heur=heur1;         //starting time interval (#of relabels)
   1.666 +    int numrelabel=0;
   1.667 +     
   1.668 +    bool what_heur=1;       
   1.669 +    //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   1.670 +
   1.671 +    bool end=false;     
   1.672 +    //Needed for 'bound decrease', true means no active nodes are above bound b.
   1.673 +
   1.674 +    int k=n-2;  //bound on the highest level under n containing a node
   1.675 +    int b=k;    //bound on the highest level under n of an active node
   1.676 +      
   1.677 +    VecStack active(n);
   1.678 +      
   1.679 +    NNMap left(*g, INVALID);
   1.680 +    NNMap right(*g, INVALID);
   1.681 +    VecNode level_list(n,INVALID);
   1.682 +    //List of the nodes in level i<n, set to n.
   1.683 +
   1.684 +    NodeIt v;
   1.685 +    for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   1.686 +    //setting each node to level n
   1.687 +      
   1.688 +    switch ( fe ) {
   1.689 +    case PRE_FLOW:
   1.690 +      {
   1.691 +	//counting the excess
   1.692 +	NodeIt v;
   1.693 +	for(g->first(v); g->valid(v); g->next(v)) {
   1.694 +	  Num exc=0;
   1.695 +	  
   1.696 +	  InEdgeIt e;
   1.697 +	  for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
   1.698 +	  OutEdgeIt f;
   1.699 +	  for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
   1.700 +	    
   1.701 +	  excess.set(v,exc);	  
   1.702 +	    
   1.703 +	  //putting the active nodes into the stack
   1.704 +	  int lev=level[v];
   1.705 +	  if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
   1.706 +	}
   1.707 +	break;
   1.708 +      }
   1.709 +    case GEN_FLOW:
   1.710 +      {
   1.711 +	//Counting the excess of t
   1.712 +	Num exc=0;
   1.713 +	  
   1.714 +	InEdgeIt e;
   1.715 +	for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   1.716 +	OutEdgeIt f;
   1.717 +	for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
   1.718 +	  
   1.719 +	excess.set(t,exc);	
   1.720 +	  
   1.721 +	break;
   1.722 +      }
   1.723 +    default:
   1.724 +      break;
   1.725 +    }
   1.726 +      
   1.727 +    preflowPreproc( fe, active, level_list, left, right );
   1.728 +    //End of preprocessing 
   1.729 +      
   1.730 +      
   1.731 +    //Push/relabel on the highest level active nodes.
   1.732 +    while ( true ) {
   1.733 +      if ( b == 0 ) {
   1.734 +	if ( !what_heur && !end && k > 0 ) {
   1.735 +	  b=k;
   1.736 +	  end=true;
   1.737 +	} else break;
   1.738 +      }
   1.739 +	
   1.740 +      if ( active[b].empty() ) --b; 
   1.741 +      else {
   1.742 +	end=false;  
   1.743 +	Node w=active[b].top();
   1.744 +	active[b].pop();
   1.745 +	int newlevel=push(w,active);
   1.746 +	if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 
   1.747 +				     left, right, b, k, what_heur);
   1.748 +	  
   1.749 +	++numrelabel; 
   1.750 +	if ( numrelabel >= heur ) {
   1.751 +	  numrelabel=0;
   1.752 +	  if ( what_heur ) {
   1.753 +	    what_heur=0;
   1.754 +	    heur=heur0;
   1.755 +	    end=false;
   1.756 +	  } else {
   1.757 +	    what_heur=1;
   1.758 +	    heur=heur1;
   1.759 +	    b=k; 
   1.760 +	  }
   1.761 +	}
   1.762 +      } 
   1.763 +    } 
   1.764 +  }
   1.765 +
   1.766 +
   1.767 +
   1.768 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   1.769 +  void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase2() 
   1.770 +  {
   1.771 +      
   1.772 +    int k=n-2;  //bound on the highest level under n containing a node
   1.773 +    int b=k;    //bound on the highest level under n of an active node
   1.774 +      
   1.775 +    VecStack active(n);
   1.776 +    level.set(s,0);
   1.777 +    std::queue<Node> bfs_queue;
   1.778 +    bfs_queue.push(s);
   1.779 +	    
   1.780 +    while (!bfs_queue.empty()) {
   1.781 +	
   1.782 +      Node v=bfs_queue.front();	
   1.783 +      bfs_queue.pop();
   1.784 +      int l=level[v]+1;
   1.785 +	      
   1.786 +      InEdgeIt e;
   1.787 +      for(g->first(e,v); g->valid(e); g->next(e)) {
   1.788 +	if ( (*capacity)[e] <= (*flow)[e] ) continue;
   1.789 +	Node u=g->tail(e);
   1.790 +	if ( level[u] >= n ) { 
   1.791 +	  bfs_queue.push(u);
   1.792 +	  level.set(u, l);
   1.793 +	  if ( excess[u] > 0 ) active[l].push(u);
   1.794 +	}
   1.795 +      }
   1.796 +	
   1.797 +      OutEdgeIt f;
   1.798 +      for(g->first(f,v); g->valid(f); g->next(f)) {
   1.799 +	if ( 0 >= (*flow)[f] ) continue;
   1.800 +	Node u=g->head(f);
   1.801 +	if ( level[u] >= n ) { 
   1.802 +	  bfs_queue.push(u);
   1.803 +	  level.set(u, l);
   1.804 +	  if ( excess[u] > 0 ) active[l].push(u);
   1.805 +	}
   1.806 +      }
   1.807 +    }
   1.808 +    b=n-2;
   1.809 +
   1.810 +    while ( true ) {
   1.811 +	
   1.812 +      if ( b == 0 ) break;
   1.813 +
   1.814 +      if ( active[b].empty() ) --b; 
   1.815 +      else {
   1.816 +	Node w=active[b].top();
   1.817 +	active[b].pop();
   1.818 +	int newlevel=push(w,active);	  
   1.819 +
   1.820 +	//relabel
   1.821 +	if ( excess[w] > 0 ) {
   1.822 +	  level.set(w,++newlevel);
   1.823 +	  active[newlevel].push(w);
   1.824 +	  b=newlevel;
   1.825 +	}
   1.826 +      }  // if stack[b] is nonempty
   1.827 +    } // while(true)
   1.828 +  }
   1.829 +
   1.830 +
   1.831 +
   1.832 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   1.833 +  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() 
   1.834 +  {
   1.835 +    ResGW res_graph(*g, *capacity, *flow);
   1.836 +    bool _augment=false;
   1.837 +      
   1.838 +    //ReachedMap level(res_graph);
   1.839 +    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   1.840 +    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   1.841 +    bfs.pushAndSetReached(s);
   1.842 +	
   1.843 +    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); 
   1.844 +    pred.set(s, INVALID);
   1.845 +      
   1.846 +    typename ResGW::template NodeMap<Num> free(res_graph);
   1.847 +	
   1.848 +    //searching for augmenting path
   1.849 +    while ( !bfs.finished() ) { 
   1.850 +      ResGWOutEdgeIt e=bfs;
   1.851 +      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   1.852 +	Node v=res_graph.tail(e);
   1.853 +	Node w=res_graph.head(e);
   1.854 +	pred.set(w, e);
   1.855 +	if (res_graph.valid(pred[v])) {
   1.856 +	  free.set(w, std::min(free[v], res_graph.resCap(e)));
   1.857 +	} else {
   1.858 +	  free.set(w, res_graph.resCap(e)); 
   1.859 +	}
   1.860 +	if (res_graph.head(e)==t) { _augment=true; break; }
   1.861 +      }
   1.862 +	
   1.863 +      ++bfs;
   1.864 +    } //end of searching augmenting path
   1.865 +
   1.866 +    if (_augment) {
   1.867 +      Node n=t;
   1.868 +      Num augment_value=free[t];
   1.869 +      while (res_graph.valid(pred[n])) { 
   1.870 +	ResGWEdge e=pred[n];
   1.871 +	res_graph.augment(e, augment_value); 
   1.872 +	n=res_graph.tail(e);
   1.873 +      }
   1.874 +    }
   1.875 +
   1.876 +    return _augment;
   1.877 +  }
   1.878 +
   1.879 +
   1.880 +
   1.881 +
   1.882 +
   1.883 +
   1.884 +
   1.885 +
   1.886 +
   1.887 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   1.888 +  template<typename MutableGraph> 
   1.889 +  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() 
   1.890 +  {      
   1.891 +    typedef MutableGraph MG;
   1.892 +    bool _augment=false;
   1.893 +
   1.894 +    ResGW res_graph(*g, *capacity, *flow);
   1.895 +
   1.896 +    //bfs for distances on the residual graph
   1.897 +    //ReachedMap level(res_graph);
   1.898 +    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   1.899 +    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   1.900 +    bfs.pushAndSetReached(s);
   1.901 +    typename ResGW::template NodeMap<int> 
   1.902 +      dist(res_graph); //filled up with 0's
   1.903 +
   1.904 +    //F will contain the physical copy of the residual graph
   1.905 +    //with the set of edges which are on shortest paths
   1.906 +    MG F;
   1.907 +    typename ResGW::template NodeMap<typename MG::Node> 
   1.908 +      res_graph_to_F(res_graph);
   1.909 +    {
   1.910 +      typename ResGW::NodeIt n;
   1.911 +      for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
   1.912 +	res_graph_to_F.set(n, F.addNode());
   1.913 +      }
   1.914 +    }
   1.915 +
   1.916 +    typename MG::Node sF=res_graph_to_F[s];
   1.917 +    typename MG::Node tF=res_graph_to_F[t];
   1.918 +    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
   1.919 +    typename MG::template EdgeMap<Num> residual_capacity(F);
   1.920 +
   1.921 +    while ( !bfs.finished() ) { 
   1.922 +      ResGWOutEdgeIt e=bfs;
   1.923 +      if (res_graph.valid(e)) {
   1.924 +	if (bfs.isBNodeNewlyReached()) {
   1.925 +	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   1.926 +	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   1.927 +	  original_edge.update();
   1.928 +	  original_edge.set(f, e);
   1.929 +	  residual_capacity.update();
   1.930 +	  residual_capacity.set(f, res_graph.resCap(e));
   1.931 +	} else {
   1.932 +	  if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
   1.933 +	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   1.934 +	    original_edge.update();
   1.935 +	    original_edge.set(f, e);
   1.936 +	    residual_capacity.update();
   1.937 +	    residual_capacity.set(f, res_graph.resCap(e));
   1.938 +	  }
   1.939 +	}
   1.940 +      }
   1.941 +      ++bfs;
   1.942 +    } //computing distances from s in the residual graph
   1.943 +
   1.944 +    bool __augment=true;
   1.945 +
   1.946 +    while (__augment) {
   1.947 +      __augment=false;
   1.948 +      //computing blocking flow with dfs
   1.949 +      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
   1.950 +      typename MG::template NodeMap<typename MG::Edge> pred(F);
   1.951 +      pred.set(sF, INVALID);
   1.952 +      //invalid iterators for sources
   1.953 +
   1.954 +      typename MG::template NodeMap<Num> free(F);
   1.955 +
   1.956 +      dfs.pushAndSetReached(sF);      
   1.957 +      while (!dfs.finished()) {
   1.958 +	++dfs;
   1.959 +	if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
   1.960 +	  if (dfs.isBNodeNewlyReached()) {
   1.961 +	    typename MG::Node v=F.aNode(dfs);
   1.962 +	    typename MG::Node w=F.bNode(dfs);
   1.963 +	    pred.set(w, dfs);
   1.964 +	    if (F.valid(pred[v])) {
   1.965 +	      free.set(w, std::min(free[v], residual_capacity[dfs]));
   1.966 +	    } else {
   1.967 +	      free.set(w, residual_capacity[dfs]); 
   1.968 +	    }
   1.969 +	    if (w==tF) { 
   1.970 +	      __augment=true; 
   1.971 +	      _augment=true;
   1.972 +	      break; 
   1.973 +	    }
   1.974 +	      
   1.975 +	  } else {
   1.976 +	    F.erase(/*typename MG::OutEdgeIt*/(dfs));
   1.977 +	  }
   1.978 +	} 
   1.979 +      }
   1.980 +
   1.981 +      if (__augment) {
   1.982 +	typename MG::Node n=tF;
   1.983 +	Num augment_value=free[tF];
   1.984 +	while (F.valid(pred[n])) { 
   1.985 +	  typename MG::Edge e=pred[n];
   1.986 +	  res_graph.augment(original_edge[e], augment_value); 
   1.987 +	  n=F.tail(e);
   1.988 +	  if (residual_capacity[e]==augment_value) 
   1.989 +	    F.erase(e); 
   1.990 +	  else 
   1.991 +	    residual_capacity.set(e, residual_capacity[e]-augment_value);
   1.992 +	}
   1.993 +      }
   1.994 +	
   1.995 +    }
   1.996 +            
   1.997 +    return _augment;
   1.998 +  }
   1.999 +
  1.1000 +
  1.1001 +
  1.1002 +
  1.1003 +
  1.1004 +
  1.1005 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  1.1006 +  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() 
  1.1007 +  {
  1.1008 +    bool _augment=false;
  1.1009 +
  1.1010 +    ResGW res_graph(*g, *capacity, *flow);
  1.1011 +      
  1.1012 +    //ReachedMap level(res_graph);
  1.1013 +    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
  1.1014 +    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
  1.1015 +
  1.1016 +    bfs.pushAndSetReached(s);
  1.1017 +    DistanceMap<ResGW> dist(res_graph);
  1.1018 +    while ( !bfs.finished() ) { 
  1.1019 +      ResGWOutEdgeIt e=bfs;
  1.1020 +      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
  1.1021 +	dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
  1.1022 +      }
  1.1023 +      ++bfs;
  1.1024 +    } //computing distances from s in the residual graph
  1.1025 +
  1.1026 +      //Subgraph containing the edges on some shortest paths
  1.1027 +    ConstMap<typename ResGW::Node, bool> true_map(true);
  1.1028 +    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, 
  1.1029 +      DistanceMap<ResGW> > FilterResGW;
  1.1030 +    FilterResGW filter_res_graph(res_graph, true_map, dist);
  1.1031 +
  1.1032 +    //Subgraph, which is able to delete edges which are already 
  1.1033 +    //met by the dfs
  1.1034 +    typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> 
  1.1035 +      first_out_edges(filter_res_graph);
  1.1036 +    typename FilterResGW::NodeIt v;
  1.1037 +    for(filter_res_graph.first(v); filter_res_graph.valid(v); 
  1.1038 +	filter_res_graph.next(v)) 
  1.1039 +      {
  1.1040 + 	typename FilterResGW::OutEdgeIt e;
  1.1041 + 	filter_res_graph.first(e, v);
  1.1042 + 	first_out_edges.set(v, e);
  1.1043 +      }
  1.1044 +    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
  1.1045 +      template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
  1.1046 +    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
  1.1047 +
  1.1048 +    bool __augment=true;
  1.1049 +
  1.1050 +    while (__augment) {
  1.1051 +
  1.1052 +      __augment=false;
  1.1053 +      //computing blocking flow with dfs
  1.1054 +      DfsIterator< ErasingResGW, 
  1.1055 +	typename ErasingResGW::template NodeMap<bool> > 
  1.1056 +	dfs(erasing_res_graph);
  1.1057 +      typename ErasingResGW::
  1.1058 +	template NodeMap<typename ErasingResGW::OutEdgeIt> 
  1.1059 +	pred(erasing_res_graph); 
  1.1060 +      pred.set(s, INVALID);
  1.1061 +      //invalid iterators for sources
  1.1062 +
  1.1063 +      typename ErasingResGW::template NodeMap<Num> 
  1.1064 +	free1(erasing_res_graph);
  1.1065 +
  1.1066 +      dfs.pushAndSetReached(
  1.1067 +			    typename ErasingResGW::Node(
  1.1068 +							typename FilterResGW::Node(
  1.1069 +										   typename ResGW::Node(s)
  1.1070 +										   )
  1.1071 +							)
  1.1072 +			    );
  1.1073 +      while (!dfs.finished()) {
  1.1074 +	++dfs;
  1.1075 +	if (erasing_res_graph.valid(
  1.1076 +				    typename ErasingResGW::OutEdgeIt(dfs))) 
  1.1077 + 	  { 
  1.1078 +  	    if (dfs.isBNodeNewlyReached()) {
  1.1079 +	  
  1.1080 + 	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
  1.1081 + 	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
  1.1082 +
  1.1083 + 	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
  1.1084 + 	      if (erasing_res_graph.valid(pred[v])) {
  1.1085 + 		free1.set(w, std::min(free1[v], res_graph.resCap(
  1.1086 +								 typename ErasingResGW::OutEdgeIt(dfs))));
  1.1087 + 	      } else {
  1.1088 + 		free1.set(w, res_graph.resCap(
  1.1089 +					      typename ErasingResGW::OutEdgeIt(dfs))); 
  1.1090 + 	      }
  1.1091 +	      
  1.1092 + 	      if (w==t) { 
  1.1093 + 		__augment=true; 
  1.1094 + 		_augment=true;
  1.1095 + 		break; 
  1.1096 + 	      }
  1.1097 + 	    } else {
  1.1098 + 	      erasing_res_graph.erase(dfs);
  1.1099 +	    }
  1.1100 +	  }
  1.1101 +      }	
  1.1102 +
  1.1103 +      if (__augment) {
  1.1104 +	typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
  1.1105 +	// 	  typename ResGW::NodeMap<Num> a(res_graph);
  1.1106 +	// 	  typename ResGW::Node b;
  1.1107 +	// 	  Num j=a[b];
  1.1108 +	// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
  1.1109 +	// 	  typename FilterResGW::Node b1;
  1.1110 +	// 	  Num j1=a1[b1];
  1.1111 +	// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
  1.1112 +	// 	  typename ErasingResGW::Node b2;
  1.1113 +	// 	  Num j2=a2[b2];
  1.1114 +	Num augment_value=free1[n];
  1.1115 +	while (erasing_res_graph.valid(pred[n])) { 
  1.1116 +	  typename ErasingResGW::OutEdgeIt e=pred[n];
  1.1117 +	  res_graph.augment(e, augment_value);
  1.1118 +	  n=erasing_res_graph.tail(e);
  1.1119 +	  if (res_graph.resCap(e)==0)
  1.1120 +	    erasing_res_graph.erase(e);
  1.1121 +	}
  1.1122 +      }
  1.1123 +      
  1.1124 +    } //while (__augment) 
  1.1125 +            
  1.1126 +    return _augment;
  1.1127 +  }
  1.1128 +
  1.1129 +
  1.1130 +
  1.1131 +  /// @}
  1.1132 +  
  1.1133 +} //END OF NAMESPACE HUGO
  1.1134 +
  1.1135 +#endif //HUGO_MAX_FLOW_H
  1.1136 +
  1.1137 +
  1.1138 +
  1.1139 +