1.1 --- a/lemon/network_simplex.h Fri Aug 09 11:07:27 2013 +0200
1.2 +++ b/lemon/network_simplex.h Sun Aug 11 15:28:12 2013 +0200
1.3 @@ -2,7 +2,7 @@
1.4 *
1.5 * This file is a part of LEMON, a generic C++ optimization library.
1.6 *
1.7 - * Copyright (C) 2003-2009
1.8 + * Copyright (C) 2003-2010
1.9 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 *
1.12 @@ -40,15 +40,17 @@
1.13 /// for finding a \ref min_cost_flow "minimum cost flow".
1.14 ///
1.15 /// \ref NetworkSimplex implements the primal Network Simplex algorithm
1.16 - /// for finding a \ref min_cost_flow "minimum cost flow".
1.17 - /// This algorithm is a specialized version of the linear programming
1.18 - /// simplex method directly for the minimum cost flow problem.
1.19 - /// It is one of the most efficient solution methods.
1.20 + /// for finding a \ref min_cost_flow "minimum cost flow"
1.21 + /// \ref amo93networkflows, \ref dantzig63linearprog,
1.22 + /// \ref kellyoneill91netsimplex.
1.23 + /// This algorithm is a highly efficient specialized version of the
1.24 + /// linear programming simplex method directly for the minimum cost
1.25 + /// flow problem.
1.26 ///
1.27 - /// In general this class is the fastest implementation available
1.28 - /// in LEMON for the minimum cost flow problem.
1.29 - /// Moreover it supports both directions of the supply/demand inequality
1.30 - /// constraints. For more information see \ref SupplyType.
1.31 + /// In general, %NetworkSimplex is the fastest implementation available
1.32 + /// in LEMON for this problem.
1.33 + /// Moreover, it supports both directions of the supply/demand inequality
1.34 + /// constraints. For more information, see \ref SupplyType.
1.35 ///
1.36 /// Most of the parameters of the problem (except for the digraph)
1.37 /// can be given using separate functions, and the algorithm can be
1.38 @@ -56,17 +58,17 @@
1.39 /// specified, then default values will be used.
1.40 ///
1.41 /// \tparam GR The digraph type the algorithm runs on.
1.42 - /// \tparam V The value type used for flow amounts, capacity bounds
1.43 - /// and supply values in the algorithm. By default it is \c int.
1.44 - /// \tparam C The value type used for costs and potentials in the
1.45 - /// algorithm. By default it is the same as \c V.
1.46 + /// \tparam V The number type used for flow amounts, capacity bounds
1.47 + /// and supply values in the algorithm. By default, it is \c int.
1.48 + /// \tparam C The number type used for costs and potentials in the
1.49 + /// algorithm. By default, it is the same as \c V.
1.50 ///
1.51 - /// \warning Both value types must be signed and all input data must
1.52 + /// \warning Both number types must be signed and all input data must
1.53 /// be integer.
1.54 ///
1.55 /// \note %NetworkSimplex provides five different pivot rule
1.56 /// implementations, from which the most efficient one is used
1.57 - /// by default. For more information see \ref PivotRule.
1.58 + /// by default. For more information, see \ref PivotRule.
1.59 template <typename GR, typename V = int, typename C = V>
1.60 class NetworkSimplex
1.61 {
1.62 @@ -95,7 +97,7 @@
1.63 /// infinite upper bound.
1.64 UNBOUNDED
1.65 };
1.66 -
1.67 +
1.68 /// \brief Constants for selecting the type of the supply constraints.
1.69 ///
1.70 /// Enum type containing constants for selecting the supply type,
1.71 @@ -113,7 +115,7 @@
1.72 /// supply/demand constraints in the definition of the problem.
1.73 LEQ
1.74 };
1.75 -
1.76 +
1.77 /// \brief Constants for selecting the pivot rule.
1.78 ///
1.79 /// Enum type containing constants for selecting the pivot rule for
1.80 @@ -122,59 +124,62 @@
1.81 /// \ref NetworkSimplex provides five different pivot rule
1.82 /// implementations that significantly affect the running time
1.83 /// of the algorithm.
1.84 - /// By default \ref BLOCK_SEARCH "Block Search" is used, which
1.85 + /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
1.86 /// proved to be the most efficient and the most robust on various
1.87 - /// test inputs according to our benchmark tests.
1.88 - /// However another pivot rule can be selected using the \ref run()
1.89 + /// test inputs.
1.90 + /// However, another pivot rule can be selected using the \ref run()
1.91 /// function with the proper parameter.
1.92 enum PivotRule {
1.93
1.94 - /// The First Eligible pivot rule.
1.95 + /// The \e First \e Eligible pivot rule.
1.96 /// The next eligible arc is selected in a wraparound fashion
1.97 /// in every iteration.
1.98 FIRST_ELIGIBLE,
1.99
1.100 - /// The Best Eligible pivot rule.
1.101 + /// The \e Best \e Eligible pivot rule.
1.102 /// The best eligible arc is selected in every iteration.
1.103 BEST_ELIGIBLE,
1.104
1.105 - /// The Block Search pivot rule.
1.106 + /// The \e Block \e Search pivot rule.
1.107 /// A specified number of arcs are examined in every iteration
1.108 /// in a wraparound fashion and the best eligible arc is selected
1.109 /// from this block.
1.110 BLOCK_SEARCH,
1.111
1.112 - /// The Candidate List pivot rule.
1.113 + /// The \e Candidate \e List pivot rule.
1.114 /// In a major iteration a candidate list is built from eligible arcs
1.115 /// in a wraparound fashion and in the following minor iterations
1.116 /// the best eligible arc is selected from this list.
1.117 CANDIDATE_LIST,
1.118
1.119 - /// The Altering Candidate List pivot rule.
1.120 + /// The \e Altering \e Candidate \e List pivot rule.
1.121 /// It is a modified version of the Candidate List method.
1.122 /// It keeps only the several best eligible arcs from the former
1.123 /// candidate list and extends this list in every iteration.
1.124 ALTERING_LIST
1.125 };
1.126 -
1.127 +
1.128 private:
1.129
1.130 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.131
1.132 - typedef std::vector<Arc> ArcVector;
1.133 - typedef std::vector<Node> NodeVector;
1.134 typedef std::vector<int> IntVector;
1.135 - typedef std::vector<bool> BoolVector;
1.136 typedef std::vector<Value> ValueVector;
1.137 typedef std::vector<Cost> CostVector;
1.138 + typedef std::vector<char> BoolVector;
1.139 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.140
1.141 // State constants for arcs
1.142 - enum ArcStateEnum {
1.143 + enum ArcState {
1.144 STATE_UPPER = -1,
1.145 STATE_TREE = 0,
1.146 STATE_LOWER = 1
1.147 };
1.148
1.149 + typedef std::vector<signed char> StateVector;
1.150 + // Note: vector<signed char> is used instead of vector<ArcState> for
1.151 + // efficiency reasons
1.152 +
1.153 private:
1.154
1.155 // Data related to the underlying digraph
1.156 @@ -194,6 +199,7 @@
1.157 IntArcMap _arc_id;
1.158 IntVector _source;
1.159 IntVector _target;
1.160 + bool _arc_mixing;
1.161
1.162 // Node and arc data
1.163 ValueVector _lower;
1.164 @@ -213,7 +219,7 @@
1.165 IntVector _last_succ;
1.166 IntVector _dirty_revs;
1.167 BoolVector _forward;
1.168 - IntVector _state;
1.169 + StateVector _state;
1.170 int _root;
1.171
1.172 // Temporary data used in the current pivot iteration
1.173 @@ -222,8 +228,10 @@
1.174 int stem, par_stem, new_stem;
1.175 Value delta;
1.176
1.177 + const Value MAX;
1.178 +
1.179 public:
1.180 -
1.181 +
1.182 /// \brief Constant for infinite upper bounds (capacities).
1.183 ///
1.184 /// Constant for infinite upper bounds (capacities).
1.185 @@ -242,7 +250,7 @@
1.186 const IntVector &_source;
1.187 const IntVector &_target;
1.188 const CostVector &_cost;
1.189 - const IntVector &_state;
1.190 + const StateVector &_state;
1.191 const CostVector &_pi;
1.192 int &_in_arc;
1.193 int _search_arc_num;
1.194 @@ -263,7 +271,7 @@
1.195 // Find next entering arc
1.196 bool findEnteringArc() {
1.197 Cost c;
1.198 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.199 + for (int e = _next_arc; e != _search_arc_num; ++e) {
1.200 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.201 if (c < 0) {
1.202 _in_arc = e;
1.203 @@ -271,7 +279,7 @@
1.204 return true;
1.205 }
1.206 }
1.207 - for (int e = 0; e < _next_arc; ++e) {
1.208 + for (int e = 0; e != _next_arc; ++e) {
1.209 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.210 if (c < 0) {
1.211 _in_arc = e;
1.212 @@ -294,7 +302,7 @@
1.213 const IntVector &_source;
1.214 const IntVector &_target;
1.215 const CostVector &_cost;
1.216 - const IntVector &_state;
1.217 + const StateVector &_state;
1.218 const CostVector &_pi;
1.219 int &_in_arc;
1.220 int _search_arc_num;
1.221 @@ -311,7 +319,7 @@
1.222 // Find next entering arc
1.223 bool findEnteringArc() {
1.224 Cost c, min = 0;
1.225 - for (int e = 0; e < _search_arc_num; ++e) {
1.226 + for (int e = 0; e != _search_arc_num; ++e) {
1.227 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.228 if (c < min) {
1.229 min = c;
1.230 @@ -333,7 +341,7 @@
1.231 const IntVector &_source;
1.232 const IntVector &_target;
1.233 const CostVector &_cost;
1.234 - const IntVector &_state;
1.235 + const StateVector &_state;
1.236 const CostVector &_pi;
1.237 int &_in_arc;
1.238 int _search_arc_num;
1.239 @@ -352,7 +360,7 @@
1.240 _next_arc(0)
1.241 {
1.242 // The main parameters of the pivot rule
1.243 - const double BLOCK_SIZE_FACTOR = 0.5;
1.244 + const double BLOCK_SIZE_FACTOR = 1.0;
1.245 const int MIN_BLOCK_SIZE = 10;
1.246
1.247 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.248 @@ -364,33 +372,32 @@
1.249 bool findEnteringArc() {
1.250 Cost c, min = 0;
1.251 int cnt = _block_size;
1.252 - int e, min_arc = _next_arc;
1.253 - for (e = _next_arc; e < _search_arc_num; ++e) {
1.254 + int e;
1.255 + for (e = _next_arc; e != _search_arc_num; ++e) {
1.256 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.257 if (c < min) {
1.258 min = c;
1.259 - min_arc = e;
1.260 + _in_arc = e;
1.261 }
1.262 if (--cnt == 0) {
1.263 - if (min < 0) break;
1.264 + if (min < 0) goto search_end;
1.265 cnt = _block_size;
1.266 }
1.267 }
1.268 - if (min == 0 || cnt > 0) {
1.269 - for (e = 0; e < _next_arc; ++e) {
1.270 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.271 - if (c < min) {
1.272 - min = c;
1.273 - min_arc = e;
1.274 - }
1.275 - if (--cnt == 0) {
1.276 - if (min < 0) break;
1.277 - cnt = _block_size;
1.278 - }
1.279 + for (e = 0; e != _next_arc; ++e) {
1.280 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.281 + if (c < min) {
1.282 + min = c;
1.283 + _in_arc = e;
1.284 + }
1.285 + if (--cnt == 0) {
1.286 + if (min < 0) goto search_end;
1.287 + cnt = _block_size;
1.288 }
1.289 }
1.290 if (min >= 0) return false;
1.291 - _in_arc = min_arc;
1.292 +
1.293 + search_end:
1.294 _next_arc = e;
1.295 return true;
1.296 }
1.297 @@ -407,7 +414,7 @@
1.298 const IntVector &_source;
1.299 const IntVector &_target;
1.300 const CostVector &_cost;
1.301 - const IntVector &_state;
1.302 + const StateVector &_state;
1.303 const CostVector &_pi;
1.304 int &_in_arc;
1.305 int _search_arc_num;
1.306 @@ -428,7 +435,7 @@
1.307 _next_arc(0)
1.308 {
1.309 // The main parameters of the pivot rule
1.310 - const double LIST_LENGTH_FACTOR = 1.0;
1.311 + const double LIST_LENGTH_FACTOR = 0.25;
1.312 const int MIN_LIST_LENGTH = 10;
1.313 const double MINOR_LIMIT_FACTOR = 0.1;
1.314 const int MIN_MINOR_LIMIT = 3;
1.315 @@ -445,7 +452,7 @@
1.316 /// Find next entering arc
1.317 bool findEnteringArc() {
1.318 Cost min, c;
1.319 - int e, min_arc = _next_arc;
1.320 + int e;
1.321 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.322 // Minor iteration: select the best eligible arc from the
1.323 // current candidate list
1.324 @@ -456,48 +463,44 @@
1.325 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.326 if (c < min) {
1.327 min = c;
1.328 - min_arc = e;
1.329 + _in_arc = e;
1.330 }
1.331 - if (c >= 0) {
1.332 + else if (c >= 0) {
1.333 _candidates[i--] = _candidates[--_curr_length];
1.334 }
1.335 }
1.336 - if (min < 0) {
1.337 - _in_arc = min_arc;
1.338 - return true;
1.339 - }
1.340 + if (min < 0) return true;
1.341 }
1.342
1.343 // Major iteration: build a new candidate list
1.344 min = 0;
1.345 _curr_length = 0;
1.346 - for (e = _next_arc; e < _search_arc_num; ++e) {
1.347 + for (e = _next_arc; e != _search_arc_num; ++e) {
1.348 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.349 if (c < 0) {
1.350 _candidates[_curr_length++] = e;
1.351 if (c < min) {
1.352 min = c;
1.353 - min_arc = e;
1.354 + _in_arc = e;
1.355 }
1.356 - if (_curr_length == _list_length) break;
1.357 + if (_curr_length == _list_length) goto search_end;
1.358 }
1.359 }
1.360 - if (_curr_length < _list_length) {
1.361 - for (e = 0; e < _next_arc; ++e) {
1.362 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.363 - if (c < 0) {
1.364 - _candidates[_curr_length++] = e;
1.365 - if (c < min) {
1.366 - min = c;
1.367 - min_arc = e;
1.368 - }
1.369 - if (_curr_length == _list_length) break;
1.370 + for (e = 0; e != _next_arc; ++e) {
1.371 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.372 + if (c < 0) {
1.373 + _candidates[_curr_length++] = e;
1.374 + if (c < min) {
1.375 + min = c;
1.376 + _in_arc = e;
1.377 }
1.378 + if (_curr_length == _list_length) goto search_end;
1.379 }
1.380 }
1.381 if (_curr_length == 0) return false;
1.382 +
1.383 + search_end:
1.384 _minor_count = 1;
1.385 - _in_arc = min_arc;
1.386 _next_arc = e;
1.387 return true;
1.388 }
1.389 @@ -514,7 +517,7 @@
1.390 const IntVector &_source;
1.391 const IntVector &_target;
1.392 const CostVector &_cost;
1.393 - const IntVector &_state;
1.394 + const StateVector &_state;
1.395 const CostVector &_pi;
1.396 int &_in_arc;
1.397 int _search_arc_num;
1.398 @@ -549,7 +552,7 @@
1.399 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.400 {
1.401 // The main parameters of the pivot rule
1.402 - const double BLOCK_SIZE_FACTOR = 1.5;
1.403 + const double BLOCK_SIZE_FACTOR = 1.0;
1.404 const int MIN_BLOCK_SIZE = 10;
1.405 const double HEAD_LENGTH_FACTOR = 0.1;
1.406 const int MIN_HEAD_LENGTH = 3;
1.407 @@ -567,7 +570,7 @@
1.408 bool findEnteringArc() {
1.409 // Check the current candidate list
1.410 int e;
1.411 - for (int i = 0; i < _curr_length; ++i) {
1.412 + for (int i = 0; i != _curr_length; ++i) {
1.413 e = _candidates[i];
1.414 _cand_cost[e] = _state[e] *
1.415 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.416 @@ -578,39 +581,35 @@
1.417
1.418 // Extend the list
1.419 int cnt = _block_size;
1.420 - int last_arc = 0;
1.421 int limit = _head_length;
1.422
1.423 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.424 + for (e = _next_arc; e != _search_arc_num; ++e) {
1.425 _cand_cost[e] = _state[e] *
1.426 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.427 if (_cand_cost[e] < 0) {
1.428 _candidates[_curr_length++] = e;
1.429 - last_arc = e;
1.430 }
1.431 if (--cnt == 0) {
1.432 - if (_curr_length > limit) break;
1.433 + if (_curr_length > limit) goto search_end;
1.434 limit = 0;
1.435 cnt = _block_size;
1.436 }
1.437 }
1.438 - if (_curr_length <= limit) {
1.439 - for (int e = 0; e < _next_arc; ++e) {
1.440 - _cand_cost[e] = _state[e] *
1.441 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.442 - if (_cand_cost[e] < 0) {
1.443 - _candidates[_curr_length++] = e;
1.444 - last_arc = e;
1.445 - }
1.446 - if (--cnt == 0) {
1.447 - if (_curr_length > limit) break;
1.448 - limit = 0;
1.449 - cnt = _block_size;
1.450 - }
1.451 + for (e = 0; e != _next_arc; ++e) {
1.452 + _cand_cost[e] = _state[e] *
1.453 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.454 + if (_cand_cost[e] < 0) {
1.455 + _candidates[_curr_length++] = e;
1.456 + }
1.457 + if (--cnt == 0) {
1.458 + if (_curr_length > limit) goto search_end;
1.459 + limit = 0;
1.460 + cnt = _block_size;
1.461 }
1.462 }
1.463 if (_curr_length == 0) return false;
1.464 - _next_arc = last_arc + 1;
1.465 +
1.466 + search_end:
1.467
1.468 // Make heap of the candidate list (approximating a partial sort)
1.469 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.470 @@ -618,6 +617,7 @@
1.471
1.472 // Pop the first element of the heap
1.473 _in_arc = _candidates[0];
1.474 + _next_arc = e;
1.475 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.476 _sort_func );
1.477 _curr_length = std::min(_head_length, _curr_length - 1);
1.478 @@ -633,69 +633,25 @@
1.479 /// The constructor of the class.
1.480 ///
1.481 /// \param graph The digraph the algorithm runs on.
1.482 - NetworkSimplex(const GR& graph) :
1.483 + /// \param arc_mixing Indicate if the arcs have to be stored in a
1.484 + /// mixed order in the internal data structure.
1.485 + /// In special cases, it could lead to better overall performance,
1.486 + /// but it is usually slower. Therefore it is disabled by default.
1.487 + NetworkSimplex(const GR& graph, bool arc_mixing = false) :
1.488 _graph(graph), _node_id(graph), _arc_id(graph),
1.489 + _arc_mixing(arc_mixing),
1.490 + MAX(std::numeric_limits<Value>::max()),
1.491 INF(std::numeric_limits<Value>::has_infinity ?
1.492 - std::numeric_limits<Value>::infinity() :
1.493 - std::numeric_limits<Value>::max())
1.494 + std::numeric_limits<Value>::infinity() : MAX)
1.495 {
1.496 - // Check the value types
1.497 + // Check the number types
1.498 LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.499 "The flow type of NetworkSimplex must be signed");
1.500 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.501 "The cost type of NetworkSimplex must be signed");
1.502 -
1.503 - // Resize vectors
1.504 - _node_num = countNodes(_graph);
1.505 - _arc_num = countArcs(_graph);
1.506 - int all_node_num = _node_num + 1;
1.507 - int max_arc_num = _arc_num + 2 * _node_num;
1.508
1.509 - _source.resize(max_arc_num);
1.510 - _target.resize(max_arc_num);
1.511 -
1.512 - _lower.resize(_arc_num);
1.513 - _upper.resize(_arc_num);
1.514 - _cap.resize(max_arc_num);
1.515 - _cost.resize(max_arc_num);
1.516 - _supply.resize(all_node_num);
1.517 - _flow.resize(max_arc_num);
1.518 - _pi.resize(all_node_num);
1.519 -
1.520 - _parent.resize(all_node_num);
1.521 - _pred.resize(all_node_num);
1.522 - _forward.resize(all_node_num);
1.523 - _thread.resize(all_node_num);
1.524 - _rev_thread.resize(all_node_num);
1.525 - _succ_num.resize(all_node_num);
1.526 - _last_succ.resize(all_node_num);
1.527 - _state.resize(max_arc_num);
1.528 -
1.529 - // Copy the graph (store the arcs in a mixed order)
1.530 - int i = 0;
1.531 - for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.532 - _node_id[n] = i;
1.533 - }
1.534 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.535 - i = 0;
1.536 - for (ArcIt a(_graph); a != INVALID; ++a) {
1.537 - _arc_id[a] = i;
1.538 - _source[i] = _node_id[_graph.source(a)];
1.539 - _target[i] = _node_id[_graph.target(a)];
1.540 - if ((i += k) >= _arc_num) i = (i % k) + 1;
1.541 - }
1.542 -
1.543 - // Initialize maps
1.544 - for (int i = 0; i != _node_num; ++i) {
1.545 - _supply[i] = 0;
1.546 - }
1.547 - for (int i = 0; i != _arc_num; ++i) {
1.548 - _lower[i] = 0;
1.549 - _upper[i] = INF;
1.550 - _cost[i] = 1;
1.551 - }
1.552 - _have_lower = false;
1.553 - _stype = GEQ;
1.554 + // Reset data structures
1.555 + reset();
1.556 }
1.557
1.558 /// \name Parameters
1.559 @@ -729,7 +685,7 @@
1.560 /// This function sets the upper bounds (capacities) on the arcs.
1.561 /// If it is not used before calling \ref run(), the upper bounds
1.562 /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.563 - /// unbounded from above on each arc).
1.564 + /// unbounded from above).
1.565 ///
1.566 /// \param map An arc map storing the upper bounds.
1.567 /// Its \c Value type must be convertible to the \c Value type
1.568 @@ -768,7 +724,6 @@
1.569 /// This function sets the supply values of the nodes.
1.570 /// If neither this function nor \ref stSupply() is used before
1.571 /// calling \ref run(), the supply of each node will be set to zero.
1.572 - /// (It makes sense only if non-zero lower bounds are given.)
1.573 ///
1.574 /// \param map A node map storing the supply values.
1.575 /// Its \c Value type must be convertible to the \c Value type
1.576 @@ -789,7 +744,6 @@
1.577 /// and the required flow value.
1.578 /// If neither this function nor \ref supplyMap() is used before
1.579 /// calling \ref run(), the supply of each node will be set to zero.
1.580 - /// (It makes sense only if non-zero lower bounds are given.)
1.581 ///
1.582 /// Using this function has the same effect as using \ref supplyMap()
1.583 /// with such a map in which \c k is assigned to \c s, \c -k is
1.584 @@ -809,14 +763,14 @@
1.585 _supply[_node_id[t]] = -k;
1.586 return *this;
1.587 }
1.588 -
1.589 +
1.590 /// \brief Set the type of the supply constraints.
1.591 ///
1.592 /// This function sets the type of the supply/demand constraints.
1.593 /// If it is not used before calling \ref run(), the \ref GEQ supply
1.594 /// type will be used.
1.595 ///
1.596 - /// For more information see \ref SupplyType.
1.597 + /// For more information, see \ref SupplyType.
1.598 ///
1.599 /// \return <tt>(*this)</tt>
1.600 NetworkSimplex& supplyType(SupplyType supply_type) {
1.601 @@ -835,7 +789,7 @@
1.602 ///
1.603 /// This function runs the algorithm.
1.604 /// The paramters can be specified using functions \ref lowerMap(),
1.605 - /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.606 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.607 /// \ref supplyType().
1.608 /// For example,
1.609 /// \code
1.610 @@ -844,15 +798,15 @@
1.611 /// .supplyMap(sup).run();
1.612 /// \endcode
1.613 ///
1.614 - /// This function can be called more than once. All the parameters
1.615 - /// that have been given are kept for the next call, unless
1.616 - /// \ref reset() is called, thus only the modified parameters
1.617 - /// have to be set again. See \ref reset() for examples.
1.618 - /// However the underlying digraph must not be modified after this
1.619 - /// class have been constructed, since it copies and extends the graph.
1.620 + /// This function can be called more than once. All the given parameters
1.621 + /// are kept for the next call, unless \ref resetParams() or \ref reset()
1.622 + /// is used, thus only the modified parameters have to be set again.
1.623 + /// If the underlying digraph was also modified after the construction
1.624 + /// of the class (or the last \ref reset() call), then the \ref reset()
1.625 + /// function must be called.
1.626 ///
1.627 /// \param pivot_rule The pivot rule that will be used during the
1.628 - /// algorithm. For more information see \ref PivotRule.
1.629 + /// algorithm. For more information, see \ref PivotRule.
1.630 ///
1.631 /// \return \c INFEASIBLE if no feasible flow exists,
1.632 /// \n \c OPTIMAL if the problem has optimal solution
1.633 @@ -863,6 +817,7 @@
1.634 /// cost and infinite upper bound.
1.635 ///
1.636 /// \see ProblemType, PivotRule
1.637 + /// \see resetParams(), reset()
1.638 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.639 if (!init()) return INFEASIBLE;
1.640 return start(pivot_rule);
1.641 @@ -874,11 +829,12 @@
1.642 /// before using functions \ref lowerMap(), \ref upperMap(),
1.643 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
1.644 ///
1.645 - /// It is useful for multiple run() calls. If this function is not
1.646 - /// used, all the parameters given before are kept for the next
1.647 - /// \ref run() call.
1.648 - /// However the underlying digraph must not be modified after this
1.649 - /// class have been constructed, since it copies and extends the graph.
1.650 + /// It is useful for multiple \ref run() calls. Basically, all the given
1.651 + /// parameters are kept for the next \ref run() call, unless
1.652 + /// \ref resetParams() or \ref reset() is used.
1.653 + /// If the underlying digraph was also modified after the construction
1.654 + /// of the class or the last \ref reset() call, then the \ref reset()
1.655 + /// function must be used, otherwise \ref resetParams() is sufficient.
1.656 ///
1.657 /// For example,
1.658 /// \code
1.659 @@ -888,20 +844,22 @@
1.660 /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.661 /// .supplyMap(sup).run();
1.662 ///
1.663 - /// // Run again with modified cost map (reset() is not called,
1.664 + /// // Run again with modified cost map (resetParams() is not called,
1.665 /// // so only the cost map have to be set again)
1.666 /// cost[e] += 100;
1.667 /// ns.costMap(cost).run();
1.668 ///
1.669 - /// // Run again from scratch using reset()
1.670 + /// // Run again from scratch using resetParams()
1.671 /// // (the lower bounds will be set to zero on all arcs)
1.672 - /// ns.reset();
1.673 + /// ns.resetParams();
1.674 /// ns.upperMap(capacity).costMap(cost)
1.675 /// .supplyMap(sup).run();
1.676 /// \endcode
1.677 ///
1.678 /// \return <tt>(*this)</tt>
1.679 - NetworkSimplex& reset() {
1.680 + ///
1.681 + /// \see reset(), run()
1.682 + NetworkSimplex& resetParams() {
1.683 for (int i = 0; i != _node_num; ++i) {
1.684 _supply[i] = 0;
1.685 }
1.686 @@ -915,6 +873,83 @@
1.687 return *this;
1.688 }
1.689
1.690 + /// \brief Reset the internal data structures and all the parameters
1.691 + /// that have been given before.
1.692 + ///
1.693 + /// This function resets the internal data structures and all the
1.694 + /// paramaters that have been given before using functions \ref lowerMap(),
1.695 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.696 + /// \ref supplyType().
1.697 + ///
1.698 + /// It is useful for multiple \ref run() calls. Basically, all the given
1.699 + /// parameters are kept for the next \ref run() call, unless
1.700 + /// \ref resetParams() or \ref reset() is used.
1.701 + /// If the underlying digraph was also modified after the construction
1.702 + /// of the class or the last \ref reset() call, then the \ref reset()
1.703 + /// function must be used, otherwise \ref resetParams() is sufficient.
1.704 + ///
1.705 + /// See \ref resetParams() for examples.
1.706 + ///
1.707 + /// \return <tt>(*this)</tt>
1.708 + ///
1.709 + /// \see resetParams(), run()
1.710 + NetworkSimplex& reset() {
1.711 + // Resize vectors
1.712 + _node_num = countNodes(_graph);
1.713 + _arc_num = countArcs(_graph);
1.714 + int all_node_num = _node_num + 1;
1.715 + int max_arc_num = _arc_num + 2 * _node_num;
1.716 +
1.717 + _source.resize(max_arc_num);
1.718 + _target.resize(max_arc_num);
1.719 +
1.720 + _lower.resize(_arc_num);
1.721 + _upper.resize(_arc_num);
1.722 + _cap.resize(max_arc_num);
1.723 + _cost.resize(max_arc_num);
1.724 + _supply.resize(all_node_num);
1.725 + _flow.resize(max_arc_num);
1.726 + _pi.resize(all_node_num);
1.727 +
1.728 + _parent.resize(all_node_num);
1.729 + _pred.resize(all_node_num);
1.730 + _forward.resize(all_node_num);
1.731 + _thread.resize(all_node_num);
1.732 + _rev_thread.resize(all_node_num);
1.733 + _succ_num.resize(all_node_num);
1.734 + _last_succ.resize(all_node_num);
1.735 + _state.resize(max_arc_num);
1.736 +
1.737 + // Copy the graph
1.738 + int i = 0;
1.739 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.740 + _node_id[n] = i;
1.741 + }
1.742 + if (_arc_mixing) {
1.743 + // Store the arcs in a mixed order
1.744 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.745 + int i = 0, j = 0;
1.746 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.747 + _arc_id[a] = i;
1.748 + _source[i] = _node_id[_graph.source(a)];
1.749 + _target[i] = _node_id[_graph.target(a)];
1.750 + if ((i += k) >= _arc_num) i = ++j;
1.751 + }
1.752 + } else {
1.753 + // Store the arcs in the original order
1.754 + int i = 0;
1.755 + for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
1.756 + _arc_id[a] = i;
1.757 + _source[i] = _node_id[_graph.source(a)];
1.758 + _target[i] = _node_id[_graph.target(a)];
1.759 + }
1.760 + }
1.761 +
1.762 + // Reset parameters
1.763 + resetParams();
1.764 + return *this;
1.765 + }
1.766 +
1.767 /// @}
1.768
1.769 /// \name Query Functions
1.770 @@ -1024,9 +1059,9 @@
1.771 for (int i = 0; i != _arc_num; ++i) {
1.772 Value c = _lower[i];
1.773 if (c >= 0) {
1.774 - _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
1.775 + _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
1.776 } else {
1.777 - _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
1.778 + _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
1.779 }
1.780 _supply[_source[i]] -= c;
1.781 _supply[_target[i]] += c;
1.782 @@ -1054,7 +1089,7 @@
1.783 _flow[i] = 0;
1.784 _state[i] = STATE_LOWER;
1.785 }
1.786 -
1.787 +
1.788 // Set data for the artificial root node
1.789 _root = _node_num;
1.790 _parent[_root] = -1;
1.791 @@ -1218,7 +1253,7 @@
1.792 for (int u = first; u != join; u = _parent[u]) {
1.793 e = _pred[u];
1.794 d = _forward[u] ?
1.795 - _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1.796 + _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
1.797 if (d < delta) {
1.798 delta = d;
1.799 u_out = u;
1.800 @@ -1228,8 +1263,8 @@
1.801 // Search the cycle along the path form the second node to the root
1.802 for (int u = second; u != join; u = _parent[u]) {
1.803 e = _pred[u];
1.804 - d = _forward[u] ?
1.805 - (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1.806 + d = _forward[u] ?
1.807 + (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
1.808 if (d <= delta) {
1.809 delta = d;
1.810 u_out = u;
1.811 @@ -1330,7 +1365,7 @@
1.812 }
1.813
1.814 // Update _rev_thread using the new _thread values
1.815 - for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1.816 + for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1.817 u = _dirty_revs[i];
1.818 _rev_thread[_thread[u]] = u;
1.819 }
1.820 @@ -1402,6 +1437,100 @@
1.821 }
1.822 }
1.823
1.824 + // Heuristic initial pivots
1.825 + bool initialPivots() {
1.826 + Value curr, total = 0;
1.827 + std::vector<Node> supply_nodes, demand_nodes;
1.828 + for (NodeIt u(_graph); u != INVALID; ++u) {
1.829 + curr = _supply[_node_id[u]];
1.830 + if (curr > 0) {
1.831 + total += curr;
1.832 + supply_nodes.push_back(u);
1.833 + }
1.834 + else if (curr < 0) {
1.835 + demand_nodes.push_back(u);
1.836 + }
1.837 + }
1.838 + if (_sum_supply > 0) total -= _sum_supply;
1.839 + if (total <= 0) return true;
1.840 +
1.841 + IntVector arc_vector;
1.842 + if (_sum_supply >= 0) {
1.843 + if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
1.844 + // Perform a reverse graph search from the sink to the source
1.845 + typename GR::template NodeMap<bool> reached(_graph, false);
1.846 + Node s = supply_nodes[0], t = demand_nodes[0];
1.847 + std::vector<Node> stack;
1.848 + reached[t] = true;
1.849 + stack.push_back(t);
1.850 + while (!stack.empty()) {
1.851 + Node u, v = stack.back();
1.852 + stack.pop_back();
1.853 + if (v == s) break;
1.854 + for (InArcIt a(_graph, v); a != INVALID; ++a) {
1.855 + if (reached[u = _graph.source(a)]) continue;
1.856 + int j = _arc_id[a];
1.857 + if (_cap[j] >= total) {
1.858 + arc_vector.push_back(j);
1.859 + reached[u] = true;
1.860 + stack.push_back(u);
1.861 + }
1.862 + }
1.863 + }
1.864 + } else {
1.865 + // Find the min. cost incomming arc for each demand node
1.866 + for (int i = 0; i != int(demand_nodes.size()); ++i) {
1.867 + Node v = demand_nodes[i];
1.868 + Cost c, min_cost = std::numeric_limits<Cost>::max();
1.869 + Arc min_arc = INVALID;
1.870 + for (InArcIt a(_graph, v); a != INVALID; ++a) {
1.871 + c = _cost[_arc_id[a]];
1.872 + if (c < min_cost) {
1.873 + min_cost = c;
1.874 + min_arc = a;
1.875 + }
1.876 + }
1.877 + if (min_arc != INVALID) {
1.878 + arc_vector.push_back(_arc_id[min_arc]);
1.879 + }
1.880 + }
1.881 + }
1.882 + } else {
1.883 + // Find the min. cost outgoing arc for each supply node
1.884 + for (int i = 0; i != int(supply_nodes.size()); ++i) {
1.885 + Node u = supply_nodes[i];
1.886 + Cost c, min_cost = std::numeric_limits<Cost>::max();
1.887 + Arc min_arc = INVALID;
1.888 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.889 + c = _cost[_arc_id[a]];
1.890 + if (c < min_cost) {
1.891 + min_cost = c;
1.892 + min_arc = a;
1.893 + }
1.894 + }
1.895 + if (min_arc != INVALID) {
1.896 + arc_vector.push_back(_arc_id[min_arc]);
1.897 + }
1.898 + }
1.899 + }
1.900 +
1.901 + // Perform heuristic initial pivots
1.902 + for (int i = 0; i != int(arc_vector.size()); ++i) {
1.903 + in_arc = arc_vector[i];
1.904 + if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
1.905 + _pi[_target[in_arc]]) >= 0) continue;
1.906 + findJoinNode();
1.907 + bool change = findLeavingArc();
1.908 + if (delta >= MAX) return false;
1.909 + changeFlow(change);
1.910 + if (change) {
1.911 + updateTreeStructure();
1.912 + updatePotential();
1.913 + }
1.914 + }
1.915 + return true;
1.916 + }
1.917 +
1.918 // Execute the algorithm
1.919 ProblemType start(PivotRule pivot_rule) {
1.920 // Select the pivot rule implementation
1.921 @@ -1424,18 +1553,21 @@
1.922 ProblemType start() {
1.923 PivotRuleImpl pivot(*this);
1.924
1.925 + // Perform heuristic initial pivots
1.926 + if (!initialPivots()) return UNBOUNDED;
1.927 +
1.928 // Execute the Network Simplex algorithm
1.929 while (pivot.findEnteringArc()) {
1.930 findJoinNode();
1.931 bool change = findLeavingArc();
1.932 - if (delta >= INF) return UNBOUNDED;
1.933 + if (delta >= MAX) return UNBOUNDED;
1.934 changeFlow(change);
1.935 if (change) {
1.936 updateTreeStructure();
1.937 updatePotential();
1.938 }
1.939 }
1.940 -
1.941 +
1.942 // Check feasibility
1.943 for (int e = _search_arc_num; e != _all_arc_num; ++e) {
1.944 if (_flow[e] != 0) return INFEASIBLE;
1.945 @@ -1452,7 +1584,7 @@
1.946 }
1.947 }
1.948 }
1.949 -
1.950 +
1.951 // Shift potentials to meet the requirements of the GEQ/LEQ type
1.952 // optimality conditions
1.953 if (_sum_supply == 0) {