Changes in lemon/cost_scaling.h [1041:f112c18bc304:1049:a07b6b27fe69] in lemon
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/cost_scaling.h
r1041 r1049 98 98 /// "preflow pushrelabel" algorithm for the maximum flow problem. 99 99 /// 100 /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest 101 /// implementations available in LEMON for this problem. 102 /// 100 103 /// Most of the parameters of the problem (except for the digraph) 101 104 /// can be given using separate functions, and the algorithm can be … … 114 117 /// consider to use the named template parameters instead. 115 118 /// 116 /// \warning Both number types must be signed and all input data must 119 /// \warning Both \c V and \c C must be signed number types. 120 /// \warning All input data (capacities, supply values, and costs) must 117 121 /// be integer. 118 /// \warning This algorithm does not support negative costs for such119 /// arcs that haveinfinite upper bound.122 /// \warning This algorithm does not support negative costs for 123 /// arcs having infinite upper bound. 120 124 /// 121 125 /// \note %CostScaling provides three different internal methods, … … 179 183 /// relabel operation. 180 184 /// By default, the so called \ref PARTIAL_AUGMENT 181 /// "Partial AugmentRelabel" method is used, which provedto be185 /// "Partial AugmentRelabel" method is used, which turned out to be 182 186 /// the most efficient and the most robust on various test inputs. 183 187 /// However, the other methods can be selected using the \ref run() … … 234 238 }; 235 239 236 typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;237 240 typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap; 238 241 … … 285 288 int _max_rank; 286 289 287 // Data for a StaticDigraph structure288 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 295 290 public: 296 291 … … 339 334 CostScaling(const GR& graph) : 340 335 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), 341 _cost_map(_cost_vec), _pi_map(_pi),342 336 INF(std::numeric_limits<Value>::has_infinity ? 343 337 std::numeric_limits<Value>::infinity() : … … 448 442 /// 449 443 /// Using this function has the same effect as using \ref supplyMap() 450 /// with sucha map in which \c k is assigned to \c s, \c k is444 /// with a map in which \c k is assigned to \c s, \c k is 451 445 /// assigned to \c t and all other nodes have zero supply value. 452 446 /// … … 494 488 /// \param method The internal method that will be used in the 495 489 /// algorithm. For more information, see \ref Method. 496 /// \param factor The cost scaling factor. It must be larger than one.490 /// \param factor The cost scaling factor. It must be at least two. 497 491 /// 498 492 /// \return \c INFEASIBLE if no feasible flow exists, … … 508 502 /// \see ProblemType, Method 509 503 /// \see resetParams(), reset() 510 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) { 504 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) { 505 LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2"); 511 506 _alpha = factor; 512 507 ProblemType pt = init(); … … 572 567 } 573 568 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. 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 /// 585 585 /// \return <tt>(*this)</tt> 586 /// 587 /// \see resetParams(), run() 586 588 CostScaling& reset() { 587 589 // Resize vectors … … 608 610 _excess.resize(_res_node_num); 609 611 _next_out.resize(_res_node_num); 610 611 _arc_vec.reserve(_res_arc_num);612 _cost_vec.reserve(_res_arc_num);613 612 614 613 // Copy the graph … … 887 886 } 888 887 889 return OPTIMAL;890 }891 892 // Execute the algorithm and transform the results893 void start(Method method) {894 // Maximum path length for partial augment895 const int MAX_PATH_LENGTH = 4;896 897 888 // Initialize data structures for buckets 898 889 _max_rank = _alpha * _res_node_num; … … 902 893 _rank.resize(_res_node_num + 1); 903 894 904 // Execute the algorithm 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 905 902 switch (method) { 906 903 case PUSH: … … 911 908 break; 912 909 case PARTIAL_AUGMENT: 913 startAugment(MAX_PA TH_LENGTH);910 startAugment(MAX_PARTIAL_PATH_LENGTH); 914 911 break; 915 912 } 916 913 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(); 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 } 933 971 934 972 // Handle nonzero lower bounds … … 948 986 LargeCost pi_u = _pi[u]; 949 987 for (int a = _first_out[u]; a != last_out; ++a) { 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; 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 } 957 997 } 958 998 } … … 970 1010 } 971 1011 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 BellmanFord 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; 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 epsilonoptimal 1064 if (top_rank == 0) { 1065 return true; 1066 } 1067 1068 // Process buckets in topdown 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); 996 1235 } 997 1236 998 1237 // Global potential update heuristic 999 1238 void globalUpdate() { 1000 int bucket_end = _root + 1;1239 const int bucket_end = _root + 1; 1001 1240 1002 1241 // Initialize buckets … … 1005 1244 } 1006 1245 Value total_excess = 0; 1246 int b0 = bucket_end; 1007 1247 for (int i = 0; i != _res_node_num; ++i) { 1008 1248 if (_excess[i] < 0) { 1009 1249 _rank[i] = 0; 1010 _bucket_next[i] = _buckets[0];1011 _bucket_prev[ _buckets[0]] = i;1012 _buckets[0]= i;1250 _bucket_next[i] = b0; 1251 _bucket_prev[b0] = i; 1252 b0 = i; 1013 1253 } else { 1014 1254 total_excess += _excess[i]; … … 1017 1257 } 1018 1258 if (total_excess == 0) return; 1259 _buckets[0] = b0; 1019 1260 1020 1261 // Search the buckets … … 1038 1279 LargeCost nrc = (_cost[ra] + _pi[v]  pi_u) / _epsilon; 1039 1280 int new_rank_v = old_rank_v; 1040 if (nrc < LargeCost(_max_rank)) 1041 new_rank_v = r + 1 + int(nrc); 1281 if (nrc < LargeCost(_max_rank)) { 1282 new_rank_v = r + 1 + static_cast<int>(nrc); 1283 } 1042 1284 1043 1285 // Change the rank of v … … 1051 1293 _buckets[old_rank_v] = _bucket_next[v]; 1052 1294 } else { 1053 _bucket_next[_bucket_prev[v]] = _bucket_next[v]; 1054 _bucket_prev[_bucket_next[v]] = _bucket_prev[v]; 1295 int pv = _bucket_prev[v], nv = _bucket_next[v]; 1296 _bucket_next[pv] = nv; 1297 _bucket_prev[nv] = pv; 1055 1298 } 1056 1299 } 1057 1300 1058 // Insert v to its new bucket 1059 _bucket_next[v] = _buckets[new_rank_v]; 1060 _bucket_prev[_buckets[new_rank_v]] = v; 1301 // Insert v into its new bucket 1302 int nv = _buckets[new_rank_v]; 1303 _bucket_next[v] = nv; 1304 _bucket_prev[nv] = v; 1061 1305 _buckets[new_rank_v] = v; 1062 1306 } … … 1087 1331 void startAugment(int max_length) { 1088 1332 // Paramters for heuristics 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 * 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 * 1093 1336 (_res_node_num + _sup_node_num * _sup_node_num)); 1094 int next_update_limit = global_update_freq; 1095 1337 int next_global_update_limit = global_update_skip; 1338 1339 // Perform cost scaling phases 1340 IntVector path; 1341 BoolVector path_arc(_res_arc_num, false); 1096 1342 int relabel_cnt = 0; 1097 1098 // Perform cost scaling phases 1099 std::vector<int> path; 1343 int eps_phase_cnt = 0; 1100 1344 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1101 1345 1 : _epsilon / _alpha ) 1102 1346 { 1103 // Early termination heuristic 1104 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1105 if (earlyTermination()) break; 1347 ++eps_phase_cnt; 1348 1349 // Price refinement heuristic 1350 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1351 if (priceRefinement()) continue; 1106 1352 } 1107 1353 … … 1120 1366 1121 1367 // Find an augmenting path from the start node 1122 path.clear();1123 1368 int tip = start; 1124 while ( _excess[tip] >= 0 && int(path.size()) < max_length) {1369 while (int(path.size()) < max_length && _excess[tip] >= 0) { 1125 1370 int u; 1126 LargeCost min_red_cost, rc, pi_tip = _pi[tip]; 1371 LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max(); 1372 LargeCost pi_tip = _pi[tip]; 1127 1373 int last_out = _first_out[tip+1]; 1128 1374 for (int a = _next_out[tip]; a != last_out; ++a) { 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; 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 } 1135 1391 } 1136 1392 } 1137 1393 1138 1394 // Relabel tip node 1139 min_red_cost = std::numeric_limits<LargeCost>::max();1140 1395 if (tip != start) { 1141 1396 int ra = _reverse[path.back()]; 1142 min_red_cost = _cost[ra] + pi_tip  _pi[_target[ra]]; 1397 min_red_cost = 1398 std::min(min_red_cost, _cost[ra] + pi_tip  _pi[_target[ra]]); 1143 1399 } 1400 last_out = _next_out[tip]; 1144 1401 for (int a = _first_out[tip]; a != last_out; ++a) { 1145 rc = _cost[a] + pi_tip  _pi[_target[a]]; 1146 if (_res_cap[a] > 0 && rc < min_red_cost) { 1147 min_red_cost = rc; 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 } 1148 1407 } 1149 1408 } … … 1154 1413 // Step back 1155 1414 if (tip != start) { 1156 tip = _source[path.back()]; 1415 int pa = path.back(); 1416 path_arc[pa] = false; 1417 tip = _source[pa]; 1157 1418 path.pop_back(); 1158 1419 } … … 1162 1423 1163 1424 // Augment along the found path (as much flow as possible) 1425 augment: 1164 1426 Value delta; 1165 1427 int pa, u, v = start; … … 1168 1430 u = v; 1169 1431 v = _target[pa]; 1432 path_arc[pa] = false; 1170 1433 delta = std::min(_res_cap[pa], _excess[u]); 1171 1434 _res_cap[pa] = delta; … … 1173 1436 _excess[u] = delta; 1174 1437 _excess[v] += delta; 1175 if (_excess[v] > 0 && _excess[v] <= delta) 1438 if (_excess[v] > 0 && _excess[v] <= delta) { 1176 1439 _active_nodes.push_back(v); 1177 } 1440 } 1441 } 1442 path.clear(); 1178 1443 1179 1444 // Global update heuristic 1180 if (relabel_cnt >= next_ update_limit) {1445 if (relabel_cnt >= next_global_update_limit) { 1181 1446 globalUpdate(); 1182 next_update_limit += global_update_freq; 1183 } 1184 } 1185 } 1447 next_global_update_limit += global_update_skip; 1448 } 1449 } 1450 1451 } 1452 1186 1453 } 1187 1454 … … 1189 1456 void startPush() { 1190 1457 // Paramters for heuristics 1191 const int EARLY_TERM_EPSILON_LIMIT = 1000;1458 const int PRICE_REFINEMENT_LIMIT = 2; 1192 1459 const double GLOBAL_UPDATE_FACTOR = 2.0; 1193 1460 1194 const int global_update_ freq = int(GLOBAL_UPDATE_FACTOR *1461 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * 1195 1462 (_res_node_num + _sup_node_num * _sup_node_num)); 1196 int next_update_limit = global_update_freq; 1197 1198 int relabel_cnt = 0; 1463 int next_global_update_limit = global_update_skip; 1199 1464 1200 1465 // Perform cost scaling phases 1201 1466 BoolVector hyper(_res_node_num, false); 1202 1467 LargeCostVector hyper_cost(_res_node_num); 1468 int relabel_cnt = 0; 1469 int eps_phase_cnt = 0; 1203 1470 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1204 1471 1 : _epsilon / _alpha ) 1205 1472 { 1206 // Early termination heuristic 1207 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1208 if (earlyTermination()) break; 1473 ++eps_phase_cnt; 1474 1475 // Price refinement heuristic 1476 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1477 if (priceRefinement()) continue; 1209 1478 } 1210 1479 … … 1278 1547 std::numeric_limits<LargeCost>::max(); 1279 1548 for (int a = _first_out[n]; a != last_out; ++a) { 1280 rc = _cost[a] + pi_n  _pi[_target[a]]; 1281 if (_res_cap[a] > 0 && rc < min_red_cost) { 1282 min_red_cost = rc; 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 } 1283 1554 } 1284 1555 } … … 1298 1569 1299 1570 // Global update heuristic 1300 if (relabel_cnt >= next_ update_limit) {1571 if (relabel_cnt >= next_global_update_limit) { 1301 1572 globalUpdate(); 1302 1573 for (int u = 0; u != _res_node_num; ++u) 1303 1574 hyper[u] = false; 1304 next_ update_limit += global_update_freq;1575 next_global_update_limit += global_update_skip; 1305 1576 } 1306 1577 }
Note: See TracChangeset
for help on using the changeset viewer.