lemon/hao_orlin.h
changeset 2530 f86f7e4eb2ba
parent 2391 14a343be7a5a
child 2553 bfced05fa852
     1.1 --- a/lemon/hao_orlin.h	Fri Nov 30 09:22:38 2007 +0000
     1.2 +++ b/lemon/hao_orlin.h	Tue Dec 04 10:55:27 2007 +0000
     1.3 @@ -19,13 +19,9 @@
     1.4  #ifndef LEMON_HAO_ORLIN_H
     1.5  #define LEMON_HAO_ORLIN_H
     1.6  
     1.7 -#include <cassert>
     1.8 - 
     1.9 -
    1.10 -
    1.11  #include <vector>
    1.12 -#include <queue>
    1.13  #include <list>
    1.14 +#include <ext/hash_set>
    1.15  #include <limits>
    1.16  
    1.17  #include <lemon/maps.h>
    1.18 @@ -37,7 +33,7 @@
    1.19  /// \ingroup min_cut
    1.20  /// \brief Implementation of the Hao-Orlin algorithm.
    1.21  ///
    1.22 -/// Implementation of the HaoOrlin algorithms class for testing network 
    1.23 +/// Implementation of the Hao-Orlin algorithm class for testing network 
    1.24  /// reliability.
    1.25  
    1.26  namespace lemon {
    1.27 @@ -46,24 +42,25 @@
    1.28    ///
    1.29    /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
    1.30    ///
    1.31 -  /// Hao-Orlin calculates a minimum cut in a directed graph 
    1.32 -  /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
    1.33 -  /// of two phases: in the first phase it determines a minimum cut
    1.34 -  /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
    1.35 -  /// with \f$ source \in X \f$ and minimal out-degree) and in the
    1.36 -  /// second phase it determines a minimum cut with \f$ source \f$ on the
    1.37 -  /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
    1.38 -  /// and minimal out-degree). Obviously, the smaller of these two
    1.39 -  /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
    1.40 -  /// modified push-relabel preflow algorithm and our implementation
    1.41 -  /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
    1.42 -  /// highest-label rule). The purpose of such an algorithm is testing
    1.43 -  /// network reliability. For an undirected graph with \f$ n \f$
    1.44 -  /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
    1.45 -  /// and Ibaraki which solves the undirected problem in 
    1.46 -  /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut 
    1.47 -  /// algorithm
    1.48 -  /// class.
    1.49 +  /// Hao-Orlin calculates a minimum cut in a directed graph
    1.50 +  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
    1.51 +  /// consists of two phases: in the first phase it determines a
    1.52 +  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
    1.53 +  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
    1.54 +  /// out-degree) and in the second phase it determines a minimum cut
    1.55 +  /// with \f$ source \f$ on the sink-side (i.e. a set 
    1.56 +  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
    1.57 +  /// out-degree). Obviously, the smaller of these two cuts will be a
    1.58 +  /// minimum cut of \f$ D \f$. The algorithm is a modified
    1.59 +  /// push-relabel preflow algorithm and our implementation calculates
    1.60 +  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
    1.61 +  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
    1.62 +  /// purpose of such algorithm is testing network reliability. For an
    1.63 +  /// undirected graph you can run just the first phase of the
    1.64 +  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
    1.65 +  /// which solves the undirected problem in 
    1.66 +  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
    1.67 +  /// NagamochiIbaraki algorithm class.
    1.68    ///
    1.69    /// \param _Graph is the graph type of the algorithm.
    1.70    /// \param _CapacityMap is an edge map of capacities which should
    1.71 @@ -80,7 +77,7 @@
    1.72              typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
    1.73  #endif
    1.74    class HaoOrlin {
    1.75 -  protected:
    1.76 +  private:
    1.77  
    1.78      typedef _Graph Graph;
    1.79      typedef _CapacityMap CapacityMap;
    1.80 @@ -88,44 +85,29 @@
    1.81  
    1.82      typedef typename CapacityMap::Value Value;
    1.83  
    1.84 +    GRAPH_TYPEDEFS(typename Graph);
    1.85      
    1.86 -    typedef typename Graph::Node Node;
    1.87 -    typedef typename Graph::NodeIt NodeIt;
    1.88 -    typedef typename Graph::EdgeIt EdgeIt;
    1.89 -    typedef typename Graph::OutEdgeIt OutEdgeIt;
    1.90 -    typedef typename Graph::InEdgeIt InEdgeIt;
    1.91 -
    1.92 -    const Graph* _graph;
    1.93 -
    1.94 +    const Graph& _graph;
    1.95      const CapacityMap* _capacity;
    1.96  
    1.97      typedef typename Graph::template EdgeMap<Value> FlowMap;
    1.98 +    FlowMap* _flow;
    1.99  
   1.100 -    FlowMap* _preflow;
   1.101 +    Node _source;
   1.102  
   1.103 -    Node _source, _target;
   1.104      int _node_num;
   1.105  
   1.106 -    typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 
   1.107 -                            FlowMap, Tolerance> OutResGraph;
   1.108 -    typedef typename OutResGraph::Edge OutResEdge;
   1.109 +    // Bucketing structure
   1.110 +    std::vector<Node> _first, _last;
   1.111 +    typename Graph::template NodeMap<Node>* _next;
   1.112 +    typename Graph::template NodeMap<Node>* _prev;    
   1.113 +    typename Graph::template NodeMap<bool>* _active;
   1.114 +    typename Graph::template NodeMap<int>* _bucket;
   1.115      
   1.116 -    OutResGraph* _out_res_graph;
   1.117 +    std::vector<bool> _dormant;
   1.118  
   1.119 -    typedef RevGraphAdaptor<const Graph> RevGraph;
   1.120 -    RevGraph* _rev_graph;
   1.121 -
   1.122 -    typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 
   1.123 -                            FlowMap, Tolerance> InResGraph;
   1.124 -    typedef typename InResGraph::Edge InResEdge;
   1.125 -    
   1.126 -    InResGraph* _in_res_graph;
   1.127 -
   1.128 -    typedef IterableBoolMap<Graph, Node> WakeMap;
   1.129 -    WakeMap* _wake;
   1.130 -
   1.131 -    typedef typename Graph::template NodeMap<int> DistMap;
   1.132 -    DistMap* _dist;  
   1.133 +    std::list<std::list<int> > _sets;
   1.134 +    std::list<int>::iterator _highest;
   1.135      
   1.136      typedef typename Graph::template NodeMap<Value> ExcessMap;
   1.137      ExcessMap* _excess;
   1.138 @@ -133,15 +115,6 @@
   1.139      typedef typename Graph::template NodeMap<bool> SourceSetMap;
   1.140      SourceSetMap* _source_set;
   1.141  
   1.142 -    std::vector<int> _level_size;
   1.143 -
   1.144 -    int _highest_active;
   1.145 -    std::vector<std::list<Node> > _active_nodes;
   1.146 -
   1.147 -    int _dormant_max;
   1.148 -    std::vector<std::list<Node> > _dormant;
   1.149 -
   1.150 -
   1.151      Value _min_cut;
   1.152  
   1.153      typedef typename Graph::template NodeMap<bool> MinCutMap;
   1.154 @@ -156,267 +129,695 @@
   1.155      /// Constructor of the algorithm class. 
   1.156      HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
   1.157               const Tolerance& tolerance = Tolerance()) :
   1.158 -      _graph(&graph), _capacity(&capacity), 
   1.159 -      _preflow(0), _source(), _target(), 
   1.160 -      _out_res_graph(0), _rev_graph(0), _in_res_graph(0),
   1.161 -      _wake(0),_dist(0), _excess(0), _source_set(0), 
   1.162 -      _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 
   1.163 -      _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
   1.164 +      _graph(graph), _capacity(&capacity), _flow(0), _source(),
   1.165 +      _node_num(), _first(), _last(), _next(0), _prev(0), 
   1.166 +      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
   1.167 +      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0), 
   1.168 +      _tolerance(tolerance) {}
   1.169  
   1.170      ~HaoOrlin() {
   1.171        if (_min_cut_map) {
   1.172          delete _min_cut_map;
   1.173        } 
   1.174 -      if (_in_res_graph) {
   1.175 -        delete _in_res_graph;
   1.176 -      }
   1.177 -      if (_rev_graph) {
   1.178 -        delete _rev_graph;
   1.179 -      }
   1.180 -      if (_out_res_graph) {
   1.181 -        delete _out_res_graph;
   1.182 -      }
   1.183        if (_source_set) {
   1.184          delete _source_set;
   1.185        }
   1.186        if (_excess) {
   1.187          delete _excess;
   1.188        }
   1.189 -      if (_dist) {
   1.190 -        delete _dist;
   1.191 +      if (_next) {
   1.192 +	delete _next;
   1.193        }
   1.194 -      if (_wake) {
   1.195 -        delete _wake;
   1.196 +      if (_prev) {
   1.197 +	delete _prev;
   1.198        }
   1.199 -      if (_preflow) {
   1.200 -        delete _preflow;
   1.201 +      if (_active) {
   1.202 +	delete _active;
   1.203 +      }
   1.204 +      if (_bucket) {
   1.205 +	delete _bucket;
   1.206 +      }
   1.207 +      if (_flow) {
   1.208 +        delete _flow;
   1.209        }
   1.210      }
   1.211      
   1.212    private:
   1.213 +
   1.214 +    void activate(const Node& i) {
   1.215 +      _active->set(i, true);
   1.216 +
   1.217 +      int bucket = (*_bucket)[i];
   1.218 +
   1.219 +      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;	    
   1.220 +      //unlace
   1.221 +      _next->set((*_prev)[i], (*_next)[i]);
   1.222 +      if ((*_next)[i] != INVALID) {
   1.223 +	_prev->set((*_next)[i], (*_prev)[i]);
   1.224 +      } else {
   1.225 +	_last[bucket] = (*_prev)[i];
   1.226 +      }
   1.227 +      //lace
   1.228 +      _next->set(i, _first[bucket]);
   1.229 +      _prev->set(_first[bucket], i);
   1.230 +      _prev->set(i, INVALID);
   1.231 +      _first[bucket] = i;
   1.232 +    }
   1.233 +
   1.234 +    void deactivate(const Node& i) {
   1.235 +      _active->set(i, false);
   1.236 +      int bucket = (*_bucket)[i];
   1.237 +
   1.238 +      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
   1.239 +      
   1.240 +      //unlace
   1.241 +      _prev->set((*_next)[i], (*_prev)[i]);
   1.242 +      if ((*_prev)[i] != INVALID) {
   1.243 +	_next->set((*_prev)[i], (*_next)[i]);
   1.244 +      } else {
   1.245 +	_first[bucket] = (*_next)[i];
   1.246 +      }
   1.247 +      //lace
   1.248 +      _prev->set(i, _last[bucket]);
   1.249 +      _next->set(_last[bucket], i);
   1.250 +      _next->set(i, INVALID);
   1.251 +      _last[bucket] = i;
   1.252 +    }
   1.253 +
   1.254 +    void addItem(const Node& i, int bucket) {
   1.255 +      (*_bucket)[i] = bucket;
   1.256 +      if (_last[bucket] != INVALID) {
   1.257 +	_prev->set(i, _last[bucket]);
   1.258 +	_next->set(_last[bucket], i);
   1.259 +	_next->set(i, INVALID);
   1.260 +	_last[bucket] = i;
   1.261 +      } else {
   1.262 +	_prev->set(i, INVALID);
   1.263 +	_first[bucket] = i;
   1.264 +	_next->set(i, INVALID);
   1.265 +	_last[bucket] = i;
   1.266 +      }
   1.267 +    }
   1.268      
   1.269 -    template <typename ResGraph>
   1.270 -    void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
   1.271 -      typedef typename ResGraph::Edge ResEdge;
   1.272 -      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   1.273 +    void findMinCutOut() {
   1.274  
   1.275 -      for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
   1.276 -        (*_preflow)[it] = 0;      
   1.277 -      }
   1.278 -      for (NodeIt it(*_graph); it != INVALID; ++it) {
   1.279 -        (*_wake)[it] = true;
   1.280 -        (*_dist)[it] = 1;
   1.281 -        (*_excess)[it] = 0;
   1.282 -        (*_source_set)[it] = false;
   1.283 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.284 +	_excess->set(n, 0);
   1.285        }
   1.286  
   1.287 -      _dormant[0].push_front(_source);
   1.288 -      (*_source_set)[_source] = true;
   1.289 -      _dormant_max = 0;
   1.290 -      (*_wake)[_source] = false;
   1.291 -
   1.292 -      _level_size[0] = 1;
   1.293 -      _level_size[1] = _node_num - 1;
   1.294 -
   1.295 -      _target = target;
   1.296 -      (*_dist)[target] = 0;
   1.297 -
   1.298 -      for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
   1.299 -        Value delta = res_graph.rescap(it);
   1.300 -        (*_excess)[_source] -= delta;
   1.301 -        res_graph.augment(it, delta);
   1.302 -        Node a = res_graph.target(it);
   1.303 -        if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
   1.304 -          _active_nodes[(*_dist)[a]].push_front(a);
   1.305 -          if (_highest_active < (*_dist)[a]) {
   1.306 -            _highest_active = (*_dist)[a];
   1.307 -          }
   1.308 -        }
   1.309 -        (*_excess)[a] += delta;
   1.310 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.311 +	_flow->set(e, 0);
   1.312        }
   1.313  
   1.314 +      int bucket_num = 1;
   1.315 +      
   1.316 +      {
   1.317 +	typename Graph::template NodeMap<bool> reached(_graph, false);
   1.318 +	
   1.319 +	reached.set(_source, true);
   1.320  
   1.321 -      do {
   1.322 -	Node n;
   1.323 -	while ((n = findActiveNode()) != INVALID) {
   1.324 -	  for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   1.325 -            Node a = res_graph.target(e);
   1.326 -            if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
   1.327 -	    Value delta = res_graph.rescap(e);
   1.328 -	    if (_tolerance.positive((*_excess)[n] - delta)) {
   1.329 -              (*_excess)[n] -= delta;
   1.330 -	    } else {
   1.331 -	      delta = (*_excess)[n];
   1.332 -              (*_excess)[n] = 0;
   1.333 +	bool first_set = true;
   1.334 +
   1.335 +	for (NodeIt t(_graph); t != INVALID; ++t) {
   1.336 +	  if (reached[t]) continue;
   1.337 +	  _sets.push_front(std::list<int>());
   1.338 +	  _sets.front().push_front(bucket_num);
   1.339 +	  _dormant[bucket_num] = !first_set;
   1.340 +
   1.341 +	  _bucket->set(t, bucket_num);
   1.342 +	  _first[bucket_num] = _last[bucket_num] = t;
   1.343 +	  _next->set(t, INVALID);
   1.344 +	  _prev->set(t, INVALID);
   1.345 +
   1.346 +	  ++bucket_num;
   1.347 +	  
   1.348 +	  std::vector<Node> queue;
   1.349 +	  queue.push_back(t);
   1.350 +	  reached.set(t, true);
   1.351 +	  
   1.352 +	  while (!queue.empty()) {
   1.353 +	    _sets.front().push_front(bucket_num);
   1.354 +	    _dormant[bucket_num] = !first_set;
   1.355 +	    _first[bucket_num] = _last[bucket_num] = INVALID;
   1.356 +	    
   1.357 +	    std::vector<Node> nqueue;
   1.358 +	    for (int i = 0; i < int(queue.size()); ++i) {
   1.359 +	      Node n = queue[i];
   1.360 +	      for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
   1.361 +		Node u = _graph.source(e);
   1.362 +		if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
   1.363 +		  reached.set(u, true);
   1.364 +		  addItem(u, bucket_num);
   1.365 +		  nqueue.push_back(u);
   1.366 +		}
   1.367 +	      }
   1.368  	    }
   1.369 -	    res_graph.augment(e, delta);
   1.370 -	    if ((*_excess)[a] == 0 && a != _target) {
   1.371 -	      _active_nodes[(*_dist)[a]].push_front(a);
   1.372 -	    }
   1.373 -	    (*_excess)[a] += delta;
   1.374 -            if ((*_excess)[n] == 0) break;
   1.375 +	    queue.swap(nqueue);
   1.376 +	    ++bucket_num;
   1.377  	  }
   1.378 -	  if ((*_excess)[n] != 0) {
   1.379 -	    relabel(n, res_graph);
   1.380 -          }
   1.381 +	  _sets.front().pop_front();
   1.382 +	  --bucket_num;
   1.383 +	  first_set = false;
   1.384  	}
   1.385  
   1.386 -	Value current_value = cutValue(out);
   1.387 - 	if (_min_cut > current_value){
   1.388 -          if (out) {
   1.389 -            for (NodeIt it(*_graph); it != INVALID; ++it) {
   1.390 -              _min_cut_map->set(it, !(*_wake)[it]);
   1.391 -            } 
   1.392 -          } else {
   1.393 -            for (NodeIt it(*_graph); it != INVALID; ++it) {
   1.394 -              _min_cut_map->set(it, (*_wake)[it]);
   1.395 -            } 
   1.396 -          }
   1.397 -
   1.398 -	  _min_cut = current_value;
   1.399 - 	}
   1.400 -
   1.401 -      } while (selectNewSink(res_graph));
   1.402 -    }
   1.403 -
   1.404 -    template <typename ResGraph>
   1.405 -    void relabel(const Node& n, ResGraph& res_graph) {
   1.406 -      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   1.407 -
   1.408 -      int k = (*_dist)[n];
   1.409 -      if (_level_size[k] == 1) {
   1.410 -	++_dormant_max;
   1.411 -	for (NodeIt it(*_graph); it != INVALID; ++it) {
   1.412 -	  if ((*_wake)[it] && (*_dist)[it] >= k) {
   1.413 -	    (*_wake)[it] = false;
   1.414 -	    _dormant[_dormant_max].push_front(it);
   1.415 -	    --_level_size[(*_dist)[it]];
   1.416 +	_bucket->set(_source, 0);
   1.417 +	_dormant[0] = true;
   1.418 +      }
   1.419 +      _source_set->set(_source, true);	  
   1.420 +	  
   1.421 +      Node target = _last[_sets.back().back()];
   1.422 +      {
   1.423 +	for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
   1.424 +	  if (_tolerance.positive((*_capacity)[e])) {
   1.425 +	    Node u = _graph.target(e);
   1.426 +	    _flow->set(e, (*_capacity)[e]);
   1.427 +	    _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
   1.428 +	    if (!(*_active)[u] && u != _source) {
   1.429 +	      activate(u);
   1.430 +	    }
   1.431  	  }
   1.432  	}
   1.433 -        --_highest_active;
   1.434 -      } else {	
   1.435 -        int new_dist = _node_num;
   1.436 -        for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   1.437 -          Node t = res_graph.target(e);
   1.438 -          if ((*_wake)[t] && new_dist > (*_dist)[t]) {
   1.439 -            new_dist = (*_dist)[t];
   1.440 -          }
   1.441 -        }
   1.442 -        if (new_dist == _node_num) {
   1.443 -	  ++_dormant_max;
   1.444 -	  (*_wake)[n] = false;
   1.445 -	  _dormant[_dormant_max].push_front(n);
   1.446 -	  --_level_size[(*_dist)[n]];
   1.447 -	} else {	    
   1.448 -	  --_level_size[(*_dist)[n]];
   1.449 -	  (*_dist)[n] = new_dist + 1;
   1.450 -	  _highest_active = (*_dist)[n];
   1.451 -	  _active_nodes[_highest_active].push_front(n);
   1.452 -	  ++_level_size[(*_dist)[n]];
   1.453 +	if ((*_active)[target]) {
   1.454 +	  deactivate(target);
   1.455 +	}
   1.456 +	
   1.457 +	_highest = _sets.back().begin();
   1.458 +	while (_highest != _sets.back().end() && 
   1.459 +	       !(*_active)[_first[*_highest]]) {
   1.460 +	  ++_highest;
   1.461 +	}
   1.462 +      }
   1.463 +
   1.464 +
   1.465 +      while (true) {
   1.466 +	while (_highest != _sets.back().end()) {
   1.467 +	  Node n = _first[*_highest];
   1.468 +	  Value excess = (*_excess)[n];
   1.469 +	  int next_bucket = _node_num;
   1.470 +
   1.471 +	  int under_bucket;
   1.472 +	  if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   1.473 +	    under_bucket = -1;
   1.474 +	  } else {
   1.475 +	    under_bucket = *(++std::list<int>::iterator(_highest));
   1.476 +	  }
   1.477 +
   1.478 +	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
   1.479 +	    Node v = _graph.target(e);
   1.480 +	    if (_dormant[(*_bucket)[v]]) continue;
   1.481 +	    Value rem = (*_capacity)[e] - (*_flow)[e];
   1.482 +	    if (!_tolerance.positive(rem)) continue;
   1.483 +	    if ((*_bucket)[v] == under_bucket) {
   1.484 +	      if (!(*_active)[v] && v != target) {
   1.485 +		activate(v);
   1.486 +	      }
   1.487 +	      if (!_tolerance.less(rem, excess)) {
   1.488 +		_flow->set(e, (*_flow)[e] + excess);
   1.489 +		_excess->set(v, (*_excess)[v] + excess);
   1.490 +		excess = 0;
   1.491 +		goto no_more_push;
   1.492 +	      } else {
   1.493 +		excess -= rem;
   1.494 +		_excess->set(v, (*_excess)[v] + rem);
   1.495 +		_flow->set(e, (*_capacity)[e]);
   1.496 +	      }
   1.497 +	    } else if (next_bucket > (*_bucket)[v]) {
   1.498 +	      next_bucket = (*_bucket)[v];
   1.499 +	    }
   1.500 +	  }
   1.501 +
   1.502 +	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
   1.503 +	    Node v = _graph.source(e);
   1.504 +	    if (_dormant[(*_bucket)[v]]) continue;
   1.505 +	    Value rem = (*_flow)[e];
   1.506 +	    if (!_tolerance.positive(rem)) continue;
   1.507 +	    if ((*_bucket)[v] == under_bucket) {
   1.508 +	      if (!(*_active)[v] && v != target) {
   1.509 +		activate(v);
   1.510 +	      }
   1.511 +	      if (!_tolerance.less(rem, excess)) {
   1.512 +		_flow->set(e, (*_flow)[e] - excess);
   1.513 +		_excess->set(v, (*_excess)[v] + excess);
   1.514 +		excess = 0;
   1.515 +		goto no_more_push;
   1.516 +	      } else {
   1.517 +		excess -= rem;
   1.518 +		_excess->set(v, (*_excess)[v] + rem);
   1.519 +		_flow->set(e, 0);
   1.520 +	      }
   1.521 +	    } else if (next_bucket > (*_bucket)[v]) {
   1.522 +	      next_bucket = (*_bucket)[v];
   1.523 +	    }
   1.524 +	  }
   1.525 +	  
   1.526 +	no_more_push:
   1.527 +	  
   1.528 +	  _excess->set(n, excess);
   1.529 +	  
   1.530 +	  if (excess != 0) {
   1.531 +	    if ((*_next)[n] == INVALID) {
   1.532 +	      typename std::list<std::list<int> >::iterator new_set = 
   1.533 +		_sets.insert(--_sets.end(), std::list<int>());
   1.534 +	      new_set->splice(new_set->end(), _sets.back(), 
   1.535 +			      _sets.back().begin(), ++_highest);
   1.536 +	      for (std::list<int>::iterator it = new_set->begin();
   1.537 +		   it != new_set->end(); ++it) {
   1.538 +		_dormant[*it] = true;
   1.539 +	      }
   1.540 +	      while (_highest != _sets.back().end() && 
   1.541 +		     !(*_active)[_first[*_highest]]) {
   1.542 +		++_highest;
   1.543 +	      }
   1.544 +	    } else if (next_bucket == _node_num) {
   1.545 +	      _first[(*_bucket)[n]] = (*_next)[n];
   1.546 +	      _prev->set((*_next)[n], INVALID);
   1.547 +	      
   1.548 +	      std::list<std::list<int> >::iterator new_set = 
   1.549 +		_sets.insert(--_sets.end(), std::list<int>());
   1.550 +	      
   1.551 +	      new_set->push_front(bucket_num);
   1.552 +	      _bucket->set(n, bucket_num);
   1.553 +	      _first[bucket_num] = _last[bucket_num] = n;
   1.554 +	      _next->set(n, INVALID);
   1.555 +	      _prev->set(n, INVALID);
   1.556 +	      _dormant[bucket_num] = true;
   1.557 +	      ++bucket_num;
   1.558 +
   1.559 +	      while (_highest != _sets.back().end() && 
   1.560 +		     !(*_active)[_first[*_highest]]) {
   1.561 +		++_highest;
   1.562 +	      }
   1.563 +	    } else {
   1.564 +	      _first[*_highest] = (*_next)[n];
   1.565 +	      _prev->set((*_next)[n], INVALID);
   1.566 +	      
   1.567 +	      while (next_bucket != *_highest) {
   1.568 +		--_highest;
   1.569 +	      }
   1.570 +
   1.571 +	      if (_highest == _sets.back().begin()) {
   1.572 +		_sets.back().push_front(bucket_num);
   1.573 +		_dormant[bucket_num] = false;
   1.574 +		_first[bucket_num] = _last[bucket_num] = INVALID;
   1.575 +		++bucket_num;
   1.576 +	      }
   1.577 +	      --_highest;
   1.578 +
   1.579 +	      _bucket->set(n, *_highest);
   1.580 +	      _next->set(n, _first[*_highest]);
   1.581 +	      if (_first[*_highest] != INVALID) {
   1.582 +		_prev->set(_first[*_highest], n);
   1.583 +	      } else {
   1.584 +		_last[*_highest] = n;
   1.585 +	      }
   1.586 +	      _first[*_highest] = n;	      
   1.587 +	    }
   1.588 +	  } else {
   1.589 +
   1.590 +	    deactivate(n);
   1.591 +	    if (!(*_active)[_first[*_highest]]) {
   1.592 +	      ++_highest;
   1.593 +	      if (_highest != _sets.back().end() && 
   1.594 +		  !(*_active)[_first[*_highest]]) {
   1.595 +		_highest = _sets.back().end();
   1.596 +	      }
   1.597 +	    }
   1.598 +	  }
   1.599 +	}
   1.600 +
   1.601 +	if ((*_excess)[target] < _min_cut) {
   1.602 +	  _min_cut = (*_excess)[target];
   1.603 +	  for (NodeIt i(_graph); i != INVALID; ++i) {
   1.604 +	    _min_cut_map->set(i, true);
   1.605 +	  }
   1.606 +	  for (std::list<int>::iterator it = _sets.back().begin();
   1.607 +	       it != _sets.back().end(); ++it) {
   1.608 +	    Node n = _first[*it];
   1.609 +	    while (n != INVALID) {
   1.610 +	      _min_cut_map->set(n, false);
   1.611 +	      n = (*_next)[n];
   1.612 +	    }
   1.613 +	  }
   1.614 +	}
   1.615 +
   1.616 +	{
   1.617 +	  Node new_target;
   1.618 +	  if ((*_prev)[target] != INVALID) {
   1.619 +	    _last[(*_bucket)[target]] = (*_prev)[target];
   1.620 +	    _next->set((*_prev)[target], INVALID);
   1.621 +	    new_target = (*_prev)[target];
   1.622 +	  } else {
   1.623 +	    _sets.back().pop_back();
   1.624 +	    if (_sets.back().empty()) {
   1.625 +	      _sets.pop_back();
   1.626 +	      if (_sets.empty())
   1.627 +		break;
   1.628 +	      for (std::list<int>::iterator it = _sets.back().begin();
   1.629 +		   it != _sets.back().end(); ++it) {
   1.630 +		_dormant[*it] = false;
   1.631 +	      }
   1.632 +	    }
   1.633 +	    new_target = _last[_sets.back().back()];
   1.634 +	  }
   1.635 +
   1.636 +	  _bucket->set(target, 0);
   1.637 +
   1.638 +	  _source_set->set(target, true);	  
   1.639 +	  for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
   1.640 +	    Value rem = (*_capacity)[e] - (*_flow)[e];
   1.641 +	    if (!_tolerance.positive(rem)) continue;
   1.642 +	    Node v = _graph.target(e);
   1.643 +	    if (!(*_active)[v] && !(*_source_set)[v]) {
   1.644 +	      activate(v);
   1.645 +	    }
   1.646 +	    _excess->set(v, (*_excess)[v] + rem);
   1.647 +	    _flow->set(e, (*_capacity)[e]);
   1.648 +	  }
   1.649 +	  
   1.650 +	  for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
   1.651 +	    Value rem = (*_flow)[e];
   1.652 +	    if (!_tolerance.positive(rem)) continue;
   1.653 +	    Node v = _graph.source(e);
   1.654 +	    if (!(*_active)[v] && !(*_source_set)[v]) {
   1.655 +	      activate(v);
   1.656 +	    }
   1.657 +	    _excess->set(v, (*_excess)[v] + rem);
   1.658 +	    _flow->set(e, 0);
   1.659 +	  }
   1.660 +	  
   1.661 +	  target = new_target;
   1.662 +	  if ((*_active)[target]) {
   1.663 +	    deactivate(target);
   1.664 +	  }
   1.665 +
   1.666 +	  _highest = _sets.back().begin();
   1.667 +	  while (_highest != _sets.back().end() && 
   1.668 +		 !(*_active)[_first[*_highest]]) {
   1.669 +	    ++_highest;
   1.670 +	  }
   1.671 +	}
   1.672 +      }
   1.673 +    }    
   1.674 +
   1.675 +    void findMinCutIn() {
   1.676 +
   1.677 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.678 +	_excess->set(n, 0);
   1.679 +      }
   1.680 +
   1.681 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.682 +	_flow->set(e, 0);
   1.683 +      }
   1.684 +
   1.685 +      int bucket_num = 1;
   1.686 +      
   1.687 +      {
   1.688 +	typename Graph::template NodeMap<bool> reached(_graph, false);
   1.689 +	
   1.690 +	reached.set(_source, true);
   1.691 +
   1.692 +	bool first_set = true;
   1.693 +
   1.694 +	for (NodeIt t(_graph); t != INVALID; ++t) {
   1.695 +	  if (reached[t]) continue;
   1.696 +	  _sets.push_front(std::list<int>());
   1.697 +	  _sets.front().push_front(bucket_num);
   1.698 +	  _dormant[bucket_num] = !first_set;
   1.699 +
   1.700 +	  _bucket->set(t, bucket_num);
   1.701 +	  _first[bucket_num] = _last[bucket_num] = t;
   1.702 +	  _next->set(t, INVALID);
   1.703 +	  _prev->set(t, INVALID);
   1.704 +
   1.705 +	  ++bucket_num;
   1.706 +	  
   1.707 +	  std::vector<Node> queue;
   1.708 +	  queue.push_back(t);
   1.709 +	  reached.set(t, true);
   1.710 +	  
   1.711 +	  while (!queue.empty()) {
   1.712 +	    _sets.front().push_front(bucket_num);
   1.713 +	    _dormant[bucket_num] = !first_set;
   1.714 +	    _first[bucket_num] = _last[bucket_num] = INVALID;
   1.715 +	    
   1.716 +	    std::vector<Node> nqueue;
   1.717 +	    for (int i = 0; i < int(queue.size()); ++i) {
   1.718 +	      Node n = queue[i];
   1.719 +	      for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
   1.720 +		Node u = _graph.target(e);
   1.721 +		if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
   1.722 +		  reached.set(u, true);
   1.723 +		  addItem(u, bucket_num);
   1.724 +		  nqueue.push_back(u);
   1.725 +		}
   1.726 +	      }
   1.727 +	    }
   1.728 +	    queue.swap(nqueue);
   1.729 +	    ++bucket_num;
   1.730 +	  }
   1.731 +	  _sets.front().pop_front();
   1.732 +	  --bucket_num;
   1.733 +	  first_set = false;
   1.734 +	}
   1.735 +
   1.736 +	_bucket->set(_source, 0);
   1.737 +	_dormant[0] = true;
   1.738 +      }
   1.739 +      _source_set->set(_source, true);	  
   1.740 +	  
   1.741 +      Node target = _last[_sets.back().back()];
   1.742 +      {
   1.743 +	for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
   1.744 +	  if (_tolerance.positive((*_capacity)[e])) {
   1.745 +	    Node u = _graph.source(e);
   1.746 +	    _flow->set(e, (*_capacity)[e]);
   1.747 +	    _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
   1.748 +	    if (!(*_active)[u] && u != _source) {
   1.749 +	      activate(u);
   1.750 +	    }
   1.751 +	  }
   1.752 +	}
   1.753 +	if ((*_active)[target]) {
   1.754 +	  deactivate(target);
   1.755 +	}
   1.756 +	
   1.757 +	_highest = _sets.back().begin();
   1.758 +	while (_highest != _sets.back().end() && 
   1.759 +	       !(*_active)[_first[*_highest]]) {
   1.760 +	  ++_highest;
   1.761 +	}
   1.762 +      }
   1.763 +
   1.764 +
   1.765 +      while (true) {
   1.766 +	while (_highest != _sets.back().end()) {
   1.767 +	  Node n = _first[*_highest];
   1.768 +	  Value excess = (*_excess)[n];
   1.769 +	  int next_bucket = _node_num;
   1.770 +
   1.771 +	  int under_bucket;
   1.772 +	  if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   1.773 +	    under_bucket = -1;
   1.774 +	  } else {
   1.775 +	    under_bucket = *(++std::list<int>::iterator(_highest));
   1.776 +	  }
   1.777 +
   1.778 +	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
   1.779 +	    Node v = _graph.source(e);
   1.780 +	    if (_dormant[(*_bucket)[v]]) continue;
   1.781 +	    Value rem = (*_capacity)[e] - (*_flow)[e];
   1.782 +	    if (!_tolerance.positive(rem)) continue;
   1.783 +	    if ((*_bucket)[v] == under_bucket) {
   1.784 +	      if (!(*_active)[v] && v != target) {
   1.785 +		activate(v);
   1.786 +	      }
   1.787 +	      if (!_tolerance.less(rem, excess)) {
   1.788 +		_flow->set(e, (*_flow)[e] + excess);
   1.789 +		_excess->set(v, (*_excess)[v] + excess);
   1.790 +		excess = 0;
   1.791 +		goto no_more_push;
   1.792 +	      } else {
   1.793 +		excess -= rem;
   1.794 +		_excess->set(v, (*_excess)[v] + rem);
   1.795 +		_flow->set(e, (*_capacity)[e]);
   1.796 +	      }
   1.797 +	    } else if (next_bucket > (*_bucket)[v]) {
   1.798 +	      next_bucket = (*_bucket)[v];
   1.799 +	    }
   1.800 +	  }
   1.801 +
   1.802 +	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
   1.803 +	    Node v = _graph.target(e);
   1.804 +	    if (_dormant[(*_bucket)[v]]) continue;
   1.805 +	    Value rem = (*_flow)[e];
   1.806 +	    if (!_tolerance.positive(rem)) continue;
   1.807 +	    if ((*_bucket)[v] == under_bucket) {
   1.808 +	      if (!(*_active)[v] && v != target) {
   1.809 +		activate(v);
   1.810 +	      }
   1.811 +	      if (!_tolerance.less(rem, excess)) {
   1.812 +		_flow->set(e, (*_flow)[e] - excess);
   1.813 +		_excess->set(v, (*_excess)[v] + excess);
   1.814 +		excess = 0;
   1.815 +		goto no_more_push;
   1.816 +	      } else {
   1.817 +		excess -= rem;
   1.818 +		_excess->set(v, (*_excess)[v] + rem);
   1.819 +		_flow->set(e, 0);
   1.820 +	      }
   1.821 +	    } else if (next_bucket > (*_bucket)[v]) {
   1.822 +	      next_bucket = (*_bucket)[v];
   1.823 +	    }
   1.824 +	  }
   1.825 +	  
   1.826 +	no_more_push:
   1.827 +	  
   1.828 +	  _excess->set(n, excess);
   1.829 +	  
   1.830 +	  if (excess != 0) {
   1.831 +	    if ((*_next)[n] == INVALID) {
   1.832 +	      typename std::list<std::list<int> >::iterator new_set = 
   1.833 +		_sets.insert(--_sets.end(), std::list<int>());
   1.834 +	      new_set->splice(new_set->end(), _sets.back(), 
   1.835 +			      _sets.back().begin(), ++_highest);
   1.836 +	      for (std::list<int>::iterator it = new_set->begin();
   1.837 +		   it != new_set->end(); ++it) {
   1.838 +		_dormant[*it] = true;
   1.839 +	      }
   1.840 +	      while (_highest != _sets.back().end() && 
   1.841 +		     !(*_active)[_first[*_highest]]) {
   1.842 +		++_highest;
   1.843 +	      }
   1.844 +	    } else if (next_bucket == _node_num) {
   1.845 +	      _first[(*_bucket)[n]] = (*_next)[n];
   1.846 +	      _prev->set((*_next)[n], INVALID);
   1.847 +	      
   1.848 +	      std::list<std::list<int> >::iterator new_set = 
   1.849 +		_sets.insert(--_sets.end(), std::list<int>());
   1.850 +	      
   1.851 +	      new_set->push_front(bucket_num);
   1.852 +	      _bucket->set(n, bucket_num);
   1.853 +	      _first[bucket_num] = _last[bucket_num] = n;
   1.854 +	      _next->set(n, INVALID);
   1.855 +	      _prev->set(n, INVALID);
   1.856 +	      _dormant[bucket_num] = true;
   1.857 +	      ++bucket_num;
   1.858 +
   1.859 +	      while (_highest != _sets.back().end() && 
   1.860 +		     !(*_active)[_first[*_highest]]) {
   1.861 +		++_highest;
   1.862 +	      }
   1.863 +	    } else {
   1.864 +	      _first[*_highest] = (*_next)[n];
   1.865 +	      _prev->set((*_next)[n], INVALID);
   1.866 +
   1.867 +	      while (next_bucket != *_highest) {
   1.868 +		--_highest;
   1.869 +	      }
   1.870 +	      if (_highest == _sets.back().begin()) {
   1.871 +		_sets.back().push_front(bucket_num);
   1.872 +		_dormant[bucket_num] = false;
   1.873 +		_first[bucket_num] = _last[bucket_num] = INVALID;
   1.874 +		++bucket_num;
   1.875 +	      }
   1.876 +	      --_highest;
   1.877 +
   1.878 +	      _bucket->set(n, *_highest);
   1.879 +	      _next->set(n, _first[*_highest]);
   1.880 +	      if (_first[*_highest] != INVALID) {
   1.881 +		_prev->set(_first[*_highest], n);
   1.882 +	      } else {
   1.883 +		_last[*_highest] = n;
   1.884 +	      }
   1.885 +	      _first[*_highest] = n;	      
   1.886 +	    }
   1.887 +	  } else {
   1.888 +
   1.889 +	    deactivate(n);
   1.890 +	    if (!(*_active)[_first[*_highest]]) {
   1.891 +	      ++_highest;
   1.892 +	      if (_highest != _sets.back().end() && 
   1.893 +		  !(*_active)[_first[*_highest]]) {
   1.894 +		_highest = _sets.back().end();
   1.895 +	      }
   1.896 +	    }
   1.897 +	  }
   1.898 +	}
   1.899 +
   1.900 +	if ((*_excess)[target] < _min_cut) {
   1.901 +	  _min_cut = (*_excess)[target];
   1.902 +	  for (NodeIt i(_graph); i != INVALID; ++i) {
   1.903 +	    _min_cut_map->set(i, false);
   1.904 +	  }
   1.905 +	  for (std::list<int>::iterator it = _sets.back().begin();
   1.906 +	       it != _sets.back().end(); ++it) {
   1.907 +	    Node n = _first[*it];
   1.908 +	    while (n != INVALID) {
   1.909 +	      _min_cut_map->set(n, true);
   1.910 +	      n = (*_next)[n];
   1.911 +	    }
   1.912 +	  }
   1.913 +	}
   1.914 +
   1.915 +	{
   1.916 +	  Node new_target;
   1.917 +	  if ((*_prev)[target] != INVALID) {
   1.918 +	    _last[(*_bucket)[target]] = (*_prev)[target];
   1.919 +	    _next->set((*_prev)[target], INVALID);
   1.920 +	    new_target = (*_prev)[target];
   1.921 +	  } else {
   1.922 +	    _sets.back().pop_back();
   1.923 +	    if (_sets.back().empty()) {
   1.924 +	      _sets.pop_back();
   1.925 +	      if (_sets.empty())
   1.926 +		break;
   1.927 +	      for (std::list<int>::iterator it = _sets.back().begin();
   1.928 +		   it != _sets.back().end(); ++it) {
   1.929 +		_dormant[*it] = false;
   1.930 +	      }
   1.931 +	    }
   1.932 +	    new_target = _last[_sets.back().back()];
   1.933 +	  }
   1.934 +
   1.935 +	  _bucket->set(target, 0);
   1.936 +
   1.937 +	  _source_set->set(target, true);	  
   1.938 +	  for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
   1.939 +	    Value rem = (*_capacity)[e] - (*_flow)[e];
   1.940 +	    if (!_tolerance.positive(rem)) continue;
   1.941 +	    Node v = _graph.source(e);
   1.942 +	    if (!(*_active)[v] && !(*_source_set)[v]) {
   1.943 +	      activate(v);
   1.944 +	    }
   1.945 +	    _excess->set(v, (*_excess)[v] + rem);
   1.946 +	    _flow->set(e, (*_capacity)[e]);
   1.947 +	  }
   1.948 +	  
   1.949 +	  for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
   1.950 +	    Value rem = (*_flow)[e];
   1.951 +	    if (!_tolerance.positive(rem)) continue;
   1.952 +	    Node v = _graph.target(e);
   1.953 +	    if (!(*_active)[v] && !(*_source_set)[v]) {
   1.954 +	      activate(v);
   1.955 +	    }
   1.956 +	    _excess->set(v, (*_excess)[v] + rem);
   1.957 +	    _flow->set(e, 0);
   1.958 +	  }
   1.959 +	  
   1.960 +	  target = new_target;
   1.961 +	  if ((*_active)[target]) {
   1.962 +	    deactivate(target);
   1.963 +	  }
   1.964 +
   1.965 +	  _highest = _sets.back().begin();
   1.966 +	  while (_highest != _sets.back().end() && 
   1.967 +		 !(*_active)[_first[*_highest]]) {
   1.968 +	    ++_highest;
   1.969 +	  }
   1.970  	}
   1.971        }
   1.972      }
   1.973  
   1.974 -    template <typename ResGraph>
   1.975 -    bool selectNewSink(ResGraph& res_graph) {
   1.976 -      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   1.977 -
   1.978 -      Node old_target = _target;
   1.979 -      (*_wake)[_target] = false;
   1.980 -      --_level_size[(*_dist)[_target]];
   1.981 -      _dormant[0].push_front(_target);
   1.982 -      (*_source_set)[_target] = true;
   1.983 -      if (int(_dormant[0].size()) == _node_num){
   1.984 -        _dormant[0].clear();
   1.985 -	return false;
   1.986 -      }
   1.987 -
   1.988 -      bool wake_was_empty = false;
   1.989 -
   1.990 -      if(_wake->trueNum() == 0) {
   1.991 -	while (!_dormant[_dormant_max].empty()){
   1.992 -	  (*_wake)[_dormant[_dormant_max].front()] = true;
   1.993 -	  ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
   1.994 -	  _dormant[_dormant_max].pop_front();
   1.995 -	}
   1.996 -	--_dormant_max;
   1.997 -	wake_was_empty = true;
   1.998 -      }
   1.999 -
  1.1000 -      int min_dist = std::numeric_limits<int>::max();
  1.1001 -      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
  1.1002 -	if (min_dist > (*_dist)[it]){
  1.1003 -	  _target = it;
  1.1004 -	  min_dist = (*_dist)[it];
  1.1005 -	}
  1.1006 -      }
  1.1007 -
  1.1008 -      if (wake_was_empty){
  1.1009 -	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
  1.1010 -          if ((*_excess)[it] != 0 && it != _target) {
  1.1011 -            _active_nodes[(*_dist)[it]].push_front(it);
  1.1012 -            if (_highest_active < (*_dist)[it]) {
  1.1013 -              _highest_active = (*_dist)[it];		    
  1.1014 -            }
  1.1015 -	  }
  1.1016 -	}
  1.1017 -      }
  1.1018 -
  1.1019 -      Node n = old_target;
  1.1020 -      for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
  1.1021 -        Node a = res_graph.target(e);
  1.1022 -	if (!(*_wake)[a]) continue;
  1.1023 -        Value delta = res_graph.rescap(e);
  1.1024 -        res_graph.augment(e, delta);
  1.1025 -        (*_excess)[n] -= delta;
  1.1026 -        if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
  1.1027 -          _active_nodes[(*_dist)[a]].push_front(a);
  1.1028 -          if (_highest_active < (*_dist)[a]) {
  1.1029 -            _highest_active = (*_dist)[a];
  1.1030 -          }
  1.1031 -        }
  1.1032 -        (*_excess)[a] += delta;
  1.1033 -      }
  1.1034 -      
  1.1035 -      return true;
  1.1036 -    }
  1.1037 -
  1.1038 -    Node findActiveNode() {
  1.1039 -      while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
  1.1040 -	--_highest_active;
  1.1041 -      }
  1.1042 -      if( _highest_active > 0) {
  1.1043 -       	Node n = _active_nodes[_highest_active].front();
  1.1044 -	_active_nodes[_highest_active].pop_front();
  1.1045 -	return n;
  1.1046 -      } else {
  1.1047 -	return INVALID;
  1.1048 -      }
  1.1049 -    }
  1.1050 -
  1.1051 -    Value cutValue(bool out) {
  1.1052 -      Value value = 0;
  1.1053 -      if (out) {
  1.1054 -        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
  1.1055 -          for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
  1.1056 -            if (!(*_wake)[_graph->source(e)]){
  1.1057 -              value += (*_capacity)[e];
  1.1058 -            }	
  1.1059 -          }
  1.1060 -        }
  1.1061 -      } else {
  1.1062 -        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
  1.1063 -          for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
  1.1064 -            if (!(*_wake)[_graph->target(e)]){
  1.1065 -              value += (*_capacity)[e];
  1.1066 -            }	
  1.1067 -          }
  1.1068 -        }
  1.1069 -      }
  1.1070 -      return value;
  1.1071 -    }
  1.1072 -
  1.1073 -
  1.1074    public:
  1.1075  
  1.1076      /// \name Execution control
  1.1077 @@ -435,7 +836,7 @@
  1.1078      /// the maps, residual graph adaptors and some bucket structures
  1.1079      /// for the algorithm. 
  1.1080      void init() {
  1.1081 -      init(NodeIt(*_graph));
  1.1082 +      init(NodeIt(_graph));
  1.1083      }
  1.1084  
  1.1085      /// \brief Initializes the internal data structures.
  1.1086 @@ -446,38 +847,37 @@
  1.1087      /// algorithm's source.
  1.1088      void init(const Node& source) {
  1.1089        _source = source;
  1.1090 -      _node_num = countNodes(*_graph);
  1.1091 +      
  1.1092 +      _node_num = countNodes(_graph);
  1.1093 +      
  1.1094 +      _first.resize(_node_num);
  1.1095 +      _last.resize(_node_num);
  1.1096  
  1.1097        _dormant.resize(_node_num);
  1.1098 -      _level_size.resize(_node_num, 0);
  1.1099 -      _active_nodes.resize(_node_num);
  1.1100  
  1.1101 -      if (!_preflow) {
  1.1102 -        _preflow = new FlowMap(*_graph);
  1.1103 +      if (!_flow) {
  1.1104 +	_flow = new FlowMap(_graph);
  1.1105        }
  1.1106 -      if (!_wake) {
  1.1107 -        _wake = new WakeMap(*_graph);
  1.1108 +      if (!_next) {
  1.1109 +	_next = new typename Graph::template NodeMap<Node>(_graph);
  1.1110        }
  1.1111 -      if (!_dist) {
  1.1112 -        _dist = new DistMap(*_graph);
  1.1113 +      if (!_prev) {
  1.1114 +	_prev = new typename Graph::template NodeMap<Node>(_graph);
  1.1115 +      }
  1.1116 +      if (!_active) {
  1.1117 +	_active = new typename Graph::template NodeMap<bool>(_graph);
  1.1118 +      }
  1.1119 +      if (!_bucket) {
  1.1120 +	_bucket = new typename Graph::template NodeMap<int>(_graph);
  1.1121        }
  1.1122        if (!_excess) {
  1.1123 -        _excess = new ExcessMap(*_graph);
  1.1124 +	_excess = new ExcessMap(_graph);
  1.1125        }
  1.1126        if (!_source_set) {
  1.1127 -        _source_set = new SourceSetMap(*_graph);
  1.1128 -      }
  1.1129 -      if (!_out_res_graph) {
  1.1130 -        _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
  1.1131 -      }
  1.1132 -      if (!_rev_graph) {
  1.1133 -        _rev_graph = new RevGraph(*_graph);
  1.1134 -      }
  1.1135 -      if (!_in_res_graph) {
  1.1136 -        _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
  1.1137 +	_source_set = new SourceSetMap(_graph);
  1.1138        }
  1.1139        if (!_min_cut_map) {
  1.1140 -        _min_cut_map = new MinCutMap(*_graph);
  1.1141 +	_min_cut_map = new MinCutMap(_graph);
  1.1142        }
  1.1143  
  1.1144        _min_cut = std::numeric_limits<Value>::max();
  1.1145 @@ -487,56 +887,23 @@
  1.1146      /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1147      /// source-side.
  1.1148      ///
  1.1149 -    /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1150 -    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
  1.1151 -    ///  and minimal out-degree).
  1.1152 +    /// Calculates a minimum cut with \f$ source \f$ on the
  1.1153 +    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
  1.1154 +    /// \in X \f$ and minimal out-degree).
  1.1155      void calculateOut() {
  1.1156 -      for (NodeIt it(*_graph); it != INVALID; ++it) {
  1.1157 -        if (it != _source) {
  1.1158 -          calculateOut(it);
  1.1159 -          return;
  1.1160 -        }
  1.1161 -      }
  1.1162 +      findMinCutOut();
  1.1163      }
  1.1164  
  1.1165      /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1166 -    /// source-side.
  1.1167 +    /// target-side.
  1.1168      ///
  1.1169 -    /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1170 -    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
  1.1171 -    ///  and minimal out-degree). The \c target is the initial target
  1.1172 -    /// for the push-relabel algorithm.
  1.1173 -    void calculateOut(const Node& target) {
  1.1174 -      findMinCut(target, true, *_out_res_graph);
  1.1175 +    /// Calculates a minimum cut with \f$ source \f$ on the
  1.1176 +    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
  1.1177 +    /// \in X \f$ and minimal out-degree).
  1.1178 +    void calculateIn() {
  1.1179 +      findMinCutIn();
  1.1180      }
  1.1181  
  1.1182 -    /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1183 -    /// sink-side.
  1.1184 -    ///
  1.1185 -    /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1186 -    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with 
  1.1187 -    /// \f$ source \notin X \f$
  1.1188 -    /// and minimal out-degree).
  1.1189 -    void calculateIn() {
  1.1190 -      for (NodeIt it(*_graph); it != INVALID; ++it) {
  1.1191 -        if (it != _source) {
  1.1192 -          calculateIn(it);
  1.1193 -          return;
  1.1194 -        }
  1.1195 -      }
  1.1196 -    }
  1.1197 -
  1.1198 -    /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1199 -    /// sink-side.
  1.1200 -    ///
  1.1201 -    /// \brief Calculates a minimum cut with \f$ source \f$ on the
  1.1202 -    /// sink-side (i.e. a set \f$ X\subsetneq V 
  1.1203 -    /// \f$ with \f$ source \notin X \f$ and minimal out-degree).  
  1.1204 -    /// The \c target is the initial
  1.1205 -    /// target for the push-relabel algorithm.
  1.1206 -    void calculateIn(const Node& target) {
  1.1207 -      findMinCut(target, false, *_in_res_graph);
  1.1208 -    }
  1.1209  
  1.1210      /// \brief Runs the algorithm.
  1.1211      ///
  1.1212 @@ -545,13 +912,8 @@
  1.1213      /// and \ref calculateIn().
  1.1214      void run() {
  1.1215        init();
  1.1216 -      for (NodeIt it(*_graph); it != INVALID; ++it) {
  1.1217 -        if (it != _source) {
  1.1218 -          calculateOut(it);
  1.1219 -          calculateIn(it);
  1.1220 -          return;
  1.1221 -        }
  1.1222 -      }
  1.1223 +      calculateOut();
  1.1224 +      calculateIn();
  1.1225      }
  1.1226  
  1.1227      /// \brief Runs the algorithm.
  1.1228 @@ -561,23 +923,8 @@
  1.1229      /// calculateOut() and \ref calculateIn().
  1.1230      void run(const Node& s) {
  1.1231        init(s);
  1.1232 -      for (NodeIt it(*_graph); it != INVALID; ++it) {
  1.1233 -        if (it != _source) {
  1.1234 -          calculateOut(it);
  1.1235 -          calculateIn(it);
  1.1236 -          return;
  1.1237 -        }
  1.1238 -      }
  1.1239 -    }
  1.1240 -
  1.1241 -    /// \brief Runs the algorithm.
  1.1242 -    ///
  1.1243 -    /// Runs the algorithm. It just calls the \ref init() and then
  1.1244 -    /// \ref calculateOut() and \ref calculateIn().
  1.1245 -    void run(const Node& s, const Node& t) {
  1.1246 -      init(s); 
  1.1247 -      calculateOut(t);
  1.1248 -      calculateIn(t);
  1.1249 +      calculateOut();
  1.1250 +      calculateIn();
  1.1251      }
  1.1252  
  1.1253      /// @}
  1.1254 @@ -608,7 +955,7 @@
  1.1255      /// bool-valued node-map.
  1.1256      template <typename NodeMap>
  1.1257      Value minCut(NodeMap& nodeMap) const {
  1.1258 -      for (NodeIt it(*_graph); it != INVALID; ++it) {
  1.1259 +      for (NodeIt it(_graph); it != INVALID; ++it) {
  1.1260  	nodeMap.set(it, (*_min_cut_map)[it]);
  1.1261        }
  1.1262        return minCut();