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