1.1 --- a/lemon/network_simplex.h Thu Nov 05 10:01:02 2009 +0100
1.2 +++ b/lemon/network_simplex.h Thu Nov 05 10:23:16 2009 +0100
1.3 @@ -40,7 +40,9 @@
1.4 /// for finding a \ref min_cost_flow "minimum cost flow".
1.5 ///
1.6 /// \ref NetworkSimplex implements the primal Network Simplex algorithm
1.7 - /// for finding a \ref min_cost_flow "minimum cost flow".
1.8 + /// for finding a \ref min_cost_flow "minimum cost flow"
1.9 + /// \ref amo93networkflows, \ref dantzig63linearprog,
1.10 + /// \ref kellyoneill91netsimplex.
1.11 /// This algorithm is a specialized version of the linear programming
1.12 /// simplex method directly for the minimum cost flow problem.
1.13 /// It is one of the most efficient solution methods.
1.14 @@ -161,8 +163,6 @@
1.15
1.16 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.17
1.18 - typedef std::vector<Arc> ArcVector;
1.19 - typedef std::vector<Node> NodeVector;
1.20 typedef std::vector<int> IntVector;
1.21 typedef std::vector<bool> BoolVector;
1.22 typedef std::vector<Value> ValueVector;
1.23 @@ -364,33 +364,32 @@
1.24 bool findEnteringArc() {
1.25 Cost c, min = 0;
1.26 int cnt = _block_size;
1.27 - int e, min_arc = _next_arc;
1.28 + int e;
1.29 for (e = _next_arc; e < _search_arc_num; ++e) {
1.30 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.31 if (c < min) {
1.32 min = c;
1.33 - min_arc = e;
1.34 + _in_arc = e;
1.35 }
1.36 if (--cnt == 0) {
1.37 - if (min < 0) break;
1.38 + if (min < 0) goto search_end;
1.39 cnt = _block_size;
1.40 }
1.41 }
1.42 - if (min == 0 || cnt > 0) {
1.43 - for (e = 0; e < _next_arc; ++e) {
1.44 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.45 - if (c < min) {
1.46 - min = c;
1.47 - min_arc = e;
1.48 - }
1.49 - if (--cnt == 0) {
1.50 - if (min < 0) break;
1.51 - cnt = _block_size;
1.52 - }
1.53 + for (e = 0; e < _next_arc; ++e) {
1.54 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.55 + if (c < min) {
1.56 + min = c;
1.57 + _in_arc = e;
1.58 + }
1.59 + if (--cnt == 0) {
1.60 + if (min < 0) goto search_end;
1.61 + cnt = _block_size;
1.62 }
1.63 }
1.64 if (min >= 0) return false;
1.65 - _in_arc = min_arc;
1.66 +
1.67 + search_end:
1.68 _next_arc = e;
1.69 return true;
1.70 }
1.71 @@ -428,7 +427,7 @@
1.72 _next_arc(0)
1.73 {
1.74 // The main parameters of the pivot rule
1.75 - const double LIST_LENGTH_FACTOR = 1.0;
1.76 + const double LIST_LENGTH_FACTOR = 0.25;
1.77 const int MIN_LIST_LENGTH = 10;
1.78 const double MINOR_LIMIT_FACTOR = 0.1;
1.79 const int MIN_MINOR_LIMIT = 3;
1.80 @@ -445,7 +444,7 @@
1.81 /// Find next entering arc
1.82 bool findEnteringArc() {
1.83 Cost min, c;
1.84 - int e, min_arc = _next_arc;
1.85 + int e;
1.86 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.87 // Minor iteration: select the best eligible arc from the
1.88 // current candidate list
1.89 @@ -456,16 +455,13 @@
1.90 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.91 if (c < min) {
1.92 min = c;
1.93 - min_arc = e;
1.94 + _in_arc = e;
1.95 }
1.96 - if (c >= 0) {
1.97 + else if (c >= 0) {
1.98 _candidates[i--] = _candidates[--_curr_length];
1.99 }
1.100 }
1.101 - if (min < 0) {
1.102 - _in_arc = min_arc;
1.103 - return true;
1.104 - }
1.105 + if (min < 0) return true;
1.106 }
1.107
1.108 // Major iteration: build a new candidate list
1.109 @@ -477,27 +473,26 @@
1.110 _candidates[_curr_length++] = e;
1.111 if (c < min) {
1.112 min = c;
1.113 - min_arc = e;
1.114 + _in_arc = e;
1.115 }
1.116 - if (_curr_length == _list_length) break;
1.117 + if (_curr_length == _list_length) goto search_end;
1.118 }
1.119 }
1.120 - if (_curr_length < _list_length) {
1.121 - for (e = 0; e < _next_arc; ++e) {
1.122 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.123 - if (c < 0) {
1.124 - _candidates[_curr_length++] = e;
1.125 - if (c < min) {
1.126 - min = c;
1.127 - min_arc = e;
1.128 - }
1.129 - if (_curr_length == _list_length) break;
1.130 + for (e = 0; e < _next_arc; ++e) {
1.131 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.132 + if (c < 0) {
1.133 + _candidates[_curr_length++] = e;
1.134 + if (c < min) {
1.135 + min = c;
1.136 + _in_arc = e;
1.137 }
1.138 + if (_curr_length == _list_length) goto search_end;
1.139 }
1.140 }
1.141 if (_curr_length == 0) return false;
1.142 +
1.143 + search_end:
1.144 _minor_count = 1;
1.145 - _in_arc = min_arc;
1.146 _next_arc = e;
1.147 return true;
1.148 }
1.149 @@ -549,7 +544,7 @@
1.150 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.151 {
1.152 // The main parameters of the pivot rule
1.153 - const double BLOCK_SIZE_FACTOR = 1.5;
1.154 + const double BLOCK_SIZE_FACTOR = 1.0;
1.155 const int MIN_BLOCK_SIZE = 10;
1.156 const double HEAD_LENGTH_FACTOR = 0.1;
1.157 const int MIN_HEAD_LENGTH = 3;
1.158 @@ -578,39 +573,35 @@
1.159
1.160 // Extend the list
1.161 int cnt = _block_size;
1.162 - int last_arc = 0;
1.163 int limit = _head_length;
1.164
1.165 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.166 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.167 _cand_cost[e] = _state[e] *
1.168 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.169 if (_cand_cost[e] < 0) {
1.170 _candidates[_curr_length++] = e;
1.171 - last_arc = e;
1.172 }
1.173 if (--cnt == 0) {
1.174 - if (_curr_length > limit) break;
1.175 + if (_curr_length > limit) goto search_end;
1.176 limit = 0;
1.177 cnt = _block_size;
1.178 }
1.179 }
1.180 - if (_curr_length <= limit) {
1.181 - for (int e = 0; e < _next_arc; ++e) {
1.182 - _cand_cost[e] = _state[e] *
1.183 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.184 - if (_cand_cost[e] < 0) {
1.185 - _candidates[_curr_length++] = e;
1.186 - last_arc = e;
1.187 - }
1.188 - if (--cnt == 0) {
1.189 - if (_curr_length > limit) break;
1.190 - limit = 0;
1.191 - cnt = _block_size;
1.192 - }
1.193 + for (e = 0; e < _next_arc; ++e) {
1.194 + _cand_cost[e] = _state[e] *
1.195 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.196 + if (_cand_cost[e] < 0) {
1.197 + _candidates[_curr_length++] = e;
1.198 + }
1.199 + if (--cnt == 0) {
1.200 + if (_curr_length > limit) goto search_end;
1.201 + limit = 0;
1.202 + cnt = _block_size;
1.203 }
1.204 }
1.205 if (_curr_length == 0) return false;
1.206 - _next_arc = last_arc + 1;
1.207 +
1.208 + search_end:
1.209
1.210 // Make heap of the candidate list (approximating a partial sort)
1.211 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.212 @@ -618,6 +609,7 @@
1.213
1.214 // Pop the first element of the heap
1.215 _in_arc = _candidates[0];
1.216 + _next_arc = e;
1.217 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.218 _sort_func );
1.219 _curr_length = std::min(_head_length, _curr_length - 1);
1.220 @@ -633,7 +625,11 @@
1.221 /// The constructor of the class.
1.222 ///
1.223 /// \param graph The digraph the algorithm runs on.
1.224 - NetworkSimplex(const GR& graph) :
1.225 + /// \param arc_mixing Indicate if the arcs have to be stored in a
1.226 + /// mixed order in the internal data structure.
1.227 + /// In special cases, it could lead to better overall performance,
1.228 + /// but it is usually slower. Therefore it is disabled by default.
1.229 + NetworkSimplex(const GR& graph, bool arc_mixing = false) :
1.230 _graph(graph), _node_id(graph), _arc_id(graph),
1.231 INF(std::numeric_limits<Value>::has_infinity ?
1.232 std::numeric_limits<Value>::infinity() :
1.233 @@ -671,31 +667,33 @@
1.234 _last_succ.resize(all_node_num);
1.235 _state.resize(max_arc_num);
1.236
1.237 - // Copy the graph (store the arcs in a mixed order)
1.238 + // Copy the graph
1.239 int i = 0;
1.240 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.241 _node_id[n] = i;
1.242 }
1.243 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.244 - i = 0;
1.245 - for (ArcIt a(_graph); a != INVALID; ++a) {
1.246 - _arc_id[a] = i;
1.247 - _source[i] = _node_id[_graph.source(a)];
1.248 - _target[i] = _node_id[_graph.target(a)];
1.249 - if ((i += k) >= _arc_num) i = (i % k) + 1;
1.250 + if (arc_mixing) {
1.251 + // Store the arcs in a mixed order
1.252 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.253 + int i = 0, j = 0;
1.254 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.255 + _arc_id[a] = i;
1.256 + _source[i] = _node_id[_graph.source(a)];
1.257 + _target[i] = _node_id[_graph.target(a)];
1.258 + if ((i += k) >= _arc_num) i = ++j;
1.259 + }
1.260 + } else {
1.261 + // Store the arcs in the original order
1.262 + int i = 0;
1.263 + for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
1.264 + _arc_id[a] = i;
1.265 + _source[i] = _node_id[_graph.source(a)];
1.266 + _target[i] = _node_id[_graph.target(a)];
1.267 + }
1.268 }
1.269
1.270 - // Initialize maps
1.271 - for (int i = 0; i != _node_num; ++i) {
1.272 - _supply[i] = 0;
1.273 - }
1.274 - for (int i = 0; i != _arc_num; ++i) {
1.275 - _lower[i] = 0;
1.276 - _upper[i] = INF;
1.277 - _cost[i] = 1;
1.278 - }
1.279 - _have_lower = false;
1.280 - _stype = GEQ;
1.281 + // Reset parameters
1.282 + reset();
1.283 }
1.284
1.285 /// \name Parameters
1.286 @@ -768,7 +766,6 @@
1.287 /// This function sets the supply values of the nodes.
1.288 /// If neither this function nor \ref stSupply() is used before
1.289 /// calling \ref run(), the supply of each node will be set to zero.
1.290 - /// (It makes sense only if non-zero lower bounds are given.)
1.291 ///
1.292 /// \param map A node map storing the supply values.
1.293 /// Its \c Value type must be convertible to the \c Value type
1.294 @@ -789,7 +786,6 @@
1.295 /// and the required flow value.
1.296 /// If neither this function nor \ref supplyMap() is used before
1.297 /// calling \ref run(), the supply of each node will be set to zero.
1.298 - /// (It makes sense only if non-zero lower bounds are given.)
1.299 ///
1.300 /// Using this function has the same effect as using \ref supplyMap()
1.301 /// with such a map in which \c k is assigned to \c s, \c -k is