# HG changeset patch
# User deba
# Date 1168550559 0
# Node ID 03c71d75499073498e82274596a199692c0fb08d
# Parent  c329fe995b40838fe7a7182f3b60c7a0adf47df4
Make Hao-Orlin epsilon-safe

diff -r c329fe995b40 -r 03c71d754990 lemon/graph_adaptor.h
--- a/lemon/graph_adaptor.h	Thu Jan 11 21:20:57 2007 +0000
+++ b/lemon/graph_adaptor.h	Thu Jan 11 21:22:39 2007 +0000
@@ -1450,7 +1450,7 @@
     void setFlow(const FlowMap& _flow) { flow = &_flow; }
 
     bool operator[](const typename Graph::Edge& e) const {
-      return tolerance.less((*flow)[e], (*capacity)[e]);
+      return tolerance.positive((*capacity)[e] - (*flow)[e]);
     }
   };
 
@@ -1473,7 +1473,7 @@
     void setCapacity(const CapacityMap& _capacity) { capacity = &_capacity; }
     void setFlow(const FlowMap& _flow) { flow = &_flow; }
     bool operator[](const typename Graph::Edge& e) const {
-      return tolerance.less(0, Number((*flow)[e]));
+      return tolerance.positive((*flow)[e]);
     }
   };
 
diff -r c329fe995b40 -r 03c71d754990 lemon/hao_orlin.h
--- a/lemon/hao_orlin.h	Thu Jan 11 21:20:57 2007 +0000
+++ b/lemon/hao_orlin.h	Thu Jan 11 21:22:39 2007 +0000
@@ -19,15 +19,19 @@
 #ifndef LEMON_HAO_ORLIN_H
 #define LEMON_HAO_ORLIN_H
 
+#include <cassert>
+ 
+
+
 #include <vector>
 #include <queue>
+#include <list>
 #include <limits>
 
 #include <lemon/maps.h>
 #include <lemon/graph_utils.h>
 #include <lemon/graph_adaptor.h>
 #include <lemon/iterable_maps.h>
- 
 
 /// \file
 /// \ingroup flowalgs
@@ -108,9 +112,6 @@
     
     OutResGraph* _out_res_graph;
 
-    typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap;
-    OutCurrentEdgeMap* _out_current_edge;  
-
     typedef RevGraphAdaptor<const Graph> RevGraph;
     RevGraph* _rev_graph;
 
@@ -120,10 +121,6 @@
     
     InResGraph* _in_res_graph;
 
-    typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap;
-    InCurrentEdgeMap* _in_current_edge;  
-
-
     typedef IterableBoolMap<Graph, Node> WakeMap;
     WakeMap* _wake;
 
@@ -161,8 +158,7 @@
              const Tolerance& tolerance = Tolerance()) :
       _graph(&graph), _capacity(&capacity), 
       _preflow(0), _source(), _target(), 
-      _out_res_graph(0), _out_current_edge(0),
-      _rev_graph(0), _in_res_graph(0), _in_current_edge(0),
+      _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) {}
@@ -171,18 +167,12 @@
       if (_min_cut_map) {
         delete _min_cut_map;
       } 
-      if (_in_current_edge) {
-        delete _in_current_edge;
-      }
       if (_in_res_graph) {
         delete _in_res_graph;
       }
       if (_rev_graph) {
         delete _rev_graph;
       }
-      if (_out_current_edge) {
-        delete _out_current_edge;
-      }
       if (_out_res_graph) {
         delete _out_res_graph;
       }
@@ -205,9 +195,9 @@
     
   private:
     
-    template <typename ResGraph, typename EdgeMap>
-    void findMinCut(const Node& target, bool out, 
-                    ResGraph& res_graph, EdgeMap& current_edge) {
+    template <typename ResGraph>
+    void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
+      typedef typename Graph::Node Node;
       typedef typename ResGraph::Edge ResEdge;
       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
 
@@ -219,28 +209,6 @@
         (*_dist)[it] = 1;
         (*_excess)[it] = 0;
         (*_source_set)[it] = false;
-
-        res_graph.firstOut(current_edge[it], it);
-      }
-
-      _target = target;
-      (*_dist)[target] = 0;
-
-      for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
-        Value delta = res_graph.rescap(it);
-        if (!_tolerance.positive(delta)) continue;
-        
-        (*_excess)[res_graph.source(it)] -= delta;
-        res_graph.augment(it, delta);
-        Node a = res_graph.target(it);
-        if (!_tolerance.positive((*_excess)[a]) && 
-            (*_wake)[a] && a != _target) {
-          _active_nodes[(*_dist)[a]].push_front(a);
-          if (_highest_active < (*_dist)[a]) {
-            _highest_active = (*_dist)[a];
-          }
-        }
-        (*_excess)[a] += delta;
       }
 
       _dormant[0].push_front(_source);
@@ -251,31 +219,46 @@
       _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) {
-	  ResEdge e;
-	  while (_tolerance.positive((*_excess)[n]) && 
-                 (e = findAdmissibleEdge(n, res_graph, current_edge)) 
-                 != INVALID){
-	    Value delta;
-	    if ((*_excess)[n] < res_graph.rescap(e)) {
+	  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;
+	    } else {
 	      delta = (*_excess)[n];
-	    } else {
-	      delta = res_graph.rescap(e);
-	      res_graph.nextOut(current_edge[n]);
+              (*_excess)[n] = 0;
 	    }
-            if (!_tolerance.positive(delta)) continue;
 	    res_graph.augment(e, delta);
-	    (*_excess)[res_graph.source(e)] -= delta;
-	    Node a = res_graph.target(e);
-	    if (!_tolerance.positive((*_excess)[a]) && a != _target) {
+	    if ((*_excess)[a] == 0 && a != _target) {
 	      _active_nodes[(*_dist)[a]].push_front(a);
 	    }
 	    (*_excess)[a] += delta;
+            if ((*_excess)[n] == 0) break;
 	  }
-	  if (_tolerance.positive((*_excess)[n])) {
-	    relabel(n, res_graph, current_edge);
+	  if ((*_excess)[n] != 0) {
+	    relabel(n, res_graph);
           }
 	}
 
@@ -297,8 +280,8 @@
       } while (selectNewSink(res_graph));
     }
 
-    template <typename ResGraph, typename EdgeMap>
-    void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) {
+    template <typename ResGraph>
+    void relabel(const Node& n, ResGraph& res_graph) {
       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
 
       int k = (*_dist)[n];
@@ -311,7 +294,7 @@
 	    --_level_size[(*_dist)[it]];
 	  }
 	}
-	--_highest_active;
+        --_highest_active;
       } else {	
         int new_dist = _node_num;
         for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
@@ -331,7 +314,6 @@
 	  _highest_active = (*_dist)[n];
 	  _active_nodes[_highest_active].push_front(n);
 	  ++_level_size[(*_dist)[n]];
-	  res_graph.firstOut(current_edge[n], n);
 	}
       }
     }
@@ -372,38 +354,34 @@
 
       if (wake_was_empty){
 	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
-          if (_tolerance.positive((*_excess)[it])) {
-	    if ((*_wake)[it] && it != _target) {
-	      _active_nodes[(*_dist)[it]].push_front(it);
-            }
-	    if (_highest_active < (*_dist)[it]) {
-	      _highest_active = (*_dist)[it];		    
+          if ((*_excess)[it] != 0 && it != _target) {
+            _active_nodes[(*_dist)[it]].push_front(it);
+            if (_highest_active < (*_dist)[it]) {
+              _highest_active = (*_dist)[it];		    
             }
 	  }
 	}
       }
 
-      for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){
-	if (!(*_source_set)[res_graph.target(e)]) {
-          Value delta = res_graph.rescap(e);
-          if (!_tolerance.positive(delta)) continue;
-          res_graph.augment(e, delta);
-          (*_excess)[res_graph.source(e)] -= delta;
-          Node a = res_graph.target(e);
-          if (!_tolerance.positive((*_excess)[a]) && 
-              (*_wake)[a] && a != _target) {
-            _active_nodes[(*_dist)[a]].push_front(a);
-            if (_highest_active < (*_dist)[a]) {
-              _highest_active = (*_dist)[a];
-            }
+      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;
-	}
+        }
+        (*_excess)[a] += delta;
       }
       
       return true;
     }
-    
+
     Node findActiveNode() {
       while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
 	--_highest_active;
@@ -417,25 +395,6 @@
       }
     }
 
-    template <typename ResGraph, typename EdgeMap>
-    typename ResGraph::Edge findAdmissibleEdge(const Node& n, 
-                                               ResGraph& res_graph, 
-                                               EdgeMap& current_edge) {
-      typedef typename ResGraph::Edge ResEdge;
-      ResEdge e = current_edge[n];
-      while (e != INVALID && 
-             ((*_dist)[n] <= (*_dist)[res_graph.target(e)] || 
-              !(*_wake)[res_graph.target(e)])) {
-	res_graph.nextOut(e);
-      }
-      if (e != INVALID) {
-	current_edge[n] = e;	
-	return e;
-      } else {
-	return INVALID;
-      }
-    }
-
     Value cutValue(bool out) {
       Value value = 0;
       if (out) {
@@ -512,18 +471,12 @@
       if (!_out_res_graph) {
         _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
       }
-      if (!_out_current_edge) {
-        _out_current_edge = new OutCurrentEdgeMap(*_graph);
-      }
       if (!_rev_graph) {
         _rev_graph = new RevGraph(*_graph);
       }
       if (!_in_res_graph) {
         _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
       }
-      if (!_in_current_edge) {
-        _in_current_edge = new InCurrentEdgeMap(*_graph);
-      }
       if (!_min_cut_map) {
         _min_cut_map = new MinCutMap(*_graph);
       }
@@ -555,7 +508,7 @@
     ///  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, *_out_current_edge);
+      findMinCut(target, true, *_out_res_graph);
     }
 
     /// \brief Calculates a minimum cut with \f$ source \f$ on the
@@ -583,7 +536,7 @@
     /// The \c target is the initial
     /// target for the push-relabel algorithm.
     void calculateIn(const Node& target) {
-      findMinCut(target, false, *_in_res_graph, *_in_current_edge);
+      findMinCut(target, false, *_in_res_graph);
     }
 
     /// \brief Runs the algorithm.