1.1 --- a/lemon/network_simplex.h Tue Jun 22 16:13:00 2010 +0200
1.2 +++ b/lemon/network_simplex.h Sun Aug 22 23:54:10 2010 +0200
1.3 @@ -166,8 +166,9 @@
1.4 typedef std::vector<int> IntVector;
1.5 typedef std::vector<Value> ValueVector;
1.6 typedef std::vector<Cost> CostVector;
1.7 - typedef std::vector<char> BoolVector;
1.8 - // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.9 + typedef std::vector<signed char> CharVector;
1.10 + // Note: vector<signed char> is used instead of vector<ArcState> and
1.11 + // vector<ArcDirection> for efficiency reasons
1.12
1.13 // State constants for arcs
1.14 enum ArcState {
1.15 @@ -176,9 +177,11 @@
1.16 STATE_LOWER = 1
1.17 };
1.18
1.19 - typedef std::vector<signed char> StateVector;
1.20 - // Note: vector<signed char> is used instead of vector<ArcState> for
1.21 - // efficiency reasons
1.22 + // Direction constants for tree arcs
1.23 + enum ArcDirection {
1.24 + DIR_DOWN = -1,
1.25 + DIR_UP = 1
1.26 + };
1.27
1.28 private:
1.29
1.30 @@ -217,15 +220,13 @@
1.31 IntVector _rev_thread;
1.32 IntVector _succ_num;
1.33 IntVector _last_succ;
1.34 + CharVector _pred_dir;
1.35 + CharVector _state;
1.36 IntVector _dirty_revs;
1.37 - BoolVector _forward;
1.38 - StateVector _state;
1.39 int _root;
1.40
1.41 // Temporary data used in the current pivot iteration
1.42 int in_arc, join, u_in, v_in, u_out, v_out;
1.43 - int first, second, right, last;
1.44 - int stem, par_stem, new_stem;
1.45 Value delta;
1.46
1.47 const Value MAX;
1.48 @@ -250,7 +251,7 @@
1.49 const IntVector &_source;
1.50 const IntVector &_target;
1.51 const CostVector &_cost;
1.52 - const StateVector &_state;
1.53 + const CharVector &_state;
1.54 const CostVector &_pi;
1.55 int &_in_arc;
1.56 int _search_arc_num;
1.57 @@ -302,7 +303,7 @@
1.58 const IntVector &_source;
1.59 const IntVector &_target;
1.60 const CostVector &_cost;
1.61 - const StateVector &_state;
1.62 + const CharVector &_state;
1.63 const CostVector &_pi;
1.64 int &_in_arc;
1.65 int _search_arc_num;
1.66 @@ -341,7 +342,7 @@
1.67 const IntVector &_source;
1.68 const IntVector &_target;
1.69 const CostVector &_cost;
1.70 - const StateVector &_state;
1.71 + const CharVector &_state;
1.72 const CostVector &_pi;
1.73 int &_in_arc;
1.74 int _search_arc_num;
1.75 @@ -414,7 +415,7 @@
1.76 const IntVector &_source;
1.77 const IntVector &_target;
1.78 const CostVector &_cost;
1.79 - const StateVector &_state;
1.80 + const CharVector &_state;
1.81 const CostVector &_pi;
1.82 int &_in_arc;
1.83 int _search_arc_num;
1.84 @@ -517,7 +518,7 @@
1.85 const IntVector &_source;
1.86 const IntVector &_target;
1.87 const CostVector &_cost;
1.88 - const StateVector &_state;
1.89 + const CharVector &_state;
1.90 const CostVector &_pi;
1.91 int &_in_arc;
1.92 int _search_arc_num;
1.93 @@ -570,11 +571,13 @@
1.94 bool findEnteringArc() {
1.95 // Check the current candidate list
1.96 int e;
1.97 + Cost c;
1.98 for (int i = 0; i != _curr_length; ++i) {
1.99 e = _candidates[i];
1.100 - _cand_cost[e] = _state[e] *
1.101 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.102 - if (_cand_cost[e] >= 0) {
1.103 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.104 + if (c < 0) {
1.105 + _cand_cost[e] = c;
1.106 + } else {
1.107 _candidates[i--] = _candidates[--_curr_length];
1.108 }
1.109 }
1.110 @@ -584,9 +587,9 @@
1.111 int limit = _head_length;
1.112
1.113 for (e = _next_arc; e != _search_arc_num; ++e) {
1.114 - _cand_cost[e] = _state[e] *
1.115 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.116 - if (_cand_cost[e] < 0) {
1.117 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.118 + if (c < 0) {
1.119 + _cand_cost[e] = c;
1.120 _candidates[_curr_length++] = e;
1.121 }
1.122 if (--cnt == 0) {
1.123 @@ -913,7 +916,7 @@
1.124
1.125 _parent.resize(all_node_num);
1.126 _pred.resize(all_node_num);
1.127 - _forward.resize(all_node_num);
1.128 + _pred_dir.resize(all_node_num);
1.129 _thread.resize(all_node_num);
1.130 _rev_thread.resize(all_node_num);
1.131 _succ_num.resize(all_node_num);
1.132 @@ -1116,14 +1119,14 @@
1.133 _cap[e] = INF;
1.134 _state[e] = STATE_TREE;
1.135 if (_supply[u] >= 0) {
1.136 - _forward[u] = true;
1.137 + _pred_dir[u] = DIR_UP;
1.138 _pi[u] = 0;
1.139 _source[e] = u;
1.140 _target[e] = _root;
1.141 _flow[e] = _supply[u];
1.142 _cost[e] = 0;
1.143 } else {
1.144 - _forward[u] = false;
1.145 + _pred_dir[u] = DIR_DOWN;
1.146 _pi[u] = ART_COST;
1.147 _source[e] = _root;
1.148 _target[e] = u;
1.149 @@ -1143,7 +1146,7 @@
1.150 _succ_num[u] = 1;
1.151 _last_succ[u] = u;
1.152 if (_supply[u] >= 0) {
1.153 - _forward[u] = true;
1.154 + _pred_dir[u] = DIR_UP;
1.155 _pi[u] = 0;
1.156 _pred[u] = e;
1.157 _source[e] = u;
1.158 @@ -1153,7 +1156,7 @@
1.159 _cost[e] = 0;
1.160 _state[e] = STATE_TREE;
1.161 } else {
1.162 - _forward[u] = false;
1.163 + _pred_dir[u] = DIR_DOWN;
1.164 _pi[u] = ART_COST;
1.165 _pred[u] = f;
1.166 _source[f] = _root;
1.167 @@ -1184,7 +1187,7 @@
1.168 _succ_num[u] = 1;
1.169 _last_succ[u] = u;
1.170 if (_supply[u] <= 0) {
1.171 - _forward[u] = false;
1.172 + _pred_dir[u] = DIR_DOWN;
1.173 _pi[u] = 0;
1.174 _pred[u] = e;
1.175 _source[e] = _root;
1.176 @@ -1194,7 +1197,7 @@
1.177 _cost[e] = 0;
1.178 _state[e] = STATE_TREE;
1.179 } else {
1.180 - _forward[u] = true;
1.181 + _pred_dir[u] = DIR_UP;
1.182 _pi[u] = -ART_COST;
1.183 _pred[u] = f;
1.184 _source[f] = u;
1.185 @@ -1237,6 +1240,7 @@
1.186 bool findLeavingArc() {
1.187 // Initialize first and second nodes according to the direction
1.188 // of the cycle
1.189 + int first, second;
1.190 if (_state[in_arc] == STATE_LOWER) {
1.191 first = _source[in_arc];
1.192 second = _target[in_arc];
1.193 @@ -1246,25 +1250,32 @@
1.194 }
1.195 delta = _cap[in_arc];
1.196 int result = 0;
1.197 - Value d;
1.198 + Value c, d;
1.199 int e;
1.200
1.201 - // Search the cycle along the path form the first node to the root
1.202 + // Search the cycle form the first node to the join node
1.203 for (int u = first; u != join; u = _parent[u]) {
1.204 e = _pred[u];
1.205 - d = _forward[u] ?
1.206 - _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
1.207 + d = _flow[e];
1.208 + if (_pred_dir[u] == DIR_DOWN) {
1.209 + c = _cap[e];
1.210 + d = c >= MAX ? INF : c - d;
1.211 + }
1.212 if (d < delta) {
1.213 delta = d;
1.214 u_out = u;
1.215 result = 1;
1.216 }
1.217 }
1.218 - // Search the cycle along the path form the second node to the root
1.219 +
1.220 + // Search the cycle form the second node to the join node
1.221 for (int u = second; u != join; u = _parent[u]) {
1.222 e = _pred[u];
1.223 - d = _forward[u] ?
1.224 - (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
1.225 + d = _flow[e];
1.226 + if (_pred_dir[u] == DIR_UP) {
1.227 + c = _cap[e];
1.228 + d = c >= MAX ? INF : c - d;
1.229 + }
1.230 if (d <= delta) {
1.231 delta = d;
1.232 u_out = u;
1.233 @@ -1289,10 +1300,10 @@
1.234 Value val = _state[in_arc] * delta;
1.235 _flow[in_arc] += val;
1.236 for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1.237 - _flow[_pred[u]] += _forward[u] ? -val : val;
1.238 + _flow[_pred[u]] -= _pred_dir[u] * val;
1.239 }
1.240 for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1.241 - _flow[_pred[u]] += _forward[u] ? val : -val;
1.242 + _flow[_pred[u]] += _pred_dir[u] * val;
1.243 }
1.244 }
1.245 // Update the state of the entering and leaving arcs
1.246 @@ -1307,130 +1318,134 @@
1.247
1.248 // Update the tree structure
1.249 void updateTreeStructure() {
1.250 - int u, w;
1.251 int old_rev_thread = _rev_thread[u_out];
1.252 int old_succ_num = _succ_num[u_out];
1.253 int old_last_succ = _last_succ[u_out];
1.254 v_out = _parent[u_out];
1.255
1.256 - u = _last_succ[u_in]; // the last successor of u_in
1.257 - right = _thread[u]; // the node after it
1.258 + // Check if u_in and u_out coincide
1.259 + if (u_in == u_out) {
1.260 + // Update _parent, _pred, _pred_dir
1.261 + _parent[u_in] = v_in;
1.262 + _pred[u_in] = in_arc;
1.263 + _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1.264
1.265 - // Handle the case when old_rev_thread equals to v_in
1.266 - // (it also means that join and v_out coincide)
1.267 - if (old_rev_thread == v_in) {
1.268 - last = _thread[_last_succ[u_out]];
1.269 + // Update _thread and _rev_thread
1.270 + if (_thread[v_in] != u_out) {
1.271 + int after = _thread[old_last_succ];
1.272 + _thread[old_rev_thread] = after;
1.273 + _rev_thread[after] = old_rev_thread;
1.274 + after = _thread[v_in];
1.275 + _thread[v_in] = u_out;
1.276 + _rev_thread[u_out] = v_in;
1.277 + _thread[old_last_succ] = after;
1.278 + _rev_thread[after] = old_last_succ;
1.279 + }
1.280 } else {
1.281 - last = _thread[v_in];
1.282 - }
1.283 + // Handle the case when old_rev_thread equals to v_in
1.284 + // (it also means that join and v_out coincide)
1.285 + int thread_continue = old_rev_thread == v_in ?
1.286 + _thread[old_last_succ] : _thread[v_in];
1.287
1.288 - // Update _thread and _parent along the stem nodes (i.e. the nodes
1.289 - // between u_in and u_out, whose parent have to be changed)
1.290 - _thread[v_in] = stem = u_in;
1.291 - _dirty_revs.clear();
1.292 - _dirty_revs.push_back(v_in);
1.293 - par_stem = v_in;
1.294 - while (stem != u_out) {
1.295 - // Insert the next stem node into the thread list
1.296 - new_stem = _parent[stem];
1.297 - _thread[u] = new_stem;
1.298 - _dirty_revs.push_back(u);
1.299 + // Update _thread and _parent along the stem nodes (i.e. the nodes
1.300 + // between u_in and u_out, whose parent have to be changed)
1.301 + int stem = u_in; // the current stem node
1.302 + int par_stem = v_in; // the new parent of stem
1.303 + int next_stem; // the next stem node
1.304 + int last = _last_succ[u_in]; // the last successor of stem
1.305 + int before, after = _thread[last];
1.306 + _thread[v_in] = u_in;
1.307 + _dirty_revs.clear();
1.308 + _dirty_revs.push_back(v_in);
1.309 + while (stem != u_out) {
1.310 + // Insert the next stem node into the thread list
1.311 + next_stem = _parent[stem];
1.312 + _thread[last] = next_stem;
1.313 + _dirty_revs.push_back(last);
1.314
1.315 - // Remove the subtree of stem from the thread list
1.316 - w = _rev_thread[stem];
1.317 - _thread[w] = right;
1.318 - _rev_thread[right] = w;
1.319 + // Remove the subtree of stem from the thread list
1.320 + before = _rev_thread[stem];
1.321 + _thread[before] = after;
1.322 + _rev_thread[after] = before;
1.323
1.324 - // Change the parent node and shift stem nodes
1.325 - _parent[stem] = par_stem;
1.326 - par_stem = stem;
1.327 - stem = new_stem;
1.328 + // Change the parent node and shift stem nodes
1.329 + _parent[stem] = par_stem;
1.330 + par_stem = stem;
1.331 + stem = next_stem;
1.332
1.333 - // Update u and right
1.334 - u = _last_succ[stem] == _last_succ[par_stem] ?
1.335 - _rev_thread[par_stem] : _last_succ[stem];
1.336 - right = _thread[u];
1.337 - }
1.338 - _parent[u_out] = par_stem;
1.339 - _thread[u] = last;
1.340 - _rev_thread[last] = u;
1.341 - _last_succ[u_out] = u;
1.342 + // Update last and after
1.343 + last = _last_succ[stem] == _last_succ[par_stem] ?
1.344 + _rev_thread[par_stem] : _last_succ[stem];
1.345 + after = _thread[last];
1.346 + }
1.347 + _parent[u_out] = par_stem;
1.348 + _thread[last] = thread_continue;
1.349 + _rev_thread[thread_continue] = last;
1.350 + _last_succ[u_out] = last;
1.351
1.352 - // Remove the subtree of u_out from the thread list except for
1.353 - // the case when old_rev_thread equals to v_in
1.354 - // (it also means that join and v_out coincide)
1.355 - if (old_rev_thread != v_in) {
1.356 - _thread[old_rev_thread] = right;
1.357 - _rev_thread[right] = old_rev_thread;
1.358 - }
1.359 + // Remove the subtree of u_out from the thread list except for
1.360 + // the case when old_rev_thread equals to v_in
1.361 + if (old_rev_thread != v_in) {
1.362 + _thread[old_rev_thread] = after;
1.363 + _rev_thread[after] = old_rev_thread;
1.364 + }
1.365
1.366 - // Update _rev_thread using the new _thread values
1.367 - for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1.368 - u = _dirty_revs[i];
1.369 - _rev_thread[_thread[u]] = u;
1.370 - }
1.371 + // Update _rev_thread using the new _thread values
1.372 + for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1.373 + int u = _dirty_revs[i];
1.374 + _rev_thread[_thread[u]] = u;
1.375 + }
1.376
1.377 - // Update _pred, _forward, _last_succ and _succ_num for the
1.378 - // stem nodes from u_out to u_in
1.379 - int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1.380 - u = u_out;
1.381 - while (u != u_in) {
1.382 - w = _parent[u];
1.383 - _pred[u] = _pred[w];
1.384 - _forward[u] = !_forward[w];
1.385 - tmp_sc += _succ_num[u] - _succ_num[w];
1.386 - _succ_num[u] = tmp_sc;
1.387 - _last_succ[w] = tmp_ls;
1.388 - u = w;
1.389 - }
1.390 - _pred[u_in] = in_arc;
1.391 - _forward[u_in] = (u_in == _source[in_arc]);
1.392 - _succ_num[u_in] = old_succ_num;
1.393 -
1.394 - // Set limits for updating _last_succ form v_in and v_out
1.395 - // towards the root
1.396 - int up_limit_in = -1;
1.397 - int up_limit_out = -1;
1.398 - if (_last_succ[join] == v_in) {
1.399 - up_limit_out = join;
1.400 - } else {
1.401 - up_limit_in = join;
1.402 + // Update _pred, _pred_dir, _last_succ and _succ_num for the
1.403 + // stem nodes from u_out to u_in
1.404 + int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1.405 + for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
1.406 + _pred[u] = _pred[p];
1.407 + _pred_dir[u] = -_pred_dir[p];
1.408 + tmp_sc += _succ_num[u] - _succ_num[p];
1.409 + _succ_num[u] = tmp_sc;
1.410 + _last_succ[p] = tmp_ls;
1.411 + }
1.412 + _pred[u_in] = in_arc;
1.413 + _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1.414 + _succ_num[u_in] = old_succ_num;
1.415 }
1.416
1.417 // Update _last_succ from v_in towards the root
1.418 - for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1.419 - u = _parent[u]) {
1.420 - _last_succ[u] = _last_succ[u_out];
1.421 + int up_limit_out = _last_succ[join] == v_in ? join : -1;
1.422 + int last_succ_out = _last_succ[u_out];
1.423 + for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
1.424 + _last_succ[u] = last_succ_out;
1.425 }
1.426 +
1.427 // Update _last_succ from v_out towards the root
1.428 if (join != old_rev_thread && v_in != old_rev_thread) {
1.429 - for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.430 + for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.431 u = _parent[u]) {
1.432 _last_succ[u] = old_rev_thread;
1.433 }
1.434 - } else {
1.435 - for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.436 + }
1.437 + else if (last_succ_out != old_last_succ) {
1.438 + for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.439 u = _parent[u]) {
1.440 - _last_succ[u] = _last_succ[u_out];
1.441 + _last_succ[u] = last_succ_out;
1.442 }
1.443 }
1.444
1.445 // Update _succ_num from v_in to join
1.446 - for (u = v_in; u != join; u = _parent[u]) {
1.447 + for (int u = v_in; u != join; u = _parent[u]) {
1.448 _succ_num[u] += old_succ_num;
1.449 }
1.450 // Update _succ_num from v_out to join
1.451 - for (u = v_out; u != join; u = _parent[u]) {
1.452 + for (int u = v_out; u != join; u = _parent[u]) {
1.453 _succ_num[u] -= old_succ_num;
1.454 }
1.455 }
1.456
1.457 - // Update potentials
1.458 + // Update potentials in the subtree that has been moved
1.459 void updatePotential() {
1.460 - Cost sigma = _forward[u_in] ?
1.461 - _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1.462 - _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1.463 - // Update potentials in the subtree, which has been moved
1.464 + Cost sigma = _pi[v_in] - _pi[u_in] -
1.465 + _pred_dir[u_in] * _cost[in_arc];
1.466 int end = _thread[_last_succ[u_in]];
1.467 for (int u = u_in; u != end; u = _thread[u]) {
1.468 _pi[u] += sigma;