1.1 --- a/lemon/cost_scaling.h Tue Mar 15 19:16:20 2011 +0100
1.2 +++ b/lemon/cost_scaling.h Tue Mar 15 19:32:21 2011 +0100
1.3 @@ -980,30 +980,229 @@
1.4 }
1.5 }
1.6
1.7 - // Early termination heuristic
1.8 - bool earlyTermination() {
1.9 - const double EARLY_TERM_FACTOR = 3.0;
1.10 + // Price (potential) refinement heuristic
1.11 + bool priceRefinement() {
1.12
1.13 - // Build a static residual graph
1.14 - _arc_vec.clear();
1.15 - _cost_vec.clear();
1.16 - for (int j = 0; j != _res_arc_num; ++j) {
1.17 - if (_res_cap[j] > 0) {
1.18 - _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.19 - _cost_vec.push_back(_cost[j] + 1);
1.20 + // Stack for stroing the topological order
1.21 + IntVector stack(_res_node_num);
1.22 + int stack_top;
1.23 +
1.24 + // Perform phases
1.25 + while (topologicalSort(stack, stack_top)) {
1.26 +
1.27 + // Compute node ranks in the acyclic admissible network and
1.28 + // store the nodes in buckets
1.29 + for (int i = 0; i != _res_node_num; ++i) {
1.30 + _rank[i] = 0;
1.31 }
1.32 + const int bucket_end = _root + 1;
1.33 + for (int r = 0; r != _max_rank; ++r) {
1.34 + _buckets[r] = bucket_end;
1.35 + }
1.36 + int top_rank = 0;
1.37 + for ( ; stack_top >= 0; --stack_top) {
1.38 + int u = stack[stack_top], v;
1.39 + int rank_u = _rank[u];
1.40 +
1.41 + LargeCost rc, pi_u = _pi[u];
1.42 + int last_out = _first_out[u+1];
1.43 + for (int a = _first_out[u]; a != last_out; ++a) {
1.44 + if (_res_cap[a] > 0) {
1.45 + v = _target[a];
1.46 + rc = _cost[a] + pi_u - _pi[v];
1.47 + if (rc < 0) {
1.48 + LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
1.49 + if (nrc < LargeCost(_max_rank)) {
1.50 + int new_rank_v = rank_u + static_cast<int>(nrc);
1.51 + if (new_rank_v > _rank[v]) {
1.52 + _rank[v] = new_rank_v;
1.53 + }
1.54 + }
1.55 + }
1.56 + }
1.57 + }
1.58 +
1.59 + if (rank_u > 0) {
1.60 + top_rank = std::max(top_rank, rank_u);
1.61 + int bfirst = _buckets[rank_u];
1.62 + _bucket_next[u] = bfirst;
1.63 + _bucket_prev[bfirst] = u;
1.64 + _buckets[rank_u] = u;
1.65 + }
1.66 + }
1.67 +
1.68 + // Check if the current flow is epsilon-optimal
1.69 + if (top_rank == 0) {
1.70 + return true;
1.71 + }
1.72 +
1.73 + // Process buckets in top-down order
1.74 + for (int rank = top_rank; rank > 0; --rank) {
1.75 + while (_buckets[rank] != bucket_end) {
1.76 + // Remove the first node from the current bucket
1.77 + int u = _buckets[rank];
1.78 + _buckets[rank] = _bucket_next[u];
1.79 +
1.80 + // Search the outgoing arcs of u
1.81 + LargeCost rc, pi_u = _pi[u];
1.82 + int last_out = _first_out[u+1];
1.83 + int v, old_rank_v, new_rank_v;
1.84 + for (int a = _first_out[u]; a != last_out; ++a) {
1.85 + if (_res_cap[a] > 0) {
1.86 + v = _target[a];
1.87 + old_rank_v = _rank[v];
1.88 +
1.89 + if (old_rank_v < rank) {
1.90 +
1.91 + // Compute the new rank of node v
1.92 + rc = _cost[a] + pi_u - _pi[v];
1.93 + if (rc < 0) {
1.94 + new_rank_v = rank;
1.95 + } else {
1.96 + LargeCost nrc = rc / _epsilon;
1.97 + new_rank_v = 0;
1.98 + if (nrc < LargeCost(_max_rank)) {
1.99 + new_rank_v = rank - 1 - static_cast<int>(nrc);
1.100 + }
1.101 + }
1.102 +
1.103 + // Change the rank of node v
1.104 + if (new_rank_v > old_rank_v) {
1.105 + _rank[v] = new_rank_v;
1.106 +
1.107 + // Remove v from its old bucket
1.108 + if (old_rank_v > 0) {
1.109 + if (_buckets[old_rank_v] == v) {
1.110 + _buckets[old_rank_v] = _bucket_next[v];
1.111 + } else {
1.112 + int pv = _bucket_prev[v], nv = _bucket_next[v];
1.113 + _bucket_next[pv] = nv;
1.114 + _bucket_prev[nv] = pv;
1.115 + }
1.116 + }
1.117 +
1.118 + // Insert v into its new bucket
1.119 + int nv = _buckets[new_rank_v];
1.120 + _bucket_next[v] = nv;
1.121 + _bucket_prev[nv] = v;
1.122 + _buckets[new_rank_v] = v;
1.123 + }
1.124 + }
1.125 + }
1.126 + }
1.127 +
1.128 + // Refine potential of node u
1.129 + _pi[u] -= rank * _epsilon;
1.130 + }
1.131 + }
1.132 +
1.133 }
1.134 - _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.135
1.136 - // Run Bellman-Ford algorithm to check if the current flow is optimal
1.137 - BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.138 - bf.init(0);
1.139 - bool done = false;
1.140 - int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
1.141 - for (int i = 0; i < K && !done; ++i) {
1.142 - done = bf.processNextWeakRound();
1.143 + return false;
1.144 + }
1.145 +
1.146 + // Find and cancel cycles in the admissible network and
1.147 + // determine topological order using DFS
1.148 + bool topologicalSort(IntVector &stack, int &stack_top) {
1.149 + const int MAX_CYCLE_CANCEL = 1;
1.150 +
1.151 + BoolVector reached(_res_node_num, false);
1.152 + BoolVector processed(_res_node_num, false);
1.153 + IntVector pred(_res_node_num);
1.154 + for (int i = 0; i != _res_node_num; ++i) {
1.155 + _next_out[i] = _first_out[i];
1.156 }
1.157 - return done;
1.158 + stack_top = -1;
1.159 +
1.160 + int cycle_cnt = 0;
1.161 + for (int start = 0; start != _res_node_num; ++start) {
1.162 + if (reached[start]) continue;
1.163 +
1.164 + // Start DFS search from this start node
1.165 + pred[start] = -1;
1.166 + int tip = start, v;
1.167 + while (true) {
1.168 + // Check the outgoing arcs of the current tip node
1.169 + reached[tip] = true;
1.170 + LargeCost pi_tip = _pi[tip];
1.171 + int a, last_out = _first_out[tip+1];
1.172 + for (a = _next_out[tip]; a != last_out; ++a) {
1.173 + if (_res_cap[a] > 0) {
1.174 + v = _target[a];
1.175 + if (_cost[a] + pi_tip - _pi[v] < 0) {
1.176 + if (!reached[v]) {
1.177 + // A new node is reached
1.178 + reached[v] = true;
1.179 + pred[v] = tip;
1.180 + _next_out[tip] = a;
1.181 + tip = v;
1.182 + a = _next_out[tip];
1.183 + last_out = _first_out[tip+1];
1.184 + break;
1.185 + }
1.186 + else if (!processed[v]) {
1.187 + // A cycle is found
1.188 + ++cycle_cnt;
1.189 + _next_out[tip] = a;
1.190 +
1.191 + // Find the minimum residual capacity along the cycle
1.192 + Value d, delta = _res_cap[a];
1.193 + int u, delta_node = tip;
1.194 + for (u = tip; u != v; ) {
1.195 + u = pred[u];
1.196 + d = _res_cap[_next_out[u]];
1.197 + if (d <= delta) {
1.198 + delta = d;
1.199 + delta_node = u;
1.200 + }
1.201 + }
1.202 +
1.203 + // Augment along the cycle
1.204 + _res_cap[a] -= delta;
1.205 + _res_cap[_reverse[a]] += delta;
1.206 + for (u = tip; u != v; ) {
1.207 + u = pred[u];
1.208 + int ca = _next_out[u];
1.209 + _res_cap[ca] -= delta;
1.210 + _res_cap[_reverse[ca]] += delta;
1.211 + }
1.212 +
1.213 + // Check the maximum number of cycle canceling
1.214 + if (cycle_cnt >= MAX_CYCLE_CANCEL) {
1.215 + return false;
1.216 + }
1.217 +
1.218 + // Roll back search to delta_node
1.219 + if (delta_node != tip) {
1.220 + for (u = tip; u != delta_node; u = pred[u]) {
1.221 + reached[u] = false;
1.222 + }
1.223 + tip = delta_node;
1.224 + a = _next_out[tip] + 1;
1.225 + last_out = _first_out[tip+1];
1.226 + break;
1.227 + }
1.228 + }
1.229 + }
1.230 + }
1.231 + }
1.232 +
1.233 + // Step back to the previous node
1.234 + if (a == last_out) {
1.235 + processed[tip] = true;
1.236 + stack[++stack_top] = tip;
1.237 + tip = pred[tip];
1.238 + if (tip < 0) {
1.239 + // Finish DFS from the current start node
1.240 + break;
1.241 + }
1.242 + ++_next_out[tip];
1.243 + }
1.244 + }
1.245 +
1.246 + }
1.247 +
1.248 + return (cycle_cnt == 0);
1.249 }
1.250
1.251 // Global potential update heuristic
1.252 @@ -1102,7 +1301,7 @@
1.253 /// Execute the algorithm performing augment and relabel operations
1.254 void startAugment(int max_length) {
1.255 // Paramters for heuristics
1.256 - const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.257 + const int PRICE_REFINEMENT_LIMIT = 2;
1.258 const double GLOBAL_UPDATE_FACTOR = 1.0;
1.259 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1.260 (_res_node_num + _sup_node_num * _sup_node_num));
1.261 @@ -1112,12 +1311,15 @@
1.262 IntVector path;
1.263 BoolVector path_arc(_res_arc_num, false);
1.264 int relabel_cnt = 0;
1.265 + int eps_phase_cnt = 0;
1.266 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.267 1 : _epsilon / _alpha )
1.268 {
1.269 - // Early termination heuristic
1.270 - if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.271 - if (earlyTermination()) break;
1.272 + ++eps_phase_cnt;
1.273 +
1.274 + // Price refinement heuristic
1.275 + if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
1.276 + if (priceRefinement()) continue;
1.277 }
1.278
1.279 // Initialize current phase
1.280 @@ -1224,7 +1426,7 @@
1.281 /// Execute the algorithm performing push and relabel operations
1.282 void startPush() {
1.283 // Paramters for heuristics
1.284 - const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.285 + const int PRICE_REFINEMENT_LIMIT = 2;
1.286 const double GLOBAL_UPDATE_FACTOR = 2.0;
1.287
1.288 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1.289 @@ -1235,12 +1437,15 @@
1.290 BoolVector hyper(_res_node_num, false);
1.291 LargeCostVector hyper_cost(_res_node_num);
1.292 int relabel_cnt = 0;
1.293 + int eps_phase_cnt = 0;
1.294 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.295 1 : _epsilon / _alpha )
1.296 {
1.297 - // Early termination heuristic
1.298 - if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.299 - if (earlyTermination()) break;
1.300 + ++eps_phase_cnt;
1.301 +
1.302 + // Price refinement heuristic
1.303 + if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
1.304 + if (priceRefinement()) continue;
1.305 }
1.306
1.307 // Initialize current phase