1.1 --- a/lemon/network_simplex.h Mon Jul 16 16:21:40 2018 +0200
1.2 +++ b/lemon/network_simplex.h Wed Oct 17 19:14:07 2018 +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-2010
1.8 + * Copyright (C) 2003-2013
1.9 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 *
1.12 @@ -41,16 +41,17 @@
1.13 ///
1.14 /// \ref NetworkSimplex implements the primal Network Simplex algorithm
1.15 /// for finding a \ref min_cost_flow "minimum cost flow"
1.16 - /// \ref amo93networkflows, \ref dantzig63linearprog,
1.17 - /// \ref kellyoneill91netsimplex.
1.18 + /// \cite amo93networkflows, \cite dantzig63linearprog,
1.19 + /// \cite kellyoneill91netsimplex.
1.20 /// This algorithm is a highly efficient specialized version of the
1.21 /// linear programming simplex method directly for the minimum cost
1.22 /// flow problem.
1.23 ///
1.24 - /// In general, %NetworkSimplex is the fastest implementation available
1.25 - /// in LEMON for this problem.
1.26 - /// Moreover, it supports both directions of the supply/demand inequality
1.27 - /// constraints. For more information, see \ref SupplyType.
1.28 + /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
1.29 + /// implementations available in LEMON for solving this problem.
1.30 + /// (For more information, see \ref min_cost_flow_algs "the module page".)
1.31 + /// Furthermore, this class supports both directions of the supply/demand
1.32 + /// inequality constraints. For more information, see \ref SupplyType.
1.33 ///
1.34 /// Most of the parameters of the problem (except for the digraph)
1.35 /// can be given using separate functions, and the algorithm can be
1.36 @@ -63,7 +64,8 @@
1.37 /// \tparam C The number type used for costs and potentials in the
1.38 /// algorithm. By default, it is the same as \c V.
1.39 ///
1.40 - /// \warning Both number types must be signed and all input data must
1.41 + /// \warning Both \c V and \c C must be signed number types.
1.42 + /// \warning All input data (capacities, supply values, and costs) must
1.43 /// be integer.
1.44 ///
1.45 /// \note %NetworkSimplex provides five different pivot rule
1.46 @@ -121,14 +123,17 @@
1.47 /// Enum type containing constants for selecting the pivot rule for
1.48 /// the \ref run() function.
1.49 ///
1.50 - /// \ref NetworkSimplex provides five different pivot rule
1.51 - /// implementations that significantly affect the running time
1.52 + /// \ref NetworkSimplex provides five different implementations for
1.53 + /// the pivot strategy that significantly affects the running time
1.54 /// of the algorithm.
1.55 - /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
1.56 - /// proved to be the most efficient and the most robust on various
1.57 - /// test inputs.
1.58 - /// However, another pivot rule can be selected using the \ref run()
1.59 - /// function with the proper parameter.
1.60 + /// According to experimental tests conducted on various problem
1.61 + /// instances, \ref BLOCK_SEARCH "Block Search" and
1.62 + /// \ref ALTERING_LIST "Altering Candidate List" rules turned out
1.63 + /// to be the most efficient.
1.64 + /// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that
1.65 + /// seemed to be slightly more robust, it is used by default.
1.66 + /// However, another pivot rule can easily be selected using the
1.67 + /// \ref run() function with the proper parameter.
1.68 enum PivotRule {
1.69
1.70 /// The \e First \e Eligible pivot rule.
1.71 @@ -154,7 +159,7 @@
1.72
1.73 /// The \e Altering \e Candidate \e List pivot rule.
1.74 /// It is a modified version of the Candidate List method.
1.75 - /// It keeps only the several best eligible arcs from the former
1.76 + /// It keeps only a few of the best eligible arcs from the former
1.77 /// candidate list and extends this list in every iteration.
1.78 ALTERING_LIST
1.79 };
1.80 @@ -166,8 +171,9 @@
1.81 typedef std::vector<int> IntVector;
1.82 typedef std::vector<Value> ValueVector;
1.83 typedef std::vector<Cost> CostVector;
1.84 - typedef std::vector<char> BoolVector;
1.85 - // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.86 + typedef std::vector<signed char> CharVector;
1.87 + // Note: vector<signed char> is used instead of vector<ArcState> and
1.88 + // vector<ArcDirection> for efficiency reasons
1.89
1.90 // State constants for arcs
1.91 enum ArcState {
1.92 @@ -176,9 +182,11 @@
1.93 STATE_LOWER = 1
1.94 };
1.95
1.96 - typedef std::vector<signed char> StateVector;
1.97 - // Note: vector<signed char> is used instead of vector<ArcState> for
1.98 - // efficiency reasons
1.99 + // Direction constants for tree arcs
1.100 + enum ArcDirection {
1.101 + DIR_DOWN = -1,
1.102 + DIR_UP = 1
1.103 + };
1.104
1.105 private:
1.106
1.107 @@ -190,7 +198,7 @@
1.108 int _search_arc_num;
1.109
1.110 // Parameters of the problem
1.111 - bool _have_lower;
1.112 + bool _has_lower;
1.113 SupplyType _stype;
1.114 Value _sum_supply;
1.115
1.116 @@ -217,15 +225,13 @@
1.117 IntVector _rev_thread;
1.118 IntVector _succ_num;
1.119 IntVector _last_succ;
1.120 + CharVector _pred_dir;
1.121 + CharVector _state;
1.122 IntVector _dirty_revs;
1.123 - BoolVector _forward;
1.124 - StateVector _state;
1.125 int _root;
1.126
1.127 // Temporary data used in the current pivot iteration
1.128 int in_arc, join, u_in, v_in, u_out, v_out;
1.129 - int first, second, right, last;
1.130 - int stem, par_stem, new_stem;
1.131 Value delta;
1.132
1.133 const Value MAX;
1.134 @@ -250,7 +256,7 @@
1.135 const IntVector &_source;
1.136 const IntVector &_target;
1.137 const CostVector &_cost;
1.138 - const StateVector &_state;
1.139 + const CharVector &_state;
1.140 const CostVector &_pi;
1.141 int &_in_arc;
1.142 int _search_arc_num;
1.143 @@ -302,7 +308,7 @@
1.144 const IntVector &_source;
1.145 const IntVector &_target;
1.146 const CostVector &_cost;
1.147 - const StateVector &_state;
1.148 + const CharVector &_state;
1.149 const CostVector &_pi;
1.150 int &_in_arc;
1.151 int _search_arc_num;
1.152 @@ -341,7 +347,7 @@
1.153 const IntVector &_source;
1.154 const IntVector &_target;
1.155 const CostVector &_cost;
1.156 - const StateVector &_state;
1.157 + const CharVector &_state;
1.158 const CostVector &_pi;
1.159 int &_in_arc;
1.160 int _search_arc_num;
1.161 @@ -414,7 +420,7 @@
1.162 const IntVector &_source;
1.163 const IntVector &_target;
1.164 const CostVector &_cost;
1.165 - const StateVector &_state;
1.166 + const CharVector &_state;
1.167 const CostVector &_pi;
1.168 int &_in_arc;
1.169 int _search_arc_num;
1.170 @@ -517,7 +523,7 @@
1.171 const IntVector &_source;
1.172 const IntVector &_target;
1.173 const CostVector &_cost;
1.174 - const StateVector &_state;
1.175 + const CharVector &_state;
1.176 const CostVector &_pi;
1.177 int &_in_arc;
1.178 int _search_arc_num;
1.179 @@ -536,7 +542,7 @@
1.180 public:
1.181 SortFunc(const CostVector &map) : _map(map) {}
1.182 bool operator()(int left, int right) {
1.183 - return _map[left] > _map[right];
1.184 + return _map[left] < _map[right];
1.185 }
1.186 };
1.187
1.188 @@ -554,7 +560,7 @@
1.189 // The main parameters of the pivot rule
1.190 const double BLOCK_SIZE_FACTOR = 1.0;
1.191 const int MIN_BLOCK_SIZE = 10;
1.192 - const double HEAD_LENGTH_FACTOR = 0.1;
1.193 + const double HEAD_LENGTH_FACTOR = 0.01;
1.194 const int MIN_HEAD_LENGTH = 3;
1.195
1.196 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.197 @@ -570,11 +576,13 @@
1.198 bool findEnteringArc() {
1.199 // Check the current candidate list
1.200 int e;
1.201 + Cost c;
1.202 for (int i = 0; i != _curr_length; ++i) {
1.203 e = _candidates[i];
1.204 - _cand_cost[e] = _state[e] *
1.205 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.206 - if (_cand_cost[e] >= 0) {
1.207 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.208 + if (c < 0) {
1.209 + _cand_cost[e] = c;
1.210 + } else {
1.211 _candidates[i--] = _candidates[--_curr_length];
1.212 }
1.213 }
1.214 @@ -584,9 +592,9 @@
1.215 int limit = _head_length;
1.216
1.217 for (e = _next_arc; e != _search_arc_num; ++e) {
1.218 - _cand_cost[e] = _state[e] *
1.219 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.220 - if (_cand_cost[e] < 0) {
1.221 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.222 + if (c < 0) {
1.223 + _cand_cost[e] = c;
1.224 _candidates[_curr_length++] = e;
1.225 }
1.226 if (--cnt == 0) {
1.227 @@ -596,9 +604,9 @@
1.228 }
1.229 }
1.230 for (e = 0; e != _next_arc; ++e) {
1.231 - _cand_cost[e] = _state[e] *
1.232 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.233 - if (_cand_cost[e] < 0) {
1.234 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.235 + if (c < 0) {
1.236 + _cand_cost[e] = c;
1.237 _candidates[_curr_length++] = e;
1.238 }
1.239 if (--cnt == 0) {
1.240 @@ -611,16 +619,16 @@
1.241
1.242 search_end:
1.243
1.244 - // Make heap of the candidate list (approximating a partial sort)
1.245 - make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.246 - _sort_func );
1.247 + // Perform partial sort operation on the candidate list
1.248 + int new_length = std::min(_head_length + 1, _curr_length);
1.249 + std::partial_sort(_candidates.begin(), _candidates.begin() + new_length,
1.250 + _candidates.begin() + _curr_length, _sort_func);
1.251
1.252 - // Pop the first element of the heap
1.253 + // Select the entering arc and remove it from the list
1.254 _in_arc = _candidates[0];
1.255 _next_arc = e;
1.256 - pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.257 - _sort_func );
1.258 - _curr_length = std::min(_head_length, _curr_length - 1);
1.259 + _candidates[0] = _candidates[new_length - 1];
1.260 + _curr_length = new_length - 1;
1.261 return true;
1.262 }
1.263
1.264 @@ -633,11 +641,12 @@
1.265 /// The constructor of the class.
1.266 ///
1.267 /// \param graph The digraph the algorithm runs on.
1.268 - /// \param arc_mixing Indicate if the arcs have to be stored in a
1.269 + /// \param arc_mixing Indicate if the arcs will be stored in a
1.270 /// mixed order in the internal data structure.
1.271 - /// In special cases, it could lead to better overall performance,
1.272 - /// but it is usually slower. Therefore it is disabled by default.
1.273 - NetworkSimplex(const GR& graph, bool arc_mixing = false) :
1.274 + /// In general, it leads to similar performance as using the original
1.275 + /// arc order, but it makes the algorithm more robust and in special
1.276 + /// cases, even significantly faster. Therefore, it is enabled by default.
1.277 + NetworkSimplex(const GR& graph, bool arc_mixing = true) :
1.278 _graph(graph), _node_id(graph), _arc_id(graph),
1.279 _arc_mixing(arc_mixing),
1.280 MAX(std::numeric_limits<Value>::max()),
1.281 @@ -673,7 +682,7 @@
1.282 /// \return <tt>(*this)</tt>
1.283 template <typename LowerMap>
1.284 NetworkSimplex& lowerMap(const LowerMap& map) {
1.285 - _have_lower = true;
1.286 + _has_lower = true;
1.287 for (ArcIt a(_graph); a != INVALID; ++a) {
1.288 _lower[_arc_id[a]] = map[a];
1.289 }
1.290 @@ -730,6 +739,8 @@
1.291 /// of the algorithm.
1.292 ///
1.293 /// \return <tt>(*this)</tt>
1.294 + ///
1.295 + /// \sa supplyType()
1.296 template<typename SupplyMap>
1.297 NetworkSimplex& supplyMap(const SupplyMap& map) {
1.298 for (NodeIt n(_graph); n != INVALID; ++n) {
1.299 @@ -746,7 +757,7 @@
1.300 /// calling \ref run(), the supply of each node will be set to zero.
1.301 ///
1.302 /// Using this function has the same effect as using \ref supplyMap()
1.303 - /// with such a map in which \c k is assigned to \c s, \c -k is
1.304 + /// with a map in which \c k is assigned to \c s, \c -k is
1.305 /// assigned to \c t and all other nodes have zero supply value.
1.306 ///
1.307 /// \param s The source node.
1.308 @@ -868,7 +879,7 @@
1.309 _upper[i] = INF;
1.310 _cost[i] = 1;
1.311 }
1.312 - _have_lower = false;
1.313 + _has_lower = false;
1.314 _stype = GEQ;
1.315 return *this;
1.316 }
1.317 @@ -913,7 +924,7 @@
1.318
1.319 _parent.resize(all_node_num);
1.320 _pred.resize(all_node_num);
1.321 - _forward.resize(all_node_num);
1.322 + _pred_dir.resize(all_node_num);
1.323 _thread.resize(all_node_num);
1.324 _rev_thread.resize(all_node_num);
1.325 _succ_num.resize(all_node_num);
1.326 @@ -925,15 +936,15 @@
1.327 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.328 _node_id[n] = i;
1.329 }
1.330 - if (_arc_mixing) {
1.331 + if (_arc_mixing && _node_num > 1) {
1.332 // Store the arcs in a mixed order
1.333 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.334 + const int skip = std::max(_arc_num / _node_num, 3);
1.335 int i = 0, j = 0;
1.336 for (ArcIt a(_graph); a != INVALID; ++a) {
1.337 _arc_id[a] = i;
1.338 _source[i] = _node_id[_graph.source(a)];
1.339 _target[i] = _node_id[_graph.target(a)];
1.340 - if ((i += k) >= _arc_num) i = ++j;
1.341 + if ((i += skip) >= _arc_num) i = ++j;
1.342 }
1.343 } else {
1.344 // Store the arcs in the original order
1.345 @@ -962,7 +973,7 @@
1.346 /// \brief Return the total cost of the found flow.
1.347 ///
1.348 /// This function returns the total cost of the found flow.
1.349 - /// Its complexity is O(e).
1.350 + /// Its complexity is O(m).
1.351 ///
1.352 /// \note The return type of the function can be specified as a
1.353 /// template parameter. For example,
1.354 @@ -999,7 +1010,8 @@
1.355 return _flow[_arc_id[a]];
1.356 }
1.357
1.358 - /// \brief Return the flow map (the primal solution).
1.359 + /// \brief Copy the flow values (the primal solution) into the
1.360 + /// given map.
1.361 ///
1.362 /// This function copies the flow value on each arc into the given
1.363 /// map. The \c Value type of the algorithm must be convertible to
1.364 @@ -1023,7 +1035,8 @@
1.365 return _pi[_node_id[n]];
1.366 }
1.367
1.368 - /// \brief Return the potential map (the dual solution).
1.369 + /// \brief Copy the potential values (the dual solution) into the
1.370 + /// given map.
1.371 ///
1.372 /// This function copies the potential (dual value) of each node
1.373 /// into the given map.
1.374 @@ -1054,8 +1067,12 @@
1.375 if ( !((_stype == GEQ && _sum_supply <= 0) ||
1.376 (_stype == LEQ && _sum_supply >= 0)) ) return false;
1.377
1.378 + // Check lower and upper bounds
1.379 + LEMON_DEBUG(checkBoundMaps(),
1.380 + "Upper bounds must be greater or equal to the lower bounds");
1.381 +
1.382 // Remove non-zero lower bounds
1.383 - if (_have_lower) {
1.384 + if (_has_lower) {
1.385 for (int i = 0; i != _arc_num; ++i) {
1.386 Value c = _lower[i];
1.387 if (c >= 0) {
1.388 @@ -1116,14 +1133,14 @@
1.389 _cap[e] = INF;
1.390 _state[e] = STATE_TREE;
1.391 if (_supply[u] >= 0) {
1.392 - _forward[u] = true;
1.393 + _pred_dir[u] = DIR_UP;
1.394 _pi[u] = 0;
1.395 _source[e] = u;
1.396 _target[e] = _root;
1.397 _flow[e] = _supply[u];
1.398 _cost[e] = 0;
1.399 } else {
1.400 - _forward[u] = false;
1.401 + _pred_dir[u] = DIR_DOWN;
1.402 _pi[u] = ART_COST;
1.403 _source[e] = _root;
1.404 _target[e] = u;
1.405 @@ -1143,7 +1160,7 @@
1.406 _succ_num[u] = 1;
1.407 _last_succ[u] = u;
1.408 if (_supply[u] >= 0) {
1.409 - _forward[u] = true;
1.410 + _pred_dir[u] = DIR_UP;
1.411 _pi[u] = 0;
1.412 _pred[u] = e;
1.413 _source[e] = u;
1.414 @@ -1153,7 +1170,7 @@
1.415 _cost[e] = 0;
1.416 _state[e] = STATE_TREE;
1.417 } else {
1.418 - _forward[u] = false;
1.419 + _pred_dir[u] = DIR_DOWN;
1.420 _pi[u] = ART_COST;
1.421 _pred[u] = f;
1.422 _source[f] = _root;
1.423 @@ -1184,7 +1201,7 @@
1.424 _succ_num[u] = 1;
1.425 _last_succ[u] = u;
1.426 if (_supply[u] <= 0) {
1.427 - _forward[u] = false;
1.428 + _pred_dir[u] = DIR_DOWN;
1.429 _pi[u] = 0;
1.430 _pred[u] = e;
1.431 _source[e] = _root;
1.432 @@ -1194,7 +1211,7 @@
1.433 _cost[e] = 0;
1.434 _state[e] = STATE_TREE;
1.435 } else {
1.436 - _forward[u] = true;
1.437 + _pred_dir[u] = DIR_UP;
1.438 _pi[u] = -ART_COST;
1.439 _pred[u] = f;
1.440 _source[f] = u;
1.441 @@ -1218,6 +1235,15 @@
1.442 return true;
1.443 }
1.444
1.445 + // Check if the upper bound is greater than or equal to the lower bound
1.446 + // on each arc.
1.447 + bool checkBoundMaps() {
1.448 + for (int j = 0; j != _arc_num; ++j) {
1.449 + if (_upper[j] < _lower[j]) return false;
1.450 + }
1.451 + return true;
1.452 + }
1.453 +
1.454 // Find the join node
1.455 void findJoinNode() {
1.456 int u = _source[in_arc];
1.457 @@ -1237,6 +1263,7 @@
1.458 bool findLeavingArc() {
1.459 // Initialize first and second nodes according to the direction
1.460 // of the cycle
1.461 + int first, second;
1.462 if (_state[in_arc] == STATE_LOWER) {
1.463 first = _source[in_arc];
1.464 second = _target[in_arc];
1.465 @@ -1246,25 +1273,32 @@
1.466 }
1.467 delta = _cap[in_arc];
1.468 int result = 0;
1.469 - Value d;
1.470 + Value c, d;
1.471 int e;
1.472
1.473 - // Search the cycle along the path form the first node to the root
1.474 + // Search the cycle form the first node to the join node
1.475 for (int u = first; u != join; u = _parent[u]) {
1.476 e = _pred[u];
1.477 - d = _forward[u] ?
1.478 - _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
1.479 + d = _flow[e];
1.480 + if (_pred_dir[u] == DIR_DOWN) {
1.481 + c = _cap[e];
1.482 + d = c >= MAX ? INF : c - d;
1.483 + }
1.484 if (d < delta) {
1.485 delta = d;
1.486 u_out = u;
1.487 result = 1;
1.488 }
1.489 }
1.490 - // Search the cycle along the path form the second node to the root
1.491 +
1.492 + // Search the cycle form the second node to the join node
1.493 for (int u = second; u != join; u = _parent[u]) {
1.494 e = _pred[u];
1.495 - d = _forward[u] ?
1.496 - (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
1.497 + d = _flow[e];
1.498 + if (_pred_dir[u] == DIR_UP) {
1.499 + c = _cap[e];
1.500 + d = c >= MAX ? INF : c - d;
1.501 + }
1.502 if (d <= delta) {
1.503 delta = d;
1.504 u_out = u;
1.505 @@ -1289,10 +1323,10 @@
1.506 Value val = _state[in_arc] * delta;
1.507 _flow[in_arc] += val;
1.508 for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1.509 - _flow[_pred[u]] += _forward[u] ? -val : val;
1.510 + _flow[_pred[u]] -= _pred_dir[u] * val;
1.511 }
1.512 for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1.513 - _flow[_pred[u]] += _forward[u] ? val : -val;
1.514 + _flow[_pred[u]] += _pred_dir[u] * val;
1.515 }
1.516 }
1.517 // Update the state of the entering and leaving arcs
1.518 @@ -1307,130 +1341,134 @@
1.519
1.520 // Update the tree structure
1.521 void updateTreeStructure() {
1.522 - int u, w;
1.523 int old_rev_thread = _rev_thread[u_out];
1.524 int old_succ_num = _succ_num[u_out];
1.525 int old_last_succ = _last_succ[u_out];
1.526 v_out = _parent[u_out];
1.527
1.528 - u = _last_succ[u_in]; // the last successor of u_in
1.529 - right = _thread[u]; // the node after it
1.530 + // Check if u_in and u_out coincide
1.531 + if (u_in == u_out) {
1.532 + // Update _parent, _pred, _pred_dir
1.533 + _parent[u_in] = v_in;
1.534 + _pred[u_in] = in_arc;
1.535 + _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1.536
1.537 - // Handle the case when old_rev_thread equals to v_in
1.538 - // (it also means that join and v_out coincide)
1.539 - if (old_rev_thread == v_in) {
1.540 - last = _thread[_last_succ[u_out]];
1.541 + // Update _thread and _rev_thread
1.542 + if (_thread[v_in] != u_out) {
1.543 + int after = _thread[old_last_succ];
1.544 + _thread[old_rev_thread] = after;
1.545 + _rev_thread[after] = old_rev_thread;
1.546 + after = _thread[v_in];
1.547 + _thread[v_in] = u_out;
1.548 + _rev_thread[u_out] = v_in;
1.549 + _thread[old_last_succ] = after;
1.550 + _rev_thread[after] = old_last_succ;
1.551 + }
1.552 } else {
1.553 - last = _thread[v_in];
1.554 - }
1.555 + // Handle the case when old_rev_thread equals to v_in
1.556 + // (it also means that join and v_out coincide)
1.557 + int thread_continue = old_rev_thread == v_in ?
1.558 + _thread[old_last_succ] : _thread[v_in];
1.559
1.560 - // Update _thread and _parent along the stem nodes (i.e. the nodes
1.561 - // between u_in and u_out, whose parent have to be changed)
1.562 - _thread[v_in] = stem = u_in;
1.563 - _dirty_revs.clear();
1.564 - _dirty_revs.push_back(v_in);
1.565 - par_stem = v_in;
1.566 - while (stem != u_out) {
1.567 - // Insert the next stem node into the thread list
1.568 - new_stem = _parent[stem];
1.569 - _thread[u] = new_stem;
1.570 - _dirty_revs.push_back(u);
1.571 + // Update _thread and _parent along the stem nodes (i.e. the nodes
1.572 + // between u_in and u_out, whose parent have to be changed)
1.573 + int stem = u_in; // the current stem node
1.574 + int par_stem = v_in; // the new parent of stem
1.575 + int next_stem; // the next stem node
1.576 + int last = _last_succ[u_in]; // the last successor of stem
1.577 + int before, after = _thread[last];
1.578 + _thread[v_in] = u_in;
1.579 + _dirty_revs.clear();
1.580 + _dirty_revs.push_back(v_in);
1.581 + while (stem != u_out) {
1.582 + // Insert the next stem node into the thread list
1.583 + next_stem = _parent[stem];
1.584 + _thread[last] = next_stem;
1.585 + _dirty_revs.push_back(last);
1.586
1.587 - // Remove the subtree of stem from the thread list
1.588 - w = _rev_thread[stem];
1.589 - _thread[w] = right;
1.590 - _rev_thread[right] = w;
1.591 + // Remove the subtree of stem from the thread list
1.592 + before = _rev_thread[stem];
1.593 + _thread[before] = after;
1.594 + _rev_thread[after] = before;
1.595
1.596 - // Change the parent node and shift stem nodes
1.597 - _parent[stem] = par_stem;
1.598 - par_stem = stem;
1.599 - stem = new_stem;
1.600 + // Change the parent node and shift stem nodes
1.601 + _parent[stem] = par_stem;
1.602 + par_stem = stem;
1.603 + stem = next_stem;
1.604
1.605 - // Update u and right
1.606 - u = _last_succ[stem] == _last_succ[par_stem] ?
1.607 - _rev_thread[par_stem] : _last_succ[stem];
1.608 - right = _thread[u];
1.609 - }
1.610 - _parent[u_out] = par_stem;
1.611 - _thread[u] = last;
1.612 - _rev_thread[last] = u;
1.613 - _last_succ[u_out] = u;
1.614 + // Update last and after
1.615 + last = _last_succ[stem] == _last_succ[par_stem] ?
1.616 + _rev_thread[par_stem] : _last_succ[stem];
1.617 + after = _thread[last];
1.618 + }
1.619 + _parent[u_out] = par_stem;
1.620 + _thread[last] = thread_continue;
1.621 + _rev_thread[thread_continue] = last;
1.622 + _last_succ[u_out] = last;
1.623
1.624 - // Remove the subtree of u_out from the thread list except for
1.625 - // the case when old_rev_thread equals to v_in
1.626 - // (it also means that join and v_out coincide)
1.627 - if (old_rev_thread != v_in) {
1.628 - _thread[old_rev_thread] = right;
1.629 - _rev_thread[right] = old_rev_thread;
1.630 - }
1.631 + // Remove the subtree of u_out from the thread list except for
1.632 + // the case when old_rev_thread equals to v_in
1.633 + if (old_rev_thread != v_in) {
1.634 + _thread[old_rev_thread] = after;
1.635 + _rev_thread[after] = old_rev_thread;
1.636 + }
1.637
1.638 - // Update _rev_thread using the new _thread values
1.639 - for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1.640 - u = _dirty_revs[i];
1.641 - _rev_thread[_thread[u]] = u;
1.642 - }
1.643 + // Update _rev_thread using the new _thread values
1.644 + for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1.645 + int u = _dirty_revs[i];
1.646 + _rev_thread[_thread[u]] = u;
1.647 + }
1.648
1.649 - // Update _pred, _forward, _last_succ and _succ_num for the
1.650 - // stem nodes from u_out to u_in
1.651 - int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1.652 - u = u_out;
1.653 - while (u != u_in) {
1.654 - w = _parent[u];
1.655 - _pred[u] = _pred[w];
1.656 - _forward[u] = !_forward[w];
1.657 - tmp_sc += _succ_num[u] - _succ_num[w];
1.658 - _succ_num[u] = tmp_sc;
1.659 - _last_succ[w] = tmp_ls;
1.660 - u = w;
1.661 - }
1.662 - _pred[u_in] = in_arc;
1.663 - _forward[u_in] = (u_in == _source[in_arc]);
1.664 - _succ_num[u_in] = old_succ_num;
1.665 -
1.666 - // Set limits for updating _last_succ form v_in and v_out
1.667 - // towards the root
1.668 - int up_limit_in = -1;
1.669 - int up_limit_out = -1;
1.670 - if (_last_succ[join] == v_in) {
1.671 - up_limit_out = join;
1.672 - } else {
1.673 - up_limit_in = join;
1.674 + // Update _pred, _pred_dir, _last_succ and _succ_num for the
1.675 + // stem nodes from u_out to u_in
1.676 + int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1.677 + for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
1.678 + _pred[u] = _pred[p];
1.679 + _pred_dir[u] = -_pred_dir[p];
1.680 + tmp_sc += _succ_num[u] - _succ_num[p];
1.681 + _succ_num[u] = tmp_sc;
1.682 + _last_succ[p] = tmp_ls;
1.683 + }
1.684 + _pred[u_in] = in_arc;
1.685 + _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1.686 + _succ_num[u_in] = old_succ_num;
1.687 }
1.688
1.689 // Update _last_succ from v_in towards the root
1.690 - for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1.691 - u = _parent[u]) {
1.692 - _last_succ[u] = _last_succ[u_out];
1.693 + int up_limit_out = _last_succ[join] == v_in ? join : -1;
1.694 + int last_succ_out = _last_succ[u_out];
1.695 + for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
1.696 + _last_succ[u] = last_succ_out;
1.697 }
1.698 +
1.699 // Update _last_succ from v_out towards the root
1.700 if (join != old_rev_thread && v_in != old_rev_thread) {
1.701 - for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.702 + for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.703 u = _parent[u]) {
1.704 _last_succ[u] = old_rev_thread;
1.705 }
1.706 - } else {
1.707 - for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.708 + }
1.709 + else if (last_succ_out != old_last_succ) {
1.710 + for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.711 u = _parent[u]) {
1.712 - _last_succ[u] = _last_succ[u_out];
1.713 + _last_succ[u] = last_succ_out;
1.714 }
1.715 }
1.716
1.717 // Update _succ_num from v_in to join
1.718 - for (u = v_in; u != join; u = _parent[u]) {
1.719 + for (int u = v_in; u != join; u = _parent[u]) {
1.720 _succ_num[u] += old_succ_num;
1.721 }
1.722 // Update _succ_num from v_out to join
1.723 - for (u = v_out; u != join; u = _parent[u]) {
1.724 + for (int u = v_out; u != join; u = _parent[u]) {
1.725 _succ_num[u] -= old_succ_num;
1.726 }
1.727 }
1.728
1.729 - // Update potentials
1.730 + // Update potentials in the subtree that has been moved
1.731 void updatePotential() {
1.732 - Cost sigma = _forward[u_in] ?
1.733 - _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1.734 - _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1.735 - // Update potentials in the subtree, which has been moved
1.736 + Cost sigma = _pi[v_in] - _pi[u_in] -
1.737 + _pred_dir[u_in] * _cost[in_arc];
1.738 int end = _thread[_last_succ[u_in]];
1.739 for (int u = u_in; u != end; u = _thread[u]) {
1.740 _pi[u] += sigma;
1.741 @@ -1478,7 +1516,7 @@
1.742 }
1.743 }
1.744 } else {
1.745 - // Find the min. cost incomming arc for each demand node
1.746 + // Find the min. cost incoming arc for each demand node
1.747 for (int i = 0; i != int(demand_nodes.size()); ++i) {
1.748 Node v = demand_nodes[i];
1.749 Cost c, min_cost = std::numeric_limits<Cost>::max();
1.750 @@ -1574,7 +1612,7 @@
1.751 }
1.752
1.753 // Transform the solution and the supply map to the original form
1.754 - if (_have_lower) {
1.755 + if (_has_lower) {
1.756 for (int i = 0; i != _arc_num; ++i) {
1.757 Value c = _lower[i];
1.758 if (c != 0) {