Use XTI implementation instead of ATI in NetworkSimplex (#234)
authorPeter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Mar 2009 00:18:25 +0100
changeset 6048c3112a66878
parent 603 425cc8328c0e
child 605 5232721b3f14
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).
lemon/network_simplex.h
     1.1 --- a/lemon/network_simplex.h	Mon Mar 23 23:54:42 2009 +0100
     1.2 +++ b/lemon/network_simplex.h	Tue Mar 24 00:18:25 2009 +0100
     1.3 @@ -154,10 +154,13 @@
     1.4      CostVector _pi;
     1.5  
     1.6      // Data for storing the spanning tree structure
     1.7 -    IntVector _depth;
     1.8      IntVector _parent;
     1.9      IntVector _pred;
    1.10      IntVector _thread;
    1.11 +    IntVector _rev_thread;
    1.12 +    IntVector _succ_num;
    1.13 +    IntVector _last_succ;
    1.14 +    IntVector _dirty_revs;
    1.15      BoolVector _forward;
    1.16      IntVector _state;
    1.17      int _root;
    1.18 @@ -850,11 +853,13 @@
    1.19        _flow.resize(all_arc_num, 0);
    1.20        _pi.resize(all_node_num, 0);
    1.21  
    1.22 -      _depth.resize(all_node_num);
    1.23        _parent.resize(all_node_num);
    1.24        _pred.resize(all_node_num);
    1.25        _forward.resize(all_node_num);
    1.26        _thread.resize(all_node_num);
    1.27 +      _rev_thread.resize(all_node_num);
    1.28 +      _succ_num.resize(all_node_num);
    1.29 +      _last_succ.resize(all_node_num);
    1.30        _state.resize(all_arc_num, STATE_LOWER);
    1.31  
    1.32        // Initialize node related data
    1.33 @@ -881,10 +886,12 @@
    1.34  
    1.35        // Set data for the artificial root node
    1.36        _root = _node_num;
    1.37 -      _depth[_root] = 0;
    1.38        _parent[_root] = -1;
    1.39        _pred[_root] = -1;
    1.40        _thread[_root] = 0;
    1.41 +      _rev_thread[0] = _root;
    1.42 +      _succ_num[_root] = all_node_num;
    1.43 +      _last_succ[_root] = _root - 1;
    1.44        _supply[_root] = 0;
    1.45        _pi[_root] = 0;
    1.46  
    1.47 @@ -922,7 +929,9 @@
    1.48        Capacity max_cap = std::numeric_limits<Capacity>::max();
    1.49        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
    1.50          _thread[u] = u + 1;
    1.51 -        _depth[u] = 1;
    1.52 +        _rev_thread[u + 1] = u;
    1.53 +        _succ_num[u] = 1;
    1.54 +        _last_succ[u] = u;
    1.55          _parent[u] = _root;
    1.56          _pred[u] = e;
    1.57          if (_supply[u] >= 0) {
    1.58 @@ -946,11 +955,12 @@
    1.59      void findJoinNode() {
    1.60        int u = _source[in_arc];
    1.61        int v = _target[in_arc];
    1.62 -      while (_depth[u] > _depth[v]) u = _parent[u];
    1.63 -      while (_depth[v] > _depth[u]) v = _parent[v];
    1.64        while (u != v) {
    1.65 -        u = _parent[u];
    1.66 -        v = _parent[v];
    1.67 +        if (_succ_num[u] < _succ_num[v]) {
    1.68 +          u = _parent[u];
    1.69 +        } else {
    1.70 +          v = _parent[v];
    1.71 +        }
    1.72        }
    1.73        join = u;
    1.74      }
    1.75 @@ -1026,88 +1036,147 @@
    1.76        }
    1.77      }
    1.78  
    1.79 -    // Update _thread and _parent vectors
    1.80 -    void updateThreadParent() {
    1.81 -      int u;
    1.82 +    // Update the tree structure
    1.83 +    void updateTreeStructure() {
    1.84 +      int u, w;
    1.85 +      int old_rev_thread = _rev_thread[u_out];
    1.86 +      int old_succ_num = _succ_num[u_out];
    1.87 +      int old_last_succ = _last_succ[u_out];
    1.88        v_out = _parent[u_out];
    1.89  
    1.90 -      // Handle the case when join and v_out coincide
    1.91 -      bool par_first = false;
    1.92 -      if (join == v_out) {
    1.93 -        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
    1.94 -        if (u == v_in) {
    1.95 -          par_first = true;
    1.96 -          while (_thread[u] != u_out) u = _thread[u];
    1.97 -          first = u;
    1.98 -        }
    1.99 +      u = _last_succ[u_in];  // the last successor of u_in
   1.100 +      right = _thread[u];    // the node after it
   1.101 +
   1.102 +      // Handle the case when old_rev_thread equals to v_in
   1.103 +      // (it also means that join and v_out coincide)
   1.104 +      if (old_rev_thread == v_in) {
   1.105 +        last = _thread[_last_succ[u_out]];
   1.106 +      } else {
   1.107 +        last = _thread[v_in];
   1.108        }
   1.109  
   1.110 -      // Find the last successor of u_in (u) and the node after it (right)
   1.111 -      // according to the thread index
   1.112 -      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
   1.113 -      right = _thread[u];
   1.114 -      if (_thread[v_in] == u_out) {
   1.115 -        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
   1.116 -        if (last == u_out) last = _thread[last];
   1.117 -      }
   1.118 -      else last = _thread[v_in];
   1.119 -
   1.120 -      // Update stem nodes
   1.121 +      // Update _thread and _parent along the stem nodes (i.e. the nodes
   1.122 +      // between u_in and u_out, whose parent have to be changed)
   1.123        _thread[v_in] = stem = u_in;
   1.124 +      _dirty_revs.clear();
   1.125 +      _dirty_revs.push_back(v_in);
   1.126        par_stem = v_in;
   1.127        while (stem != u_out) {
   1.128 -        _thread[u] = new_stem = _parent[stem];
   1.129 +        // Insert the next stem node into the thread list
   1.130 +        new_stem = _parent[stem];
   1.131 +        _thread[u] = new_stem;
   1.132 +        _dirty_revs.push_back(u);
   1.133  
   1.134 -        // Find the node just before the stem node (u) according to
   1.135 -        // the original thread index
   1.136 -        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
   1.137 -        _thread[u] = right;
   1.138 +        // Remove the subtree of stem from the thread list
   1.139 +        w = _rev_thread[stem];
   1.140 +        _thread[w] = right;
   1.141 +        _rev_thread[right] = w;
   1.142  
   1.143 -        // Change the parent node of stem and shift stem and par_stem nodes
   1.144 +        // Change the parent node and shift stem nodes
   1.145          _parent[stem] = par_stem;
   1.146          par_stem = stem;
   1.147          stem = new_stem;
   1.148  
   1.149 -        // Find the last successor of stem (u) and the node after it (right)
   1.150 -        // according to the thread index
   1.151 -        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
   1.152 +        // Update u and right
   1.153 +        u = _last_succ[stem] == _last_succ[par_stem] ?
   1.154 +          _rev_thread[par_stem] : _last_succ[stem];
   1.155          right = _thread[u];
   1.156        }
   1.157        _parent[u_out] = par_stem;
   1.158        _thread[u] = last;
   1.159 +      _rev_thread[last] = u;
   1.160 +      _last_succ[u_out] = u;
   1.161  
   1.162 -      if (join == v_out && par_first) {
   1.163 -        if (first != v_in) _thread[first] = right;
   1.164 +      // Remove the subtree of u_out from the thread list except for
   1.165 +      // the case when old_rev_thread equals to v_in
   1.166 +      // (it also means that join and v_out coincide)
   1.167 +      if (old_rev_thread != v_in) {
   1.168 +        _thread[old_rev_thread] = right;
   1.169 +        _rev_thread[right] = old_rev_thread;
   1.170 +      }
   1.171 +
   1.172 +      // Update _rev_thread using the new _thread values
   1.173 +      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
   1.174 +        u = _dirty_revs[i];
   1.175 +        _rev_thread[_thread[u]] = u;
   1.176 +      }
   1.177 +
   1.178 +      // Update _pred, _forward, _last_succ and _succ_num for the
   1.179 +      // stem nodes from u_out to u_in
   1.180 +      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
   1.181 +      u = u_out;
   1.182 +      while (u != u_in) {
   1.183 +        w = _parent[u];
   1.184 +        _pred[u] = _pred[w];
   1.185 +        _forward[u] = !_forward[w];
   1.186 +        tmp_sc += _succ_num[u] - _succ_num[w];
   1.187 +        _succ_num[u] = tmp_sc;
   1.188 +        _last_succ[w] = tmp_ls;
   1.189 +        u = w;
   1.190 +      }
   1.191 +      _pred[u_in] = in_arc;
   1.192 +      _forward[u_in] = (u_in == _source[in_arc]);
   1.193 +      _succ_num[u_in] = old_succ_num;
   1.194 +
   1.195 +      // Set limits for updating _last_succ form v_in and v_out
   1.196 +      // towards the root
   1.197 +      int up_limit_in = -1;
   1.198 +      int up_limit_out = -1;
   1.199 +      if (_last_succ[join] == v_in) {
   1.200 +        up_limit_out = join;
   1.201        } else {
   1.202 -        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
   1.203 -        _thread[u] = right;
   1.204 +        up_limit_in = join;
   1.205 +      }
   1.206 +
   1.207 +      // Update _last_succ from v_in towards the root
   1.208 +      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
   1.209 +           u = _parent[u]) {
   1.210 +        _last_succ[u] = _last_succ[u_out];
   1.211 +      }
   1.212 +      // Update _last_succ from v_out towards the root
   1.213 +      if (join != old_rev_thread && v_in != old_rev_thread) {
   1.214 +        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
   1.215 +             u = _parent[u]) {
   1.216 +          _last_succ[u] = old_rev_thread;
   1.217 +        }
   1.218 +      } else {
   1.219 +        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
   1.220 +             u = _parent[u]) {
   1.221 +          _last_succ[u] = _last_succ[u_out];
   1.222 +        }
   1.223 +      }
   1.224 +
   1.225 +      // Update _succ_num from v_in to join
   1.226 +      for (u = v_in; u != join; u = _parent[u]) {
   1.227 +        _succ_num[u] += old_succ_num;
   1.228 +      }
   1.229 +      // Update _succ_num from v_out to join
   1.230 +      for (u = v_out; u != join; u = _parent[u]) {
   1.231 +        _succ_num[u] -= old_succ_num;
   1.232        }
   1.233      }
   1.234  
   1.235 -    // Update _pred and _forward vectors
   1.236 -    void updatePredArc() {
   1.237 -      int u = u_out, v;
   1.238 -      while (u != u_in) {
   1.239 -        v = _parent[u];
   1.240 -        _pred[u] = _pred[v];
   1.241 -        _forward[u] = !_forward[v];
   1.242 -        u = v;
   1.243 -      }
   1.244 -      _pred[u_in] = in_arc;
   1.245 -      _forward[u_in] = (u_in == _source[in_arc]);
   1.246 -    }
   1.247 -
   1.248 -    // Update _depth and _potential vectors
   1.249 -    void updateDepthPotential() {
   1.250 -      _depth[u_in] = _depth[v_in] + 1;
   1.251 +    // Update potentials
   1.252 +    void updatePotential() {
   1.253        Cost sigma = _forward[u_in] ?
   1.254          _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
   1.255          _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
   1.256 -      _pi[u_in] += sigma;
   1.257 -      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
   1.258 -        _depth[u] = _depth[_parent[u]] + 1;
   1.259 -        if (_depth[u] <= _depth[u_in]) break;
   1.260 -        _pi[u] += sigma;
   1.261 +      if (_succ_num[u_in] > _node_num / 2) {
   1.262 +        // Update in the upper subtree (which contains the root)
   1.263 +        int before = _rev_thread[u_in];
   1.264 +        int after = _thread[_last_succ[u_in]];
   1.265 +        _thread[before] = after;
   1.266 +        _pi[_root] -= sigma;
   1.267 +        for (int u = _thread[_root]; u != _root; u = _thread[u]) {
   1.268 +          _pi[u] -= sigma;
   1.269 +        }
   1.270 +        _thread[before] = u_in;
   1.271 +      } else {
   1.272 +        // Update in the lower subtree (which has been moved)
   1.273 +        int end = _thread[_last_succ[u_in]];
   1.274 +        for (int u = u_in; u != end; u = _thread[u]) {
   1.275 +          _pi[u] += sigma;
   1.276 +        }
   1.277        }
   1.278      }
   1.279  
   1.280 @@ -1139,9 +1208,8 @@
   1.281          bool change = findLeavingArc();
   1.282          changeFlow(change);
   1.283          if (change) {
   1.284 -          updateThreadParent();
   1.285 -          updatePredArc();
   1.286 -          updateDepthPotential();
   1.287 +          updateTreeStructure();
   1.288 +          updatePotential();
   1.289          }
   1.290        }
   1.291