1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/network_simplex.h Thu Dec 10 17:05:35 2009 +0100
1.3 @@ -0,0 +1,1489 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2009
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_NETWORK_SIMPLEX_H
1.23 +#define LEMON_NETWORK_SIMPLEX_H
1.24 +
1.25 +/// \ingroup min_cost_flow_algs
1.26 +///
1.27 +/// \file
1.28 +/// \brief Network Simplex algorithm for finding a minimum cost flow.
1.29 +
1.30 +#include <vector>
1.31 +#include <limits>
1.32 +#include <algorithm>
1.33 +
1.34 +#include <lemon/core.h>
1.35 +#include <lemon/math.h>
1.36 +
1.37 +namespace lemon {
1.38 +
1.39 + /// \addtogroup min_cost_flow_algs
1.40 + /// @{
1.41 +
1.42 + /// \brief Implementation of the primal Network Simplex algorithm
1.43 + /// for finding a \ref min_cost_flow "minimum cost flow".
1.44 + ///
1.45 + /// \ref NetworkSimplex implements the primal Network Simplex algorithm
1.46 + /// for finding a \ref min_cost_flow "minimum cost flow".
1.47 + /// This algorithm is a specialized version of the linear programming
1.48 + /// simplex method directly for the minimum cost flow problem.
1.49 + /// It is one of the most efficient solution methods.
1.50 + ///
1.51 + /// In general this class is the fastest implementation available
1.52 + /// in LEMON for the minimum cost flow problem.
1.53 + /// Moreover it supports both directions of the supply/demand inequality
1.54 + /// constraints. For more information see \ref SupplyType.
1.55 + ///
1.56 + /// Most of the parameters of the problem (except for the digraph)
1.57 + /// can be given using separate functions, and the algorithm can be
1.58 + /// executed using the \ref run() function. If some parameters are not
1.59 + /// specified, then default values will be used.
1.60 + ///
1.61 + /// \tparam GR The digraph type the algorithm runs on.
1.62 + /// \tparam V The value type used for flow amounts, capacity bounds
1.63 + /// and supply values in the algorithm. By default it is \c int.
1.64 + /// \tparam C The value type used for costs and potentials in the
1.65 + /// algorithm. By default it is the same as \c V.
1.66 + ///
1.67 + /// \warning Both value types must be signed and all input data must
1.68 + /// be integer.
1.69 + ///
1.70 + /// \note %NetworkSimplex provides five different pivot rule
1.71 + /// implementations, from which the most efficient one is used
1.72 + /// by default. For more information see \ref PivotRule.
1.73 + template <typename GR, typename V = int, typename C = V>
1.74 + class NetworkSimplex
1.75 + {
1.76 + public:
1.77 +
1.78 + /// The type of the flow amounts, capacity bounds and supply values
1.79 + typedef V Value;
1.80 + /// The type of the arc costs
1.81 + typedef C Cost;
1.82 +
1.83 + public:
1.84 +
1.85 + /// \brief Problem type constants for the \c run() function.
1.86 + ///
1.87 + /// Enum type containing the problem type constants that can be
1.88 + /// returned by the \ref run() function of the algorithm.
1.89 + enum ProblemType {
1.90 + /// The problem has no feasible solution (flow).
1.91 + INFEASIBLE,
1.92 + /// The problem has optimal solution (i.e. it is feasible and
1.93 + /// bounded), and the algorithm has found optimal flow and node
1.94 + /// potentials (primal and dual solutions).
1.95 + OPTIMAL,
1.96 + /// The objective function of the problem is unbounded, i.e.
1.97 + /// there is a directed cycle having negative total cost and
1.98 + /// infinite upper bound.
1.99 + UNBOUNDED
1.100 + };
1.101 +
1.102 + /// \brief Constants for selecting the type of the supply constraints.
1.103 + ///
1.104 + /// Enum type containing constants for selecting the supply type,
1.105 + /// i.e. the direction of the inequalities in the supply/demand
1.106 + /// constraints of the \ref min_cost_flow "minimum cost flow problem".
1.107 + ///
1.108 + /// The default supply type is \c GEQ, the \c LEQ type can be
1.109 + /// selected using \ref supplyType().
1.110 + /// The equality form is a special case of both supply types.
1.111 + enum SupplyType {
1.112 + /// This option means that there are <em>"greater or equal"</em>
1.113 + /// supply/demand constraints in the definition of the problem.
1.114 + GEQ,
1.115 + /// This option means that there are <em>"less or equal"</em>
1.116 + /// supply/demand constraints in the definition of the problem.
1.117 + LEQ
1.118 + };
1.119 +
1.120 + /// \brief Constants for selecting the pivot rule.
1.121 + ///
1.122 + /// Enum type containing constants for selecting the pivot rule for
1.123 + /// the \ref run() function.
1.124 + ///
1.125 + /// \ref NetworkSimplex provides five different pivot rule
1.126 + /// implementations that significantly affect the running time
1.127 + /// of the algorithm.
1.128 + /// By default \ref BLOCK_SEARCH "Block Search" is used, which
1.129 + /// proved to be the most efficient and the most robust on various
1.130 + /// test inputs according to our benchmark tests.
1.131 + /// However another pivot rule can be selected using the \ref run()
1.132 + /// function with the proper parameter.
1.133 + enum PivotRule {
1.134 +
1.135 + /// The First Eligible pivot rule.
1.136 + /// The next eligible arc is selected in a wraparound fashion
1.137 + /// in every iteration.
1.138 + FIRST_ELIGIBLE,
1.139 +
1.140 + /// The Best Eligible pivot rule.
1.141 + /// The best eligible arc is selected in every iteration.
1.142 + BEST_ELIGIBLE,
1.143 +
1.144 + /// The Block Search pivot rule.
1.145 + /// A specified number of arcs are examined in every iteration
1.146 + /// in a wraparound fashion and the best eligible arc is selected
1.147 + /// from this block.
1.148 + BLOCK_SEARCH,
1.149 +
1.150 + /// The Candidate List pivot rule.
1.151 + /// In a major iteration a candidate list is built from eligible arcs
1.152 + /// in a wraparound fashion and in the following minor iterations
1.153 + /// the best eligible arc is selected from this list.
1.154 + CANDIDATE_LIST,
1.155 +
1.156 + /// The Altering Candidate List pivot rule.
1.157 + /// It is a modified version of the Candidate List method.
1.158 + /// It keeps only the several best eligible arcs from the former
1.159 + /// candidate list and extends this list in every iteration.
1.160 + ALTERING_LIST
1.161 + };
1.162 +
1.163 + private:
1.164 +
1.165 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.166 +
1.167 + typedef std::vector<Arc> ArcVector;
1.168 + typedef std::vector<Node> NodeVector;
1.169 + typedef std::vector<int> IntVector;
1.170 + typedef std::vector<bool> BoolVector;
1.171 + typedef std::vector<Value> ValueVector;
1.172 + typedef std::vector<Cost> CostVector;
1.173 +
1.174 + // State constants for arcs
1.175 + enum ArcStateEnum {
1.176 + STATE_UPPER = -1,
1.177 + STATE_TREE = 0,
1.178 + STATE_LOWER = 1
1.179 + };
1.180 +
1.181 + private:
1.182 +
1.183 + // Data related to the underlying digraph
1.184 + const GR &_graph;
1.185 + int _node_num;
1.186 + int _arc_num;
1.187 + int _all_arc_num;
1.188 + int _search_arc_num;
1.189 +
1.190 + // Parameters of the problem
1.191 + bool _have_lower;
1.192 + SupplyType _stype;
1.193 + Value _sum_supply;
1.194 +
1.195 + // Data structures for storing the digraph
1.196 + IntNodeMap _node_id;
1.197 + IntArcMap _arc_id;
1.198 + IntVector _source;
1.199 + IntVector _target;
1.200 +
1.201 + // Node and arc data
1.202 + ValueVector _lower;
1.203 + ValueVector _upper;
1.204 + ValueVector _cap;
1.205 + CostVector _cost;
1.206 + ValueVector _supply;
1.207 + ValueVector _flow;
1.208 + CostVector _pi;
1.209 +
1.210 + // Data for storing the spanning tree structure
1.211 + IntVector _parent;
1.212 + IntVector _pred;
1.213 + IntVector _thread;
1.214 + IntVector _rev_thread;
1.215 + IntVector _succ_num;
1.216 + IntVector _last_succ;
1.217 + IntVector _dirty_revs;
1.218 + BoolVector _forward;
1.219 + IntVector _state;
1.220 + int _root;
1.221 +
1.222 + // Temporary data used in the current pivot iteration
1.223 + int in_arc, join, u_in, v_in, u_out, v_out;
1.224 + int first, second, right, last;
1.225 + int stem, par_stem, new_stem;
1.226 + Value delta;
1.227 +
1.228 + public:
1.229 +
1.230 + /// \brief Constant for infinite upper bounds (capacities).
1.231 + ///
1.232 + /// Constant for infinite upper bounds (capacities).
1.233 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.234 + /// \c std::numeric_limits<Value>::max() otherwise.
1.235 + const Value INF;
1.236 +
1.237 + private:
1.238 +
1.239 + // Implementation of the First Eligible pivot rule
1.240 + class FirstEligiblePivotRule
1.241 + {
1.242 + private:
1.243 +
1.244 + // References to the NetworkSimplex class
1.245 + const IntVector &_source;
1.246 + const IntVector &_target;
1.247 + const CostVector &_cost;
1.248 + const IntVector &_state;
1.249 + const CostVector &_pi;
1.250 + int &_in_arc;
1.251 + int _search_arc_num;
1.252 +
1.253 + // Pivot rule data
1.254 + int _next_arc;
1.255 +
1.256 + public:
1.257 +
1.258 + // Constructor
1.259 + FirstEligiblePivotRule(NetworkSimplex &ns) :
1.260 + _source(ns._source), _target(ns._target),
1.261 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.262 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.263 + _next_arc(0)
1.264 + {}
1.265 +
1.266 + // Find next entering arc
1.267 + bool findEnteringArc() {
1.268 + Cost c;
1.269 + for (int e = _next_arc; e < _search_arc_num; ++e) {
1.270 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.271 + if (c < 0) {
1.272 + _in_arc = e;
1.273 + _next_arc = e + 1;
1.274 + return true;
1.275 + }
1.276 + }
1.277 + for (int e = 0; e < _next_arc; ++e) {
1.278 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.279 + if (c < 0) {
1.280 + _in_arc = e;
1.281 + _next_arc = e + 1;
1.282 + return true;
1.283 + }
1.284 + }
1.285 + return false;
1.286 + }
1.287 +
1.288 + }; //class FirstEligiblePivotRule
1.289 +
1.290 +
1.291 + // Implementation of the Best Eligible pivot rule
1.292 + class BestEligiblePivotRule
1.293 + {
1.294 + private:
1.295 +
1.296 + // References to the NetworkSimplex class
1.297 + const IntVector &_source;
1.298 + const IntVector &_target;
1.299 + const CostVector &_cost;
1.300 + const IntVector &_state;
1.301 + const CostVector &_pi;
1.302 + int &_in_arc;
1.303 + int _search_arc_num;
1.304 +
1.305 + public:
1.306 +
1.307 + // Constructor
1.308 + BestEligiblePivotRule(NetworkSimplex &ns) :
1.309 + _source(ns._source), _target(ns._target),
1.310 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.311 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
1.312 + {}
1.313 +
1.314 + // Find next entering arc
1.315 + bool findEnteringArc() {
1.316 + Cost c, min = 0;
1.317 + for (int e = 0; e < _search_arc_num; ++e) {
1.318 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.319 + if (c < min) {
1.320 + min = c;
1.321 + _in_arc = e;
1.322 + }
1.323 + }
1.324 + return min < 0;
1.325 + }
1.326 +
1.327 + }; //class BestEligiblePivotRule
1.328 +
1.329 +
1.330 + // Implementation of the Block Search pivot rule
1.331 + class BlockSearchPivotRule
1.332 + {
1.333 + private:
1.334 +
1.335 + // References to the NetworkSimplex class
1.336 + const IntVector &_source;
1.337 + const IntVector &_target;
1.338 + const CostVector &_cost;
1.339 + const IntVector &_state;
1.340 + const CostVector &_pi;
1.341 + int &_in_arc;
1.342 + int _search_arc_num;
1.343 +
1.344 + // Pivot rule data
1.345 + int _block_size;
1.346 + int _next_arc;
1.347 +
1.348 + public:
1.349 +
1.350 + // Constructor
1.351 + BlockSearchPivotRule(NetworkSimplex &ns) :
1.352 + _source(ns._source), _target(ns._target),
1.353 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.354 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.355 + _next_arc(0)
1.356 + {
1.357 + // The main parameters of the pivot rule
1.358 + const double BLOCK_SIZE_FACTOR = 0.5;
1.359 + const int MIN_BLOCK_SIZE = 10;
1.360 +
1.361 + _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.362 + std::sqrt(double(_search_arc_num))),
1.363 + MIN_BLOCK_SIZE );
1.364 + }
1.365 +
1.366 + // Find next entering arc
1.367 + bool findEnteringArc() {
1.368 + Cost c, min = 0;
1.369 + int cnt = _block_size;
1.370 + int e, min_arc = _next_arc;
1.371 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.372 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.373 + if (c < min) {
1.374 + min = c;
1.375 + min_arc = e;
1.376 + }
1.377 + if (--cnt == 0) {
1.378 + if (min < 0) break;
1.379 + cnt = _block_size;
1.380 + }
1.381 + }
1.382 + if (min == 0 || cnt > 0) {
1.383 + for (e = 0; e < _next_arc; ++e) {
1.384 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.385 + if (c < min) {
1.386 + min = c;
1.387 + min_arc = e;
1.388 + }
1.389 + if (--cnt == 0) {
1.390 + if (min < 0) break;
1.391 + cnt = _block_size;
1.392 + }
1.393 + }
1.394 + }
1.395 + if (min >= 0) return false;
1.396 + _in_arc = min_arc;
1.397 + _next_arc = e;
1.398 + return true;
1.399 + }
1.400 +
1.401 + }; //class BlockSearchPivotRule
1.402 +
1.403 +
1.404 + // Implementation of the Candidate List pivot rule
1.405 + class CandidateListPivotRule
1.406 + {
1.407 + private:
1.408 +
1.409 + // References to the NetworkSimplex class
1.410 + const IntVector &_source;
1.411 + const IntVector &_target;
1.412 + const CostVector &_cost;
1.413 + const IntVector &_state;
1.414 + const CostVector &_pi;
1.415 + int &_in_arc;
1.416 + int _search_arc_num;
1.417 +
1.418 + // Pivot rule data
1.419 + IntVector _candidates;
1.420 + int _list_length, _minor_limit;
1.421 + int _curr_length, _minor_count;
1.422 + int _next_arc;
1.423 +
1.424 + public:
1.425 +
1.426 + /// Constructor
1.427 + CandidateListPivotRule(NetworkSimplex &ns) :
1.428 + _source(ns._source), _target(ns._target),
1.429 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.430 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.431 + _next_arc(0)
1.432 + {
1.433 + // The main parameters of the pivot rule
1.434 + const double LIST_LENGTH_FACTOR = 1.0;
1.435 + const int MIN_LIST_LENGTH = 10;
1.436 + const double MINOR_LIMIT_FACTOR = 0.1;
1.437 + const int MIN_MINOR_LIMIT = 3;
1.438 +
1.439 + _list_length = std::max( int(LIST_LENGTH_FACTOR *
1.440 + std::sqrt(double(_search_arc_num))),
1.441 + MIN_LIST_LENGTH );
1.442 + _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
1.443 + MIN_MINOR_LIMIT );
1.444 + _curr_length = _minor_count = 0;
1.445 + _candidates.resize(_list_length);
1.446 + }
1.447 +
1.448 + /// Find next entering arc
1.449 + bool findEnteringArc() {
1.450 + Cost min, c;
1.451 + int e, min_arc = _next_arc;
1.452 + if (_curr_length > 0 && _minor_count < _minor_limit) {
1.453 + // Minor iteration: select the best eligible arc from the
1.454 + // current candidate list
1.455 + ++_minor_count;
1.456 + min = 0;
1.457 + for (int i = 0; i < _curr_length; ++i) {
1.458 + e = _candidates[i];
1.459 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.460 + if (c < min) {
1.461 + min = c;
1.462 + min_arc = e;
1.463 + }
1.464 + if (c >= 0) {
1.465 + _candidates[i--] = _candidates[--_curr_length];
1.466 + }
1.467 + }
1.468 + if (min < 0) {
1.469 + _in_arc = min_arc;
1.470 + return true;
1.471 + }
1.472 + }
1.473 +
1.474 + // Major iteration: build a new candidate list
1.475 + min = 0;
1.476 + _curr_length = 0;
1.477 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.478 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.479 + if (c < 0) {
1.480 + _candidates[_curr_length++] = e;
1.481 + if (c < min) {
1.482 + min = c;
1.483 + min_arc = e;
1.484 + }
1.485 + if (_curr_length == _list_length) break;
1.486 + }
1.487 + }
1.488 + if (_curr_length < _list_length) {
1.489 + for (e = 0; e < _next_arc; ++e) {
1.490 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.491 + if (c < 0) {
1.492 + _candidates[_curr_length++] = e;
1.493 + if (c < min) {
1.494 + min = c;
1.495 + min_arc = e;
1.496 + }
1.497 + if (_curr_length == _list_length) break;
1.498 + }
1.499 + }
1.500 + }
1.501 + if (_curr_length == 0) return false;
1.502 + _minor_count = 1;
1.503 + _in_arc = min_arc;
1.504 + _next_arc = e;
1.505 + return true;
1.506 + }
1.507 +
1.508 + }; //class CandidateListPivotRule
1.509 +
1.510 +
1.511 + // Implementation of the Altering Candidate List pivot rule
1.512 + class AlteringListPivotRule
1.513 + {
1.514 + private:
1.515 +
1.516 + // References to the NetworkSimplex class
1.517 + const IntVector &_source;
1.518 + const IntVector &_target;
1.519 + const CostVector &_cost;
1.520 + const IntVector &_state;
1.521 + const CostVector &_pi;
1.522 + int &_in_arc;
1.523 + int _search_arc_num;
1.524 +
1.525 + // Pivot rule data
1.526 + int _block_size, _head_length, _curr_length;
1.527 + int _next_arc;
1.528 + IntVector _candidates;
1.529 + CostVector _cand_cost;
1.530 +
1.531 + // Functor class to compare arcs during sort of the candidate list
1.532 + class SortFunc
1.533 + {
1.534 + private:
1.535 + const CostVector &_map;
1.536 + public:
1.537 + SortFunc(const CostVector &map) : _map(map) {}
1.538 + bool operator()(int left, int right) {
1.539 + return _map[left] > _map[right];
1.540 + }
1.541 + };
1.542 +
1.543 + SortFunc _sort_func;
1.544 +
1.545 + public:
1.546 +
1.547 + // Constructor
1.548 + AlteringListPivotRule(NetworkSimplex &ns) :
1.549 + _source(ns._source), _target(ns._target),
1.550 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.551 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
1.552 + _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.553 + {
1.554 + // The main parameters of the pivot rule
1.555 + const double BLOCK_SIZE_FACTOR = 1.5;
1.556 + const int MIN_BLOCK_SIZE = 10;
1.557 + const double HEAD_LENGTH_FACTOR = 0.1;
1.558 + const int MIN_HEAD_LENGTH = 3;
1.559 +
1.560 + _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.561 + std::sqrt(double(_search_arc_num))),
1.562 + MIN_BLOCK_SIZE );
1.563 + _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
1.564 + MIN_HEAD_LENGTH );
1.565 + _candidates.resize(_head_length + _block_size);
1.566 + _curr_length = 0;
1.567 + }
1.568 +
1.569 + // Find next entering arc
1.570 + bool findEnteringArc() {
1.571 + // Check the current candidate list
1.572 + int e;
1.573 + for (int i = 0; i < _curr_length; ++i) {
1.574 + e = _candidates[i];
1.575 + _cand_cost[e] = _state[e] *
1.576 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.577 + if (_cand_cost[e] >= 0) {
1.578 + _candidates[i--] = _candidates[--_curr_length];
1.579 + }
1.580 + }
1.581 +
1.582 + // Extend the list
1.583 + int cnt = _block_size;
1.584 + int last_arc = 0;
1.585 + int limit = _head_length;
1.586 +
1.587 + for (int e = _next_arc; e < _search_arc_num; ++e) {
1.588 + _cand_cost[e] = _state[e] *
1.589 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.590 + if (_cand_cost[e] < 0) {
1.591 + _candidates[_curr_length++] = e;
1.592 + last_arc = e;
1.593 + }
1.594 + if (--cnt == 0) {
1.595 + if (_curr_length > limit) break;
1.596 + limit = 0;
1.597 + cnt = _block_size;
1.598 + }
1.599 + }
1.600 + if (_curr_length <= limit) {
1.601 + for (int e = 0; e < _next_arc; ++e) {
1.602 + _cand_cost[e] = _state[e] *
1.603 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.604 + if (_cand_cost[e] < 0) {
1.605 + _candidates[_curr_length++] = e;
1.606 + last_arc = e;
1.607 + }
1.608 + if (--cnt == 0) {
1.609 + if (_curr_length > limit) break;
1.610 + limit = 0;
1.611 + cnt = _block_size;
1.612 + }
1.613 + }
1.614 + }
1.615 + if (_curr_length == 0) return false;
1.616 + _next_arc = last_arc + 1;
1.617 +
1.618 + // Make heap of the candidate list (approximating a partial sort)
1.619 + make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.620 + _sort_func );
1.621 +
1.622 + // Pop the first element of the heap
1.623 + _in_arc = _candidates[0];
1.624 + pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.625 + _sort_func );
1.626 + _curr_length = std::min(_head_length, _curr_length - 1);
1.627 + return true;
1.628 + }
1.629 +
1.630 + }; //class AlteringListPivotRule
1.631 +
1.632 + public:
1.633 +
1.634 + /// \brief Constructor.
1.635 + ///
1.636 + /// The constructor of the class.
1.637 + ///
1.638 + /// \param graph The digraph the algorithm runs on.
1.639 + NetworkSimplex(const GR& graph) :
1.640 + _graph(graph), _node_id(graph), _arc_id(graph),
1.641 + INF(std::numeric_limits<Value>::has_infinity ?
1.642 + std::numeric_limits<Value>::infinity() :
1.643 + std::numeric_limits<Value>::max())
1.644 + {
1.645 + // Check the value types
1.646 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.647 + "The flow type of NetworkSimplex must be signed");
1.648 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.649 + "The cost type of NetworkSimplex must be signed");
1.650 +
1.651 + // Resize vectors
1.652 + _node_num = countNodes(_graph);
1.653 + _arc_num = countArcs(_graph);
1.654 + int all_node_num = _node_num + 1;
1.655 + int max_arc_num = _arc_num + 2 * _node_num;
1.656 +
1.657 + _source.resize(max_arc_num);
1.658 + _target.resize(max_arc_num);
1.659 +
1.660 + _lower.resize(_arc_num);
1.661 + _upper.resize(_arc_num);
1.662 + _cap.resize(max_arc_num);
1.663 + _cost.resize(max_arc_num);
1.664 + _supply.resize(all_node_num);
1.665 + _flow.resize(max_arc_num);
1.666 + _pi.resize(all_node_num);
1.667 +
1.668 + _parent.resize(all_node_num);
1.669 + _pred.resize(all_node_num);
1.670 + _forward.resize(all_node_num);
1.671 + _thread.resize(all_node_num);
1.672 + _rev_thread.resize(all_node_num);
1.673 + _succ_num.resize(all_node_num);
1.674 + _last_succ.resize(all_node_num);
1.675 + _state.resize(max_arc_num);
1.676 +
1.677 + // Copy the graph (store the arcs in a mixed order)
1.678 + int i = 0;
1.679 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.680 + _node_id[n] = i;
1.681 + }
1.682 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.683 + i = 0;
1.684 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.685 + _arc_id[a] = i;
1.686 + _source[i] = _node_id[_graph.source(a)];
1.687 + _target[i] = _node_id[_graph.target(a)];
1.688 + if ((i += k) >= _arc_num) i = (i % k) + 1;
1.689 + }
1.690 +
1.691 + // Initialize maps
1.692 + for (int i = 0; i != _node_num; ++i) {
1.693 + _supply[i] = 0;
1.694 + }
1.695 + for (int i = 0; i != _arc_num; ++i) {
1.696 + _lower[i] = 0;
1.697 + _upper[i] = INF;
1.698 + _cost[i] = 1;
1.699 + }
1.700 + _have_lower = false;
1.701 + _stype = GEQ;
1.702 + }
1.703 +
1.704 + /// \name Parameters
1.705 + /// The parameters of the algorithm can be specified using these
1.706 + /// functions.
1.707 +
1.708 + /// @{
1.709 +
1.710 + /// \brief Set the lower bounds on the arcs.
1.711 + ///
1.712 + /// This function sets the lower bounds on the arcs.
1.713 + /// If it is not used before calling \ref run(), the lower bounds
1.714 + /// will be set to zero on all arcs.
1.715 + ///
1.716 + /// \param map An arc map storing the lower bounds.
1.717 + /// Its \c Value type must be convertible to the \c Value type
1.718 + /// of the algorithm.
1.719 + ///
1.720 + /// \return <tt>(*this)</tt>
1.721 + template <typename LowerMap>
1.722 + NetworkSimplex& lowerMap(const LowerMap& map) {
1.723 + _have_lower = true;
1.724 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.725 + _lower[_arc_id[a]] = map[a];
1.726 + }
1.727 + return *this;
1.728 + }
1.729 +
1.730 + /// \brief Set the upper bounds (capacities) on the arcs.
1.731 + ///
1.732 + /// This function sets the upper bounds (capacities) on the arcs.
1.733 + /// If it is not used before calling \ref run(), the upper bounds
1.734 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.735 + /// unbounded from above on each arc).
1.736 + ///
1.737 + /// \param map An arc map storing the upper bounds.
1.738 + /// Its \c Value type must be convertible to the \c Value type
1.739 + /// of the algorithm.
1.740 + ///
1.741 + /// \return <tt>(*this)</tt>
1.742 + template<typename UpperMap>
1.743 + NetworkSimplex& upperMap(const UpperMap& map) {
1.744 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.745 + _upper[_arc_id[a]] = map[a];
1.746 + }
1.747 + return *this;
1.748 + }
1.749 +
1.750 + /// \brief Set the costs of the arcs.
1.751 + ///
1.752 + /// This function sets the costs of the arcs.
1.753 + /// If it is not used before calling \ref run(), the costs
1.754 + /// will be set to \c 1 on all arcs.
1.755 + ///
1.756 + /// \param map An arc map storing the costs.
1.757 + /// Its \c Value type must be convertible to the \c Cost type
1.758 + /// of the algorithm.
1.759 + ///
1.760 + /// \return <tt>(*this)</tt>
1.761 + template<typename CostMap>
1.762 + NetworkSimplex& costMap(const CostMap& map) {
1.763 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.764 + _cost[_arc_id[a]] = map[a];
1.765 + }
1.766 + return *this;
1.767 + }
1.768 +
1.769 + /// \brief Set the supply values of the nodes.
1.770 + ///
1.771 + /// This function sets the supply values of the nodes.
1.772 + /// If neither this function nor \ref stSupply() is used before
1.773 + /// calling \ref run(), the supply of each node will be set to zero.
1.774 + /// (It makes sense only if non-zero lower bounds are given.)
1.775 + ///
1.776 + /// \param map A node map storing the supply values.
1.777 + /// Its \c Value type must be convertible to the \c Value type
1.778 + /// of the algorithm.
1.779 + ///
1.780 + /// \return <tt>(*this)</tt>
1.781 + template<typename SupplyMap>
1.782 + NetworkSimplex& supplyMap(const SupplyMap& map) {
1.783 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.784 + _supply[_node_id[n]] = map[n];
1.785 + }
1.786 + return *this;
1.787 + }
1.788 +
1.789 + /// \brief Set single source and target nodes and a supply value.
1.790 + ///
1.791 + /// This function sets a single source node and a single target node
1.792 + /// and the required flow value.
1.793 + /// If neither this function nor \ref supplyMap() is used before
1.794 + /// calling \ref run(), the supply of each node will be set to zero.
1.795 + /// (It makes sense only if non-zero lower bounds are given.)
1.796 + ///
1.797 + /// Using this function has the same effect as using \ref supplyMap()
1.798 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.799 + /// assigned to \c t and all other nodes have zero supply value.
1.800 + ///
1.801 + /// \param s The source node.
1.802 + /// \param t The target node.
1.803 + /// \param k The required amount of flow from node \c s to node \c t
1.804 + /// (i.e. the supply of \c s and the demand of \c t).
1.805 + ///
1.806 + /// \return <tt>(*this)</tt>
1.807 + NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
1.808 + for (int i = 0; i != _node_num; ++i) {
1.809 + _supply[i] = 0;
1.810 + }
1.811 + _supply[_node_id[s]] = k;
1.812 + _supply[_node_id[t]] = -k;
1.813 + return *this;
1.814 + }
1.815 +
1.816 + /// \brief Set the type of the supply constraints.
1.817 + ///
1.818 + /// This function sets the type of the supply/demand constraints.
1.819 + /// If it is not used before calling \ref run(), the \ref GEQ supply
1.820 + /// type will be used.
1.821 + ///
1.822 + /// For more information see \ref SupplyType.
1.823 + ///
1.824 + /// \return <tt>(*this)</tt>
1.825 + NetworkSimplex& supplyType(SupplyType supply_type) {
1.826 + _stype = supply_type;
1.827 + return *this;
1.828 + }
1.829 +
1.830 + /// @}
1.831 +
1.832 + /// \name Execution Control
1.833 + /// The algorithm can be executed using \ref run().
1.834 +
1.835 + /// @{
1.836 +
1.837 + /// \brief Run the algorithm.
1.838 + ///
1.839 + /// This function runs the algorithm.
1.840 + /// The paramters can be specified using functions \ref lowerMap(),
1.841 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.842 + /// \ref supplyType().
1.843 + /// For example,
1.844 + /// \code
1.845 + /// NetworkSimplex<ListDigraph> ns(graph);
1.846 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.847 + /// .supplyMap(sup).run();
1.848 + /// \endcode
1.849 + ///
1.850 + /// This function can be called more than once. All the parameters
1.851 + /// that have been given are kept for the next call, unless
1.852 + /// \ref reset() is called, thus only the modified parameters
1.853 + /// have to be set again. See \ref reset() for examples.
1.854 + /// However the underlying digraph must not be modified after this
1.855 + /// class have been constructed, since it copies and extends the graph.
1.856 + ///
1.857 + /// \param pivot_rule The pivot rule that will be used during the
1.858 + /// algorithm. For more information see \ref PivotRule.
1.859 + ///
1.860 + /// \return \c INFEASIBLE if no feasible flow exists,
1.861 + /// \n \c OPTIMAL if the problem has optimal solution
1.862 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.863 + /// optimal flow and node potentials (primal and dual solutions),
1.864 + /// \n \c UNBOUNDED if the objective function of the problem is
1.865 + /// unbounded, i.e. there is a directed cycle having negative total
1.866 + /// cost and infinite upper bound.
1.867 + ///
1.868 + /// \see ProblemType, PivotRule
1.869 + ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.870 + if (!init()) return INFEASIBLE;
1.871 + return start(pivot_rule);
1.872 + }
1.873 +
1.874 + /// \brief Reset all the parameters that have been given before.
1.875 + ///
1.876 + /// This function resets all the paramaters that have been given
1.877 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.878 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
1.879 + ///
1.880 + /// It is useful for multiple run() calls. If this function is not
1.881 + /// used, all the parameters given before are kept for the next
1.882 + /// \ref run() call.
1.883 + /// However the underlying digraph must not be modified after this
1.884 + /// class have been constructed, since it copies and extends the graph.
1.885 + ///
1.886 + /// For example,
1.887 + /// \code
1.888 + /// NetworkSimplex<ListDigraph> ns(graph);
1.889 + ///
1.890 + /// // First run
1.891 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.892 + /// .supplyMap(sup).run();
1.893 + ///
1.894 + /// // Run again with modified cost map (reset() is not called,
1.895 + /// // so only the cost map have to be set again)
1.896 + /// cost[e] += 100;
1.897 + /// ns.costMap(cost).run();
1.898 + ///
1.899 + /// // Run again from scratch using reset()
1.900 + /// // (the lower bounds will be set to zero on all arcs)
1.901 + /// ns.reset();
1.902 + /// ns.upperMap(capacity).costMap(cost)
1.903 + /// .supplyMap(sup).run();
1.904 + /// \endcode
1.905 + ///
1.906 + /// \return <tt>(*this)</tt>
1.907 + NetworkSimplex& reset() {
1.908 + for (int i = 0; i != _node_num; ++i) {
1.909 + _supply[i] = 0;
1.910 + }
1.911 + for (int i = 0; i != _arc_num; ++i) {
1.912 + _lower[i] = 0;
1.913 + _upper[i] = INF;
1.914 + _cost[i] = 1;
1.915 + }
1.916 + _have_lower = false;
1.917 + _stype = GEQ;
1.918 + return *this;
1.919 + }
1.920 +
1.921 + /// @}
1.922 +
1.923 + /// \name Query Functions
1.924 + /// The results of the algorithm can be obtained using these
1.925 + /// functions.\n
1.926 + /// The \ref run() function must be called before using them.
1.927 +
1.928 + /// @{
1.929 +
1.930 + /// \brief Return the total cost of the found flow.
1.931 + ///
1.932 + /// This function returns the total cost of the found flow.
1.933 + /// Its complexity is O(e).
1.934 + ///
1.935 + /// \note The return type of the function can be specified as a
1.936 + /// template parameter. For example,
1.937 + /// \code
1.938 + /// ns.totalCost<double>();
1.939 + /// \endcode
1.940 + /// It is useful if the total cost cannot be stored in the \c Cost
1.941 + /// type of the algorithm, which is the default return type of the
1.942 + /// function.
1.943 + ///
1.944 + /// \pre \ref run() must be called before using this function.
1.945 + template <typename Number>
1.946 + Number totalCost() const {
1.947 + Number c = 0;
1.948 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.949 + int i = _arc_id[a];
1.950 + c += Number(_flow[i]) * Number(_cost[i]);
1.951 + }
1.952 + return c;
1.953 + }
1.954 +
1.955 +#ifndef DOXYGEN
1.956 + Cost totalCost() const {
1.957 + return totalCost<Cost>();
1.958 + }
1.959 +#endif
1.960 +
1.961 + /// \brief Return the flow on the given arc.
1.962 + ///
1.963 + /// This function returns the flow on the given arc.
1.964 + ///
1.965 + /// \pre \ref run() must be called before using this function.
1.966 + Value flow(const Arc& a) const {
1.967 + return _flow[_arc_id[a]];
1.968 + }
1.969 +
1.970 + /// \brief Return the flow map (the primal solution).
1.971 + ///
1.972 + /// This function copies the flow value on each arc into the given
1.973 + /// map. The \c Value type of the algorithm must be convertible to
1.974 + /// the \c Value type of the map.
1.975 + ///
1.976 + /// \pre \ref run() must be called before using this function.
1.977 + template <typename FlowMap>
1.978 + void flowMap(FlowMap &map) const {
1.979 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.980 + map.set(a, _flow[_arc_id[a]]);
1.981 + }
1.982 + }
1.983 +
1.984 + /// \brief Return the potential (dual value) of the given node.
1.985 + ///
1.986 + /// This function returns the potential (dual value) of the
1.987 + /// given node.
1.988 + ///
1.989 + /// \pre \ref run() must be called before using this function.
1.990 + Cost potential(const Node& n) const {
1.991 + return _pi[_node_id[n]];
1.992 + }
1.993 +
1.994 + /// \brief Return the potential map (the dual solution).
1.995 + ///
1.996 + /// This function copies the potential (dual value) of each node
1.997 + /// into the given map.
1.998 + /// The \c Cost type of the algorithm must be convertible to the
1.999 + /// \c Value type of the map.
1.1000 + ///
1.1001 + /// \pre \ref run() must be called before using this function.
1.1002 + template <typename PotentialMap>
1.1003 + void potentialMap(PotentialMap &map) const {
1.1004 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1005 + map.set(n, _pi[_node_id[n]]);
1.1006 + }
1.1007 + }
1.1008 +
1.1009 + /// @}
1.1010 +
1.1011 + private:
1.1012 +
1.1013 + // Initialize internal data structures
1.1014 + bool init() {
1.1015 + if (_node_num == 0) return false;
1.1016 +
1.1017 + // Check the sum of supply values
1.1018 + _sum_supply = 0;
1.1019 + for (int i = 0; i != _node_num; ++i) {
1.1020 + _sum_supply += _supply[i];
1.1021 + }
1.1022 + if ( !((_stype == GEQ && _sum_supply <= 0) ||
1.1023 + (_stype == LEQ && _sum_supply >= 0)) ) return false;
1.1024 +
1.1025 + // Remove non-zero lower bounds
1.1026 + if (_have_lower) {
1.1027 + for (int i = 0; i != _arc_num; ++i) {
1.1028 + Value c = _lower[i];
1.1029 + if (c >= 0) {
1.1030 + _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
1.1031 + } else {
1.1032 + _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
1.1033 + }
1.1034 + _supply[_source[i]] -= c;
1.1035 + _supply[_target[i]] += c;
1.1036 + }
1.1037 + } else {
1.1038 + for (int i = 0; i != _arc_num; ++i) {
1.1039 + _cap[i] = _upper[i];
1.1040 + }
1.1041 + }
1.1042 +
1.1043 + // Initialize artifical cost
1.1044 + Cost ART_COST;
1.1045 + if (std::numeric_limits<Cost>::is_exact) {
1.1046 + ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
1.1047 + } else {
1.1048 + ART_COST = std::numeric_limits<Cost>::min();
1.1049 + for (int i = 0; i != _arc_num; ++i) {
1.1050 + if (_cost[i] > ART_COST) ART_COST = _cost[i];
1.1051 + }
1.1052 + ART_COST = (ART_COST + 1) * _node_num;
1.1053 + }
1.1054 +
1.1055 + // Initialize arc maps
1.1056 + for (int i = 0; i != _arc_num; ++i) {
1.1057 + _flow[i] = 0;
1.1058 + _state[i] = STATE_LOWER;
1.1059 + }
1.1060 +
1.1061 + // Set data for the artificial root node
1.1062 + _root = _node_num;
1.1063 + _parent[_root] = -1;
1.1064 + _pred[_root] = -1;
1.1065 + _thread[_root] = 0;
1.1066 + _rev_thread[0] = _root;
1.1067 + _succ_num[_root] = _node_num + 1;
1.1068 + _last_succ[_root] = _root - 1;
1.1069 + _supply[_root] = -_sum_supply;
1.1070 + _pi[_root] = 0;
1.1071 +
1.1072 + // Add artificial arcs and initialize the spanning tree data structure
1.1073 + if (_sum_supply == 0) {
1.1074 + // EQ supply constraints
1.1075 + _search_arc_num = _arc_num;
1.1076 + _all_arc_num = _arc_num + _node_num;
1.1077 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1078 + _parent[u] = _root;
1.1079 + _pred[u] = e;
1.1080 + _thread[u] = u + 1;
1.1081 + _rev_thread[u + 1] = u;
1.1082 + _succ_num[u] = 1;
1.1083 + _last_succ[u] = u;
1.1084 + _cap[e] = INF;
1.1085 + _state[e] = STATE_TREE;
1.1086 + if (_supply[u] >= 0) {
1.1087 + _forward[u] = true;
1.1088 + _pi[u] = 0;
1.1089 + _source[e] = u;
1.1090 + _target[e] = _root;
1.1091 + _flow[e] = _supply[u];
1.1092 + _cost[e] = 0;
1.1093 + } else {
1.1094 + _forward[u] = false;
1.1095 + _pi[u] = ART_COST;
1.1096 + _source[e] = _root;
1.1097 + _target[e] = u;
1.1098 + _flow[e] = -_supply[u];
1.1099 + _cost[e] = ART_COST;
1.1100 + }
1.1101 + }
1.1102 + }
1.1103 + else if (_sum_supply > 0) {
1.1104 + // LEQ supply constraints
1.1105 + _search_arc_num = _arc_num + _node_num;
1.1106 + int f = _arc_num + _node_num;
1.1107 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1108 + _parent[u] = _root;
1.1109 + _thread[u] = u + 1;
1.1110 + _rev_thread[u + 1] = u;
1.1111 + _succ_num[u] = 1;
1.1112 + _last_succ[u] = u;
1.1113 + if (_supply[u] >= 0) {
1.1114 + _forward[u] = true;
1.1115 + _pi[u] = 0;
1.1116 + _pred[u] = e;
1.1117 + _source[e] = u;
1.1118 + _target[e] = _root;
1.1119 + _cap[e] = INF;
1.1120 + _flow[e] = _supply[u];
1.1121 + _cost[e] = 0;
1.1122 + _state[e] = STATE_TREE;
1.1123 + } else {
1.1124 + _forward[u] = false;
1.1125 + _pi[u] = ART_COST;
1.1126 + _pred[u] = f;
1.1127 + _source[f] = _root;
1.1128 + _target[f] = u;
1.1129 + _cap[f] = INF;
1.1130 + _flow[f] = -_supply[u];
1.1131 + _cost[f] = ART_COST;
1.1132 + _state[f] = STATE_TREE;
1.1133 + _source[e] = u;
1.1134 + _target[e] = _root;
1.1135 + _cap[e] = INF;
1.1136 + _flow[e] = 0;
1.1137 + _cost[e] = 0;
1.1138 + _state[e] = STATE_LOWER;
1.1139 + ++f;
1.1140 + }
1.1141 + }
1.1142 + _all_arc_num = f;
1.1143 + }
1.1144 + else {
1.1145 + // GEQ supply constraints
1.1146 + _search_arc_num = _arc_num + _node_num;
1.1147 + int f = _arc_num + _node_num;
1.1148 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1149 + _parent[u] = _root;
1.1150 + _thread[u] = u + 1;
1.1151 + _rev_thread[u + 1] = u;
1.1152 + _succ_num[u] = 1;
1.1153 + _last_succ[u] = u;
1.1154 + if (_supply[u] <= 0) {
1.1155 + _forward[u] = false;
1.1156 + _pi[u] = 0;
1.1157 + _pred[u] = e;
1.1158 + _source[e] = _root;
1.1159 + _target[e] = u;
1.1160 + _cap[e] = INF;
1.1161 + _flow[e] = -_supply[u];
1.1162 + _cost[e] = 0;
1.1163 + _state[e] = STATE_TREE;
1.1164 + } else {
1.1165 + _forward[u] = true;
1.1166 + _pi[u] = -ART_COST;
1.1167 + _pred[u] = f;
1.1168 + _source[f] = u;
1.1169 + _target[f] = _root;
1.1170 + _cap[f] = INF;
1.1171 + _flow[f] = _supply[u];
1.1172 + _state[f] = STATE_TREE;
1.1173 + _cost[f] = ART_COST;
1.1174 + _source[e] = _root;
1.1175 + _target[e] = u;
1.1176 + _cap[e] = INF;
1.1177 + _flow[e] = 0;
1.1178 + _cost[e] = 0;
1.1179 + _state[e] = STATE_LOWER;
1.1180 + ++f;
1.1181 + }
1.1182 + }
1.1183 + _all_arc_num = f;
1.1184 + }
1.1185 +
1.1186 + return true;
1.1187 + }
1.1188 +
1.1189 + // Find the join node
1.1190 + void findJoinNode() {
1.1191 + int u = _source[in_arc];
1.1192 + int v = _target[in_arc];
1.1193 + while (u != v) {
1.1194 + if (_succ_num[u] < _succ_num[v]) {
1.1195 + u = _parent[u];
1.1196 + } else {
1.1197 + v = _parent[v];
1.1198 + }
1.1199 + }
1.1200 + join = u;
1.1201 + }
1.1202 +
1.1203 + // Find the leaving arc of the cycle and returns true if the
1.1204 + // leaving arc is not the same as the entering arc
1.1205 + bool findLeavingArc() {
1.1206 + // Initialize first and second nodes according to the direction
1.1207 + // of the cycle
1.1208 + if (_state[in_arc] == STATE_LOWER) {
1.1209 + first = _source[in_arc];
1.1210 + second = _target[in_arc];
1.1211 + } else {
1.1212 + first = _target[in_arc];
1.1213 + second = _source[in_arc];
1.1214 + }
1.1215 + delta = _cap[in_arc];
1.1216 + int result = 0;
1.1217 + Value d;
1.1218 + int e;
1.1219 +
1.1220 + // Search the cycle along the path form the first node to the root
1.1221 + for (int u = first; u != join; u = _parent[u]) {
1.1222 + e = _pred[u];
1.1223 + d = _forward[u] ?
1.1224 + _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1.1225 + if (d < delta) {
1.1226 + delta = d;
1.1227 + u_out = u;
1.1228 + result = 1;
1.1229 + }
1.1230 + }
1.1231 + // Search the cycle along the path form the second node to the root
1.1232 + for (int u = second; u != join; u = _parent[u]) {
1.1233 + e = _pred[u];
1.1234 + d = _forward[u] ?
1.1235 + (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1.1236 + if (d <= delta) {
1.1237 + delta = d;
1.1238 + u_out = u;
1.1239 + result = 2;
1.1240 + }
1.1241 + }
1.1242 +
1.1243 + if (result == 1) {
1.1244 + u_in = first;
1.1245 + v_in = second;
1.1246 + } else {
1.1247 + u_in = second;
1.1248 + v_in = first;
1.1249 + }
1.1250 + return result != 0;
1.1251 + }
1.1252 +
1.1253 + // Change _flow and _state vectors
1.1254 + void changeFlow(bool change) {
1.1255 + // Augment along the cycle
1.1256 + if (delta > 0) {
1.1257 + Value val = _state[in_arc] * delta;
1.1258 + _flow[in_arc] += val;
1.1259 + for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1.1260 + _flow[_pred[u]] += _forward[u] ? -val : val;
1.1261 + }
1.1262 + for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1.1263 + _flow[_pred[u]] += _forward[u] ? val : -val;
1.1264 + }
1.1265 + }
1.1266 + // Update the state of the entering and leaving arcs
1.1267 + if (change) {
1.1268 + _state[in_arc] = STATE_TREE;
1.1269 + _state[_pred[u_out]] =
1.1270 + (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1.1271 + } else {
1.1272 + _state[in_arc] = -_state[in_arc];
1.1273 + }
1.1274 + }
1.1275 +
1.1276 + // Update the tree structure
1.1277 + void updateTreeStructure() {
1.1278 + int u, w;
1.1279 + int old_rev_thread = _rev_thread[u_out];
1.1280 + int old_succ_num = _succ_num[u_out];
1.1281 + int old_last_succ = _last_succ[u_out];
1.1282 + v_out = _parent[u_out];
1.1283 +
1.1284 + u = _last_succ[u_in]; // the last successor of u_in
1.1285 + right = _thread[u]; // the node after it
1.1286 +
1.1287 + // Handle the case when old_rev_thread equals to v_in
1.1288 + // (it also means that join and v_out coincide)
1.1289 + if (old_rev_thread == v_in) {
1.1290 + last = _thread[_last_succ[u_out]];
1.1291 + } else {
1.1292 + last = _thread[v_in];
1.1293 + }
1.1294 +
1.1295 + // Update _thread and _parent along the stem nodes (i.e. the nodes
1.1296 + // between u_in and u_out, whose parent have to be changed)
1.1297 + _thread[v_in] = stem = u_in;
1.1298 + _dirty_revs.clear();
1.1299 + _dirty_revs.push_back(v_in);
1.1300 + par_stem = v_in;
1.1301 + while (stem != u_out) {
1.1302 + // Insert the next stem node into the thread list
1.1303 + new_stem = _parent[stem];
1.1304 + _thread[u] = new_stem;
1.1305 + _dirty_revs.push_back(u);
1.1306 +
1.1307 + // Remove the subtree of stem from the thread list
1.1308 + w = _rev_thread[stem];
1.1309 + _thread[w] = right;
1.1310 + _rev_thread[right] = w;
1.1311 +
1.1312 + // Change the parent node and shift stem nodes
1.1313 + _parent[stem] = par_stem;
1.1314 + par_stem = stem;
1.1315 + stem = new_stem;
1.1316 +
1.1317 + // Update u and right
1.1318 + u = _last_succ[stem] == _last_succ[par_stem] ?
1.1319 + _rev_thread[par_stem] : _last_succ[stem];
1.1320 + right = _thread[u];
1.1321 + }
1.1322 + _parent[u_out] = par_stem;
1.1323 + _thread[u] = last;
1.1324 + _rev_thread[last] = u;
1.1325 + _last_succ[u_out] = u;
1.1326 +
1.1327 + // Remove the subtree of u_out from the thread list except for
1.1328 + // the case when old_rev_thread equals to v_in
1.1329 + // (it also means that join and v_out coincide)
1.1330 + if (old_rev_thread != v_in) {
1.1331 + _thread[old_rev_thread] = right;
1.1332 + _rev_thread[right] = old_rev_thread;
1.1333 + }
1.1334 +
1.1335 + // Update _rev_thread using the new _thread values
1.1336 + for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1.1337 + u = _dirty_revs[i];
1.1338 + _rev_thread[_thread[u]] = u;
1.1339 + }
1.1340 +
1.1341 + // Update _pred, _forward, _last_succ and _succ_num for the
1.1342 + // stem nodes from u_out to u_in
1.1343 + int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1.1344 + u = u_out;
1.1345 + while (u != u_in) {
1.1346 + w = _parent[u];
1.1347 + _pred[u] = _pred[w];
1.1348 + _forward[u] = !_forward[w];
1.1349 + tmp_sc += _succ_num[u] - _succ_num[w];
1.1350 + _succ_num[u] = tmp_sc;
1.1351 + _last_succ[w] = tmp_ls;
1.1352 + u = w;
1.1353 + }
1.1354 + _pred[u_in] = in_arc;
1.1355 + _forward[u_in] = (u_in == _source[in_arc]);
1.1356 + _succ_num[u_in] = old_succ_num;
1.1357 +
1.1358 + // Set limits for updating _last_succ form v_in and v_out
1.1359 + // towards the root
1.1360 + int up_limit_in = -1;
1.1361 + int up_limit_out = -1;
1.1362 + if (_last_succ[join] == v_in) {
1.1363 + up_limit_out = join;
1.1364 + } else {
1.1365 + up_limit_in = join;
1.1366 + }
1.1367 +
1.1368 + // Update _last_succ from v_in towards the root
1.1369 + for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1.1370 + u = _parent[u]) {
1.1371 + _last_succ[u] = _last_succ[u_out];
1.1372 + }
1.1373 + // Update _last_succ from v_out towards the root
1.1374 + if (join != old_rev_thread && v_in != old_rev_thread) {
1.1375 + for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.1376 + u = _parent[u]) {
1.1377 + _last_succ[u] = old_rev_thread;
1.1378 + }
1.1379 + } else {
1.1380 + for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.1381 + u = _parent[u]) {
1.1382 + _last_succ[u] = _last_succ[u_out];
1.1383 + }
1.1384 + }
1.1385 +
1.1386 + // Update _succ_num from v_in to join
1.1387 + for (u = v_in; u != join; u = _parent[u]) {
1.1388 + _succ_num[u] += old_succ_num;
1.1389 + }
1.1390 + // Update _succ_num from v_out to join
1.1391 + for (u = v_out; u != join; u = _parent[u]) {
1.1392 + _succ_num[u] -= old_succ_num;
1.1393 + }
1.1394 + }
1.1395 +
1.1396 + // Update potentials
1.1397 + void updatePotential() {
1.1398 + Cost sigma = _forward[u_in] ?
1.1399 + _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1.1400 + _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1.1401 + // Update potentials in the subtree, which has been moved
1.1402 + int end = _thread[_last_succ[u_in]];
1.1403 + for (int u = u_in; u != end; u = _thread[u]) {
1.1404 + _pi[u] += sigma;
1.1405 + }
1.1406 + }
1.1407 +
1.1408 + // Execute the algorithm
1.1409 + ProblemType start(PivotRule pivot_rule) {
1.1410 + // Select the pivot rule implementation
1.1411 + switch (pivot_rule) {
1.1412 + case FIRST_ELIGIBLE:
1.1413 + return start<FirstEligiblePivotRule>();
1.1414 + case BEST_ELIGIBLE:
1.1415 + return start<BestEligiblePivotRule>();
1.1416 + case BLOCK_SEARCH:
1.1417 + return start<BlockSearchPivotRule>();
1.1418 + case CANDIDATE_LIST:
1.1419 + return start<CandidateListPivotRule>();
1.1420 + case ALTERING_LIST:
1.1421 + return start<AlteringListPivotRule>();
1.1422 + }
1.1423 + return INFEASIBLE; // avoid warning
1.1424 + }
1.1425 +
1.1426 + template <typename PivotRuleImpl>
1.1427 + ProblemType start() {
1.1428 + PivotRuleImpl pivot(*this);
1.1429 +
1.1430 + // Execute the Network Simplex algorithm
1.1431 + while (pivot.findEnteringArc()) {
1.1432 + findJoinNode();
1.1433 + bool change = findLeavingArc();
1.1434 + if (delta >= INF) return UNBOUNDED;
1.1435 + changeFlow(change);
1.1436 + if (change) {
1.1437 + updateTreeStructure();
1.1438 + updatePotential();
1.1439 + }
1.1440 + }
1.1441 +
1.1442 + // Check feasibility
1.1443 + for (int e = _search_arc_num; e != _all_arc_num; ++e) {
1.1444 + if (_flow[e] != 0) return INFEASIBLE;
1.1445 + }
1.1446 +
1.1447 + // Transform the solution and the supply map to the original form
1.1448 + if (_have_lower) {
1.1449 + for (int i = 0; i != _arc_num; ++i) {
1.1450 + Value c = _lower[i];
1.1451 + if (c != 0) {
1.1452 + _flow[i] += c;
1.1453 + _supply[_source[i]] += c;
1.1454 + _supply[_target[i]] -= c;
1.1455 + }
1.1456 + }
1.1457 + }
1.1458 +
1.1459 + // Shift potentials to meet the requirements of the GEQ/LEQ type
1.1460 + // optimality conditions
1.1461 + if (_sum_supply == 0) {
1.1462 + if (_stype == GEQ) {
1.1463 + Cost max_pot = std::numeric_limits<Cost>::min();
1.1464 + for (int i = 0; i != _node_num; ++i) {
1.1465 + if (_pi[i] > max_pot) max_pot = _pi[i];
1.1466 + }
1.1467 + if (max_pot > 0) {
1.1468 + for (int i = 0; i != _node_num; ++i)
1.1469 + _pi[i] -= max_pot;
1.1470 + }
1.1471 + } else {
1.1472 + Cost min_pot = std::numeric_limits<Cost>::max();
1.1473 + for (int i = 0; i != _node_num; ++i) {
1.1474 + if (_pi[i] < min_pot) min_pot = _pi[i];
1.1475 + }
1.1476 + if (min_pot < 0) {
1.1477 + for (int i = 0; i != _node_num; ++i)
1.1478 + _pi[i] -= min_pot;
1.1479 + }
1.1480 + }
1.1481 + }
1.1482 +
1.1483 + return OPTIMAL;
1.1484 + }
1.1485 +
1.1486 + }; //class NetworkSimplex
1.1487 +
1.1488 + ///@}
1.1489 +
1.1490 +} //namespace lemon
1.1491 +
1.1492 +#endif //LEMON_NETWORK_SIMPLEX_H