1.1 --- a/lemon/network_simplex.h Sun Oct 04 10:15:32 2009 +0200
1.2 +++ b/lemon/network_simplex.h Wed Dec 09 11:14:06 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 @@ -48,7 +50,7 @@
1.15 /// In general this class is the fastest implementation available
1.16 /// in LEMON for the minimum cost flow problem.
1.17 /// Moreover it supports both directions of the supply/demand inequality
1.18 - /// constraints. For more information see \ref SupplyType.
1.19 + /// constraints. For more information, see \ref SupplyType.
1.20 ///
1.21 /// Most of the parameters of the problem (except for the digraph)
1.22 /// can be given using separate functions, and the algorithm can be
1.23 @@ -57,16 +59,16 @@
1.24 ///
1.25 /// \tparam GR The digraph type the algorithm runs on.
1.26 /// \tparam V The value type used for flow amounts, capacity bounds
1.27 - /// and supply values in the algorithm. By default it is \c int.
1.28 + /// and supply values in the algorithm. By default, it is \c int.
1.29 /// \tparam C The value type used for costs and potentials in the
1.30 - /// algorithm. By default it is the same as \c V.
1.31 + /// algorithm. By default, it is the same as \c V.
1.32 ///
1.33 /// \warning Both value types must be signed and all input data must
1.34 /// be integer.
1.35 ///
1.36 /// \note %NetworkSimplex provides five different pivot rule
1.37 /// implementations, from which the most efficient one is used
1.38 - /// by default. For more information see \ref PivotRule.
1.39 + /// by default. For more information, see \ref PivotRule.
1.40 template <typename GR, typename V = int, typename C = V>
1.41 class NetworkSimplex
1.42 {
1.43 @@ -122,35 +124,35 @@
1.44 /// \ref NetworkSimplex provides five different pivot rule
1.45 /// implementations that significantly affect the running time
1.46 /// of the algorithm.
1.47 - /// By default \ref BLOCK_SEARCH "Block Search" is used, which
1.48 + /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
1.49 /// proved to be the most efficient and the most robust on various
1.50 /// test inputs according to our benchmark tests.
1.51 - /// However another pivot rule can be selected using the \ref run()
1.52 + /// However, another pivot rule can be selected using the \ref run()
1.53 /// function with the proper parameter.
1.54 enum PivotRule {
1.55
1.56 - /// The First Eligible pivot rule.
1.57 + /// The \e First \e Eligible pivot rule.
1.58 /// The next eligible arc is selected in a wraparound fashion
1.59 /// in every iteration.
1.60 FIRST_ELIGIBLE,
1.61
1.62 - /// The Best Eligible pivot rule.
1.63 + /// The \e Best \e Eligible pivot rule.
1.64 /// The best eligible arc is selected in every iteration.
1.65 BEST_ELIGIBLE,
1.66
1.67 - /// The Block Search pivot rule.
1.68 + /// The \e Block \e Search pivot rule.
1.69 /// A specified number of arcs are examined in every iteration
1.70 /// in a wraparound fashion and the best eligible arc is selected
1.71 /// from this block.
1.72 BLOCK_SEARCH,
1.73
1.74 - /// The Candidate List pivot rule.
1.75 + /// The \e Candidate \e List pivot rule.
1.76 /// In a major iteration a candidate list is built from eligible arcs
1.77 /// in a wraparound fashion and in the following minor iterations
1.78 /// the best eligible arc is selected from this list.
1.79 CANDIDATE_LIST,
1.80
1.81 - /// The Altering Candidate List pivot rule.
1.82 + /// The \e Altering \e Candidate \e List pivot rule.
1.83 /// It is a modified version of the Candidate List method.
1.84 /// It keeps only the several best eligible arcs from the former
1.85 /// candidate list and extends this list in every iteration.
1.86 @@ -161,8 +163,6 @@
1.87
1.88 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.89
1.90 - typedef std::vector<Arc> ArcVector;
1.91 - typedef std::vector<Node> NodeVector;
1.92 typedef std::vector<int> IntVector;
1.93 typedef std::vector<bool> BoolVector;
1.94 typedef std::vector<Value> ValueVector;
1.95 @@ -364,33 +364,32 @@
1.96 bool findEnteringArc() {
1.97 Cost c, min = 0;
1.98 int cnt = _block_size;
1.99 - int e, min_arc = _next_arc;
1.100 + int e;
1.101 for (e = _next_arc; e < _search_arc_num; ++e) {
1.102 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.103 if (c < min) {
1.104 min = c;
1.105 - min_arc = e;
1.106 + _in_arc = e;
1.107 }
1.108 if (--cnt == 0) {
1.109 - if (min < 0) break;
1.110 + if (min < 0) goto search_end;
1.111 cnt = _block_size;
1.112 }
1.113 }
1.114 - if (min == 0 || cnt > 0) {
1.115 - for (e = 0; e < _next_arc; ++e) {
1.116 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.117 - if (c < min) {
1.118 - min = c;
1.119 - min_arc = e;
1.120 - }
1.121 - if (--cnt == 0) {
1.122 - if (min < 0) break;
1.123 - cnt = _block_size;
1.124 - }
1.125 + for (e = 0; e < _next_arc; ++e) {
1.126 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.127 + if (c < min) {
1.128 + min = c;
1.129 + _in_arc = e;
1.130 + }
1.131 + if (--cnt == 0) {
1.132 + if (min < 0) goto search_end;
1.133 + cnt = _block_size;
1.134 }
1.135 }
1.136 if (min >= 0) return false;
1.137 - _in_arc = min_arc;
1.138 +
1.139 + search_end:
1.140 _next_arc = e;
1.141 return true;
1.142 }
1.143 @@ -428,7 +427,7 @@
1.144 _next_arc(0)
1.145 {
1.146 // The main parameters of the pivot rule
1.147 - const double LIST_LENGTH_FACTOR = 1.0;
1.148 + const double LIST_LENGTH_FACTOR = 0.25;
1.149 const int MIN_LIST_LENGTH = 10;
1.150 const double MINOR_LIMIT_FACTOR = 0.1;
1.151 const int MIN_MINOR_LIMIT = 3;
1.152 @@ -445,7 +444,7 @@
1.153 /// Find next entering arc
1.154 bool findEnteringArc() {
1.155 Cost min, c;
1.156 - int e, min_arc = _next_arc;
1.157 + int e;
1.158 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.159 // Minor iteration: select the best eligible arc from the
1.160 // current candidate list
1.161 @@ -456,16 +455,13 @@
1.162 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.163 if (c < min) {
1.164 min = c;
1.165 - min_arc = e;
1.166 + _in_arc = e;
1.167 }
1.168 - if (c >= 0) {
1.169 + else if (c >= 0) {
1.170 _candidates[i--] = _candidates[--_curr_length];
1.171 }
1.172 }
1.173 - if (min < 0) {
1.174 - _in_arc = min_arc;
1.175 - return true;
1.176 - }
1.177 + if (min < 0) return true;
1.178 }
1.179
1.180 // Major iteration: build a new candidate list
1.181 @@ -477,27 +473,26 @@
1.182 _candidates[_curr_length++] = e;
1.183 if (c < min) {
1.184 min = c;
1.185 - min_arc = e;
1.186 + _in_arc = e;
1.187 }
1.188 - if (_curr_length == _list_length) break;
1.189 + if (_curr_length == _list_length) goto search_end;
1.190 }
1.191 }
1.192 - if (_curr_length < _list_length) {
1.193 - for (e = 0; e < _next_arc; ++e) {
1.194 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.195 - if (c < 0) {
1.196 - _candidates[_curr_length++] = e;
1.197 - if (c < min) {
1.198 - min = c;
1.199 - min_arc = e;
1.200 - }
1.201 - if (_curr_length == _list_length) break;
1.202 + for (e = 0; e < _next_arc; ++e) {
1.203 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.204 + if (c < 0) {
1.205 + _candidates[_curr_length++] = e;
1.206 + if (c < min) {
1.207 + min = c;
1.208 + _in_arc = e;
1.209 }
1.210 + if (_curr_length == _list_length) goto search_end;
1.211 }
1.212 }
1.213 if (_curr_length == 0) return false;
1.214 +
1.215 + search_end:
1.216 _minor_count = 1;
1.217 - _in_arc = min_arc;
1.218 _next_arc = e;
1.219 return true;
1.220 }
1.221 @@ -549,7 +544,7 @@
1.222 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.223 {
1.224 // The main parameters of the pivot rule
1.225 - const double BLOCK_SIZE_FACTOR = 1.5;
1.226 + const double BLOCK_SIZE_FACTOR = 1.0;
1.227 const int MIN_BLOCK_SIZE = 10;
1.228 const double HEAD_LENGTH_FACTOR = 0.1;
1.229 const int MIN_HEAD_LENGTH = 3;
1.230 @@ -578,39 +573,35 @@
1.231
1.232 // Extend the list
1.233 int cnt = _block_size;
1.234 - int last_arc = 0;
1.235 int limit = _head_length;
1.236
1.237 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.238 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.239 _cand_cost[e] = _state[e] *
1.240 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.241 if (_cand_cost[e] < 0) {
1.242 _candidates[_curr_length++] = e;
1.243 - last_arc = e;
1.244 }
1.245 if (--cnt == 0) {
1.246 - if (_curr_length > limit) break;
1.247 + if (_curr_length > limit) goto search_end;
1.248 limit = 0;
1.249 cnt = _block_size;
1.250 }
1.251 }
1.252 - if (_curr_length <= limit) {
1.253 - for (int e = 0; e < _next_arc; ++e) {
1.254 - _cand_cost[e] = _state[e] *
1.255 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.256 - if (_cand_cost[e] < 0) {
1.257 - _candidates[_curr_length++] = e;
1.258 - last_arc = e;
1.259 - }
1.260 - if (--cnt == 0) {
1.261 - if (_curr_length > limit) break;
1.262 - limit = 0;
1.263 - cnt = _block_size;
1.264 - }
1.265 + for (e = 0; e < _next_arc; ++e) {
1.266 + _cand_cost[e] = _state[e] *
1.267 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.268 + if (_cand_cost[e] < 0) {
1.269 + _candidates[_curr_length++] = e;
1.270 + }
1.271 + if (--cnt == 0) {
1.272 + if (_curr_length > limit) goto search_end;
1.273 + limit = 0;
1.274 + cnt = _block_size;
1.275 }
1.276 }
1.277 if (_curr_length == 0) return false;
1.278 - _next_arc = last_arc + 1;
1.279 +
1.280 + search_end:
1.281
1.282 // Make heap of the candidate list (approximating a partial sort)
1.283 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.284 @@ -618,6 +609,7 @@
1.285
1.286 // Pop the first element of the heap
1.287 _in_arc = _candidates[0];
1.288 + _next_arc = e;
1.289 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.290 _sort_func );
1.291 _curr_length = std::min(_head_length, _curr_length - 1);
1.292 @@ -633,7 +625,11 @@
1.293 /// The constructor of the class.
1.294 ///
1.295 /// \param graph The digraph the algorithm runs on.
1.296 - NetworkSimplex(const GR& graph) :
1.297 + /// \param arc_mixing Indicate if the arcs have to be stored in a
1.298 + /// mixed order in the internal data structure.
1.299 + /// In special cases, it could lead to better overall performance,
1.300 + /// but it is usually slower. Therefore it is disabled by default.
1.301 + NetworkSimplex(const GR& graph, bool arc_mixing = false) :
1.302 _graph(graph), _node_id(graph), _arc_id(graph),
1.303 INF(std::numeric_limits<Value>::has_infinity ?
1.304 std::numeric_limits<Value>::infinity() :
1.305 @@ -671,31 +667,33 @@
1.306 _last_succ.resize(all_node_num);
1.307 _state.resize(max_arc_num);
1.308
1.309 - // Copy the graph (store the arcs in a mixed order)
1.310 + // Copy the graph
1.311 int i = 0;
1.312 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.313 _node_id[n] = i;
1.314 }
1.315 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.316 - i = 0;
1.317 - for (ArcIt a(_graph); a != INVALID; ++a) {
1.318 - _arc_id[a] = i;
1.319 - _source[i] = _node_id[_graph.source(a)];
1.320 - _target[i] = _node_id[_graph.target(a)];
1.321 - if ((i += k) >= _arc_num) i = (i % k) + 1;
1.322 + if (arc_mixing) {
1.323 + // Store the arcs in a mixed order
1.324 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.325 + int i = 0, j = 0;
1.326 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.327 + _arc_id[a] = i;
1.328 + _source[i] = _node_id[_graph.source(a)];
1.329 + _target[i] = _node_id[_graph.target(a)];
1.330 + if ((i += k) >= _arc_num) i = ++j;
1.331 + }
1.332 + } else {
1.333 + // Store the arcs in the original order
1.334 + int i = 0;
1.335 + for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
1.336 + _arc_id[a] = i;
1.337 + _source[i] = _node_id[_graph.source(a)];
1.338 + _target[i] = _node_id[_graph.target(a)];
1.339 + }
1.340 }
1.341
1.342 - // Initialize maps
1.343 - for (int i = 0; i != _node_num; ++i) {
1.344 - _supply[i] = 0;
1.345 - }
1.346 - for (int i = 0; i != _arc_num; ++i) {
1.347 - _lower[i] = 0;
1.348 - _upper[i] = INF;
1.349 - _cost[i] = 1;
1.350 - }
1.351 - _have_lower = false;
1.352 - _stype = GEQ;
1.353 + // Reset parameters
1.354 + reset();
1.355 }
1.356
1.357 /// \name Parameters
1.358 @@ -768,7 +766,6 @@
1.359 /// This function sets the supply values of the nodes.
1.360 /// If neither this function nor \ref stSupply() is used before
1.361 /// calling \ref run(), the supply of each node will be set to zero.
1.362 - /// (It makes sense only if non-zero lower bounds are given.)
1.363 ///
1.364 /// \param map A node map storing the supply values.
1.365 /// Its \c Value type must be convertible to the \c Value type
1.366 @@ -789,7 +786,6 @@
1.367 /// and the required flow value.
1.368 /// If neither this function nor \ref supplyMap() is used before
1.369 /// calling \ref run(), the supply of each node will be set to zero.
1.370 - /// (It makes sense only if non-zero lower bounds are given.)
1.371 ///
1.372 /// Using this function has the same effect as using \ref supplyMap()
1.373 /// with such a map in which \c k is assigned to \c s, \c -k is
1.374 @@ -816,7 +812,7 @@
1.375 /// If it is not used before calling \ref run(), the \ref GEQ supply
1.376 /// type will be used.
1.377 ///
1.378 - /// For more information see \ref SupplyType.
1.379 + /// For more information, see \ref SupplyType.
1.380 ///
1.381 /// \return <tt>(*this)</tt>
1.382 NetworkSimplex& supplyType(SupplyType supply_type) {
1.383 @@ -848,11 +844,11 @@
1.384 /// that have been given are kept for the next call, unless
1.385 /// \ref reset() is called, thus only the modified parameters
1.386 /// have to be set again. See \ref reset() for examples.
1.387 - /// However the underlying digraph must not be modified after this
1.388 + /// However, the underlying digraph must not be modified after this
1.389 /// class have been constructed, since it copies and extends the graph.
1.390 ///
1.391 /// \param pivot_rule The pivot rule that will be used during the
1.392 - /// algorithm. For more information see \ref PivotRule.
1.393 + /// algorithm. For more information, see \ref PivotRule.
1.394 ///
1.395 /// \return \c INFEASIBLE if no feasible flow exists,
1.396 /// \n \c OPTIMAL if the problem has optimal solution
1.397 @@ -877,7 +873,7 @@
1.398 /// It is useful for multiple run() calls. If this function is not
1.399 /// used, all the parameters given before are kept for the next
1.400 /// \ref run() call.
1.401 - /// However the underlying digraph must not be modified after this
1.402 + /// However, the underlying digraph must not be modified after this
1.403 /// class have been constructed, since it copies and extends the graph.
1.404 ///
1.405 /// For example,