Port CapacityScaling from SVN -r3524 (#180)
authorPeter Kovacs <kpeter@inf.elte.hu>
Thu, 12 Nov 2009 23:17:34 +0100
changeset 871d3e32a777d0b
parent 870 4db8d5ccd26b
child 872 fa6f37d7a25b
Port CapacityScaling from SVN -r3524 (#180)
lemon/Makefile.am
lemon/capacity_scaling.h
     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 = &map;
   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 = &map;
   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