52 /// \tparam F The value type used for flow amounts, capacity bounds |
52 /// \tparam F The value type used for flow amounts, capacity bounds |
53 /// and supply values in the algorithm. By default it is \c int. |
53 /// and supply values in the algorithm. By default it is \c int. |
54 /// \tparam C The value type used for costs and potentials in the |
54 /// \tparam C The value type used for costs and potentials in the |
55 /// algorithm. By default it is the same as \c F. |
55 /// algorithm. By default it is the same as \c F. |
56 /// |
56 /// |
57 /// \warning Both value types must be signed integer types. |
57 /// \warning Both value types must be signed and all input data must |
|
58 /// be integer. |
58 /// |
59 /// |
59 /// \note %NetworkSimplex provides five different pivot rule |
60 /// \note %NetworkSimplex provides five different pivot rule |
60 /// implementations. For more information see \ref PivotRule. |
61 /// implementations. For more information see \ref PivotRule. |
61 template <typename GR, typename F = int, typename C = F> |
62 template <typename GR, typename F = int, typename C = F> |
62 class NetworkSimplex |
63 class NetworkSimplex |
1042 _arc_ref[i] = e; |
1043 _arc_ref[i] = e; |
1043 if ((i += k) >= _arc_num) i = (i % k) + 1; |
1044 if ((i += k) >= _arc_num) i = (i % k) + 1; |
1044 } |
1045 } |
1045 |
1046 |
1046 // Initialize arc maps |
1047 // Initialize arc maps |
1047 Flow max_cap = std::numeric_limits<Flow>::max(); |
1048 Flow inf_cap = |
1048 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
1049 std::numeric_limits<Flow>::has_infinity ? |
|
1050 std::numeric_limits<Flow>::infinity() : |
|
1051 std::numeric_limits<Flow>::max(); |
1049 if (_pupper && _pcost) { |
1052 if (_pupper && _pcost) { |
1050 for (int i = 0; i != _arc_num; ++i) { |
1053 for (int i = 0; i != _arc_num; ++i) { |
1051 Arc e = _arc_ref[i]; |
1054 Arc e = _arc_ref[i]; |
1052 _source[i] = _node_id[_graph.source(e)]; |
1055 _source[i] = _node_id[_graph.source(e)]; |
1053 _target[i] = _node_id[_graph.target(e)]; |
1056 _target[i] = _node_id[_graph.target(e)]; |
1067 if (_pupper) { |
1070 if (_pupper) { |
1068 for (int i = 0; i != _arc_num; ++i) |
1071 for (int i = 0; i != _arc_num; ++i) |
1069 _cap[i] = (*_pupper)[_arc_ref[i]]; |
1072 _cap[i] = (*_pupper)[_arc_ref[i]]; |
1070 } else { |
1073 } else { |
1071 for (int i = 0; i != _arc_num; ++i) |
1074 for (int i = 0; i != _arc_num; ++i) |
1072 _cap[i] = max_cap; |
1075 _cap[i] = inf_cap; |
1073 } |
1076 } |
1074 if (_pcost) { |
1077 if (_pcost) { |
1075 for (int i = 0; i != _arc_num; ++i) |
1078 for (int i = 0; i != _arc_num; ++i) |
1076 _cost[i] = (*_pcost)[_arc_ref[i]]; |
1079 _cost[i] = (*_pcost)[_arc_ref[i]]; |
1077 } else { |
1080 } else { |
1078 for (int i = 0; i != _arc_num; ++i) |
1081 for (int i = 0; i != _arc_num; ++i) |
1079 _cost[i] = 1; |
1082 _cost[i] = 1; |
1080 } |
1083 } |
|
1084 } |
|
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; |
1081 } |
1096 } |
1082 |
1097 |
1083 // Remove non-zero lower bounds |
1098 // Remove non-zero lower bounds |
1084 if (_plower) { |
1099 if (_plower) { |
1085 for (int i = 0; i != _arc_num; ++i) { |
1100 for (int i = 0; i != _arc_num; ++i) { |
1098 _rev_thread[u + 1] = u; |
1113 _rev_thread[u + 1] = u; |
1099 _succ_num[u] = 1; |
1114 _succ_num[u] = 1; |
1100 _last_succ[u] = u; |
1115 _last_succ[u] = u; |
1101 _parent[u] = _root; |
1116 _parent[u] = _root; |
1102 _pred[u] = e; |
1117 _pred[u] = e; |
1103 _cost[e] = max_cost; |
1118 _cost[e] = art_cost; |
1104 _cap[e] = max_cap; |
1119 _cap[e] = inf_cap; |
1105 _state[e] = STATE_TREE; |
1120 _state[e] = STATE_TREE; |
1106 if (_supply[u] >= 0) { |
1121 if (_supply[u] >= 0) { |
1107 _flow[e] = _supply[u]; |
1122 _flow[e] = _supply[u]; |
1108 _forward[u] = true; |
1123 _forward[u] = true; |
1109 _pi[u] = -max_cost; |
1124 _pi[u] = -art_cost; |
1110 } else { |
1125 } else { |
1111 _flow[e] = -_supply[u]; |
1126 _flow[e] = -_supply[u]; |
1112 _forward[u] = false; |
1127 _forward[u] = false; |
1113 _pi[u] = max_cost; |
1128 _pi[u] = art_cost; |
1114 } |
1129 } |
1115 } |
1130 } |
1116 |
1131 |
1117 return true; |
1132 return true; |
1118 } |
1133 } |
1325 // Update potentials |
1340 // Update potentials |
1326 void updatePotential() { |
1341 void updatePotential() { |
1327 Cost sigma = _forward[u_in] ? |
1342 Cost sigma = _forward[u_in] ? |
1328 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1343 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1329 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1344 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1330 if (_succ_num[u_in] > _node_num / 2) { |
1345 // Update potentials in the subtree, which has been moved |
1331 // Update in the upper subtree (which contains the root) |
1346 int end = _thread[_last_succ[u_in]]; |
1332 int before = _rev_thread[u_in]; |
1347 for (int u = u_in; u != end; u = _thread[u]) { |
1333 int after = _thread[_last_succ[u_in]]; |
1348 _pi[u] += sigma; |
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) |
|
1342 int end = _thread[_last_succ[u_in]]; |
|
1343 for (int u = u_in; u != end; u = _thread[u]) { |
|
1344 _pi[u] += sigma; |
|
1345 } |
|
1346 } |
1349 } |
1347 } |
1350 } |
1348 |
1351 |
1349 // Execute the algorithm |
1352 // Execute the algorithm |
1350 bool start(PivotRule pivot_rule) { |
1353 bool start(PivotRule pivot_rule) { |