Changeset 2530:f86f7e4eb2ba in lemon-0.x

Timestamp:

12/04/07 11:55:27 (15 years ago)

Author:

Balazs Dezso

Branch:

default

Phase:

public

Convert:

svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3407

Message:

Reimplementation of Hao-Orlin algorithm
Little modifictaion in NagamochiIbaraki?
More docs for minimum cut algorithms

Files:

• doc/groups.dox (modified) (1 diff)
• lemon/hao_orlin.h (modified) (15 diffs)
• lemon/nagamochi_ibaraki.h (modified) (4 diffs)

doc/groups.dox

@@ -259 +259 @@
 
 /**
-@defgroup min_cut Minimum Cut algorithms
-@ingroup algs
-\brief This group describes the algorithms
-for finding minimum cut in graphs.
-
-This group describes the algorithms
-for finding minimum cut in graphs.
+@defgroup min_cut Minimum Cut algorithms 
+@ingroup algs 
+
+\brief This group describes the algorithms for finding minimum cut in
+graphs.
+
+This group describes the algorithms for finding minimum cut in graphs.
+
+The minimum cut problem is to find a non-empty and non-complete
+\f$X\f$ subset of the vertices with minimum overall capacity on
+outgoing arcs. Formally, there is \f$G=(V,A)\f$ directed graph, an
+\f$c_a:A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
+cut is the solution of the next optimization problem:
+
+\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}\sum_{uv\in A, u\in X, v\not\in X}c_{uv}\f]
+
+The lemon contains several algorithms related to minimum cut problems:
+
+- \ref lemon::HaoOrlin "Hao-Orlin algorithm" for calculate minimum cut
+  in directed graphs  
+- \ref lemon::NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
+  calculate minimum cut in undirected graphs
+- \ref lemon::GomoryHuTree "Gomory-Hu tree computation" for calculate all
+  pairs minimum cut in undirected graphs
+
+If you want to find minimum cut just between two distinict nodes,
+please see the \ref max_flow "Maximum Flow page".
+
 */
 

lemon/hao_orlin.h

@@ -20 +20 @@
 #define LEMON_HAO_ORLIN_H
 
-#include <cassert>
- 
-
-
 #include <vector>
-#include <queue>
 #include <list>
+#include <ext/hash_set>
 #include <limits>
 
@@ -38 +34 @@
 /// \brief Implementation of the Hao-Orlin algorithm.
 ///
-/// Implementation of the HaoOrlin algorithms class for testing network 
+/// Implementation of the Hao-Orlin algorithm class for testing network 
 /// reliability.
 
@@ -47 +43 @@
   /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
   ///
-  /// Hao-Orlin calculates a minimum cut in a directed graph 
-  /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
-  /// of two phases: in the first phase it determines a minimum cut
-  /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
-  /// with \f$ source \in X \f$ and minimal out-degree) and in the
-  /// second phase it determines a minimum cut with \f$ source \f$ on the
-  /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
-  /// and minimal out-degree). Obviously, the smaller of these two
-  /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
-  /// modified push-relabel preflow algorithm and our implementation
-  /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
-  /// highest-label rule). The purpose of such an algorithm is testing
-  /// network reliability. For an undirected graph with \f$ n \f$
-  /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
-  /// and Ibaraki which solves the undirected problem in 
-  /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut 
-  /// algorithm
-  /// class.
+  /// Hao-Orlin calculates a minimum cut in a directed graph
+  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
+  /// consists of two phases: in the first phase it determines a
+  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
+  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
+  /// out-degree) and in the second phase it determines a minimum cut
+  /// with \f$ source \f$ on the sink-side (i.e. a set 
+  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
+  /// out-degree). Obviously, the smaller of these two cuts will be a
+  /// minimum cut of \f$ D \f$. The algorithm is a modified
+  /// push-relabel preflow algorithm and our implementation calculates
+  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
+  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
+  /// purpose of such algorithm is testing network reliability. For an
+  /// undirected graph you can run just the first phase of the
+  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
+  /// which solves the undirected problem in 
+  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
+  /// NagamochiIbaraki algorithm class.
   ///
   /// \param _Graph is the graph type of the algorithm.
@@ -81 +78 @@
 #endif
   class HaoOrlin {
-  protected:
+  private:
 
     typedef _Graph Graph;
@@ -89 +86 @@
     typedef typename CapacityMap::Value Value;
 
+    GRAPH_TYPEDEFS(typename Graph);
     
-    typedef typename Graph::Node Node;
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::EdgeIt EdgeIt;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
-    typedef typename Graph::InEdgeIt InEdgeIt;
-
-    const Graph* _graph;
-
+    const Graph& _graph;
     const CapacityMap* _capacity;
 
     typedef typename Graph::template EdgeMap<Value> FlowMap;
-
-    FlowMap* _preflow;
-
-    Node _source, _target;
+    FlowMap* _flow;
+
+    Node _source;
+
     int _node_num;
 
-    typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 
-                            FlowMap, Tolerance> OutResGraph;
-    typedef typename OutResGraph::Edge OutResEdge;
+    // Bucketing structure
+    std::vector<Node> _first, _last;
+    typename Graph::template NodeMap<Node>* _next;
+    typename Graph::template NodeMap<Node>* _prev;    
+    typename Graph::template NodeMap<bool>* _active;
+    typename Graph::template NodeMap<int>* _bucket;
     
-    OutResGraph* _out_res_graph;
-
-    typedef RevGraphAdaptor<const Graph> RevGraph;
-    RevGraph* _rev_graph;
-
-    typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 
-                            FlowMap, Tolerance> InResGraph;
-    typedef typename InResGraph::Edge InResEdge;
-    
-    InResGraph* _in_res_graph;
-
-    typedef IterableBoolMap<Graph, Node> WakeMap;
-    WakeMap* _wake;
-
-    typedef typename Graph::template NodeMap<int> DistMap;
-    DistMap* _dist;  
+    std::vector<bool> _dormant;
+
+    std::list<std::list<int> > _sets;
+    std::list<int>::iterator _highest;
     
     typedef typename Graph::template NodeMap<Value> ExcessMap;
@@ -133 +115 @@
     typedef typename Graph::template NodeMap<bool> SourceSetMap;
     SourceSetMap* _source_set;
-
-    std::vector<int> _level_size;
-
-    int _highest_active;
-    std::vector<std::list<Node> > _active_nodes;
-
-    int _dormant_max;
-    std::vector<std::list<Node> > _dormant;
-
 
     Value _min_cut;
@@ -157 +130 @@
     HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
              const Tolerance& tolerance = Tolerance()) :
-      _graph(&graph), _capacity(&capacity), 
-      _preflow(0), _source(), _target(), 
-      _out_res_graph(0), _rev_graph(0), _in_res_graph(0),
-      _wake(0),_dist(0), _excess(0), _source_set(0), 
-      _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 
-      _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
+      _graph(graph), _capacity(&capacity), _flow(0), _source(),
+      _node_num(), _first(), _last(), _next(0), _prev(0), 
+      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
+      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0), 
+      _tolerance(tolerance) {}
 
     ~HaoOrlin() {
@@ -168 +140 @@
         delete _min_cut_map;
       } 
-      if (_in_res_graph) {
-        delete _in_res_graph;
-      }
-      if (_rev_graph) {
-        delete _rev_graph;
-      }
-      if (_out_res_graph) {
-        delete _out_res_graph;
-      }
       if (_source_set) {
         delete _source_set;
@@ -183 +146 @@
         delete _excess;
       }
-      if (_dist) {
-        delete _dist;
-      }
-      if (_wake) {
-        delete _wake;
-      }
-      if (_preflow) {
-        delete _preflow;
+      if (_next) {
+        delete _next;
+      }
+      if (_prev) {
+        delete _prev;
+      }
+      if (_active) {
+        delete _active;
+      }
+      if (_bucket) {
+        delete _bucket;
+      }
+      if (_flow) {
+        delete _flow;
       }
     }
     
   private:
+
+    void activate(const Node& i) {
+      _active->set(i, true);
+
+      int bucket = (*_bucket)[i];
+
+      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;        
+      //unlace
+      _next->set((*_prev)[i], (*_next)[i]);
+      if ((*_next)[i] != INVALID) {
+        _prev->set((*_next)[i], (*_prev)[i]);
+      } else {
+        _last[bucket] = (*_prev)[i];
+      }
+      //lace
+      _next->set(i, _first[bucket]);
+      _prev->set(_first[bucket], i);
+      _prev->set(i, INVALID);
+      _first[bucket] = i;
+    }
+
+    void deactivate(const Node& i) {
+      _active->set(i, false);
+      int bucket = (*_bucket)[i];
+
+      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
+      
+      //unlace
+      _prev->set((*_next)[i], (*_prev)[i]);
+      if ((*_prev)[i] != INVALID) {
+        _next->set((*_prev)[i], (*_next)[i]);
+      } else {
+        _first[bucket] = (*_next)[i];
+      }
+      //lace
+      _prev->set(i, _last[bucket]);
+      _next->set(_last[bucket], i);
+      _next->set(i, INVALID);
+      _last[bucket] = i;
+    }
+
+    void addItem(const Node& i, int bucket) {
+      (*_bucket)[i] = bucket;
+      if (_last[bucket] != INVALID) {
+        _prev->set(i, _last[bucket]);
+        _next->set(_last[bucket], i);
+        _next->set(i, INVALID);
+        _last[bucket] = i;
+      } else {
+        _prev->set(i, INVALID);
+        _first[bucket] = i;
+        _next->set(i, INVALID);
+        _last[bucket] = i;
+      }
+    }
     
-    template <typename ResGraph>
-    void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
-      typedef typename ResGraph::Edge ResEdge;
-      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
-
-      for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
-        (*_preflow)[it] = 0;      
-      }
-      for (NodeIt it(*_graph); it != INVALID; ++it) {
-        (*_wake)[it] = true;
-        (*_dist)[it] = 1;
-        (*_excess)[it] = 0;
-        (*_source_set)[it] = false;
-      }
-
-      _dormant[0].push_front(_source);
-      (*_source_set)[_source] = true;
-      _dormant_max = 0;
-      (*_wake)[_source] = false;
-
-      _level_size[0] = 1;
-      _level_size[1] = _node_num - 1;
-
-      _target = target;
-      (*_dist)[target] = 0;
-
-      for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
-        Value delta = res_graph.rescap(it);
-        (*_excess)[_source] -= delta;
-        res_graph.augment(it, delta);
-        Node a = res_graph.target(it);
-        if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
-          _active_nodes[(*_dist)[a]].push_front(a);
-          if (_highest_active < (*_dist)[a]) {
-            _highest_active = (*_dist)[a];
-          }
-        }
-        (*_excess)[a] += delta;
-      }
-
-
-      do {
-        Node n;
-        while ((n = findActiveNode()) != INVALID) {
-          for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
-            Node a = res_graph.target(e);
-            if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
-            Value delta = res_graph.rescap(e);
-            if (_tolerance.positive((*_excess)[n] - delta)) {
-              (*_excess)[n] -= delta;
+    void findMinCutOut() {
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _excess->set(n, 0);
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow->set(e, 0);
+      }
+
+      int bucket_num = 1;
+      
+      {
+        typename Graph::template NodeMap<bool> reached(_graph, false);
+        
+        reached.set(_source, true);
+
+        bool first_set = true;
+
+        for (NodeIt t(_graph); t != INVALID; ++t) {
+          if (reached[t]) continue;
+          _sets.push_front(std::list<int>());
+          _sets.front().push_front(bucket_num);
+          _dormant[bucket_num] = !first_set;
+
+          _bucket->set(t, bucket_num);
+          _first[bucket_num] = _last[bucket_num] = t;
+          _next->set(t, INVALID);
+          _prev->set(t, INVALID);
+
+          ++bucket_num;
+          
+          std::vector<Node> queue;
+          queue.push_back(t);
+          reached.set(t, true);
+          
+          while (!queue.empty()) {
+            _sets.front().push_front(bucket_num);
+            _dormant[bucket_num] = !first_set;
+            _first[bucket_num] = _last[bucket_num] = INVALID;
+            
+            std::vector<Node> nqueue;
+            for (int i = 0; i < int(queue.size()); ++i) {
+              Node n = queue[i];
+              for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+                Node u = _graph.source(e);
+                if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+                  reached.set(u, true);
+                  addItem(u, bucket_num);
+                  nqueue.push_back(u);
+                }
+              }
+            }
+            queue.swap(nqueue);
+            ++bucket_num;
+          }
+          _sets.front().pop_front();
+          --bucket_num;
+          first_set = false;
+        }
+
+        _bucket->set(_source, 0);
+        _dormant[0] = true;
+      }
+      _source_set->set(_source, true);    
+          
+      Node target = _last[_sets.back().back()];
+      {
+        for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+          if (_tolerance.positive((*_capacity)[e])) {
+            Node u = _graph.target(e);
+            _flow->set(e, (*_capacity)[e]);
+            _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+            if (!(*_active)[u] && u != _source) {
+              activate(u);
+            }
+          }
+        }
+        if ((*_active)[target]) {
+          deactivate(target);
+        }
+        
+        _highest = _sets.back().begin();
+        while (_highest != _sets.back().end() && 
+               !(*_active)[_first[*_highest]]) {
+          ++_highest;
+        }
+      }
+
+
+      while (true) {
+        while (_highest != _sets.back().end()) {
+          Node n = _first[*_highest];
+          Value excess = (*_excess)[n];
+          int next_bucket = _node_num;
+
+          int under_bucket;
+          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
+            under_bucket = -1;
+          } else {
+            under_bucket = *(++std::list<int>::iterator(_highest));
+          }
+
+          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node v = _graph.target(e);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] + excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, (*_capacity)[e]);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node v = _graph.source(e);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] - excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, 0);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+          
+        no_more_push:
+          
+          _excess->set(n, excess);
+          
+          if (excess != 0) {
+            if ((*_next)[n] == INVALID) {
+              typename std::list<std::list<int> >::iterator new_set = 
+                _sets.insert(--_sets.end(), std::list<int>());
+              new_set->splice(new_set->end(), _sets.back(), 
+                              _sets.back().begin(), ++_highest);
+              for (std::list<int>::iterator it = new_set->begin();
+                   it != new_set->end(); ++it) {
+                _dormant[*it] = true;
+              }
+              while (_highest != _sets.back().end() && 
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else if (next_bucket == _node_num) {
+              _first[(*_bucket)[n]] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+              
+              std::list<std::list<int> >::iterator new_set = 
+                _sets.insert(--_sets.end(), std::list<int>());
+              
+              new_set->push_front(bucket_num);
+              _bucket->set(n, bucket_num);
+              _first[bucket_num] = _last[bucket_num] = n;
+              _next->set(n, INVALID);
+              _prev->set(n, INVALID);
+              _dormant[bucket_num] = true;
+              ++bucket_num;
+
+              while (_highest != _sets.back().end() && 
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
             } else {
-              delta = (*_excess)[n];
-              (*_excess)[n] = 0;
-            }
-            res_graph.augment(e, delta);
-            if ((*_excess)[a] == 0 && a != _target) {
-              _active_nodes[(*_dist)[a]].push_front(a);
-            }
-            (*_excess)[a] += delta;
-            if ((*_excess)[n] == 0) break;
-          }
-          if ((*_excess)[n] != 0) {
-            relabel(n, res_graph);
-          }
-        }
-
-        Value current_value = cutValue(out);
-        if (_min_cut > current_value){
-          if (out) {
-            for (NodeIt it(*_graph); it != INVALID; ++it) {
-              _min_cut_map->set(it, !(*_wake)[it]);
-            } 
-          } else {
-            for (NodeIt it(*_graph); it != INVALID; ++it) {
-              _min_cut_map->set(it, (*_wake)[it]);
-            } 
-          }
-
-          _min_cut = current_value;
-        }
-
-      } while (selectNewSink(res_graph));
-    }
-
-    template <typename ResGraph>
-    void relabel(const Node& n, ResGraph& res_graph) {
-      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
-
-      int k = (*_dist)[n];
-      if (_level_size[k] == 1) {
-        ++_dormant_max;
-        for (NodeIt it(*_graph); it != INVALID; ++it) {
-          if ((*_wake)[it] && (*_dist)[it] >= k) {
-            (*_wake)[it] = false;
-            _dormant[_dormant_max].push_front(it);
-            --_level_size[(*_dist)[it]];
-          }
-        }
-        --_highest_active;
-      } else {  
-        int new_dist = _node_num;
-        for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
-          Node t = res_graph.target(e);
-          if ((*_wake)[t] && new_dist > (*_dist)[t]) {
-            new_dist = (*_dist)[t];
-          }
-        }
-        if (new_dist == _node_num) {
-          ++_dormant_max;
-          (*_wake)[n] = false;
-          _dormant[_dormant_max].push_front(n);
-          --_level_size[(*_dist)[n]];
-        } else {            
-          --_level_size[(*_dist)[n]];
-          (*_dist)[n] = new_dist + 1;
-          _highest_active = (*_dist)[n];
-          _active_nodes[_highest_active].push_front(n);
-          ++_level_size[(*_dist)[n]];
-        }
-      }
-    }
-
-    template <typename ResGraph>
-    bool selectNewSink(ResGraph& res_graph) {
-      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
-
-      Node old_target = _target;
-      (*_wake)[_target] = false;
-      --_level_size[(*_dist)[_target]];
-      _dormant[0].push_front(_target);
-      (*_source_set)[_target] = true;
-      if (int(_dormant[0].size()) == _node_num){
-        _dormant[0].clear();
-        return false;
-      }
-
-      bool wake_was_empty = false;
-
-      if(_wake->trueNum() == 0) {
-        while (!_dormant[_dormant_max].empty()){
-          (*_wake)[_dormant[_dormant_max].front()] = true;
-          ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
-          _dormant[_dormant_max].pop_front();
-        }
-        --_dormant_max;
-        wake_was_empty = true;
-      }
-
-      int min_dist = std::numeric_limits<int>::max();
-      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
-        if (min_dist > (*_dist)[it]){
-          _target = it;
-          min_dist = (*_dist)[it];
-        }
-      }
-
-      if (wake_was_empty){
-        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
-          if ((*_excess)[it] != 0 && it != _target) {
-            _active_nodes[(*_dist)[it]].push_front(it);
-            if (_highest_active < (*_dist)[it]) {
-              _highest_active = (*_dist)[it];               
-            }
-          }
-        }
-      }
-
-      Node n = old_target;
-      for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
-        Node a = res_graph.target(e);
-        if (!(*_wake)[a]) continue;
-        Value delta = res_graph.rescap(e);
-        res_graph.augment(e, delta);
-        (*_excess)[n] -= delta;
-        if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
-          _active_nodes[(*_dist)[a]].push_front(a);
-          if (_highest_active < (*_dist)[a]) {
-            _highest_active = (*_dist)[a];
-          }
-        }
-        (*_excess)[a] += delta;
-      }
+              _first[*_highest] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+              
+              while (next_bucket != *_highest) {
+                --_highest;
+              }
+
+              if (_highest == _sets.back().begin()) {
+                _sets.back().push_front(bucket_num);
+                _dormant[bucket_num] = false;
+                _first[bucket_num] = _last[bucket_num] = INVALID;
+                ++bucket_num;
+              }
+              --_highest;
+
+              _bucket->set(n, *_highest);
+              _next->set(n, _first[*_highest]);
+              if (_first[*_highest] != INVALID) {
+                _prev->set(_first[*_highest], n);
+              } else {
+                _last[*_highest] = n;
+              }
+              _first[*_highest] = n;          
+            }
+          } else {
+
+            deactivate(n);
+            if (!(*_active)[_first[*_highest]]) {
+              ++_highest;
+              if (_highest != _sets.back().end() && 
+                  !(*_active)[_first[*_highest]]) {
+                _highest = _sets.back().end();
+              }
+            }
+          }
+        }
+
+        if ((*_excess)[target] < _min_cut) {
+          _min_cut = (*_excess)[target];
+          for (NodeIt i(_graph); i != INVALID; ++i) {
+            _min_cut_map->set(i, true);
+          }
+          for (std::list<int>::iterator it = _sets.back().begin();
+               it != _sets.back().end(); ++it) {
+            Node n = _first[*it];
+            while (n != INVALID) {
+              _min_cut_map->set(n, false);
+              n = (*_next)[n];
+            }
+          }
+        }
+
+        {
+          Node new_target;
+          if ((*_prev)[target] != INVALID) {
+            _last[(*_bucket)[target]] = (*_prev)[target];
+            _next->set((*_prev)[target], INVALID);
+            new_target = (*_prev)[target];
+          } else {
+            _sets.back().pop_back();
+            if (_sets.back().empty()) {
+              _sets.pop_back();
+              if (_sets.empty())
+                break;
+              for (std::list<int>::iterator it = _sets.back().begin();
+                   it != _sets.back().end(); ++it) {
+                _dormant[*it] = false;
+              }
+            }
+            new_target = _last[_sets.back().back()];
+          }
+
+          _bucket->set(target, 0);
+
+          _source_set->set(target, true);         
+          for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.target(e);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(e, (*_capacity)[e]);
+          }
+          
+          for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
+            Value rem = (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.source(e);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(e, 0);
+          }
+          
+          target = new_target;
+          if ((*_active)[target]) {
+            deactivate(target);
+          }
+
+          _highest = _sets.back().begin();
+          while (_highest != _sets.back().end() && 
+                 !(*_active)[_first[*_highest]]) {
+            ++_highest;
+          }
+        }
+      }
+    }    
+
+    void findMinCutIn() {
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _excess->set(n, 0);
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow->set(e, 0);
+      }
+
+      int bucket_num = 1;
       
-      return true;
-    }
-
-    Node findActiveNode() {
-      while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
-        --_highest_active;
-      }
-      if( _highest_active > 0) {
-        Node n = _active_nodes[_highest_active].front();
-        _active_nodes[_highest_active].pop_front();
-        return n;
-      } else {
-        return INVALID;
-      }
-    }
-
-    Value cutValue(bool out) {
-      Value value = 0;
-      if (out) {
-        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
-          for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
-            if (!(*_wake)[_graph->source(e)]){
-              value += (*_capacity)[e];
-            }   
-          }
-        }
-      } else {
-        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
-          for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
-            if (!(*_wake)[_graph->target(e)]){
-              value += (*_capacity)[e];
-            }   
-          }
-        }
-      }
-      return value;
-    }
-
+      {
+        typename Graph::template NodeMap<bool> reached(_graph, false);
+        
+        reached.set(_source, true);
+
+        bool first_set = true;
+
+        for (NodeIt t(_graph); t != INVALID; ++t) {
+          if (reached[t]) continue;
+          _sets.push_front(std::list<int>());
+          _sets.front().push_front(bucket_num);
+          _dormant[bucket_num] = !first_set;
+
+          _bucket->set(t, bucket_num);
+          _first[bucket_num] = _last[bucket_num] = t;
+          _next->set(t, INVALID);
+          _prev->set(t, INVALID);
+
+          ++bucket_num;
+          
+          std::vector<Node> queue;
+          queue.push_back(t);
+          reached.set(t, true);
+          
+          while (!queue.empty()) {
+            _sets.front().push_front(bucket_num);
+            _dormant[bucket_num] = !first_set;
+            _first[bucket_num] = _last[bucket_num] = INVALID;
+            
+            std::vector<Node> nqueue;
+            for (int i = 0; i < int(queue.size()); ++i) {
+              Node n = queue[i];
+              for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+                Node u = _graph.target(e);
+                if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+                  reached.set(u, true);
+                  addItem(u, bucket_num);
+                  nqueue.push_back(u);
+                }
+              }
+            }
+            queue.swap(nqueue);
+            ++bucket_num;
+          }
+          _sets.front().pop_front();
+          --bucket_num;
+          first_set = false;
+        }
+
+        _bucket->set(_source, 0);
+        _dormant[0] = true;
+      }
+      _source_set->set(_source, true);    
+          
+      Node target = _last[_sets.back().back()];
+      {
+        for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
+          if (_tolerance.positive((*_capacity)[e])) {
+            Node u = _graph.source(e);
+            _flow->set(e, (*_capacity)[e]);
+            _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+            if (!(*_active)[u] && u != _source) {
+              activate(u);
+            }
+          }
+        }
+        if ((*_active)[target]) {
+          deactivate(target);
+        }
+        
+        _highest = _sets.back().begin();
+        while (_highest != _sets.back().end() && 
+               !(*_active)[_first[*_highest]]) {
+          ++_highest;
+        }
+      }
+
+
+      while (true) {
+        while (_highest != _sets.back().end()) {
+          Node n = _first[*_highest];
+          Value excess = (*_excess)[n];
+          int next_bucket = _node_num;
+
+          int under_bucket;
+          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
+            under_bucket = -1;
+          } else {
+            under_bucket = *(++std::list<int>::iterator(_highest));
+          }
+
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node v = _graph.source(e);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] + excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, (*_capacity)[e]);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+
+          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node v = _graph.target(e);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] - excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, 0);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+          
+        no_more_push:
+          
+          _excess->set(n, excess);
+          
+          if (excess != 0) {
+            if ((*_next)[n] == INVALID) {
+              typename std::list<std::list<int> >::iterator new_set = 
+                _sets.insert(--_sets.end(), std::list<int>());
+              new_set->splice(new_set->end(), _sets.back(), 
+                              _sets.back().begin(), ++_highest);
+              for (std::list<int>::iterator it = new_set->begin();
+                   it != new_set->end(); ++it) {
+                _dormant[*it] = true;
+              }
+              while (_highest != _sets.back().end() && 
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else if (next_bucket == _node_num) {
+              _first[(*_bucket)[n]] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+              
+              std::list<std::list<int> >::iterator new_set = 
+                _sets.insert(--_sets.end(), std::list<int>());
+              
+              new_set->push_front(bucket_num);
+              _bucket->set(n, bucket_num);
+              _first[bucket_num] = _last[bucket_num] = n;
+              _next->set(n, INVALID);
+              _prev->set(n, INVALID);
+              _dormant[bucket_num] = true;
+              ++bucket_num;
+
+              while (_highest != _sets.back().end() && 
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else {
+              _first[*_highest] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+
+              while (next_bucket != *_highest) {
+                --_highest;
+              }
+              if (_highest == _sets.back().begin()) {
+                _sets.back().push_front(bucket_num);
+                _dormant[bucket_num] = false;
+                _first[bucket_num] = _last[bucket_num] = INVALID;
+                ++bucket_num;
+              }
+              --_highest;
+
+              _bucket->set(n, *_highest);
+              _next->set(n, _first[*_highest]);
+              if (_first[*_highest] != INVALID) {
+                _prev->set(_first[*_highest], n);
+              } else {
+                _last[*_highest] = n;
+              }
+              _first[*_highest] = n;          
+            }
+          } else {
+
+            deactivate(n);
+            if (!(*_active)[_first[*_highest]]) {
+              ++_highest;
+              if (_highest != _sets.back().end() && 
+                  !(*_active)[_first[*_highest]]) {
+                _highest = _sets.back().end();
+              }
+            }
+          }
+        }
+
+        if ((*_excess)[target] < _min_cut) {
+          _min_cut = (*_excess)[target];
+          for (NodeIt i(_graph); i != INVALID; ++i) {
+            _min_cut_map->set(i, false);
+          }
+          for (std::list<int>::iterator it = _sets.back().begin();
+               it != _sets.back().end(); ++it) {
+            Node n = _first[*it];
+            while (n != INVALID) {
+              _min_cut_map->set(n, true);
+              n = (*_next)[n];
+            }
+          }
+        }
+
+        {
+          Node new_target;
+          if ((*_prev)[target] != INVALID) {
+            _last[(*_bucket)[target]] = (*_prev)[target];
+            _next->set((*_prev)[target], INVALID);
+            new_target = (*_prev)[target];
+          } else {
+            _sets.back().pop_back();
+            if (_sets.back().empty()) {
+              _sets.pop_back();
+              if (_sets.empty())
+                break;
+              for (std::list<int>::iterator it = _sets.back().begin();
+                   it != _sets.back().end(); ++it) {
+                _dormant[*it] = false;
+              }
+            }
+            new_target = _last[_sets.back().back()];
+          }
+
+          _bucket->set(target, 0);
+
+          _source_set->set(target, true);         
+          for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.source(e);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(e, (*_capacity)[e]);
+          }
+          
+          for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
+            Value rem = (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.target(e);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(e, 0);
+          }
+          
+          target = new_target;
+          if ((*_active)[target]) {
+            deactivate(target);
+          }
+
+          _highest = _sets.back().begin();
+          while (_highest != _sets.back().end() && 
+                 !(*_active)[_first[*_highest]]) {
+            ++_highest;
+          }
+        }
+      }
+    }
 
   public:
@@ -436 +837 @@
     /// for the algorithm. 
     void init() {
-      init(NodeIt(*_graph));
+      init(NodeIt(_graph));
     }
 
@@ -447 +848 @@
     void init(const Node& source) {
       _source = source;
-      _node_num = countNodes(*_graph);
+      
+      _node_num = countNodes(_graph);
+      
+      _first.resize(_node_num);
+      _last.resize(_node_num);
 
       _dormant.resize(_node_num);
-      _level_size.resize(_node_num, 0);
-      _active_nodes.resize(_node_num);
-
-      if (!_preflow) {
-        _preflow = new FlowMap(*_graph);
-      }
-      if (!_wake) {
-        _wake = new WakeMap(*_graph);
-      }
-      if (!_dist) {
-        _dist = new DistMap(*_graph);
+
+      if (!_flow) {
+        _flow = new FlowMap(_graph);
+      }
+      if (!_next) {
+        _next = new typename Graph::template NodeMap<Node>(_graph);
+      }
+      if (!_prev) {
+        _prev = new typename Graph::template NodeMap<Node>(_graph);
+      }
+      if (!_active) {
+        _active = new typename Graph::template NodeMap<bool>(_graph);
+      }
+      if (!_bucket) {
+        _bucket = new typename Graph::template NodeMap<int>(_graph);
       }
       if (!_excess) {
-        _excess = new ExcessMap(*_graph);
+        _excess = new ExcessMap(_graph);
       }
       if (!_source_set) {
-        _source_set = new SourceSetMap(*_graph);
-      }
-      if (!_out_res_graph) {
-        _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
-      }
-      if (!_rev_graph) {
-        _rev_graph = new RevGraph(*_graph);
-      }
-      if (!_in_res_graph) {
-        _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
+        _source_set = new SourceSetMap(_graph);
       }
       if (!_min_cut_map) {
-        _min_cut_map = new MinCutMap(*_graph);
+        _min_cut_map = new MinCutMap(_graph);
       }
 
@@ -488 +888 @@
     /// source-side.
     ///
+    /// Calculates a minimum cut with \f$ source \f$ on the
+    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
+    /// \in X \f$ and minimal out-degree).
+    void calculateOut() {
+      findMinCutOut();
+    }
+
     /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
-    ///  and minimal out-degree).
-    void calculateOut() {
-      for (NodeIt it(*_graph); it != INVALID; ++it) {
-        if (it != _source) {
-          calculateOut(it);
-          return;
-        }
-      }
-    }
-
-    /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// source-side.
+    /// target-side.
     ///
-    /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
-    ///  and minimal out-degree). The \c target is the initial target
-    /// for the push-relabel algorithm.
-    void calculateOut(const Node& target) {
-      findMinCut(target, true, *_out_res_graph);
-    }
-
-    /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// sink-side.
-    ///
-    /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with 
-    /// \f$ source \notin X \f$
-    /// and minimal out-degree).
+    /// Calculates a minimum cut with \f$ source \f$ on the
+    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
+    /// \in X \f$ and minimal out-degree).
     void calculateIn() {
-      for (NodeIt it(*_graph); it != INVALID; ++it) {
-        if (it != _source) {
-          calculateIn(it);
-          return;
-        }
-      }
-    }
-
-    /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// sink-side.
-    ///
-    /// \brief Calculates a minimum cut with \f$ source \f$ on the
-    /// sink-side (i.e. a set \f$ X\subsetneq V 
-    /// \f$ with \f$ source \notin X \f$ and minimal out-degree).  
-    /// The \c target is the initial
-    /// target for the push-relabel algorithm.
-    void calculateIn(const Node& target) {
-      findMinCut(target, false, *_in_res_graph);
-    }
+      findMinCutIn();
+    }
+
 
     /// \brief Runs the algorithm.
@@ -546 +913 @@
     void run() {
       init();
-      for (NodeIt it(*_graph); it != INVALID; ++it) {
-        if (it != _source) {
-          calculateOut(it);
-          calculateIn(it);
-          return;
-        }
-      }
+      calculateOut();
+      calculateIn();
     }
 
@@ -562 +924 @@
     void run(const Node& s) {
       init(s);
-      for (NodeIt it(*_graph); it != INVALID; ++it) {
-        if (it != _source) {
-          calculateOut(it);
-          calculateIn(it);
-          return;
-        }
-      }
-    }
-
-    /// \brief Runs the algorithm.
-    ///
-    /// Runs the algorithm. It just calls the \ref init() and then
-    /// \ref calculateOut() and \ref calculateIn().
-    void run(const Node& s, const Node& t) {
-      init(s); 
-      calculateOut(t);
-      calculateIn(t);
+      calculateOut();
+      calculateIn();
     }
 
@@ -609 +956 @@
     template <typename NodeMap>
     Value minCut(NodeMap& nodeMap) const {
-      for (NodeIt it(*_graph); it != INVALID; ++it) {
+      for (NodeIt it(_graph); it != INVALID; ++it) {
         nodeMap.set(it, (*_min_cut_map)[it]);
       }

lemon/nagamochi_ibaraki.h

@@ -848 +848 @@
   ///
   /// The complexity of the algorithm is \f$ O(ne\log(n)) \f$ but with
-  /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$. 
-  /// When capacity map is neutral then it uses BucketHeap which
-  /// results \f$ O(ne) \f$ time complexity.
+  /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$.
+  /// When unit capacity minimum cut is computed then it uses
+  /// BucketHeap which results \f$ O(ne) \f$ time complexity.
   ///
   /// \warning The value type of the capacity map should be able to hold
@@ -907 +907 @@
     ///@{
 
-    struct DefNeutralCapacityTraits : public Traits {
+    struct DefUnitCapacityTraits : public Traits {
       typedef ConstMap<typename Graph::UEdge, Const<int, 1> > CapacityMap;
       static CapacityMap *createCapacityMap(const Graph&) {
@@ -918 +918 @@
     /// \ref named-templ-param "Named parameter" for setting 
     /// the capacity type to constMap<UEdge, int, 1>()
-    struct DefNeutralCapacity
+    struct DefUnitCapacity
       : public NagamochiIbaraki<Graph, CapacityMap, 
-                                DefNeutralCapacityTraits> { 
+                                DefUnitCapacityTraits> { 
       typedef NagamochiIbaraki<Graph, CapacityMap, 
-                               DefNeutralCapacityTraits> Create;
+                               DefUnitCapacityTraits> Create;
     };
 
@@ -1135 +1135 @@
     /// This constructor can be used only when the Traits class
     /// defines how can we instantiate a local capacity map.
-    /// If the DefNeutralCapacity used the algorithm automatically
+    /// If the DefUnitCapacity used the algorithm automatically
     /// construct the capacity map.
     ///