# HG changeset patch
# User kpeter
# Date 1244078346 0
# Node ID 61bf01404ede7afad822fe0bca76a335cdc252a0
# Parent  bafe730864daa80c0ec2fa4f4923be964f64736a
Various improvements for NS pivot rules

diff -r bafe730864da -r 61bf01404ede lemon/network_simplex.h
--- a/lemon/network_simplex.h	Mon Jun 01 16:53:59 2009 +0000
+++ b/lemon/network_simplex.h	Thu Jun 04 01:19:06 2009 +0000
@@ -241,10 +241,10 @@
       BlockSearchPivotRule(NetworkSimplex &ns) :
         _edge(ns._edge), _source(ns._source), _target(ns._target),
         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-        _in_edge(ns._in_edge), _edge_num(ns._edge_num + ns._node_num), _next_edge(0)
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
       {
         // The main parameters of the pivot rule
-        const double BLOCK_SIZE_FACTOR = 2.0;
+        const double BLOCK_SIZE_FACTOR = 0.5;
         const int MIN_BLOCK_SIZE = 10;
 
         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
@@ -255,33 +255,32 @@
       bool findEnteringEdge() {
         Cost c, min = 0;
         int cnt = _block_size;
-        int e, min_edge = _next_edge;
+        int e;
         for (e = _next_edge; e < _edge_num; ++e) {
           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (c < min) {
             min = c;
-            min_edge = e;
+            _in_edge = e;
           }
           if (--cnt == 0) {
-            if (min < 0) break;
+            if (min < 0) goto search_end;
             cnt = _block_size;
           }
         }
-        if (min == 0 || cnt > 0) {
-          for (e = 0; e < _next_edge; ++e) {
-            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-            if (c < min) {
-              min = c;
-              min_edge = e;
-            }
-            if (--cnt == 0) {
-              if (min < 0) break;
-              cnt = _block_size;
-            }
+        for (e = 0; e < _next_edge; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < min) {
+            min = c;
+            _in_edge = e;
+          }
+          if (--cnt == 0) {
+            if (min < 0) goto search_end;
+            cnt = _block_size;
           }
         }
         if (min >= 0) return false;
-        _in_edge = min_edge;
+
+      search_end:
         _next_edge = e;
         return true;
       }
@@ -325,7 +324,7 @@
         _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
       {
         // The main parameters of the pivot rule
-        const double LIST_LENGTH_FACTOR = 1.0;
+        const double LIST_LENGTH_FACTOR = 0.25;
         const int MIN_LIST_LENGTH = 10;
         const double MINOR_LIMIT_FACTOR = 0.1;
         const int MIN_MINOR_LIMIT = 3;
@@ -341,7 +340,7 @@
       /// Find next entering edge
       bool findEnteringEdge() {
         Cost min, c;
-        int e, min_edge = _next_edge;
+        int e;
         if (_curr_length > 0 && _minor_count < _minor_limit) {
           // Minor iteration: select the best eligible edge from the
           // current candidate list
@@ -352,16 +351,13 @@
             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
             if (c < min) {
               min = c;
-              min_edge = e;
+              _in_edge = e;
             }
-            if (c >= 0) {
+            else if (c >= 0) {
               _candidates[i--] = _candidates[--_curr_length];
             }
           }
-          if (min < 0) {
-            _in_edge = min_edge;
-            return true;
-          }
+          if (min < 0) return true;
         }
 
         // Major iteration: build a new candidate list
@@ -373,27 +369,26 @@
             _candidates[_curr_length++] = e;
             if (c < min) {
               min = c;
-              min_edge = e;
+              _in_edge = e;
             }
-            if (_curr_length == _list_length) break;
+            if (_curr_length == _list_length) goto search_end;
           }
         }
-        if (_curr_length < _list_length) {
-          for (e = 0; e < _next_edge; ++e) {
-            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-            if (c < 0) {
-              _candidates[_curr_length++] = e;
-              if (c < min) {
-                min = c;
-                min_edge = e;
-              }
-              if (_curr_length == _list_length) break;
+        for (e = 0; e < _next_edge; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _candidates[_curr_length++] = e;
+            if (c < min) {
+              min = c;
+              _in_edge = e;
             }
+            if (_curr_length == _list_length) goto search_end;
           }
         }
         if (_curr_length == 0) return false;
+      
+      search_end:        
         _minor_count = 1;
-        _in_edge = min_edge;
         _next_edge = e;
         return true;
       }
@@ -451,7 +446,7 @@
          _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
       {
         // The main parameters of the pivot rule
-        const double BLOCK_SIZE_FACTOR = 1.5;
+        const double BLOCK_SIZE_FACTOR = 1.0;
         const int MIN_BLOCK_SIZE = 10;
         const double HEAD_LENGTH_FACTOR = 0.1;
         const int MIN_HEAD_LENGTH = 3;
@@ -479,39 +474,35 @@
 
         // Extend the list
         int cnt = _block_size;
-        int last_edge = 0;
         int limit = _head_length;
         
-        for (int e = _next_edge; e < _edge_num; ++e) {
+        for (e = _next_edge; e < _edge_num; ++e) {
           _cand_cost[e] = _state[e] *
             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (_cand_cost[e] < 0) {
             _candidates[_curr_length++] = e;
-            last_edge = e;
           }
           if (--cnt == 0) {
-            if (_curr_length > limit) break;
+            if (_curr_length > limit) goto search_end;
             limit = 0;
             cnt = _block_size;
           }
         }
-        if (_curr_length <= limit) {
-          for (int e = 0; e < _next_edge; ++e) {
-            _cand_cost[e] = _state[e] *
-              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-            if (_cand_cost[e] < 0) {
-              _candidates[_curr_length++] = e;
-              last_edge = e;
-            }
-            if (--cnt == 0) {
-              if (_curr_length > limit) break;
-              limit = 0;
-              cnt = _block_size;
-            }
+        for (e = 0; e < _next_edge; ++e) {
+          _cand_cost[e] = _state[e] *
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (_cand_cost[e] < 0) {
+            _candidates[_curr_length++] = e;
+          }
+          if (--cnt == 0) {
+            if (_curr_length > limit) goto search_end;
+            limit = 0;
+            cnt = _block_size;
           }
         }
         if (_curr_length == 0) return false;
-        _next_edge = last_edge + 1;
+        
+      search_end:
 
         // Make heap of the candidate list (approximating a partial sort)
         make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
@@ -519,6 +510,7 @@
 
         // Pop the first element of the heap
         _in_edge = _candidates[0];
+        _next_edge = e;
         pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
                   _sort_func );
         _curr_length = std::min(_head_length, _curr_length - 1);
@@ -931,12 +923,10 @@
       _supply[_root] = 0;
       _pi[_root] = 0;
       
-      // Store the edges in a mixed order
-      int k = std::max(int(sqrt(_edge_num)), 10);
+      // Store the edges
       int i = 0;
       for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
-        _edge[i] = e;
-        if ((i += k) >= _edge_num) i = (i % k) + 1;
+        _edge[i++] = e;
       }
 
       // Initialize edge maps
@@ -961,7 +951,7 @@
       }
 
       // Add artificial edges and initialize the spanning tree data structure
-      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
+      Cost max_cost = std::numeric_limits<Cost>::max() / 2 + 1;
       Capacity max_cap = std::numeric_limits<Capacity>::max();
       for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
         _parent[u] = _root;