1.1 --- a/lemon/capacity_scaling.h Tue Feb 09 23:29:51 2010 +0100
1.2 +++ b/lemon/capacity_scaling.h Sat Feb 20 18:39:03 2010 +0100
1.3 @@ -134,9 +134,10 @@
1.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.5
1.6 typedef std::vector<int> IntVector;
1.7 - typedef std::vector<char> BoolVector;
1.8 typedef std::vector<Value> ValueVector;
1.9 typedef std::vector<Cost> CostVector;
1.10 + typedef std::vector<char> BoolVector;
1.11 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.12
1.13 private:
1.14
1.15 @@ -764,15 +765,15 @@
1.16 // Initialize delta value
1.17 if (_factor > 1) {
1.18 // With scaling
1.19 - Value max_sup = 0, max_dem = 0;
1.20 - for (int i = 0; i != _node_num; ++i) {
1.21 + Value max_sup = 0, max_dem = 0, max_cap = 0;
1.22 + for (int i = 0; i != _root; ++i) {
1.23 Value ex = _excess[i];
1.24 if ( ex > max_sup) max_sup = ex;
1.25 if (-ex > max_dem) max_dem = -ex;
1.26 - }
1.27 - Value max_cap = 0;
1.28 - for (int j = 0; j != _res_arc_num; ++j) {
1.29 - if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
1.30 + int last_out = _first_out[i+1] - 1;
1.31 + for (int j = _first_out[i]; j != last_out; ++j) {
1.32 + if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
1.33 + }
1.34 }
1.35 max_sup = std::min(std::min(max_sup, max_dem), max_cap);
1.36 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
2.1 --- a/lemon/cost_scaling.h Tue Feb 09 23:29:51 2010 +0100
2.2 +++ b/lemon/cost_scaling.h Sat Feb 20 18:39:03 2010 +0100
2.3 @@ -197,10 +197,11 @@
2.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
2.5
2.6 typedef std::vector<int> IntVector;
2.7 - typedef std::vector<char> BoolVector;
2.8 typedef std::vector<Value> ValueVector;
2.9 typedef std::vector<Cost> CostVector;
2.10 typedef std::vector<LargeCost> LargeCostVector;
2.11 + typedef std::vector<char> BoolVector;
2.12 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
2.13
2.14 private:
2.15
2.16 @@ -244,6 +245,7 @@
2.17 // Parameters of the problem
2.18 bool _have_lower;
2.19 Value _sum_supply;
2.20 + int _sup_node_num;
2.21
2.22 // Data structures for storing the digraph
2.23 IntNodeMap _node_id;
2.24 @@ -272,6 +274,12 @@
2.25 LargeCost _epsilon;
2.26 int _alpha;
2.27
2.28 + IntVector _buckets;
2.29 + IntVector _bucket_next;
2.30 + IntVector _bucket_prev;
2.31 + IntVector _rank;
2.32 + int _max_rank;
2.33 +
2.34 // Data for a StaticDigraph structure
2.35 typedef std::pair<int, int> IntPair;
2.36 StaticDigraph _sgr;
2.37 @@ -802,6 +810,11 @@
2.38 }
2.39 }
2.40
2.41 + _sup_node_num = 0;
2.42 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.43 + if (sup[n] > 0) ++_sup_node_num;
2.44 + }
2.45 +
2.46 // Find a feasible flow using Circulation
2.47 Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
2.48 circ(_graph, low, cap, sup);
2.49 @@ -836,7 +849,7 @@
2.50 }
2.51 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
2.52 int ra = _reverse[a];
2.53 - _res_cap[a] = 1;
2.54 + _res_cap[a] = 0;
2.55 _res_cap[ra] = 0;
2.56 _cost[a] = 0;
2.57 _cost[ra] = 0;
2.58 @@ -850,7 +863,14 @@
2.59 void start(Method method) {
2.60 // Maximum path length for partial augment
2.61 const int MAX_PATH_LENGTH = 4;
2.62 -
2.63 +
2.64 + // Initialize data structures for buckets
2.65 + _max_rank = _alpha * _res_node_num;
2.66 + _buckets.resize(_max_rank);
2.67 + _bucket_next.resize(_res_node_num + 1);
2.68 + _bucket_prev.resize(_res_node_num + 1);
2.69 + _rank.resize(_res_node_num + 1);
2.70 +
2.71 // Execute the algorithm
2.72 switch (method) {
2.73 case PUSH:
2.74 @@ -889,63 +909,175 @@
2.75 }
2.76 }
2.77 }
2.78 +
2.79 + // Initialize a cost scaling phase
2.80 + void initPhase() {
2.81 + // Saturate arcs not satisfying the optimality condition
2.82 + for (int u = 0; u != _res_node_num; ++u) {
2.83 + int last_out = _first_out[u+1];
2.84 + LargeCost pi_u = _pi[u];
2.85 + for (int a = _first_out[u]; a != last_out; ++a) {
2.86 + int v = _target[a];
2.87 + if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
2.88 + Value delta = _res_cap[a];
2.89 + _excess[u] -= delta;
2.90 + _excess[v] += delta;
2.91 + _res_cap[a] = 0;
2.92 + _res_cap[_reverse[a]] += delta;
2.93 + }
2.94 + }
2.95 + }
2.96 +
2.97 + // Find active nodes (i.e. nodes with positive excess)
2.98 + for (int u = 0; u != _res_node_num; ++u) {
2.99 + if (_excess[u] > 0) _active_nodes.push_back(u);
2.100 + }
2.101 +
2.102 + // Initialize the next arcs
2.103 + for (int u = 0; u != _res_node_num; ++u) {
2.104 + _next_out[u] = _first_out[u];
2.105 + }
2.106 + }
2.107 +
2.108 + // Early termination heuristic
2.109 + bool earlyTermination() {
2.110 + const double EARLY_TERM_FACTOR = 3.0;
2.111 +
2.112 + // Build a static residual graph
2.113 + _arc_vec.clear();
2.114 + _cost_vec.clear();
2.115 + for (int j = 0; j != _res_arc_num; ++j) {
2.116 + if (_res_cap[j] > 0) {
2.117 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
2.118 + _cost_vec.push_back(_cost[j] + 1);
2.119 + }
2.120 + }
2.121 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
2.122 +
2.123 + // Run Bellman-Ford algorithm to check if the current flow is optimal
2.124 + BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
2.125 + bf.init(0);
2.126 + bool done = false;
2.127 + int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
2.128 + for (int i = 0; i < K && !done; ++i) {
2.129 + done = bf.processNextWeakRound();
2.130 + }
2.131 + return done;
2.132 + }
2.133 +
2.134 + // Global potential update heuristic
2.135 + void globalUpdate() {
2.136 + int bucket_end = _root + 1;
2.137 +
2.138 + // Initialize buckets
2.139 + for (int r = 0; r != _max_rank; ++r) {
2.140 + _buckets[r] = bucket_end;
2.141 + }
2.142 + Value total_excess = 0;
2.143 + for (int i = 0; i != _res_node_num; ++i) {
2.144 + if (_excess[i] < 0) {
2.145 + _rank[i] = 0;
2.146 + _bucket_next[i] = _buckets[0];
2.147 + _bucket_prev[_buckets[0]] = i;
2.148 + _buckets[0] = i;
2.149 + } else {
2.150 + total_excess += _excess[i];
2.151 + _rank[i] = _max_rank;
2.152 + }
2.153 + }
2.154 + if (total_excess == 0) return;
2.155 +
2.156 + // Search the buckets
2.157 + int r = 0;
2.158 + for ( ; r != _max_rank; ++r) {
2.159 + while (_buckets[r] != bucket_end) {
2.160 + // Remove the first node from the current bucket
2.161 + int u = _buckets[r];
2.162 + _buckets[r] = _bucket_next[u];
2.163 +
2.164 + // Search the incomming arcs of u
2.165 + LargeCost pi_u = _pi[u];
2.166 + int last_out = _first_out[u+1];
2.167 + for (int a = _first_out[u]; a != last_out; ++a) {
2.168 + int ra = _reverse[a];
2.169 + if (_res_cap[ra] > 0) {
2.170 + int v = _source[ra];
2.171 + int old_rank_v = _rank[v];
2.172 + if (r < old_rank_v) {
2.173 + // Compute the new rank of v
2.174 + LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
2.175 + int new_rank_v = old_rank_v;
2.176 + if (nrc < LargeCost(_max_rank))
2.177 + new_rank_v = r + 1 + int(nrc);
2.178 +
2.179 + // Change the rank of v
2.180 + if (new_rank_v < old_rank_v) {
2.181 + _rank[v] = new_rank_v;
2.182 + _next_out[v] = _first_out[v];
2.183 +
2.184 + // Remove v from its old bucket
2.185 + if (old_rank_v < _max_rank) {
2.186 + if (_buckets[old_rank_v] == v) {
2.187 + _buckets[old_rank_v] = _bucket_next[v];
2.188 + } else {
2.189 + _bucket_next[_bucket_prev[v]] = _bucket_next[v];
2.190 + _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
2.191 + }
2.192 + }
2.193 +
2.194 + // Insert v to its new bucket
2.195 + _bucket_next[v] = _buckets[new_rank_v];
2.196 + _bucket_prev[_buckets[new_rank_v]] = v;
2.197 + _buckets[new_rank_v] = v;
2.198 + }
2.199 + }
2.200 + }
2.201 + }
2.202 +
2.203 + // Finish search if there are no more active nodes
2.204 + if (_excess[u] > 0) {
2.205 + total_excess -= _excess[u];
2.206 + if (total_excess <= 0) break;
2.207 + }
2.208 + }
2.209 + if (total_excess <= 0) break;
2.210 + }
2.211 +
2.212 + // Relabel nodes
2.213 + for (int u = 0; u != _res_node_num; ++u) {
2.214 + int k = std::min(_rank[u], r);
2.215 + if (k > 0) {
2.216 + _pi[u] -= _epsilon * k;
2.217 + _next_out[u] = _first_out[u];
2.218 + }
2.219 + }
2.220 + }
2.221
2.222 /// Execute the algorithm performing augment and relabel operations
2.223 void startAugment(int max_length = std::numeric_limits<int>::max()) {
2.224 // Paramters for heuristics
2.225 - const int BF_HEURISTIC_EPSILON_BOUND = 1000;
2.226 - const int BF_HEURISTIC_BOUND_FACTOR = 3;
2.227 + const int EARLY_TERM_EPSILON_LIMIT = 1000;
2.228 + const double GLOBAL_UPDATE_FACTOR = 3.0;
2.229
2.230 + const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
2.231 + (_res_node_num + _sup_node_num * _sup_node_num));
2.232 + int next_update_limit = global_update_freq;
2.233 +
2.234 + int relabel_cnt = 0;
2.235 +
2.236 // Perform cost scaling phases
2.237 - IntVector pred_arc(_res_node_num);
2.238 - std::vector<int> path_nodes;
2.239 + std::vector<int> path;
2.240 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
2.241 1 : _epsilon / _alpha )
2.242 {
2.243 - // "Early Termination" heuristic: use Bellman-Ford algorithm
2.244 - // to check if the current flow is optimal
2.245 - if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
2.246 - _arc_vec.clear();
2.247 - _cost_vec.clear();
2.248 - for (int j = 0; j != _res_arc_num; ++j) {
2.249 - if (_res_cap[j] > 0) {
2.250 - _arc_vec.push_back(IntPair(_source[j], _target[j]));
2.251 - _cost_vec.push_back(_cost[j] + 1);
2.252 - }
2.253 - }
2.254 - _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
2.255 -
2.256 - BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
2.257 - bf.init(0);
2.258 - bool done = false;
2.259 - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
2.260 - for (int i = 0; i < K && !done; ++i)
2.261 - done = bf.processNextWeakRound();
2.262 - if (done) break;
2.263 - }
2.264 -
2.265 - // Saturate arcs not satisfying the optimality condition
2.266 - for (int a = 0; a != _res_arc_num; ++a) {
2.267 - if (_res_cap[a] > 0 &&
2.268 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
2.269 - Value delta = _res_cap[a];
2.270 - _excess[_source[a]] -= delta;
2.271 - _excess[_target[a]] += delta;
2.272 - _res_cap[a] = 0;
2.273 - _res_cap[_reverse[a]] += delta;
2.274 - }
2.275 + // Early termination heuristic
2.276 + if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
2.277 + if (earlyTermination()) break;
2.278 }
2.279
2.280 - // Find active nodes (i.e. nodes with positive excess)
2.281 - for (int u = 0; u != _res_node_num; ++u) {
2.282 - if (_excess[u] > 0) _active_nodes.push_back(u);
2.283 - }
2.284 -
2.285 - // Initialize the next arcs
2.286 - for (int u = 0; u != _res_node_num; ++u) {
2.287 - _next_out[u] = _first_out[u];
2.288 - }
2.289 -
2.290 + // Initialize current phase
2.291 + initPhase();
2.292 +
2.293 // Perform partial augment and relabel operations
2.294 while (true) {
2.295 // Select an active node (FIFO selection)
2.296 @@ -955,46 +1087,44 @@
2.297 }
2.298 if (_active_nodes.size() == 0) break;
2.299 int start = _active_nodes.front();
2.300 - path_nodes.clear();
2.301 - path_nodes.push_back(start);
2.302
2.303 // Find an augmenting path from the start node
2.304 + path.clear();
2.305 int tip = start;
2.306 - while (_excess[tip] >= 0 &&
2.307 - int(path_nodes.size()) <= max_length) {
2.308 + while (_excess[tip] >= 0 && int(path.size()) < max_length) {
2.309 int u;
2.310 - LargeCost min_red_cost, rc;
2.311 - int last_out = _sum_supply < 0 ?
2.312 - _first_out[tip+1] : _first_out[tip+1] - 1;
2.313 + LargeCost min_red_cost, rc, pi_tip = _pi[tip];
2.314 + int last_out = _first_out[tip+1];
2.315 for (int a = _next_out[tip]; a != last_out; ++a) {
2.316 - if (_res_cap[a] > 0 &&
2.317 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
2.318 - u = _target[a];
2.319 - pred_arc[u] = a;
2.320 + u = _target[a];
2.321 + if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
2.322 + path.push_back(a);
2.323 _next_out[tip] = a;
2.324 tip = u;
2.325 - path_nodes.push_back(tip);
2.326 goto next_step;
2.327 }
2.328 }
2.329
2.330 // Relabel tip node
2.331 - min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
2.332 + min_red_cost = std::numeric_limits<LargeCost>::max();
2.333 + if (tip != start) {
2.334 + int ra = _reverse[path.back()];
2.335 + min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
2.336 + }
2.337 for (int a = _first_out[tip]; a != last_out; ++a) {
2.338 - rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
2.339 + rc = _cost[a] + pi_tip - _pi[_target[a]];
2.340 if (_res_cap[a] > 0 && rc < min_red_cost) {
2.341 min_red_cost = rc;
2.342 }
2.343 }
2.344 _pi[tip] -= min_red_cost + _epsilon;
2.345 -
2.346 - // Reset the next arc of tip
2.347 _next_out[tip] = _first_out[tip];
2.348 + ++relabel_cnt;
2.349
2.350 // Step back
2.351 if (tip != start) {
2.352 - path_nodes.pop_back();
2.353 - tip = path_nodes.back();
2.354 + tip = _source[path.back()];
2.355 + path.pop_back();
2.356 }
2.357
2.358 next_step: ;
2.359 @@ -1002,11 +1132,11 @@
2.360
2.361 // Augment along the found path (as much flow as possible)
2.362 Value delta;
2.363 - int u, v = path_nodes.front(), pa;
2.364 - for (int i = 1; i < int(path_nodes.size()); ++i) {
2.365 + int pa, u, v = start;
2.366 + for (int i = 0; i != int(path.size()); ++i) {
2.367 + pa = path[i];
2.368 u = v;
2.369 - v = path_nodes[i];
2.370 - pa = pred_arc[v];
2.371 + v = _target[pa];
2.372 delta = std::min(_res_cap[pa], _excess[u]);
2.373 _res_cap[pa] -= delta;
2.374 _res_cap[_reverse[pa]] += delta;
2.375 @@ -1015,6 +1145,12 @@
2.376 if (_excess[v] > 0 && _excess[v] <= delta)
2.377 _active_nodes.push_back(v);
2.378 }
2.379 +
2.380 + // Global update heuristic
2.381 + if (relabel_cnt >= next_update_limit) {
2.382 + globalUpdate();
2.383 + next_update_limit += global_update_freq;
2.384 + }
2.385 }
2.386 }
2.387 }
2.388 @@ -1022,98 +1158,70 @@
2.389 /// Execute the algorithm performing push and relabel operations
2.390 void startPush() {
2.391 // Paramters for heuristics
2.392 - const int BF_HEURISTIC_EPSILON_BOUND = 1000;
2.393 - const int BF_HEURISTIC_BOUND_FACTOR = 3;
2.394 + const int EARLY_TERM_EPSILON_LIMIT = 1000;
2.395 + const double GLOBAL_UPDATE_FACTOR = 2.0;
2.396
2.397 + const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
2.398 + (_res_node_num + _sup_node_num * _sup_node_num));
2.399 + int next_update_limit = global_update_freq;
2.400 +
2.401 + int relabel_cnt = 0;
2.402 +
2.403 // Perform cost scaling phases
2.404 BoolVector hyper(_res_node_num, false);
2.405 + LargeCostVector hyper_cost(_res_node_num);
2.406 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
2.407 1 : _epsilon / _alpha )
2.408 {
2.409 - // "Early Termination" heuristic: use Bellman-Ford algorithm
2.410 - // to check if the current flow is optimal
2.411 - if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
2.412 - _arc_vec.clear();
2.413 - _cost_vec.clear();
2.414 - for (int j = 0; j != _res_arc_num; ++j) {
2.415 - if (_res_cap[j] > 0) {
2.416 - _arc_vec.push_back(IntPair(_source[j], _target[j]));
2.417 - _cost_vec.push_back(_cost[j] + 1);
2.418 - }
2.419 - }
2.420 - _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
2.421 -
2.422 - BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
2.423 - bf.init(0);
2.424 - bool done = false;
2.425 - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
2.426 - for (int i = 0; i < K && !done; ++i)
2.427 - done = bf.processNextWeakRound();
2.428 - if (done) break;
2.429 + // Early termination heuristic
2.430 + if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
2.431 + if (earlyTermination()) break;
2.432 }
2.433 -
2.434 - // Saturate arcs not satisfying the optimality condition
2.435 - for (int a = 0; a != _res_arc_num; ++a) {
2.436 - if (_res_cap[a] > 0 &&
2.437 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
2.438 - Value delta = _res_cap[a];
2.439 - _excess[_source[a]] -= delta;
2.440 - _excess[_target[a]] += delta;
2.441 - _res_cap[a] = 0;
2.442 - _res_cap[_reverse[a]] += delta;
2.443 - }
2.444 - }
2.445 -
2.446 - // Find active nodes (i.e. nodes with positive excess)
2.447 - for (int u = 0; u != _res_node_num; ++u) {
2.448 - if (_excess[u] > 0) _active_nodes.push_back(u);
2.449 - }
2.450 -
2.451 - // Initialize the next arcs
2.452 - for (int u = 0; u != _res_node_num; ++u) {
2.453 - _next_out[u] = _first_out[u];
2.454 - }
2.455 +
2.456 + // Initialize current phase
2.457 + initPhase();
2.458
2.459 // Perform push and relabel operations
2.460 while (_active_nodes.size() > 0) {
2.461 - LargeCost min_red_cost, rc;
2.462 + LargeCost min_red_cost, rc, pi_n;
2.463 Value delta;
2.464 int n, t, a, last_out = _res_arc_num;
2.465
2.466 + next_node:
2.467 // Select an active node (FIFO selection)
2.468 - next_node:
2.469 n = _active_nodes.front();
2.470 - last_out = _sum_supply < 0 ?
2.471 - _first_out[n+1] : _first_out[n+1] - 1;
2.472 -
2.473 + last_out = _first_out[n+1];
2.474 + pi_n = _pi[n];
2.475 +
2.476 // Perform push operations if there are admissible arcs
2.477 if (_excess[n] > 0) {
2.478 for (a = _next_out[n]; a != last_out; ++a) {
2.479 if (_res_cap[a] > 0 &&
2.480 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
2.481 + _cost[a] + pi_n - _pi[_target[a]] < 0) {
2.482 delta = std::min(_res_cap[a], _excess[n]);
2.483 t = _target[a];
2.484
2.485 // Push-look-ahead heuristic
2.486 Value ahead = -_excess[t];
2.487 - int last_out_t = _sum_supply < 0 ?
2.488 - _first_out[t+1] : _first_out[t+1] - 1;
2.489 + int last_out_t = _first_out[t+1];
2.490 + LargeCost pi_t = _pi[t];
2.491 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
2.492 if (_res_cap[ta] > 0 &&
2.493 - _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
2.494 + _cost[ta] + pi_t - _pi[_target[ta]] < 0)
2.495 ahead += _res_cap[ta];
2.496 if (ahead >= delta) break;
2.497 }
2.498 if (ahead < 0) ahead = 0;
2.499
2.500 // Push flow along the arc
2.501 - if (ahead < delta) {
2.502 + if (ahead < delta && !hyper[t]) {
2.503 _res_cap[a] -= ahead;
2.504 _res_cap[_reverse[a]] += ahead;
2.505 _excess[n] -= ahead;
2.506 _excess[t] += ahead;
2.507 _active_nodes.push_front(t);
2.508 hyper[t] = true;
2.509 + hyper_cost[t] = _cost[a] + pi_n - pi_t;
2.510 _next_out[n] = a;
2.511 goto next_node;
2.512 } else {
2.513 @@ -1136,18 +1244,18 @@
2.514
2.515 // Relabel the node if it is still active (or hyper)
2.516 if (_excess[n] > 0 || hyper[n]) {
2.517 - min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
2.518 + min_red_cost = hyper[n] ? -hyper_cost[n] :
2.519 + std::numeric_limits<LargeCost>::max();
2.520 for (int a = _first_out[n]; a != last_out; ++a) {
2.521 - rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
2.522 + rc = _cost[a] + pi_n - _pi[_target[a]];
2.523 if (_res_cap[a] > 0 && rc < min_red_cost) {
2.524 min_red_cost = rc;
2.525 }
2.526 }
2.527 _pi[n] -= min_red_cost + _epsilon;
2.528 + _next_out[n] = _first_out[n];
2.529 hyper[n] = false;
2.530 -
2.531 - // Reset the next arc
2.532 - _next_out[n] = _first_out[n];
2.533 + ++relabel_cnt;
2.534 }
2.535
2.536 // Remove nodes that are not active nor hyper
2.537 @@ -1157,6 +1265,14 @@
2.538 !hyper[_active_nodes.front()] ) {
2.539 _active_nodes.pop_front();
2.540 }
2.541 +
2.542 + // Global update heuristic
2.543 + if (relabel_cnt >= next_update_limit) {
2.544 + globalUpdate();
2.545 + for (int u = 0; u != _res_node_num; ++u)
2.546 + hyper[u] = false;
2.547 + next_update_limit += global_update_freq;
2.548 + }
2.549 }
2.550 }
2.551 }
3.1 --- a/lemon/cycle_canceling.h Tue Feb 09 23:29:51 2010 +0100
3.2 +++ b/lemon/cycle_canceling.h Sat Feb 20 18:39:03 2010 +0100
3.3 @@ -144,10 +144,11 @@
3.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
3.5
3.6 typedef std::vector<int> IntVector;
3.7 - typedef std::vector<char> CharVector;
3.8 typedef std::vector<double> DoubleVector;
3.9 typedef std::vector<Value> ValueVector;
3.10 typedef std::vector<Cost> CostVector;
3.11 + typedef std::vector<char> BoolVector;
3.12 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
3.13
3.14 private:
3.15
3.16 @@ -198,7 +199,7 @@
3.17 IntArcMap _arc_idf;
3.18 IntArcMap _arc_idb;
3.19 IntVector _first_out;
3.20 - CharVector _forward;
3.21 + BoolVector _forward;
3.22 IntVector _source;
3.23 IntVector _target;
3.24 IntVector _reverse;
3.25 @@ -933,8 +934,8 @@
3.26 // Contruct auxiliary data vectors
3.27 DoubleVector pi(_res_node_num, 0.0);
3.28 IntVector level(_res_node_num);
3.29 - CharVector reached(_res_node_num);
3.30 - CharVector processed(_res_node_num);
3.31 + BoolVector reached(_res_node_num);
3.32 + BoolVector processed(_res_node_num);
3.33 IntVector pred_node(_res_node_num);
3.34 IntVector pred_arc(_res_node_num);
3.35 std::vector<int> stack(_res_node_num);
4.1 --- a/lemon/network_simplex.h Tue Feb 09 23:29:51 2010 +0100
4.2 +++ b/lemon/network_simplex.h Sat Feb 20 18:39:03 2010 +0100
4.3 @@ -164,9 +164,10 @@
4.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
4.5
4.6 typedef std::vector<int> IntVector;
4.7 - typedef std::vector<char> CharVector;
4.8 typedef std::vector<Value> ValueVector;
4.9 typedef std::vector<Cost> CostVector;
4.10 + typedef std::vector<char> BoolVector;
4.11 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
4.12
4.13 // State constants for arcs
4.14 enum ArcStateEnum {
4.15 @@ -212,8 +213,8 @@
4.16 IntVector _succ_num;
4.17 IntVector _last_succ;
4.18 IntVector _dirty_revs;
4.19 - CharVector _forward;
4.20 - CharVector _state;
4.21 + BoolVector _forward;
4.22 + BoolVector _state;
4.23 int _root;
4.24
4.25 // Temporary data used in the current pivot iteration
4.26 @@ -244,7 +245,7 @@
4.27 const IntVector &_source;
4.28 const IntVector &_target;
4.29 const CostVector &_cost;
4.30 - const CharVector &_state;
4.31 + const BoolVector &_state;
4.32 const CostVector &_pi;
4.33 int &_in_arc;
4.34 int _search_arc_num;
4.35 @@ -265,7 +266,7 @@
4.36 // Find next entering arc
4.37 bool findEnteringArc() {
4.38 Cost c;
4.39 - for (int e = _next_arc; e < _search_arc_num; ++e) {
4.40 + for (int e = _next_arc; e != _search_arc_num; ++e) {
4.41 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.42 if (c < 0) {
4.43 _in_arc = e;
4.44 @@ -273,7 +274,7 @@
4.45 return true;
4.46 }
4.47 }
4.48 - for (int e = 0; e < _next_arc; ++e) {
4.49 + for (int e = 0; e != _next_arc; ++e) {
4.50 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.51 if (c < 0) {
4.52 _in_arc = e;
4.53 @@ -296,7 +297,7 @@
4.54 const IntVector &_source;
4.55 const IntVector &_target;
4.56 const CostVector &_cost;
4.57 - const CharVector &_state;
4.58 + const BoolVector &_state;
4.59 const CostVector &_pi;
4.60 int &_in_arc;
4.61 int _search_arc_num;
4.62 @@ -313,7 +314,7 @@
4.63 // Find next entering arc
4.64 bool findEnteringArc() {
4.65 Cost c, min = 0;
4.66 - for (int e = 0; e < _search_arc_num; ++e) {
4.67 + for (int e = 0; e != _search_arc_num; ++e) {
4.68 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.69 if (c < min) {
4.70 min = c;
4.71 @@ -335,7 +336,7 @@
4.72 const IntVector &_source;
4.73 const IntVector &_target;
4.74 const CostVector &_cost;
4.75 - const CharVector &_state;
4.76 + const BoolVector &_state;
4.77 const CostVector &_pi;
4.78 int &_in_arc;
4.79 int _search_arc_num;
4.80 @@ -354,7 +355,7 @@
4.81 _next_arc(0)
4.82 {
4.83 // The main parameters of the pivot rule
4.84 - const double BLOCK_SIZE_FACTOR = 0.5;
4.85 + const double BLOCK_SIZE_FACTOR = 1.0;
4.86 const int MIN_BLOCK_SIZE = 10;
4.87
4.88 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
4.89 @@ -367,7 +368,7 @@
4.90 Cost c, min = 0;
4.91 int cnt = _block_size;
4.92 int e;
4.93 - for (e = _next_arc; e < _search_arc_num; ++e) {
4.94 + for (e = _next_arc; e != _search_arc_num; ++e) {
4.95 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.96 if (c < min) {
4.97 min = c;
4.98 @@ -378,7 +379,7 @@
4.99 cnt = _block_size;
4.100 }
4.101 }
4.102 - for (e = 0; e < _next_arc; ++e) {
4.103 + for (e = 0; e != _next_arc; ++e) {
4.104 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.105 if (c < min) {
4.106 min = c;
4.107 @@ -408,7 +409,7 @@
4.108 const IntVector &_source;
4.109 const IntVector &_target;
4.110 const CostVector &_cost;
4.111 - const CharVector &_state;
4.112 + const BoolVector &_state;
4.113 const CostVector &_pi;
4.114 int &_in_arc;
4.115 int _search_arc_num;
4.116 @@ -469,7 +470,7 @@
4.117 // Major iteration: build a new candidate list
4.118 min = 0;
4.119 _curr_length = 0;
4.120 - for (e = _next_arc; e < _search_arc_num; ++e) {
4.121 + for (e = _next_arc; e != _search_arc_num; ++e) {
4.122 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.123 if (c < 0) {
4.124 _candidates[_curr_length++] = e;
4.125 @@ -480,7 +481,7 @@
4.126 if (_curr_length == _list_length) goto search_end;
4.127 }
4.128 }
4.129 - for (e = 0; e < _next_arc; ++e) {
4.130 + for (e = 0; e != _next_arc; ++e) {
4.131 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.132 if (c < 0) {
4.133 _candidates[_curr_length++] = e;
4.134 @@ -511,7 +512,7 @@
4.135 const IntVector &_source;
4.136 const IntVector &_target;
4.137 const CostVector &_cost;
4.138 - const CharVector &_state;
4.139 + const BoolVector &_state;
4.140 const CostVector &_pi;
4.141 int &_in_arc;
4.142 int _search_arc_num;
4.143 @@ -564,7 +565,7 @@
4.144 bool findEnteringArc() {
4.145 // Check the current candidate list
4.146 int e;
4.147 - for (int i = 0; i < _curr_length; ++i) {
4.148 + for (int i = 0; i != _curr_length; ++i) {
4.149 e = _candidates[i];
4.150 _cand_cost[e] = _state[e] *
4.151 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.152 @@ -577,7 +578,7 @@
4.153 int cnt = _block_size;
4.154 int limit = _head_length;
4.155
4.156 - for (e = _next_arc; e < _search_arc_num; ++e) {
4.157 + for (e = _next_arc; e != _search_arc_num; ++e) {
4.158 _cand_cost[e] = _state[e] *
4.159 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.160 if (_cand_cost[e] < 0) {
4.161 @@ -589,7 +590,7 @@
4.162 cnt = _block_size;
4.163 }
4.164 }
4.165 - for (e = 0; e < _next_arc; ++e) {
4.166 + for (e = 0; e != _next_arc; ++e) {
4.167 _cand_cost[e] = _state[e] *
4.168 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.169 if (_cand_cost[e] < 0) {
4.170 @@ -1328,7 +1329,7 @@
4.171 }
4.172
4.173 // Update _rev_thread using the new _thread values
4.174 - for (int i = 0; i < int(_dirty_revs.size()); ++i) {
4.175 + for (int i = 0; i != int(_dirty_revs.size()); ++i) {
4.176 u = _dirty_revs[i];
4.177 _rev_thread[_thread[u]] = u;
4.178 }
4.179 @@ -1400,6 +1401,100 @@
4.180 }
4.181 }
4.182
4.183 + // Heuristic initial pivots
4.184 + bool initialPivots() {
4.185 + Value curr, total = 0;
4.186 + std::vector<Node> supply_nodes, demand_nodes;
4.187 + for (NodeIt u(_graph); u != INVALID; ++u) {
4.188 + curr = _supply[_node_id[u]];
4.189 + if (curr > 0) {
4.190 + total += curr;
4.191 + supply_nodes.push_back(u);
4.192 + }
4.193 + else if (curr < 0) {
4.194 + demand_nodes.push_back(u);
4.195 + }
4.196 + }
4.197 + if (_sum_supply > 0) total -= _sum_supply;
4.198 + if (total <= 0) return true;
4.199 +
4.200 + IntVector arc_vector;
4.201 + if (_sum_supply >= 0) {
4.202 + if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
4.203 + // Perform a reverse graph search from the sink to the source
4.204 + typename GR::template NodeMap<bool> reached(_graph, false);
4.205 + Node s = supply_nodes[0], t = demand_nodes[0];
4.206 + std::vector<Node> stack;
4.207 + reached[t] = true;
4.208 + stack.push_back(t);
4.209 + while (!stack.empty()) {
4.210 + Node u, v = stack.back();
4.211 + stack.pop_back();
4.212 + if (v == s) break;
4.213 + for (InArcIt a(_graph, v); a != INVALID; ++a) {
4.214 + if (reached[u = _graph.source(a)]) continue;
4.215 + int j = _arc_id[a];
4.216 + if (_cap[j] >= total) {
4.217 + arc_vector.push_back(j);
4.218 + reached[u] = true;
4.219 + stack.push_back(u);
4.220 + }
4.221 + }
4.222 + }
4.223 + } else {
4.224 + // Find the min. cost incomming arc for each demand node
4.225 + for (int i = 0; i != int(demand_nodes.size()); ++i) {
4.226 + Node v = demand_nodes[i];
4.227 + Cost c, min_cost = std::numeric_limits<Cost>::max();
4.228 + Arc min_arc = INVALID;
4.229 + for (InArcIt a(_graph, v); a != INVALID; ++a) {
4.230 + c = _cost[_arc_id[a]];
4.231 + if (c < min_cost) {
4.232 + min_cost = c;
4.233 + min_arc = a;
4.234 + }
4.235 + }
4.236 + if (min_arc != INVALID) {
4.237 + arc_vector.push_back(_arc_id[min_arc]);
4.238 + }
4.239 + }
4.240 + }
4.241 + } else {
4.242 + // Find the min. cost outgoing arc for each supply node
4.243 + for (int i = 0; i != int(supply_nodes.size()); ++i) {
4.244 + Node u = supply_nodes[i];
4.245 + Cost c, min_cost = std::numeric_limits<Cost>::max();
4.246 + Arc min_arc = INVALID;
4.247 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
4.248 + c = _cost[_arc_id[a]];
4.249 + if (c < min_cost) {
4.250 + min_cost = c;
4.251 + min_arc = a;
4.252 + }
4.253 + }
4.254 + if (min_arc != INVALID) {
4.255 + arc_vector.push_back(_arc_id[min_arc]);
4.256 + }
4.257 + }
4.258 + }
4.259 +
4.260 + // Perform heuristic initial pivots
4.261 + for (int i = 0; i != int(arc_vector.size()); ++i) {
4.262 + in_arc = arc_vector[i];
4.263 + if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
4.264 + _pi[_target[in_arc]]) >= 0) continue;
4.265 + findJoinNode();
4.266 + bool change = findLeavingArc();
4.267 + if (delta >= MAX) return false;
4.268 + changeFlow(change);
4.269 + if (change) {
4.270 + updateTreeStructure();
4.271 + updatePotential();
4.272 + }
4.273 + }
4.274 + return true;
4.275 + }
4.276 +
4.277 // Execute the algorithm
4.278 ProblemType start(PivotRule pivot_rule) {
4.279 // Select the pivot rule implementation
4.280 @@ -1422,6 +1517,9 @@
4.281 ProblemType start() {
4.282 PivotRuleImpl pivot(*this);
4.283
4.284 + // Perform heuristic initial pivots
4.285 + if (!initialPivots()) return UNBOUNDED;
4.286 +
4.287 // Execute the Network Simplex algorithm
4.288 while (pivot.findEnteringArc()) {
4.289 findJoinNode();