max_flow.h (wich doesn't use STL
authoralpar
Thu, 22 Jul 2004 14:19:23 +0000
changeset 726835ebe1b3250
parent 725 9fa4045571cd
child 727 aada518af30f
max_flow.h (wich doesn't use STL
- max_flow.h (which doesn't use STL stack) moved to /src/hugo
- for_each_macros.h was removed from max_flow.h
- (blocking) flow augmenting stuffs was removed.
src/hugo/max_flow.h
src/work/jacint/max_flow_no_stack.h
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/hugo/max_flow.h	Thu Jul 22 14:19:23 2004 +0000
     1.3 @@ -0,0 +1,892 @@
     1.4 +// -*- C++ -*-
     1.5 +#ifndef HUGO_MAX_FLOW_NO_STACK_H
     1.6 +#define HUGO_MAX_FLOW_NO_STACK_H
     1.7 +
     1.8 +#include <vector>
     1.9 +#include <queue>
    1.10 +//#include <stack>
    1.11 +
    1.12 +#include <hugo/graph_wrapper.h>
    1.13 +#include <hugo/invalid.h>
    1.14 +#include <hugo/maps.h>
    1.15 +
    1.16 +/// \file
    1.17 +/// \brief The same as max_flow.h, but without using stl stack for the active nodes. Only for test.
    1.18 +/// \ingroup galgs
    1.19 +
    1.20 +namespace hugo {
    1.21 +
    1.22 +  /// \addtogroup galgs
    1.23 +  /// @{                                                                                                                                        
    1.24 +  ///Maximum flow algorithms class.
    1.25 +
    1.26 +  ///This class provides various algorithms for finding a flow of
    1.27 +  ///maximum value in a directed graph. The \e source node, the \e
    1.28 +  ///target node, the \e capacity of the edges and the \e starting \e
    1.29 +  ///flow value of the edges should be passed to the algorithm through the
    1.30 +  ///constructor. It is possible to change these quantities using the
    1.31 +  ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
    1.32 +  ///\ref resetFlow. Before any subsequent runs of any algorithm of
    1.33 +  ///the class \ref resetFlow should be called. 
    1.34 +
    1.35 +  ///After running an algorithm of the class, the actual flow value 
    1.36 +  ///can be obtained by calling \ref flowValue(). The minimum
    1.37 +  ///value cut can be written into a \c node map of \c bools by
    1.38 +  ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
    1.39 +  ///the inclusionwise minimum and maximum of the minimum value
    1.40 +  ///cuts, resp.)                                                                                                                               
    1.41 +  ///\param Graph The directed graph type the algorithm runs on.
    1.42 +  ///\param Num The number type of the capacities and the flow values.
    1.43 +  ///\param CapMap The capacity map type.
    1.44 +  ///\param FlowMap The flow map type.                                                                                                           
    1.45 +  ///\author Marton Makai, Jacint Szabo 
    1.46 +  template <typename Graph, typename Num,
    1.47 +	    typename CapMap=typename Graph::template EdgeMap<Num>,
    1.48 +            typename FlowMap=typename Graph::template EdgeMap<Num> >
    1.49 +  class MaxFlow {
    1.50 +  protected:
    1.51 +    typedef typename Graph::Node Node;
    1.52 +    typedef typename Graph::NodeIt NodeIt;
    1.53 +    typedef typename Graph::EdgeIt EdgeIt;
    1.54 +    typedef typename Graph::OutEdgeIt OutEdgeIt;
    1.55 +    typedef typename Graph::InEdgeIt InEdgeIt;
    1.56 +
    1.57 +    //    typedef typename std::vector<std::stack<Node> > VecStack;
    1.58 +    typedef typename std::vector<Node> VecFirst;
    1.59 +    typedef typename Graph::template NodeMap<Node> NNMap;
    1.60 +    typedef typename std::vector<Node> VecNode;
    1.61 +
    1.62 +    const Graph* g;
    1.63 +    Node s;
    1.64 +    Node t;
    1.65 +    const CapMap* capacity;
    1.66 +    FlowMap* flow;
    1.67 +    int n;      //the number of nodes of G
    1.68 +    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    1.69 +    //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    1.70 +    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    1.71 +    typedef typename ResGW::Edge ResGWEdge;
    1.72 +    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    1.73 +    typedef typename Graph::template NodeMap<int> ReachedMap;
    1.74 +
    1.75 +
    1.76 +    //level works as a bool map in augmenting path algorithms and is
    1.77 +    //used by bfs for storing reached information.  In preflow, it
    1.78 +    //shows the levels of nodes.     
    1.79 +    ReachedMap level;
    1.80 +
    1.81 +    //excess is needed only in preflow
    1.82 +    typename Graph::template NodeMap<Num> excess;
    1.83 +
    1.84 +    // constants used for heuristics
    1.85 +    static const int H0=20;
    1.86 +    static const int H1=1;
    1.87 +
    1.88 +  public:
    1.89 +
    1.90 +    ///Indicates the property of the starting flow.
    1.91 +
    1.92 +    ///Indicates the property of the starting flow. The meanings are as follows:
    1.93 +    ///- \c ZERO_FLOW: constant zero flow
    1.94 +    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
    1.95 +    ///the sum of the out-flows in every node except the \e source and
    1.96 +    ///the \e target.
    1.97 +    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at 
    1.98 +    ///least the sum of the out-flows in every node except the \e source.
    1.99 +    ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be 
   1.100 +    ///set to the constant zero flow in the beginning of the algorithm in this case.
   1.101 +    enum FlowEnum{
   1.102 +      ZERO_FLOW,
   1.103 +      GEN_FLOW,
   1.104 +      PRE_FLOW,
   1.105 +      NO_FLOW
   1.106 +    };
   1.107 +
   1.108 +    enum StatusEnum {
   1.109 +      AFTER_NOTHING,
   1.110 +      AFTER_AUGMENTING,
   1.111 +      AFTER_FAST_AUGMENTING, 
   1.112 +      AFTER_PRE_FLOW_PHASE_1,      
   1.113 +      AFTER_PRE_FLOW_PHASE_2
   1.114 +    };
   1.115 +
   1.116 +    /// Don not needle this flag only if necessary.
   1.117 +    StatusEnum status;
   1.118 +
   1.119 +//     int number_of_augmentations;
   1.120 +
   1.121 +
   1.122 +//     template<typename IntMap>
   1.123 +//     class TrickyReachedMap {
   1.124 +//     protected:
   1.125 +//       IntMap* map;
   1.126 +//       int* number_of_augmentations;
   1.127 +//     public:
   1.128 +//       TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : 
   1.129 +// 	map(&_map), number_of_augmentations(&_number_of_augmentations) { }
   1.130 +//       void set(const Node& n, bool b) {
   1.131 +// 	if (b)
   1.132 +// 	  map->set(n, *number_of_augmentations);
   1.133 +// 	else 
   1.134 +// 	  map->set(n, *number_of_augmentations-1);
   1.135 +//       }
   1.136 +//       bool operator[](const Node& n) const { 
   1.137 +// 	return (*map)[n]==*number_of_augmentations; 
   1.138 +//       }
   1.139 +//     };
   1.140 +    
   1.141 +    ///Constructor
   1.142 +
   1.143 +    ///\todo Document, please.
   1.144 +    ///
   1.145 +    MaxFlow(const Graph& _G, Node _s, Node _t,
   1.146 +		   const CapMap& _capacity, FlowMap& _flow) :
   1.147 +      g(&_G), s(_s), t(_t), capacity(&_capacity),
   1.148 +      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0), 
   1.149 +      status(AFTER_NOTHING) { }
   1.150 +
   1.151 +    ///Runs a maximum flow algorithm.
   1.152 +
   1.153 +    ///Runs a preflow algorithm, which is the fastest maximum flow
   1.154 +    ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
   1.155 +    ///\pre The starting flow must be
   1.156 +    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   1.157 +    /// - an arbitary flow if \c fe is \c GEN_FLOW,
   1.158 +    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
   1.159 +    /// - any map if \c fe is NO_FLOW.
   1.160 +    void run(FlowEnum fe=ZERO_FLOW) {
   1.161 +      preflow(fe);
   1.162 +    }
   1.163 +
   1.164 +                                                                              
   1.165 +    ///Runs a preflow algorithm.  
   1.166 +
   1.167 +    ///Runs a preflow algorithm. The preflow algorithms provide the
   1.168 +    ///fastest way to compute a maximum flow in a directed graph.
   1.169 +    ///\pre The starting flow must be
   1.170 +    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   1.171 +    /// - an arbitary flow if \c fe is \c GEN_FLOW,
   1.172 +    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
   1.173 +    /// - any map if \c fe is NO_FLOW.
   1.174 +    ///
   1.175 +    ///\todo NO_FLOW should be the default flow.
   1.176 +    void preflow(FlowEnum fe) {
   1.177 +      preflowPhase1(fe);
   1.178 +      preflowPhase2();
   1.179 +    }
   1.180 +    // Heuristics:
   1.181 +    //   2 phase
   1.182 +    //   gap
   1.183 +    //   list 'level_list' on the nodes on level i implemented by hand
   1.184 +    //   stack 'active' on the active nodes on level i                                                                                    
   1.185 +    //   runs heuristic 'highest label' for H1*n relabels
   1.186 +    //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
   1.187 +    //   Parameters H0 and H1 are initialized to 20 and 1.
   1.188 +
   1.189 +    ///Runs the first phase of the preflow algorithm.
   1.190 +
   1.191 +    ///The preflow algorithm consists of two phases, this method runs the
   1.192 +    ///first phase. After the first phase the maximum flow value and a
   1.193 +    ///minimum value cut can already be computed, though a maximum flow
   1.194 +    ///is net yet obtained. So after calling this method \ref flowValue
   1.195 +    ///and \ref actMinCut gives proper results.
   1.196 +    ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
   1.197 +    ///give minimum value cuts unless calling \ref preflowPhase2.
   1.198 +    ///\pre The starting flow must be
   1.199 +    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   1.200 +    /// - an arbitary flow if \c fe is \c GEN_FLOW,
   1.201 +    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
   1.202 +    /// - any map if \c fe is NO_FLOW.
   1.203 +    void preflowPhase1(FlowEnum fe)
   1.204 +    {
   1.205 +
   1.206 +      int heur0=(int)(H0*n);  //time while running 'bound decrease'
   1.207 +      int heur1=(int)(H1*n);  //time while running 'highest label'
   1.208 +      int heur=heur1;         //starting time interval (#of relabels)
   1.209 +      int numrelabel=0;
   1.210 +
   1.211 +      bool what_heur=1;
   1.212 +      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   1.213 +
   1.214 +      bool end=false;
   1.215 +      //Needed for 'bound decrease', true means no active nodes are above bound
   1.216 +      //b.
   1.217 +
   1.218 +      int k=n-2;  //bound on the highest level under n containing a node
   1.219 +      int b=k;    //bound on the highest level under n of an active node
   1.220 +
   1.221 +      VecFirst first(n, INVALID);
   1.222 +      NNMap next(*g, INVALID); //maybe INVALID is not needed
   1.223 +      //    VecStack active(n);
   1.224 +
   1.225 +      NNMap left(*g, INVALID);
   1.226 +      NNMap right(*g, INVALID);
   1.227 +      VecNode level_list(n,INVALID);
   1.228 +      //List of the nodes in level i<n, set to n.
   1.229 +
   1.230 +      NodeIt v;
   1.231 +      for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   1.232 +      //setting each node to level n
   1.233 +
   1.234 +      if ( fe == NO_FLOW ) {
   1.235 +	EdgeIt e;
   1.236 +	for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
   1.237 +      }
   1.238 +
   1.239 +      switch (fe) { //computing the excess
   1.240 +      case PRE_FLOW:
   1.241 +	{
   1.242 +	  NodeIt v;
   1.243 +	  for(g->first(v); g->valid(v); g->next(v)) {
   1.244 +	    Num exc=0;
   1.245 +
   1.246 +	    InEdgeIt e;
   1.247 +	    for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
   1.248 +	    OutEdgeIt f;
   1.249 +	    for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
   1.250 +
   1.251 +	    excess.set(v,exc);
   1.252 +
   1.253 +	    //putting the active nodes into the stack
   1.254 +	    int lev=level[v];
   1.255 +	    if ( exc > 0 && lev < n && v != t ) 
   1.256 +	      {
   1.257 +		next.set(v,first[lev]);
   1.258 +		first[lev]=v;
   1.259 +	      }
   1.260 +	    //	  active[lev].push(v);
   1.261 +	  }
   1.262 +	  break;
   1.263 +	}
   1.264 +      case GEN_FLOW:
   1.265 +	{
   1.266 +	  NodeIt v;
   1.267 +	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
   1.268 +
   1.269 +	  Num exc=0;
   1.270 +	  InEdgeIt e;
   1.271 +	  for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   1.272 +	  OutEdgeIt f;
   1.273 +	  for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
   1.274 +	  excess.set(t,exc);
   1.275 +	  break;
   1.276 +	}
   1.277 +      case ZERO_FLOW:
   1.278 +      case NO_FLOW:
   1.279 +	{
   1.280 +	  NodeIt v;
   1.281 +	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
   1.282 +	  break;
   1.283 +	}
   1.284 +      }
   1.285 +
   1.286 +      preflowPreproc(fe, next, first,/*active*/ level_list, left, right);
   1.287 +      //End of preprocessing
   1.288 +
   1.289 +
   1.290 +      //Push/relabel on the highest level active nodes.
   1.291 +      while ( true ) {
   1.292 +	if ( b == 0 ) {
   1.293 +	  if ( !what_heur && !end && k > 0 ) {
   1.294 +	    b=k;
   1.295 +	    end=true;
   1.296 +	  } else break;
   1.297 +	}
   1.298 +
   1.299 +	if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
   1.300 +	else {
   1.301 +	  end=false;
   1.302 +	  Node w=first[b];
   1.303 +	  first[b]=next[w];
   1.304 +	  /*	Node w=active[b].top();
   1.305 +		active[b].pop();*/
   1.306 +	  int newlevel=push(w,/*active*/next, first);
   1.307 +	  if ( excess[w] > 0 ) relabel(w, newlevel, /*active*/next, first, level_list,
   1.308 +				       left, right, b, k, what_heur);
   1.309 +
   1.310 +	  ++numrelabel;
   1.311 +	  if ( numrelabel >= heur ) {
   1.312 +	    numrelabel=0;
   1.313 +	    if ( what_heur ) {
   1.314 +	      what_heur=0;
   1.315 +	      heur=heur0;
   1.316 +	      end=false;
   1.317 +	    } else {
   1.318 +	      what_heur=1;
   1.319 +	      heur=heur1;
   1.320 +	      b=k;
   1.321 +	    }
   1.322 +	  }
   1.323 +	}
   1.324 +      }
   1.325 +
   1.326 +      status=AFTER_PRE_FLOW_PHASE_1;
   1.327 +    }
   1.328 +
   1.329 +
   1.330 +    ///Runs the second phase of the preflow algorithm.
   1.331 +
   1.332 +    ///The preflow algorithm consists of two phases, this method runs
   1.333 +    ///the second phase. After calling \ref preflowPhase1 and then
   1.334 +    ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
   1.335 +    ///\ref minMinCut and \ref maxMinCut give proper results.
   1.336 +    ///\pre \ref preflowPhase1 must be called before.
   1.337 +    void preflowPhase2()
   1.338 +    {
   1.339 +
   1.340 +      int k=n-2;  //bound on the highest level under n containing a node
   1.341 +      int b=k;    //bound on the highest level under n of an active node
   1.342 +
   1.343 +    
   1.344 +      VecFirst first(n, INVALID);
   1.345 +      NNMap next(*g, INVALID); //maybe INVALID is not needed
   1.346 +      //    VecStack active(n);
   1.347 +      level.set(s,0);
   1.348 +      std::queue<Node> bfs_queue;
   1.349 +      bfs_queue.push(s);
   1.350 +
   1.351 +      while (!bfs_queue.empty()) {
   1.352 +
   1.353 +	Node v=bfs_queue.front();
   1.354 +	bfs_queue.pop();
   1.355 +	int l=level[v]+1;
   1.356 +
   1.357 +	InEdgeIt e;
   1.358 +	for(g->first(e,v); g->valid(e); g->next(e)) {
   1.359 +	  if ( (*capacity)[e] <= (*flow)[e] ) continue;
   1.360 +	  Node u=g->tail(e);
   1.361 +	  if ( level[u] >= n ) {
   1.362 +	    bfs_queue.push(u);
   1.363 +	    level.set(u, l);
   1.364 +	    if ( excess[u] > 0 ) {
   1.365 +	      next.set(u,first[l]);
   1.366 +	      first[l]=u;
   1.367 +	      //active[l].push(u);
   1.368 +	    }
   1.369 +	  }
   1.370 +	}
   1.371 +
   1.372 +	OutEdgeIt f;
   1.373 +	for(g->first(f,v); g->valid(f); g->next(f)) {
   1.374 +	  if ( 0 >= (*flow)[f] ) continue;
   1.375 +	  Node u=g->head(f);
   1.376 +	  if ( level[u] >= n ) {
   1.377 +	    bfs_queue.push(u);
   1.378 +	    level.set(u, l);
   1.379 +	    if ( excess[u] > 0 ) {
   1.380 +	      next.set(u,first[l]);
   1.381 +	      first[l]=u;
   1.382 +	      //active[l].push(u);
   1.383 +	    }
   1.384 +	  }
   1.385 +	}
   1.386 +      }
   1.387 +      b=n-2;
   1.388 +
   1.389 +      while ( true ) {
   1.390 +
   1.391 +	if ( b == 0 ) break;
   1.392 +
   1.393 +	if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
   1.394 +	else {
   1.395 +
   1.396 +	  Node w=first[b];
   1.397 +	  first[b]=next[w];
   1.398 +	  /*	Node w=active[b].top();
   1.399 +		active[b].pop();*/
   1.400 +	  int newlevel=push(w,next, first/*active*/);
   1.401 +
   1.402 +	  //relabel
   1.403 +	  if ( excess[w] > 0 ) {
   1.404 +	    level.set(w,++newlevel);
   1.405 +	    next.set(w,first[newlevel]);
   1.406 +	    first[newlevel]=w;
   1.407 +	    //active[newlevel].push(w);
   1.408 +	    b=newlevel;
   1.409 +	  }
   1.410 +	}  // if stack[b] is nonempty
   1.411 +      } // while(true)
   1.412 +
   1.413 +      status=AFTER_PRE_FLOW_PHASE_2;
   1.414 +    }
   1.415 +
   1.416 +
   1.417 +    /// Returns the maximum value of a flow.
   1.418 +
   1.419 +    /// Returns the maximum value of a flow, by counting the 
   1.420 +    /// over-flow of the target node \ref t.
   1.421 +    /// It can be called already after running \ref preflowPhase1.
   1.422 +    Num flowValue() const {
   1.423 +      Num a=0;
   1.424 +      for(InEdgeIt e(*g,t);g->valid(e);G.next(e)) a+=(*flow)[e];
   1.425 +      for(OutEdgeIt e(*g,t);g->valid(e);G.next(e)) a-=(*flow)[e];
   1.426 +
   1.427 +      //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
   1.428 +    }
   1.429 +
   1.430 +    ///Returns a minimum value cut after calling \ref preflowPhase1.
   1.431 +
   1.432 +    ///After the first phase of the preflow algorithm the maximum flow
   1.433 +    ///value and a minimum value cut can already be computed. This
   1.434 +    ///method can be called after running \ref preflowPhase1 for
   1.435 +    ///obtaining a minimum value cut.
   1.436 +    /// \warning Gives proper result only right after calling \ref
   1.437 +    /// preflowPhase1.
   1.438 +    /// \todo We have to make some status variable which shows the
   1.439 +    /// actual state
   1.440 +    /// of the class. This enables us to determine which methods are valid
   1.441 +    /// for MinCut computation
   1.442 +    template<typename _CutMap>
   1.443 +    void actMinCut(_CutMap& M) const {
   1.444 +      NodeIt v;
   1.445 +      switch (status) {
   1.446 +      case AFTER_PRE_FLOW_PHASE_1:
   1.447 +	for(g->first(v); g->valid(v); g->next(v)) {
   1.448 +	  if (level[v] < n) {
   1.449 +	    M.set(v, false);
   1.450 +	  } else {
   1.451 +	    M.set(v, true);
   1.452 +	  }
   1.453 +	}
   1.454 +	break;
   1.455 +      case AFTER_PRE_FLOW_PHASE_2:
   1.456 +      case AFTER_NOTHING:
   1.457 +	minMinCut(M);
   1.458 +	break;
   1.459 +      case AFTER_AUGMENTING:
   1.460 +	for(g->first(v); g->valid(v); g->next(v)) {
   1.461 +	  if (level[v]) {
   1.462 +	    M.set(v, true);
   1.463 +	  } else {
   1.464 +	    M.set(v, false);
   1.465 +	  }
   1.466 +	}
   1.467 +	break;
   1.468 +      case AFTER_FAST_AUGMENTING:
   1.469 +	for(g->first(v); g->valid(v); g->next(v)) {
   1.470 +	  if (level[v]==number_of_augmentations) {
   1.471 +	    M.set(v, true);
   1.472 +	  } else {
   1.473 +	    M.set(v, false);
   1.474 +	  }
   1.475 +	}
   1.476 +	break;
   1.477 +      }
   1.478 +    }
   1.479 +
   1.480 +    ///Returns the inclusionwise minimum of the minimum value cuts.
   1.481 +
   1.482 +    ///Sets \c M to the characteristic vector of the minimum value cut
   1.483 +    ///which is inclusionwise minimum. It is computed by processing
   1.484 +    ///a bfs from the source node \c s in the residual graph.
   1.485 +    ///\pre M should be a node map of bools initialized to false.
   1.486 +    ///\pre \c flow must be a maximum flow.
   1.487 +    template<typename _CutMap>
   1.488 +    void minMinCut(_CutMap& M) const {
   1.489 +      std::queue<Node> queue;
   1.490 +
   1.491 +      M.set(s,true);
   1.492 +      queue.push(s);
   1.493 +
   1.494 +      while (!queue.empty()) {
   1.495 +        Node w=queue.front();
   1.496 +	queue.pop();
   1.497 +
   1.498 +	OutEdgeIt e;
   1.499 +	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   1.500 +	  Node v=g->head(e);
   1.501 +	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   1.502 +	    queue.push(v);
   1.503 +	    M.set(v, true);
   1.504 +	  }
   1.505 +	}
   1.506 +
   1.507 +	InEdgeIt f;
   1.508 +	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   1.509 +	  Node v=g->tail(f);
   1.510 +	  if (!M[v] && (*flow)[f] > 0 ) {
   1.511 +	    queue.push(v);
   1.512 +	    M.set(v, true);
   1.513 +	  }
   1.514 +	}
   1.515 +      }
   1.516 +    }
   1.517 +
   1.518 +    ///Returns the inclusionwise maximum of the minimum value cuts.
   1.519 +
   1.520 +    ///Sets \c M to the characteristic vector of the minimum value cut
   1.521 +    ///which is inclusionwise maximum. It is computed by processing a
   1.522 +    ///backward bfs from the target node \c t in the residual graph.
   1.523 +    ///\pre M should be a node map of bools initialized to false.
   1.524 +    ///\pre \c flow must be a maximum flow. 
   1.525 +    template<typename _CutMap>
   1.526 +    void maxMinCut(_CutMap& M) const {
   1.527 +
   1.528 +      NodeIt v;
   1.529 +      for(g->first(v) ; g->valid(v); g->next(v)) {
   1.530 +	M.set(v, true);
   1.531 +      }
   1.532 +
   1.533 +      std::queue<Node> queue;
   1.534 +
   1.535 +      M.set(t,false);
   1.536 +      queue.push(t);
   1.537 +
   1.538 +      while (!queue.empty()) {
   1.539 +        Node w=queue.front();
   1.540 +	queue.pop();
   1.541 +
   1.542 +	InEdgeIt e;
   1.543 +	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   1.544 +	  Node v=g->tail(e);
   1.545 +	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   1.546 +	    queue.push(v);
   1.547 +	    M.set(v, false);
   1.548 +	  }
   1.549 +	}
   1.550 +
   1.551 +	OutEdgeIt f;
   1.552 +	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   1.553 +	  Node v=g->head(f);
   1.554 +	  if (M[v] && (*flow)[f] > 0 ) {
   1.555 +	    queue.push(v);
   1.556 +	    M.set(v, false);
   1.557 +	  }
   1.558 +	}
   1.559 +      }
   1.560 +    }
   1.561 +
   1.562 +    ///Returns a minimum value cut.
   1.563 +
   1.564 +    ///Sets \c M to the characteristic vector of a minimum value cut.
   1.565 +    ///\pre M should be a node map of bools initialized to false.
   1.566 +    ///\pre \c flow must be a maximum flow.    
   1.567 +    template<typename CutMap>
   1.568 +    void minCut(CutMap& M) const { minMinCut(M); }
   1.569 +
   1.570 +    ///Resets the source node to \c _s.
   1.571 +
   1.572 +    ///Resets the source node to \c _s.
   1.573 +    /// 
   1.574 +    void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
   1.575 +
   1.576 +    ///Resets the target node to \c _t.
   1.577 +
   1.578 +    ///Resets the target node to \c _t.
   1.579 +    ///
   1.580 +    void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
   1.581 +
   1.582 +    /// Resets the edge map of the capacities to _cap.
   1.583 +
   1.584 +    /// Resets the edge map of the capacities to _cap.
   1.585 +    /// 
   1.586 +    void resetCap(const CapMap& _cap)
   1.587 +    { capacity=&_cap; status=AFTER_NOTHING; }
   1.588 +
   1.589 +    /// Resets the edge map of the flows to _flow.
   1.590 +
   1.591 +    /// Resets the edge map of the flows to _flow.
   1.592 +    /// 
   1.593 +    void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
   1.594 +
   1.595 +
   1.596 +  private:
   1.597 +
   1.598 +    int push(Node w, NNMap& next, VecFirst& first) {
   1.599 +
   1.600 +      int lev=level[w];
   1.601 +      Num exc=excess[w];
   1.602 +      int newlevel=n;       //bound on the next level of w
   1.603 +
   1.604 +      OutEdgeIt e;
   1.605 +      for(g->first(e,w); g->valid(e); g->next(e)) {
   1.606 +
   1.607 +	if ( (*flow)[e] >= (*capacity)[e] ) continue;
   1.608 +	Node v=g->head(e);
   1.609 +
   1.610 +	if( lev > level[v] ) { //Push is allowed now
   1.611 +
   1.612 +	  if ( excess[v]<=0 && v!=t && v!=s ) {
   1.613 +	    next.set(v,first[level[v]]);
   1.614 +	    first[level[v]]=v;
   1.615 +	    //	    int lev_v=level[v];
   1.616 +	    //active[lev_v].push(v);
   1.617 +	  }
   1.618 +
   1.619 +	  Num cap=(*capacity)[e];
   1.620 +	  Num flo=(*flow)[e];
   1.621 +	  Num remcap=cap-flo;
   1.622 +
   1.623 +	  if ( remcap >= exc ) { //A nonsaturating push.
   1.624 +
   1.625 +	    flow->set(e, flo+exc);
   1.626 +	    excess.set(v, excess[v]+exc);
   1.627 +	    exc=0;
   1.628 +	    break;
   1.629 +
   1.630 +	  } else { //A saturating push.
   1.631 +	    flow->set(e, cap);
   1.632 +	    excess.set(v, excess[v]+remcap);
   1.633 +	    exc-=remcap;
   1.634 +	  }
   1.635 +	} else if ( newlevel > level[v] ) newlevel = level[v];
   1.636 +      } //for out edges wv
   1.637 +
   1.638 +      if ( exc > 0 ) {
   1.639 +	InEdgeIt e;
   1.640 +	for(g->first(e,w); g->valid(e); g->next(e)) {
   1.641 +
   1.642 +	  if( (*flow)[e] <= 0 ) continue;
   1.643 +	  Node v=g->tail(e);
   1.644 +
   1.645 +	  if( lev > level[v] ) { //Push is allowed now
   1.646 +
   1.647 +	    if ( excess[v]<=0 && v!=t && v!=s ) {
   1.648 +	      next.set(v,first[level[v]]);
   1.649 +	      first[level[v]]=v;
   1.650 +	      //int lev_v=level[v];
   1.651 +	      //active[lev_v].push(v);
   1.652 +	    }
   1.653 +
   1.654 +	    Num flo=(*flow)[e];
   1.655 +
   1.656 +	    if ( flo >= exc ) { //A nonsaturating push.
   1.657 +
   1.658 +	      flow->set(e, flo-exc);
   1.659 +	      excess.set(v, excess[v]+exc);
   1.660 +	      exc=0;
   1.661 +	      break;
   1.662 +	    } else {  //A saturating push.
   1.663 +
   1.664 +	      excess.set(v, excess[v]+flo);
   1.665 +	      exc-=flo;
   1.666 +	      flow->set(e,0);
   1.667 +	    }
   1.668 +	  } else if ( newlevel > level[v] ) newlevel = level[v];
   1.669 +	} //for in edges vw
   1.670 +
   1.671 +      } // if w still has excess after the out edge for cycle
   1.672 +
   1.673 +      excess.set(w, exc);
   1.674 +
   1.675 +      return newlevel;
   1.676 +    }
   1.677 +
   1.678 +
   1.679 +    void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
   1.680 +			VecNode& level_list, NNMap& left, NNMap& right)
   1.681 +    {
   1.682 +      std::queue<Node> bfs_queue;
   1.683 +
   1.684 +      switch (fe) {
   1.685 +      case NO_FLOW:   //flow is already set to const zero in this case
   1.686 +      case ZERO_FLOW:
   1.687 +	{
   1.688 +	  //Reverse_bfs from t, to find the starting level.
   1.689 +	  level.set(t,0);
   1.690 +	  bfs_queue.push(t);
   1.691 +
   1.692 +	  while (!bfs_queue.empty()) {
   1.693 +
   1.694 +	    Node v=bfs_queue.front();
   1.695 +	    bfs_queue.pop();
   1.696 +	    int l=level[v]+1;
   1.697 +
   1.698 +	    InEdgeIt e;
   1.699 +	    for(g->first(e,v); g->valid(e); g->next(e)) {
   1.700 +	      Node w=g->tail(e);
   1.701 +	      if ( level[w] == n && w != s ) {
   1.702 +		bfs_queue.push(w);
   1.703 +		Node z=level_list[l];
   1.704 +		if ( g->valid(z) ) left.set(z,w);
   1.705 +		right.set(w,z);
   1.706 +		level_list[l]=w;
   1.707 +		level.set(w, l);
   1.708 +	      }
   1.709 +	    }
   1.710 +	  }
   1.711 +
   1.712 +	  //the starting flow
   1.713 +	  OutEdgeIt e;
   1.714 +	  for(g->first(e,s); g->valid(e); g->next(e))
   1.715 +	    {
   1.716 +	      Num c=(*capacity)[e];
   1.717 +	      if ( c <= 0 ) continue;
   1.718 +	      Node w=g->head(e);
   1.719 +	      if ( level[w] < n ) {
   1.720 +		if ( excess[w] <= 0 && w!=t ) 
   1.721 +		  {
   1.722 +		    next.set(w,first[level[w]]);
   1.723 +		    first[level[w]]=w;
   1.724 +		    //active[level[w]].push(w);
   1.725 +		  }
   1.726 +		flow->set(e, c);
   1.727 +		excess.set(w, excess[w]+c);
   1.728 +	      }
   1.729 +	    }
   1.730 +	  break;
   1.731 +	}
   1.732 +
   1.733 +      case GEN_FLOW:
   1.734 +      case PRE_FLOW:
   1.735 +	{
   1.736 +	  //Reverse_bfs from t in the residual graph,
   1.737 +	  //to find the starting level.
   1.738 +	  level.set(t,0);
   1.739 +	  bfs_queue.push(t);
   1.740 +
   1.741 +	  while (!bfs_queue.empty()) {
   1.742 +
   1.743 +	    Node v=bfs_queue.front();
   1.744 +	    bfs_queue.pop();
   1.745 +	    int l=level[v]+1;
   1.746 +
   1.747 +	    InEdgeIt e;
   1.748 +	    for(g->first(e,v); g->valid(e); g->next(e)) {
   1.749 +	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
   1.750 +	      Node w=g->tail(e);
   1.751 +	      if ( level[w] == n && w != s ) {
   1.752 +		bfs_queue.push(w);
   1.753 +		Node z=level_list[l];
   1.754 +		if ( g->valid(z) ) left.set(z,w);
   1.755 +		right.set(w,z);
   1.756 +		level_list[l]=w;
   1.757 +		level.set(w, l);
   1.758 +	      }
   1.759 +	    }
   1.760 +
   1.761 +	    OutEdgeIt f;
   1.762 +	    for(g->first(f,v); g->valid(f); g->next(f)) {
   1.763 +	      if ( 0 >= (*flow)[f] ) continue;
   1.764 +	      Node w=g->head(f);
   1.765 +	      if ( level[w] == n && w != s ) {
   1.766 +		bfs_queue.push(w);
   1.767 +		Node z=level_list[l];
   1.768 +		if ( g->valid(z) ) left.set(z,w);
   1.769 +		right.set(w,z);
   1.770 +		level_list[l]=w;
   1.771 +		level.set(w, l);
   1.772 +	      }
   1.773 +	    }
   1.774 +	  }
   1.775 +
   1.776 +
   1.777 +	  //the starting flow
   1.778 +	  OutEdgeIt e;
   1.779 +	  for(g->first(e,s); g->valid(e); g->next(e))
   1.780 +	    {
   1.781 +	      Num rem=(*capacity)[e]-(*flow)[e];
   1.782 +	      if ( rem <= 0 ) continue;
   1.783 +	      Node w=g->head(e);
   1.784 +	      if ( level[w] < n ) {
   1.785 +		if ( excess[w] <= 0 && w!=t )
   1.786 +		  {
   1.787 +		    next.set(w,first[level[w]]);
   1.788 +		    first[level[w]]=w;
   1.789 +		    //active[level[w]].push(w);
   1.790 +		  }   
   1.791 +		flow->set(e, (*capacity)[e]);
   1.792 +		excess.set(w, excess[w]+rem);
   1.793 +	      }
   1.794 +	    }
   1.795 +
   1.796 +	  InEdgeIt f;
   1.797 +	  for(g->first(f,s); g->valid(f); g->next(f))
   1.798 +	    {
   1.799 +	      if ( (*flow)[f] <= 0 ) continue;
   1.800 +	      Node w=g->tail(f);
   1.801 +	      if ( level[w] < n ) {
   1.802 +		if ( excess[w] <= 0 && w!=t )
   1.803 +		  {
   1.804 +		    next.set(w,first[level[w]]);
   1.805 +		    first[level[w]]=w;
   1.806 +		    //active[level[w]].push(w);
   1.807 +		  }   
   1.808 +		excess.set(w, excess[w]+(*flow)[f]);
   1.809 +		flow->set(f, 0);
   1.810 +	      }
   1.811 +	    }
   1.812 +	  break;
   1.813 +	} //case PRE_FLOW
   1.814 +      }
   1.815 +    } //preflowPreproc
   1.816 +
   1.817 +
   1.818 +
   1.819 +    void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
   1.820 +		 VecNode& level_list, NNMap& left,
   1.821 +		 NNMap& right, int& b, int& k, bool what_heur )
   1.822 +    {
   1.823 +
   1.824 +      Num lev=level[w];
   1.825 +
   1.826 +      Node right_n=right[w];
   1.827 +      Node left_n=left[w];
   1.828 +
   1.829 +      //unlacing starts
   1.830 +      if ( g->valid(right_n) ) {
   1.831 +	if ( g->valid(left_n) ) {
   1.832 +	  right.set(left_n, right_n);
   1.833 +	  left.set(right_n, left_n);
   1.834 +	} else {
   1.835 +	  level_list[lev]=right_n;
   1.836 +	  left.set(right_n, INVALID);
   1.837 +	}
   1.838 +      } else {
   1.839 +	if ( g->valid(left_n) ) {
   1.840 +	  right.set(left_n, INVALID);
   1.841 +	} else {
   1.842 +	  level_list[lev]=INVALID;
   1.843 +	}
   1.844 +      }
   1.845 +      //unlacing ends
   1.846 +
   1.847 +      if ( !g->valid(level_list[lev]) ) {
   1.848 +
   1.849 +	//gapping starts
   1.850 +	for (int i=lev; i!=k ; ) {
   1.851 +	  Node v=level_list[++i];
   1.852 +	  while ( g->valid(v) ) {
   1.853 +	    level.set(v,n);
   1.854 +	    v=right[v];
   1.855 +	  }
   1.856 +	  level_list[i]=INVALID;
   1.857 +	  if ( !what_heur ) first[i]=INVALID;
   1.858 +	  /*{
   1.859 +	    while ( !active[i].empty() ) {
   1.860 +	    active[i].pop();    //FIXME: ezt szebben kene
   1.861 +	    }
   1.862 +	    }*/
   1.863 +	}
   1.864 +
   1.865 +	level.set(w,n);
   1.866 +	b=lev-1;
   1.867 +	k=b;
   1.868 +	//gapping ends
   1.869 +
   1.870 +      } else {
   1.871 +
   1.872 +	if ( newlevel == n ) level.set(w,n);
   1.873 +	else {
   1.874 +	  level.set(w,++newlevel);
   1.875 +	  next.set(w,first[newlevel]);
   1.876 +	  first[newlevel]=w;
   1.877 +	  //	  active[newlevel].push(w);
   1.878 +	  if ( what_heur ) b=newlevel;
   1.879 +	  if ( k < newlevel ) ++k;      //now k=newlevel
   1.880 +	  Node z=level_list[newlevel];
   1.881 +	  if ( g->valid(z) ) left.set(z,w);
   1.882 +	  right.set(w,z);
   1.883 +	  left.set(w,INVALID);
   1.884 +	  level_list[newlevel]=w;
   1.885 +	}
   1.886 +      }
   1.887 +    } //relabel
   1.888 +  };  //class MaxFlow
   1.889 +} //namespace hugo
   1.890 +
   1.891 +#endif //HUGO_MAX_FLOW_H
   1.892 +
   1.893 +
   1.894 +
   1.895 +
     2.1 --- a/src/work/jacint/max_flow_no_stack.h	Thu Jul 22 14:09:21 2004 +0000
     2.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.3 @@ -1,1318 +0,0 @@
     2.4 -// -*- C++ -*-
     2.5 -#ifndef HUGO_MAX_FLOW_NO_STACK_H
     2.6 -#define HUGO_MAX_FLOW_NO_STACK_H
     2.7 -
     2.8 -#include <vector>
     2.9 -#include <queue>
    2.10 -//#include <stack>
    2.11 -
    2.12 -#include <hugo/graph_wrapper.h>
    2.13 -#include <bfs_dfs.h>
    2.14 -#include <hugo/invalid.h>
    2.15 -#include <hugo/maps.h>
    2.16 -#include <hugo/for_each_macros.h>
    2.17 -
    2.18 -/// \file
    2.19 -/// \brief The same as max_flow.h, but without using stl stack for the active nodes. Only for test.
    2.20 -/// \ingroup galgs
    2.21 -
    2.22 -namespace hugo {
    2.23 -
    2.24 -  /// \addtogroup galgs
    2.25 -  /// @{                                                                                                                                        
    2.26 -  ///Maximum flow algorithms class.
    2.27 -
    2.28 -  ///This class provides various algorithms for finding a flow of
    2.29 -  ///maximum value in a directed graph. The \e source node, the \e
    2.30 -  ///target node, the \e capacity of the edges and the \e starting \e
    2.31 -  ///flow value of the edges should be passed to the algorithm through the
    2.32 -  ///constructor. It is possible to change these quantities using the
    2.33 -  ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
    2.34 -  ///\ref resetFlow. Before any subsequent runs of any algorithm of
    2.35 -  ///the class \ref resetFlow should be called. 
    2.36 -
    2.37 -  ///After running an algorithm of the class, the actual flow value 
    2.38 -  ///can be obtained by calling \ref flowValue(). The minimum
    2.39 -  ///value cut can be written into a \c node map of \c bools by
    2.40 -  ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
    2.41 -  ///the inclusionwise minimum and maximum of the minimum value
    2.42 -  ///cuts, resp.)                                                                                                                               
    2.43 -  ///\param Graph The directed graph type the algorithm runs on.
    2.44 -  ///\param Num The number type of the capacities and the flow values.
    2.45 -  ///\param CapMap The capacity map type.
    2.46 -  ///\param FlowMap The flow map type.                                                                                                           
    2.47 -  ///\author Marton Makai, Jacint Szabo 
    2.48 -  template <typename Graph, typename Num,
    2.49 -	    typename CapMap=typename Graph::template EdgeMap<Num>,
    2.50 -            typename FlowMap=typename Graph::template EdgeMap<Num> >
    2.51 -  class MaxFlowNoStack {
    2.52 -  protected:
    2.53 -    typedef typename Graph::Node Node;
    2.54 -    typedef typename Graph::NodeIt NodeIt;
    2.55 -    typedef typename Graph::EdgeIt EdgeIt;
    2.56 -    typedef typename Graph::OutEdgeIt OutEdgeIt;
    2.57 -    typedef typename Graph::InEdgeIt InEdgeIt;
    2.58 -
    2.59 -    //    typedef typename std::vector<std::stack<Node> > VecStack;
    2.60 -    typedef typename std::vector<Node> VecFirst;
    2.61 -    typedef typename Graph::template NodeMap<Node> NNMap;
    2.62 -    typedef typename std::vector<Node> VecNode;
    2.63 -
    2.64 -    const Graph* g;
    2.65 -    Node s;
    2.66 -    Node t;
    2.67 -    const CapMap* capacity;
    2.68 -    FlowMap* flow;
    2.69 -    int n;      //the number of nodes of G
    2.70 -    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    2.71 -    //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    2.72 -    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    2.73 -    typedef typename ResGW::Edge ResGWEdge;
    2.74 -    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    2.75 -    typedef typename Graph::template NodeMap<int> ReachedMap;
    2.76 -
    2.77 -
    2.78 -    //level works as a bool map in augmenting path algorithms and is
    2.79 -    //used by bfs for storing reached information.  In preflow, it
    2.80 -    //shows the levels of nodes.     
    2.81 -    ReachedMap level;
    2.82 -
    2.83 -    //excess is needed only in preflow
    2.84 -    typename Graph::template NodeMap<Num> excess;
    2.85 -
    2.86 -    //fixme    
    2.87 -//   protected:
    2.88 -    //     MaxFlow() { }
    2.89 -    //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
    2.90 -    // 	     FlowMap& _flow)
    2.91 -    //       {
    2.92 -    // 	g=&_G;
    2.93 -    // 	s=_s;
    2.94 -    // 	t=_t;
    2.95 -    // 	capacity=&_capacity;
    2.96 -    // 	flow=&_flow;
    2.97 -    // 	n=_G.nodeNum;
    2.98 -    // 	level.set (_G); //kellene vmi ilyesmi fv
    2.99 -    // 	excess(_G,0); //itt is
   2.100 -    //       }
   2.101 -
   2.102 -    // constants used for heuristics
   2.103 -    static const int H0=20;
   2.104 -    static const int H1=1;
   2.105 -
   2.106 -  public:
   2.107 -
   2.108 -    ///Indicates the property of the starting flow.
   2.109 -
   2.110 -    ///Indicates the property of the starting flow. The meanings are as follows:
   2.111 -    ///- \c ZERO_FLOW: constant zero flow
   2.112 -    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
   2.113 -    ///the sum of the out-flows in every node except the \e source and
   2.114 -    ///the \e target.
   2.115 -    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at 
   2.116 -    ///least the sum of the out-flows in every node except the \e source.
   2.117 -    ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be 
   2.118 -    ///set to the constant zero flow in the beginning of the algorithm in this case.
   2.119 -    enum FlowEnum{
   2.120 -      ZERO_FLOW,
   2.121 -      GEN_FLOW,
   2.122 -      PRE_FLOW,
   2.123 -      NO_FLOW
   2.124 -    };
   2.125 -
   2.126 -    enum StatusEnum {
   2.127 -      AFTER_NOTHING,
   2.128 -      AFTER_AUGMENTING,
   2.129 -      AFTER_FAST_AUGMENTING, 
   2.130 -      AFTER_PRE_FLOW_PHASE_1,      
   2.131 -      AFTER_PRE_FLOW_PHASE_2
   2.132 -    };
   2.133 -
   2.134 -    /// Don not needle this flag only if necessary.
   2.135 -    StatusEnum status;
   2.136 -    int number_of_augmentations;
   2.137 -
   2.138 -
   2.139 -    template<typename IntMap>
   2.140 -    class TrickyReachedMap {
   2.141 -    protected:
   2.142 -      IntMap* map;
   2.143 -      int* number_of_augmentations;
   2.144 -    public:
   2.145 -      TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : 
   2.146 -	map(&_map), number_of_augmentations(&_number_of_augmentations) { }
   2.147 -      void set(const Node& n, bool b) {
   2.148 -	if (b)
   2.149 -	  map->set(n, *number_of_augmentations);
   2.150 -	else 
   2.151 -	  map->set(n, *number_of_augmentations-1);
   2.152 -      }
   2.153 -      bool operator[](const Node& n) const { 
   2.154 -	return (*map)[n]==*number_of_augmentations; 
   2.155 -      }
   2.156 -    };
   2.157 -    
   2.158 -    ///Constructor
   2.159 -
   2.160 -    ///\todo Document, please.
   2.161 -    ///
   2.162 -    MaxFlowNoStack(const Graph& _G, Node _s, Node _t,
   2.163 -		   const CapMap& _capacity, FlowMap& _flow) :
   2.164 -      g(&_G), s(_s), t(_t), capacity(&_capacity),
   2.165 -      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0), 
   2.166 -      status(AFTER_NOTHING), number_of_augmentations(0) { }
   2.167 -
   2.168 -    ///Runs a maximum flow algorithm.
   2.169 -
   2.170 -    ///Runs a preflow algorithm, which is the fastest maximum flow
   2.171 -    ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
   2.172 -    ///\pre The starting flow must be
   2.173 -    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   2.174 -    /// - an arbitary flow if \c fe is \c GEN_FLOW,
   2.175 -    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
   2.176 -    /// - any map if \c fe is NO_FLOW.
   2.177 -    void run(FlowEnum fe=ZERO_FLOW) {
   2.178 -      preflow(fe);
   2.179 -    }
   2.180 -
   2.181 -                                                                              
   2.182 -    ///Runs a preflow algorithm.  
   2.183 -
   2.184 -    ///Runs a preflow algorithm. The preflow algorithms provide the
   2.185 -    ///fastest way to compute a maximum flow in a directed graph.
   2.186 -    ///\pre The starting flow must be
   2.187 -    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   2.188 -    /// - an arbitary flow if \c fe is \c GEN_FLOW,
   2.189 -    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
   2.190 -    /// - any map if \c fe is NO_FLOW.
   2.191 -    ///
   2.192 -    ///\todo NO_FLOW should be the default flow.
   2.193 -    void preflow(FlowEnum fe) {
   2.194 -      preflowPhase1(fe);
   2.195 -      preflowPhase2();
   2.196 -    }
   2.197 -    // Heuristics:
   2.198 -    //   2 phase
   2.199 -    //   gap
   2.200 -    //   list 'level_list' on the nodes on level i implemented by hand
   2.201 -    //   stack 'active' on the active nodes on level i                                                                                    
   2.202 -    //   runs heuristic 'highest label' for H1*n relabels
   2.203 -    //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
   2.204 -    //   Parameters H0 and H1 are initialized to 20 and 1.
   2.205 -
   2.206 -    ///Runs the first phase of the preflow algorithm.
   2.207 -
   2.208 -    ///The preflow algorithm consists of two phases, this method runs the
   2.209 -    ///first phase. After the first phase the maximum flow value and a
   2.210 -    ///minimum value cut can already be computed, though a maximum flow
   2.211 -    ///is net yet obtained. So after calling this method \ref flowValue
   2.212 -    ///and \ref actMinCut gives proper results.
   2.213 -    ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
   2.214 -    ///give minimum value cuts unless calling \ref preflowPhase2.
   2.215 -    ///\pre The starting flow must be
   2.216 -    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   2.217 -    /// - an arbitary flow if \c fe is \c GEN_FLOW,
   2.218 -    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
   2.219 -    /// - any map if \c fe is NO_FLOW.
   2.220 -    void preflowPhase1(FlowEnum fe);
   2.221 -
   2.222 -    ///Runs the second phase of the preflow algorithm.
   2.223 -
   2.224 -    ///The preflow algorithm consists of two phases, this method runs
   2.225 -    ///the second phase. After calling \ref preflowPhase1 and then
   2.226 -    ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
   2.227 -    ///\ref minMinCut and \ref maxMinCut give proper results.
   2.228 -    ///\pre \ref preflowPhase1 must be called before.
   2.229 -    void preflowPhase2();
   2.230 -
   2.231 -    /// Starting from a flow, this method searches for an augmenting path
   2.232 -    /// according to the Edmonds-Karp algorithm
   2.233 -    /// and augments the flow on if any.
   2.234 -    /// The return value shows if the augmentation was succesful.
   2.235 -    bool augmentOnShortestPath();
   2.236 -    bool augmentOnShortestPath2();
   2.237 -
   2.238 -    /// Starting from a flow, this method searches for an augmenting blocking
   2.239 -    /// flow according to Dinits' algorithm and augments the flow on if any.
   2.240 -    /// The blocking flow is computed in a physically constructed
   2.241 -    /// residual graph of type \c Mutablegraph.
   2.242 -    /// The return value show sif the augmentation was succesful.
   2.243 -    template<typename MutableGraph> bool augmentOnBlockingFlow();
   2.244 -
   2.245 -    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
   2.246 -    /// residual graph is not constructed physically.
   2.247 -    /// The return value shows if the augmentation was succesful.
   2.248 -    bool augmentOnBlockingFlow2();
   2.249 -
   2.250 -    /// Returns the maximum value of a flow.
   2.251 -
   2.252 -    /// Returns the maximum value of a flow, by counting the 
   2.253 -    /// over-flow of the target node \ref t.
   2.254 -    /// It can be called already after running \ref preflowPhase1.
   2.255 -    Num flowValue() const {
   2.256 -      Num a=0;
   2.257 -      FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e];
   2.258 -      FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e];
   2.259 -      return a;
   2.260 -      //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
   2.261 -    }
   2.262 -
   2.263 -    ///Returns a minimum value cut after calling \ref preflowPhase1.
   2.264 -
   2.265 -    ///After the first phase of the preflow algorithm the maximum flow
   2.266 -    ///value and a minimum value cut can already be computed. This
   2.267 -    ///method can be called after running \ref preflowPhase1 for
   2.268 -    ///obtaining a minimum value cut.
   2.269 -    /// \warning Gives proper result only right after calling \ref
   2.270 -    /// preflowPhase1.
   2.271 -    /// \todo We have to make some status variable which shows the
   2.272 -    /// actual state
   2.273 -    /// of the class. This enables us to determine which methods are valid
   2.274 -    /// for MinCut computation
   2.275 -    template<typename _CutMap>
   2.276 -    void actMinCut(_CutMap& M) const {
   2.277 -      NodeIt v;
   2.278 -      switch (status) {
   2.279 -      case AFTER_PRE_FLOW_PHASE_1:
   2.280 -	for(g->first(v); g->valid(v); g->next(v)) {
   2.281 -	  if (level[v] < n) {
   2.282 -	    M.set(v, false);
   2.283 -	  } else {
   2.284 -	    M.set(v, true);
   2.285 -	  }
   2.286 -	}
   2.287 -	break;
   2.288 -      case AFTER_PRE_FLOW_PHASE_2:
   2.289 -      case AFTER_NOTHING:
   2.290 -	minMinCut(M);
   2.291 -	break;
   2.292 -      case AFTER_AUGMENTING:
   2.293 -	for(g->first(v); g->valid(v); g->next(v)) {
   2.294 -	  if (level[v]) {
   2.295 -	    M.set(v, true);
   2.296 -	  } else {
   2.297 -	    M.set(v, false);
   2.298 -	  }
   2.299 -	}
   2.300 -	break;
   2.301 -      case AFTER_FAST_AUGMENTING:
   2.302 -	for(g->first(v); g->valid(v); g->next(v)) {
   2.303 -	  if (level[v]==number_of_augmentations) {
   2.304 -	    M.set(v, true);
   2.305 -	  } else {
   2.306 -	    M.set(v, false);
   2.307 -	  }
   2.308 -	}
   2.309 -	break;
   2.310 -      }
   2.311 -    }
   2.312 -
   2.313 -    ///Returns the inclusionwise minimum of the minimum value cuts.
   2.314 -
   2.315 -    ///Sets \c M to the characteristic vector of the minimum value cut
   2.316 -    ///which is inclusionwise minimum. It is computed by processing
   2.317 -    ///a bfs from the source node \c s in the residual graph.
   2.318 -    ///\pre M should be a node map of bools initialized to false.
   2.319 -    ///\pre \c flow must be a maximum flow.
   2.320 -    template<typename _CutMap>
   2.321 -    void minMinCut(_CutMap& M) const {
   2.322 -      std::queue<Node> queue;
   2.323 -
   2.324 -      M.set(s,true);
   2.325 -      queue.push(s);
   2.326 -
   2.327 -      while (!queue.empty()) {
   2.328 -        Node w=queue.front();
   2.329 -	queue.pop();
   2.330 -
   2.331 -	OutEdgeIt e;
   2.332 -	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   2.333 -	  Node v=g->head(e);
   2.334 -	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   2.335 -	    queue.push(v);
   2.336 -	    M.set(v, true);
   2.337 -	  }
   2.338 -	}
   2.339 -
   2.340 -	InEdgeIt f;
   2.341 -	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   2.342 -	  Node v=g->tail(f);
   2.343 -	  if (!M[v] && (*flow)[f] > 0 ) {
   2.344 -	    queue.push(v);
   2.345 -	    M.set(v, true);
   2.346 -	  }
   2.347 -	}
   2.348 -      }
   2.349 -    }
   2.350 -
   2.351 -    ///Returns the inclusionwise maximum of the minimum value cuts.
   2.352 -
   2.353 -    ///Sets \c M to the characteristic vector of the minimum value cut
   2.354 -    ///which is inclusionwise maximum. It is computed by processing a
   2.355 -    ///backward bfs from the target node \c t in the residual graph.
   2.356 -    ///\pre M should be a node map of bools initialized to false.
   2.357 -    ///\pre \c flow must be a maximum flow. 
   2.358 -    template<typename _CutMap>
   2.359 -    void maxMinCut(_CutMap& M) const {
   2.360 -
   2.361 -      NodeIt v;
   2.362 -      for(g->first(v) ; g->valid(v); g->next(v)) {
   2.363 -	M.set(v, true);
   2.364 -      }
   2.365 -
   2.366 -      std::queue<Node> queue;
   2.367 -
   2.368 -      M.set(t,false);
   2.369 -      queue.push(t);
   2.370 -
   2.371 -      while (!queue.empty()) {
   2.372 -        Node w=queue.front();
   2.373 -	queue.pop();
   2.374 -
   2.375 -	InEdgeIt e;
   2.376 -	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   2.377 -	  Node v=g->tail(e);
   2.378 -	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   2.379 -	    queue.push(v);
   2.380 -	    M.set(v, false);
   2.381 -	  }
   2.382 -	}
   2.383 -
   2.384 -	OutEdgeIt f;
   2.385 -	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   2.386 -	  Node v=g->head(f);
   2.387 -	  if (M[v] && (*flow)[f] > 0 ) {
   2.388 -	    queue.push(v);
   2.389 -	    M.set(v, false);
   2.390 -	  }
   2.391 -	}
   2.392 -      }
   2.393 -    }
   2.394 -
   2.395 -    ///Returns a minimum value cut.
   2.396 -
   2.397 -    ///Sets \c M to the characteristic vector of a minimum value cut.
   2.398 -    ///\pre M should be a node map of bools initialized to false.
   2.399 -    ///\pre \c flow must be a maximum flow.    
   2.400 -    template<typename CutMap>
   2.401 -    void minCut(CutMap& M) const { minMinCut(M); }
   2.402 -
   2.403 -    ///Resets the source node to \c _s.
   2.404 -
   2.405 -    ///Resets the source node to \c _s.
   2.406 -    /// 
   2.407 -    void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
   2.408 -
   2.409 -    ///Resets the target node to \c _t.
   2.410 -
   2.411 -    ///Resets the target node to \c _t.
   2.412 -    ///
   2.413 -    void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
   2.414 -
   2.415 -    /// Resets the edge map of the capacities to _cap.
   2.416 -
   2.417 -    /// Resets the edge map of the capacities to _cap.
   2.418 -    /// 
   2.419 -    void resetCap(const CapMap& _cap)
   2.420 -    { capacity=&_cap; status=AFTER_NOTHING; }
   2.421 -
   2.422 -    /// Resets the edge map of the flows to _flow.
   2.423 -
   2.424 -    /// Resets the edge map of the flows to _flow.
   2.425 -    /// 
   2.426 -    void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
   2.427 -
   2.428 -
   2.429 -  private:
   2.430 -
   2.431 -    int push(Node w, NNMap& next, VecFirst& first) {
   2.432 -
   2.433 -      int lev=level[w];
   2.434 -      Num exc=excess[w];
   2.435 -      int newlevel=n;       //bound on the next level of w
   2.436 -
   2.437 -      OutEdgeIt e;
   2.438 -      for(g->first(e,w); g->valid(e); g->next(e)) {
   2.439 -
   2.440 -	if ( (*flow)[e] >= (*capacity)[e] ) continue;
   2.441 -	Node v=g->head(e);
   2.442 -
   2.443 -	if( lev > level[v] ) { //Push is allowed now
   2.444 -
   2.445 -	  if ( excess[v]<=0 && v!=t && v!=s ) {
   2.446 -	    next.set(v,first[level[v]]);
   2.447 -	    first[level[v]]=v;
   2.448 -	    //	    int lev_v=level[v];
   2.449 -	    //active[lev_v].push(v);
   2.450 -	  }
   2.451 -
   2.452 -	  Num cap=(*capacity)[e];
   2.453 -	  Num flo=(*flow)[e];
   2.454 -	  Num remcap=cap-flo;
   2.455 -
   2.456 -	  if ( remcap >= exc ) { //A nonsaturating push.
   2.457 -
   2.458 -	    flow->set(e, flo+exc);
   2.459 -	    excess.set(v, excess[v]+exc);
   2.460 -	    exc=0;
   2.461 -	    break;
   2.462 -
   2.463 -	  } else { //A saturating push.
   2.464 -	    flow->set(e, cap);
   2.465 -	    excess.set(v, excess[v]+remcap);
   2.466 -	    exc-=remcap;
   2.467 -	  }
   2.468 -	} else if ( newlevel > level[v] ) newlevel = level[v];
   2.469 -      } //for out edges wv
   2.470 -
   2.471 -      if ( exc > 0 ) {
   2.472 -	InEdgeIt e;
   2.473 -	for(g->first(e,w); g->valid(e); g->next(e)) {
   2.474 -
   2.475 -	  if( (*flow)[e] <= 0 ) continue;
   2.476 -	  Node v=g->tail(e);
   2.477 -
   2.478 -	  if( lev > level[v] ) { //Push is allowed now
   2.479 -
   2.480 -	    if ( excess[v]<=0 && v!=t && v!=s ) {
   2.481 -	      next.set(v,first[level[v]]);
   2.482 -	      first[level[v]]=v;
   2.483 -	      //int lev_v=level[v];
   2.484 -	      //active[lev_v].push(v);
   2.485 -	    }
   2.486 -
   2.487 -	    Num flo=(*flow)[e];
   2.488 -
   2.489 -	    if ( flo >= exc ) { //A nonsaturating push.
   2.490 -
   2.491 -	      flow->set(e, flo-exc);
   2.492 -	      excess.set(v, excess[v]+exc);
   2.493 -	      exc=0;
   2.494 -	      break;
   2.495 -	    } else {  //A saturating push.
   2.496 -
   2.497 -	      excess.set(v, excess[v]+flo);
   2.498 -	      exc-=flo;
   2.499 -	      flow->set(e,0);
   2.500 -	    }
   2.501 -	  } else if ( newlevel > level[v] ) newlevel = level[v];
   2.502 -	} //for in edges vw
   2.503 -
   2.504 -      } // if w still has excess after the out edge for cycle
   2.505 -
   2.506 -      excess.set(w, exc);
   2.507 -
   2.508 -      return newlevel;
   2.509 -    }
   2.510 -
   2.511 -
   2.512 -    void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
   2.513 -			VecNode& level_list, NNMap& left, NNMap& right)
   2.514 -    {
   2.515 -      std::queue<Node> bfs_queue;
   2.516 -
   2.517 -      switch (fe) {
   2.518 -      case NO_FLOW:   //flow is already set to const zero in this case
   2.519 -      case ZERO_FLOW:
   2.520 -	{
   2.521 -	  //Reverse_bfs from t, to find the starting level.
   2.522 -	  level.set(t,0);
   2.523 -	  bfs_queue.push(t);
   2.524 -
   2.525 -	  while (!bfs_queue.empty()) {
   2.526 -
   2.527 -	    Node v=bfs_queue.front();
   2.528 -	    bfs_queue.pop();
   2.529 -	    int l=level[v]+1;
   2.530 -
   2.531 -	    InEdgeIt e;
   2.532 -	    for(g->first(e,v); g->valid(e); g->next(e)) {
   2.533 -	      Node w=g->tail(e);
   2.534 -	      if ( level[w] == n && w != s ) {
   2.535 -		bfs_queue.push(w);
   2.536 -		Node z=level_list[l];
   2.537 -		if ( g->valid(z) ) left.set(z,w);
   2.538 -		right.set(w,z);
   2.539 -		level_list[l]=w;
   2.540 -		level.set(w, l);
   2.541 -	      }
   2.542 -	    }
   2.543 -	  }
   2.544 -
   2.545 -	  //the starting flow
   2.546 -	  OutEdgeIt e;
   2.547 -	  for(g->first(e,s); g->valid(e); g->next(e))
   2.548 -	    {
   2.549 -	      Num c=(*capacity)[e];
   2.550 -	      if ( c <= 0 ) continue;
   2.551 -	      Node w=g->head(e);
   2.552 -	      if ( level[w] < n ) {
   2.553 -		if ( excess[w] <= 0 && w!=t ) 
   2.554 -		  {
   2.555 -		    next.set(w,first[level[w]]);
   2.556 -		    first[level[w]]=w;
   2.557 -		    //active[level[w]].push(w);
   2.558 -		  }
   2.559 -		flow->set(e, c);
   2.560 -		excess.set(w, excess[w]+c);
   2.561 -	      }
   2.562 -	    }
   2.563 -	  break;
   2.564 -	}
   2.565 -
   2.566 -      case GEN_FLOW:
   2.567 -      case PRE_FLOW:
   2.568 -	{
   2.569 -	  //Reverse_bfs from t in the residual graph,
   2.570 -	  //to find the starting level.
   2.571 -	  level.set(t,0);
   2.572 -	  bfs_queue.push(t);
   2.573 -
   2.574 -	  while (!bfs_queue.empty()) {
   2.575 -
   2.576 -	    Node v=bfs_queue.front();
   2.577 -	    bfs_queue.pop();
   2.578 -	    int l=level[v]+1;
   2.579 -
   2.580 -	    InEdgeIt e;
   2.581 -	    for(g->first(e,v); g->valid(e); g->next(e)) {
   2.582 -	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
   2.583 -	      Node w=g->tail(e);
   2.584 -	      if ( level[w] == n && w != s ) {
   2.585 -		bfs_queue.push(w);
   2.586 -		Node z=level_list[l];
   2.587 -		if ( g->valid(z) ) left.set(z,w);
   2.588 -		right.set(w,z);
   2.589 -		level_list[l]=w;
   2.590 -		level.set(w, l);
   2.591 -	      }
   2.592 -	    }
   2.593 -
   2.594 -	    OutEdgeIt f;
   2.595 -	    for(g->first(f,v); g->valid(f); g->next(f)) {
   2.596 -	      if ( 0 >= (*flow)[f] ) continue;
   2.597 -	      Node w=g->head(f);
   2.598 -	      if ( level[w] == n && w != s ) {
   2.599 -		bfs_queue.push(w);
   2.600 -		Node z=level_list[l];
   2.601 -		if ( g->valid(z) ) left.set(z,w);
   2.602 -		right.set(w,z);
   2.603 -		level_list[l]=w;
   2.604 -		level.set(w, l);
   2.605 -	      }
   2.606 -	    }
   2.607 -	  }
   2.608 -
   2.609 -
   2.610 -	  //the starting flow
   2.611 -	  OutEdgeIt e;
   2.612 -	  for(g->first(e,s); g->valid(e); g->next(e))
   2.613 -	    {
   2.614 -	      Num rem=(*capacity)[e]-(*flow)[e];
   2.615 -	      if ( rem <= 0 ) continue;
   2.616 -	      Node w=g->head(e);
   2.617 -	      if ( level[w] < n ) {
   2.618 -		if ( excess[w] <= 0 && w!=t )
   2.619 -		  {
   2.620 -		    next.set(w,first[level[w]]);
   2.621 -		    first[level[w]]=w;
   2.622 -		    //active[level[w]].push(w);
   2.623 -		  }   
   2.624 -		flow->set(e, (*capacity)[e]);
   2.625 -		excess.set(w, excess[w]+rem);
   2.626 -	      }
   2.627 -	    }
   2.628 -
   2.629 -	  InEdgeIt f;
   2.630 -	  for(g->first(f,s); g->valid(f); g->next(f))
   2.631 -	    {
   2.632 -	      if ( (*flow)[f] <= 0 ) continue;
   2.633 -	      Node w=g->tail(f);
   2.634 -	      if ( level[w] < n ) {
   2.635 -		if ( excess[w] <= 0 && w!=t )
   2.636 -		  {
   2.637 -		    next.set(w,first[level[w]]);
   2.638 -		    first[level[w]]=w;
   2.639 -		    //active[level[w]].push(w);
   2.640 -		  }   
   2.641 -		excess.set(w, excess[w]+(*flow)[f]);
   2.642 -		flow->set(f, 0);
   2.643 -	      }
   2.644 -	    }
   2.645 -	  break;
   2.646 -	} //case PRE_FLOW
   2.647 -      }
   2.648 -    } //preflowPreproc
   2.649 -
   2.650 -
   2.651 -
   2.652 -    void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
   2.653 -		 VecNode& level_list, NNMap& left,
   2.654 -		 NNMap& right, int& b, int& k, bool what_heur )
   2.655 -    {
   2.656 -
   2.657 -      Num lev=level[w];
   2.658 -
   2.659 -      Node right_n=right[w];
   2.660 -      Node left_n=left[w];
   2.661 -
   2.662 -      //unlacing starts
   2.663 -      if ( g->valid(right_n) ) {
   2.664 -	if ( g->valid(left_n) ) {
   2.665 -	  right.set(left_n, right_n);
   2.666 -	  left.set(right_n, left_n);
   2.667 -	} else {
   2.668 -	  level_list[lev]=right_n;
   2.669 -	  left.set(right_n, INVALID);
   2.670 -	}
   2.671 -      } else {
   2.672 -	if ( g->valid(left_n) ) {
   2.673 -	  right.set(left_n, INVALID);
   2.674 -	} else {
   2.675 -	  level_list[lev]=INVALID;
   2.676 -	}
   2.677 -      }
   2.678 -      //unlacing ends
   2.679 -
   2.680 -      if ( !g->valid(level_list[lev]) ) {
   2.681 -
   2.682 -	//gapping starts
   2.683 -	for (int i=lev; i!=k ; ) {
   2.684 -	  Node v=level_list[++i];
   2.685 -	  while ( g->valid(v) ) {
   2.686 -	    level.set(v,n);
   2.687 -	    v=right[v];
   2.688 -	  }
   2.689 -	  level_list[i]=INVALID;
   2.690 -	  if ( !what_heur ) first[i]=INVALID;
   2.691 -	  /*{
   2.692 -	    while ( !active[i].empty() ) {
   2.693 -	    active[i].pop();    //FIXME: ezt szebben kene
   2.694 -	    }
   2.695 -	    }*/
   2.696 -	}
   2.697 -
   2.698 -	level.set(w,n);
   2.699 -	b=lev-1;
   2.700 -	k=b;
   2.701 -	//gapping ends
   2.702 -
   2.703 -      } else {
   2.704 -
   2.705 -	if ( newlevel == n ) level.set(w,n);
   2.706 -	else {
   2.707 -	  level.set(w,++newlevel);
   2.708 -	  next.set(w,first[newlevel]);
   2.709 -	  first[newlevel]=w;
   2.710 -	  //	  active[newlevel].push(w);
   2.711 -	  if ( what_heur ) b=newlevel;
   2.712 -	  if ( k < newlevel ) ++k;      //now k=newlevel
   2.713 -	  Node z=level_list[newlevel];
   2.714 -	  if ( g->valid(z) ) left.set(z,w);
   2.715 -	  right.set(w,z);
   2.716 -	  left.set(w,INVALID);
   2.717 -	  level_list[newlevel]=w;
   2.718 -	}
   2.719 -      }
   2.720 -
   2.721 -    } //relabel
   2.722 -
   2.723 -
   2.724 -    template<typename MapGraphWrapper>
   2.725 -    class DistanceMap {
   2.726 -    protected:
   2.727 -      const MapGraphWrapper* g;
   2.728 -      typename MapGraphWrapper::template NodeMap<int> dist;
   2.729 -    public:
   2.730 -      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
   2.731 -      void set(const typename MapGraphWrapper::Node& n, int a) {
   2.732 -	dist.set(n, a);
   2.733 -      }
   2.734 -      int operator[](const typename MapGraphWrapper::Node& n) const { 
   2.735 -	return dist[n]; 
   2.736 -      }
   2.737 -      //       int get(const typename MapGraphWrapper::Node& n) const {
   2.738 -      // 	return dist[n]; }
   2.739 -      //       bool get(const typename MapGraphWrapper::Edge& e) const {
   2.740 -      // 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
   2.741 -      bool operator[](const typename MapGraphWrapper::Edge& e) const {
   2.742 -	return (dist[g->tail(e)]<dist[g->head(e)]);
   2.743 -      }
   2.744 -    };
   2.745 -
   2.746 -  };
   2.747 -
   2.748 -
   2.749 -  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   2.750 -  void MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::preflowPhase1(FlowEnum fe)
   2.751 -  {
   2.752 -
   2.753 -    int heur0=(int)(H0*n);  //time while running 'bound decrease'
   2.754 -    int heur1=(int)(H1*n);  //time while running 'highest label'
   2.755 -    int heur=heur1;         //starting time interval (#of relabels)
   2.756 -    int numrelabel=0;
   2.757 -
   2.758 -    bool what_heur=1;
   2.759 -    //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   2.760 -
   2.761 -    bool end=false;
   2.762 -    //Needed for 'bound decrease', true means no active nodes are above bound
   2.763 -    //b.
   2.764 -
   2.765 -    int k=n-2;  //bound on the highest level under n containing a node
   2.766 -    int b=k;    //bound on the highest level under n of an active node
   2.767 -
   2.768 -    VecFirst first(n, INVALID);
   2.769 -    NNMap next(*g, INVALID); //maybe INVALID is not needed
   2.770 -    //    VecStack active(n);
   2.771 -
   2.772 -    NNMap left(*g, INVALID);
   2.773 -    NNMap right(*g, INVALID);
   2.774 -    VecNode level_list(n,INVALID);
   2.775 -    //List of the nodes in level i<n, set to n.
   2.776 -
   2.777 -    NodeIt v;
   2.778 -    for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   2.779 -    //setting each node to level n
   2.780 -
   2.781 -    if ( fe == NO_FLOW ) {
   2.782 -      EdgeIt e;
   2.783 -      for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
   2.784 -    }
   2.785 -
   2.786 -    switch (fe) { //computing the excess
   2.787 -    case PRE_FLOW:
   2.788 -      {
   2.789 -	NodeIt v;
   2.790 -	for(g->first(v); g->valid(v); g->next(v)) {
   2.791 -	  Num exc=0;
   2.792 -
   2.793 -	  InEdgeIt e;
   2.794 -	  for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
   2.795 -	  OutEdgeIt f;
   2.796 -	  for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
   2.797 -
   2.798 -	  excess.set(v,exc);
   2.799 -
   2.800 -	  //putting the active nodes into the stack
   2.801 -	  int lev=level[v];
   2.802 -	  if ( exc > 0 && lev < n && v != t ) 
   2.803 -	    {
   2.804 -	      next.set(v,first[lev]);
   2.805 -	      first[lev]=v;
   2.806 -	    }
   2.807 -	  //	  active[lev].push(v);
   2.808 -	}
   2.809 -	break;
   2.810 -      }
   2.811 -    case GEN_FLOW:
   2.812 -      {
   2.813 -	NodeIt v;
   2.814 -	for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
   2.815 -
   2.816 -	Num exc=0;
   2.817 -	InEdgeIt e;
   2.818 -	for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   2.819 -	OutEdgeIt f;
   2.820 -	for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
   2.821 -	excess.set(t,exc);
   2.822 -	break;
   2.823 -      }
   2.824 -    case ZERO_FLOW:
   2.825 -    case NO_FLOW:
   2.826 -      {
   2.827 -	NodeIt v;
   2.828 -        for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
   2.829 -	break;
   2.830 -      }
   2.831 -    }
   2.832 -
   2.833 -    preflowPreproc(fe, next, first,/*active*/ level_list, left, right);
   2.834 -    //End of preprocessing
   2.835 -
   2.836 -
   2.837 -    //Push/relabel on the highest level active nodes.
   2.838 -    while ( true ) {
   2.839 -      if ( b == 0 ) {
   2.840 -	if ( !what_heur && !end && k > 0 ) {
   2.841 -	  b=k;
   2.842 -	  end=true;
   2.843 -	} else break;
   2.844 -      }
   2.845 -
   2.846 -      if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
   2.847 -      else {
   2.848 -	end=false;
   2.849 -	Node w=first[b];
   2.850 -	first[b]=next[w];
   2.851 -	/*	Node w=active[b].top();
   2.852 -		active[b].pop();*/
   2.853 -	int newlevel=push(w,/*active*/next, first);
   2.854 -	if ( excess[w] > 0 ) relabel(w, newlevel, /*active*/next, first, level_list,
   2.855 -				     left, right, b, k, what_heur);
   2.856 -
   2.857 -	++numrelabel;
   2.858 -	if ( numrelabel >= heur ) {
   2.859 -	  numrelabel=0;
   2.860 -	  if ( what_heur ) {
   2.861 -	    what_heur=0;
   2.862 -	    heur=heur0;
   2.863 -	    end=false;
   2.864 -	  } else {
   2.865 -	    what_heur=1;
   2.866 -	    heur=heur1;
   2.867 -	    b=k;
   2.868 -	  }
   2.869 -	}
   2.870 -      }
   2.871 -    }
   2.872 -
   2.873 -    status=AFTER_PRE_FLOW_PHASE_1;
   2.874 -  }
   2.875 -
   2.876 -
   2.877 -
   2.878 -  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   2.879 -  void MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::preflowPhase2()
   2.880 -  {
   2.881 -
   2.882 -    int k=n-2;  //bound on the highest level under n containing a node
   2.883 -    int b=k;    //bound on the highest level under n of an active node
   2.884 -
   2.885 -    
   2.886 -    VecFirst first(n, INVALID);
   2.887 -    NNMap next(*g, INVALID); //maybe INVALID is not needed
   2.888 -    //    VecStack active(n);
   2.889 -    level.set(s,0);
   2.890 -    std::queue<Node> bfs_queue;
   2.891 -    bfs_queue.push(s);
   2.892 -
   2.893 -    while (!bfs_queue.empty()) {
   2.894 -
   2.895 -      Node v=bfs_queue.front();
   2.896 -      bfs_queue.pop();
   2.897 -      int l=level[v]+1;
   2.898 -
   2.899 -      InEdgeIt e;
   2.900 -      for(g->first(e,v); g->valid(e); g->next(e)) {
   2.901 -	if ( (*capacity)[e] <= (*flow)[e] ) continue;
   2.902 -	Node u=g->tail(e);
   2.903 -	if ( level[u] >= n ) {
   2.904 -	  bfs_queue.push(u);
   2.905 -	  level.set(u, l);
   2.906 -	  if ( excess[u] > 0 ) {
   2.907 -	    next.set(u,first[l]);
   2.908 -	    first[l]=u;
   2.909 -	    //active[l].push(u);
   2.910 -	  }
   2.911 -	}
   2.912 -      }
   2.913 -
   2.914 -      OutEdgeIt f;
   2.915 -      for(g->first(f,v); g->valid(f); g->next(f)) {
   2.916 -	if ( 0 >= (*flow)[f] ) continue;
   2.917 -	Node u=g->head(f);
   2.918 -	if ( level[u] >= n ) {
   2.919 -	  bfs_queue.push(u);
   2.920 -	  level.set(u, l);
   2.921 -	  if ( excess[u] > 0 ) {
   2.922 -	    next.set(u,first[l]);
   2.923 -	    first[l]=u;
   2.924 -	    //active[l].push(u);
   2.925 -	  }
   2.926 -	}
   2.927 -      }
   2.928 -    }
   2.929 -    b=n-2;
   2.930 -
   2.931 -    while ( true ) {
   2.932 -
   2.933 -      if ( b == 0 ) break;
   2.934 -
   2.935 -      if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
   2.936 -      else {
   2.937 -
   2.938 -	Node w=first[b];
   2.939 -	first[b]=next[w];
   2.940 -	/*	Node w=active[b].top();
   2.941 -		active[b].pop();*/
   2.942 -	int newlevel=push(w,next, first/*active*/);
   2.943 -
   2.944 -	//relabel
   2.945 -	if ( excess[w] > 0 ) {
   2.946 -	  level.set(w,++newlevel);
   2.947 -	  next.set(w,first[newlevel]);
   2.948 -	  first[newlevel]=w;
   2.949 -	  //active[newlevel].push(w);
   2.950 -	  b=newlevel;
   2.951 -	}
   2.952 -      }  // if stack[b] is nonempty
   2.953 -    } // while(true)
   2.954 -
   2.955 -    status=AFTER_PRE_FLOW_PHASE_2;
   2.956 -  }
   2.957 -
   2.958 -
   2.959 -
   2.960 -  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   2.961 -  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
   2.962 -  {
   2.963 -    ResGW res_graph(*g, *capacity, *flow);
   2.964 -    bool _augment=false;
   2.965 -
   2.966 -    //ReachedMap level(res_graph);
   2.967 -    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   2.968 -    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   2.969 -    bfs.pushAndSetReached(s);
   2.970 -
   2.971 -    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
   2.972 -    pred.set(s, INVALID);
   2.973 -
   2.974 -    typename ResGW::template NodeMap<Num> free(res_graph);
   2.975 -
   2.976 -    //searching for augmenting path
   2.977 -    while ( !bfs.finished() ) {
   2.978 -      ResGWOutEdgeIt e=bfs;
   2.979 -      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   2.980 -	Node v=res_graph.tail(e);
   2.981 -	Node w=res_graph.head(e);
   2.982 -	pred.set(w, e);
   2.983 -	if (res_graph.valid(pred[v])) {
   2.984 -	  free.set(w, std::min(free[v], res_graph.resCap(e)));
   2.985 -	} else {
   2.986 -	  free.set(w, res_graph.resCap(e));
   2.987 -	}
   2.988 -	if (res_graph.head(e)==t) { _augment=true; break; }
   2.989 -      }
   2.990 -
   2.991 -      ++bfs;
   2.992 -    } //end of searching augmenting path
   2.993 -
   2.994 -    if (_augment) {
   2.995 -      Node n=t;
   2.996 -      Num augment_value=free[t];
   2.997 -      while (res_graph.valid(pred[n])) {
   2.998 -	ResGWEdge e=pred[n];
   2.999 -	res_graph.augment(e, augment_value);
  2.1000 -	n=res_graph.tail(e);
  2.1001 -      }
  2.1002 -    }
  2.1003 -
  2.1004 -    status=AFTER_AUGMENTING;
  2.1005 -    return _augment;
  2.1006 -  }
  2.1007 -
  2.1008 -
  2.1009 -  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  2.1010 -  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
  2.1011 -  {
  2.1012 -    ResGW res_graph(*g, *capacity, *flow);
  2.1013 -    bool _augment=false;
  2.1014 -
  2.1015 -    if (status!=AFTER_FAST_AUGMENTING) {
  2.1016 -      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); 
  2.1017 -      number_of_augmentations=1;
  2.1018 -    } else {
  2.1019 -      ++number_of_augmentations;
  2.1020 -    }
  2.1021 -    TrickyReachedMap<ReachedMap> 
  2.1022 -      tricky_reached_map(level, number_of_augmentations);
  2.1023 -    //ReachedMap level(res_graph);
  2.1024 -//    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
  2.1025 -    BfsIterator<ResGW, TrickyReachedMap<ReachedMap> > 
  2.1026 -      bfs(res_graph, tricky_reached_map);
  2.1027 -    bfs.pushAndSetReached(s);
  2.1028 -
  2.1029 -    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
  2.1030 -    pred.set(s, INVALID);
  2.1031 -
  2.1032 -    typename ResGW::template NodeMap<Num> free(res_graph);
  2.1033 -
  2.1034 -    //searching for augmenting path
  2.1035 -    while ( !bfs.finished() ) {
  2.1036 -      ResGWOutEdgeIt e=bfs;
  2.1037 -      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
  2.1038 -	Node v=res_graph.tail(e);
  2.1039 -	Node w=res_graph.head(e);
  2.1040 -	pred.set(w, e);
  2.1041 -	if (res_graph.valid(pred[v])) {
  2.1042 -	  free.set(w, std::min(free[v], res_graph.resCap(e)));
  2.1043 -	} else {
  2.1044 -	  free.set(w, res_graph.resCap(e));
  2.1045 -	}
  2.1046 -	if (res_graph.head(e)==t) { _augment=true; break; }
  2.1047 -      }
  2.1048 -
  2.1049 -      ++bfs;
  2.1050 -    } //end of searching augmenting path
  2.1051 -
  2.1052 -    if (_augment) {
  2.1053 -      Node n=t;
  2.1054 -      Num augment_value=free[t];
  2.1055 -      while (res_graph.valid(pred[n])) {
  2.1056 -	ResGWEdge e=pred[n];
  2.1057 -	res_graph.augment(e, augment_value);
  2.1058 -	n=res_graph.tail(e);
  2.1059 -      }
  2.1060 -    }
  2.1061 -
  2.1062 -    status=AFTER_FAST_AUGMENTING;
  2.1063 -    return _augment;
  2.1064 -  }
  2.1065 -
  2.1066 -
  2.1067 -  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  2.1068 -  template<typename MutableGraph>
  2.1069 -  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
  2.1070 -  {
  2.1071 -    typedef MutableGraph MG;
  2.1072 -    bool _augment=false;
  2.1073 -
  2.1074 -    ResGW res_graph(*g, *capacity, *flow);
  2.1075 -
  2.1076 -    //bfs for distances on the residual graph
  2.1077 -    //ReachedMap level(res_graph);
  2.1078 -    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
  2.1079 -    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
  2.1080 -    bfs.pushAndSetReached(s);
  2.1081 -    typename ResGW::template NodeMap<int>
  2.1082 -      dist(res_graph); //filled up with 0's
  2.1083 -
  2.1084 -    //F will contain the physical copy of the residual graph
  2.1085 -    //with the set of edges which are on shortest paths
  2.1086 -    MG F;
  2.1087 -    typename ResGW::template NodeMap<typename MG::Node>
  2.1088 -      res_graph_to_F(res_graph);
  2.1089 -    {
  2.1090 -      typename ResGW::NodeIt n;
  2.1091 -      for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
  2.1092 -	res_graph_to_F.set(n, F.addNode());
  2.1093 -      }
  2.1094 -    }
  2.1095 -
  2.1096 -    typename MG::Node sF=res_graph_to_F[s];
  2.1097 -    typename MG::Node tF=res_graph_to_F[t];
  2.1098 -    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
  2.1099 -    typename MG::template EdgeMap<Num> residual_capacity(F);
  2.1100 -
  2.1101 -    while ( !bfs.finished() ) {
  2.1102 -      ResGWOutEdgeIt e=bfs;
  2.1103 -      if (res_graph.valid(e)) {
  2.1104 -	if (bfs.isBNodeNewlyReached()) {
  2.1105 -	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
  2.1106 -	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
  2.1107 -					res_graph_to_F[res_graph.head(e)]);
  2.1108 -	  original_edge.update();
  2.1109 -	  original_edge.set(f, e);
  2.1110 -	  residual_capacity.update();
  2.1111 -	  residual_capacity.set(f, res_graph.resCap(e));
  2.1112 -	} else {
  2.1113 -	  if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
  2.1114 -	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
  2.1115 -					  res_graph_to_F[res_graph.head(e)]);
  2.1116 -	    original_edge.update();
  2.1117 -	    original_edge.set(f, e);
  2.1118 -	    residual_capacity.update();
  2.1119 -	    residual_capacity.set(f, res_graph.resCap(e));
  2.1120 -	  }
  2.1121 -	}
  2.1122 -      }
  2.1123 -      ++bfs;
  2.1124 -    } //computing distances from s in the residual graph
  2.1125 -
  2.1126 -    bool __augment=true;
  2.1127 -
  2.1128 -    while (__augment) {
  2.1129 -      __augment=false;
  2.1130 -      //computing blocking flow with dfs
  2.1131 -      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
  2.1132 -      typename MG::template NodeMap<typename MG::Edge> pred(F);
  2.1133 -      pred.set(sF, INVALID);
  2.1134 -      //invalid iterators for sources
  2.1135 -
  2.1136 -      typename MG::template NodeMap<Num> free(F);
  2.1137 -
  2.1138 -      dfs.pushAndSetReached(sF);
  2.1139 -      while (!dfs.finished()) {
  2.1140 -	++dfs;
  2.1141 -	if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
  2.1142 -	  if (dfs.isBNodeNewlyReached()) {
  2.1143 -	    typename MG::Node v=F.aNode(dfs);
  2.1144 -	    typename MG::Node w=F.bNode(dfs);
  2.1145 -	    pred.set(w, dfs);
  2.1146 -	    if (F.valid(pred[v])) {
  2.1147 -	      free.set(w, std::min(free[v], residual_capacity[dfs]));
  2.1148 -	    } else {
  2.1149 -	      free.set(w, residual_capacity[dfs]);
  2.1150 -	    }
  2.1151 -	    if (w==tF) {
  2.1152 -	      __augment=true;
  2.1153 -	      _augment=true;
  2.1154 -	      break;
  2.1155 -	    }
  2.1156 -
  2.1157 -	  } else {
  2.1158 -	    F.erase(/*typename MG::OutEdgeIt*/(dfs));
  2.1159 -	  }
  2.1160 -	}
  2.1161 -      }
  2.1162 -
  2.1163 -      if (__augment) {
  2.1164 -	typename MG::Node n=tF;
  2.1165 -	Num augment_value=free[tF];
  2.1166 -	while (F.valid(pred[n])) {
  2.1167 -	  typename MG::Edge e=pred[n];
  2.1168 -	  res_graph.augment(original_edge[e], augment_value);
  2.1169 -	  n=F.tail(e);
  2.1170 -	  if (residual_capacity[e]==augment_value)
  2.1171 -	    F.erase(e);
  2.1172 -	  else
  2.1173 -	    residual_capacity.set(e, residual_capacity[e]-augment_value);
  2.1174 -	}
  2.1175 -      }
  2.1176 -
  2.1177 -    }
  2.1178 -
  2.1179 -    status=AFTER_AUGMENTING;
  2.1180 -    return _augment;
  2.1181 -  }
  2.1182 -
  2.1183 -
  2.1184 -
  2.1185 -
  2.1186 -  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  2.1187 -  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
  2.1188 -  {
  2.1189 -    bool _augment=false;
  2.1190 -
  2.1191 -    ResGW res_graph(*g, *capacity, *flow);
  2.1192 -
  2.1193 -    //ReachedMap level(res_graph);
  2.1194 -    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
  2.1195 -    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
  2.1196 -
  2.1197 -    bfs.pushAndSetReached(s);
  2.1198 -    DistanceMap<ResGW> dist(res_graph);
  2.1199 -    while ( !bfs.finished() ) {
  2.1200 -      ResGWOutEdgeIt e=bfs;
  2.1201 -      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
  2.1202 -	dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
  2.1203 -      }
  2.1204 -      ++bfs;
  2.1205 -    } //computing distances from s in the residual graph
  2.1206 -
  2.1207 -      //Subgraph containing the edges on some shortest paths
  2.1208 -    ConstMap<typename ResGW::Node, bool> true_map(true);
  2.1209 -    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
  2.1210 -      DistanceMap<ResGW> > FilterResGW;
  2.1211 -    FilterResGW filter_res_graph(res_graph, true_map, dist);
  2.1212 -
  2.1213 -    //Subgraph, which is able to delete edges which are already
  2.1214 -    //met by the dfs
  2.1215 -    typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
  2.1216 -      first_out_edges(filter_res_graph);
  2.1217 -    typename FilterResGW::NodeIt v;
  2.1218 -    for(filter_res_graph.first(v); filter_res_graph.valid(v);
  2.1219 -	filter_res_graph.next(v))
  2.1220 -      {
  2.1221 - 	typename FilterResGW::OutEdgeIt e;
  2.1222 - 	filter_res_graph.first(e, v);
  2.1223 - 	first_out_edges.set(v, e);
  2.1224 -      }
  2.1225 -    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
  2.1226 -      template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
  2.1227 -    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
  2.1228 -
  2.1229 -    bool __augment=true;
  2.1230 -
  2.1231 -    while (__augment) {
  2.1232 -
  2.1233 -      __augment=false;
  2.1234 -      //computing blocking flow with dfs
  2.1235 -      DfsIterator< ErasingResGW,
  2.1236 -	typename ErasingResGW::template NodeMap<bool> >
  2.1237 -	dfs(erasing_res_graph);
  2.1238 -      typename ErasingResGW::
  2.1239 -	template NodeMap<typename ErasingResGW::OutEdgeIt>
  2.1240 -	pred(erasing_res_graph);
  2.1241 -      pred.set(s, INVALID);
  2.1242 -      //invalid iterators for sources
  2.1243 -
  2.1244 -      typename ErasingResGW::template NodeMap<Num>
  2.1245 -	free1(erasing_res_graph);
  2.1246 -
  2.1247 -      dfs.pushAndSetReached
  2.1248 -	///\bug hugo 0.2
  2.1249 -	(typename ErasingResGW::Node
  2.1250 -	 (typename FilterResGW::Node
  2.1251 -	  (typename ResGW::Node(s)
  2.1252 -	   )
  2.1253 -	  )
  2.1254 -	 );
  2.1255 -      while (!dfs.finished()) {
  2.1256 -	++dfs;
  2.1257 -	if (erasing_res_graph.valid(typename ErasingResGW::OutEdgeIt(dfs)))
  2.1258 - 	  {
  2.1259 -  	    if (dfs.isBNodeNewlyReached()) {
  2.1260 -
  2.1261 - 	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
  2.1262 - 	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
  2.1263 -
  2.1264 - 	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
  2.1265 - 	      if (erasing_res_graph.valid(pred[v])) {
  2.1266 - 		free1.set
  2.1267 -		  (w, std::min(free1[v], res_graph.resCap
  2.1268 -			       (typename ErasingResGW::OutEdgeIt(dfs))));
  2.1269 - 	      } else {
  2.1270 - 		free1.set
  2.1271 -		  (w, res_graph.resCap
  2.1272 -		   (typename ErasingResGW::OutEdgeIt(dfs)));
  2.1273 - 	      }
  2.1274 -
  2.1275 - 	      if (w==t) {
  2.1276 - 		__augment=true;
  2.1277 - 		_augment=true;
  2.1278 - 		break;
  2.1279 - 	      }
  2.1280 - 	    } else {
  2.1281 - 	      erasing_res_graph.erase(dfs);
  2.1282 -	    }
  2.1283 -	  }
  2.1284 -      }
  2.1285 -
  2.1286 -      if (__augment) {
  2.1287 -	typename ErasingResGW::Node
  2.1288 -	  n=typename FilterResGW::Node(typename ResGW::Node(t));
  2.1289 -	// 	  typename ResGW::NodeMap<Num> a(res_graph);
  2.1290 -	// 	  typename ResGW::Node b;
  2.1291 -	// 	  Num j=a[b];
  2.1292 -	// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
  2.1293 -	// 	  typename FilterResGW::Node b1;
  2.1294 -	// 	  Num j1=a1[b1];
  2.1295 -	// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
  2.1296 -	// 	  typename ErasingResGW::Node b2;
  2.1297 -	// 	  Num j2=a2[b2];
  2.1298 -	Num augment_value=free1[n];
  2.1299 -	while (erasing_res_graph.valid(pred[n])) {
  2.1300 -	  typename ErasingResGW::OutEdgeIt e=pred[n];
  2.1301 -	  res_graph.augment(e, augment_value);
  2.1302 -	  n=erasing_res_graph.tail(e);
  2.1303 -	  if (res_graph.resCap(e)==0)
  2.1304 -	    erasing_res_graph.erase(e);
  2.1305 -	}
  2.1306 -      }
  2.1307 -
  2.1308 -    } //while (__augment)
  2.1309 -
  2.1310 -    status=AFTER_AUGMENTING;
  2.1311 -    return _augment;
  2.1312 -  }
  2.1313 -
  2.1314 -
  2.1315 -} //namespace hugo
  2.1316 -
  2.1317 -#endif //HUGO_MAX_FLOW_H
  2.1318 -
  2.1319 -
  2.1320 -
  2.1321 -