Changes in lemon/cost_scaling.h [1049:a07b6b27fe69:1041:f112c18bc304] in lemon
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/cost_scaling.h
r1049 r1041 98 98 /// "preflow push-relabel" algorithm for the maximum flow problem. 99 99 /// 100 /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest101 /// implementations available in LEMON for this problem.102 ///103 100 /// Most of the parameters of the problem (except for the digraph) 104 101 /// can be given using separate functions, and the algorithm can be … … 117 114 /// consider to use the named template parameters instead. 118 115 /// 119 /// \warning Both \c V and \c C must be signed number types. 120 /// \warning All input data (capacities, supply values, and costs) must 116 /// \warning Both number types must be signed and all input data must 121 117 /// be integer. 122 /// \warning This algorithm does not support negative costs for 123 /// arcs havinginfinite upper bound.118 /// \warning This algorithm does not support negative costs for such 119 /// arcs that have infinite upper bound. 124 120 /// 125 121 /// \note %CostScaling provides three different internal methods, … … 183 179 /// relabel operation. 184 180 /// By default, the so called \ref PARTIAL_AUGMENT 185 /// "Partial Augment-Relabel" method is used, which turned outto be181 /// "Partial Augment-Relabel" method is used, which proved to be 186 182 /// the most efficient and the most robust on various test inputs. 187 183 /// However, the other methods can be selected using the \ref run() … … 238 234 }; 239 235 236 typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap; 240 237 typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap; 241 238 … … 288 285 int _max_rank; 289 286 287 // Data for a StaticDigraph structure 288 typedef std::pair<int, int> IntPair; 289 StaticDigraph _sgr; 290 std::vector<IntPair> _arc_vec; 291 std::vector<LargeCost> _cost_vec; 292 LargeCostArcMap _cost_map; 293 LargeCostNodeMap _pi_map; 294 290 295 public: 291 296 … … 334 339 CostScaling(const GR& graph) : 335 340 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), 341 _cost_map(_cost_vec), _pi_map(_pi), 336 342 INF(std::numeric_limits<Value>::has_infinity ? 337 343 std::numeric_limits<Value>::infinity() : … … 442 448 /// 443 449 /// Using this function has the same effect as using \ref supplyMap() 444 /// with a map in which \c k is assigned to \c s, \c -k is450 /// with such a map in which \c k is assigned to \c s, \c -k is 445 451 /// assigned to \c t and all other nodes have zero supply value. 446 452 /// … … 488 494 /// \param method The internal method that will be used in the 489 495 /// algorithm. For more information, see \ref Method. 490 /// \param factor The cost scaling factor. It must be at least two.496 /// \param factor The cost scaling factor. It must be larger than one. 491 497 /// 492 498 /// \return \c INFEASIBLE if no feasible flow exists, … … 502 508 /// \see ProblemType, Method 503 509 /// \see resetParams(), reset() 504 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) { 505 LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2"); 510 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) { 506 511 _alpha = factor; 507 512 ProblemType pt = init(); … … 567 572 } 568 573 569 /// \brief Reset the internal data structures and all the parameters 570 /// that have been given before. 571 /// 572 /// This function resets the internal data structures and all the 573 /// paramaters that have been given before using functions \ref lowerMap(), 574 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). 575 /// 576 /// It is useful for multiple \ref run() calls. By default, all the given 577 /// parameters are kept for the next \ref run() call, unless 578 /// \ref resetParams() or \ref reset() is used. 579 /// If the underlying digraph was also modified after the construction 580 /// of the class or the last \ref reset() call, then the \ref reset() 581 /// function must be used, otherwise \ref resetParams() is sufficient. 582 /// 583 /// See \ref resetParams() for examples. 584 /// 574 /// \brief Reset all the parameters that have been given before. 575 /// 576 /// This function resets all the paramaters that have been given 577 /// before using functions \ref lowerMap(), \ref upperMap(), 578 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). 579 /// 580 /// It is useful for multiple run() calls. If this function is not 581 /// used, all the parameters given before are kept for the next 582 /// \ref run() call. 583 /// However, the underlying digraph must not be modified after this 584 /// class have been constructed, since it copies and extends the graph. 585 585 /// \return <tt>(*this)</tt> 586 ///587 /// \see resetParams(), run()588 586 CostScaling& reset() { 589 587 // Resize vectors … … 610 608 _excess.resize(_res_node_num); 611 609 _next_out.resize(_res_node_num); 610 611 _arc_vec.reserve(_res_arc_num); 612 _cost_vec.reserve(_res_arc_num); 612 613 613 614 // Copy the graph … … 886 887 } 887 888 889 return OPTIMAL; 890 } 891 892 // Execute the algorithm and transform the results 893 void start(Method method) { 894 // Maximum path length for partial augment 895 const int MAX_PATH_LENGTH = 4; 896 888 897 // Initialize data structures for buckets 889 898 _max_rank = _alpha * _res_node_num; … … 893 902 _rank.resize(_res_node_num + 1); 894 903 895 return OPTIMAL; 896 } 897 898 // Execute the algorithm and transform the results 899 void start(Method method) { 900 const int MAX_PARTIAL_PATH_LENGTH = 4; 901 904 // Execute the algorithm 902 905 switch (method) { 903 906 case PUSH: … … 908 911 break; 909 912 case PARTIAL_AUGMENT: 910 startAugment(MAX_PA RTIAL_PATH_LENGTH);913 startAugment(MAX_PATH_LENGTH); 911 914 break; 912 915 } 913 916 914 // Compute node potentials (dual solution) 915 for (int i = 0; i != _res_node_num; ++i) { 916 _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha)); 917 } 918 bool optimal = true; 919 for (int i = 0; optimal && i != _res_node_num; ++i) { 920 LargeCost pi_i = _pi[i]; 921 int last_out = _first_out[i+1]; 922 for (int j = _first_out[i]; j != last_out; ++j) { 923 if (_res_cap[j] > 0 && _scost[j] + pi_i - _pi[_target[j]] < 0) { 924 optimal = false; 925 break; 926 } 927 } 928 } 929 930 if (!optimal) { 931 // Compute node potentials for the original costs with BellmanFord 932 // (if it is necessary) 933 typedef std::pair<int, int> IntPair; 934 StaticDigraph sgr; 935 std::vector<IntPair> arc_vec; 936 std::vector<LargeCost> cost_vec; 937 LargeCostArcMap cost_map(cost_vec); 938 939 arc_vec.clear(); 940 cost_vec.clear(); 941 for (int j = 0; j != _res_arc_num; ++j) { 942 if (_res_cap[j] > 0) { 943 int u = _source[j], v = _target[j]; 944 arc_vec.push_back(IntPair(u, v)); 945 cost_vec.push_back(_scost[j] + _pi[u] - _pi[v]); 946 } 947 } 948 sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end()); 949 950 typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create 951 bf(sgr, cost_map); 952 bf.init(0); 953 bf.start(); 954 955 for (int i = 0; i != _res_node_num; ++i) { 956 _pi[i] += bf.dist(sgr.node(i)); 957 } 958 } 959 960 // Shift potentials to meet the requirements of the GEQ type 961 // optimality conditions 962 LargeCost max_pot = _pi[_root]; 963 for (int i = 0; i != _res_node_num; ++i) { 964 if (_pi[i] > max_pot) max_pot = _pi[i]; 965 } 966 if (max_pot != 0) { 967 for (int i = 0; i != _res_node_num; ++i) { 968 _pi[i] -= max_pot; 969 } 970 } 917 // Compute node potentials for the original costs 918 _arc_vec.clear(); 919 _cost_vec.clear(); 920 for (int j = 0; j != _res_arc_num; ++j) { 921 if (_res_cap[j] > 0) { 922 _arc_vec.push_back(IntPair(_source[j], _target[j])); 923 _cost_vec.push_back(_scost[j]); 924 } 925 } 926 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 927 928 typename BellmanFord<StaticDigraph, LargeCostArcMap> 929 ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map); 930 bf.distMap(_pi_map); 931 bf.init(0); 932 bf.start(); 971 933 972 934 // Handle non-zero lower bounds … … 986 948 LargeCost pi_u = _pi[u]; 987 949 for (int a = _first_out[u]; a != last_out; ++a) { 988 Value delta = _res_cap[a]; 989 if (delta > 0) { 990 int v = _target[a]; 991 if (_cost[a] + pi_u - _pi[v] < 0) { 992 _excess[u] -= delta; 993 _excess[v] += delta; 994 _res_cap[a] = 0; 995 _res_cap[_reverse[a]] += delta; 996 } 950 int v = _target[a]; 951 if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) { 952 Value delta = _res_cap[a]; 953 _excess[u] -= delta; 954 _excess[v] += delta; 955 _res_cap[a] = 0; 956 _res_cap[_reverse[a]] += delta; 997 957 } 998 958 } … … 1010 970 } 1011 971 1012 // Price (potential) refinement heuristic 1013 bool priceRefinement() { 1014 1015 // Stack for stroing the topological order 1016 IntVector stack(_res_node_num); 1017 int stack_top; 1018 1019 // Perform phases 1020 while (topologicalSort(stack, stack_top)) { 1021 1022 // Compute node ranks in the acyclic admissible network and 1023 // store the nodes in buckets 1024 for (int i = 0; i != _res_node_num; ++i) { 1025 _rank[i] = 0; 1026 } 1027 const int bucket_end = _root + 1; 1028 for (int r = 0; r != _max_rank; ++r) { 1029 _buckets[r] = bucket_end; 1030 } 1031 int top_rank = 0; 1032 for ( ; stack_top >= 0; --stack_top) { 1033 int u = stack[stack_top], v; 1034 int rank_u = _rank[u]; 1035 1036 LargeCost rc, pi_u = _pi[u]; 1037 int last_out = _first_out[u+1]; 1038 for (int a = _first_out[u]; a != last_out; ++a) { 1039 if (_res_cap[a] > 0) { 1040 v = _target[a]; 1041 rc = _cost[a] + pi_u - _pi[v]; 1042 if (rc < 0) { 1043 LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon); 1044 if (nrc < LargeCost(_max_rank)) { 1045 int new_rank_v = rank_u + static_cast<int>(nrc); 1046 if (new_rank_v > _rank[v]) { 1047 _rank[v] = new_rank_v; 1048 } 1049 } 1050 } 1051 } 1052 } 1053 1054 if (rank_u > 0) { 1055 top_rank = std::max(top_rank, rank_u); 1056 int bfirst = _buckets[rank_u]; 1057 _bucket_next[u] = bfirst; 1058 _bucket_prev[bfirst] = u; 1059 _buckets[rank_u] = u; 1060 } 1061 } 1062 1063 // Check if the current flow is epsilon-optimal 1064 if (top_rank == 0) { 1065 return true; 1066 } 1067 1068 // Process buckets in top-down order 1069 for (int rank = top_rank; rank > 0; --rank) { 1070 while (_buckets[rank] != bucket_end) { 1071 // Remove the first node from the current bucket 1072 int u = _buckets[rank]; 1073 _buckets[rank] = _bucket_next[u]; 1074 1075 // Search the outgoing arcs of u 1076 LargeCost rc, pi_u = _pi[u]; 1077 int last_out = _first_out[u+1]; 1078 int v, old_rank_v, new_rank_v; 1079 for (int a = _first_out[u]; a != last_out; ++a) { 1080 if (_res_cap[a] > 0) { 1081 v = _target[a]; 1082 old_rank_v = _rank[v]; 1083 1084 if (old_rank_v < rank) { 1085 1086 // Compute the new rank of node v 1087 rc = _cost[a] + pi_u - _pi[v]; 1088 if (rc < 0) { 1089 new_rank_v = rank; 1090 } else { 1091 LargeCost nrc = rc / _epsilon; 1092 new_rank_v = 0; 1093 if (nrc < LargeCost(_max_rank)) { 1094 new_rank_v = rank - 1 - static_cast<int>(nrc); 1095 } 1096 } 1097 1098 // Change the rank of node v 1099 if (new_rank_v > old_rank_v) { 1100 _rank[v] = new_rank_v; 1101 1102 // Remove v from its old bucket 1103 if (old_rank_v > 0) { 1104 if (_buckets[old_rank_v] == v) { 1105 _buckets[old_rank_v] = _bucket_next[v]; 1106 } else { 1107 int pv = _bucket_prev[v], nv = _bucket_next[v]; 1108 _bucket_next[pv] = nv; 1109 _bucket_prev[nv] = pv; 1110 } 1111 } 1112 1113 // Insert v into its new bucket 1114 int nv = _buckets[new_rank_v]; 1115 _bucket_next[v] = nv; 1116 _bucket_prev[nv] = v; 1117 _buckets[new_rank_v] = v; 1118 } 1119 } 1120 } 1121 } 1122 1123 // Refine potential of node u 1124 _pi[u] -= rank * _epsilon; 1125 } 1126 } 1127 1128 } 1129 1130 return false; 1131 } 1132 1133 // Find and cancel cycles in the admissible network and 1134 // determine topological order using DFS 1135 bool topologicalSort(IntVector &stack, int &stack_top) { 1136 const int MAX_CYCLE_CANCEL = 1; 1137 1138 BoolVector reached(_res_node_num, false); 1139 BoolVector processed(_res_node_num, false); 1140 IntVector pred(_res_node_num); 1141 for (int i = 0; i != _res_node_num; ++i) { 1142 _next_out[i] = _first_out[i]; 1143 } 1144 stack_top = -1; 1145 1146 int cycle_cnt = 0; 1147 for (int start = 0; start != _res_node_num; ++start) { 1148 if (reached[start]) continue; 1149 1150 // Start DFS search from this start node 1151 pred[start] = -1; 1152 int tip = start, v; 1153 while (true) { 1154 // Check the outgoing arcs of the current tip node 1155 reached[tip] = true; 1156 LargeCost pi_tip = _pi[tip]; 1157 int a, last_out = _first_out[tip+1]; 1158 for (a = _next_out[tip]; a != last_out; ++a) { 1159 if (_res_cap[a] > 0) { 1160 v = _target[a]; 1161 if (_cost[a] + pi_tip - _pi[v] < 0) { 1162 if (!reached[v]) { 1163 // A new node is reached 1164 reached[v] = true; 1165 pred[v] = tip; 1166 _next_out[tip] = a; 1167 tip = v; 1168 a = _next_out[tip]; 1169 last_out = _first_out[tip+1]; 1170 break; 1171 } 1172 else if (!processed[v]) { 1173 // A cycle is found 1174 ++cycle_cnt; 1175 _next_out[tip] = a; 1176 1177 // Find the minimum residual capacity along the cycle 1178 Value d, delta = _res_cap[a]; 1179 int u, delta_node = tip; 1180 for (u = tip; u != v; ) { 1181 u = pred[u]; 1182 d = _res_cap[_next_out[u]]; 1183 if (d <= delta) { 1184 delta = d; 1185 delta_node = u; 1186 } 1187 } 1188 1189 // Augment along the cycle 1190 _res_cap[a] -= delta; 1191 _res_cap[_reverse[a]] += delta; 1192 for (u = tip; u != v; ) { 1193 u = pred[u]; 1194 int ca = _next_out[u]; 1195 _res_cap[ca] -= delta; 1196 _res_cap[_reverse[ca]] += delta; 1197 } 1198 1199 // Check the maximum number of cycle canceling 1200 if (cycle_cnt >= MAX_CYCLE_CANCEL) { 1201 return false; 1202 } 1203 1204 // Roll back search to delta_node 1205 if (delta_node != tip) { 1206 for (u = tip; u != delta_node; u = pred[u]) { 1207 reached[u] = false; 1208 } 1209 tip = delta_node; 1210 a = _next_out[tip] + 1; 1211 last_out = _first_out[tip+1]; 1212 break; 1213 } 1214 } 1215 } 1216 } 1217 } 1218 1219 // Step back to the previous node 1220 if (a == last_out) { 1221 processed[tip] = true; 1222 stack[++stack_top] = tip; 1223 tip = pred[tip]; 1224 if (tip < 0) { 1225 // Finish DFS from the current start node 1226 break; 1227 } 1228 ++_next_out[tip]; 1229 } 1230 } 1231 1232 } 1233 1234 return (cycle_cnt == 0); 972 // Early termination heuristic 973 bool earlyTermination() { 974 const double EARLY_TERM_FACTOR = 3.0; 975 976 // Build a static residual graph 977 _arc_vec.clear(); 978 _cost_vec.clear(); 979 for (int j = 0; j != _res_arc_num; ++j) { 980 if (_res_cap[j] > 0) { 981 _arc_vec.push_back(IntPair(_source[j], _target[j])); 982 _cost_vec.push_back(_cost[j] + 1); 983 } 984 } 985 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 986 987 // Run Bellman-Ford algorithm to check if the current flow is optimal 988 BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); 989 bf.init(0); 990 bool done = false; 991 int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num))); 992 for (int i = 0; i < K && !done; ++i) { 993 done = bf.processNextWeakRound(); 994 } 995 return done; 1235 996 } 1236 997 1237 998 // Global potential update heuristic 1238 999 void globalUpdate() { 1239 constint bucket_end = _root + 1;1000 int bucket_end = _root + 1; 1240 1001 1241 1002 // Initialize buckets … … 1244 1005 } 1245 1006 Value total_excess = 0; 1246 int b0 = bucket_end;1247 1007 for (int i = 0; i != _res_node_num; ++i) { 1248 1008 if (_excess[i] < 0) { 1249 1009 _rank[i] = 0; 1250 _bucket_next[i] = b0;1251 _bucket_prev[ b0] = i;1252 b0= i;1010 _bucket_next[i] = _buckets[0]; 1011 _bucket_prev[_buckets[0]] = i; 1012 _buckets[0] = i; 1253 1013 } else { 1254 1014 total_excess += _excess[i]; … … 1257 1017 } 1258 1018 if (total_excess == 0) return; 1259 _buckets[0] = b0;1260 1019 1261 1020 // Search the buckets … … 1279 1038 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon; 1280 1039 int new_rank_v = old_rank_v; 1281 if (nrc < LargeCost(_max_rank)) { 1282 new_rank_v = r + 1 + static_cast<int>(nrc); 1283 } 1040 if (nrc < LargeCost(_max_rank)) 1041 new_rank_v = r + 1 + int(nrc); 1284 1042 1285 1043 // Change the rank of v … … 1293 1051 _buckets[old_rank_v] = _bucket_next[v]; 1294 1052 } else { 1295 int pv = _bucket_prev[v], nv = _bucket_next[v]; 1296 _bucket_next[pv] = nv; 1297 _bucket_prev[nv] = pv; 1053 _bucket_next[_bucket_prev[v]] = _bucket_next[v]; 1054 _bucket_prev[_bucket_next[v]] = _bucket_prev[v]; 1298 1055 } 1299 1056 } 1300 1057 1301 // Insert v into its new bucket 1302 int nv = _buckets[new_rank_v]; 1303 _bucket_next[v] = nv; 1304 _bucket_prev[nv] = v; 1058 // Insert v to its new bucket 1059 _bucket_next[v] = _buckets[new_rank_v]; 1060 _bucket_prev[_buckets[new_rank_v]] = v; 1305 1061 _buckets[new_rank_v] = v; 1306 1062 } … … 1331 1087 void startAugment(int max_length) { 1332 1088 // Paramters for heuristics 1333 const int PRICE_REFINEMENT_LIMIT = 2; 1334 const double GLOBAL_UPDATE_FACTOR = 1.0; 1335 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * 1089 const int EARLY_TERM_EPSILON_LIMIT = 1000; 1090 const double GLOBAL_UPDATE_FACTOR = 3.0; 1091 1092 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * 1336 1093 (_res_node_num + _sup_node_num * _sup_node_num)); 1337 int next_global_update_limit = global_update_skip; 1094 int next_update_limit = global_update_freq; 1095 1096 int relabel_cnt = 0; 1338 1097 1339 1098 // Perform cost scaling phases 1340 IntVector path; 1341 BoolVector path_arc(_res_arc_num, false); 1342 int relabel_cnt = 0; 1343 int eps_phase_cnt = 0; 1099 std::vector<int> path; 1344 1100 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1345 1101 1 : _epsilon / _alpha ) 1346 1102 { 1347 ++eps_phase_cnt; 1348 1349 // Price refinement heuristic 1350 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1351 if (priceRefinement()) continue; 1103 // Early termination heuristic 1104 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1105 if (earlyTermination()) break; 1352 1106 } 1353 1107 … … 1366 1120 1367 1121 // Find an augmenting path from the start node 1122 path.clear(); 1368 1123 int tip = start; 1369 while ( int(path.size()) < max_length && _excess[tip] >= 0) {1124 while (_excess[tip] >= 0 && int(path.size()) < max_length) { 1370 1125 int u; 1371 LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max(); 1372 LargeCost pi_tip = _pi[tip]; 1126 LargeCost min_red_cost, rc, pi_tip = _pi[tip]; 1373 1127 int last_out = _first_out[tip+1]; 1374 1128 for (int a = _next_out[tip]; a != last_out; ++a) { 1375 if (_res_cap[a] > 0) { 1376 u = _target[a]; 1377 rc = _cost[a] + pi_tip - _pi[u]; 1378 if (rc < 0) { 1379 path.push_back(a); 1380 _next_out[tip] = a; 1381 if (path_arc[a]) { 1382 goto augment; // a cycle is found, stop path search 1383 } 1384 tip = u; 1385 path_arc[a] = true; 1386 goto next_step; 1387 } 1388 else if (rc < min_red_cost) { 1389 min_red_cost = rc; 1390 } 1129 u = _target[a]; 1130 if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) { 1131 path.push_back(a); 1132 _next_out[tip] = a; 1133 tip = u; 1134 goto next_step; 1391 1135 } 1392 1136 } 1393 1137 1394 1138 // Relabel tip node 1139 min_red_cost = std::numeric_limits<LargeCost>::max(); 1395 1140 if (tip != start) { 1396 1141 int ra = _reverse[path.back()]; 1397 min_red_cost = 1398 std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]); 1142 min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]]; 1399 1143 } 1400 last_out = _next_out[tip];1401 1144 for (int a = _first_out[tip]; a != last_out; ++a) { 1402 if (_res_cap[a] > 0) { 1403 rc = _cost[a] + pi_tip - _pi[_target[a]]; 1404 if (rc < min_red_cost) { 1405 min_red_cost = rc; 1406 } 1145 rc = _cost[a] + pi_tip - _pi[_target[a]]; 1146 if (_res_cap[a] > 0 && rc < min_red_cost) { 1147 min_red_cost = rc; 1407 1148 } 1408 1149 } … … 1413 1154 // Step back 1414 1155 if (tip != start) { 1415 int pa = path.back(); 1416 path_arc[pa] = false; 1417 tip = _source[pa]; 1156 tip = _source[path.back()]; 1418 1157 path.pop_back(); 1419 1158 } … … 1423 1162 1424 1163 // Augment along the found path (as much flow as possible) 1425 augment:1426 1164 Value delta; 1427 1165 int pa, u, v = start; … … 1430 1168 u = v; 1431 1169 v = _target[pa]; 1432 path_arc[pa] = false;1433 1170 delta = std::min(_res_cap[pa], _excess[u]); 1434 1171 _res_cap[pa] -= delta; … … 1436 1173 _excess[u] -= delta; 1437 1174 _excess[v] += delta; 1438 if (_excess[v] > 0 && _excess[v] <= delta) {1175 if (_excess[v] > 0 && _excess[v] <= delta) 1439 1176 _active_nodes.push_back(v); 1440 }1441 1177 } 1442 path.clear();1443 1178 1444 1179 // Global update heuristic 1445 if (relabel_cnt >= next_ global_update_limit) {1180 if (relabel_cnt >= next_update_limit) { 1446 1181 globalUpdate(); 1447 next_ global_update_limit += global_update_skip;1182 next_update_limit += global_update_freq; 1448 1183 } 1449 1184 } 1450 1451 } 1452 1185 } 1453 1186 } 1454 1187 … … 1456 1189 void startPush() { 1457 1190 // Paramters for heuristics 1458 const int PRICE_REFINEMENT_LIMIT = 2;1191 const int EARLY_TERM_EPSILON_LIMIT = 1000; 1459 1192 const double GLOBAL_UPDATE_FACTOR = 2.0; 1460 1193 1461 const int global_update_ skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *1194 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * 1462 1195 (_res_node_num + _sup_node_num * _sup_node_num)); 1463 int next_global_update_limit = global_update_skip; 1196 int next_update_limit = global_update_freq; 1197 1198 int relabel_cnt = 0; 1464 1199 1465 1200 // Perform cost scaling phases 1466 1201 BoolVector hyper(_res_node_num, false); 1467 1202 LargeCostVector hyper_cost(_res_node_num); 1468 int relabel_cnt = 0;1469 int eps_phase_cnt = 0;1470 1203 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1471 1204 1 : _epsilon / _alpha ) 1472 1205 { 1473 ++eps_phase_cnt; 1474 1475 // Price refinement heuristic 1476 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1477 if (priceRefinement()) continue; 1206 // Early termination heuristic 1207 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1208 if (earlyTermination()) break; 1478 1209 } 1479 1210 … … 1547 1278 std::numeric_limits<LargeCost>::max(); 1548 1279 for (int a = _first_out[n]; a != last_out; ++a) { 1549 if (_res_cap[a] > 0) { 1550 rc = _cost[a] + pi_n - _pi[_target[a]]; 1551 if (rc < min_red_cost) { 1552 min_red_cost = rc; 1553 } 1280 rc = _cost[a] + pi_n - _pi[_target[a]]; 1281 if (_res_cap[a] > 0 && rc < min_red_cost) { 1282 min_red_cost = rc; 1554 1283 } 1555 1284 } … … 1569 1298 1570 1299 // Global update heuristic 1571 if (relabel_cnt >= next_ global_update_limit) {1300 if (relabel_cnt >= next_update_limit) { 1572 1301 globalUpdate(); 1573 1302 for (int u = 0; u != _res_node_num; ++u) 1574 1303 hyper[u] = false; 1575 next_ global_update_limit += global_update_skip;1304 next_update_limit += global_update_freq; 1576 1305 } 1577 1306 }
Note: See TracChangeset
for help on using the changeset viewer.