[Lemon-commits] Peter Kovacs: Support real types + numerical sta...

[Lemon HG]

details:   http://lemon.cs.elte.hu/hg/lemon/rev/6ac5d9ae1d3d
changeset: 640:6ac5d9ae1d3d
user:      Peter Kovacs <kpeter [at] inf.elte.hu>
date:      Fri Apr 03 18:59:15 2009 +0200
description:
	Support real types + numerical stability fix in NS (#254)

	 - Real types are supported by appropriate inicialization.
	 - A feature of the XTI spanning tree structure is removed to ensure
	numerical stability (could cause problems using integer types).
	The node potentials are updated always on the lower subtree, in
	order to prevent overflow problems. The former method isn't
	notably faster during to our tests.

diffstat:

 lemon/network_simplex.h |  51 +++++++++++++++++++++++++++------------------------
 1 files changed, 27 insertions(+), 24 deletions(-)

diffs (103 lines):

diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
--- a/lemon/network_simplex.h
+++ b/lemon/network_simplex.h
@@ -54,7 +54,8 @@
   /// \tparam C The value type used for costs and potentials in the
   /// algorithm. By default it is the same as \c F.
   ///
-  /// \warning Both value types must be signed integer types.
+  /// \warning Both value types must be signed and all input data must
+  /// be integer.
   ///
   /// \note %NetworkSimplex provides five different pivot rule
   /// implementations. For more information see \ref PivotRule.
@@ -1044,8 +1045,10 @@
       }
 
       // Initialize arc maps
-      Flow max_cap = std::numeric_limits<Flow>::max();
-      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
+      Flow inf_cap =
+        std::numeric_limits<Flow>::has_infinity ?
+        std::numeric_limits<Flow>::infinity() :
+        std::numeric_limits<Flow>::max();
       if (_pupper && _pcost) {
         for (int i = 0; i != _arc_num; ++i) {
           Arc e = _arc_ref[i];
@@ -1069,7 +1072,7 @@
             _cap[i] = (*_pupper)[_arc_ref[i]];
         } else {
           for (int i = 0; i != _arc_num; ++i)
-            _cap[i] = max_cap;
+            _cap[i] = inf_cap;
         }
         if (_pcost) {
           for (int i = 0; i != _arc_num; ++i)
@@ -1079,6 +1082,18 @@
             _cost[i] = 1;
         }
       }
+      
+      // Initialize artifical cost
+      Cost art_cost;
+      if (std::numeric_limits<Cost>::is_exact) {
+        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
+      } else {
+        art_cost = std::numeric_limits<Cost>::min();
+        for (int i = 0; i != _arc_num; ++i) {
+          if (_cost[i] > art_cost) art_cost = _cost[i];
+        }
+        art_cost = (art_cost + 1) * _node_num;
+      }
 
       // Remove non-zero lower bounds
       if (_plower) {
@@ -1100,17 +1115,17 @@
         _last_succ[u] = u;
         _parent[u] = _root;
         _pred[u] = e;
-        _cost[e] = max_cost;
-        _cap[e] = max_cap;
+        _cost[e] = art_cost;
+        _cap[e] = inf_cap;
         _state[e] = STATE_TREE;
         if (_supply[u] >= 0) {
           _flow[e] = _supply[u];
           _forward[u] = true;
-          _pi[u] = -max_cost;
+          _pi[u] = -art_cost;
         } else {
           _flow[e] = -_supply[u];
           _forward[u] = false;
-          _pi[u] = max_cost;
+          _pi[u] = art_cost;
         }
       }
 
@@ -1327,22 +1342,10 @@
       Cost sigma = _forward[u_in] ?
         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
-      if (_succ_num[u_in] > _node_num / 2) {
-        // Update in the upper subtree (which contains the root)
-        int before = _rev_thread[u_in];
-        int after = _thread[_last_succ[u_in]];
-        _thread[before] = after;
-        _pi[_root] -= sigma;
-        for (int u = _thread[_root]; u != _root; u = _thread[u]) {
-          _pi[u] -= sigma;
-        }
-        _thread[before] = u_in;
-      } else {
-        // Update in the lower subtree (which has been moved)
-        int end = _thread[_last_succ[u_in]];
-        for (int u = u_in; u != end; u = _thread[u]) {
-          _pi[u] += sigma;
-        }
+      // Update potentials in the subtree, which has been moved
+      int end = _thread[_last_succ[u_in]];
+      for (int u = u_in; u != end; u = _thread[u]) {
+        _pi[u] += sigma;
       }
     }