lemon/hao_orlin.h
changeset 2530 f86f7e4eb2ba
parent 2391 14a343be7a5a
child 2553 bfced05fa852
equal deleted inserted replaced
8:b3d5ecb478ed 9:fe253f918d60
    17  */
    17  */
    18 
    18 
    19 #ifndef LEMON_HAO_ORLIN_H
    19 #ifndef LEMON_HAO_ORLIN_H
    20 #define LEMON_HAO_ORLIN_H
    20 #define LEMON_HAO_ORLIN_H
    21 
    21 
    22 #include <cassert>
       
    23  
       
    24 
       
    25 
       
    26 #include <vector>
    22 #include <vector>
    27 #include <queue>
       
    28 #include <list>
    23 #include <list>
       
    24 #include <ext/hash_set>
    29 #include <limits>
    25 #include <limits>
    30 
    26 
    31 #include <lemon/maps.h>
    27 #include <lemon/maps.h>
    32 #include <lemon/graph_utils.h>
    28 #include <lemon/graph_utils.h>
    33 #include <lemon/graph_adaptor.h>
    29 #include <lemon/graph_adaptor.h>
    35 
    31 
    36 /// \file
    32 /// \file
    37 /// \ingroup min_cut
    33 /// \ingroup min_cut
    38 /// \brief Implementation of the Hao-Orlin algorithm.
    34 /// \brief Implementation of the Hao-Orlin algorithm.
    39 ///
    35 ///
    40 /// Implementation of the HaoOrlin algorithms class for testing network 
    36 /// Implementation of the Hao-Orlin algorithm class for testing network 
    41 /// reliability.
    37 /// reliability.
    42 
    38 
    43 namespace lemon {
    39 namespace lemon {
    44 
    40 
    45   /// \ingroup min_cut
    41   /// \ingroup min_cut
    46   ///
    42   ///
    47   /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
    43   /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
    48   ///
    44   ///
    49   /// Hao-Orlin calculates a minimum cut in a directed graph 
    45   /// Hao-Orlin calculates a minimum cut in a directed graph
    50   /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
    46   /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
    51   /// of two phases: in the first phase it determines a minimum cut
    47   /// consists of two phases: in the first phase it determines a
    52   /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
    48   /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
    53   /// with \f$ source \in X \f$ and minimal out-degree) and in the
    49   /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
    54   /// second phase it determines a minimum cut with \f$ source \f$ on the
    50   /// out-degree) and in the second phase it determines a minimum cut
    55   /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
    51   /// with \f$ source \f$ on the sink-side (i.e. a set 
    56   /// and minimal out-degree). Obviously, the smaller of these two
    52   /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
    57   /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
    53   /// out-degree). Obviously, the smaller of these two cuts will be a
    58   /// modified push-relabel preflow algorithm and our implementation
    54   /// minimum cut of \f$ D \f$. The algorithm is a modified
    59   /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
    55   /// push-relabel preflow algorithm and our implementation calculates
    60   /// highest-label rule). The purpose of such an algorithm is testing
    56   /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
    61   /// network reliability. For an undirected graph with \f$ n \f$
    57   /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
    62   /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
    58   /// purpose of such algorithm is testing network reliability. For an
    63   /// and Ibaraki which solves the undirected problem in 
    59   /// undirected graph you can run just the first phase of the
    64   /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut 
    60   /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
    65   /// algorithm
    61   /// which solves the undirected problem in 
    66   /// class.
    62   /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
       
    63   /// NagamochiIbaraki algorithm class.
    67   ///
    64   ///
    68   /// \param _Graph is the graph type of the algorithm.
    65   /// \param _Graph is the graph type of the algorithm.
    69   /// \param _CapacityMap is an edge map of capacities which should
    66   /// \param _CapacityMap is an edge map of capacities which should
    70   /// be any numreric type. The default type is _Graph::EdgeMap<int>.
    67   /// be any numreric type. The default type is _Graph::EdgeMap<int>.
    71   /// \param _Tolerance is the handler of the inexact computation. The
    68   /// \param _Tolerance is the handler of the inexact computation. The
    78   template <typename _Graph,
    75   template <typename _Graph,
    79 	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
    76 	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
    80             typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
    77             typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
    81 #endif
    78 #endif
    82   class HaoOrlin {
    79   class HaoOrlin {
    83   protected:
    80   private:
    84 
    81 
    85     typedef _Graph Graph;
    82     typedef _Graph Graph;
    86     typedef _CapacityMap CapacityMap;
    83     typedef _CapacityMap CapacityMap;
    87     typedef _Tolerance Tolerance;
    84     typedef _Tolerance Tolerance;
    88 
    85 
    89     typedef typename CapacityMap::Value Value;
    86     typedef typename CapacityMap::Value Value;
    90 
    87 
       
    88     GRAPH_TYPEDEFS(typename Graph);
    91     
    89     
    92     typedef typename Graph::Node Node;
    90     const Graph& _graph;
    93     typedef typename Graph::NodeIt NodeIt;
       
    94     typedef typename Graph::EdgeIt EdgeIt;
       
    95     typedef typename Graph::OutEdgeIt OutEdgeIt;
       
    96     typedef typename Graph::InEdgeIt InEdgeIt;
       
    97 
       
    98     const Graph* _graph;
       
    99 
       
   100     const CapacityMap* _capacity;
    91     const CapacityMap* _capacity;
   101 
    92 
   102     typedef typename Graph::template EdgeMap<Value> FlowMap;
    93     typedef typename Graph::template EdgeMap<Value> FlowMap;
   103 
    94     FlowMap* _flow;
   104     FlowMap* _preflow;
    95 
   105 
    96     Node _source;
   106     Node _source, _target;
    97 
   107     int _node_num;
    98     int _node_num;
   108 
    99 
   109     typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 
   100     // Bucketing structure
   110                             FlowMap, Tolerance> OutResGraph;
   101     std::vector<Node> _first, _last;
   111     typedef typename OutResGraph::Edge OutResEdge;
   102     typename Graph::template NodeMap<Node>* _next;
       
   103     typename Graph::template NodeMap<Node>* _prev;    
       
   104     typename Graph::template NodeMap<bool>* _active;
       
   105     typename Graph::template NodeMap<int>* _bucket;
   112     
   106     
   113     OutResGraph* _out_res_graph;
   107     std::vector<bool> _dormant;
   114 
   108 
   115     typedef RevGraphAdaptor<const Graph> RevGraph;
   109     std::list<std::list<int> > _sets;
   116     RevGraph* _rev_graph;
   110     std::list<int>::iterator _highest;
   117 
       
   118     typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 
       
   119                             FlowMap, Tolerance> InResGraph;
       
   120     typedef typename InResGraph::Edge InResEdge;
       
   121     
       
   122     InResGraph* _in_res_graph;
       
   123 
       
   124     typedef IterableBoolMap<Graph, Node> WakeMap;
       
   125     WakeMap* _wake;
       
   126 
       
   127     typedef typename Graph::template NodeMap<int> DistMap;
       
   128     DistMap* _dist;  
       
   129     
   111     
   130     typedef typename Graph::template NodeMap<Value> ExcessMap;
   112     typedef typename Graph::template NodeMap<Value> ExcessMap;
   131     ExcessMap* _excess;
   113     ExcessMap* _excess;
   132 
   114 
   133     typedef typename Graph::template NodeMap<bool> SourceSetMap;
   115     typedef typename Graph::template NodeMap<bool> SourceSetMap;
   134     SourceSetMap* _source_set;
   116     SourceSetMap* _source_set;
   135 
       
   136     std::vector<int> _level_size;
       
   137 
       
   138     int _highest_active;
       
   139     std::vector<std::list<Node> > _active_nodes;
       
   140 
       
   141     int _dormant_max;
       
   142     std::vector<std::list<Node> > _dormant;
       
   143 
       
   144 
   117 
   145     Value _min_cut;
   118     Value _min_cut;
   146 
   119 
   147     typedef typename Graph::template NodeMap<bool> MinCutMap;
   120     typedef typename Graph::template NodeMap<bool> MinCutMap;
   148     MinCutMap* _min_cut_map;
   121     MinCutMap* _min_cut_map;
   154     /// \brief Constructor
   127     /// \brief Constructor
   155     ///
   128     ///
   156     /// Constructor of the algorithm class. 
   129     /// Constructor of the algorithm class. 
   157     HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
   130     HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
   158              const Tolerance& tolerance = Tolerance()) :
   131              const Tolerance& tolerance = Tolerance()) :
   159       _graph(&graph), _capacity(&capacity), 
   132       _graph(graph), _capacity(&capacity), _flow(0), _source(),
   160       _preflow(0), _source(), _target(), 
   133       _node_num(), _first(), _last(), _next(0), _prev(0), 
   161       _out_res_graph(0), _rev_graph(0), _in_res_graph(0),
   134       _active(0), _bucket(0), _dormant(), _sets(), _highest(),
   162       _wake(0),_dist(0), _excess(0), _source_set(0), 
   135       _excess(0), _source_set(0), _min_cut(), _min_cut_map(0), 
   163       _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 
   136       _tolerance(tolerance) {}
   164       _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
       
   165 
   137 
   166     ~HaoOrlin() {
   138     ~HaoOrlin() {
   167       if (_min_cut_map) {
   139       if (_min_cut_map) {
   168         delete _min_cut_map;
   140         delete _min_cut_map;
   169       } 
   141       } 
   170       if (_in_res_graph) {
       
   171         delete _in_res_graph;
       
   172       }
       
   173       if (_rev_graph) {
       
   174         delete _rev_graph;
       
   175       }
       
   176       if (_out_res_graph) {
       
   177         delete _out_res_graph;
       
   178       }
       
   179       if (_source_set) {
   142       if (_source_set) {
   180         delete _source_set;
   143         delete _source_set;
   181       }
   144       }
   182       if (_excess) {
   145       if (_excess) {
   183         delete _excess;
   146         delete _excess;
   184       }
   147       }
   185       if (_dist) {
   148       if (_next) {
   186         delete _dist;
   149 	delete _next;
   187       }
   150       }
   188       if (_wake) {
   151       if (_prev) {
   189         delete _wake;
   152 	delete _prev;
   190       }
   153       }
   191       if (_preflow) {
   154       if (_active) {
   192         delete _preflow;
   155 	delete _active;
       
   156       }
       
   157       if (_bucket) {
       
   158 	delete _bucket;
       
   159       }
       
   160       if (_flow) {
       
   161         delete _flow;
   193       }
   162       }
   194     }
   163     }
   195     
   164     
   196   private:
   165   private:
       
   166 
       
   167     void activate(const Node& i) {
       
   168       _active->set(i, true);
       
   169 
       
   170       int bucket = (*_bucket)[i];
       
   171 
       
   172       if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;	    
       
   173       //unlace
       
   174       _next->set((*_prev)[i], (*_next)[i]);
       
   175       if ((*_next)[i] != INVALID) {
       
   176 	_prev->set((*_next)[i], (*_prev)[i]);
       
   177       } else {
       
   178 	_last[bucket] = (*_prev)[i];
       
   179       }
       
   180       //lace
       
   181       _next->set(i, _first[bucket]);
       
   182       _prev->set(_first[bucket], i);
       
   183       _prev->set(i, INVALID);
       
   184       _first[bucket] = i;
       
   185     }
       
   186 
       
   187     void deactivate(const Node& i) {
       
   188       _active->set(i, false);
       
   189       int bucket = (*_bucket)[i];
       
   190 
       
   191       if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
       
   192       
       
   193       //unlace
       
   194       _prev->set((*_next)[i], (*_prev)[i]);
       
   195       if ((*_prev)[i] != INVALID) {
       
   196 	_next->set((*_prev)[i], (*_next)[i]);
       
   197       } else {
       
   198 	_first[bucket] = (*_next)[i];
       
   199       }
       
   200       //lace
       
   201       _prev->set(i, _last[bucket]);
       
   202       _next->set(_last[bucket], i);
       
   203       _next->set(i, INVALID);
       
   204       _last[bucket] = i;
       
   205     }
       
   206 
       
   207     void addItem(const Node& i, int bucket) {
       
   208       (*_bucket)[i] = bucket;
       
   209       if (_last[bucket] != INVALID) {
       
   210 	_prev->set(i, _last[bucket]);
       
   211 	_next->set(_last[bucket], i);
       
   212 	_next->set(i, INVALID);
       
   213 	_last[bucket] = i;
       
   214       } else {
       
   215 	_prev->set(i, INVALID);
       
   216 	_first[bucket] = i;
       
   217 	_next->set(i, INVALID);
       
   218 	_last[bucket] = i;
       
   219       }
       
   220     }
   197     
   221     
   198     template <typename ResGraph>
   222     void findMinCutOut() {
   199     void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
   223 
   200       typedef typename ResGraph::Edge ResEdge;
   224       for (NodeIt n(_graph); n != INVALID; ++n) {
   201       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   225 	_excess->set(n, 0);
   202 
   226       }
   203       for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
   227 
   204         (*_preflow)[it] = 0;      
   228       for (EdgeIt e(_graph); e != INVALID; ++e) {
   205       }
   229 	_flow->set(e, 0);
   206       for (NodeIt it(*_graph); it != INVALID; ++it) {
   230       }
   207         (*_wake)[it] = true;
   231 
   208         (*_dist)[it] = 1;
   232       int bucket_num = 1;
   209         (*_excess)[it] = 0;
   233       
   210         (*_source_set)[it] = false;
   234       {
   211       }
   235 	typename Graph::template NodeMap<bool> reached(_graph, false);
   212 
   236 	
   213       _dormant[0].push_front(_source);
   237 	reached.set(_source, true);
   214       (*_source_set)[_source] = true;
   238 
   215       _dormant_max = 0;
   239 	bool first_set = true;
   216       (*_wake)[_source] = false;
   240 
   217 
   241 	for (NodeIt t(_graph); t != INVALID; ++t) {
   218       _level_size[0] = 1;
   242 	  if (reached[t]) continue;
   219       _level_size[1] = _node_num - 1;
   243 	  _sets.push_front(std::list<int>());
   220 
   244 	  _sets.front().push_front(bucket_num);
   221       _target = target;
   245 	  _dormant[bucket_num] = !first_set;
   222       (*_dist)[target] = 0;
   246 
   223 
   247 	  _bucket->set(t, bucket_num);
   224       for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
   248 	  _first[bucket_num] = _last[bucket_num] = t;
   225         Value delta = res_graph.rescap(it);
   249 	  _next->set(t, INVALID);
   226         (*_excess)[_source] -= delta;
   250 	  _prev->set(t, INVALID);
   227         res_graph.augment(it, delta);
   251 
   228         Node a = res_graph.target(it);
   252 	  ++bucket_num;
   229         if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
   253 	  
   230           _active_nodes[(*_dist)[a]].push_front(a);
   254 	  std::vector<Node> queue;
   231           if (_highest_active < (*_dist)[a]) {
   255 	  queue.push_back(t);
   232             _highest_active = (*_dist)[a];
   256 	  reached.set(t, true);
   233           }
   257 	  
   234         }
   258 	  while (!queue.empty()) {
   235         (*_excess)[a] += delta;
   259 	    _sets.front().push_front(bucket_num);
   236       }
   260 	    _dormant[bucket_num] = !first_set;
   237 
   261 	    _first[bucket_num] = _last[bucket_num] = INVALID;
   238 
   262 	    
   239       do {
   263 	    std::vector<Node> nqueue;
   240 	Node n;
   264 	    for (int i = 0; i < int(queue.size()); ++i) {
   241 	while ((n = findActiveNode()) != INVALID) {
   265 	      Node n = queue[i];
   242 	  for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   266 	      for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
   243             Node a = res_graph.target(e);
   267 		Node u = _graph.source(e);
   244             if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
   268 		if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
   245 	    Value delta = res_graph.rescap(e);
   269 		  reached.set(u, true);
   246 	    if (_tolerance.positive((*_excess)[n] - delta)) {
   270 		  addItem(u, bucket_num);
   247               (*_excess)[n] -= delta;
   271 		  nqueue.push_back(u);
       
   272 		}
       
   273 	      }
       
   274 	    }
       
   275 	    queue.swap(nqueue);
       
   276 	    ++bucket_num;
       
   277 	  }
       
   278 	  _sets.front().pop_front();
       
   279 	  --bucket_num;
       
   280 	  first_set = false;
       
   281 	}
       
   282 
       
   283 	_bucket->set(_source, 0);
       
   284 	_dormant[0] = true;
       
   285       }
       
   286       _source_set->set(_source, true);	  
       
   287 	  
       
   288       Node target = _last[_sets.back().back()];
       
   289       {
       
   290 	for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
       
   291 	  if (_tolerance.positive((*_capacity)[e])) {
       
   292 	    Node u = _graph.target(e);
       
   293 	    _flow->set(e, (*_capacity)[e]);
       
   294 	    _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
       
   295 	    if (!(*_active)[u] && u != _source) {
       
   296 	      activate(u);
       
   297 	    }
       
   298 	  }
       
   299 	}
       
   300 	if ((*_active)[target]) {
       
   301 	  deactivate(target);
       
   302 	}
       
   303 	
       
   304 	_highest = _sets.back().begin();
       
   305 	while (_highest != _sets.back().end() && 
       
   306 	       !(*_active)[_first[*_highest]]) {
       
   307 	  ++_highest;
       
   308 	}
       
   309       }
       
   310 
       
   311 
       
   312       while (true) {
       
   313 	while (_highest != _sets.back().end()) {
       
   314 	  Node n = _first[*_highest];
       
   315 	  Value excess = (*_excess)[n];
       
   316 	  int next_bucket = _node_num;
       
   317 
       
   318 	  int under_bucket;
       
   319 	  if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
       
   320 	    under_bucket = -1;
       
   321 	  } else {
       
   322 	    under_bucket = *(++std::list<int>::iterator(_highest));
       
   323 	  }
       
   324 
       
   325 	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
       
   326 	    Node v = _graph.target(e);
       
   327 	    if (_dormant[(*_bucket)[v]]) continue;
       
   328 	    Value rem = (*_capacity)[e] - (*_flow)[e];
       
   329 	    if (!_tolerance.positive(rem)) continue;
       
   330 	    if ((*_bucket)[v] == under_bucket) {
       
   331 	      if (!(*_active)[v] && v != target) {
       
   332 		activate(v);
       
   333 	      }
       
   334 	      if (!_tolerance.less(rem, excess)) {
       
   335 		_flow->set(e, (*_flow)[e] + excess);
       
   336 		_excess->set(v, (*_excess)[v] + excess);
       
   337 		excess = 0;
       
   338 		goto no_more_push;
       
   339 	      } else {
       
   340 		excess -= rem;
       
   341 		_excess->set(v, (*_excess)[v] + rem);
       
   342 		_flow->set(e, (*_capacity)[e]);
       
   343 	      }
       
   344 	    } else if (next_bucket > (*_bucket)[v]) {
       
   345 	      next_bucket = (*_bucket)[v];
       
   346 	    }
       
   347 	  }
       
   348 
       
   349 	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
       
   350 	    Node v = _graph.source(e);
       
   351 	    if (_dormant[(*_bucket)[v]]) continue;
       
   352 	    Value rem = (*_flow)[e];
       
   353 	    if (!_tolerance.positive(rem)) continue;
       
   354 	    if ((*_bucket)[v] == under_bucket) {
       
   355 	      if (!(*_active)[v] && v != target) {
       
   356 		activate(v);
       
   357 	      }
       
   358 	      if (!_tolerance.less(rem, excess)) {
       
   359 		_flow->set(e, (*_flow)[e] - excess);
       
   360 		_excess->set(v, (*_excess)[v] + excess);
       
   361 		excess = 0;
       
   362 		goto no_more_push;
       
   363 	      } else {
       
   364 		excess -= rem;
       
   365 		_excess->set(v, (*_excess)[v] + rem);
       
   366 		_flow->set(e, 0);
       
   367 	      }
       
   368 	    } else if (next_bucket > (*_bucket)[v]) {
       
   369 	      next_bucket = (*_bucket)[v];
       
   370 	    }
       
   371 	  }
       
   372 	  
       
   373 	no_more_push:
       
   374 	  
       
   375 	  _excess->set(n, excess);
       
   376 	  
       
   377 	  if (excess != 0) {
       
   378 	    if ((*_next)[n] == INVALID) {
       
   379 	      typename std::list<std::list<int> >::iterator new_set = 
       
   380 		_sets.insert(--_sets.end(), std::list<int>());
       
   381 	      new_set->splice(new_set->end(), _sets.back(), 
       
   382 			      _sets.back().begin(), ++_highest);
       
   383 	      for (std::list<int>::iterator it = new_set->begin();
       
   384 		   it != new_set->end(); ++it) {
       
   385 		_dormant[*it] = true;
       
   386 	      }
       
   387 	      while (_highest != _sets.back().end() && 
       
   388 		     !(*_active)[_first[*_highest]]) {
       
   389 		++_highest;
       
   390 	      }
       
   391 	    } else if (next_bucket == _node_num) {
       
   392 	      _first[(*_bucket)[n]] = (*_next)[n];
       
   393 	      _prev->set((*_next)[n], INVALID);
       
   394 	      
       
   395 	      std::list<std::list<int> >::iterator new_set = 
       
   396 		_sets.insert(--_sets.end(), std::list<int>());
       
   397 	      
       
   398 	      new_set->push_front(bucket_num);
       
   399 	      _bucket->set(n, bucket_num);
       
   400 	      _first[bucket_num] = _last[bucket_num] = n;
       
   401 	      _next->set(n, INVALID);
       
   402 	      _prev->set(n, INVALID);
       
   403 	      _dormant[bucket_num] = true;
       
   404 	      ++bucket_num;
       
   405 
       
   406 	      while (_highest != _sets.back().end() && 
       
   407 		     !(*_active)[_first[*_highest]]) {
       
   408 		++_highest;
       
   409 	      }
   248 	    } else {
   410 	    } else {
   249 	      delta = (*_excess)[n];
   411 	      _first[*_highest] = (*_next)[n];
   250               (*_excess)[n] = 0;
   412 	      _prev->set((*_next)[n], INVALID);
   251 	    }
   413 	      
   252 	    res_graph.augment(e, delta);
   414 	      while (next_bucket != *_highest) {
   253 	    if ((*_excess)[a] == 0 && a != _target) {
   415 		--_highest;
   254 	      _active_nodes[(*_dist)[a]].push_front(a);
   416 	      }
   255 	    }
   417 
   256 	    (*_excess)[a] += delta;
   418 	      if (_highest == _sets.back().begin()) {
   257             if ((*_excess)[n] == 0) break;
   419 		_sets.back().push_front(bucket_num);
   258 	  }
   420 		_dormant[bucket_num] = false;
   259 	  if ((*_excess)[n] != 0) {
   421 		_first[bucket_num] = _last[bucket_num] = INVALID;
   260 	    relabel(n, res_graph);
   422 		++bucket_num;
   261           }
   423 	      }
   262 	}
   424 	      --_highest;
   263 
   425 
   264 	Value current_value = cutValue(out);
   426 	      _bucket->set(n, *_highest);
   265  	if (_min_cut > current_value){
   427 	      _next->set(n, _first[*_highest]);
   266           if (out) {
   428 	      if (_first[*_highest] != INVALID) {
   267             for (NodeIt it(*_graph); it != INVALID; ++it) {
   429 		_prev->set(_first[*_highest], n);
   268               _min_cut_map->set(it, !(*_wake)[it]);
   430 	      } else {
   269             } 
   431 		_last[*_highest] = n;
   270           } else {
   432 	      }
   271             for (NodeIt it(*_graph); it != INVALID; ++it) {
   433 	      _first[*_highest] = n;	      
   272               _min_cut_map->set(it, (*_wake)[it]);
   434 	    }
   273             } 
   435 	  } else {
   274           }
   436 
   275 
   437 	    deactivate(n);
   276 	  _min_cut = current_value;
   438 	    if (!(*_active)[_first[*_highest]]) {
   277  	}
   439 	      ++_highest;
   278 
   440 	      if (_highest != _sets.back().end() && 
   279       } while (selectNewSink(res_graph));
   441 		  !(*_active)[_first[*_highest]]) {
   280     }
   442 		_highest = _sets.back().end();
   281 
   443 	      }
   282     template <typename ResGraph>
   444 	    }
   283     void relabel(const Node& n, ResGraph& res_graph) {
   445 	  }
   284       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   446 	}
   285 
   447 
   286       int k = (*_dist)[n];
   448 	if ((*_excess)[target] < _min_cut) {
   287       if (_level_size[k] == 1) {
   449 	  _min_cut = (*_excess)[target];
   288 	++_dormant_max;
   450 	  for (NodeIt i(_graph); i != INVALID; ++i) {
   289 	for (NodeIt it(*_graph); it != INVALID; ++it) {
   451 	    _min_cut_map->set(i, true);
   290 	  if ((*_wake)[it] && (*_dist)[it] >= k) {
   452 	  }
   291 	    (*_wake)[it] = false;
   453 	  for (std::list<int>::iterator it = _sets.back().begin();
   292 	    _dormant[_dormant_max].push_front(it);
   454 	       it != _sets.back().end(); ++it) {
   293 	    --_level_size[(*_dist)[it]];
   455 	    Node n = _first[*it];
   294 	  }
   456 	    while (n != INVALID) {
   295 	}
   457 	      _min_cut_map->set(n, false);
   296         --_highest_active;
   458 	      n = (*_next)[n];
   297       } else {	
   459 	    }
   298         int new_dist = _node_num;
   460 	  }
   299         for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   461 	}
   300           Node t = res_graph.target(e);
   462 
   301           if ((*_wake)[t] && new_dist > (*_dist)[t]) {
   463 	{
   302             new_dist = (*_dist)[t];
   464 	  Node new_target;
   303           }
   465 	  if ((*_prev)[target] != INVALID) {
   304         }
   466 	    _last[(*_bucket)[target]] = (*_prev)[target];
   305         if (new_dist == _node_num) {
   467 	    _next->set((*_prev)[target], INVALID);
   306 	  ++_dormant_max;
   468 	    new_target = (*_prev)[target];
   307 	  (*_wake)[n] = false;
   469 	  } else {
   308 	  _dormant[_dormant_max].push_front(n);
   470 	    _sets.back().pop_back();
   309 	  --_level_size[(*_dist)[n]];
   471 	    if (_sets.back().empty()) {
   310 	} else {	    
   472 	      _sets.pop_back();
   311 	  --_level_size[(*_dist)[n]];
   473 	      if (_sets.empty())
   312 	  (*_dist)[n] = new_dist + 1;
   474 		break;
   313 	  _highest_active = (*_dist)[n];
   475 	      for (std::list<int>::iterator it = _sets.back().begin();
   314 	  _active_nodes[_highest_active].push_front(n);
   476 		   it != _sets.back().end(); ++it) {
   315 	  ++_level_size[(*_dist)[n]];
   477 		_dormant[*it] = false;
   316 	}
   478 	      }
   317       }
   479 	    }
   318     }
   480 	    new_target = _last[_sets.back().back()];
   319 
   481 	  }
   320     template <typename ResGraph>
   482 
   321     bool selectNewSink(ResGraph& res_graph) {
   483 	  _bucket->set(target, 0);
   322       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   484 
   323 
   485 	  _source_set->set(target, true);	  
   324       Node old_target = _target;
   486 	  for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
   325       (*_wake)[_target] = false;
   487 	    Value rem = (*_capacity)[e] - (*_flow)[e];
   326       --_level_size[(*_dist)[_target]];
   488 	    if (!_tolerance.positive(rem)) continue;
   327       _dormant[0].push_front(_target);
   489 	    Node v = _graph.target(e);
   328       (*_source_set)[_target] = true;
   490 	    if (!(*_active)[v] && !(*_source_set)[v]) {
   329       if (int(_dormant[0].size()) == _node_num){
   491 	      activate(v);
   330         _dormant[0].clear();
   492 	    }
   331 	return false;
   493 	    _excess->set(v, (*_excess)[v] + rem);
   332       }
   494 	    _flow->set(e, (*_capacity)[e]);
   333 
   495 	  }
   334       bool wake_was_empty = false;
   496 	  
   335 
   497 	  for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
   336       if(_wake->trueNum() == 0) {
   498 	    Value rem = (*_flow)[e];
   337 	while (!_dormant[_dormant_max].empty()){
   499 	    if (!_tolerance.positive(rem)) continue;
   338 	  (*_wake)[_dormant[_dormant_max].front()] = true;
   500 	    Node v = _graph.source(e);
   339 	  ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
   501 	    if (!(*_active)[v] && !(*_source_set)[v]) {
   340 	  _dormant[_dormant_max].pop_front();
   502 	      activate(v);
   341 	}
   503 	    }
   342 	--_dormant_max;
   504 	    _excess->set(v, (*_excess)[v] + rem);
   343 	wake_was_empty = true;
   505 	    _flow->set(e, 0);
   344       }
   506 	  }
   345 
   507 	  
   346       int min_dist = std::numeric_limits<int>::max();
   508 	  target = new_target;
   347       for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
   509 	  if ((*_active)[target]) {
   348 	if (min_dist > (*_dist)[it]){
   510 	    deactivate(target);
   349 	  _target = it;
   511 	  }
   350 	  min_dist = (*_dist)[it];
   512 
   351 	}
   513 	  _highest = _sets.back().begin();
   352       }
   514 	  while (_highest != _sets.back().end() && 
   353 
   515 		 !(*_active)[_first[*_highest]]) {
   354       if (wake_was_empty){
   516 	    ++_highest;
   355 	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
   517 	  }
   356           if ((*_excess)[it] != 0 && it != _target) {
   518 	}
   357             _active_nodes[(*_dist)[it]].push_front(it);
   519       }
   358             if (_highest_active < (*_dist)[it]) {
   520     }    
   359               _highest_active = (*_dist)[it];		    
   521 
   360             }
   522     void findMinCutIn() {
   361 	  }
   523 
   362 	}
   524       for (NodeIt n(_graph); n != INVALID; ++n) {
   363       }
   525 	_excess->set(n, 0);
   364 
   526       }
   365       Node n = old_target;
   527 
   366       for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
   528       for (EdgeIt e(_graph); e != INVALID; ++e) {
   367         Node a = res_graph.target(e);
   529 	_flow->set(e, 0);
   368 	if (!(*_wake)[a]) continue;
   530       }
   369         Value delta = res_graph.rescap(e);
   531 
   370         res_graph.augment(e, delta);
   532       int bucket_num = 1;
   371         (*_excess)[n] -= delta;
       
   372         if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
       
   373           _active_nodes[(*_dist)[a]].push_front(a);
       
   374           if (_highest_active < (*_dist)[a]) {
       
   375             _highest_active = (*_dist)[a];
       
   376           }
       
   377         }
       
   378         (*_excess)[a] += delta;
       
   379       }
       
   380       
   533       
   381       return true;
   534       {
   382     }
   535 	typename Graph::template NodeMap<bool> reached(_graph, false);
   383 
   536 	
   384     Node findActiveNode() {
   537 	reached.set(_source, true);
   385       while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
   538 
   386 	--_highest_active;
   539 	bool first_set = true;
   387       }
   540 
   388       if( _highest_active > 0) {
   541 	for (NodeIt t(_graph); t != INVALID; ++t) {
   389        	Node n = _active_nodes[_highest_active].front();
   542 	  if (reached[t]) continue;
   390 	_active_nodes[_highest_active].pop_front();
   543 	  _sets.push_front(std::list<int>());
   391 	return n;
   544 	  _sets.front().push_front(bucket_num);
   392       } else {
   545 	  _dormant[bucket_num] = !first_set;
   393 	return INVALID;
   546 
   394       }
   547 	  _bucket->set(t, bucket_num);
   395     }
   548 	  _first[bucket_num] = _last[bucket_num] = t;
   396 
   549 	  _next->set(t, INVALID);
   397     Value cutValue(bool out) {
   550 	  _prev->set(t, INVALID);
   398       Value value = 0;
   551 
   399       if (out) {
   552 	  ++bucket_num;
   400         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
   553 	  
   401           for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
   554 	  std::vector<Node> queue;
   402             if (!(*_wake)[_graph->source(e)]){
   555 	  queue.push_back(t);
   403               value += (*_capacity)[e];
   556 	  reached.set(t, true);
   404             }	
   557 	  
   405           }
   558 	  while (!queue.empty()) {
   406         }
   559 	    _sets.front().push_front(bucket_num);
   407       } else {
   560 	    _dormant[bucket_num] = !first_set;
   408         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
   561 	    _first[bucket_num] = _last[bucket_num] = INVALID;
   409           for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
   562 	    
   410             if (!(*_wake)[_graph->target(e)]){
   563 	    std::vector<Node> nqueue;
   411               value += (*_capacity)[e];
   564 	    for (int i = 0; i < int(queue.size()); ++i) {
   412             }	
   565 	      Node n = queue[i];
   413           }
   566 	      for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
   414         }
   567 		Node u = _graph.target(e);
   415       }
   568 		if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
   416       return value;
   569 		  reached.set(u, true);
   417     }
   570 		  addItem(u, bucket_num);
   418 
   571 		  nqueue.push_back(u);
       
   572 		}
       
   573 	      }
       
   574 	    }
       
   575 	    queue.swap(nqueue);
       
   576 	    ++bucket_num;
       
   577 	  }
       
   578 	  _sets.front().pop_front();
       
   579 	  --bucket_num;
       
   580 	  first_set = false;
       
   581 	}
       
   582 
       
   583 	_bucket->set(_source, 0);
       
   584 	_dormant[0] = true;
       
   585       }
       
   586       _source_set->set(_source, true);	  
       
   587 	  
       
   588       Node target = _last[_sets.back().back()];
       
   589       {
       
   590 	for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
       
   591 	  if (_tolerance.positive((*_capacity)[e])) {
       
   592 	    Node u = _graph.source(e);
       
   593 	    _flow->set(e, (*_capacity)[e]);
       
   594 	    _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
       
   595 	    if (!(*_active)[u] && u != _source) {
       
   596 	      activate(u);
       
   597 	    }
       
   598 	  }
       
   599 	}
       
   600 	if ((*_active)[target]) {
       
   601 	  deactivate(target);
       
   602 	}
       
   603 	
       
   604 	_highest = _sets.back().begin();
       
   605 	while (_highest != _sets.back().end() && 
       
   606 	       !(*_active)[_first[*_highest]]) {
       
   607 	  ++_highest;
       
   608 	}
       
   609       }
       
   610 
       
   611 
       
   612       while (true) {
       
   613 	while (_highest != _sets.back().end()) {
       
   614 	  Node n = _first[*_highest];
       
   615 	  Value excess = (*_excess)[n];
       
   616 	  int next_bucket = _node_num;
       
   617 
       
   618 	  int under_bucket;
       
   619 	  if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
       
   620 	    under_bucket = -1;
       
   621 	  } else {
       
   622 	    under_bucket = *(++std::list<int>::iterator(_highest));
       
   623 	  }
       
   624 
       
   625 	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
       
   626 	    Node v = _graph.source(e);
       
   627 	    if (_dormant[(*_bucket)[v]]) continue;
       
   628 	    Value rem = (*_capacity)[e] - (*_flow)[e];
       
   629 	    if (!_tolerance.positive(rem)) continue;
       
   630 	    if ((*_bucket)[v] == under_bucket) {
       
   631 	      if (!(*_active)[v] && v != target) {
       
   632 		activate(v);
       
   633 	      }
       
   634 	      if (!_tolerance.less(rem, excess)) {
       
   635 		_flow->set(e, (*_flow)[e] + excess);
       
   636 		_excess->set(v, (*_excess)[v] + excess);
       
   637 		excess = 0;
       
   638 		goto no_more_push;
       
   639 	      } else {
       
   640 		excess -= rem;
       
   641 		_excess->set(v, (*_excess)[v] + rem);
       
   642 		_flow->set(e, (*_capacity)[e]);
       
   643 	      }
       
   644 	    } else if (next_bucket > (*_bucket)[v]) {
       
   645 	      next_bucket = (*_bucket)[v];
       
   646 	    }
       
   647 	  }
       
   648 
       
   649 	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
       
   650 	    Node v = _graph.target(e);
       
   651 	    if (_dormant[(*_bucket)[v]]) continue;
       
   652 	    Value rem = (*_flow)[e];
       
   653 	    if (!_tolerance.positive(rem)) continue;
       
   654 	    if ((*_bucket)[v] == under_bucket) {
       
   655 	      if (!(*_active)[v] && v != target) {
       
   656 		activate(v);
       
   657 	      }
       
   658 	      if (!_tolerance.less(rem, excess)) {
       
   659 		_flow->set(e, (*_flow)[e] - excess);
       
   660 		_excess->set(v, (*_excess)[v] + excess);
       
   661 		excess = 0;
       
   662 		goto no_more_push;
       
   663 	      } else {
       
   664 		excess -= rem;
       
   665 		_excess->set(v, (*_excess)[v] + rem);
       
   666 		_flow->set(e, 0);
       
   667 	      }
       
   668 	    } else if (next_bucket > (*_bucket)[v]) {
       
   669 	      next_bucket = (*_bucket)[v];
       
   670 	    }
       
   671 	  }
       
   672 	  
       
   673 	no_more_push:
       
   674 	  
       
   675 	  _excess->set(n, excess);
       
   676 	  
       
   677 	  if (excess != 0) {
       
   678 	    if ((*_next)[n] == INVALID) {
       
   679 	      typename std::list<std::list<int> >::iterator new_set = 
       
   680 		_sets.insert(--_sets.end(), std::list<int>());
       
   681 	      new_set->splice(new_set->end(), _sets.back(), 
       
   682 			      _sets.back().begin(), ++_highest);
       
   683 	      for (std::list<int>::iterator it = new_set->begin();
       
   684 		   it != new_set->end(); ++it) {
       
   685 		_dormant[*it] = true;
       
   686 	      }
       
   687 	      while (_highest != _sets.back().end() && 
       
   688 		     !(*_active)[_first[*_highest]]) {
       
   689 		++_highest;
       
   690 	      }
       
   691 	    } else if (next_bucket == _node_num) {
       
   692 	      _first[(*_bucket)[n]] = (*_next)[n];
       
   693 	      _prev->set((*_next)[n], INVALID);
       
   694 	      
       
   695 	      std::list<std::list<int> >::iterator new_set = 
       
   696 		_sets.insert(--_sets.end(), std::list<int>());
       
   697 	      
       
   698 	      new_set->push_front(bucket_num);
       
   699 	      _bucket->set(n, bucket_num);
       
   700 	      _first[bucket_num] = _last[bucket_num] = n;
       
   701 	      _next->set(n, INVALID);
       
   702 	      _prev->set(n, INVALID);
       
   703 	      _dormant[bucket_num] = true;
       
   704 	      ++bucket_num;
       
   705 
       
   706 	      while (_highest != _sets.back().end() && 
       
   707 		     !(*_active)[_first[*_highest]]) {
       
   708 		++_highest;
       
   709 	      }
       
   710 	    } else {
       
   711 	      _first[*_highest] = (*_next)[n];
       
   712 	      _prev->set((*_next)[n], INVALID);
       
   713 
       
   714 	      while (next_bucket != *_highest) {
       
   715 		--_highest;
       
   716 	      }
       
   717 	      if (_highest == _sets.back().begin()) {
       
   718 		_sets.back().push_front(bucket_num);
       
   719 		_dormant[bucket_num] = false;
       
   720 		_first[bucket_num] = _last[bucket_num] = INVALID;
       
   721 		++bucket_num;
       
   722 	      }
       
   723 	      --_highest;
       
   724 
       
   725 	      _bucket->set(n, *_highest);
       
   726 	      _next->set(n, _first[*_highest]);
       
   727 	      if (_first[*_highest] != INVALID) {
       
   728 		_prev->set(_first[*_highest], n);
       
   729 	      } else {
       
   730 		_last[*_highest] = n;
       
   731 	      }
       
   732 	      _first[*_highest] = n;	      
       
   733 	    }
       
   734 	  } else {
       
   735 
       
   736 	    deactivate(n);
       
   737 	    if (!(*_active)[_first[*_highest]]) {
       
   738 	      ++_highest;
       
   739 	      if (_highest != _sets.back().end() && 
       
   740 		  !(*_active)[_first[*_highest]]) {
       
   741 		_highest = _sets.back().end();
       
   742 	      }
       
   743 	    }
       
   744 	  }
       
   745 	}
       
   746 
       
   747 	if ((*_excess)[target] < _min_cut) {
       
   748 	  _min_cut = (*_excess)[target];
       
   749 	  for (NodeIt i(_graph); i != INVALID; ++i) {
       
   750 	    _min_cut_map->set(i, false);
       
   751 	  }
       
   752 	  for (std::list<int>::iterator it = _sets.back().begin();
       
   753 	       it != _sets.back().end(); ++it) {
       
   754 	    Node n = _first[*it];
       
   755 	    while (n != INVALID) {
       
   756 	      _min_cut_map->set(n, true);
       
   757 	      n = (*_next)[n];
       
   758 	    }
       
   759 	  }
       
   760 	}
       
   761 
       
   762 	{
       
   763 	  Node new_target;
       
   764 	  if ((*_prev)[target] != INVALID) {
       
   765 	    _last[(*_bucket)[target]] = (*_prev)[target];
       
   766 	    _next->set((*_prev)[target], INVALID);
       
   767 	    new_target = (*_prev)[target];
       
   768 	  } else {
       
   769 	    _sets.back().pop_back();
       
   770 	    if (_sets.back().empty()) {
       
   771 	      _sets.pop_back();
       
   772 	      if (_sets.empty())
       
   773 		break;
       
   774 	      for (std::list<int>::iterator it = _sets.back().begin();
       
   775 		   it != _sets.back().end(); ++it) {
       
   776 		_dormant[*it] = false;
       
   777 	      }
       
   778 	    }
       
   779 	    new_target = _last[_sets.back().back()];
       
   780 	  }
       
   781 
       
   782 	  _bucket->set(target, 0);
       
   783 
       
   784 	  _source_set->set(target, true);	  
       
   785 	  for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
       
   786 	    Value rem = (*_capacity)[e] - (*_flow)[e];
       
   787 	    if (!_tolerance.positive(rem)) continue;
       
   788 	    Node v = _graph.source(e);
       
   789 	    if (!(*_active)[v] && !(*_source_set)[v]) {
       
   790 	      activate(v);
       
   791 	    }
       
   792 	    _excess->set(v, (*_excess)[v] + rem);
       
   793 	    _flow->set(e, (*_capacity)[e]);
       
   794 	  }
       
   795 	  
       
   796 	  for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
       
   797 	    Value rem = (*_flow)[e];
       
   798 	    if (!_tolerance.positive(rem)) continue;
       
   799 	    Node v = _graph.target(e);
       
   800 	    if (!(*_active)[v] && !(*_source_set)[v]) {
       
   801 	      activate(v);
       
   802 	    }
       
   803 	    _excess->set(v, (*_excess)[v] + rem);
       
   804 	    _flow->set(e, 0);
       
   805 	  }
       
   806 	  
       
   807 	  target = new_target;
       
   808 	  if ((*_active)[target]) {
       
   809 	    deactivate(target);
       
   810 	  }
       
   811 
       
   812 	  _highest = _sets.back().begin();
       
   813 	  while (_highest != _sets.back().end() && 
       
   814 		 !(*_active)[_first[*_highest]]) {
       
   815 	    ++_highest;
       
   816 	  }
       
   817 	}
       
   818       }
       
   819     }
   419 
   820 
   420   public:
   821   public:
   421 
   822 
   422     /// \name Execution control
   823     /// \name Execution control
   423     /// The simplest way to execute the algorithm is to use
   824     /// The simplest way to execute the algorithm is to use
   433     ///
   834     ///
   434     /// Initializes the internal data structures. It creates
   835     /// Initializes the internal data structures. It creates
   435     /// the maps, residual graph adaptors and some bucket structures
   836     /// the maps, residual graph adaptors and some bucket structures
   436     /// for the algorithm. 
   837     /// for the algorithm. 
   437     void init() {
   838     void init() {
   438       init(NodeIt(*_graph));
   839       init(NodeIt(_graph));
   439     }
   840     }
   440 
   841 
   441     /// \brief Initializes the internal data structures.
   842     /// \brief Initializes the internal data structures.
   442     ///
   843     ///
   443     /// Initializes the internal data structures. It creates
   844     /// Initializes the internal data structures. It creates
   444     /// the maps, residual graph adaptor and some bucket structures
   845     /// the maps, residual graph adaptor and some bucket structures
   445     /// for the algorithm. Node \c source  is used as the push-relabel
   846     /// for the algorithm. Node \c source  is used as the push-relabel
   446     /// algorithm's source.
   847     /// algorithm's source.
   447     void init(const Node& source) {
   848     void init(const Node& source) {
   448       _source = source;
   849       _source = source;
   449       _node_num = countNodes(*_graph);
   850       
       
   851       _node_num = countNodes(_graph);
       
   852       
       
   853       _first.resize(_node_num);
       
   854       _last.resize(_node_num);
   450 
   855 
   451       _dormant.resize(_node_num);
   856       _dormant.resize(_node_num);
   452       _level_size.resize(_node_num, 0);
   857 
   453       _active_nodes.resize(_node_num);
   858       if (!_flow) {
   454 
   859 	_flow = new FlowMap(_graph);
   455       if (!_preflow) {
   860       }
   456         _preflow = new FlowMap(*_graph);
   861       if (!_next) {
   457       }
   862 	_next = new typename Graph::template NodeMap<Node>(_graph);
   458       if (!_wake) {
   863       }
   459         _wake = new WakeMap(*_graph);
   864       if (!_prev) {
   460       }
   865 	_prev = new typename Graph::template NodeMap<Node>(_graph);
   461       if (!_dist) {
   866       }
   462         _dist = new DistMap(*_graph);
   867       if (!_active) {
       
   868 	_active = new typename Graph::template NodeMap<bool>(_graph);
       
   869       }
       
   870       if (!_bucket) {
       
   871 	_bucket = new typename Graph::template NodeMap<int>(_graph);
   463       }
   872       }
   464       if (!_excess) {
   873       if (!_excess) {
   465         _excess = new ExcessMap(*_graph);
   874 	_excess = new ExcessMap(_graph);
   466       }
   875       }
   467       if (!_source_set) {
   876       if (!_source_set) {
   468         _source_set = new SourceSetMap(*_graph);
   877 	_source_set = new SourceSetMap(_graph);
   469       }
       
   470       if (!_out_res_graph) {
       
   471         _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
       
   472       }
       
   473       if (!_rev_graph) {
       
   474         _rev_graph = new RevGraph(*_graph);
       
   475       }
       
   476       if (!_in_res_graph) {
       
   477         _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
       
   478       }
   878       }
   479       if (!_min_cut_map) {
   879       if (!_min_cut_map) {
   480         _min_cut_map = new MinCutMap(*_graph);
   880 	_min_cut_map = new MinCutMap(_graph);
   481       }
   881       }
   482 
   882 
   483       _min_cut = std::numeric_limits<Value>::max();
   883       _min_cut = std::numeric_limits<Value>::max();
   484     }
   884     }
   485 
   885 
   486 
   886 
   487     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   887     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   488     /// source-side.
   888     /// source-side.
   489     ///
   889     ///
       
   890     /// Calculates a minimum cut with \f$ source \f$ on the
       
   891     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
       
   892     /// \in X \f$ and minimal out-degree).
       
   893     void calculateOut() {
       
   894       findMinCutOut();
       
   895     }
       
   896 
   490     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   897     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   491     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
   898     /// target-side.
   492     ///  and minimal out-degree).
       
   493     void calculateOut() {
       
   494       for (NodeIt it(*_graph); it != INVALID; ++it) {
       
   495         if (it != _source) {
       
   496           calculateOut(it);
       
   497           return;
       
   498         }
       
   499       }
       
   500     }
       
   501 
       
   502     /// \brief Calculates a minimum cut with \f$ source \f$ on the
       
   503     /// source-side.
       
   504     ///
   899     ///
   505     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   900     /// Calculates a minimum cut with \f$ source \f$ on the
   506     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
   901     /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
   507     ///  and minimal out-degree). The \c target is the initial target
   902     /// \in X \f$ and minimal out-degree).
   508     /// for the push-relabel algorithm.
       
   509     void calculateOut(const Node& target) {
       
   510       findMinCut(target, true, *_out_res_graph);
       
   511     }
       
   512 
       
   513     /// \brief Calculates a minimum cut with \f$ source \f$ on the
       
   514     /// sink-side.
       
   515     ///
       
   516     /// \brief Calculates a minimum cut with \f$ source \f$ on the
       
   517     /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with 
       
   518     /// \f$ source \notin X \f$
       
   519     /// and minimal out-degree).
       
   520     void calculateIn() {
   903     void calculateIn() {
   521       for (NodeIt it(*_graph); it != INVALID; ++it) {
   904       findMinCutIn();
   522         if (it != _source) {
   905     }
   523           calculateIn(it);
   906 
   524           return;
       
   525         }
       
   526       }
       
   527     }
       
   528 
       
   529     /// \brief Calculates a minimum cut with \f$ source \f$ on the
       
   530     /// sink-side.
       
   531     ///
       
   532     /// \brief Calculates a minimum cut with \f$ source \f$ on the
       
   533     /// sink-side (i.e. a set \f$ X\subsetneq V 
       
   534     /// \f$ with \f$ source \notin X \f$ and minimal out-degree).  
       
   535     /// The \c target is the initial
       
   536     /// target for the push-relabel algorithm.
       
   537     void calculateIn(const Node& target) {
       
   538       findMinCut(target, false, *_in_res_graph);
       
   539     }
       
   540 
   907 
   541     /// \brief Runs the algorithm.
   908     /// \brief Runs the algorithm.
   542     ///
   909     ///
   543     /// Runs the algorithm. It finds nodes \c source and \c target
   910     /// Runs the algorithm. It finds nodes \c source and \c target
   544     /// arbitrarily and then calls \ref init(), \ref calculateOut()
   911     /// arbitrarily and then calls \ref init(), \ref calculateOut()
   545     /// and \ref calculateIn().
   912     /// and \ref calculateIn().
   546     void run() {
   913     void run() {
   547       init();
   914       init();
   548       for (NodeIt it(*_graph); it != INVALID; ++it) {
   915       calculateOut();
   549         if (it != _source) {
   916       calculateIn();
   550           calculateOut(it);
       
   551           calculateIn(it);
       
   552           return;
       
   553         }
       
   554       }
       
   555     }
   917     }
   556 
   918 
   557     /// \brief Runs the algorithm.
   919     /// \brief Runs the algorithm.
   558     ///
   920     ///
   559     /// Runs the algorithm. It uses the given \c source node, finds a
   921     /// Runs the algorithm. It uses the given \c source node, finds a
   560     /// proper \c target and then calls the \ref init(), \ref
   922     /// proper \c target and then calls the \ref init(), \ref
   561     /// calculateOut() and \ref calculateIn().
   923     /// calculateOut() and \ref calculateIn().
   562     void run(const Node& s) {
   924     void run(const Node& s) {
   563       init(s);
   925       init(s);
   564       for (NodeIt it(*_graph); it != INVALID; ++it) {
   926       calculateOut();
   565         if (it != _source) {
   927       calculateIn();
   566           calculateOut(it);
       
   567           calculateIn(it);
       
   568           return;
       
   569         }
       
   570       }
       
   571     }
       
   572 
       
   573     /// \brief Runs the algorithm.
       
   574     ///
       
   575     /// Runs the algorithm. It just calls the \ref init() and then
       
   576     /// \ref calculateOut() and \ref calculateIn().
       
   577     void run(const Node& s, const Node& t) {
       
   578       init(s); 
       
   579       calculateOut(t);
       
   580       calculateIn(t);
       
   581     }
   928     }
   582 
   929 
   583     /// @}
   930     /// @}
   584     
   931     
   585     /// \name Query Functions 
   932     /// \name Query Functions 
   606     /// with minimal out-degree (i.e. \c nodeMap will be true exactly
   953     /// with minimal out-degree (i.e. \c nodeMap will be true exactly
   607     /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
   954     /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
   608     /// bool-valued node-map.
   955     /// bool-valued node-map.
   609     template <typename NodeMap>
   956     template <typename NodeMap>
   610     Value minCut(NodeMap& nodeMap) const {
   957     Value minCut(NodeMap& nodeMap) const {
   611       for (NodeIt it(*_graph); it != INVALID; ++it) {
   958       for (NodeIt it(_graph); it != INVALID; ++it) {
   612 	nodeMap.set(it, (*_min_cut_map)[it]);
   959 	nodeMap.set(it, (*_min_cut_map)[it]);
   613       }
   960       }
   614       return minCut();
   961       return minCut();
   615     }
   962     }
   616 
   963