1.1 --- a/lemon/cost_scaling.h Tue Mar 15 17:59:57 2011 +0100
1.2 +++ b/lemon/cost_scaling.h Tue Mar 15 19:16:20 2011 +0100
1.3 @@ -1103,16 +1103,15 @@
1.4 void startAugment(int max_length) {
1.5 // Paramters for heuristics
1.6 const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.7 - const double GLOBAL_UPDATE_FACTOR = 3.0;
1.8 -
1.9 - const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.10 + const double GLOBAL_UPDATE_FACTOR = 1.0;
1.11 + const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1.12 (_res_node_num + _sup_node_num * _sup_node_num));
1.13 - int next_update_limit = global_update_freq;
1.14 -
1.15 - int relabel_cnt = 0;
1.16 + int next_global_update_limit = global_update_skip;
1.17
1.18 // Perform cost scaling phases
1.19 - std::vector<int> path;
1.20 + IntVector path;
1.21 + BoolVector path_arc(_res_arc_num, false);
1.22 + int relabel_cnt = 0;
1.23 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.24 1 : _epsilon / _alpha )
1.25 {
1.26 @@ -1135,32 +1134,45 @@
1.27 int start = _active_nodes.front();
1.28
1.29 // Find an augmenting path from the start node
1.30 - path.clear();
1.31 int tip = start;
1.32 - while (_excess[tip] >= 0 && int(path.size()) < max_length) {
1.33 + while (int(path.size()) < max_length && _excess[tip] >= 0) {
1.34 int u;
1.35 - LargeCost min_red_cost, rc, pi_tip = _pi[tip];
1.36 + LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
1.37 + LargeCost pi_tip = _pi[tip];
1.38 int last_out = _first_out[tip+1];
1.39 for (int a = _next_out[tip]; a != last_out; ++a) {
1.40 - u = _target[a];
1.41 - if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
1.42 - path.push_back(a);
1.43 - _next_out[tip] = a;
1.44 - tip = u;
1.45 - goto next_step;
1.46 + if (_res_cap[a] > 0) {
1.47 + u = _target[a];
1.48 + rc = _cost[a] + pi_tip - _pi[u];
1.49 + if (rc < 0) {
1.50 + path.push_back(a);
1.51 + _next_out[tip] = a;
1.52 + if (path_arc[a]) {
1.53 + goto augment; // a cycle is found, stop path search
1.54 + }
1.55 + tip = u;
1.56 + path_arc[a] = true;
1.57 + goto next_step;
1.58 + }
1.59 + else if (rc < min_red_cost) {
1.60 + min_red_cost = rc;
1.61 + }
1.62 }
1.63 }
1.64
1.65 // Relabel tip node
1.66 - min_red_cost = std::numeric_limits<LargeCost>::max();
1.67 if (tip != start) {
1.68 int ra = _reverse[path.back()];
1.69 - min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
1.70 + min_red_cost =
1.71 + std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
1.72 }
1.73 + last_out = _next_out[tip];
1.74 for (int a = _first_out[tip]; a != last_out; ++a) {
1.75 - rc = _cost[a] + pi_tip - _pi[_target[a]];
1.76 - if (_res_cap[a] > 0 && rc < min_red_cost) {
1.77 - min_red_cost = rc;
1.78 + if (_res_cap[a] > 0) {
1.79 + rc = _cost[a] + pi_tip - _pi[_target[a]];
1.80 + if (rc < min_red_cost) {
1.81 + min_red_cost = rc;
1.82 + }
1.83 }
1.84 }
1.85 _pi[tip] -= min_red_cost + _epsilon;
1.86 @@ -1169,7 +1181,9 @@
1.87
1.88 // Step back
1.89 if (tip != start) {
1.90 - tip = _source[path.back()];
1.91 + int pa = path.back();
1.92 + path_arc[pa] = false;
1.93 + tip = _source[pa];
1.94 path.pop_back();
1.95 }
1.96
1.97 @@ -1177,28 +1191,34 @@
1.98 }
1.99
1.100 // Augment along the found path (as much flow as possible)
1.101 + augment:
1.102 Value delta;
1.103 int pa, u, v = start;
1.104 for (int i = 0; i != int(path.size()); ++i) {
1.105 pa = path[i];
1.106 u = v;
1.107 v = _target[pa];
1.108 + path_arc[pa] = false;
1.109 delta = std::min(_res_cap[pa], _excess[u]);
1.110 _res_cap[pa] -= delta;
1.111 _res_cap[_reverse[pa]] += delta;
1.112 _excess[u] -= delta;
1.113 _excess[v] += delta;
1.114 - if (_excess[v] > 0 && _excess[v] <= delta)
1.115 + if (_excess[v] > 0 && _excess[v] <= delta) {
1.116 _active_nodes.push_back(v);
1.117 + }
1.118 }
1.119 + path.clear();
1.120
1.121 // Global update heuristic
1.122 - if (relabel_cnt >= next_update_limit) {
1.123 + if (relabel_cnt >= next_global_update_limit) {
1.124 globalUpdate();
1.125 - next_update_limit += global_update_freq;
1.126 + next_global_update_limit += global_update_skip;
1.127 }
1.128 }
1.129 +
1.130 }
1.131 +
1.132 }
1.133
1.134 /// Execute the algorithm performing push and relabel operations
1.135 @@ -1207,15 +1227,14 @@
1.136 const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.137 const double GLOBAL_UPDATE_FACTOR = 2.0;
1.138
1.139 - const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.140 + const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1.141 (_res_node_num + _sup_node_num * _sup_node_num));
1.142 - int next_update_limit = global_update_freq;
1.143 -
1.144 - int relabel_cnt = 0;
1.145 + int next_global_update_limit = global_update_skip;
1.146
1.147 // Perform cost scaling phases
1.148 BoolVector hyper(_res_node_num, false);
1.149 LargeCostVector hyper_cost(_res_node_num);
1.150 + int relabel_cnt = 0;
1.151 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.152 1 : _epsilon / _alpha )
1.153 {
1.154 @@ -1293,9 +1312,11 @@
1.155 min_red_cost = hyper[n] ? -hyper_cost[n] :
1.156 std::numeric_limits<LargeCost>::max();
1.157 for (int a = _first_out[n]; a != last_out; ++a) {
1.158 - rc = _cost[a] + pi_n - _pi[_target[a]];
1.159 - if (_res_cap[a] > 0 && rc < min_red_cost) {
1.160 - min_red_cost = rc;
1.161 + if (_res_cap[a] > 0) {
1.162 + rc = _cost[a] + pi_n - _pi[_target[a]];
1.163 + if (rc < min_red_cost) {
1.164 + min_red_cost = rc;
1.165 + }
1.166 }
1.167 }
1.168 _pi[n] -= min_red_cost + _epsilon;
1.169 @@ -1313,11 +1334,11 @@
1.170 }
1.171
1.172 // Global update heuristic
1.173 - if (relabel_cnt >= next_update_limit) {
1.174 + if (relabel_cnt >= next_global_update_limit) {
1.175 globalUpdate();
1.176 for (int u = 0; u != _res_node_num; ++u)
1.177 hyper[u] = false;
1.178 - next_update_limit += global_update_freq;
1.179 + next_global_update_limit += global_update_skip;
1.180 }
1.181 }
1.182 }