1.1 --- a/lemon/network_simplex.h Mon May 11 16:38:21 2009 +0100
1.2 +++ b/lemon/network_simplex.h Tue May 12 12:06:40 2009 +0200
1.3 @@ -19,7 +19,7 @@
1.4 #ifndef LEMON_NETWORK_SIMPLEX_H
1.5 #define LEMON_NETWORK_SIMPLEX_H
1.6
1.7 -/// \ingroup min_cost_flow
1.8 +/// \ingroup min_cost_flow_algs
1.9 ///
1.10 /// \file
1.11 /// \brief Network Simplex algorithm for finding a minimum cost flow.
1.12 @@ -33,7 +33,7 @@
1.13
1.14 namespace lemon {
1.15
1.16 - /// \addtogroup min_cost_flow
1.17 + /// \addtogroup min_cost_flow_algs
1.18 /// @{
1.19
1.20 /// \brief Implementation of the primal Network Simplex algorithm
1.21 @@ -102,50 +102,16 @@
1.22 /// i.e. the direction of the inequalities in the supply/demand
1.23 /// constraints of the \ref min_cost_flow "minimum cost flow problem".
1.24 ///
1.25 - /// The default supply type is \c GEQ, since this form is supported
1.26 - /// by other minimum cost flow algorithms and the \ref Circulation
1.27 - /// algorithm, as well.
1.28 - /// The \c LEQ problem type can be selected using the \ref supplyType()
1.29 - /// function.
1.30 - ///
1.31 - /// Note that the equality form is a special case of both supply types.
1.32 + /// The default supply type is \c GEQ, the \c LEQ type can be
1.33 + /// selected using \ref supplyType().
1.34 + /// The equality form is a special case of both supply types.
1.35 enum SupplyType {
1.36 -
1.37 /// This option means that there are <em>"greater or equal"</em>
1.38 - /// supply/demand constraints in the definition, i.e. the exact
1.39 - /// formulation of the problem is the following.
1.40 - /**
1.41 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.42 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.43 - sup(u) \quad \forall u\in V \f]
1.44 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.45 - */
1.46 - /// It means that the total demand must be greater or equal to the
1.47 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.48 - /// negative) and all the supplies have to be carried out from
1.49 - /// the supply nodes, but there could be demands that are not
1.50 - /// satisfied.
1.51 + /// supply/demand constraints in the definition of the problem.
1.52 GEQ,
1.53 - /// It is just an alias for the \c GEQ option.
1.54 - CARRY_SUPPLIES = GEQ,
1.55 -
1.56 /// This option means that there are <em>"less or equal"</em>
1.57 - /// supply/demand constraints in the definition, i.e. the exact
1.58 - /// formulation of the problem is the following.
1.59 - /**
1.60 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.61 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
1.62 - sup(u) \quad \forall u\in V \f]
1.63 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.64 - */
1.65 - /// It means that the total demand must be less or equal to the
1.66 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.67 - /// positive) and all the demands have to be satisfied, but there
1.68 - /// could be supplies that are not carried out from the supply
1.69 - /// nodes.
1.70 - LEQ,
1.71 - /// It is just an alias for the \c LEQ option.
1.72 - SATISFY_DEMANDS = LEQ
1.73 + /// supply/demand constraints in the definition of the problem.
1.74 + LEQ
1.75 };
1.76
1.77 /// \brief Constants for selecting the pivot rule.
1.78 @@ -215,6 +181,8 @@
1.79 const GR &_graph;
1.80 int _node_num;
1.81 int _arc_num;
1.82 + int _all_arc_num;
1.83 + int _search_arc_num;
1.84
1.85 // Parameters of the problem
1.86 bool _have_lower;
1.87 @@ -277,7 +245,7 @@
1.88 const IntVector &_state;
1.89 const CostVector &_pi;
1.90 int &_in_arc;
1.91 - int _arc_num;
1.92 + int _search_arc_num;
1.93
1.94 // Pivot rule data
1.95 int _next_arc;
1.96 @@ -288,13 +256,14 @@
1.97 FirstEligiblePivotRule(NetworkSimplex &ns) :
1.98 _source(ns._source), _target(ns._target),
1.99 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.100 - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
1.101 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.102 + _next_arc(0)
1.103 {}
1.104
1.105 // Find next entering arc
1.106 bool findEnteringArc() {
1.107 Cost c;
1.108 - for (int e = _next_arc; e < _arc_num; ++e) {
1.109 + for (int e = _next_arc; e < _search_arc_num; ++e) {
1.110 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.111 if (c < 0) {
1.112 _in_arc = e;
1.113 @@ -328,7 +297,7 @@
1.114 const IntVector &_state;
1.115 const CostVector &_pi;
1.116 int &_in_arc;
1.117 - int _arc_num;
1.118 + int _search_arc_num;
1.119
1.120 public:
1.121
1.122 @@ -336,13 +305,13 @@
1.123 BestEligiblePivotRule(NetworkSimplex &ns) :
1.124 _source(ns._source), _target(ns._target),
1.125 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.126 - _in_arc(ns.in_arc), _arc_num(ns._arc_num)
1.127 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
1.128 {}
1.129
1.130 // Find next entering arc
1.131 bool findEnteringArc() {
1.132 Cost c, min = 0;
1.133 - for (int e = 0; e < _arc_num; ++e) {
1.134 + for (int e = 0; e < _search_arc_num; ++e) {
1.135 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.136 if (c < min) {
1.137 min = c;
1.138 @@ -367,7 +336,7 @@
1.139 const IntVector &_state;
1.140 const CostVector &_pi;
1.141 int &_in_arc;
1.142 - int _arc_num;
1.143 + int _search_arc_num;
1.144
1.145 // Pivot rule data
1.146 int _block_size;
1.147 @@ -379,14 +348,15 @@
1.148 BlockSearchPivotRule(NetworkSimplex &ns) :
1.149 _source(ns._source), _target(ns._target),
1.150 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.151 - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
1.152 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.153 + _next_arc(0)
1.154 {
1.155 // The main parameters of the pivot rule
1.156 - const double BLOCK_SIZE_FACTOR = 2.0;
1.157 + const double BLOCK_SIZE_FACTOR = 0.5;
1.158 const int MIN_BLOCK_SIZE = 10;
1.159
1.160 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.161 - std::sqrt(double(_arc_num))),
1.162 + std::sqrt(double(_search_arc_num))),
1.163 MIN_BLOCK_SIZE );
1.164 }
1.165
1.166 @@ -395,7 +365,7 @@
1.167 Cost c, min = 0;
1.168 int cnt = _block_size;
1.169 int e, min_arc = _next_arc;
1.170 - for (e = _next_arc; e < _arc_num; ++e) {
1.171 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.172 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.173 if (c < min) {
1.174 min = c;
1.175 @@ -440,7 +410,7 @@
1.176 const IntVector &_state;
1.177 const CostVector &_pi;
1.178 int &_in_arc;
1.179 - int _arc_num;
1.180 + int _search_arc_num;
1.181
1.182 // Pivot rule data
1.183 IntVector _candidates;
1.184 @@ -454,7 +424,8 @@
1.185 CandidateListPivotRule(NetworkSimplex &ns) :
1.186 _source(ns._source), _target(ns._target),
1.187 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.188 - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
1.189 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.190 + _next_arc(0)
1.191 {
1.192 // The main parameters of the pivot rule
1.193 const double LIST_LENGTH_FACTOR = 1.0;
1.194 @@ -463,7 +434,7 @@
1.195 const int MIN_MINOR_LIMIT = 3;
1.196
1.197 _list_length = std::max( int(LIST_LENGTH_FACTOR *
1.198 - std::sqrt(double(_arc_num))),
1.199 + std::sqrt(double(_search_arc_num))),
1.200 MIN_LIST_LENGTH );
1.201 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
1.202 MIN_MINOR_LIMIT );
1.203 @@ -500,7 +471,7 @@
1.204 // Major iteration: build a new candidate list
1.205 min = 0;
1.206 _curr_length = 0;
1.207 - for (e = _next_arc; e < _arc_num; ++e) {
1.208 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.209 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.210 if (c < 0) {
1.211 _candidates[_curr_length++] = e;
1.212 @@ -546,7 +517,7 @@
1.213 const IntVector &_state;
1.214 const CostVector &_pi;
1.215 int &_in_arc;
1.216 - int _arc_num;
1.217 + int _search_arc_num;
1.218
1.219 // Pivot rule data
1.220 int _block_size, _head_length, _curr_length;
1.221 @@ -574,8 +545,8 @@
1.222 AlteringListPivotRule(NetworkSimplex &ns) :
1.223 _source(ns._source), _target(ns._target),
1.224 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.225 - _in_arc(ns.in_arc), _arc_num(ns._arc_num),
1.226 - _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
1.227 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.228 + _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.229 {
1.230 // The main parameters of the pivot rule
1.231 const double BLOCK_SIZE_FACTOR = 1.5;
1.232 @@ -584,7 +555,7 @@
1.233 const int MIN_HEAD_LENGTH = 3;
1.234
1.235 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.236 - std::sqrt(double(_arc_num))),
1.237 + std::sqrt(double(_search_arc_num))),
1.238 MIN_BLOCK_SIZE );
1.239 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
1.240 MIN_HEAD_LENGTH );
1.241 @@ -610,7 +581,7 @@
1.242 int last_arc = 0;
1.243 int limit = _head_length;
1.244
1.245 - for (int e = _next_arc; e < _arc_num; ++e) {
1.246 + for (int e = _next_arc; e < _search_arc_num; ++e) {
1.247 _cand_cost[e] = _state[e] *
1.248 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.249 if (_cand_cost[e] < 0) {
1.250 @@ -678,17 +649,17 @@
1.251 _node_num = countNodes(_graph);
1.252 _arc_num = countArcs(_graph);
1.253 int all_node_num = _node_num + 1;
1.254 - int all_arc_num = _arc_num + _node_num;
1.255 + int max_arc_num = _arc_num + 2 * _node_num;
1.256
1.257 - _source.resize(all_arc_num);
1.258 - _target.resize(all_arc_num);
1.259 + _source.resize(max_arc_num);
1.260 + _target.resize(max_arc_num);
1.261
1.262 - _lower.resize(all_arc_num);
1.263 - _upper.resize(all_arc_num);
1.264 - _cap.resize(all_arc_num);
1.265 - _cost.resize(all_arc_num);
1.266 + _lower.resize(_arc_num);
1.267 + _upper.resize(_arc_num);
1.268 + _cap.resize(max_arc_num);
1.269 + _cost.resize(max_arc_num);
1.270 _supply.resize(all_node_num);
1.271 - _flow.resize(all_arc_num);
1.272 + _flow.resize(max_arc_num);
1.273 _pi.resize(all_node_num);
1.274
1.275 _parent.resize(all_node_num);
1.276 @@ -698,7 +669,7 @@
1.277 _rev_thread.resize(all_node_num);
1.278 _succ_num.resize(all_node_num);
1.279 _last_succ.resize(all_node_num);
1.280 - _state.resize(all_arc_num);
1.281 + _state.resize(max_arc_num);
1.282
1.283 // Copy the graph (store the arcs in a mixed order)
1.284 int i = 0;
1.285 @@ -1069,7 +1040,7 @@
1.286 // Initialize artifical cost
1.287 Cost ART_COST;
1.288 if (std::numeric_limits<Cost>::is_exact) {
1.289 - ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
1.290 + ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
1.291 } else {
1.292 ART_COST = std::numeric_limits<Cost>::min();
1.293 for (int i = 0; i != _arc_num; ++i) {
1.294 @@ -1093,29 +1064,121 @@
1.295 _succ_num[_root] = _node_num + 1;
1.296 _last_succ[_root] = _root - 1;
1.297 _supply[_root] = -_sum_supply;
1.298 - _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
1.299 + _pi[_root] = 0;
1.300
1.301 // Add artificial arcs and initialize the spanning tree data structure
1.302 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.303 - _parent[u] = _root;
1.304 - _pred[u] = e;
1.305 - _thread[u] = u + 1;
1.306 - _rev_thread[u + 1] = u;
1.307 - _succ_num[u] = 1;
1.308 - _last_succ[u] = u;
1.309 - _cost[e] = ART_COST;
1.310 - _cap[e] = INF;
1.311 - _state[e] = STATE_TREE;
1.312 - if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
1.313 - _flow[e] = _supply[u];
1.314 - _forward[u] = true;
1.315 - _pi[u] = -ART_COST + _pi[_root];
1.316 - } else {
1.317 - _flow[e] = -_supply[u];
1.318 - _forward[u] = false;
1.319 - _pi[u] = ART_COST + _pi[_root];
1.320 + if (_sum_supply == 0) {
1.321 + // EQ supply constraints
1.322 + _search_arc_num = _arc_num;
1.323 + _all_arc_num = _arc_num + _node_num;
1.324 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.325 + _parent[u] = _root;
1.326 + _pred[u] = e;
1.327 + _thread[u] = u + 1;
1.328 + _rev_thread[u + 1] = u;
1.329 + _succ_num[u] = 1;
1.330 + _last_succ[u] = u;
1.331 + _cap[e] = INF;
1.332 + _state[e] = STATE_TREE;
1.333 + if (_supply[u] >= 0) {
1.334 + _forward[u] = true;
1.335 + _pi[u] = 0;
1.336 + _source[e] = u;
1.337 + _target[e] = _root;
1.338 + _flow[e] = _supply[u];
1.339 + _cost[e] = 0;
1.340 + } else {
1.341 + _forward[u] = false;
1.342 + _pi[u] = ART_COST;
1.343 + _source[e] = _root;
1.344 + _target[e] = u;
1.345 + _flow[e] = -_supply[u];
1.346 + _cost[e] = ART_COST;
1.347 + }
1.348 }
1.349 }
1.350 + else if (_sum_supply > 0) {
1.351 + // LEQ supply constraints
1.352 + _search_arc_num = _arc_num + _node_num;
1.353 + int f = _arc_num + _node_num;
1.354 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.355 + _parent[u] = _root;
1.356 + _thread[u] = u + 1;
1.357 + _rev_thread[u + 1] = u;
1.358 + _succ_num[u] = 1;
1.359 + _last_succ[u] = u;
1.360 + if (_supply[u] >= 0) {
1.361 + _forward[u] = true;
1.362 + _pi[u] = 0;
1.363 + _pred[u] = e;
1.364 + _source[e] = u;
1.365 + _target[e] = _root;
1.366 + _cap[e] = INF;
1.367 + _flow[e] = _supply[u];
1.368 + _cost[e] = 0;
1.369 + _state[e] = STATE_TREE;
1.370 + } else {
1.371 + _forward[u] = false;
1.372 + _pi[u] = ART_COST;
1.373 + _pred[u] = f;
1.374 + _source[f] = _root;
1.375 + _target[f] = u;
1.376 + _cap[f] = INF;
1.377 + _flow[f] = -_supply[u];
1.378 + _cost[f] = ART_COST;
1.379 + _state[f] = STATE_TREE;
1.380 + _source[e] = u;
1.381 + _target[e] = _root;
1.382 + _cap[e] = INF;
1.383 + _flow[e] = 0;
1.384 + _cost[e] = 0;
1.385 + _state[e] = STATE_LOWER;
1.386 + ++f;
1.387 + }
1.388 + }
1.389 + _all_arc_num = f;
1.390 + }
1.391 + else {
1.392 + // GEQ supply constraints
1.393 + _search_arc_num = _arc_num + _node_num;
1.394 + int f = _arc_num + _node_num;
1.395 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.396 + _parent[u] = _root;
1.397 + _thread[u] = u + 1;
1.398 + _rev_thread[u + 1] = u;
1.399 + _succ_num[u] = 1;
1.400 + _last_succ[u] = u;
1.401 + if (_supply[u] <= 0) {
1.402 + _forward[u] = false;
1.403 + _pi[u] = 0;
1.404 + _pred[u] = e;
1.405 + _source[e] = _root;
1.406 + _target[e] = u;
1.407 + _cap[e] = INF;
1.408 + _flow[e] = -_supply[u];
1.409 + _cost[e] = 0;
1.410 + _state[e] = STATE_TREE;
1.411 + } else {
1.412 + _forward[u] = true;
1.413 + _pi[u] = -ART_COST;
1.414 + _pred[u] = f;
1.415 + _source[f] = u;
1.416 + _target[f] = _root;
1.417 + _cap[f] = INF;
1.418 + _flow[f] = _supply[u];
1.419 + _state[f] = STATE_TREE;
1.420 + _cost[f] = ART_COST;
1.421 + _source[e] = _root;
1.422 + _target[e] = u;
1.423 + _cap[e] = INF;
1.424 + _flow[e] = 0;
1.425 + _cost[e] = 0;
1.426 + _state[e] = STATE_LOWER;
1.427 + ++f;
1.428 + }
1.429 + }
1.430 + _all_arc_num = f;
1.431 + }
1.432
1.433 return true;
1.434 }
1.435 @@ -1374,20 +1437,8 @@
1.436 }
1.437
1.438 // Check feasibility
1.439 - if (_sum_supply < 0) {
1.440 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.441 - if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1.442 - }
1.443 - }
1.444 - else if (_sum_supply > 0) {
1.445 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.446 - if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1.447 - }
1.448 - }
1.449 - else {
1.450 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.451 - if (_flow[e] != 0) return INFEASIBLE;
1.452 - }
1.453 + for (int e = _search_arc_num; e != _all_arc_num; ++e) {
1.454 + if (_flow[e] != 0) return INFEASIBLE;
1.455 }
1.456
1.457 // Transform the solution and the supply map to the original form
1.458 @@ -1401,6 +1452,30 @@
1.459 }
1.460 }
1.461 }
1.462 +
1.463 + // Shift potentials to meet the requirements of the GEQ/LEQ type
1.464 + // optimality conditions
1.465 + if (_sum_supply == 0) {
1.466 + if (_stype == GEQ) {
1.467 + Cost max_pot = std::numeric_limits<Cost>::min();
1.468 + for (int i = 0; i != _node_num; ++i) {
1.469 + if (_pi[i] > max_pot) max_pot = _pi[i];
1.470 + }
1.471 + if (max_pot > 0) {
1.472 + for (int i = 0; i != _node_num; ++i)
1.473 + _pi[i] -= max_pot;
1.474 + }
1.475 + } else {
1.476 + Cost min_pot = std::numeric_limits<Cost>::max();
1.477 + for (int i = 0; i != _node_num; ++i) {
1.478 + if (_pi[i] < min_pot) min_pot = _pi[i];
1.479 + }
1.480 + if (min_pot < 0) {
1.481 + for (int i = 0; i != _node_num; ++i)
1.482 + _pi[i] -= min_pot;
1.483 + }
1.484 + }
1.485 + }
1.486
1.487 return OPTIMAL;
1.488 }