1.1 --- a/lemon/cost_scaling.h Sun Feb 21 18:55:30 2010 +0100
1.2 +++ b/lemon/cost_scaling.h Fri Feb 26 14:00:20 2010 +0100
1.3 @@ -201,10 +201,11 @@
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<LargeCost> LargeCostVector;
1.11 + typedef std::vector<char> BoolVector;
1.12 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.13
1.14 private:
1.15
1.16 @@ -248,6 +249,7 @@
1.17 // Parameters of the problem
1.18 bool _have_lower;
1.19 Value _sum_supply;
1.20 + int _sup_node_num;
1.21
1.22 // Data structures for storing the digraph
1.23 IntNodeMap _node_id;
1.24 @@ -276,6 +278,12 @@
1.25 LargeCost _epsilon;
1.26 int _alpha;
1.27
1.28 + IntVector _buckets;
1.29 + IntVector _bucket_next;
1.30 + IntVector _bucket_prev;
1.31 + IntVector _rank;
1.32 + int _max_rank;
1.33 +
1.34 // Data for a StaticDigraph structure
1.35 typedef std::pair<int, int> IntPair;
1.36 StaticDigraph _sgr;
1.37 @@ -828,6 +836,11 @@
1.38 }
1.39 }
1.40
1.41 + _sup_node_num = 0;
1.42 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.43 + if (sup[n] > 0) ++_sup_node_num;
1.44 + }
1.45 +
1.46 // Find a feasible flow using Circulation
1.47 Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
1.48 circ(_graph, low, cap, sup);
1.49 @@ -862,7 +875,7 @@
1.50 }
1.51 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.52 int ra = _reverse[a];
1.53 - _res_cap[a] = 1;
1.54 + _res_cap[a] = 0;
1.55 _res_cap[ra] = 0;
1.56 _cost[a] = 0;
1.57 _cost[ra] = 0;
1.58 @@ -876,7 +889,14 @@
1.59 void start(Method method) {
1.60 // Maximum path length for partial augment
1.61 const int MAX_PATH_LENGTH = 4;
1.62 -
1.63 +
1.64 + // Initialize data structures for buckets
1.65 + _max_rank = _alpha * _res_node_num;
1.66 + _buckets.resize(_max_rank);
1.67 + _bucket_next.resize(_res_node_num + 1);
1.68 + _bucket_prev.resize(_res_node_num + 1);
1.69 + _rank.resize(_res_node_num + 1);
1.70 +
1.71 // Execute the algorithm
1.72 switch (method) {
1.73 case PUSH:
1.74 @@ -915,63 +935,175 @@
1.75 }
1.76 }
1.77 }
1.78 +
1.79 + // Initialize a cost scaling phase
1.80 + void initPhase() {
1.81 + // Saturate arcs not satisfying the optimality condition
1.82 + for (int u = 0; u != _res_node_num; ++u) {
1.83 + int last_out = _first_out[u+1];
1.84 + LargeCost pi_u = _pi[u];
1.85 + for (int a = _first_out[u]; a != last_out; ++a) {
1.86 + int v = _target[a];
1.87 + if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
1.88 + Value delta = _res_cap[a];
1.89 + _excess[u] -= delta;
1.90 + _excess[v] += delta;
1.91 + _res_cap[a] = 0;
1.92 + _res_cap[_reverse[a]] += delta;
1.93 + }
1.94 + }
1.95 + }
1.96 +
1.97 + // Find active nodes (i.e. nodes with positive excess)
1.98 + for (int u = 0; u != _res_node_num; ++u) {
1.99 + if (_excess[u] > 0) _active_nodes.push_back(u);
1.100 + }
1.101 +
1.102 + // Initialize the next arcs
1.103 + for (int u = 0; u != _res_node_num; ++u) {
1.104 + _next_out[u] = _first_out[u];
1.105 + }
1.106 + }
1.107 +
1.108 + // Early termination heuristic
1.109 + bool earlyTermination() {
1.110 + const double EARLY_TERM_FACTOR = 3.0;
1.111 +
1.112 + // Build a static residual graph
1.113 + _arc_vec.clear();
1.114 + _cost_vec.clear();
1.115 + for (int j = 0; j != _res_arc_num; ++j) {
1.116 + if (_res_cap[j] > 0) {
1.117 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.118 + _cost_vec.push_back(_cost[j] + 1);
1.119 + }
1.120 + }
1.121 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.122 +
1.123 + // Run Bellman-Ford algorithm to check if the current flow is optimal
1.124 + BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.125 + bf.init(0);
1.126 + bool done = false;
1.127 + int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
1.128 + for (int i = 0; i < K && !done; ++i) {
1.129 + done = bf.processNextWeakRound();
1.130 + }
1.131 + return done;
1.132 + }
1.133 +
1.134 + // Global potential update heuristic
1.135 + void globalUpdate() {
1.136 + int bucket_end = _root + 1;
1.137 +
1.138 + // Initialize buckets
1.139 + for (int r = 0; r != _max_rank; ++r) {
1.140 + _buckets[r] = bucket_end;
1.141 + }
1.142 + Value total_excess = 0;
1.143 + for (int i = 0; i != _res_node_num; ++i) {
1.144 + if (_excess[i] < 0) {
1.145 + _rank[i] = 0;
1.146 + _bucket_next[i] = _buckets[0];
1.147 + _bucket_prev[_buckets[0]] = i;
1.148 + _buckets[0] = i;
1.149 + } else {
1.150 + total_excess += _excess[i];
1.151 + _rank[i] = _max_rank;
1.152 + }
1.153 + }
1.154 + if (total_excess == 0) return;
1.155 +
1.156 + // Search the buckets
1.157 + int r = 0;
1.158 + for ( ; r != _max_rank; ++r) {
1.159 + while (_buckets[r] != bucket_end) {
1.160 + // Remove the first node from the current bucket
1.161 + int u = _buckets[r];
1.162 + _buckets[r] = _bucket_next[u];
1.163 +
1.164 + // Search the incomming arcs of u
1.165 + LargeCost pi_u = _pi[u];
1.166 + int last_out = _first_out[u+1];
1.167 + for (int a = _first_out[u]; a != last_out; ++a) {
1.168 + int ra = _reverse[a];
1.169 + if (_res_cap[ra] > 0) {
1.170 + int v = _source[ra];
1.171 + int old_rank_v = _rank[v];
1.172 + if (r < old_rank_v) {
1.173 + // Compute the new rank of v
1.174 + LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
1.175 + int new_rank_v = old_rank_v;
1.176 + if (nrc < LargeCost(_max_rank))
1.177 + new_rank_v = r + 1 + int(nrc);
1.178 +
1.179 + // Change the rank of v
1.180 + if (new_rank_v < old_rank_v) {
1.181 + _rank[v] = new_rank_v;
1.182 + _next_out[v] = _first_out[v];
1.183 +
1.184 + // Remove v from its old bucket
1.185 + if (old_rank_v < _max_rank) {
1.186 + if (_buckets[old_rank_v] == v) {
1.187 + _buckets[old_rank_v] = _bucket_next[v];
1.188 + } else {
1.189 + _bucket_next[_bucket_prev[v]] = _bucket_next[v];
1.190 + _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
1.191 + }
1.192 + }
1.193 +
1.194 + // Insert v to its new bucket
1.195 + _bucket_next[v] = _buckets[new_rank_v];
1.196 + _bucket_prev[_buckets[new_rank_v]] = v;
1.197 + _buckets[new_rank_v] = v;
1.198 + }
1.199 + }
1.200 + }
1.201 + }
1.202 +
1.203 + // Finish search if there are no more active nodes
1.204 + if (_excess[u] > 0) {
1.205 + total_excess -= _excess[u];
1.206 + if (total_excess <= 0) break;
1.207 + }
1.208 + }
1.209 + if (total_excess <= 0) break;
1.210 + }
1.211 +
1.212 + // Relabel nodes
1.213 + for (int u = 0; u != _res_node_num; ++u) {
1.214 + int k = std::min(_rank[u], r);
1.215 + if (k > 0) {
1.216 + _pi[u] -= _epsilon * k;
1.217 + _next_out[u] = _first_out[u];
1.218 + }
1.219 + }
1.220 + }
1.221
1.222 /// Execute the algorithm performing augment and relabel operations
1.223 void startAugment(int max_length = std::numeric_limits<int>::max()) {
1.224 // Paramters for heuristics
1.225 - const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.226 - const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.227 + const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.228 + const double GLOBAL_UPDATE_FACTOR = 3.0;
1.229
1.230 + const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.231 + (_res_node_num + _sup_node_num * _sup_node_num));
1.232 + int next_update_limit = global_update_freq;
1.233 +
1.234 + int relabel_cnt = 0;
1.235 +
1.236 // Perform cost scaling phases
1.237 - IntVector pred_arc(_res_node_num);
1.238 - std::vector<int> path_nodes;
1.239 + std::vector<int> path;
1.240 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.241 1 : _epsilon / _alpha )
1.242 {
1.243 - // "Early Termination" heuristic: use Bellman-Ford algorithm
1.244 - // to check if the current flow is optimal
1.245 - if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.246 - _arc_vec.clear();
1.247 - _cost_vec.clear();
1.248 - for (int j = 0; j != _res_arc_num; ++j) {
1.249 - if (_res_cap[j] > 0) {
1.250 - _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.251 - _cost_vec.push_back(_cost[j] + 1);
1.252 - }
1.253 - }
1.254 - _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.255 -
1.256 - BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.257 - bf.init(0);
1.258 - bool done = false;
1.259 - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1.260 - for (int i = 0; i < K && !done; ++i)
1.261 - done = bf.processNextWeakRound();
1.262 - if (done) break;
1.263 - }
1.264 -
1.265 - // Saturate arcs not satisfying the optimality condition
1.266 - for (int a = 0; a != _res_arc_num; ++a) {
1.267 - if (_res_cap[a] > 0 &&
1.268 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.269 - Value delta = _res_cap[a];
1.270 - _excess[_source[a]] -= delta;
1.271 - _excess[_target[a]] += delta;
1.272 - _res_cap[a] = 0;
1.273 - _res_cap[_reverse[a]] += delta;
1.274 - }
1.275 + // Early termination heuristic
1.276 + if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.277 + if (earlyTermination()) break;
1.278 }
1.279
1.280 - // Find active nodes (i.e. nodes with positive excess)
1.281 - for (int u = 0; u != _res_node_num; ++u) {
1.282 - if (_excess[u] > 0) _active_nodes.push_back(u);
1.283 - }
1.284 -
1.285 - // Initialize the next arcs
1.286 - for (int u = 0; u != _res_node_num; ++u) {
1.287 - _next_out[u] = _first_out[u];
1.288 - }
1.289 -
1.290 + // Initialize current phase
1.291 + initPhase();
1.292 +
1.293 // Perform partial augment and relabel operations
1.294 while (true) {
1.295 // Select an active node (FIFO selection)
1.296 @@ -981,46 +1113,44 @@
1.297 }
1.298 if (_active_nodes.size() == 0) break;
1.299 int start = _active_nodes.front();
1.300 - path_nodes.clear();
1.301 - path_nodes.push_back(start);
1.302
1.303 // Find an augmenting path from the start node
1.304 + path.clear();
1.305 int tip = start;
1.306 - while (_excess[tip] >= 0 &&
1.307 - int(path_nodes.size()) <= max_length) {
1.308 + while (_excess[tip] >= 0 && int(path.size()) < max_length) {
1.309 int u;
1.310 - LargeCost min_red_cost, rc;
1.311 - int last_out = _sum_supply < 0 ?
1.312 - _first_out[tip+1] : _first_out[tip+1] - 1;
1.313 + LargeCost min_red_cost, rc, pi_tip = _pi[tip];
1.314 + int last_out = _first_out[tip+1];
1.315 for (int a = _next_out[tip]; a != last_out; ++a) {
1.316 - if (_res_cap[a] > 0 &&
1.317 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.318 - u = _target[a];
1.319 - pred_arc[u] = a;
1.320 + u = _target[a];
1.321 + if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
1.322 + path.push_back(a);
1.323 _next_out[tip] = a;
1.324 tip = u;
1.325 - path_nodes.push_back(tip);
1.326 goto next_step;
1.327 }
1.328 }
1.329
1.330 // Relabel tip node
1.331 - min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1.332 + min_red_cost = std::numeric_limits<LargeCost>::max();
1.333 + if (tip != start) {
1.334 + int ra = _reverse[path.back()];
1.335 + min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
1.336 + }
1.337 for (int a = _first_out[tip]; a != last_out; ++a) {
1.338 - rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1.339 + rc = _cost[a] + pi_tip - _pi[_target[a]];
1.340 if (_res_cap[a] > 0 && rc < min_red_cost) {
1.341 min_red_cost = rc;
1.342 }
1.343 }
1.344 _pi[tip] -= min_red_cost + _epsilon;
1.345 -
1.346 - // Reset the next arc of tip
1.347 _next_out[tip] = _first_out[tip];
1.348 + ++relabel_cnt;
1.349
1.350 // Step back
1.351 if (tip != start) {
1.352 - path_nodes.pop_back();
1.353 - tip = path_nodes.back();
1.354 + tip = _source[path.back()];
1.355 + path.pop_back();
1.356 }
1.357
1.358 next_step: ;
1.359 @@ -1028,11 +1158,11 @@
1.360
1.361 // Augment along the found path (as much flow as possible)
1.362 Value delta;
1.363 - int u, v = path_nodes.front(), pa;
1.364 - for (int i = 1; i < int(path_nodes.size()); ++i) {
1.365 + int pa, u, v = start;
1.366 + for (int i = 0; i != int(path.size()); ++i) {
1.367 + pa = path[i];
1.368 u = v;
1.369 - v = path_nodes[i];
1.370 - pa = pred_arc[v];
1.371 + v = _target[pa];
1.372 delta = std::min(_res_cap[pa], _excess[u]);
1.373 _res_cap[pa] -= delta;
1.374 _res_cap[_reverse[pa]] += delta;
1.375 @@ -1041,6 +1171,12 @@
1.376 if (_excess[v] > 0 && _excess[v] <= delta)
1.377 _active_nodes.push_back(v);
1.378 }
1.379 +
1.380 + // Global update heuristic
1.381 + if (relabel_cnt >= next_update_limit) {
1.382 + globalUpdate();
1.383 + next_update_limit += global_update_freq;
1.384 + }
1.385 }
1.386 }
1.387 }
1.388 @@ -1048,98 +1184,70 @@
1.389 /// Execute the algorithm performing push and relabel operations
1.390 void startPush() {
1.391 // Paramters for heuristics
1.392 - const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.393 - const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.394 + const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.395 + const double GLOBAL_UPDATE_FACTOR = 2.0;
1.396
1.397 + const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.398 + (_res_node_num + _sup_node_num * _sup_node_num));
1.399 + int next_update_limit = global_update_freq;
1.400 +
1.401 + int relabel_cnt = 0;
1.402 +
1.403 // Perform cost scaling phases
1.404 BoolVector hyper(_res_node_num, false);
1.405 + LargeCostVector hyper_cost(_res_node_num);
1.406 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.407 1 : _epsilon / _alpha )
1.408 {
1.409 - // "Early Termination" heuristic: use Bellman-Ford algorithm
1.410 - // to check if the current flow is optimal
1.411 - if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.412 - _arc_vec.clear();
1.413 - _cost_vec.clear();
1.414 - for (int j = 0; j != _res_arc_num; ++j) {
1.415 - if (_res_cap[j] > 0) {
1.416 - _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.417 - _cost_vec.push_back(_cost[j] + 1);
1.418 - }
1.419 - }
1.420 - _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.421 -
1.422 - BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.423 - bf.init(0);
1.424 - bool done = false;
1.425 - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1.426 - for (int i = 0; i < K && !done; ++i)
1.427 - done = bf.processNextWeakRound();
1.428 - if (done) break;
1.429 + // Early termination heuristic
1.430 + if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.431 + if (earlyTermination()) break;
1.432 }
1.433 -
1.434 - // Saturate arcs not satisfying the optimality condition
1.435 - for (int a = 0; a != _res_arc_num; ++a) {
1.436 - if (_res_cap[a] > 0 &&
1.437 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.438 - Value delta = _res_cap[a];
1.439 - _excess[_source[a]] -= delta;
1.440 - _excess[_target[a]] += delta;
1.441 - _res_cap[a] = 0;
1.442 - _res_cap[_reverse[a]] += delta;
1.443 - }
1.444 - }
1.445 -
1.446 - // Find active nodes (i.e. nodes with positive excess)
1.447 - for (int u = 0; u != _res_node_num; ++u) {
1.448 - if (_excess[u] > 0) _active_nodes.push_back(u);
1.449 - }
1.450 -
1.451 - // Initialize the next arcs
1.452 - for (int u = 0; u != _res_node_num; ++u) {
1.453 - _next_out[u] = _first_out[u];
1.454 - }
1.455 +
1.456 + // Initialize current phase
1.457 + initPhase();
1.458
1.459 // Perform push and relabel operations
1.460 while (_active_nodes.size() > 0) {
1.461 - LargeCost min_red_cost, rc;
1.462 + LargeCost min_red_cost, rc, pi_n;
1.463 Value delta;
1.464 int n, t, a, last_out = _res_arc_num;
1.465
1.466 + next_node:
1.467 // Select an active node (FIFO selection)
1.468 - next_node:
1.469 n = _active_nodes.front();
1.470 - last_out = _sum_supply < 0 ?
1.471 - _first_out[n+1] : _first_out[n+1] - 1;
1.472 -
1.473 + last_out = _first_out[n+1];
1.474 + pi_n = _pi[n];
1.475 +
1.476 // Perform push operations if there are admissible arcs
1.477 if (_excess[n] > 0) {
1.478 for (a = _next_out[n]; a != last_out; ++a) {
1.479 if (_res_cap[a] > 0 &&
1.480 - _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.481 + _cost[a] + pi_n - _pi[_target[a]] < 0) {
1.482 delta = std::min(_res_cap[a], _excess[n]);
1.483 t = _target[a];
1.484
1.485 // Push-look-ahead heuristic
1.486 Value ahead = -_excess[t];
1.487 - int last_out_t = _sum_supply < 0 ?
1.488 - _first_out[t+1] : _first_out[t+1] - 1;
1.489 + int last_out_t = _first_out[t+1];
1.490 + LargeCost pi_t = _pi[t];
1.491 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1.492 if (_res_cap[ta] > 0 &&
1.493 - _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
1.494 + _cost[ta] + pi_t - _pi[_target[ta]] < 0)
1.495 ahead += _res_cap[ta];
1.496 if (ahead >= delta) break;
1.497 }
1.498 if (ahead < 0) ahead = 0;
1.499
1.500 // Push flow along the arc
1.501 - if (ahead < delta) {
1.502 + if (ahead < delta && !hyper[t]) {
1.503 _res_cap[a] -= ahead;
1.504 _res_cap[_reverse[a]] += ahead;
1.505 _excess[n] -= ahead;
1.506 _excess[t] += ahead;
1.507 _active_nodes.push_front(t);
1.508 hyper[t] = true;
1.509 + hyper_cost[t] = _cost[a] + pi_n - pi_t;
1.510 _next_out[n] = a;
1.511 goto next_node;
1.512 } else {
1.513 @@ -1162,18 +1270,18 @@
1.514
1.515 // Relabel the node if it is still active (or hyper)
1.516 if (_excess[n] > 0 || hyper[n]) {
1.517 - min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1.518 + min_red_cost = hyper[n] ? -hyper_cost[n] :
1.519 + std::numeric_limits<LargeCost>::max();
1.520 for (int a = _first_out[n]; a != last_out; ++a) {
1.521 - rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1.522 + rc = _cost[a] + pi_n - _pi[_target[a]];
1.523 if (_res_cap[a] > 0 && rc < min_red_cost) {
1.524 min_red_cost = rc;
1.525 }
1.526 }
1.527 _pi[n] -= min_red_cost + _epsilon;
1.528 + _next_out[n] = _first_out[n];
1.529 hyper[n] = false;
1.530 -
1.531 - // Reset the next arc
1.532 - _next_out[n] = _first_out[n];
1.533 + ++relabel_cnt;
1.534 }
1.535
1.536 // Remove nodes that are not active nor hyper
1.537 @@ -1183,6 +1291,14 @@
1.538 !hyper[_active_nodes.front()] ) {
1.539 _active_nodes.pop_front();
1.540 }
1.541 +
1.542 + // Global update heuristic
1.543 + if (relabel_cnt >= next_update_limit) {
1.544 + globalUpdate();
1.545 + for (int u = 0; u != _res_node_num; ++u)
1.546 + hyper[u] = false;
1.547 + next_update_limit += global_update_freq;
1.548 + }
1.549 }
1.550 }
1.551 }