r655:6ac5d9ae1d3d
Fri, 03 Apr 2009 18:59:15
[Not Reviewed]
0 comments (0 inline)
0
1
0
| ... | ... |
@@ -45,25 +45,26 @@ |
| 45 | 45 |
/// simplex method directly for the minimum cost flow problem. |
| 46 | 46 |
/// It is one of the most efficient solution methods. |
| 47 | 47 |
/// |
| 48 | 48 |
/// In general this class is the fastest implementation available |
| 49 | 49 |
/// in LEMON for the minimum cost flow problem. |
| 50 | 50 |
/// |
| 51 | 51 |
/// \tparam GR The digraph type the algorithm runs on. |
| 52 | 52 |
/// \tparam F The value type used for flow amounts, capacity bounds |
| 53 | 53 |
/// and supply values in the algorithm. By default it is \c int. |
| 54 | 54 |
/// \tparam C The value type used for costs and potentials in the |
| 55 | 55 |
/// algorithm. By default it is the same as \c F. |
| 56 | 56 |
/// |
| 57 |
/// \warning Both value types must be signed |
|
| 57 |
/// \warning Both value types must be signed and all input data must |
|
| 58 |
/// be integer. |
|
| 58 | 59 |
/// |
| 59 | 60 |
/// \note %NetworkSimplex provides five different pivot rule |
| 60 | 61 |
/// implementations. For more information see \ref PivotRule. |
| 61 | 62 |
template <typename GR, typename F = int, typename C = F> |
| 62 | 63 |
class NetworkSimplex |
| 63 | 64 |
{
|
| 64 | 65 |
public: |
| 65 | 66 |
|
| 66 | 67 |
/// The flow type of the algorithm |
| 67 | 68 |
typedef F Flow; |
| 68 | 69 |
/// The cost type of the algorithm |
| 69 | 70 |
typedef C Cost; |
| ... | ... |
@@ -1035,91 +1036,105 @@ |
| 1035 | 1036 |
_supply[_root] = 0; |
| 1036 | 1037 |
_pi[_root] = 0; |
| 1037 | 1038 |
|
| 1038 | 1039 |
// Store the arcs in a mixed order |
| 1039 | 1040 |
int k = std::max(int(sqrt(_arc_num)), 10); |
| 1040 | 1041 |
int i = 0; |
| 1041 | 1042 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
| 1042 | 1043 |
_arc_ref[i] = e; |
| 1043 | 1044 |
if ((i += k) >= _arc_num) i = (i % k) + 1; |
| 1044 | 1045 |
} |
| 1045 | 1046 |
|
| 1046 | 1047 |
// Initialize arc maps |
| 1047 |
Flow max_cap = std::numeric_limits<Flow>::max(); |
|
| 1048 |
Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
|
| 1048 |
Flow inf_cap = |
|
| 1049 |
std::numeric_limits<Flow>::has_infinity ? |
|
| 1050 |
std::numeric_limits<Flow>::infinity() : |
|
| 1051 |
std::numeric_limits<Flow>::max(); |
|
| 1049 | 1052 |
if (_pupper && _pcost) {
|
| 1050 | 1053 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1051 | 1054 |
Arc e = _arc_ref[i]; |
| 1052 | 1055 |
_source[i] = _node_id[_graph.source(e)]; |
| 1053 | 1056 |
_target[i] = _node_id[_graph.target(e)]; |
| 1054 | 1057 |
_cap[i] = (*_pupper)[e]; |
| 1055 | 1058 |
_cost[i] = (*_pcost)[e]; |
| 1056 | 1059 |
_flow[i] = 0; |
| 1057 | 1060 |
_state[i] = STATE_LOWER; |
| 1058 | 1061 |
} |
| 1059 | 1062 |
} else {
|
| 1060 | 1063 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1061 | 1064 |
Arc e = _arc_ref[i]; |
| 1062 | 1065 |
_source[i] = _node_id[_graph.source(e)]; |
| 1063 | 1066 |
_target[i] = _node_id[_graph.target(e)]; |
| 1064 | 1067 |
_flow[i] = 0; |
| 1065 | 1068 |
_state[i] = STATE_LOWER; |
| 1066 | 1069 |
} |
| 1067 | 1070 |
if (_pupper) {
|
| 1068 | 1071 |
for (int i = 0; i != _arc_num; ++i) |
| 1069 | 1072 |
_cap[i] = (*_pupper)[_arc_ref[i]]; |
| 1070 | 1073 |
} else {
|
| 1071 | 1074 |
for (int i = 0; i != _arc_num; ++i) |
| 1072 |
_cap[i] = |
|
| 1075 |
_cap[i] = inf_cap; |
|
| 1073 | 1076 |
} |
| 1074 | 1077 |
if (_pcost) {
|
| 1075 | 1078 |
for (int i = 0; i != _arc_num; ++i) |
| 1076 | 1079 |
_cost[i] = (*_pcost)[_arc_ref[i]]; |
| 1077 | 1080 |
} else {
|
| 1078 | 1081 |
for (int i = 0; i != _arc_num; ++i) |
| 1079 | 1082 |
_cost[i] = 1; |
| 1080 | 1083 |
} |
| 1081 | 1084 |
} |
| 1082 | 1085 |
|
| 1086 |
// Initialize artifical cost |
|
| 1087 |
Cost art_cost; |
|
| 1088 |
if (std::numeric_limits<Cost>::is_exact) {
|
|
| 1089 |
art_cost = std::numeric_limits<Cost>::max() / 4 + 1; |
|
| 1090 |
} else {
|
|
| 1091 |
art_cost = std::numeric_limits<Cost>::min(); |
|
| 1092 |
for (int i = 0; i != _arc_num; ++i) {
|
|
| 1093 |
if (_cost[i] > art_cost) art_cost = _cost[i]; |
|
| 1094 |
} |
|
| 1095 |
art_cost = (art_cost + 1) * _node_num; |
|
| 1096 |
} |
|
| 1097 |
|
|
| 1083 | 1098 |
// Remove non-zero lower bounds |
| 1084 | 1099 |
if (_plower) {
|
| 1085 | 1100 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1086 | 1101 |
Flow c = (*_plower)[_arc_ref[i]]; |
| 1087 | 1102 |
if (c != 0) {
|
| 1088 | 1103 |
_cap[i] -= c; |
| 1089 | 1104 |
_supply[_source[i]] -= c; |
| 1090 | 1105 |
_supply[_target[i]] += c; |
| 1091 | 1106 |
} |
| 1092 | 1107 |
} |
| 1093 | 1108 |
} |
| 1094 | 1109 |
|
| 1095 | 1110 |
// Add artificial arcs and initialize the spanning tree data structure |
| 1096 | 1111 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
| 1097 | 1112 |
_thread[u] = u + 1; |
| 1098 | 1113 |
_rev_thread[u + 1] = u; |
| 1099 | 1114 |
_succ_num[u] = 1; |
| 1100 | 1115 |
_last_succ[u] = u; |
| 1101 | 1116 |
_parent[u] = _root; |
| 1102 | 1117 |
_pred[u] = e; |
| 1103 |
_cost[e] = max_cost; |
|
| 1104 |
_cap[e] = max_cap; |
|
| 1118 |
_cost[e] = art_cost; |
|
| 1119 |
_cap[e] = inf_cap; |
|
| 1105 | 1120 |
_state[e] = STATE_TREE; |
| 1106 | 1121 |
if (_supply[u] >= 0) {
|
| 1107 | 1122 |
_flow[e] = _supply[u]; |
| 1108 | 1123 |
_forward[u] = true; |
| 1109 |
_pi[u] = - |
|
| 1124 |
_pi[u] = -art_cost; |
|
| 1110 | 1125 |
} else {
|
| 1111 | 1126 |
_flow[e] = -_supply[u]; |
| 1112 | 1127 |
_forward[u] = false; |
| 1113 |
_pi[u] = |
|
| 1128 |
_pi[u] = art_cost; |
|
| 1114 | 1129 |
} |
| 1115 | 1130 |
} |
| 1116 | 1131 |
|
| 1117 | 1132 |
return true; |
| 1118 | 1133 |
} |
| 1119 | 1134 |
|
| 1120 | 1135 |
// Find the join node |
| 1121 | 1136 |
void findJoinNode() {
|
| 1122 | 1137 |
int u = _source[in_arc]; |
| 1123 | 1138 |
int v = _target[in_arc]; |
| 1124 | 1139 |
while (u != v) {
|
| 1125 | 1140 |
if (_succ_num[u] < _succ_num[v]) {
|
| ... | ... |
@@ -1318,42 +1333,30 @@ |
| 1318 | 1333 |
} |
| 1319 | 1334 |
// Update _succ_num from v_out to join |
| 1320 | 1335 |
for (u = v_out; u != join; u = _parent[u]) {
|
| 1321 | 1336 |
_succ_num[u] -= old_succ_num; |
| 1322 | 1337 |
} |
| 1323 | 1338 |
} |
| 1324 | 1339 |
|
| 1325 | 1340 |
// Update potentials |
| 1326 | 1341 |
void updatePotential() {
|
| 1327 | 1342 |
Cost sigma = _forward[u_in] ? |
| 1328 | 1343 |
_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
| 1329 | 1344 |
_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
| 1330 |
if (_succ_num[u_in] > _node_num / 2) {
|
|
| 1331 |
// Update in the upper subtree (which contains the root) |
|
| 1332 |
int before = _rev_thread[u_in]; |
|
| 1333 |
int after = _thread[_last_succ[u_in]]; |
|
| 1334 |
_thread[before] = after; |
|
| 1335 |
_pi[_root] -= sigma; |
|
| 1336 |
for (int u = _thread[_root]; u != _root; u = _thread[u]) {
|
|
| 1337 |
_pi[u] -= sigma; |
|
| 1338 |
} |
|
| 1339 |
_thread[before] = u_in; |
|
| 1340 |
} else {
|
|
| 1341 |
// Update in the lower subtree (which has been moved) |
|
| 1345 |
// Update potentials in the subtree, which has been moved |
|
| 1342 | 1346 |
int end = _thread[_last_succ[u_in]]; |
| 1343 | 1347 |
for (int u = u_in; u != end; u = _thread[u]) {
|
| 1344 | 1348 |
_pi[u] += sigma; |
| 1345 | 1349 |
} |
| 1346 | 1350 |
} |
| 1347 |
} |
|
| 1348 | 1351 |
|
| 1349 | 1352 |
// Execute the algorithm |
| 1350 | 1353 |
bool start(PivotRule pivot_rule) {
|
| 1351 | 1354 |
// Select the pivot rule implementation |
| 1352 | 1355 |
switch (pivot_rule) {
|
| 1353 | 1356 |
case FIRST_ELIGIBLE: |
| 1354 | 1357 |
return start<FirstEligiblePivotRule>(); |
| 1355 | 1358 |
case BEST_ELIGIBLE: |
| 1356 | 1359 |
return start<BestEligiblePivotRule>(); |
| 1357 | 1360 |
case BLOCK_SEARCH: |
| 1358 | 1361 |
return start<BlockSearchPivotRule>(); |
| 1359 | 1362 |
case CANDIDATE_LIST: |
0 comments (0 inline)