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