Changes in lemon/network_simplex.h [1026:9312d6c89d02:978:cbf32bf95954] in lemon
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r1026 r978 48 48 /// flow problem. 49 49 /// 50 /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest51 /// i mplementations available in LEMON for this problem.52 /// Furthermore, this class supports both directions of the supply/demand53 /// inequalityconstraints. For more information, see \ref SupplyType.50 /// In general, %NetworkSimplex is the fastest implementation available 51 /// in LEMON for this problem. 52 /// Moreover, it supports both directions of the supply/demand inequality 53 /// constraints. For more information, see \ref SupplyType. 54 54 /// 55 55 /// Most of the parameters of the problem (except for the digraph) … … 64 64 /// algorithm. By default, it is the same as \c V. 65 65 /// 66 /// \warning Both \c V and \c C must be signed number types. 67 /// \warning All input data (capacities, supply values, and costs) must 66 /// \warning Both number types must be signed and all input data must 68 67 /// be integer. 69 68 /// … … 127 126 /// of the algorithm. 128 127 /// By default, \ref BLOCK_SEARCH "Block Search" is used, which 129 /// turend outto be the most efficient and the most robust on various128 /// proved to be the most efficient and the most robust on various 130 129 /// test inputs. 131 130 /// However, another pivot rule can be selected using the \ref run() … … 168 167 typedef std::vector<Value> ValueVector; 169 168 typedef std::vector<Cost> CostVector; 170 typedef std::vector<signed char> CharVector; 171 // Note: vector<signed char> is used instead of vector<ArcState> and 172 // vector<ArcDirection> for efficiency reasons 169 typedef std::vector<char> BoolVector; 170 // Note: vector<char> is used instead of vector<bool> for efficiency reasons 173 171 174 172 // State constants for arcs … … 179 177 }; 180 178 181 // Direction constants for tree arcs 182 enum ArcDirection { 183 DIR_DOWN = 1, 184 DIR_UP = 1 185 }; 179 typedef std::vector<signed char> StateVector; 180 // Note: vector<signed char> is used instead of vector<ArcState> for 181 // efficiency reasons 186 182 187 183 private: … … 222 218 IntVector _succ_num; 223 219 IntVector _last_succ; 224 CharVector _pred_dir;225 CharVector _state;226 220 IntVector _dirty_revs; 221 BoolVector _forward; 222 StateVector _state; 227 223 int _root; 228 224 229 225 // Temporary data used in the current pivot iteration 230 226 int in_arc, join, u_in, v_in, u_out, v_out; 227 int first, second, right, last; 228 int stem, par_stem, new_stem; 231 229 Value delta; 232 230 … … 253 251 const IntVector &_target; 254 252 const CostVector &_cost; 255 const CharVector &_state;253 const StateVector &_state; 256 254 const CostVector &_pi; 257 255 int &_in_arc; … … 305 303 const IntVector &_target; 306 304 const CostVector &_cost; 307 const CharVector &_state;305 const StateVector &_state; 308 306 const CostVector &_pi; 309 307 int &_in_arc; … … 344 342 const IntVector &_target; 345 343 const CostVector &_cost; 346 const CharVector &_state;344 const StateVector &_state; 347 345 const CostVector &_pi; 348 346 int &_in_arc; … … 417 415 const IntVector &_target; 418 416 const CostVector &_cost; 419 const CharVector &_state;417 const StateVector &_state; 420 418 const CostVector &_pi; 421 419 int &_in_arc; … … 520 518 const IntVector &_target; 521 519 const CostVector &_cost; 522 const CharVector &_state;520 const StateVector &_state; 523 521 const CostVector &_pi; 524 522 int &_in_arc; … … 573 571 // Check the current candidate list 574 572 int e; 575 Cost c;576 573 for (int i = 0; i != _curr_length; ++i) { 577 574 e = _candidates[i]; 578 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 579 if (c < 0) { 580 _cand_cost[e] = c; 581 } else { 575 _cand_cost[e] = _state[e] * 576 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 577 if (_cand_cost[e] >= 0) { 582 578 _candidates[i] = _candidates[_curr_length]; 583 579 } … … 589 585 590 586 for (e = _next_arc; e != _search_arc_num; ++e) { 591 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]);592 if (c < 0) {593 _cand_cost[e] = c;587 _cand_cost[e] = _state[e] * 588 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 589 if (_cand_cost[e] < 0) { 594 590 _candidates[_curr_length++] = e; 595 591 } … … 638 634 /// 639 635 /// \param graph The digraph the algorithm runs on. 640 /// \param arc_mixing Indicate if the arcs willbe stored in a636 /// \param arc_mixing Indicate if the arcs have to be stored in a 641 637 /// mixed order in the internal data structure. 642 /// In general, it leads to similar performance as using the original 643 /// arc order, but it makes the algorithm more robust and in special 644 /// cases, even significantly faster. Therefore, it is enabled by default. 645 NetworkSimplex(const GR& graph, bool arc_mixing = true) : 638 /// In special cases, it could lead to better overall performance, 639 /// but it is usually slower. Therefore it is disabled by default. 640 NetworkSimplex(const GR& graph, bool arc_mixing = false) : 646 641 _graph(graph), _node_id(graph), _arc_id(graph), 647 642 _arc_mixing(arc_mixing), … … 736 731 /// 737 732 /// \return <tt>(*this)</tt> 738 ///739 /// \sa supplyType()740 733 template<typename SupplyMap> 741 734 NetworkSimplex& supplyMap(const SupplyMap& map) { … … 754 747 /// 755 748 /// Using this function has the same effect as using \ref supplyMap() 756 /// with a map in which \c k is assigned to \c s, \c k is749 /// with such a map in which \c k is assigned to \c s, \c k is 757 750 /// assigned to \c t and all other nodes have zero supply value. 758 751 /// … … 921 914 _parent.resize(all_node_num); 922 915 _pred.resize(all_node_num); 923 _ pred_dir.resize(all_node_num);916 _forward.resize(all_node_num); 924 917 _thread.resize(all_node_num); 925 918 _rev_thread.resize(all_node_num); … … 935 928 if (_arc_mixing) { 936 929 // Store the arcs in a mixed order 937 const int skip = std::max(_arc_num / _node_num, 3);930 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 938 931 int i = 0, j = 0; 939 932 for (ArcIt a(_graph); a != INVALID; ++a) { … … 941 934 _source[i] = _node_id[_graph.source(a)]; 942 935 _target[i] = _node_id[_graph.target(a)]; 943 if ((i += skip) >= _arc_num) i = ++j;936 if ((i += k) >= _arc_num) i = ++j; 944 937 } 945 938 } else { … … 1124 1117 _state[e] = STATE_TREE; 1125 1118 if (_supply[u] >= 0) { 1126 _ pred_dir[u] = DIR_UP;1119 _forward[u] = true; 1127 1120 _pi[u] = 0; 1128 1121 _source[e] = u; … … 1131 1124 _cost[e] = 0; 1132 1125 } else { 1133 _ pred_dir[u] = DIR_DOWN;1126 _forward[u] = false; 1134 1127 _pi[u] = ART_COST; 1135 1128 _source[e] = _root; … … 1151 1144 _last_succ[u] = u; 1152 1145 if (_supply[u] >= 0) { 1153 _ pred_dir[u] = DIR_UP;1146 _forward[u] = true; 1154 1147 _pi[u] = 0; 1155 1148 _pred[u] = e; … … 1161 1154 _state[e] = STATE_TREE; 1162 1155 } else { 1163 _ pred_dir[u] = DIR_DOWN;1156 _forward[u] = false; 1164 1157 _pi[u] = ART_COST; 1165 1158 _pred[u] = f; … … 1192 1185 _last_succ[u] = u; 1193 1186 if (_supply[u] <= 0) { 1194 _ pred_dir[u] = DIR_DOWN;1187 _forward[u] = false; 1195 1188 _pi[u] = 0; 1196 1189 _pred[u] = e; … … 1202 1195 _state[e] = STATE_TREE; 1203 1196 } else { 1204 _ pred_dir[u] = DIR_UP;1197 _forward[u] = true; 1205 1198 _pi[u] = ART_COST; 1206 1199 _pred[u] = f; … … 1245 1238 // Initialize first and second nodes according to the direction 1246 1239 // of the cycle 1247 int first, second;1248 1240 if (_state[in_arc] == STATE_LOWER) { 1249 1241 first = _source[in_arc]; … … 1255 1247 delta = _cap[in_arc]; 1256 1248 int result = 0; 1257 Value c,d;1249 Value d; 1258 1250 int e; 1259 1251 1260 // Search the cycle form the first node to the join node1252 // Search the cycle along the path form the first node to the root 1261 1253 for (int u = first; u != join; u = _parent[u]) { 1262 1254 e = _pred[u]; 1263 d = _flow[e]; 1264 if (_pred_dir[u] == DIR_DOWN) { 1265 c = _cap[e]; 1266 d = c >= MAX ? INF : c  d; 1267 } 1255 d = _forward[u] ? 1256 _flow[e] : (_cap[e] >= MAX ? INF : _cap[e]  _flow[e]); 1268 1257 if (d < delta) { 1269 1258 delta = d; … … 1272 1261 } 1273 1262 } 1274 1275 // Search the cycle form the second node to the join node 1263 // Search the cycle along the path form the second node to the root 1276 1264 for (int u = second; u != join; u = _parent[u]) { 1277 1265 e = _pred[u]; 1278 d = _flow[e]; 1279 if (_pred_dir[u] == DIR_UP) { 1280 c = _cap[e]; 1281 d = c >= MAX ? INF : c  d; 1282 } 1266 d = _forward[u] ? 1267 (_cap[e] >= MAX ? INF : _cap[e]  _flow[e]) : _flow[e]; 1283 1268 if (d <= delta) { 1284 1269 delta = d; … … 1305 1290 _flow[in_arc] += val; 1306 1291 for (int u = _source[in_arc]; u != join; u = _parent[u]) { 1307 _flow[_pred[u]] = _pred_dir[u] *val;1292 _flow[_pred[u]] += _forward[u] ? val : val; 1308 1293 } 1309 1294 for (int u = _target[in_arc]; u != join; u = _parent[u]) { 1310 _flow[_pred[u]] += _ pred_dir[u] *val;1295 _flow[_pred[u]] += _forward[u] ? val : val; 1311 1296 } 1312 1297 } … … 1323 1308 // Update the tree structure 1324 1309 void updateTreeStructure() { 1310 int u, w; 1325 1311 int old_rev_thread = _rev_thread[u_out]; 1326 1312 int old_succ_num = _succ_num[u_out]; … … 1328 1314 v_out = _parent[u_out]; 1329 1315 1330 // Check if u_in and u_out coincide 1331 if (u_in == u_out) { 1332 // Update _parent, _pred, _pred_dir 1333 _parent[u_in] = v_in; 1334 _pred[u_in] = in_arc; 1335 _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN; 1336 1337 // Update _thread and _rev_thread 1338 if (_thread[v_in] != u_out) { 1339 int after = _thread[old_last_succ]; 1340 _thread[old_rev_thread] = after; 1341 _rev_thread[after] = old_rev_thread; 1342 after = _thread[v_in]; 1343 _thread[v_in] = u_out; 1344 _rev_thread[u_out] = v_in; 1345 _thread[old_last_succ] = after; 1346 _rev_thread[after] = old_last_succ; 1347 } 1316 u = _last_succ[u_in]; // the last successor of u_in 1317 right = _thread[u]; // the node after it 1318 1319 // Handle the case when old_rev_thread equals to v_in 1320 // (it also means that join and v_out coincide) 1321 if (old_rev_thread == v_in) { 1322 last = _thread[_last_succ[u_out]]; 1348 1323 } else { 1349 // Handle the case when old_rev_thread equals to v_in 1350 // (it also means that join and v_out coincide) 1351 int thread_continue = old_rev_thread == v_in ? 1352 _thread[old_last_succ] : _thread[v_in]; 1353 1354 // Update _thread and _parent along the stem nodes (i.e. the nodes 1355 // between u_in and u_out, whose parent have to be changed) 1356 int stem = u_in; // the current stem node 1357 int par_stem = v_in; // the new parent of stem 1358 int next_stem; // the next stem node 1359 int last = _last_succ[u_in]; // the last successor of stem 1360 int before, after = _thread[last]; 1361 _thread[v_in] = u_in; 1362 _dirty_revs.clear(); 1363 _dirty_revs.push_back(v_in); 1364 while (stem != u_out) { 1365 // Insert the next stem node into the thread list 1366 next_stem = _parent[stem]; 1367 _thread[last] = next_stem; 1368 _dirty_revs.push_back(last); 1369 1370 // Remove the subtree of stem from the thread list 1371 before = _rev_thread[stem]; 1372 _thread[before] = after; 1373 _rev_thread[after] = before; 1374 1375 // Change the parent node and shift stem nodes 1376 _parent[stem] = par_stem; 1377 par_stem = stem; 1378 stem = next_stem; 1379 1380 // Update last and after 1381 last = _last_succ[stem] == _last_succ[par_stem] ? 1382 _rev_thread[par_stem] : _last_succ[stem]; 1383 after = _thread[last]; 1384 } 1385 _parent[u_out] = par_stem; 1386 _thread[last] = thread_continue; 1387 _rev_thread[thread_continue] = last; 1388 _last_succ[u_out] = last; 1389 1390 // Remove the subtree of u_out from the thread list except for 1391 // the case when old_rev_thread equals to v_in 1392 if (old_rev_thread != v_in) { 1393 _thread[old_rev_thread] = after; 1394 _rev_thread[after] = old_rev_thread; 1395 } 1396 1397 // Update _rev_thread using the new _thread values 1398 for (int i = 0; i != int(_dirty_revs.size()); ++i) { 1399 int u = _dirty_revs[i]; 1400 _rev_thread[_thread[u]] = u; 1401 } 1402 1403 // Update _pred, _pred_dir, _last_succ and _succ_num for the 1404 // stem nodes from u_out to u_in 1405 int tmp_sc = 0, tmp_ls = _last_succ[u_out]; 1406 for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) { 1407 _pred[u] = _pred[p]; 1408 _pred_dir[u] = _pred_dir[p]; 1409 tmp_sc += _succ_num[u]  _succ_num[p]; 1410 _succ_num[u] = tmp_sc; 1411 _last_succ[p] = tmp_ls; 1412 } 1413 _pred[u_in] = in_arc; 1414 _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN; 1415 _succ_num[u_in] = old_succ_num; 1324 last = _thread[v_in]; 1325 } 1326 1327 // Update _thread and _parent along the stem nodes (i.e. the nodes 1328 // between u_in and u_out, whose parent have to be changed) 1329 _thread[v_in] = stem = u_in; 1330 _dirty_revs.clear(); 1331 _dirty_revs.push_back(v_in); 1332 par_stem = v_in; 1333 while (stem != u_out) { 1334 // Insert the next stem node into the thread list 1335 new_stem = _parent[stem]; 1336 _thread[u] = new_stem; 1337 _dirty_revs.push_back(u); 1338 1339 // Remove the subtree of stem from the thread list 1340 w = _rev_thread[stem]; 1341 _thread[w] = right; 1342 _rev_thread[right] = w; 1343 1344 // Change the parent node and shift stem nodes 1345 _parent[stem] = par_stem; 1346 par_stem = stem; 1347 stem = new_stem; 1348 1349 // Update u and right 1350 u = _last_succ[stem] == _last_succ[par_stem] ? 1351 _rev_thread[par_stem] : _last_succ[stem]; 1352 right = _thread[u]; 1353 } 1354 _parent[u_out] = par_stem; 1355 _thread[u] = last; 1356 _rev_thread[last] = u; 1357 _last_succ[u_out] = u; 1358 1359 // Remove the subtree of u_out from the thread list except for 1360 // the case when old_rev_thread equals to v_in 1361 // (it also means that join and v_out coincide) 1362 if (old_rev_thread != v_in) { 1363 _thread[old_rev_thread] = right; 1364 _rev_thread[right] = old_rev_thread; 1365 } 1366 1367 // Update _rev_thread using the new _thread values 1368 for (int i = 0; i != int(_dirty_revs.size()); ++i) { 1369 u = _dirty_revs[i]; 1370 _rev_thread[_thread[u]] = u; 1371 } 1372 1373 // Update _pred, _forward, _last_succ and _succ_num for the 1374 // stem nodes from u_out to u_in 1375 int tmp_sc = 0, tmp_ls = _last_succ[u_out]; 1376 u = u_out; 1377 while (u != u_in) { 1378 w = _parent[u]; 1379 _pred[u] = _pred[w]; 1380 _forward[u] = !_forward[w]; 1381 tmp_sc += _succ_num[u]  _succ_num[w]; 1382 _succ_num[u] = tmp_sc; 1383 _last_succ[w] = tmp_ls; 1384 u = w; 1385 } 1386 _pred[u_in] = in_arc; 1387 _forward[u_in] = (u_in == _source[in_arc]); 1388 _succ_num[u_in] = old_succ_num; 1389 1390 // Set limits for updating _last_succ form v_in and v_out 1391 // towards the root 1392 int up_limit_in = 1; 1393 int up_limit_out = 1; 1394 if (_last_succ[join] == v_in) { 1395 up_limit_out = join; 1396 } else { 1397 up_limit_in = join; 1416 1398 } 1417 1399 1418 1400 // Update _last_succ from v_in towards the root 1419 int up_limit_out = _last_succ[join] == v_in ? join : 1; 1420 int last_succ_out = _last_succ[u_out]; 1421 for (int u = v_in; u != 1 && _last_succ[u] == v_in; u = _parent[u]) { 1422 _last_succ[u] = last_succ_out; 1423 } 1424 1401 for (u = v_in; u != up_limit_in && _last_succ[u] == v_in; 1402 u = _parent[u]) { 1403 _last_succ[u] = _last_succ[u_out]; 1404 } 1425 1405 // Update _last_succ from v_out towards the root 1426 1406 if (join != old_rev_thread && v_in != old_rev_thread) { 1427 for ( intu = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;1407 for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; 1428 1408 u = _parent[u]) { 1429 1409 _last_succ[u] = old_rev_thread; 1430 1410 } 1431 } 1432 else if (last_succ_out != old_last_succ) { 1433 for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; 1411 } else { 1412 for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; 1434 1413 u = _parent[u]) { 1435 _last_succ[u] = last_succ_out;1414 _last_succ[u] = _last_succ[u_out]; 1436 1415 } 1437 1416 } 1438 1417 1439 1418 // Update _succ_num from v_in to join 1440 for ( intu = v_in; u != join; u = _parent[u]) {1419 for (u = v_in; u != join; u = _parent[u]) { 1441 1420 _succ_num[u] += old_succ_num; 1442 1421 } 1443 1422 // Update _succ_num from v_out to join 1444 for ( intu = v_out; u != join; u = _parent[u]) {1423 for (u = v_out; u != join; u = _parent[u]) { 1445 1424 _succ_num[u] = old_succ_num; 1446 1425 } 1447 1426 } 1448 1427 1449 // Update potentials in the subtree that has been moved1428 // Update potentials 1450 1429 void updatePotential() { 1451 Cost sigma = _pi[v_in]  _pi[u_in]  1452 _pred_dir[u_in] * _cost[in_arc]; 1430 Cost sigma = _forward[u_in] ? 1431 _pi[v_in]  _pi[u_in]  _cost[_pred[u_in]] : 1432 _pi[v_in]  _pi[u_in] + _cost[_pred[u_in]]; 1433 // Update potentials in the subtree, which has been moved 1453 1434 int end = _thread[_last_succ[u_in]]; 1454 1435 for (int u = u_in; u != end; u = _thread[u]) {
Note: See TracChangeset
for help on using the changeset viewer.