COIN-OR::LEMON - Graph Library

Ticket #180: 180-cas1-317a0e41f3d5.patch

File 180-cas1-317a0e41f3d5.patch, 23.7 KB (added by Peter Kovacs, 10 years ago)
  • lemon/Makefile.am

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1252577519 -7200
    # Node ID 317a0e41f3d5defc8f9aea733c56a77cd22ba1ca
    # Parent  6d5f547e5bfb8131568ad50ea406129a056ac9ed
    Port CapacityScaling from SVN -r3524 (#180)
    
    diff --git a/lemon/Makefile.am b/lemon/Makefile.am
    a b  
    6262        lemon/bin_heap.h \
    6363        lemon/binom_heap.h \
    6464        lemon/bucket_heap.h \
     65        lemon/capacity_scaling.h \
    6566        lemon/cbc.h \
    6667        lemon/circulation.h \
    6768        lemon/clp.h \
  • new file lemon/capacity_scaling.h

    diff --git a/lemon/capacity_scaling.h b/lemon/capacity_scaling.h
    new file mode 100644
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_CAPACITY_SCALING_H
     20#define LEMON_CAPACITY_SCALING_H
     21
     22/// \ingroup min_cost_flow
     23///
     24/// \file
     25/// \brief Capacity scaling algorithm for finding a minimum cost flow.
     26
     27#include <vector>
     28#include <lemon/bin_heap.h>
     29
     30namespace lemon {
     31
     32  /// \addtogroup min_cost_flow
     33  /// @{
     34
     35  /// \brief Implementation of the capacity scaling algorithm for
     36  /// finding a minimum cost flow.
     37  ///
     38  /// \ref CapacityScaling implements the capacity scaling version
     39  /// of the successive shortest path algorithm for finding a minimum
     40  /// cost flow.
     41  ///
     42  /// \tparam Digraph The digraph type the algorithm runs on.
     43  /// \tparam LowerMap The type of the lower bound map.
     44  /// \tparam CapacityMap The type of the capacity (upper bound) map.
     45  /// \tparam CostMap The type of the cost (length) map.
     46  /// \tparam SupplyMap The type of the supply map.
     47  ///
     48  /// \warning
     49  /// - Arc capacities and costs should be \e non-negative \e integers.
     50  /// - Supply values should be \e signed \e integers.
     51  /// - The value types of the maps should be convertible to each other.
     52  /// - \c CostMap::Value must be signed type.
     53  ///
     54  /// \author Peter Kovacs
     55  template < typename Digraph,
     56             typename LowerMap = typename Digraph::template ArcMap<int>,
     57             typename CapacityMap = typename Digraph::template ArcMap<int>,
     58             typename CostMap = typename Digraph::template ArcMap<int>,
     59             typename SupplyMap = typename Digraph::template NodeMap<int> >
     60  class CapacityScaling
     61  {
     62    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
     63
     64    typedef typename CapacityMap::Value Capacity;
     65    typedef typename CostMap::Value Cost;
     66    typedef typename SupplyMap::Value Supply;
     67    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
     68    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
     69    typedef typename Digraph::template NodeMap<Arc> PredMap;
     70
     71  public:
     72
     73    /// The type of the flow map.
     74    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
     75    /// The type of the potential map.
     76    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
     77
     78  private:
     79
     80    /// \brief Special implementation of the \ref Dijkstra algorithm
     81    /// for finding shortest paths in the residual network.
     82    ///
     83    /// \ref ResidualDijkstra is a special implementation of the
     84    /// \ref Dijkstra algorithm for finding shortest paths in the
     85    /// residual network of the digraph with respect to the reduced arc
     86    /// costs and modifying the node potentials according to the
     87    /// distance of the nodes.
     88    class ResidualDijkstra
     89    {
     90      typedef typename Digraph::template NodeMap<int> HeapCrossRef;
     91      typedef BinHeap<Cost, HeapCrossRef> Heap;
     92
     93    private:
     94
     95      // The digraph the algorithm runs on
     96      const Digraph &_graph;
     97
     98      // The main maps
     99      const FlowMap &_flow;
     100      const CapacityArcMap &_res_cap;
     101      const CostMap &_cost;
     102      const SupplyNodeMap &_excess;
     103      PotentialMap &_potential;
     104
     105      // The distance map
     106      PotentialMap _dist;
     107      // The pred arc map
     108      PredMap &_pred;
     109      // The processed (i.e. permanently labeled) nodes
     110      std::vector<Node> _proc_nodes;
     111
     112    public:
     113
     114      /// Constructor.
     115      ResidualDijkstra( const Digraph &digraph,
     116                        const FlowMap &flow,
     117                        const CapacityArcMap &res_cap,
     118                        const CostMap &cost,
     119                        const SupplyMap &excess,
     120                        PotentialMap &potential,
     121                        PredMap &pred ) :
     122        _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost),
     123        _excess(excess), _potential(potential), _dist(digraph),
     124        _pred(pred)
     125      {}
     126
     127      /// Run the algorithm from the given source node.
     128      Node run(Node s, Capacity delta = 1) {
     129        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
     130        Heap heap(heap_cross_ref);
     131        heap.push(s, 0);
     132        _pred[s] = INVALID;
     133        _proc_nodes.clear();
     134
     135        // Processing nodes
     136        while (!heap.empty() && _excess[heap.top()] > -delta) {
     137          Node u = heap.top(), v;
     138          Cost d = heap.prio() + _potential[u], nd;
     139          _dist[u] = heap.prio();
     140          heap.pop();
     141          _proc_nodes.push_back(u);
     142
     143          // Traversing outgoing arcs
     144          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
     145            if (_res_cap[e] >= delta) {
     146              v = _graph.target(e);
     147              switch(heap.state(v)) {
     148              case Heap::PRE_HEAP:
     149                heap.push(v, d + _cost[e] - _potential[v]);
     150                _pred[v] = e;
     151                break;
     152              case Heap::IN_HEAP:
     153                nd = d + _cost[e] - _potential[v];
     154                if (nd < heap[v]) {
     155                  heap.decrease(v, nd);
     156                  _pred[v] = e;
     157                }
     158                break;
     159              case Heap::POST_HEAP:
     160                break;
     161              }
     162            }
     163          }
     164
     165          // Traversing incoming arcs
     166          for (InArcIt e(_graph, u); e != INVALID; ++e) {
     167            if (_flow[e] >= delta) {
     168              v = _graph.source(e);
     169              switch(heap.state(v)) {
     170              case Heap::PRE_HEAP:
     171                heap.push(v, d - _cost[e] - _potential[v]);
     172                _pred[v] = e;
     173                break;
     174              case Heap::IN_HEAP:
     175                nd = d - _cost[e] - _potential[v];
     176                if (nd < heap[v]) {
     177                  heap.decrease(v, nd);
     178                  _pred[v] = e;
     179                }
     180                break;
     181              case Heap::POST_HEAP:
     182                break;
     183              }
     184            }
     185          }
     186        }
     187        if (heap.empty()) return INVALID;
     188
     189        // Updating potentials of processed nodes
     190        Node t = heap.top();
     191        Cost t_dist = heap.prio();
     192        for (int i = 0; i < int(_proc_nodes.size()); ++i)
     193          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
     194
     195        return t;
     196      }
     197
     198    }; //class ResidualDijkstra
     199
     200  private:
     201
     202    // The digraph the algorithm runs on
     203    const Digraph &_graph;
     204    // The original lower bound map
     205    const LowerMap *_lower;
     206    // The modified capacity map
     207    CapacityArcMap _capacity;
     208    // The original cost map
     209    const CostMap &_cost;
     210    // The modified supply map
     211    SupplyNodeMap _supply;
     212    bool _valid_supply;
     213
     214    // Arc map of the current flow
     215    FlowMap *_flow;
     216    bool _local_flow;
     217    // Node map of the current potentials
     218    PotentialMap *_potential;
     219    bool _local_potential;
     220
     221    // The residual capacity map
     222    CapacityArcMap _res_cap;
     223    // The excess map
     224    SupplyNodeMap _excess;
     225    // The excess nodes (i.e. nodes with positive excess)
     226    std::vector<Node> _excess_nodes;
     227    // The deficit nodes (i.e. nodes with negative excess)
     228    std::vector<Node> _deficit_nodes;
     229
     230    // The delta parameter used for capacity scaling
     231    Capacity _delta;
     232    // The maximum number of phases
     233    int _phase_num;
     234
     235    // The pred arc map
     236    PredMap _pred;
     237    // Implementation of the Dijkstra algorithm for finding augmenting
     238    // shortest paths in the residual network
     239    ResidualDijkstra *_dijkstra;
     240
     241  public:
     242
     243    /// \brief General constructor (with lower bounds).
     244    ///
     245    /// General constructor (with lower bounds).
     246    ///
     247    /// \param digraph The digraph the algorithm runs on.
     248    /// \param lower The lower bounds of the arcs.
     249    /// \param capacity The capacities (upper bounds) of the arcs.
     250    /// \param cost The cost (length) values of the arcs.
     251    /// \param supply The supply values of the nodes (signed).
     252    CapacityScaling( const Digraph &digraph,
     253                     const LowerMap &lower,
     254                     const CapacityMap &capacity,
     255                     const CostMap &cost,
     256                     const SupplyMap &supply ) :
     257      _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
     258      _supply(digraph), _flow(NULL), _local_flow(false),
     259      _potential(NULL), _local_potential(false),
     260      _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL)
     261    {
     262      Supply sum = 0;
     263      for (NodeIt n(_graph); n != INVALID; ++n) {
     264        _supply[n] = supply[n];
     265        _excess[n] = supply[n];
     266        sum += supply[n];
     267      }
     268      _valid_supply = sum == 0;
     269      for (ArcIt a(_graph); a != INVALID; ++a) {
     270        _capacity[a] = capacity[a];
     271        _res_cap[a] = capacity[a];
     272      }
     273
     274      // Remove non-zero lower bounds
     275      typename LowerMap::Value lcap;
     276      for (ArcIt e(_graph); e != INVALID; ++e) {
     277        if ((lcap = lower[e]) != 0) {
     278          _capacity[e] -= lcap;
     279          _res_cap[e] -= lcap;
     280          _supply[_graph.source(e)] -= lcap;
     281          _supply[_graph.target(e)] += lcap;
     282          _excess[_graph.source(e)] -= lcap;
     283          _excess[_graph.target(e)] += lcap;
     284        }
     285      }
     286    }
     287/*
     288    /// \brief General constructor (without lower bounds).
     289    ///
     290    /// General constructor (without lower bounds).
     291    ///
     292    /// \param digraph The digraph the algorithm runs on.
     293    /// \param capacity The capacities (upper bounds) of the arcs.
     294    /// \param cost The cost (length) values of the arcs.
     295    /// \param supply The supply values of the nodes (signed).
     296    CapacityScaling( const Digraph &digraph,
     297                     const CapacityMap &capacity,
     298                     const CostMap &cost,
     299                     const SupplyMap &supply ) :
     300      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
     301      _supply(supply), _flow(NULL), _local_flow(false),
     302      _potential(NULL), _local_potential(false),
     303      _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL)
     304    {
     305      // Check the sum of supply values
     306      Supply sum = 0;
     307      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
     308      _valid_supply = sum == 0;
     309    }
     310
     311    /// \brief Simple constructor (with lower bounds).
     312    ///
     313    /// Simple constructor (with lower bounds).
     314    ///
     315    /// \param digraph The digraph the algorithm runs on.
     316    /// \param lower The lower bounds of the arcs.
     317    /// \param capacity The capacities (upper bounds) of the arcs.
     318    /// \param cost The cost (length) values of the arcs.
     319    /// \param s The source node.
     320    /// \param t The target node.
     321    /// \param flow_value The required amount of flow from node \c s
     322    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
     323    CapacityScaling( const Digraph &digraph,
     324                     const LowerMap &lower,
     325                     const CapacityMap &capacity,
     326                     const CostMap &cost,
     327                     Node s, Node t,
     328                     Supply flow_value ) :
     329      _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
     330      _supply(digraph, 0), _flow(NULL), _local_flow(false),
     331      _potential(NULL), _local_potential(false),
     332      _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
     333    {
     334      // Remove non-zero lower bounds
     335      _supply[s] = _excess[s] =  flow_value;
     336      _supply[t] = _excess[t] = -flow_value;
     337      typename LowerMap::Value lcap;
     338      for (ArcIt e(_graph); e != INVALID; ++e) {
     339        if ((lcap = lower[e]) != 0) {
     340          _capacity[e] -= lcap;
     341          _res_cap[e] -= lcap;
     342          _supply[_graph.source(e)] -= lcap;
     343          _supply[_graph.target(e)] += lcap;
     344          _excess[_graph.source(e)] -= lcap;
     345          _excess[_graph.target(e)] += lcap;
     346        }
     347      }
     348      _valid_supply = true;
     349    }
     350
     351    /// \brief Simple constructor (without lower bounds).
     352    ///
     353    /// Simple constructor (without lower bounds).
     354    ///
     355    /// \param digraph The digraph the algorithm runs on.
     356    /// \param capacity The capacities (upper bounds) of the arcs.
     357    /// \param cost The cost (length) values of the arcs.
     358    /// \param s The source node.
     359    /// \param t The target node.
     360    /// \param flow_value The required amount of flow from node \c s
     361    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
     362    CapacityScaling( const Digraph &digraph,
     363                     const CapacityMap &capacity,
     364                     const CostMap &cost,
     365                     Node s, Node t,
     366                     Supply flow_value ) :
     367      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
     368      _supply(digraph, 0), _flow(NULL), _local_flow(false),
     369      _potential(NULL), _local_potential(false),
     370      _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
     371    {
     372      _supply[s] = _excess[s] =  flow_value;
     373      _supply[t] = _excess[t] = -flow_value;
     374      _valid_supply = true;
     375    }
     376*/
     377    /// Destructor.
     378    ~CapacityScaling() {
     379      if (_local_flow) delete _flow;
     380      if (_local_potential) delete _potential;
     381      delete _dijkstra;
     382    }
     383
     384    /// \brief Set the flow map.
     385    ///
     386    /// Set the flow map.
     387    ///
     388    /// \return \c (*this)
     389    CapacityScaling& flowMap(FlowMap &map) {
     390      if (_local_flow) {
     391        delete _flow;
     392        _local_flow = false;
     393      }
     394      _flow = &map;
     395      return *this;
     396    }
     397
     398    /// \brief Set the potential map.
     399    ///
     400    /// Set the potential map.
     401    ///
     402    /// \return \c (*this)
     403    CapacityScaling& potentialMap(PotentialMap &map) {
     404      if (_local_potential) {
     405        delete _potential;
     406        _local_potential = false;
     407      }
     408      _potential = &map;
     409      return *this;
     410    }
     411
     412    /// \name Execution control
     413
     414    /// @{
     415
     416    /// \brief Run the algorithm.
     417    ///
     418    /// This function runs the algorithm.
     419    ///
     420    /// \param scaling Enable or disable capacity scaling.
     421    /// If the maximum arc capacity and/or the amount of total supply
     422    /// is rather small, the algorithm could be slightly faster without
     423    /// scaling.
     424    ///
     425    /// \return \c true if a feasible flow can be found.
     426    bool run(bool scaling = true) {
     427      return init(scaling) && start();
     428    }
     429
     430    /// @}
     431
     432    /// \name Query Functions
     433    /// The results of the algorithm can be obtained using these
     434    /// functions.\n
     435    /// \ref lemon::CapacityScaling::run() "run()" must be called before
     436    /// using them.
     437
     438    /// @{
     439
     440    /// \brief Return a const reference to the arc map storing the
     441    /// found flow.
     442    ///
     443    /// Return a const reference to the arc map storing the found flow.
     444    ///
     445    /// \pre \ref run() must be called before using this function.
     446    const FlowMap& flowMap() const {
     447      return *_flow;
     448    }
     449
     450    /// \brief Return a const reference to the node map storing the
     451    /// found potentials (the dual solution).
     452    ///
     453    /// Return a const reference to the node map storing the found
     454    /// potentials (the dual solution).
     455    ///
     456    /// \pre \ref run() must be called before using this function.
     457    const PotentialMap& potentialMap() const {
     458      return *_potential;
     459    }
     460
     461    /// \brief Return the flow on the given arc.
     462    ///
     463    /// Return the flow on the given arc.
     464    ///
     465    /// \pre \ref run() must be called before using this function.
     466    Capacity flow(const Arc& arc) const {
     467      return (*_flow)[arc];
     468    }
     469
     470    /// \brief Return the potential of the given node.
     471    ///
     472    /// Return the potential of the given node.
     473    ///
     474    /// \pre \ref run() must be called before using this function.
     475    Cost potential(const Node& node) const {
     476      return (*_potential)[node];
     477    }
     478
     479    /// \brief Return the total cost of the found flow.
     480    ///
     481    /// Return the total cost of the found flow. The complexity of the
     482    /// function is \f$ O(e) \f$.
     483    ///
     484    /// \pre \ref run() must be called before using this function.
     485    Cost totalCost() const {
     486      Cost c = 0;
     487      for (ArcIt e(_graph); e != INVALID; ++e)
     488        c += (*_flow)[e] * _cost[e];
     489      return c;
     490    }
     491
     492    /// @}
     493
     494  private:
     495
     496    /// Initialize the algorithm.
     497    bool init(bool scaling) {
     498      if (!_valid_supply) return false;
     499
     500      // Initializing maps
     501      if (!_flow) {
     502        _flow = new FlowMap(_graph);
     503        _local_flow = true;
     504      }
     505      if (!_potential) {
     506        _potential = new PotentialMap(_graph);
     507        _local_potential = true;
     508      }
     509      for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
     510      for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
     511
     512      _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
     513                                        _excess, *_potential, _pred );
     514
     515      // Initializing delta value
     516      if (scaling) {
     517        // With scaling
     518        Supply max_sup = 0, max_dem = 0;
     519        for (NodeIt n(_graph); n != INVALID; ++n) {
     520          if ( _supply[n] > max_sup) max_sup =  _supply[n];
     521          if (-_supply[n] > max_dem) max_dem = -_supply[n];
     522        }
     523        Capacity max_cap = 0;
     524        for (ArcIt e(_graph); e != INVALID; ++e) {
     525          if (_capacity[e] > max_cap) max_cap = _capacity[e];
     526        }
     527        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
     528        _phase_num = 0;
     529        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
     530          ++_phase_num;
     531      } else {
     532        // Without scaling
     533        _delta = 1;
     534      }
     535
     536      return true;
     537    }
     538
     539    bool start() {
     540      if (_delta > 1)
     541        return startWithScaling();
     542      else
     543        return startWithoutScaling();
     544    }
     545
     546    /// Execute the capacity scaling algorithm.
     547    bool startWithScaling() {
     548      // Processing capacity scaling phases
     549      Node s, t;
     550      int phase_cnt = 0;
     551      int factor = 4;
     552      while (true) {
     553        // Saturating all arcs not satisfying the optimality condition
     554        for (ArcIt e(_graph); e != INVALID; ++e) {
     555          Node u = _graph.source(e), v = _graph.target(e);
     556          Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
     557          if (c < 0 && _res_cap[e] >= _delta) {
     558            _excess[u] -= _res_cap[e];
     559            _excess[v] += _res_cap[e];
     560            (*_flow)[e] = _capacity[e];
     561            _res_cap[e] = 0;
     562          }
     563          else if (c > 0 && (*_flow)[e] >= _delta) {
     564            _excess[u] += (*_flow)[e];
     565            _excess[v] -= (*_flow)[e];
     566            (*_flow)[e] = 0;
     567            _res_cap[e] = _capacity[e];
     568          }
     569        }
     570
     571        // Finding excess nodes and deficit nodes
     572        _excess_nodes.clear();
     573        _deficit_nodes.clear();
     574        for (NodeIt n(_graph); n != INVALID; ++n) {
     575          if (_excess[n] >=  _delta) _excess_nodes.push_back(n);
     576          if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
     577        }
     578        int next_node = 0, next_def_node = 0;
     579
     580        // Finding augmenting shortest paths
     581        while (next_node < int(_excess_nodes.size())) {
     582          // Checking deficit nodes
     583          if (_delta > 1) {
     584            bool delta_deficit = false;
     585            for ( ; next_def_node < int(_deficit_nodes.size());
     586                    ++next_def_node ) {
     587              if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
     588                delta_deficit = true;
     589                break;
     590              }
     591            }
     592            if (!delta_deficit) break;
     593          }
     594
     595          // Running Dijkstra
     596          s = _excess_nodes[next_node];
     597          if ((t = _dijkstra->run(s, _delta)) == INVALID) {
     598            if (_delta > 1) {
     599              ++next_node;
     600              continue;
     601            }
     602            return false;
     603          }
     604
     605          // Augmenting along a shortest path from s to t.
     606          Capacity d = std::min(_excess[s], -_excess[t]);
     607          Node u = t;
     608          Arc e;
     609          if (d > _delta) {
     610            while ((e = _pred[u]) != INVALID) {
     611              Capacity rc;
     612              if (u == _graph.target(e)) {
     613                rc = _res_cap[e];
     614                u = _graph.source(e);
     615              } else {
     616                rc = (*_flow)[e];
     617                u = _graph.target(e);
     618              }
     619              if (rc < d) d = rc;
     620            }
     621          }
     622          u = t;
     623          while ((e = _pred[u]) != INVALID) {
     624            if (u == _graph.target(e)) {
     625              (*_flow)[e] += d;
     626              _res_cap[e] -= d;
     627              u = _graph.source(e);
     628            } else {
     629              (*_flow)[e] -= d;
     630              _res_cap[e] += d;
     631              u = _graph.target(e);
     632            }
     633          }
     634          _excess[s] -= d;
     635          _excess[t] += d;
     636
     637          if (_excess[s] < _delta) ++next_node;
     638        }
     639
     640        if (_delta == 1) break;
     641        if (++phase_cnt > _phase_num / 4) factor = 2;
     642        _delta = _delta <= factor ? 1 : _delta / factor;
     643      }
     644
     645      // Handling non-zero lower bounds
     646      if (_lower) {
     647        for (ArcIt e(_graph); e != INVALID; ++e)
     648          (*_flow)[e] += (*_lower)[e];
     649      }
     650      return true;
     651    }
     652
     653    /// Execute the successive shortest path algorithm.
     654    bool startWithoutScaling() {
     655      // Finding excess nodes
     656      for (NodeIt n(_graph); n != INVALID; ++n)
     657        if (_excess[n] > 0) _excess_nodes.push_back(n);
     658      if (_excess_nodes.size() == 0) return true;
     659      int next_node = 0;
     660
     661      // Finding shortest paths
     662      Node s, t;
     663      while ( _excess[_excess_nodes[next_node]] > 0 ||
     664              ++next_node < int(_excess_nodes.size()) )
     665      {
     666        // Running Dijkstra
     667        s = _excess_nodes[next_node];
     668        if ((t = _dijkstra->run(s)) == INVALID) return false;
     669
     670        // Augmenting along a shortest path from s to t
     671        Capacity d = std::min(_excess[s], -_excess[t]);
     672        Node u = t;
     673        Arc e;
     674        if (d > 1) {
     675          while ((e = _pred[u]) != INVALID) {
     676            Capacity rc;
     677            if (u == _graph.target(e)) {
     678              rc = _res_cap[e];
     679              u = _graph.source(e);
     680            } else {
     681              rc = (*_flow)[e];
     682              u = _graph.target(e);
     683            }
     684            if (rc < d) d = rc;
     685          }
     686        }
     687        u = t;
     688        while ((e = _pred[u]) != INVALID) {
     689          if (u == _graph.target(e)) {
     690            (*_flow)[e] += d;
     691            _res_cap[e] -= d;
     692            u = _graph.source(e);
     693          } else {
     694            (*_flow)[e] -= d;
     695            _res_cap[e] += d;
     696            u = _graph.target(e);
     697          }
     698        }
     699        _excess[s] -= d;
     700        _excess[t] += d;
     701      }
     702
     703      // Handling non-zero lower bounds
     704      if (_lower) {
     705        for (ArcIt e(_graph); e != INVALID; ++e)
     706          (*_flow)[e] += (*_lower)[e];
     707      }
     708      return true;
     709    }
     710
     711  }; //class CapacityScaling
     712
     713  ///@}
     714
     715} //namespace lemon
     716
     717#endif //LEMON_CAPACITY_SCALING_H