1.1 --- a/lemon/Makefile.am Sun Dec 13 22:19:08 2009 +0100
1.2 +++ b/lemon/Makefile.am Thu Nov 12 23:17:34 2009 +0100
1.3 @@ -62,6 +62,7 @@
1.4 lemon/bin_heap.h \
1.5 lemon/binom_heap.h \
1.6 lemon/bucket_heap.h \
1.7 + lemon/capacity_scaling.h \
1.8 lemon/cbc.h \
1.9 lemon/circulation.h \
1.10 lemon/clp.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/capacity_scaling.h Thu Nov 12 23:17:34 2009 +0100
2.3 @@ -0,0 +1,717 @@
2.4 +/* -*- C++ -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library
2.7 + *
2.8 + * Copyright (C) 2003-2008
2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.11 + *
2.12 + * Permission to use, modify and distribute this software is granted
2.13 + * provided that this copyright notice appears in all copies. For
2.14 + * precise terms see the accompanying LICENSE file.
2.15 + *
2.16 + * This software is provided "AS IS" with no warranty of any kind,
2.17 + * express or implied, and with no claim as to its suitability for any
2.18 + * purpose.
2.19 + *
2.20 + */
2.21 +
2.22 +#ifndef LEMON_CAPACITY_SCALING_H
2.23 +#define LEMON_CAPACITY_SCALING_H
2.24 +
2.25 +/// \ingroup min_cost_flow
2.26 +///
2.27 +/// \file
2.28 +/// \brief Capacity scaling algorithm for finding a minimum cost flow.
2.29 +
2.30 +#include <vector>
2.31 +#include <lemon/bin_heap.h>
2.32 +
2.33 +namespace lemon {
2.34 +
2.35 + /// \addtogroup min_cost_flow
2.36 + /// @{
2.37 +
2.38 + /// \brief Implementation of the capacity scaling algorithm for
2.39 + /// finding a minimum cost flow.
2.40 + ///
2.41 + /// \ref CapacityScaling implements the capacity scaling version
2.42 + /// of the successive shortest path algorithm for finding a minimum
2.43 + /// cost flow.
2.44 + ///
2.45 + /// \tparam Digraph The digraph type the algorithm runs on.
2.46 + /// \tparam LowerMap The type of the lower bound map.
2.47 + /// \tparam CapacityMap The type of the capacity (upper bound) map.
2.48 + /// \tparam CostMap The type of the cost (length) map.
2.49 + /// \tparam SupplyMap The type of the supply map.
2.50 + ///
2.51 + /// \warning
2.52 + /// - Arc capacities and costs should be \e non-negative \e integers.
2.53 + /// - Supply values should be \e signed \e integers.
2.54 + /// - The value types of the maps should be convertible to each other.
2.55 + /// - \c CostMap::Value must be signed type.
2.56 + ///
2.57 + /// \author Peter Kovacs
2.58 + template < typename Digraph,
2.59 + typename LowerMap = typename Digraph::template ArcMap<int>,
2.60 + typename CapacityMap = typename Digraph::template ArcMap<int>,
2.61 + typename CostMap = typename Digraph::template ArcMap<int>,
2.62 + typename SupplyMap = typename Digraph::template NodeMap<int> >
2.63 + class CapacityScaling
2.64 + {
2.65 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
2.66 +
2.67 + typedef typename CapacityMap::Value Capacity;
2.68 + typedef typename CostMap::Value Cost;
2.69 + typedef typename SupplyMap::Value Supply;
2.70 + typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
2.71 + typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
2.72 + typedef typename Digraph::template NodeMap<Arc> PredMap;
2.73 +
2.74 + public:
2.75 +
2.76 + /// The type of the flow map.
2.77 + typedef typename Digraph::template ArcMap<Capacity> FlowMap;
2.78 + /// The type of the potential map.
2.79 + typedef typename Digraph::template NodeMap<Cost> PotentialMap;
2.80 +
2.81 + private:
2.82 +
2.83 + /// \brief Special implementation of the \ref Dijkstra algorithm
2.84 + /// for finding shortest paths in the residual network.
2.85 + ///
2.86 + /// \ref ResidualDijkstra is a special implementation of the
2.87 + /// \ref Dijkstra algorithm for finding shortest paths in the
2.88 + /// residual network of the digraph with respect to the reduced arc
2.89 + /// costs and modifying the node potentials according to the
2.90 + /// distance of the nodes.
2.91 + class ResidualDijkstra
2.92 + {
2.93 + typedef typename Digraph::template NodeMap<int> HeapCrossRef;
2.94 + typedef BinHeap<Cost, HeapCrossRef> Heap;
2.95 +
2.96 + private:
2.97 +
2.98 + // The digraph the algorithm runs on
2.99 + const Digraph &_graph;
2.100 +
2.101 + // The main maps
2.102 + const FlowMap &_flow;
2.103 + const CapacityArcMap &_res_cap;
2.104 + const CostMap &_cost;
2.105 + const SupplyNodeMap &_excess;
2.106 + PotentialMap &_potential;
2.107 +
2.108 + // The distance map
2.109 + PotentialMap _dist;
2.110 + // The pred arc map
2.111 + PredMap &_pred;
2.112 + // The processed (i.e. permanently labeled) nodes
2.113 + std::vector<Node> _proc_nodes;
2.114 +
2.115 + public:
2.116 +
2.117 + /// Constructor.
2.118 + ResidualDijkstra( const Digraph &digraph,
2.119 + const FlowMap &flow,
2.120 + const CapacityArcMap &res_cap,
2.121 + const CostMap &cost,
2.122 + const SupplyMap &excess,
2.123 + PotentialMap &potential,
2.124 + PredMap &pred ) :
2.125 + _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost),
2.126 + _excess(excess), _potential(potential), _dist(digraph),
2.127 + _pred(pred)
2.128 + {}
2.129 +
2.130 + /// Run the algorithm from the given source node.
2.131 + Node run(Node s, Capacity delta = 1) {
2.132 + HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
2.133 + Heap heap(heap_cross_ref);
2.134 + heap.push(s, 0);
2.135 + _pred[s] = INVALID;
2.136 + _proc_nodes.clear();
2.137 +
2.138 + // Processing nodes
2.139 + while (!heap.empty() && _excess[heap.top()] > -delta) {
2.140 + Node u = heap.top(), v;
2.141 + Cost d = heap.prio() + _potential[u], nd;
2.142 + _dist[u] = heap.prio();
2.143 + heap.pop();
2.144 + _proc_nodes.push_back(u);
2.145 +
2.146 + // Traversing outgoing arcs
2.147 + for (OutArcIt e(_graph, u); e != INVALID; ++e) {
2.148 + if (_res_cap[e] >= delta) {
2.149 + v = _graph.target(e);
2.150 + switch(heap.state(v)) {
2.151 + case Heap::PRE_HEAP:
2.152 + heap.push(v, d + _cost[e] - _potential[v]);
2.153 + _pred[v] = e;
2.154 + break;
2.155 + case Heap::IN_HEAP:
2.156 + nd = d + _cost[e] - _potential[v];
2.157 + if (nd < heap[v]) {
2.158 + heap.decrease(v, nd);
2.159 + _pred[v] = e;
2.160 + }
2.161 + break;
2.162 + case Heap::POST_HEAP:
2.163 + break;
2.164 + }
2.165 + }
2.166 + }
2.167 +
2.168 + // Traversing incoming arcs
2.169 + for (InArcIt e(_graph, u); e != INVALID; ++e) {
2.170 + if (_flow[e] >= delta) {
2.171 + v = _graph.source(e);
2.172 + switch(heap.state(v)) {
2.173 + case Heap::PRE_HEAP:
2.174 + heap.push(v, d - _cost[e] - _potential[v]);
2.175 + _pred[v] = e;
2.176 + break;
2.177 + case Heap::IN_HEAP:
2.178 + nd = d - _cost[e] - _potential[v];
2.179 + if (nd < heap[v]) {
2.180 + heap.decrease(v, nd);
2.181 + _pred[v] = e;
2.182 + }
2.183 + break;
2.184 + case Heap::POST_HEAP:
2.185 + break;
2.186 + }
2.187 + }
2.188 + }
2.189 + }
2.190 + if (heap.empty()) return INVALID;
2.191 +
2.192 + // Updating potentials of processed nodes
2.193 + Node t = heap.top();
2.194 + Cost t_dist = heap.prio();
2.195 + for (int i = 0; i < int(_proc_nodes.size()); ++i)
2.196 + _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
2.197 +
2.198 + return t;
2.199 + }
2.200 +
2.201 + }; //class ResidualDijkstra
2.202 +
2.203 + private:
2.204 +
2.205 + // The digraph the algorithm runs on
2.206 + const Digraph &_graph;
2.207 + // The original lower bound map
2.208 + const LowerMap *_lower;
2.209 + // The modified capacity map
2.210 + CapacityArcMap _capacity;
2.211 + // The original cost map
2.212 + const CostMap &_cost;
2.213 + // The modified supply map
2.214 + SupplyNodeMap _supply;
2.215 + bool _valid_supply;
2.216 +
2.217 + // Arc map of the current flow
2.218 + FlowMap *_flow;
2.219 + bool _local_flow;
2.220 + // Node map of the current potentials
2.221 + PotentialMap *_potential;
2.222 + bool _local_potential;
2.223 +
2.224 + // The residual capacity map
2.225 + CapacityArcMap _res_cap;
2.226 + // The excess map
2.227 + SupplyNodeMap _excess;
2.228 + // The excess nodes (i.e. nodes with positive excess)
2.229 + std::vector<Node> _excess_nodes;
2.230 + // The deficit nodes (i.e. nodes with negative excess)
2.231 + std::vector<Node> _deficit_nodes;
2.232 +
2.233 + // The delta parameter used for capacity scaling
2.234 + Capacity _delta;
2.235 + // The maximum number of phases
2.236 + int _phase_num;
2.237 +
2.238 + // The pred arc map
2.239 + PredMap _pred;
2.240 + // Implementation of the Dijkstra algorithm for finding augmenting
2.241 + // shortest paths in the residual network
2.242 + ResidualDijkstra *_dijkstra;
2.243 +
2.244 + public:
2.245 +
2.246 + /// \brief General constructor (with lower bounds).
2.247 + ///
2.248 + /// General constructor (with lower bounds).
2.249 + ///
2.250 + /// \param digraph The digraph the algorithm runs on.
2.251 + /// \param lower The lower bounds of the arcs.
2.252 + /// \param capacity The capacities (upper bounds) of the arcs.
2.253 + /// \param cost The cost (length) values of the arcs.
2.254 + /// \param supply The supply values of the nodes (signed).
2.255 + CapacityScaling( const Digraph &digraph,
2.256 + const LowerMap &lower,
2.257 + const CapacityMap &capacity,
2.258 + const CostMap &cost,
2.259 + const SupplyMap &supply ) :
2.260 + _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
2.261 + _supply(digraph), _flow(NULL), _local_flow(false),
2.262 + _potential(NULL), _local_potential(false),
2.263 + _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL)
2.264 + {
2.265 + Supply sum = 0;
2.266 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.267 + _supply[n] = supply[n];
2.268 + _excess[n] = supply[n];
2.269 + sum += supply[n];
2.270 + }
2.271 + _valid_supply = sum == 0;
2.272 + for (ArcIt a(_graph); a != INVALID; ++a) {
2.273 + _capacity[a] = capacity[a];
2.274 + _res_cap[a] = capacity[a];
2.275 + }
2.276 +
2.277 + // Remove non-zero lower bounds
2.278 + typename LowerMap::Value lcap;
2.279 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.280 + if ((lcap = lower[e]) != 0) {
2.281 + _capacity[e] -= lcap;
2.282 + _res_cap[e] -= lcap;
2.283 + _supply[_graph.source(e)] -= lcap;
2.284 + _supply[_graph.target(e)] += lcap;
2.285 + _excess[_graph.source(e)] -= lcap;
2.286 + _excess[_graph.target(e)] += lcap;
2.287 + }
2.288 + }
2.289 + }
2.290 +/*
2.291 + /// \brief General constructor (without lower bounds).
2.292 + ///
2.293 + /// General constructor (without lower bounds).
2.294 + ///
2.295 + /// \param digraph The digraph the algorithm runs on.
2.296 + /// \param capacity The capacities (upper bounds) of the arcs.
2.297 + /// \param cost The cost (length) values of the arcs.
2.298 + /// \param supply The supply values of the nodes (signed).
2.299 + CapacityScaling( const Digraph &digraph,
2.300 + const CapacityMap &capacity,
2.301 + const CostMap &cost,
2.302 + const SupplyMap &supply ) :
2.303 + _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
2.304 + _supply(supply), _flow(NULL), _local_flow(false),
2.305 + _potential(NULL), _local_potential(false),
2.306 + _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL)
2.307 + {
2.308 + // Check the sum of supply values
2.309 + Supply sum = 0;
2.310 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
2.311 + _valid_supply = sum == 0;
2.312 + }
2.313 +
2.314 + /// \brief Simple constructor (with lower bounds).
2.315 + ///
2.316 + /// Simple constructor (with lower bounds).
2.317 + ///
2.318 + /// \param digraph The digraph the algorithm runs on.
2.319 + /// \param lower The lower bounds of the arcs.
2.320 + /// \param capacity The capacities (upper bounds) of the arcs.
2.321 + /// \param cost The cost (length) values of the arcs.
2.322 + /// \param s The source node.
2.323 + /// \param t The target node.
2.324 + /// \param flow_value The required amount of flow from node \c s
2.325 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
2.326 + CapacityScaling( const Digraph &digraph,
2.327 + const LowerMap &lower,
2.328 + const CapacityMap &capacity,
2.329 + const CostMap &cost,
2.330 + Node s, Node t,
2.331 + Supply flow_value ) :
2.332 + _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
2.333 + _supply(digraph, 0), _flow(NULL), _local_flow(false),
2.334 + _potential(NULL), _local_potential(false),
2.335 + _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
2.336 + {
2.337 + // Remove non-zero lower bounds
2.338 + _supply[s] = _excess[s] = flow_value;
2.339 + _supply[t] = _excess[t] = -flow_value;
2.340 + typename LowerMap::Value lcap;
2.341 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.342 + if ((lcap = lower[e]) != 0) {
2.343 + _capacity[e] -= lcap;
2.344 + _res_cap[e] -= lcap;
2.345 + _supply[_graph.source(e)] -= lcap;
2.346 + _supply[_graph.target(e)] += lcap;
2.347 + _excess[_graph.source(e)] -= lcap;
2.348 + _excess[_graph.target(e)] += lcap;
2.349 + }
2.350 + }
2.351 + _valid_supply = true;
2.352 + }
2.353 +
2.354 + /// \brief Simple constructor (without lower bounds).
2.355 + ///
2.356 + /// Simple constructor (without lower bounds).
2.357 + ///
2.358 + /// \param digraph The digraph the algorithm runs on.
2.359 + /// \param capacity The capacities (upper bounds) of the arcs.
2.360 + /// \param cost The cost (length) values of the arcs.
2.361 + /// \param s The source node.
2.362 + /// \param t The target node.
2.363 + /// \param flow_value The required amount of flow from node \c s
2.364 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
2.365 + CapacityScaling( const Digraph &digraph,
2.366 + const CapacityMap &capacity,
2.367 + const CostMap &cost,
2.368 + Node s, Node t,
2.369 + Supply flow_value ) :
2.370 + _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
2.371 + _supply(digraph, 0), _flow(NULL), _local_flow(false),
2.372 + _potential(NULL), _local_potential(false),
2.373 + _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
2.374 + {
2.375 + _supply[s] = _excess[s] = flow_value;
2.376 + _supply[t] = _excess[t] = -flow_value;
2.377 + _valid_supply = true;
2.378 + }
2.379 +*/
2.380 + /// Destructor.
2.381 + ~CapacityScaling() {
2.382 + if (_local_flow) delete _flow;
2.383 + if (_local_potential) delete _potential;
2.384 + delete _dijkstra;
2.385 + }
2.386 +
2.387 + /// \brief Set the flow map.
2.388 + ///
2.389 + /// Set the flow map.
2.390 + ///
2.391 + /// \return \c (*this)
2.392 + CapacityScaling& flowMap(FlowMap &map) {
2.393 + if (_local_flow) {
2.394 + delete _flow;
2.395 + _local_flow = false;
2.396 + }
2.397 + _flow = ↦
2.398 + return *this;
2.399 + }
2.400 +
2.401 + /// \brief Set the potential map.
2.402 + ///
2.403 + /// Set the potential map.
2.404 + ///
2.405 + /// \return \c (*this)
2.406 + CapacityScaling& potentialMap(PotentialMap &map) {
2.407 + if (_local_potential) {
2.408 + delete _potential;
2.409 + _local_potential = false;
2.410 + }
2.411 + _potential = ↦
2.412 + return *this;
2.413 + }
2.414 +
2.415 + /// \name Execution control
2.416 +
2.417 + /// @{
2.418 +
2.419 + /// \brief Run the algorithm.
2.420 + ///
2.421 + /// This function runs the algorithm.
2.422 + ///
2.423 + /// \param scaling Enable or disable capacity scaling.
2.424 + /// If the maximum arc capacity and/or the amount of total supply
2.425 + /// is rather small, the algorithm could be slightly faster without
2.426 + /// scaling.
2.427 + ///
2.428 + /// \return \c true if a feasible flow can be found.
2.429 + bool run(bool scaling = true) {
2.430 + return init(scaling) && start();
2.431 + }
2.432 +
2.433 + /// @}
2.434 +
2.435 + /// \name Query Functions
2.436 + /// The results of the algorithm can be obtained using these
2.437 + /// functions.\n
2.438 + /// \ref lemon::CapacityScaling::run() "run()" must be called before
2.439 + /// using them.
2.440 +
2.441 + /// @{
2.442 +
2.443 + /// \brief Return a const reference to the arc map storing the
2.444 + /// found flow.
2.445 + ///
2.446 + /// Return a const reference to the arc map storing the found flow.
2.447 + ///
2.448 + /// \pre \ref run() must be called before using this function.
2.449 + const FlowMap& flowMap() const {
2.450 + return *_flow;
2.451 + }
2.452 +
2.453 + /// \brief Return a const reference to the node map storing the
2.454 + /// found potentials (the dual solution).
2.455 + ///
2.456 + /// Return a const reference to the node map storing the found
2.457 + /// potentials (the dual solution).
2.458 + ///
2.459 + /// \pre \ref run() must be called before using this function.
2.460 + const PotentialMap& potentialMap() const {
2.461 + return *_potential;
2.462 + }
2.463 +
2.464 + /// \brief Return the flow on the given arc.
2.465 + ///
2.466 + /// Return the flow on the given arc.
2.467 + ///
2.468 + /// \pre \ref run() must be called before using this function.
2.469 + Capacity flow(const Arc& arc) const {
2.470 + return (*_flow)[arc];
2.471 + }
2.472 +
2.473 + /// \brief Return the potential of the given node.
2.474 + ///
2.475 + /// Return the potential of the given node.
2.476 + ///
2.477 + /// \pre \ref run() must be called before using this function.
2.478 + Cost potential(const Node& node) const {
2.479 + return (*_potential)[node];
2.480 + }
2.481 +
2.482 + /// \brief Return the total cost of the found flow.
2.483 + ///
2.484 + /// Return the total cost of the found flow. The complexity of the
2.485 + /// function is \f$ O(e) \f$.
2.486 + ///
2.487 + /// \pre \ref run() must be called before using this function.
2.488 + Cost totalCost() const {
2.489 + Cost c = 0;
2.490 + for (ArcIt e(_graph); e != INVALID; ++e)
2.491 + c += (*_flow)[e] * _cost[e];
2.492 + return c;
2.493 + }
2.494 +
2.495 + /// @}
2.496 +
2.497 + private:
2.498 +
2.499 + /// Initialize the algorithm.
2.500 + bool init(bool scaling) {
2.501 + if (!_valid_supply) return false;
2.502 +
2.503 + // Initializing maps
2.504 + if (!_flow) {
2.505 + _flow = new FlowMap(_graph);
2.506 + _local_flow = true;
2.507 + }
2.508 + if (!_potential) {
2.509 + _potential = new PotentialMap(_graph);
2.510 + _local_potential = true;
2.511 + }
2.512 + for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
2.513 + for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
2.514 +
2.515 + _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
2.516 + _excess, *_potential, _pred );
2.517 +
2.518 + // Initializing delta value
2.519 + if (scaling) {
2.520 + // With scaling
2.521 + Supply max_sup = 0, max_dem = 0;
2.522 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.523 + if ( _supply[n] > max_sup) max_sup = _supply[n];
2.524 + if (-_supply[n] > max_dem) max_dem = -_supply[n];
2.525 + }
2.526 + Capacity max_cap = 0;
2.527 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.528 + if (_capacity[e] > max_cap) max_cap = _capacity[e];
2.529 + }
2.530 + max_sup = std::min(std::min(max_sup, max_dem), max_cap);
2.531 + _phase_num = 0;
2.532 + for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
2.533 + ++_phase_num;
2.534 + } else {
2.535 + // Without scaling
2.536 + _delta = 1;
2.537 + }
2.538 +
2.539 + return true;
2.540 + }
2.541 +
2.542 + bool start() {
2.543 + if (_delta > 1)
2.544 + return startWithScaling();
2.545 + else
2.546 + return startWithoutScaling();
2.547 + }
2.548 +
2.549 + /// Execute the capacity scaling algorithm.
2.550 + bool startWithScaling() {
2.551 + // Processing capacity scaling phases
2.552 + Node s, t;
2.553 + int phase_cnt = 0;
2.554 + int factor = 4;
2.555 + while (true) {
2.556 + // Saturating all arcs not satisfying the optimality condition
2.557 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.558 + Node u = _graph.source(e), v = _graph.target(e);
2.559 + Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
2.560 + if (c < 0 && _res_cap[e] >= _delta) {
2.561 + _excess[u] -= _res_cap[e];
2.562 + _excess[v] += _res_cap[e];
2.563 + (*_flow)[e] = _capacity[e];
2.564 + _res_cap[e] = 0;
2.565 + }
2.566 + else if (c > 0 && (*_flow)[e] >= _delta) {
2.567 + _excess[u] += (*_flow)[e];
2.568 + _excess[v] -= (*_flow)[e];
2.569 + (*_flow)[e] = 0;
2.570 + _res_cap[e] = _capacity[e];
2.571 + }
2.572 + }
2.573 +
2.574 + // Finding excess nodes and deficit nodes
2.575 + _excess_nodes.clear();
2.576 + _deficit_nodes.clear();
2.577 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.578 + if (_excess[n] >= _delta) _excess_nodes.push_back(n);
2.579 + if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
2.580 + }
2.581 + int next_node = 0, next_def_node = 0;
2.582 +
2.583 + // Finding augmenting shortest paths
2.584 + while (next_node < int(_excess_nodes.size())) {
2.585 + // Checking deficit nodes
2.586 + if (_delta > 1) {
2.587 + bool delta_deficit = false;
2.588 + for ( ; next_def_node < int(_deficit_nodes.size());
2.589 + ++next_def_node ) {
2.590 + if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
2.591 + delta_deficit = true;
2.592 + break;
2.593 + }
2.594 + }
2.595 + if (!delta_deficit) break;
2.596 + }
2.597 +
2.598 + // Running Dijkstra
2.599 + s = _excess_nodes[next_node];
2.600 + if ((t = _dijkstra->run(s, _delta)) == INVALID) {
2.601 + if (_delta > 1) {
2.602 + ++next_node;
2.603 + continue;
2.604 + }
2.605 + return false;
2.606 + }
2.607 +
2.608 + // Augmenting along a shortest path from s to t.
2.609 + Capacity d = std::min(_excess[s], -_excess[t]);
2.610 + Node u = t;
2.611 + Arc e;
2.612 + if (d > _delta) {
2.613 + while ((e = _pred[u]) != INVALID) {
2.614 + Capacity rc;
2.615 + if (u == _graph.target(e)) {
2.616 + rc = _res_cap[e];
2.617 + u = _graph.source(e);
2.618 + } else {
2.619 + rc = (*_flow)[e];
2.620 + u = _graph.target(e);
2.621 + }
2.622 + if (rc < d) d = rc;
2.623 + }
2.624 + }
2.625 + u = t;
2.626 + while ((e = _pred[u]) != INVALID) {
2.627 + if (u == _graph.target(e)) {
2.628 + (*_flow)[e] += d;
2.629 + _res_cap[e] -= d;
2.630 + u = _graph.source(e);
2.631 + } else {
2.632 + (*_flow)[e] -= d;
2.633 + _res_cap[e] += d;
2.634 + u = _graph.target(e);
2.635 + }
2.636 + }
2.637 + _excess[s] -= d;
2.638 + _excess[t] += d;
2.639 +
2.640 + if (_excess[s] < _delta) ++next_node;
2.641 + }
2.642 +
2.643 + if (_delta == 1) break;
2.644 + if (++phase_cnt > _phase_num / 4) factor = 2;
2.645 + _delta = _delta <= factor ? 1 : _delta / factor;
2.646 + }
2.647 +
2.648 + // Handling non-zero lower bounds
2.649 + if (_lower) {
2.650 + for (ArcIt e(_graph); e != INVALID; ++e)
2.651 + (*_flow)[e] += (*_lower)[e];
2.652 + }
2.653 + return true;
2.654 + }
2.655 +
2.656 + /// Execute the successive shortest path algorithm.
2.657 + bool startWithoutScaling() {
2.658 + // Finding excess nodes
2.659 + for (NodeIt n(_graph); n != INVALID; ++n)
2.660 + if (_excess[n] > 0) _excess_nodes.push_back(n);
2.661 + if (_excess_nodes.size() == 0) return true;
2.662 + int next_node = 0;
2.663 +
2.664 + // Finding shortest paths
2.665 + Node s, t;
2.666 + while ( _excess[_excess_nodes[next_node]] > 0 ||
2.667 + ++next_node < int(_excess_nodes.size()) )
2.668 + {
2.669 + // Running Dijkstra
2.670 + s = _excess_nodes[next_node];
2.671 + if ((t = _dijkstra->run(s)) == INVALID) return false;
2.672 +
2.673 + // Augmenting along a shortest path from s to t
2.674 + Capacity d = std::min(_excess[s], -_excess[t]);
2.675 + Node u = t;
2.676 + Arc e;
2.677 + if (d > 1) {
2.678 + while ((e = _pred[u]) != INVALID) {
2.679 + Capacity rc;
2.680 + if (u == _graph.target(e)) {
2.681 + rc = _res_cap[e];
2.682 + u = _graph.source(e);
2.683 + } else {
2.684 + rc = (*_flow)[e];
2.685 + u = _graph.target(e);
2.686 + }
2.687 + if (rc < d) d = rc;
2.688 + }
2.689 + }
2.690 + u = t;
2.691 + while ((e = _pred[u]) != INVALID) {
2.692 + if (u == _graph.target(e)) {
2.693 + (*_flow)[e] += d;
2.694 + _res_cap[e] -= d;
2.695 + u = _graph.source(e);
2.696 + } else {
2.697 + (*_flow)[e] -= d;
2.698 + _res_cap[e] += d;
2.699 + u = _graph.target(e);
2.700 + }
2.701 + }
2.702 + _excess[s] -= d;
2.703 + _excess[t] += d;
2.704 + }
2.705 +
2.706 + // Handling non-zero lower bounds
2.707 + if (_lower) {
2.708 + for (ArcIt e(_graph); e != INVALID; ++e)
2.709 + (*_flow)[e] += (*_lower)[e];
2.710 + }
2.711 + return true;
2.712 + }
2.713 +
2.714 + }; //class CapacityScaling
2.715 +
2.716 + ///@}
2.717 +
2.718 +} //namespace lemon
2.719 +
2.720 +#endif //LEMON_CAPACITY_SCALING_H