1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/capacity_scaling.h Thu Nov 12 23:17:34 2009 +0100
1.3 @@ -0,0 +1,717 @@
1.4 +/* -*- C++ -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library
1.7 + *
1.8 + * Copyright (C) 2003-2008
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_CAPACITY_SCALING_H
1.23 +#define LEMON_CAPACITY_SCALING_H
1.24 +
1.25 +/// \ingroup min_cost_flow
1.26 +///
1.27 +/// \file
1.28 +/// \brief Capacity scaling algorithm for finding a minimum cost flow.
1.29 +
1.30 +#include <vector>
1.31 +#include <lemon/bin_heap.h>
1.32 +
1.33 +namespace lemon {
1.34 +
1.35 + /// \addtogroup min_cost_flow
1.36 + /// @{
1.37 +
1.38 + /// \brief Implementation of the capacity scaling algorithm for
1.39 + /// finding a minimum cost flow.
1.40 + ///
1.41 + /// \ref CapacityScaling implements the capacity scaling version
1.42 + /// of the successive shortest path algorithm for finding a minimum
1.43 + /// cost flow.
1.44 + ///
1.45 + /// \tparam Digraph The digraph type the algorithm runs on.
1.46 + /// \tparam LowerMap The type of the lower bound map.
1.47 + /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.48 + /// \tparam CostMap The type of the cost (length) map.
1.49 + /// \tparam SupplyMap The type of the supply map.
1.50 + ///
1.51 + /// \warning
1.52 + /// - Arc capacities and costs should be \e non-negative \e integers.
1.53 + /// - Supply values should be \e signed \e integers.
1.54 + /// - The value types of the maps should be convertible to each other.
1.55 + /// - \c CostMap::Value must be signed type.
1.56 + ///
1.57 + /// \author Peter Kovacs
1.58 + template < typename Digraph,
1.59 + typename LowerMap = typename Digraph::template ArcMap<int>,
1.60 + typename CapacityMap = typename Digraph::template ArcMap<int>,
1.61 + typename CostMap = typename Digraph::template ArcMap<int>,
1.62 + typename SupplyMap = typename Digraph::template NodeMap<int> >
1.63 + class CapacityScaling
1.64 + {
1.65 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.66 +
1.67 + typedef typename CapacityMap::Value Capacity;
1.68 + typedef typename CostMap::Value Cost;
1.69 + typedef typename SupplyMap::Value Supply;
1.70 + typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.71 + typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.72 + typedef typename Digraph::template NodeMap<Arc> PredMap;
1.73 +
1.74 + public:
1.75 +
1.76 + /// The type of the flow map.
1.77 + typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.78 + /// The type of the potential map.
1.79 + typedef typename Digraph::template NodeMap<Cost> PotentialMap;
1.80 +
1.81 + private:
1.82 +
1.83 + /// \brief Special implementation of the \ref Dijkstra algorithm
1.84 + /// for finding shortest paths in the residual network.
1.85 + ///
1.86 + /// \ref ResidualDijkstra is a special implementation of the
1.87 + /// \ref Dijkstra algorithm for finding shortest paths in the
1.88 + /// residual network of the digraph with respect to the reduced arc
1.89 + /// costs and modifying the node potentials according to the
1.90 + /// distance of the nodes.
1.91 + class ResidualDijkstra
1.92 + {
1.93 + typedef typename Digraph::template NodeMap<int> HeapCrossRef;
1.94 + typedef BinHeap<Cost, HeapCrossRef> Heap;
1.95 +
1.96 + private:
1.97 +
1.98 + // The digraph the algorithm runs on
1.99 + const Digraph &_graph;
1.100 +
1.101 + // The main maps
1.102 + const FlowMap &_flow;
1.103 + const CapacityArcMap &_res_cap;
1.104 + const CostMap &_cost;
1.105 + const SupplyNodeMap &_excess;
1.106 + PotentialMap &_potential;
1.107 +
1.108 + // The distance map
1.109 + PotentialMap _dist;
1.110 + // The pred arc map
1.111 + PredMap &_pred;
1.112 + // The processed (i.e. permanently labeled) nodes
1.113 + std::vector<Node> _proc_nodes;
1.114 +
1.115 + public:
1.116 +
1.117 + /// Constructor.
1.118 + ResidualDijkstra( const Digraph &digraph,
1.119 + const FlowMap &flow,
1.120 + const CapacityArcMap &res_cap,
1.121 + const CostMap &cost,
1.122 + const SupplyMap &excess,
1.123 + PotentialMap &potential,
1.124 + PredMap &pred ) :
1.125 + _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost),
1.126 + _excess(excess), _potential(potential), _dist(digraph),
1.127 + _pred(pred)
1.128 + {}
1.129 +
1.130 + /// Run the algorithm from the given source node.
1.131 + Node run(Node s, Capacity delta = 1) {
1.132 + HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
1.133 + Heap heap(heap_cross_ref);
1.134 + heap.push(s, 0);
1.135 + _pred[s] = INVALID;
1.136 + _proc_nodes.clear();
1.137 +
1.138 + // Processing nodes
1.139 + while (!heap.empty() && _excess[heap.top()] > -delta) {
1.140 + Node u = heap.top(), v;
1.141 + Cost d = heap.prio() + _potential[u], nd;
1.142 + _dist[u] = heap.prio();
1.143 + heap.pop();
1.144 + _proc_nodes.push_back(u);
1.145 +
1.146 + // Traversing outgoing arcs
1.147 + for (OutArcIt e(_graph, u); e != INVALID; ++e) {
1.148 + if (_res_cap[e] >= delta) {
1.149 + v = _graph.target(e);
1.150 + switch(heap.state(v)) {
1.151 + case Heap::PRE_HEAP:
1.152 + heap.push(v, d + _cost[e] - _potential[v]);
1.153 + _pred[v] = e;
1.154 + break;
1.155 + case Heap::IN_HEAP:
1.156 + nd = d + _cost[e] - _potential[v];
1.157 + if (nd < heap[v]) {
1.158 + heap.decrease(v, nd);
1.159 + _pred[v] = e;
1.160 + }
1.161 + break;
1.162 + case Heap::POST_HEAP:
1.163 + break;
1.164 + }
1.165 + }
1.166 + }
1.167 +
1.168 + // Traversing incoming arcs
1.169 + for (InArcIt e(_graph, u); e != INVALID; ++e) {
1.170 + if (_flow[e] >= delta) {
1.171 + v = _graph.source(e);
1.172 + switch(heap.state(v)) {
1.173 + case Heap::PRE_HEAP:
1.174 + heap.push(v, d - _cost[e] - _potential[v]);
1.175 + _pred[v] = e;
1.176 + break;
1.177 + case Heap::IN_HEAP:
1.178 + nd = d - _cost[e] - _potential[v];
1.179 + if (nd < heap[v]) {
1.180 + heap.decrease(v, nd);
1.181 + _pred[v] = e;
1.182 + }
1.183 + break;
1.184 + case Heap::POST_HEAP:
1.185 + break;
1.186 + }
1.187 + }
1.188 + }
1.189 + }
1.190 + if (heap.empty()) return INVALID;
1.191 +
1.192 + // Updating potentials of processed nodes
1.193 + Node t = heap.top();
1.194 + Cost t_dist = heap.prio();
1.195 + for (int i = 0; i < int(_proc_nodes.size()); ++i)
1.196 + _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
1.197 +
1.198 + return t;
1.199 + }
1.200 +
1.201 + }; //class ResidualDijkstra
1.202 +
1.203 + private:
1.204 +
1.205 + // The digraph the algorithm runs on
1.206 + const Digraph &_graph;
1.207 + // The original lower bound map
1.208 + const LowerMap *_lower;
1.209 + // The modified capacity map
1.210 + CapacityArcMap _capacity;
1.211 + // The original cost map
1.212 + const CostMap &_cost;
1.213 + // The modified supply map
1.214 + SupplyNodeMap _supply;
1.215 + bool _valid_supply;
1.216 +
1.217 + // Arc map of the current flow
1.218 + FlowMap *_flow;
1.219 + bool _local_flow;
1.220 + // Node map of the current potentials
1.221 + PotentialMap *_potential;
1.222 + bool _local_potential;
1.223 +
1.224 + // The residual capacity map
1.225 + CapacityArcMap _res_cap;
1.226 + // The excess map
1.227 + SupplyNodeMap _excess;
1.228 + // The excess nodes (i.e. nodes with positive excess)
1.229 + std::vector<Node> _excess_nodes;
1.230 + // The deficit nodes (i.e. nodes with negative excess)
1.231 + std::vector<Node> _deficit_nodes;
1.232 +
1.233 + // The delta parameter used for capacity scaling
1.234 + Capacity _delta;
1.235 + // The maximum number of phases
1.236 + int _phase_num;
1.237 +
1.238 + // The pred arc map
1.239 + PredMap _pred;
1.240 + // Implementation of the Dijkstra algorithm for finding augmenting
1.241 + // shortest paths in the residual network
1.242 + ResidualDijkstra *_dijkstra;
1.243 +
1.244 + public:
1.245 +
1.246 + /// \brief General constructor (with lower bounds).
1.247 + ///
1.248 + /// General constructor (with lower bounds).
1.249 + ///
1.250 + /// \param digraph The digraph the algorithm runs on.
1.251 + /// \param lower The lower bounds of the arcs.
1.252 + /// \param capacity The capacities (upper bounds) of the arcs.
1.253 + /// \param cost The cost (length) values of the arcs.
1.254 + /// \param supply The supply values of the nodes (signed).
1.255 + CapacityScaling( const Digraph &digraph,
1.256 + const LowerMap &lower,
1.257 + const CapacityMap &capacity,
1.258 + const CostMap &cost,
1.259 + const SupplyMap &supply ) :
1.260 + _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
1.261 + _supply(digraph), _flow(NULL), _local_flow(false),
1.262 + _potential(NULL), _local_potential(false),
1.263 + _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL)
1.264 + {
1.265 + Supply sum = 0;
1.266 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.267 + _supply[n] = supply[n];
1.268 + _excess[n] = supply[n];
1.269 + sum += supply[n];
1.270 + }
1.271 + _valid_supply = sum == 0;
1.272 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.273 + _capacity[a] = capacity[a];
1.274 + _res_cap[a] = capacity[a];
1.275 + }
1.276 +
1.277 + // Remove non-zero lower bounds
1.278 + typename LowerMap::Value lcap;
1.279 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.280 + if ((lcap = lower[e]) != 0) {
1.281 + _capacity[e] -= lcap;
1.282 + _res_cap[e] -= lcap;
1.283 + _supply[_graph.source(e)] -= lcap;
1.284 + _supply[_graph.target(e)] += lcap;
1.285 + _excess[_graph.source(e)] -= lcap;
1.286 + _excess[_graph.target(e)] += lcap;
1.287 + }
1.288 + }
1.289 + }
1.290 +/*
1.291 + /// \brief General constructor (without lower bounds).
1.292 + ///
1.293 + /// General constructor (without lower bounds).
1.294 + ///
1.295 + /// \param digraph The digraph the algorithm runs on.
1.296 + /// \param capacity The capacities (upper bounds) of the arcs.
1.297 + /// \param cost The cost (length) values of the arcs.
1.298 + /// \param supply The supply values of the nodes (signed).
1.299 + CapacityScaling( const Digraph &digraph,
1.300 + const CapacityMap &capacity,
1.301 + const CostMap &cost,
1.302 + const SupplyMap &supply ) :
1.303 + _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.304 + _supply(supply), _flow(NULL), _local_flow(false),
1.305 + _potential(NULL), _local_potential(false),
1.306 + _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL)
1.307 + {
1.308 + // Check the sum of supply values
1.309 + Supply sum = 0;
1.310 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.311 + _valid_supply = sum == 0;
1.312 + }
1.313 +
1.314 + /// \brief Simple constructor (with lower bounds).
1.315 + ///
1.316 + /// Simple constructor (with lower bounds).
1.317 + ///
1.318 + /// \param digraph The digraph the algorithm runs on.
1.319 + /// \param lower The lower bounds of the arcs.
1.320 + /// \param capacity The capacities (upper bounds) of the arcs.
1.321 + /// \param cost The cost (length) values of the arcs.
1.322 + /// \param s The source node.
1.323 + /// \param t The target node.
1.324 + /// \param flow_value The required amount of flow from node \c s
1.325 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.326 + CapacityScaling( const Digraph &digraph,
1.327 + const LowerMap &lower,
1.328 + const CapacityMap &capacity,
1.329 + const CostMap &cost,
1.330 + Node s, Node t,
1.331 + Supply flow_value ) :
1.332 + _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
1.333 + _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.334 + _potential(NULL), _local_potential(false),
1.335 + _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
1.336 + {
1.337 + // Remove non-zero lower bounds
1.338 + _supply[s] = _excess[s] = flow_value;
1.339 + _supply[t] = _excess[t] = -flow_value;
1.340 + typename LowerMap::Value lcap;
1.341 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.342 + if ((lcap = lower[e]) != 0) {
1.343 + _capacity[e] -= lcap;
1.344 + _res_cap[e] -= lcap;
1.345 + _supply[_graph.source(e)] -= lcap;
1.346 + _supply[_graph.target(e)] += lcap;
1.347 + _excess[_graph.source(e)] -= lcap;
1.348 + _excess[_graph.target(e)] += lcap;
1.349 + }
1.350 + }
1.351 + _valid_supply = true;
1.352 + }
1.353 +
1.354 + /// \brief Simple constructor (without lower bounds).
1.355 + ///
1.356 + /// Simple constructor (without lower bounds).
1.357 + ///
1.358 + /// \param digraph The digraph the algorithm runs on.
1.359 + /// \param capacity The capacities (upper bounds) of the arcs.
1.360 + /// \param cost The cost (length) values of the arcs.
1.361 + /// \param s The source node.
1.362 + /// \param t The target node.
1.363 + /// \param flow_value The required amount of flow from node \c s
1.364 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.365 + CapacityScaling( const Digraph &digraph,
1.366 + const CapacityMap &capacity,
1.367 + const CostMap &cost,
1.368 + Node s, Node t,
1.369 + Supply flow_value ) :
1.370 + _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.371 + _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.372 + _potential(NULL), _local_potential(false),
1.373 + _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
1.374 + {
1.375 + _supply[s] = _excess[s] = flow_value;
1.376 + _supply[t] = _excess[t] = -flow_value;
1.377 + _valid_supply = true;
1.378 + }
1.379 +*/
1.380 + /// Destructor.
1.381 + ~CapacityScaling() {
1.382 + if (_local_flow) delete _flow;
1.383 + if (_local_potential) delete _potential;
1.384 + delete _dijkstra;
1.385 + }
1.386 +
1.387 + /// \brief Set the flow map.
1.388 + ///
1.389 + /// Set the flow map.
1.390 + ///
1.391 + /// \return \c (*this)
1.392 + CapacityScaling& flowMap(FlowMap &map) {
1.393 + if (_local_flow) {
1.394 + delete _flow;
1.395 + _local_flow = false;
1.396 + }
1.397 + _flow = ↦
1.398 + return *this;
1.399 + }
1.400 +
1.401 + /// \brief Set the potential map.
1.402 + ///
1.403 + /// Set the potential map.
1.404 + ///
1.405 + /// \return \c (*this)
1.406 + CapacityScaling& potentialMap(PotentialMap &map) {
1.407 + if (_local_potential) {
1.408 + delete _potential;
1.409 + _local_potential = false;
1.410 + }
1.411 + _potential = ↦
1.412 + return *this;
1.413 + }
1.414 +
1.415 + /// \name Execution control
1.416 +
1.417 + /// @{
1.418 +
1.419 + /// \brief Run the algorithm.
1.420 + ///
1.421 + /// This function runs the algorithm.
1.422 + ///
1.423 + /// \param scaling Enable or disable capacity scaling.
1.424 + /// If the maximum arc capacity and/or the amount of total supply
1.425 + /// is rather small, the algorithm could be slightly faster without
1.426 + /// scaling.
1.427 + ///
1.428 + /// \return \c true if a feasible flow can be found.
1.429 + bool run(bool scaling = true) {
1.430 + return init(scaling) && start();
1.431 + }
1.432 +
1.433 + /// @}
1.434 +
1.435 + /// \name Query Functions
1.436 + /// The results of the algorithm can be obtained using these
1.437 + /// functions.\n
1.438 + /// \ref lemon::CapacityScaling::run() "run()" must be called before
1.439 + /// using them.
1.440 +
1.441 + /// @{
1.442 +
1.443 + /// \brief Return a const reference to the arc map storing the
1.444 + /// found flow.
1.445 + ///
1.446 + /// Return a const reference to the arc map storing the found flow.
1.447 + ///
1.448 + /// \pre \ref run() must be called before using this function.
1.449 + const FlowMap& flowMap() const {
1.450 + return *_flow;
1.451 + }
1.452 +
1.453 + /// \brief Return a const reference to the node map storing the
1.454 + /// found potentials (the dual solution).
1.455 + ///
1.456 + /// Return a const reference to the node map storing the found
1.457 + /// potentials (the dual solution).
1.458 + ///
1.459 + /// \pre \ref run() must be called before using this function.
1.460 + const PotentialMap& potentialMap() const {
1.461 + return *_potential;
1.462 + }
1.463 +
1.464 + /// \brief Return the flow on the given arc.
1.465 + ///
1.466 + /// Return the flow on the given arc.
1.467 + ///
1.468 + /// \pre \ref run() must be called before using this function.
1.469 + Capacity flow(const Arc& arc) const {
1.470 + return (*_flow)[arc];
1.471 + }
1.472 +
1.473 + /// \brief Return the potential of the given node.
1.474 + ///
1.475 + /// Return the potential of the given node.
1.476 + ///
1.477 + /// \pre \ref run() must be called before using this function.
1.478 + Cost potential(const Node& node) const {
1.479 + return (*_potential)[node];
1.480 + }
1.481 +
1.482 + /// \brief Return the total cost of the found flow.
1.483 + ///
1.484 + /// Return the total cost of the found flow. The complexity of the
1.485 + /// function is \f$ O(e) \f$.
1.486 + ///
1.487 + /// \pre \ref run() must be called before using this function.
1.488 + Cost totalCost() const {
1.489 + Cost c = 0;
1.490 + for (ArcIt e(_graph); e != INVALID; ++e)
1.491 + c += (*_flow)[e] * _cost[e];
1.492 + return c;
1.493 + }
1.494 +
1.495 + /// @}
1.496 +
1.497 + private:
1.498 +
1.499 + /// Initialize the algorithm.
1.500 + bool init(bool scaling) {
1.501 + if (!_valid_supply) return false;
1.502 +
1.503 + // Initializing maps
1.504 + if (!_flow) {
1.505 + _flow = new FlowMap(_graph);
1.506 + _local_flow = true;
1.507 + }
1.508 + if (!_potential) {
1.509 + _potential = new PotentialMap(_graph);
1.510 + _local_potential = true;
1.511 + }
1.512 + for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
1.513 + for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
1.514 +
1.515 + _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
1.516 + _excess, *_potential, _pred );
1.517 +
1.518 + // Initializing delta value
1.519 + if (scaling) {
1.520 + // With scaling
1.521 + Supply max_sup = 0, max_dem = 0;
1.522 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.523 + if ( _supply[n] > max_sup) max_sup = _supply[n];
1.524 + if (-_supply[n] > max_dem) max_dem = -_supply[n];
1.525 + }
1.526 + Capacity max_cap = 0;
1.527 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.528 + if (_capacity[e] > max_cap) max_cap = _capacity[e];
1.529 + }
1.530 + max_sup = std::min(std::min(max_sup, max_dem), max_cap);
1.531 + _phase_num = 0;
1.532 + for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
1.533 + ++_phase_num;
1.534 + } else {
1.535 + // Without scaling
1.536 + _delta = 1;
1.537 + }
1.538 +
1.539 + return true;
1.540 + }
1.541 +
1.542 + bool start() {
1.543 + if (_delta > 1)
1.544 + return startWithScaling();
1.545 + else
1.546 + return startWithoutScaling();
1.547 + }
1.548 +
1.549 + /// Execute the capacity scaling algorithm.
1.550 + bool startWithScaling() {
1.551 + // Processing capacity scaling phases
1.552 + Node s, t;
1.553 + int phase_cnt = 0;
1.554 + int factor = 4;
1.555 + while (true) {
1.556 + // Saturating all arcs not satisfying the optimality condition
1.557 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.558 + Node u = _graph.source(e), v = _graph.target(e);
1.559 + Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
1.560 + if (c < 0 && _res_cap[e] >= _delta) {
1.561 + _excess[u] -= _res_cap[e];
1.562 + _excess[v] += _res_cap[e];
1.563 + (*_flow)[e] = _capacity[e];
1.564 + _res_cap[e] = 0;
1.565 + }
1.566 + else if (c > 0 && (*_flow)[e] >= _delta) {
1.567 + _excess[u] += (*_flow)[e];
1.568 + _excess[v] -= (*_flow)[e];
1.569 + (*_flow)[e] = 0;
1.570 + _res_cap[e] = _capacity[e];
1.571 + }
1.572 + }
1.573 +
1.574 + // Finding excess nodes and deficit nodes
1.575 + _excess_nodes.clear();
1.576 + _deficit_nodes.clear();
1.577 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.578 + if (_excess[n] >= _delta) _excess_nodes.push_back(n);
1.579 + if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
1.580 + }
1.581 + int next_node = 0, next_def_node = 0;
1.582 +
1.583 + // Finding augmenting shortest paths
1.584 + while (next_node < int(_excess_nodes.size())) {
1.585 + // Checking deficit nodes
1.586 + if (_delta > 1) {
1.587 + bool delta_deficit = false;
1.588 + for ( ; next_def_node < int(_deficit_nodes.size());
1.589 + ++next_def_node ) {
1.590 + if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
1.591 + delta_deficit = true;
1.592 + break;
1.593 + }
1.594 + }
1.595 + if (!delta_deficit) break;
1.596 + }
1.597 +
1.598 + // Running Dijkstra
1.599 + s = _excess_nodes[next_node];
1.600 + if ((t = _dijkstra->run(s, _delta)) == INVALID) {
1.601 + if (_delta > 1) {
1.602 + ++next_node;
1.603 + continue;
1.604 + }
1.605 + return false;
1.606 + }
1.607 +
1.608 + // Augmenting along a shortest path from s to t.
1.609 + Capacity d = std::min(_excess[s], -_excess[t]);
1.610 + Node u = t;
1.611 + Arc e;
1.612 + if (d > _delta) {
1.613 + while ((e = _pred[u]) != INVALID) {
1.614 + Capacity rc;
1.615 + if (u == _graph.target(e)) {
1.616 + rc = _res_cap[e];
1.617 + u = _graph.source(e);
1.618 + } else {
1.619 + rc = (*_flow)[e];
1.620 + u = _graph.target(e);
1.621 + }
1.622 + if (rc < d) d = rc;
1.623 + }
1.624 + }
1.625 + u = t;
1.626 + while ((e = _pred[u]) != INVALID) {
1.627 + if (u == _graph.target(e)) {
1.628 + (*_flow)[e] += d;
1.629 + _res_cap[e] -= d;
1.630 + u = _graph.source(e);
1.631 + } else {
1.632 + (*_flow)[e] -= d;
1.633 + _res_cap[e] += d;
1.634 + u = _graph.target(e);
1.635 + }
1.636 + }
1.637 + _excess[s] -= d;
1.638 + _excess[t] += d;
1.639 +
1.640 + if (_excess[s] < _delta) ++next_node;
1.641 + }
1.642 +
1.643 + if (_delta == 1) break;
1.644 + if (++phase_cnt > _phase_num / 4) factor = 2;
1.645 + _delta = _delta <= factor ? 1 : _delta / factor;
1.646 + }
1.647 +
1.648 + // Handling non-zero lower bounds
1.649 + if (_lower) {
1.650 + for (ArcIt e(_graph); e != INVALID; ++e)
1.651 + (*_flow)[e] += (*_lower)[e];
1.652 + }
1.653 + return true;
1.654 + }
1.655 +
1.656 + /// Execute the successive shortest path algorithm.
1.657 + bool startWithoutScaling() {
1.658 + // Finding excess nodes
1.659 + for (NodeIt n(_graph); n != INVALID; ++n)
1.660 + if (_excess[n] > 0) _excess_nodes.push_back(n);
1.661 + if (_excess_nodes.size() == 0) return true;
1.662 + int next_node = 0;
1.663 +
1.664 + // Finding shortest paths
1.665 + Node s, t;
1.666 + while ( _excess[_excess_nodes[next_node]] > 0 ||
1.667 + ++next_node < int(_excess_nodes.size()) )
1.668 + {
1.669 + // Running Dijkstra
1.670 + s = _excess_nodes[next_node];
1.671 + if ((t = _dijkstra->run(s)) == INVALID) return false;
1.672 +
1.673 + // Augmenting along a shortest path from s to t
1.674 + Capacity d = std::min(_excess[s], -_excess[t]);
1.675 + Node u = t;
1.676 + Arc e;
1.677 + if (d > 1) {
1.678 + while ((e = _pred[u]) != INVALID) {
1.679 + Capacity rc;
1.680 + if (u == _graph.target(e)) {
1.681 + rc = _res_cap[e];
1.682 + u = _graph.source(e);
1.683 + } else {
1.684 + rc = (*_flow)[e];
1.685 + u = _graph.target(e);
1.686 + }
1.687 + if (rc < d) d = rc;
1.688 + }
1.689 + }
1.690 + u = t;
1.691 + while ((e = _pred[u]) != INVALID) {
1.692 + if (u == _graph.target(e)) {
1.693 + (*_flow)[e] += d;
1.694 + _res_cap[e] -= d;
1.695 + u = _graph.source(e);
1.696 + } else {
1.697 + (*_flow)[e] -= d;
1.698 + _res_cap[e] += d;
1.699 + u = _graph.target(e);
1.700 + }
1.701 + }
1.702 + _excess[s] -= d;
1.703 + _excess[t] += d;
1.704 + }
1.705 +
1.706 + // Handling non-zero lower bounds
1.707 + if (_lower) {
1.708 + for (ArcIt e(_graph); e != INVALID; ++e)
1.709 + (*_flow)[e] += (*_lower)[e];
1.710 + }
1.711 + return true;
1.712 + }
1.713 +
1.714 + }; //class CapacityScaling
1.715 +
1.716 + ///@}
1.717 +
1.718 +} //namespace lemon
1.719 +
1.720 +#endif //LEMON_CAPACITY_SCALING_H