Changeset 1047:ddd3c0d3d9bf in lemon for lemon/cost_scaling.h
 Timestamp:
 03/15/11 19:32:21 (8 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/cost_scaling.h
r1046 r1047 981 981 } 982 982 983 // Early termination heuristic 984 bool earlyTermination() { 985 const double EARLY_TERM_FACTOR = 3.0; 986 987 // Build a static residual graph 988 _arc_vec.clear(); 989 _cost_vec.clear(); 990 for (int j = 0; j != _res_arc_num; ++j) { 991 if (_res_cap[j] > 0) { 992 _arc_vec.push_back(IntPair(_source[j], _target[j])); 993 _cost_vec.push_back(_cost[j] + 1); 994 } 995 } 996 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 997 998 // Run BellmanFord algorithm to check if the current flow is optimal 999 BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); 1000 bf.init(0); 1001 bool done = false; 1002 int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num))); 1003 for (int i = 0; i < K && !done; ++i) { 1004 done = bf.processNextWeakRound(); 1005 } 1006 return done; 983 // Price (potential) refinement heuristic 984 bool priceRefinement() { 985 986 // Stack for stroing the topological order 987 IntVector stack(_res_node_num); 988 int stack_top; 989 990 // Perform phases 991 while (topologicalSort(stack, stack_top)) { 992 993 // Compute node ranks in the acyclic admissible network and 994 // store the nodes in buckets 995 for (int i = 0; i != _res_node_num; ++i) { 996 _rank[i] = 0; 997 } 998 const int bucket_end = _root + 1; 999 for (int r = 0; r != _max_rank; ++r) { 1000 _buckets[r] = bucket_end; 1001 } 1002 int top_rank = 0; 1003 for ( ; stack_top >= 0; stack_top) { 1004 int u = stack[stack_top], v; 1005 int rank_u = _rank[u]; 1006 1007 LargeCost rc, pi_u = _pi[u]; 1008 int last_out = _first_out[u+1]; 1009 for (int a = _first_out[u]; a != last_out; ++a) { 1010 if (_res_cap[a] > 0) { 1011 v = _target[a]; 1012 rc = _cost[a] + pi_u  _pi[v]; 1013 if (rc < 0) { 1014 LargeCost nrc = static_cast<LargeCost>((rc  0.5) / _epsilon); 1015 if (nrc < LargeCost(_max_rank)) { 1016 int new_rank_v = rank_u + static_cast<int>(nrc); 1017 if (new_rank_v > _rank[v]) { 1018 _rank[v] = new_rank_v; 1019 } 1020 } 1021 } 1022 } 1023 } 1024 1025 if (rank_u > 0) { 1026 top_rank = std::max(top_rank, rank_u); 1027 int bfirst = _buckets[rank_u]; 1028 _bucket_next[u] = bfirst; 1029 _bucket_prev[bfirst] = u; 1030 _buckets[rank_u] = u; 1031 } 1032 } 1033 1034 // Check if the current flow is epsilonoptimal 1035 if (top_rank == 0) { 1036 return true; 1037 } 1038 1039 // Process buckets in topdown order 1040 for (int rank = top_rank; rank > 0; rank) { 1041 while (_buckets[rank] != bucket_end) { 1042 // Remove the first node from the current bucket 1043 int u = _buckets[rank]; 1044 _buckets[rank] = _bucket_next[u]; 1045 1046 // Search the outgoing arcs of u 1047 LargeCost rc, pi_u = _pi[u]; 1048 int last_out = _first_out[u+1]; 1049 int v, old_rank_v, new_rank_v; 1050 for (int a = _first_out[u]; a != last_out; ++a) { 1051 if (_res_cap[a] > 0) { 1052 v = _target[a]; 1053 old_rank_v = _rank[v]; 1054 1055 if (old_rank_v < rank) { 1056 1057 // Compute the new rank of node v 1058 rc = _cost[a] + pi_u  _pi[v]; 1059 if (rc < 0) { 1060 new_rank_v = rank; 1061 } else { 1062 LargeCost nrc = rc / _epsilon; 1063 new_rank_v = 0; 1064 if (nrc < LargeCost(_max_rank)) { 1065 new_rank_v = rank  1  static_cast<int>(nrc); 1066 } 1067 } 1068 1069 // Change the rank of node v 1070 if (new_rank_v > old_rank_v) { 1071 _rank[v] = new_rank_v; 1072 1073 // Remove v from its old bucket 1074 if (old_rank_v > 0) { 1075 if (_buckets[old_rank_v] == v) { 1076 _buckets[old_rank_v] = _bucket_next[v]; 1077 } else { 1078 int pv = _bucket_prev[v], nv = _bucket_next[v]; 1079 _bucket_next[pv] = nv; 1080 _bucket_prev[nv] = pv; 1081 } 1082 } 1083 1084 // Insert v into its new bucket 1085 int nv = _buckets[new_rank_v]; 1086 _bucket_next[v] = nv; 1087 _bucket_prev[nv] = v; 1088 _buckets[new_rank_v] = v; 1089 } 1090 } 1091 } 1092 } 1093 1094 // Refine potential of node u 1095 _pi[u] = rank * _epsilon; 1096 } 1097 } 1098 1099 } 1100 1101 return false; 1102 } 1103 1104 // Find and cancel cycles in the admissible network and 1105 // determine topological order using DFS 1106 bool topologicalSort(IntVector &stack, int &stack_top) { 1107 const int MAX_CYCLE_CANCEL = 1; 1108 1109 BoolVector reached(_res_node_num, false); 1110 BoolVector processed(_res_node_num, false); 1111 IntVector pred(_res_node_num); 1112 for (int i = 0; i != _res_node_num; ++i) { 1113 _next_out[i] = _first_out[i]; 1114 } 1115 stack_top = 1; 1116 1117 int cycle_cnt = 0; 1118 for (int start = 0; start != _res_node_num; ++start) { 1119 if (reached[start]) continue; 1120 1121 // Start DFS search from this start node 1122 pred[start] = 1; 1123 int tip = start, v; 1124 while (true) { 1125 // Check the outgoing arcs of the current tip node 1126 reached[tip] = true; 1127 LargeCost pi_tip = _pi[tip]; 1128 int a, last_out = _first_out[tip+1]; 1129 for (a = _next_out[tip]; a != last_out; ++a) { 1130 if (_res_cap[a] > 0) { 1131 v = _target[a]; 1132 if (_cost[a] + pi_tip  _pi[v] < 0) { 1133 if (!reached[v]) { 1134 // A new node is reached 1135 reached[v] = true; 1136 pred[v] = tip; 1137 _next_out[tip] = a; 1138 tip = v; 1139 a = _next_out[tip]; 1140 last_out = _first_out[tip+1]; 1141 break; 1142 } 1143 else if (!processed[v]) { 1144 // A cycle is found 1145 ++cycle_cnt; 1146 _next_out[tip] = a; 1147 1148 // Find the minimum residual capacity along the cycle 1149 Value d, delta = _res_cap[a]; 1150 int u, delta_node = tip; 1151 for (u = tip; u != v; ) { 1152 u = pred[u]; 1153 d = _res_cap[_next_out[u]]; 1154 if (d <= delta) { 1155 delta = d; 1156 delta_node = u; 1157 } 1158 } 1159 1160 // Augment along the cycle 1161 _res_cap[a] = delta; 1162 _res_cap[_reverse[a]] += delta; 1163 for (u = tip; u != v; ) { 1164 u = pred[u]; 1165 int ca = _next_out[u]; 1166 _res_cap[ca] = delta; 1167 _res_cap[_reverse[ca]] += delta; 1168 } 1169 1170 // Check the maximum number of cycle canceling 1171 if (cycle_cnt >= MAX_CYCLE_CANCEL) { 1172 return false; 1173 } 1174 1175 // Roll back search to delta_node 1176 if (delta_node != tip) { 1177 for (u = tip; u != delta_node; u = pred[u]) { 1178 reached[u] = false; 1179 } 1180 tip = delta_node; 1181 a = _next_out[tip] + 1; 1182 last_out = _first_out[tip+1]; 1183 break; 1184 } 1185 } 1186 } 1187 } 1188 } 1189 1190 // Step back to the previous node 1191 if (a == last_out) { 1192 processed[tip] = true; 1193 stack[++stack_top] = tip; 1194 tip = pred[tip]; 1195 if (tip < 0) { 1196 // Finish DFS from the current start node 1197 break; 1198 } 1199 ++_next_out[tip]; 1200 } 1201 } 1202 1203 } 1204 1205 return (cycle_cnt == 0); 1007 1206 } 1008 1207 … … 1103 1302 void startAugment(int max_length) { 1104 1303 // Paramters for heuristics 1105 const int EARLY_TERM_EPSILON_LIMIT = 1000;1304 const int PRICE_REFINEMENT_LIMIT = 2; 1106 1305 const double GLOBAL_UPDATE_FACTOR = 1.0; 1107 1306 const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR * … … 1113 1312 BoolVector path_arc(_res_arc_num, false); 1114 1313 int relabel_cnt = 0; 1314 int eps_phase_cnt = 0; 1115 1315 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1116 1316 1 : _epsilon / _alpha ) 1117 1317 { 1118 // Early termination heuristic 1119 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1120 if (earlyTermination()) break; 1318 ++eps_phase_cnt; 1319 1320 // Price refinement heuristic 1321 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1322 if (priceRefinement()) continue; 1121 1323 } 1122 1324 … … 1225 1427 void startPush() { 1226 1428 // Paramters for heuristics 1227 const int EARLY_TERM_EPSILON_LIMIT = 1000;1429 const int PRICE_REFINEMENT_LIMIT = 2; 1228 1430 const double GLOBAL_UPDATE_FACTOR = 2.0; 1229 1431 … … 1236 1438 LargeCostVector hyper_cost(_res_node_num); 1237 1439 int relabel_cnt = 0; 1440 int eps_phase_cnt = 0; 1238 1441 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1239 1442 1 : _epsilon / _alpha ) 1240 1443 { 1241 // Early termination heuristic 1242 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 1243 if (earlyTermination()) break; 1444 ++eps_phase_cnt; 1445 1446 // Price refinement heuristic 1447 if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) { 1448 if (priceRefinement()) continue; 1244 1449 } 1245 1450
Note: See TracChangeset
for help on using the changeset viewer.