# HG changeset patch
# User Peter Kovacs <kpeter@inf.elte.hu>
# Date 2009-03-24 00:18:25
# Node ID 8c3112a66878ea474302c98be9cb538208f4866f
# Parent  425cc8328c0e78e2918db0dc04b4ef5ed22976fd

Use XTI implementation instead of ATI in NetworkSimplex (#234)

XTI (eXtended Threaded Index) is an imporved version of the widely
known ATI (Augmented Threaded Index) method for storing and updating
the spanning tree structure in Network Simplex algorithms.

In the ATI data structure three indices are stored for each node:
predecessor, thread and depth. In the XTI data structure depth is
replaced by the number of successors and the last successor
(according to the thread index).

diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
--- a/lemon/network_simplex.h
+++ b/lemon/network_simplex.h
@@ -154,10 +154,13 @@
     CostVector _pi;
 
     // Data for storing the spanning tree structure
-    IntVector _depth;
     IntVector _parent;
     IntVector _pred;
     IntVector _thread;
+    IntVector _rev_thread;
+    IntVector _succ_num;
+    IntVector _last_succ;
+    IntVector _dirty_revs;
     BoolVector _forward;
     IntVector _state;
     int _root;
@@ -850,11 +853,13 @@
       _flow.resize(all_arc_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);
+      _rev_thread.resize(all_node_num);
+      _succ_num.resize(all_node_num);
+      _last_succ.resize(all_node_num);
       _state.resize(all_arc_num, STATE_LOWER);
 
       // Initialize node related data
@@ -881,10 +886,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;
 
@@ -922,7 +929,9 @@
       Capacity max_cap = std::numeric_limits<Capacity>::max();
       for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
         _thread[u] = u + 1;
-        _depth[u] = 1;
+        _rev_thread[u + 1] = u;
+        _succ_num[u] = 1;
+        _last_succ[u] = u;
         _parent[u] = _root;
         _pred[u] = e;
         if (_supply[u] >= 0) {
@@ -946,11 +955,12 @@
     void findJoinNode() {
       int u = _source[in_arc];
       int v = _target[in_arc];
-      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;
     }
@@ -1026,88 +1036,147 @@
       }
     }
 
-    // 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];
       }
 
-      // 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
+      // 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) {
-        _thread[u] = new_stem = _parent[stem];
+        // Insert the next stem node into the thread list
+        new_stem = _parent[stem];
+        _thread[u] = new_stem;
+        _dirty_revs.push_back(u);
 
-        // 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;
+        // Remove the subtree of stem from the thread list
+        w = _rev_thread[stem];
+        _thread[w] = right;
+        _rev_thread[right] = w;
 
-        // Change the parent node of stem and shift stem and par_stem nodes
+        // Change the parent node and shift 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]) ;
+        // Update u and right
+        u = _last_succ[stem] == _last_succ[par_stem] ?
+          _rev_thread[par_stem] : _last_succ[stem];
         right = _thread[u];
       }
       _parent[u_out] = par_stem;
       _thread[u] = last;
+      _rev_thread[last] = u;
+      _last_succ[u_out] = u;
 
-      if (join == v_out && par_first) {
-        if (first != v_in) _thread[first] = right;
+      // 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 _pred, _forward, _last_succ and _succ_num for the
+      // stem nodes from u_out to u_in
+      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
+      u = u_out;
+      while (u != u_in) {
+        w = _parent[u];
+        _pred[u] = _pred[w];
+        _forward[u] = !_forward[w];
+        tmp_sc += _succ_num[u] - _succ_num[w];
+        _succ_num[u] = tmp_sc;
+        _last_succ[w] = tmp_ls;
+        u = w;
+      }
+      _pred[u_in] = in_arc;
+      _forward[u_in] = (u_in == _source[in_arc]);
+      _succ_num[u_in] = old_succ_num;
+
+      // 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 {
-        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
-        _thread[u] = right;
+        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];
+        }
+      }
+
+      // 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 _pred and _forward vectors
-    void updatePredArc() {
-      int u = u_out, v;
-      while (u != u_in) {
-        v = _parent[u];
-        _pred[u] = _pred[v];
-        _forward[u] = !_forward[v];
-        u = v;
-      }
-      _pred[u_in] = in_arc;
-      _forward[u_in] = (u_in == _source[in_arc]);
-    }
-
-    // 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;
-        _pi[u] += sigma;
+      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;
+        }
       }
     }
 
@@ -1139,9 +1208,8 @@
         bool change = findLeavingArc();
         changeFlow(change);
         if (change) {
-          updateThreadParent();
-          updatePredArc();
-          updateDepthPotential();
+          updateTreeStructure();
+          updatePotential();
         }
       }