Changes in lemon/cost_scaling.h [1041:f112c18bc304:1271:fb1c7da561ce] in lemon
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/cost_scaling.h
r1041 r1271 3 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 * Copyright (C) 2003201 05 * Copyright (C) 20032013 6 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). … … 92 92 /// \ref CostScaling implements a cost scaling algorithm that performs 93 93 /// push/augment and relabel operations for finding a \ref min_cost_flow 94 /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation, 95 /// \ref goldberg97efficient, \ref bunnagel98efficient. 94 /// "minimum cost flow" \cite amo93networkflows, 95 /// \cite goldberg90approximation, 96 /// \cite goldberg97efficient, \cite bunnagel98efficient. 96 97 /// It is a highly efficient primaldual solution method, which 97 98 /// can be viewed as the generalization of the \ref Preflow 98 99 /// "preflow pushrelabel" algorithm for the maximum flow problem. 100 /// It is a polynomial algorithm, its running time complexity is 101 /// \f$O(n^2m\log(nK))\f$, where <i>K</i> denotes the maximum arc cost. 102 /// 103 /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest 104 /// implementations available in LEMON for solving this problem. 105 /// (For more information, see \ref min_cost_flow_algs "the module page".) 99 106 /// 100 107 /// Most of the parameters of the problem (except for the digraph) … … 114 121 /// consider to use the named template parameters instead. 115 122 /// 116 /// \warning Both number types must be signed and all input data must 123 /// \warning Both \c V and \c C must be signed number types. 124 /// \warning All input data (capacities, supply values, and costs) must 117 125 /// be integer. 118 /// \warning This algorithm does not support negative costs for such119 /// arcs that haveinfinite upper bound.126 /// \warning This algorithm does not support negative costs for 127 /// arcs having infinite upper bound. 120 128 /// 121 129 /// \note %CostScaling provides three different internal methods, … … 146 154 typedef typename TR::LargeCost LargeCost; 147 155 148 /// The \ref CostScalingDefaultTraits "traits class" of the algorithm 156 /// \brief The \ref lemon::CostScalingDefaultTraits "traits class" 157 /// of the algorithm 149 158 typedef TR Traits; 150 159 … … 179 188 /// relabel operation. 180 189 /// By default, the so called \ref PARTIAL_AUGMENT 181 /// "Partial AugmentRelabel" method is used, which provedto be190 /// "Partial AugmentRelabel" method is used, which turned out to be 182 191 /// the most efficient and the most robust on various test inputs. 183 192 /// However, the other methods can be selected using the \ref run() … … 206 215 typedef std::vector<LargeCost> LargeCostVector; 207 216 typedef std::vector<char> BoolVector; 208 // Note: vector<char> is used instead of vector<bool> for efficiency reasons 217 // Note: vector<char> is used instead of vector<bool> 218 // for efficiency reasons 209 219 210 220 private: … … 234 244 }; 235 245 236 typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;237 246 typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap; 238 247 … … 285 294 int _max_rank; 286 295 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 296 public: 296 297 … … 339 340 CostScaling(const GR& graph) : 340 341 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), 341 _cost_map(_cost_vec), _pi_map(_pi),342 342 INF(std::numeric_limits<Value>::has_infinity ? 343 343 std::numeric_limits<Value>::infinity() : … … 448 448 /// 449 449 /// 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 is450 /// with a map in which \c k is assigned to \c s, \c k is 451 451 /// assigned to \c t and all other nodes have zero supply value. 452 452 /// … … 494 494 /// \param method The internal method that will be used in the 495 495 /// algorithm. For more information, see \ref Method. 496 /// \param factor The cost scaling factor. It must be larger than one.496 /// \param factor The cost scaling factor. It must be at least two. 497 497 /// 498 498 /// \return \c INFEASIBLE if no feasible flow exists, … … 508 508 /// \see ProblemType, Method 509 509 /// \see resetParams(), reset() 510 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) { 510 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) { 511 LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2"); 511 512 _alpha = factor; 512 513 ProblemType pt = init(); … … 572 573 } 573 574 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. 575 /// \brief Reset the internal data structures and all the parameters 576 /// that have been given before. 577 /// 578 /// This function resets the internal data structures and all the 579 /// paramaters that have been given before using functions \ref lowerMap(), 580 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). 581 /// 582 /// It is useful for multiple \ref run() calls. By default, all the given 583 /// parameters are kept for the next \ref run() call, unless 584 /// \ref resetParams() or \ref reset() is used. 585 /// If the underlying digraph was also modified after the construction 586 /// of the class or the last \ref reset() call, then the \ref reset() 587 /// function must be used, otherwise \ref resetParams() is sufficient. 588 /// 589 /// See \ref resetParams() for examples. 590 /// 585 591 /// \return <tt>(*this)</tt> 592 /// 593 /// \see resetParams(), run() 586 594 CostScaling& reset() { 587 595 // Resize vectors … … 608 616 _excess.resize(_res_node_num); 609 617 _next_out.resize(_res_node_num); 610 611 _arc_vec.reserve(_res_arc_num);612 _cost_vec.reserve(_res_arc_num);613 618 614 619 // Copy the graph … … 668 673 /// 669 674 /// This function returns the total cost of the found flow. 670 /// Its complexity is O( e).675 /// Its complexity is O(m). 671 676 /// 672 677 /// \note The return type of the function can be specified as a … … 706 711 } 707 712 708 /// \brief Return the flow map (the primal solution). 713 /// \brief Copy the flow values (the primal solution) into the 714 /// given map. 709 715 /// 710 716 /// This function copies the flow value on each arc into the given … … 730 736 } 731 737 732 /// \brief Return the potential map (the dual solution). 738 /// \brief Copy the potential values (the dual solution) into the 739 /// given map. 733 740 /// 734 741 /// This function copies the potential (dual value) of each node … … 759 766 } 760 767 if (_sum_supply > 0) return INFEASIBLE; 768 769 // Check lower and upper bounds 770 LEMON_DEBUG(checkBoundMaps(), 771 "Upper bounds must be greater or equal to the lower bounds"); 761 772 762 773 … … 887 898 } 888 899 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 900 // Initialize data structures for buckets 898 901 _max_rank = _alpha * _res_node_num; … … 902 905 _rank.resize(_res_node_num + 1); 903 906 904 // Execute the algorithm 907 return OPTIMAL; 908 } 909 910 // Check if the upper bound is greater or equal to the lower bound 911 // on each arc. 912 bool checkBoundMaps() { 913 for (int j = 0; j != _res_arc_num; ++j) { 914 if (_upper[j] < _lower[j]) return false; 915 } 916 return true; 917 } 918 919 // Execute the algorithm and transform the results 920 void start(Method method) { 921 const int MAX_PARTIAL_PATH_LENGTH = 4; 922 905 923 switch (method) { 906 924 case PUSH: … … 911 929 break; 912 930 case PARTIAL_AUGMENT: 913 startAugment(MAX_PA TH_LENGTH);931 startAugment(MAX_PARTIAL_PATH_LENGTH); 914 932 break; 915 933 } 916 934 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(); 935 // Compute node potentials (dual solution) 936 for (int i = 0; i != _res_node_num; ++i) { 937 _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha)); 938 } 939 bool optimal = true; 940 for (int i = 0; optimal && i != _res_node_num; ++i) { 941 LargeCost pi_i = _pi[i]; 942 int last_out = _first_out[i+1]; 943 for (int j = _first_out[i]; j != last_out; ++j) { 944 if (_res_cap[j] > 0 && _scost[j] + pi_i  _pi[_target[j]] < 0) { 945 optimal = false; 946 break; 947 } 948 } 949 } 950 951 if (!optimal) { 952 // Compute node potentials for the original costs with BellmanFord 953 // (if it is necessary) 954 typedef std::pair<int, int> IntPair; 955 StaticDigraph sgr; 956 std::vector<IntPair> arc_vec; 957 std::vector<LargeCost> cost_vec; 958 LargeCostArcMap cost_map(cost_vec); 959 960 arc_vec.clear(); 961 cost_vec.clear(); 962 for (int j = 0; j != _res_arc_num; ++j) { 963 if (_res_cap[j] > 0) { 964 int u = _source[j], v = _target[j]; 965 arc_vec.push_back(IntPair(u, v)); 966 cost_vec.push_back(_scost[j] + _pi[u]  _pi[v]); 967 } 968 } 969 sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end()); 970 971 typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create 972 bf(sgr, cost_map); 973 bf.init(0); 974 bf.start(); 975 976 for (int i = 0; i != _res_node_num; ++i) { 977 _pi[i] += bf.dist(sgr.node(i)); 978 } 979 } 980 981 // Shift potentials to meet the requirements of the GEQ type 982 // optimality conditions 983 LargeCost max_pot = _pi[_root]; 984 for (int i = 0; i != _res_node_num; ++i) { 985 if (_pi[i] > max_pot) max_pot = _pi[i]; 986 } 987 if (max_pot != 0) { 988 for (int i = 0; i != _res_node_num; ++i) { 989 _pi[i] = max_pot; 990 } 991 } 933 992 934 993 // Handle nonzero lower bounds … … 948 1007 LargeCost pi_u = _pi[u]; 949 1008 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; 1009 Value delta = _res_cap[a]; 1010 if (delta > 0) { 1011 int v = _target[a]; 1012 if (_cost[a] + pi_u  _pi[v] < 0) { 1013 _excess[u] = delta; 1014 _excess[v] += delta; 1015 _res_cap[a] = 0; 1016 _res_cap[_reverse[a]] += delta; 1017 } 957 1018 } 958 1019 } … … 970 1031 } 971 1032 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; 1033 // Price (potential) refinement heuristic 1034 bool priceRefinement() { 1035 1036 // Stack for stroing the topological order 1037 IntVector stack(_res_node_num); 1038 int stack_top; 1039 1040 // Perform phases 1041 while (topologicalSort(stack, stack_top)) { 1042 1043 // Compute node ranks in the acyclic admissible network and 1044 // store the nodes in buckets 1045 for (int i = 0; i != _res_node_num; ++i) { 1046 _rank[i] = 0; 1047 } 1048 const int bucket_end = _root + 1; 1049 for (int r = 0; r != _max_rank; ++r) { 1050 _buckets[r] = bucket_end; 1051 } 1052 int top_rank = 0; 1053 for ( ; stack_top >= 0; stack_top) { 1054 int u = stack[stack_top], v; 1055 int rank_u = _rank[u]; 1056 1057 LargeCost rc, pi_u = _pi[u]; 1058 int last_out = _first_out[u+1]; 1059 for (int a = _first_out[u]; a != last_out; ++a) { 1060 if (_res_cap[a] > 0) { 1061 v = _target[a]; 1062 rc = _cost[a] + pi_u  _pi[v]; 1063 if (rc < 0) { 1064 LargeCost nrc = static_cast<LargeCost>((rc  0.5) / _epsilon); 1065 if (nrc < LargeCost(_max_rank)) { 1066 int new_rank_v = rank_u + static_cast<int>(nrc); 1067 if (new_rank_v > _rank[v]) { 1068 _rank[v] = new_rank_v; 1069 } 1070 } 1071 } 1072 } 1073 } 1074 1075 if (rank_u > 0) { 1076 top_rank = std::max(top_rank, rank_u); 1077 int bfirst = _buckets[rank_u]; 1078 _bucket_next[u] = bfirst; 1079 _bucket_prev[bfirst] = u; 1080 _buckets[rank_u] = u; 1081 } 1082 } 1083 1084 // Check if the current flow is epsilonoptimal 1085 if (top_rank == 0) { 1086 return true; 1087 } 1088 1089 // Process buckets in topdown order 1090 for (int rank = top_rank; rank > 0; rank) { 1091 while (_buckets[rank] != bucket_end) { 1092 // Remove the first node from the current bucket 1093 int u = _buckets[rank]; 1094 _buckets[rank] = _bucket_next[u]; 1095 1096 // Search the outgoing arcs of u 1097 LargeCost rc, pi_u = _pi[u]; 1098 int last_out = _first_out[u+1]; 1099 int v, old_rank_v, new_rank_v; 1100 for (int a = _first_out[u]; a != last_out; ++a) { 1101 if (_res_cap[a] > 0) { 1102 v = _target[a]; 1103 old_rank_v = _rank[v]; 1104 1105 if (old_rank_v < rank) { 1106 1107 // Compute the new rank of node v 1108 rc = _cost[a] + pi_u  _pi[v]; 1109 if (rc < 0) { 1110 new_rank_v = rank; 1111 } else { 1112 LargeCost nrc = rc / _epsilon; 1113 new_rank_v = 0; 1114 if (nrc < LargeCost(_max_rank)) { 1115 new_rank_v = rank  1  static_cast<int>(nrc); 1116 } 1117 } 1118 1119 // Change the rank of node v 1120 if (new_rank_v > old_rank_v) { 1121 _rank[v] = new_rank_v; 1122 1123 // Remove v from its old bucket 1124 if (old_rank_v > 0) { 1125 if (_buckets[old_rank_v] == v) { 1126 _buckets[old_rank_v] = _bucket_next[v]; 1127 } else { 1128 int pv = _bucket_prev[v], nv = _bucket_next[v]; 1129 _bucket_next[pv] = nv; 1130 _bucket_prev[nv] = pv; 1131 } 1132 } 1133 1134 // Insert v into its new bucket 1135 int nv = _buckets[new_rank_v]; 1136 _bucket_next[v] = nv; 1137 _bucket_prev[nv] = v; 1138 _buckets[new_rank_v] = v; 1139 } 1140 } 1141 } 1142 } 1143 1144 // Refine potential of node u 1145 _pi[u] = rank * _epsilon; 1146 } 1147 } 1148 1149 } 1150 1151 return false; 1152 } 1153 1154 // Find and cancel cycles in the admissible network and 1155 // determine topological order using DFS 1156 bool topologicalSort(IntVector &stack, int &stack_top) { 1157 const int MAX_CYCLE_CANCEL = 1; 1158 1159 BoolVector reached(_res_node_num, false); 1160 BoolVector processed(_res_node_num, false); 1161 IntVector pred(_res_node_num); 1162 for (int i = 0; i != _res_node_num; ++i) { 1163 _next_out[i] = _first_out[i]; 1164 } 1165 stack_top = 1; 1166 1167 int cycle_cnt = 0; 1168 for (int start = 0; start != _res_node_num; ++start) { 1169 if (reached[start]) continue; 1170 1171 // Start DFS search from this start node 1172 pred[start] = 1; 1173 int tip = start, v; 1174 while (true) { 1175 // Check the outgoing arcs of the current tip node 1176 reached[tip] = true; 1177 LargeCost pi_tip = _pi[tip]; 1178 int a, last_out = _first_out[tip+1]; 1179 for (a = _next_out[tip]; a != last_out; ++a) { 1180 if (_res_cap[a] > 0) { 1181 v = _target[a]; 1182 if (_cost[a] + pi_tip  _pi[v] < 0) { 1183 if (!reached[v]) { 1184 // A new node is reached 1185 reached[v] = true; 1186 pred[v] = tip; 1187 _next_out[tip] = a; 1188 tip = v; 1189 a = _next_out[tip]; 1190 last_out = _first_out[tip+1]; 1191 break; 1192 } 1193 else if (!processed[v]) { 1194 // A cycle is found 1195 ++cycle_cnt; 1196 _next_out[tip] = a; 1197 1198 // Find the minimum residual capacity along the cycle 1199 Value d, delta = _res_cap[a]; 1200 int u, delta_node = tip; 1201 for (u = tip; u != v; ) { 1202 u = pred[u]; 1203 d = _res_cap[_next_out[u]]; 1204 if (d <= delta) { 1205 delta = d; 1206 delta_node = u; 1207 } 1208 } 1209 1210 // Augment along the cycle 1211 _res_cap[a] = delta; 1212 _res_cap[_reverse[a]] += delta; 1213 for (u = tip; u != v; ) { 1214 u = pred[u]; 1215 int ca = _next_out[u]; 1216 _res_cap[ca] = delta; 1217 _res_cap[_reverse[ca]] += delta; 1218 } 1219 1220 // Check the maximum number of cycle canceling 1221 if (cycle_cnt >= MAX_CYCLE_CANCEL) { 1222 return false; 1223 } 1224 1225 // Roll back search to delta_node 1226 if (delta_node != tip) { 1227 for (u = tip; u != delta_node; u = pred[u]) { 1228 reached[u] = false; 1229 } 1230 tip = delta_node; 1231 a = _next_out[tip] + 1; 1232 last_out = _first_out[tip+1]; 1233 break; 1234 } 1235 } 1236 } 1237 } 1238 } 1239 1240 // Step back to the previous node 1241 if (a == last_out) { 1242 processed[tip] = true; 1243 stack[++stack_top] = tip; 1244 tip = pred[tip]; 1245 if (tip < 0) { 1246 // Finish DFS from the current start node 1247 break; 1248 } 1249 ++_next_out[tip]; 1250 } 1251 } 1252 1253 } 1254 1255 return (cycle_cnt == 0); 996 1256 } 997 1257 998 1258 // Global potential update heuristic 999 1259 void globalUpdate() { 1000 int bucket_end = _root + 1;1260 const int bucket_end = _root + 1; 1001 1261 1002 1262 // Initialize buckets … … 1005 1265 } 1006 1266 Value total_excess = 0; 1267 int b0 = bucket_end; 1007 1268 for (int i = 0; i != _res_node_num; ++i) { 1008 1269 if (_excess[i] < 0) { 1009 1270 _rank[i] = 0; 1010 _bucket_next[i] = _buckets[0];1011 _bucket_prev[ _buckets[0]] = i;1012 _buckets[0]= i;1271 _bucket_next[i] = b0; 1272 _bucket_prev[b0] = i; 1273 b0 = i; 1013 1274 } else { 1014 1275 total_excess += _excess[i]; … … 1017 1278 } 1018 1279 if (total_excess == 0) return; 1280 _buckets[0] = b0; 1019 1281 1020 1282 // Search the buckets … … 1026 1288 _buckets[r] = _bucket_next[u]; 1027 1289 1028 // Search the incom ming arcs of u1290 // Search the incoming arcs of u 1029 1291 LargeCost pi_u = _pi[u]; 1030 1292 int last_out = _first_out[u+1]; … … 1038 1300 LargeCost nrc = (_cost[ra] + _pi[v]  pi_u) / _epsilon; 1039 1301 int new_rank_v = old_rank_v; 1040 if (nrc < LargeCost(_max_rank)) 1041 new_rank_v = r + 1 + int(nrc); 1302 if (nrc < LargeCost(_max_rank)) { 1303 new_rank_v = r + 1 + static_cast<int>(nrc); 1304 } 1042 1305 1043 1306 // Change the rank of v … … 1051 1314 _buckets[old_rank_v] = _bucket_next[v]; 1052 1315 } else { 1053 _bucket_next[_bucket_prev[v]] = _bucket_next[v]; 1054 _bucket_prev[_bucket_next[v]] = _bucket_prev[v]; 1316 int pv = _bucket_prev[v], nv = _bucket_next[v]; 1317 _bucket_next[pv] = nv; 1318 _bucket_prev[nv] = pv; 1055 1319 } 1056 1320 } 1057 1321 1058 // Insert v to its new bucket 1059 _bucket_next[v] = _buckets[new_rank_v]; 1060 _bucket_prev[_buckets[new_rank_v]] = v; 1322 // Insert v into its new bucket 1323 int nv = _buckets[new_rank_v]; 1324 _bucket_next[v] = nv; 1325 _bucket_prev[nv] = v; 1061 1326 _buckets[new_rank_v] = v; 1062 1327 } … … 1087 1352 void startAugment(int max_length) { 1088 1353 // 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 * 1354 const int PRICE_REFINEMENT_LIMIT = 2; 1355 const double GLOBAL_UPDATE_FACTOR = 1.0; 1356 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * 1093 1357 (_res_node_num + _sup_node_num * _sup_node_num)); 1094 int next_update_limit = global_update_freq; 1095 1358 int next_global_update_limit = global_update_skip; 1359 1360 // Perform cost scaling phases 1361 IntVector path; 1362 BoolVector path_arc(_res_arc_num, false); 1096 1363 int relabel_cnt = 0; 1097 1098 // Perform cost scaling phases 1099 std::vector<int> path; 1364 int eps_phase_cnt = 0; 1100 1365 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1101 1366 1 : _epsilon / _alpha ) 1102 1367 { 1103 // Early termination heuristic 1104 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1105 if (earlyTermination()) break; 1368 ++eps_phase_cnt; 1369 1370 // Price refinement heuristic 1371 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1372 if (priceRefinement()) continue; 1106 1373 } 1107 1374 … … 1120 1387 1121 1388 // Find an augmenting path from the start node 1122 path.clear();1123 1389 int tip = start; 1124 while ( _excess[tip] >= 0 && int(path.size()) < max_length) {1390 while (int(path.size()) < max_length && _excess[tip] >= 0) { 1125 1391 int u; 1126 LargeCost min_red_cost, rc, pi_tip = _pi[tip]; 1392 LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max(); 1393 LargeCost pi_tip = _pi[tip]; 1127 1394 int last_out = _first_out[tip+1]; 1128 1395 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; 1396 if (_res_cap[a] > 0) { 1397 u = _target[a]; 1398 rc = _cost[a] + pi_tip  _pi[u]; 1399 if (rc < 0) { 1400 path.push_back(a); 1401 _next_out[tip] = a; 1402 if (path_arc[a]) { 1403 goto augment; // a cycle is found, stop path search 1404 } 1405 tip = u; 1406 path_arc[a] = true; 1407 goto next_step; 1408 } 1409 else if (rc < min_red_cost) { 1410 min_red_cost = rc; 1411 } 1135 1412 } 1136 1413 } 1137 1414 1138 1415 // Relabel tip node 1139 min_red_cost = std::numeric_limits<LargeCost>::max();1140 1416 if (tip != start) { 1141 1417 int ra = _reverse[path.back()]; 1142 min_red_cost = _cost[ra] + pi_tip  _pi[_target[ra]]; 1418 min_red_cost = 1419 std::min(min_red_cost, _cost[ra] + pi_tip  _pi[_target[ra]]); 1143 1420 } 1421 last_out = _next_out[tip]; 1144 1422 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; 1423 if (_res_cap[a] > 0) { 1424 rc = _cost[a] + pi_tip  _pi[_target[a]]; 1425 if (rc < min_red_cost) { 1426 min_red_cost = rc; 1427 } 1148 1428 } 1149 1429 } … … 1154 1434 // Step back 1155 1435 if (tip != start) { 1156 tip = _source[path.back()]; 1436 int pa = path.back(); 1437 path_arc[pa] = false; 1438 tip = _source[pa]; 1157 1439 path.pop_back(); 1158 1440 } … … 1162 1444 1163 1445 // Augment along the found path (as much flow as possible) 1446 augment: 1164 1447 Value delta; 1165 1448 int pa, u, v = start; … … 1168 1451 u = v; 1169 1452 v = _target[pa]; 1453 path_arc[pa] = false; 1170 1454 delta = std::min(_res_cap[pa], _excess[u]); 1171 1455 _res_cap[pa] = delta; … … 1173 1457 _excess[u] = delta; 1174 1458 _excess[v] += delta; 1175 if (_excess[v] > 0 && _excess[v] <= delta) 1459 if (_excess[v] > 0 && _excess[v] <= delta) { 1176 1460 _active_nodes.push_back(v); 1177 } 1461 } 1462 } 1463 path.clear(); 1178 1464 1179 1465 // Global update heuristic 1180 if (relabel_cnt >= next_ update_limit) {1466 if (relabel_cnt >= next_global_update_limit) { 1181 1467 globalUpdate(); 1182 next_update_limit += global_update_freq; 1183 } 1184 } 1185 } 1468 next_global_update_limit += global_update_skip; 1469 } 1470 } 1471 1472 } 1473 1186 1474 } 1187 1475 … … 1189 1477 void startPush() { 1190 1478 // Paramters for heuristics 1191 const int EARLY_TERM_EPSILON_LIMIT = 1000;1479 const int PRICE_REFINEMENT_LIMIT = 2; 1192 1480 const double GLOBAL_UPDATE_FACTOR = 2.0; 1193 1481 1194 const int global_update_ freq = int(GLOBAL_UPDATE_FACTOR *1482 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * 1195 1483 (_res_node_num + _sup_node_num * _sup_node_num)); 1196 int next_update_limit = global_update_freq; 1197 1198 int relabel_cnt = 0; 1484 int next_global_update_limit = global_update_skip; 1199 1485 1200 1486 // Perform cost scaling phases 1201 1487 BoolVector hyper(_res_node_num, false); 1202 1488 LargeCostVector hyper_cost(_res_node_num); 1489 int relabel_cnt = 0; 1490 int eps_phase_cnt = 0; 1203 1491 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1204 1492 1 : _epsilon / _alpha ) 1205 1493 { 1206 // Early termination heuristic 1207 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1208 if (earlyTermination()) break; 1494 ++eps_phase_cnt; 1495 1496 // Price refinement heuristic 1497 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1498 if (priceRefinement()) continue; 1209 1499 } 1210 1500 … … 1278 1568 std::numeric_limits<LargeCost>::max(); 1279 1569 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; 1570 if (_res_cap[a] > 0) { 1571 rc = _cost[a] + pi_n  _pi[_target[a]]; 1572 if (rc < min_red_cost) { 1573 min_red_cost = rc; 1574 } 1283 1575 } 1284 1576 } … … 1298 1590 1299 1591 // Global update heuristic 1300 if (relabel_cnt >= next_ update_limit) {1592 if (relabel_cnt >= next_global_update_limit) { 1301 1593 globalUpdate(); 1302 1594 for (int u = 0; u != _res_node_num; ++u) 1303 1595 hyper[u] = false; 1304 next_ update_limit += global_update_freq;1596 next_global_update_limit += global_update_skip; 1305 1597 } 1306 1598 }
Note: See TracChangeset
for help on using the changeset viewer.