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