1101 |
1101 |
1102 /// Execute the algorithm performing augment and relabel operations |
1102 /// Execute the algorithm performing augment and relabel operations |
1103 void startAugment(int max_length) { |
1103 void startAugment(int max_length) { |
1104 // Paramters for heuristics |
1104 // Paramters for heuristics |
1105 const int EARLY_TERM_EPSILON_LIMIT = 1000; |
1105 const int EARLY_TERM_EPSILON_LIMIT = 1000; |
1106 const double GLOBAL_UPDATE_FACTOR = 3.0; |
1106 const double GLOBAL_UPDATE_FACTOR = 1.0; |
1107 |
1107 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * |
1108 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * |
|
1109 (_res_node_num + _sup_node_num * _sup_node_num)); |
1108 (_res_node_num + _sup_node_num * _sup_node_num)); |
1110 int next_update_limit = global_update_freq; |
1109 int next_global_update_limit = global_update_skip; |
1111 |
1110 |
|
1111 // Perform cost scaling phases |
|
1112 IntVector path; |
|
1113 BoolVector path_arc(_res_arc_num, false); |
1112 int relabel_cnt = 0; |
1114 int relabel_cnt = 0; |
1113 |
|
1114 // Perform cost scaling phases |
|
1115 std::vector<int> path; |
|
1116 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
1115 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
1117 1 : _epsilon / _alpha ) |
1116 1 : _epsilon / _alpha ) |
1118 { |
1117 { |
1119 // Early termination heuristic |
1118 // Early termination heuristic |
1120 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { |
1119 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { |
1133 } |
1132 } |
1134 if (_active_nodes.size() == 0) break; |
1133 if (_active_nodes.size() == 0) break; |
1135 int start = _active_nodes.front(); |
1134 int start = _active_nodes.front(); |
1136 |
1135 |
1137 // Find an augmenting path from the start node |
1136 // Find an augmenting path from the start node |
1138 path.clear(); |
|
1139 int tip = start; |
1137 int tip = start; |
1140 while (_excess[tip] >= 0 && int(path.size()) < max_length) { |
1138 while (int(path.size()) < max_length && _excess[tip] >= 0) { |
1141 int u; |
1139 int u; |
1142 LargeCost min_red_cost, rc, pi_tip = _pi[tip]; |
1140 LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max(); |
|
1141 LargeCost pi_tip = _pi[tip]; |
1143 int last_out = _first_out[tip+1]; |
1142 int last_out = _first_out[tip+1]; |
1144 for (int a = _next_out[tip]; a != last_out; ++a) { |
1143 for (int a = _next_out[tip]; a != last_out; ++a) { |
1145 u = _target[a]; |
1144 if (_res_cap[a] > 0) { |
1146 if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) { |
1145 u = _target[a]; |
1147 path.push_back(a); |
1146 rc = _cost[a] + pi_tip - _pi[u]; |
1148 _next_out[tip] = a; |
1147 if (rc < 0) { |
1149 tip = u; |
1148 path.push_back(a); |
1150 goto next_step; |
1149 _next_out[tip] = a; |
|
1150 if (path_arc[a]) { |
|
1151 goto augment; // a cycle is found, stop path search |
|
1152 } |
|
1153 tip = u; |
|
1154 path_arc[a] = true; |
|
1155 goto next_step; |
|
1156 } |
|
1157 else if (rc < min_red_cost) { |
|
1158 min_red_cost = rc; |
|
1159 } |
1151 } |
1160 } |
1152 } |
1161 } |
1153 |
1162 |
1154 // Relabel tip node |
1163 // Relabel tip node |
1155 min_red_cost = std::numeric_limits<LargeCost>::max(); |
|
1156 if (tip != start) { |
1164 if (tip != start) { |
1157 int ra = _reverse[path.back()]; |
1165 int ra = _reverse[path.back()]; |
1158 min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]]; |
1166 min_red_cost = |
|
1167 std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]); |
1159 } |
1168 } |
|
1169 last_out = _next_out[tip]; |
1160 for (int a = _first_out[tip]; a != last_out; ++a) { |
1170 for (int a = _first_out[tip]; a != last_out; ++a) { |
1161 rc = _cost[a] + pi_tip - _pi[_target[a]]; |
1171 if (_res_cap[a] > 0) { |
1162 if (_res_cap[a] > 0 && rc < min_red_cost) { |
1172 rc = _cost[a] + pi_tip - _pi[_target[a]]; |
1163 min_red_cost = rc; |
1173 if (rc < min_red_cost) { |
|
1174 min_red_cost = rc; |
|
1175 } |
1164 } |
1176 } |
1165 } |
1177 } |
1166 _pi[tip] -= min_red_cost + _epsilon; |
1178 _pi[tip] -= min_red_cost + _epsilon; |
1167 _next_out[tip] = _first_out[tip]; |
1179 _next_out[tip] = _first_out[tip]; |
1168 ++relabel_cnt; |
1180 ++relabel_cnt; |
1169 |
1181 |
1170 // Step back |
1182 // Step back |
1171 if (tip != start) { |
1183 if (tip != start) { |
1172 tip = _source[path.back()]; |
1184 int pa = path.back(); |
|
1185 path_arc[pa] = false; |
|
1186 tip = _source[pa]; |
1173 path.pop_back(); |
1187 path.pop_back(); |
1174 } |
1188 } |
1175 |
1189 |
1176 next_step: ; |
1190 next_step: ; |
1177 } |
1191 } |
1178 |
1192 |
1179 // Augment along the found path (as much flow as possible) |
1193 // Augment along the found path (as much flow as possible) |
|
1194 augment: |
1180 Value delta; |
1195 Value delta; |
1181 int pa, u, v = start; |
1196 int pa, u, v = start; |
1182 for (int i = 0; i != int(path.size()); ++i) { |
1197 for (int i = 0; i != int(path.size()); ++i) { |
1183 pa = path[i]; |
1198 pa = path[i]; |
1184 u = v; |
1199 u = v; |
1185 v = _target[pa]; |
1200 v = _target[pa]; |
|
1201 path_arc[pa] = false; |
1186 delta = std::min(_res_cap[pa], _excess[u]); |
1202 delta = std::min(_res_cap[pa], _excess[u]); |
1187 _res_cap[pa] -= delta; |
1203 _res_cap[pa] -= delta; |
1188 _res_cap[_reverse[pa]] += delta; |
1204 _res_cap[_reverse[pa]] += delta; |
1189 _excess[u] -= delta; |
1205 _excess[u] -= delta; |
1190 _excess[v] += delta; |
1206 _excess[v] += delta; |
1191 if (_excess[v] > 0 && _excess[v] <= delta) |
1207 if (_excess[v] > 0 && _excess[v] <= delta) { |
1192 _active_nodes.push_back(v); |
1208 _active_nodes.push_back(v); |
|
1209 } |
1193 } |
1210 } |
|
1211 path.clear(); |
1194 |
1212 |
1195 // Global update heuristic |
1213 // Global update heuristic |
1196 if (relabel_cnt >= next_update_limit) { |
1214 if (relabel_cnt >= next_global_update_limit) { |
1197 globalUpdate(); |
1215 globalUpdate(); |
1198 next_update_limit += global_update_freq; |
1216 next_global_update_limit += global_update_skip; |
1199 } |
1217 } |
1200 } |
1218 } |
1201 } |
1219 |
|
1220 } |
|
1221 |
1202 } |
1222 } |
1203 |
1223 |
1204 /// Execute the algorithm performing push and relabel operations |
1224 /// Execute the algorithm performing push and relabel operations |
1205 void startPush() { |
1225 void startPush() { |
1206 // Paramters for heuristics |
1226 // Paramters for heuristics |
1207 const int EARLY_TERM_EPSILON_LIMIT = 1000; |
1227 const int EARLY_TERM_EPSILON_LIMIT = 1000; |
1208 const double GLOBAL_UPDATE_FACTOR = 2.0; |
1228 const double GLOBAL_UPDATE_FACTOR = 2.0; |
1209 |
1229 |
1210 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * |
1230 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * |
1211 (_res_node_num + _sup_node_num * _sup_node_num)); |
1231 (_res_node_num + _sup_node_num * _sup_node_num)); |
1212 int next_update_limit = global_update_freq; |
1232 int next_global_update_limit = global_update_skip; |
1213 |
|
1214 int relabel_cnt = 0; |
|
1215 |
1233 |
1216 // Perform cost scaling phases |
1234 // Perform cost scaling phases |
1217 BoolVector hyper(_res_node_num, false); |
1235 BoolVector hyper(_res_node_num, false); |
1218 LargeCostVector hyper_cost(_res_node_num); |
1236 LargeCostVector hyper_cost(_res_node_num); |
|
1237 int relabel_cnt = 0; |
1219 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
1238 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
1220 1 : _epsilon / _alpha ) |
1239 1 : _epsilon / _alpha ) |
1221 { |
1240 { |
1222 // Early termination heuristic |
1241 // Early termination heuristic |
1223 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { |
1242 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { |