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