Changeset 608:6ac5d9ae1d3d in lemon-main

Timestamp:

04/03/09 18:59:15 (16 years ago)

Author:

Peter Kovacs <kpeter@…>

Branch:

default

Phase:

public

Message:

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.

File:

• lemon/network_simplex.h (modified) (6 diffs)

lemon/network_simplex.h

@@ -55 +55 @@
   /// 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
@@ -1045 +1046 @@
 
       // 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) {
@@ -1070 +1073 @@
         } else {
           for (int i = 0; i != _arc_num; ++i)
-            _cap[i] = max_cap;
+            _cap[i] = inf_cap;
         }
         if (_pcost) {
@@ -1079 +1082 @@
             _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;
       }
 
@@ -1101 +1116 @@
         _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;
         }
       }
@@ -1328 +1343 @@
         _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;
       }
     }