# HG changeset patch
# User kpeter
# Date 1243870527 0
# Node ID 23570e368afab1a47d4c98ccf9a5e34351511b82
# Parent  e98bbe64cca4f68e74f2da559d49d45bceee69c4
XTI data structure for NS - backport hg commit [8c3112a66878]

diff -r e98bbe64cca4 -r 23570e368afa lemon/network_simplex.h
--- a/lemon/network_simplex.h	Fri Feb 06 21:52:34 2009 +0000
+++ b/lemon/network_simplex.h	Mon Jun 01 15:35:27 2009 +0000
@@ -241,7 +241,7 @@
       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), _next_edge(0)
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num + ns._node_num), _next_edge(0)
       {
         // The main parameters of the pivot rule
         const double BLOCK_SIZE_FACTOR = 2.0;
@@ -591,14 +591,19 @@
     // The number of edges in the original graph
     int _edge_num;
 
-    // The depth vector of the spanning tree structure
-    IntVector _depth;
     // The parent vector of the spanning tree structure
     IntVector _parent;
     // The pred_edge vector of the spanning tree structure
     IntVector _pred;
     // The thread vector of the spanning tree structure
     IntVector _thread;
+
+    IntVector _rev_thread;
+    IntVector _succ_num;
+    IntVector _last_succ;
+
+    IntVector _dirty_revs;
+    
     // The forward vector of the spanning tree structure
     BoolVector _forward;
     // The state vector
@@ -882,11 +887,13 @@
       _flow.resize(all_edge_num, 0);
       _pi.resize(all_node_num, 0);
       
-      _depth.resize(all_node_num);
       _parent.resize(all_node_num);
       _pred.resize(all_node_num);
+      _forward.resize(all_node_num);
       _thread.resize(all_node_num);
-      _forward.resize(all_node_num);
+      _rev_thread.resize(all_node_num);
+      _succ_num.resize(all_node_num);
+      _last_succ.resize(all_node_num);
       _state.resize(all_edge_num, STATE_LOWER);
       
       // Initialize node related data
@@ -915,10 +922,12 @@
       
       // Set data for the artificial root node
       _root = _node_num;
-      _depth[_root] = 0;
       _parent[_root] = -1;
       _pred[_root] = -1;
       _thread[_root] = 0;
+      _rev_thread[0] = _root;
+      _succ_num[_root] = all_node_num;
+      _last_succ[_root] = _root - 1;
       _supply[_root] = 0;
       _pi[_root] = 0;
       
@@ -955,22 +964,30 @@
       Cost max_cost = std::numeric_limits<Cost>::max() / 4;
       Capacity max_cap = std::numeric_limits<Capacity>::max();
       for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
-        _thread[u] = u + 1;
-        _depth[u] = 1;
         _parent[u] = _root;
         _pred[u] = e;
+        _thread[u] = u + 1;
+        _rev_thread[u + 1] = u;
+        _succ_num[u] = 1;
+        _last_succ[u] = u;
+        _cap[e] = max_cap;
+        _state[e] = STATE_TREE;
         if (_supply[u] >= 0) {
+          _forward[u] = true;
+          _pi[u] = 0;
+          _source[e] = u;
+          _target[e] = _root;
           _flow[e] = _supply[u];
-          _forward[u] = true;
-          _pi[u] = -max_cost;
-        } else {
-          _flow[e] = -_supply[u];
+          _cost[e] = 0;
+        }
+        else {
           _forward[u] = false;
           _pi[u] = max_cost;
+          _source[e] = _root;
+          _target[e] = u;
+          _flow[e] = -_supply[u];
+          _cost[e] = max_cost;
         }
-        _cost[e] = max_cost;
-        _cap[e] = max_cap;
-        _state[e] = STATE_TREE;
       }
 
       return true;
@@ -980,11 +997,12 @@
     void findJoinNode() {
       int u = _source[_in_edge];
       int v = _target[_in_edge];
-      while (_depth[u] > _depth[v]) u = _parent[u];
-      while (_depth[v] > _depth[u]) v = _parent[v];
       while (u != v) {
-        u = _parent[u];
-        v = _parent[v];
+        if (_succ_num[u] < _succ_num[v]) {
+          u = _parent[u];
+        } else {
+          v = _parent[v];
+        }
       }
       join = u;
     }
@@ -1059,88 +1077,144 @@
         _state[_in_edge] = -_state[_in_edge];
       }
     }
-
-    // Update _thread and _parent vectors
-    void updateThreadParent() {
-      int u;
+    
+    // Update the tree structure
+    void updateTreeStructure() {
+      int u, w;
+      int old_rev_thread = _rev_thread[u_out];
+      int old_succ_num = _succ_num[u_out];
+      int old_last_succ = _last_succ[u_out];
       v_out = _parent[u_out];
 
-      // Handle the case when join and v_out coincide
-      bool par_first = false;
-      if (join == v_out) {
-        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
-        if (u == v_in) {
-          par_first = true;
-          while (_thread[u] != u_out) u = _thread[u];
-          first = u;
+      u = _last_succ[u_in];  // the last successor of u_in
+      right = _thread[u];    // the node after it
+      
+      // Handle the case when old_rev_thread equals to v_in 
+      // (it also means that join and v_out coincide)
+      if (old_rev_thread == v_in) {
+        last = _thread[_last_succ[u_out]];
+      } else {
+        last = _thread[v_in];
+      }
+
+      // Update _thread and _parent along the stem nodes (i.e. the nodes 
+      // between u_in and u_out, whose parent have to be changed)
+      _thread[v_in] = stem = u_in;
+      _dirty_revs.clear();
+      _dirty_revs.push_back(v_in);
+      par_stem = v_in;
+      while (stem != u_out) {
+        // Insert the next stem node into the thread list
+        new_stem = _parent[stem];
+        _thread[u] = new_stem;
+        _dirty_revs.push_back(u);
+        
+        // Remove the subtree of stem from the thread list
+        w = _rev_thread[stem];
+        _thread[w] = right;
+        _rev_thread[right] = w;
+
+        // Change the parent node and shift stem nodes
+        _parent[stem] = par_stem;        
+        par_stem = stem;
+        stem = new_stem;
+
+        // Update u and right nodes
+        u = _last_succ[stem] == _last_succ[par_stem] ?
+          _rev_thread[par_stem] : _last_succ[stem];
+        right = _thread[u];
+      }
+      _parent[u_out] = par_stem;
+      _last_succ[u_out] = u;
+      _thread[u] = last;
+      _rev_thread[last] = u;
+
+      // Remove the subtree of u_out from the thread list except for 
+      // the case when old_rev_thread equals to v_in 
+      // (it also means that join and v_out coincide)
+      if (old_rev_thread != v_in) {
+        _thread[old_rev_thread] = right;
+        _rev_thread[right] = old_rev_thread;
+      }
+      
+      // Update _rev_thread using the new _thread values
+      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
+        u = _dirty_revs[i];
+        _rev_thread[_thread[u]] = u;
+      }
+      
+      // Update _last_succ for the stem nodes from u_out to u_in
+      for (u = u_out; u != u_in; u = _parent[u]) {
+        _last_succ[_parent[u]] = _last_succ[u];      
+      }
+
+      // Set limits for updating _last_succ form v_in and v_out
+      // towards the root
+      int up_limit_in = -1;
+      int up_limit_out = -1;
+      if (_last_succ[join] == v_in) {
+        up_limit_out = join;
+      } else {
+        up_limit_in = join;
+      }
+
+      // Update _last_succ from v_in towards the root
+      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in; 
+           u = _parent[u]) {
+        _last_succ[u] = _last_succ[u_out];
+      }
+      // Update _last_succ from v_out towards the root
+      if (join != old_rev_thread && v_in != old_rev_thread) {
+        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; 
+             u = _parent[u]) {
+          _last_succ[u] = old_rev_thread;
+        }
+      } else {
+        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+             u = _parent[u]) {
+          _last_succ[u] = _last_succ[u_out];
         }
       }
 
-      // Find the last successor of u_in (u) and the node after it (right)
-      // according to the thread index
-      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
-      right = _thread[u];
-      if (_thread[v_in] == u_out) {
-        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
-        if (last == u_out) last = _thread[last];
-      }
-      else last = _thread[v_in];
-
-      // Update stem nodes
-      _thread[v_in] = stem = u_in;
-      par_stem = v_in;
-      while (stem != u_out) {
-        _thread[u] = new_stem = _parent[stem];
-
-        // Find the node just before the stem node (u) according to
-        // the original thread index
-        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
-        _thread[u] = right;
-
-        // Change the parent node of stem and shift stem and par_stem nodes
-        _parent[stem] = par_stem;
-        par_stem = stem;
-        stem = new_stem;
-
-        // Find the last successor of stem (u) and the node after it (right)
-        // according to the thread index
-        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
-        right = _thread[u];
-      }
-      _parent[u_out] = par_stem;
-      _thread[u] = last;
-
-      if (join == v_out && par_first) {
-        if (first != v_in) _thread[first] = right;
-      } else {
-        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
-        _thread[u] = right;
-      }
-    }
-
-    // Update _pred and _forward vectors
-    void updatePredEdge() {
-      int u = u_out, v;
+      // Update _pred and _forward for the stem nodes from u_out to u_in
+      u = u_out;
       while (u != u_in) {
-        v = _parent[u];
-        _pred[u] = _pred[v];
-        _forward[u] = !_forward[v];
-        u = v;
+        w = _parent[u];
+        _pred[u] = _pred[w];
+        _forward[u] = !_forward[w];
+        u = w;
       }
       _pred[u_in] = _in_edge;
       _forward[u_in] = (u_in == _source[_in_edge]);
+      
+      // Update _succ_num from v_in to join
+      for (u = v_in; u != join; u = _parent[u]) {
+        _succ_num[u] += old_succ_num;
+      }
+      // Update _succ_num from v_out to join
+      for (u = v_out; u != join; u = _parent[u]) {
+        _succ_num[u] -= old_succ_num;
+      }
+      // Update _succ_num for the stem nodes from u_out to u_in
+      int tmp = 0;
+      u = u_out;
+      while (u != u_in) {
+        w = _parent[u];
+        tmp = _succ_num[u] - _succ_num[w] + tmp;
+        _succ_num[u] = tmp;
+        u = w;
+      }
+      _succ_num[u_in] = old_succ_num;
     }
 
-    // Update _depth and _potential vectors
-    void updateDepthPotential() {
-      _depth[u_in] = _depth[v_in] + 1;
+    // Update potentials
+    void updatePotential() {
       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]];
-      _pi[u_in] += sigma;
-      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
-        _depth[u] = _depth[_parent[u]] + 1;
-        if (_depth[u] <= _depth[u_in]) break;
+      // 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;
       }
     }
@@ -1166,19 +1240,18 @@
     template<class PivotRuleImplementation>
     bool start() {
       PivotRuleImplementation pivot(*this);
-
+      
       // Execute the network simplex algorithm
       while (pivot.findEnteringEdge()) {
         findJoinNode();
         bool change = findLeavingEdge();
         changeFlow(change);
         if (change) {
-          updateThreadParent();
-          updatePredEdge();
-          updateDepthPotential();
+          updateTreeStructure();
+          updatePotential();
         }
       }
-
+      
       // Check if the flow amount equals zero on all the artificial edges
       for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
         if (_flow[e] > 0) return false;
@@ -1199,7 +1272,7 @@
       for (int i = 0; i != _node_num; ++i) {
         (*_potential_result)[_node[i]] = _pi[i];
       }
-
+      
       return true;
     }