COIN-OR::LEMON - Graph Library

About

Downloads

Documentation

Development

Issue Tracker

Browse Tickets

Custom Query

Roadmap

Browse Source

Local Info

Mailing Lists

Search

Changeset 872:fa6f37d7a25b in lemon


Ignore:
Timestamp:
11/12/09 23:26:13 (14 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Rebase:
39653339343431333337306539356364643765343037643134313235643439386231623630306163
Message:

Entirely rework CapacityScaling? (#180)

  • Use the new interface similarly to NetworkSimplex?.
  • Rework the implementation using an efficient internal structure for handling the residual network. This improvement made the code much faster (up to 2-5 times faster on large graphs).
  • Handle GEQ supply type (LEQ is not supported).
  • Handle negative costs for arcs of finite capacity. (Note that this algorithm cannot handle arcs of negative cost and infinite upper bound, thus it returns UNBOUNDED if such an arc exists.)
  • Extend the documentation.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • lemon/capacity_scaling.h

    r871 r872  
    2020#define LEMON_CAPACITY_SCALING_H
    2121
    22 /// \ingroup min_cost_flow
     22/// \ingroup min_cost_flow_algs
    2323///
    2424/// \file
    25 /// \brief Capacity scaling algorithm for finding a minimum cost flow.
     25/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
    2626
    2727#include <vector>
     28#include <limits>
     29#include <lemon/core.h>
    2830#include <lemon/bin_heap.h>
    2931
    3032namespace lemon {
    3133
    32   /// \addtogroup min_cost_flow
     34  /// \addtogroup min_cost_flow_algs
    3335  /// @{
    3436
    35   /// \brief Implementation of the capacity scaling algorithm for
    36   /// finding a minimum cost flow.
     37  /// \brief Implementation of the Capacity Scaling algorithm for
     38  /// finding a \ref min_cost_flow "minimum cost flow".
    3739  ///
    3840  /// \ref CapacityScaling implements the capacity scaling version
    39   /// of the successive shortest path algorithm for finding a minimum
    40   /// cost flow.
     41  /// of the successive shortest path algorithm for finding a
     42  /// \ref min_cost_flow "minimum cost flow". It is an efficient dual
     43  /// solution method.
    4144  ///
    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.
     45  /// Most of the parameters of the problem (except for the digraph)
     46  /// can be given using separate functions, and the algorithm can be
     47  /// executed using the \ref run() function. If some parameters are not
     48  /// specified, then default values will be used.
    4749  ///
    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.
     50  /// \tparam GR The digraph type the algorithm runs on.
     51  /// \tparam V The value type used for flow amounts, capacity bounds
     52  /// and supply values in the algorithm. By default it is \c int.
     53  /// \tparam C The value type used for costs and potentials in the
     54  /// algorithm. By default it is the same as \c V.
    5355  ///
    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> >
     56  /// \warning Both value types must be signed and all input data must
     57  /// be integer.
     58  /// \warning This algorithm does not support negative costs for such
     59  /// arcs that have infinite upper bound.
     60  template <typename GR, typename V = int, typename C = V>
    6061  class CapacityScaling
    6162  {
    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 
    7163  public:
    7264
    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 
     65    /// The type of the flow amounts, capacity bounds and supply values
     66    typedef V Value;
     67    /// The type of the arc costs
     68    typedef C Cost;
     69
     70  public:
     71
     72    /// \brief Problem type constants for the \c run() function.
     73    ///
     74    /// Enum type containing the problem type constants that can be
     75    /// returned by the \ref run() function of the algorithm.
     76    enum ProblemType {
     77      /// The problem has no feasible solution (flow).
     78      INFEASIBLE,
     79      /// The problem has optimal solution (i.e. it is feasible and
     80      /// bounded), and the algorithm has found optimal flow and node
     81      /// potentials (primal and dual solutions).
     82      OPTIMAL,
     83      /// The digraph contains an arc of negative cost and infinite
     84      /// upper bound. It means that the objective function is unbounded
     85      /// on that arc, however note that it could actually be bounded
     86      /// over the feasible flows, but this algroithm cannot handle
     87      /// these cases.
     88      UNBOUNDED
     89    };
     90 
    7891  private:
    7992
    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.
     93    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     94
     95    typedef std::vector<Arc> ArcVector;
     96    typedef std::vector<Node> NodeVector;
     97    typedef std::vector<int> IntVector;
     98    typedef std::vector<bool> BoolVector;
     99    typedef std::vector<Value> ValueVector;
     100    typedef std::vector<Cost> CostVector;
     101
     102  private:
     103
     104    // Data related to the underlying digraph
     105    const GR &_graph;
     106    int _node_num;
     107    int _arc_num;
     108    int _res_arc_num;
     109    int _root;
     110
     111    // Parameters of the problem
     112    bool _have_lower;
     113    Value _sum_supply;
     114
     115    // Data structures for storing the digraph
     116    IntNodeMap _node_id;
     117    IntArcMap _arc_idf;
     118    IntArcMap _arc_idb;
     119    IntVector _first_out;
     120    BoolVector _forward;
     121    IntVector _source;
     122    IntVector _target;
     123    IntVector _reverse;
     124
     125    // Node and arc data
     126    ValueVector _lower;
     127    ValueVector _upper;
     128    CostVector _cost;
     129    ValueVector _supply;
     130
     131    ValueVector _res_cap;
     132    CostVector _pi;
     133    ValueVector _excess;
     134    IntVector _excess_nodes;
     135    IntVector _deficit_nodes;
     136
     137    Value _delta;
     138    int _phase_num;
     139    IntVector _pred;
     140
     141  public:
     142 
     143    /// \brief Constant for infinite upper bounds (capacities).
     144    ///
     145    /// Constant for infinite upper bounds (capacities).
     146    /// It is \c std::numeric_limits<Value>::infinity() if available,
     147    /// \c std::numeric_limits<Value>::max() otherwise.
     148    const Value INF;
     149
     150  private:
     151
     152    // Special implementation of the Dijkstra algorithm for finding
     153    // shortest paths in the residual network of the digraph with
     154    // respect to the reduced arc costs and modifying the node
     155    // potentials according to the found distance labels.
    88156    class ResidualDijkstra
    89157    {
    90       typedef typename Digraph::template NodeMap<int> HeapCrossRef;
     158      typedef RangeMap<int> HeapCrossRef;
    91159      typedef BinHeap<Cost, HeapCrossRef> Heap;
    92160
    93161    private:
    94162
    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 
     163      int _node_num;
     164      const IntVector &_first_out;
     165      const IntVector &_target;
     166      const CostVector &_cost;
     167      const ValueVector &_res_cap;
     168      const ValueVector &_excess;
     169      CostVector &_pi;
     170      IntVector &_pred;
     171     
     172      IntVector _proc_nodes;
     173      CostVector _dist;
     174     
    112175    public:
    113176
    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)
     177      ResidualDijkstra(CapacityScaling& cs) :
     178        _node_num(cs._node_num), _first_out(cs._first_out),
     179        _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
     180        _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
     181        _dist(cs._node_num)
    125182      {}
    126183
    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);
     184      int run(int s, Value delta = 1) {
     185        HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP);
    130186        Heap heap(heap_cross_ref);
    131187        heap.push(s, 0);
    132         _pred[s] = INVALID;
     188        _pred[s] = -1;
    133189        _proc_nodes.clear();
    134190
    135         // Processing nodes
     191        // Process nodes
    136192        while (!heap.empty() && _excess[heap.top()] > -delta) {
    137           Node u = heap.top(), v;
    138           Cost d = heap.prio() + _potential[u], nd;
     193          int u = heap.top(), v;
     194          Cost d = heap.prio() + _pi[u], dn;
    139195          _dist[u] = heap.prio();
     196          _proc_nodes.push_back(u);
    140197          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)) {
     198
     199          // Traverse outgoing residual arcs
     200          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
     201            if (_res_cap[a] < delta) continue;
     202            v = _target[a];
     203            switch (heap.state(v)) {
    148204              case Heap::PRE_HEAP:
    149                 heap.push(v, d + _cost[e] - _potential[v]);
    150                 _pred[v] = e;
     205                heap.push(v, d + _cost[a] - _pi[v]);
     206                _pred[v] = a;
    151207                break;
    152208              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;
     209                dn = d + _cost[a] - _pi[v];
     210                if (dn < heap[v]) {
     211                  heap.decrease(v, dn);
     212                  _pred[v] = a;
    157213                }
    158214                break;
    159215              case Heap::POST_HEAP:
    160216                break;
    161               }
    162217            }
    163218          }
    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;
     219        }
     220        if (heap.empty()) return -1;
     221
     222        // Update potentials of processed nodes
     223        int t = heap.top();
     224        Cost dt = heap.prio();
     225        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
     226          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
     227        }
    194228
    195229        return t;
    … …  
    198232    }; //class ResidualDijkstra
    199233
    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 
    241234  public:
    242235
    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)
     236    /// \brief Constructor.
     237    ///
     238    /// The constructor of the class.
     239    ///
     240    /// \param graph The digraph the algorithm runs on.
     241    CapacityScaling(const GR& graph) :
     242      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
     243      INF(std::numeric_limits<Value>::has_infinity ?
     244          std::numeric_limits<Value>::infinity() :
     245          std::numeric_limits<Value>::max())
    261246    {
    262       Supply sum = 0;
     247      // Check the value types
     248      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
     249        "The flow type of CapacityScaling must be signed");
     250      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
     251        "The cost type of CapacityScaling must be signed");
     252
     253      // Resize vectors
     254      _node_num = countNodes(_graph);
     255      _arc_num = countArcs(_graph);
     256      _res_arc_num = 2 * (_arc_num + _node_num);
     257      _root = _node_num;
     258      ++_node_num;
     259
     260      _first_out.resize(_node_num + 1);
     261      _forward.resize(_res_arc_num);
     262      _source.resize(_res_arc_num);
     263      _target.resize(_res_arc_num);
     264      _reverse.resize(_res_arc_num);
     265
     266      _lower.resize(_res_arc_num);
     267      _upper.resize(_res_arc_num);
     268      _cost.resize(_res_arc_num);
     269      _supply.resize(_node_num);
     270     
     271      _res_cap.resize(_res_arc_num);
     272      _pi.resize(_node_num);
     273      _excess.resize(_node_num);
     274      _pred.resize(_node_num);
     275
     276      // Copy the graph
     277      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
     278      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     279        _node_id[n] = i;
     280      }
     281      i = 0;
     282      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     283        _first_out[i] = j;
     284        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
     285          _arc_idf[a] = j;
     286          _forward[j] = true;
     287          _source[j] = i;
     288          _target[j] = _node_id[_graph.runningNode(a)];
     289        }
     290        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
     291          _arc_idb[a] = j;
     292          _forward[j] = false;
     293          _source[j] = i;
     294          _target[j] = _node_id[_graph.runningNode(a)];
     295        }
     296        _forward[j] = false;
     297        _source[j] = i;
     298        _target[j] = _root;
     299        _reverse[j] = k;
     300        _forward[k] = true;
     301        _source[k] = _root;
     302        _target[k] = i;
     303        _reverse[k] = j;
     304        ++j; ++k;
     305      }
     306      _first_out[i] = j;
     307      _first_out[_node_num] = k;
     308      for (ArcIt a(_graph); a != INVALID; ++a) {
     309        int fi = _arc_idf[a];
     310        int bi = _arc_idb[a];
     311        _reverse[fi] = bi;
     312        _reverse[bi] = fi;
     313      }
     314     
     315      // Reset parameters
     316      reset();
     317    }
     318
     319    /// \name Parameters
     320    /// The parameters of the algorithm can be specified using these
     321    /// functions.
     322
     323    /// @{
     324
     325    /// \brief Set the lower bounds on the arcs.
     326    ///
     327    /// This function sets the lower bounds on the arcs.
     328    /// If it is not used before calling \ref run(), the lower bounds
     329    /// will be set to zero on all arcs.
     330    ///
     331    /// \param map An arc map storing the lower bounds.
     332    /// Its \c Value type must be convertible to the \c Value type
     333    /// of the algorithm.
     334    ///
     335    /// \return <tt>(*this)</tt>
     336    template <typename LowerMap>
     337    CapacityScaling& lowerMap(const LowerMap& map) {
     338      _have_lower = true;
     339      for (ArcIt a(_graph); a != INVALID; ++a) {
     340        _lower[_arc_idf[a]] = map[a];
     341        _lower[_arc_idb[a]] = map[a];
     342      }
     343      return *this;
     344    }
     345
     346    /// \brief Set the upper bounds (capacities) on the arcs.
     347    ///
     348    /// This function sets the upper bounds (capacities) on the arcs.
     349    /// If it is not used before calling \ref run(), the upper bounds
     350    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     351    /// unbounded from above on each arc).
     352    ///
     353    /// \param map An arc map storing the upper bounds.
     354    /// Its \c Value type must be convertible to the \c Value type
     355    /// of the algorithm.
     356    ///
     357    /// \return <tt>(*this)</tt>
     358    template<typename UpperMap>
     359    CapacityScaling& upperMap(const UpperMap& map) {
     360      for (ArcIt a(_graph); a != INVALID; ++a) {
     361        _upper[_arc_idf[a]] = map[a];
     362      }
     363      return *this;
     364    }
     365
     366    /// \brief Set the costs of the arcs.
     367    ///
     368    /// This function sets the costs of the arcs.
     369    /// If it is not used before calling \ref run(), the costs
     370    /// will be set to \c 1 on all arcs.
     371    ///
     372    /// \param map An arc map storing the costs.
     373    /// Its \c Value type must be convertible to the \c Cost type
     374    /// of the algorithm.
     375    ///
     376    /// \return <tt>(*this)</tt>
     377    template<typename CostMap>
     378    CapacityScaling& costMap(const CostMap& map) {
     379      for (ArcIt a(_graph); a != INVALID; ++a) {
     380        _cost[_arc_idf[a]] =  map[a];
     381        _cost[_arc_idb[a]] = -map[a];
     382      }
     383      return *this;
     384    }
     385
     386    /// \brief Set the supply values of the nodes.
     387    ///
     388    /// This function sets the supply values of the nodes.
     389    /// If neither this function nor \ref stSupply() is used before
     390    /// calling \ref run(), the supply of each node will be set to zero.
     391    ///
     392    /// \param map A node map storing the supply values.
     393    /// Its \c Value type must be convertible to the \c Value type
     394    /// of the algorithm.
     395    ///
     396    /// \return <tt>(*this)</tt>
     397    template<typename SupplyMap>
     398    CapacityScaling& supplyMap(const SupplyMap& map) {
    263399      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.
     400        _supply[_node_id[n]] = map[n];
     401      }
     402      return *this;
     403    }
     404
     405    /// \brief Set single source and target nodes and a supply value.
     406    ///
     407    /// This function sets a single source node and a single target node
     408    /// and the required flow value.
     409    /// If neither this function nor \ref supplyMap() is used before
     410    /// calling \ref run(), the supply of each node will be set to zero.
     411    ///
     412    /// Using this function has the same effect as using \ref supplyMap()
     413    /// with such a map in which \c k is assigned to \c s, \c -k is
     414    /// assigned to \c t and all other nodes have zero supply value.
     415    ///
    319416    /// \param s The source node.
    320417    /// \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;
     418    /// \param k The required amount of flow from node \c s to node \c t
     419    /// (i.e. the supply of \c s and the demand of \c t).
     420    ///
     421    /// \return <tt>(*this)</tt>
     422    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
     423      for (int i = 0; i != _node_num; ++i) {
     424        _supply[i] = 0;
     425      }
     426      _supply[_node_id[s]] =  k;
     427      _supply[_node_id[t]] = -k;
    395428      return *this;
    396429    }
    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     }
     430   
     431    /// @}
    411432
    412433    /// \name Execution control
    … …  
    417438    ///
    418439    /// This function runs the algorithm.
     440    /// The paramters can be specified using functions \ref lowerMap(),
     441    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     442    /// For example,
     443    /// \code
     444    ///   CapacityScaling<ListDigraph> cs(graph);
     445    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     446    ///     .supplyMap(sup).run();
     447    /// \endcode
     448    ///
     449    /// This function can be called more than once. All the parameters
     450    /// that have been given are kept for the next call, unless
     451    /// \ref reset() is called, thus only the modified parameters
     452    /// have to be set again. See \ref reset() for examples.
     453    /// However the underlying digraph must not be modified after this
     454    /// class have been constructed, since it copies the digraph.
    419455    ///
    420456    /// \param scaling Enable or disable capacity scaling.
    421     /// If the maximum arc capacity and/or the amount of total supply
     457    /// If the maximum upper bound and/or the amount of total supply
    422458    /// is rather small, the algorithm could be slightly faster without
    423459    /// scaling.
    424460    ///
    425     /// \return \c true if a feasible flow can be found.
    426     bool run(bool scaling = true) {
    427       return init(scaling) && start();
     461    /// \return \c INFEASIBLE if no feasible flow exists,
     462    /// \n \c OPTIMAL if the problem has optimal solution
     463    /// (i.e. it is feasible and bounded), and the algorithm has found
     464    /// optimal flow and node potentials (primal and dual solutions),
     465    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
     466    /// and infinite upper bound. It means that the objective function
     467    /// is unbounded on that arc, however note that it could actually be
     468    /// bounded over the feasible flows, but this algroithm cannot handle
     469    /// these cases.
     470    ///
     471    /// \see ProblemType
     472    ProblemType run(bool scaling = true) {
     473      ProblemType pt = init(scaling);
     474      if (pt != OPTIMAL) return pt;
     475      return start();
     476    }
     477
     478    /// \brief Reset all the parameters that have been given before.
     479    ///
     480    /// This function resets all the paramaters that have been given
     481    /// before using functions \ref lowerMap(), \ref upperMap(),
     482    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     483    ///
     484    /// It is useful for multiple run() calls. If this function is not
     485    /// used, all the parameters given before are kept for the next
     486    /// \ref run() call.
     487    /// However the underlying digraph must not be modified after this
     488    /// class have been constructed, since it copies and extends the graph.
     489    ///
     490    /// For example,
     491    /// \code
     492    ///   CapacityScaling<ListDigraph> cs(graph);
     493    ///
     494    ///   // First run
     495    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     496    ///     .supplyMap(sup).run();
     497    ///
     498    ///   // Run again with modified cost map (reset() is not called,
     499    ///   // so only the cost map have to be set again)
     500    ///   cost[e] += 100;
     501    ///   cs.costMap(cost).run();
     502    ///
     503    ///   // Run again from scratch using reset()
     504    ///   // (the lower bounds will be set to zero on all arcs)
     505    ///   cs.reset();
     506    ///   cs.upperMap(capacity).costMap(cost)
     507    ///     .supplyMap(sup).run();
     508    /// \endcode
     509    ///
     510    /// \return <tt>(*this)</tt>
     511    CapacityScaling& reset() {
     512      for (int i = 0; i != _node_num; ++i) {
     513        _supply[i] = 0;
     514      }
     515      for (int j = 0; j != _res_arc_num; ++j) {
     516        _lower[j] = 0;
     517        _upper[j] = INF;
     518        _cost[j] = _forward[j] ? 1 : -1;
     519      }
     520      _have_lower = false;
     521      return *this;
    428522    }
    429523
    … …  
    433527    /// The results of the algorithm can be obtained using these
    434528    /// functions.\n
    435     /// \ref lemon::CapacityScaling::run() "run()" must be called before
    436     /// using them.
     529    /// The \ref run() function must be called before using them.
    437530
    438531    /// @{
    439532
    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.
     533    /// \brief Return the total cost of the found flow.
     534    ///
     535    /// This function returns the total cost of the found flow.
     536    /// Its complexity is O(e).
     537    ///
     538    /// \note The return type of the function can be specified as a
     539    /// template parameter. For example,
     540    /// \code
     541    ///   cs.totalCost<double>();
     542    /// \endcode
     543    /// It is useful if the total cost cannot be stored in the \c Cost
     544    /// type of the algorithm, which is the default return type of the
     545    /// function.
    444546    ///
    445547    /// \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).
     548    template <typename Number>
     549    Number totalCost() const {
     550      Number c = 0;
     551      for (ArcIt a(_graph); a != INVALID; ++a) {
     552        int i = _arc_idb[a];
     553        c += static_cast<Number>(_res_cap[i]) *
     554             (-static_cast<Number>(_cost[i]));
     555      }
     556      return c;
     557    }
     558
     559#ifndef DOXYGEN
     560    Cost totalCost() const {
     561      return totalCost<Cost>();
     562    }
     563#endif
     564
     565    /// \brief Return the flow on the given arc.
     566    ///
     567    /// This function returns the flow on the given arc.
    455568    ///
    456569    /// \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.
     570    Value flow(const Arc& a) const {
     571      return _res_cap[_arc_idb[a]];
     572    }
     573
     574    /// \brief Return the flow map (the primal solution).
     575    ///
     576    /// This function copies the flow value on each arc into the given
     577    /// map. The \c Value type of the algorithm must be convertible to
     578    /// the \c Value type of the map.
    464579    ///
    465580    /// \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.
     581    template <typename FlowMap>
     582    void flowMap(FlowMap &map) const {
     583      for (ArcIt a(_graph); a != INVALID; ++a) {
     584        map.set(a, _res_cap[_arc_idb[a]]);
     585      }
     586    }
     587
     588    /// \brief Return the potential (dual value) of the given node.
     589    ///
     590    /// This function returns the potential (dual value) of the
     591    /// given node.
    473592    ///
    474593    /// \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$.
     594    Cost potential(const Node& n) const {
     595      return _pi[_node_id[n]];
     596    }
     597
     598    /// \brief Return the potential map (the dual solution).
     599    ///
     600    /// This function copies the potential (dual value) of each node
     601    /// into the given map.
     602    /// The \c Cost type of the algorithm must be convertible to the
     603    /// \c Value type of the map.
    483604    ///
    484605    /// \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;
     606    template <typename PotentialMap>
     607    void potentialMap(PotentialMap &map) const {
     608      for (NodeIt n(_graph); n != INVALID; ++n) {
     609        map.set(n, _pi[_node_id[n]]);
     610      }
    490611    }
    491612
    … …  
    494615  private:
    495616
    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
     617    // Initialize the algorithm
     618    ProblemType init(bool scaling) {
     619      if (_node_num == 0) return INFEASIBLE;
     620
     621      // Check the sum of supply values
     622      _sum_supply = 0;
     623      for (int i = 0; i != _root; ++i) {
     624        _sum_supply += _supply[i];
     625      }
     626      if (_sum_supply > 0) return INFEASIBLE;
     627     
     628      // Initialize maps
     629      for (int i = 0; i != _root; ++i) {
     630        _pi[i] = 0;
     631        _excess[i] = _supply[i];
     632      }
     633
     634      // Remove non-zero lower bounds
     635      if (_have_lower) {
     636        for (int i = 0; i != _root; ++i) {
     637          for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
     638            if (_forward[j]) {
     639              Value c = _lower[j];
     640              if (c >= 0) {
     641                _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
     642              } else {
     643                _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
     644              }
     645              _excess[i] -= c;
     646              _excess[_target[j]] += c;
     647            } else {
     648              _res_cap[j] = 0;
     649            }
     650          }
     651        }
     652      } else {
     653        for (int j = 0; j != _res_arc_num; ++j) {
     654          _res_cap[j] = _forward[j] ? _upper[j] : 0;
     655        }
     656      }
     657
     658      // Handle negative costs
     659      for (int u = 0; u != _root; ++u) {
     660        for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
     661          Value rc = _res_cap[a];
     662          if (_cost[a] < 0 && rc > 0) {
     663            if (rc == INF) return UNBOUNDED;
     664            _excess[u] -= rc;
     665            _excess[_target[a]] += rc;
     666            _res_cap[a] = 0;
     667            _res_cap[_reverse[a]] += rc;
     668          }
     669        }
     670      }
     671     
     672      // Handle GEQ supply type
     673      if (_sum_supply < 0) {
     674        _pi[_root] = 0;
     675        _excess[_root] = -_sum_supply;
     676        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
     677          int u = _target[a];
     678          if (_excess[u] < 0) {
     679            _res_cap[a] = -_excess[u] + 1;
     680          } else {
     681            _res_cap[a] = 1;
     682          }
     683          _res_cap[_reverse[a]] = 0;
     684          _cost[a] = 0;
     685          _cost[_reverse[a]] = 0;
     686        }
     687      } else {
     688        _pi[_root] = 0;
     689        _excess[_root] = 0;
     690        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
     691          _res_cap[a] = 1;
     692          _res_cap[_reverse[a]] = 0;
     693          _cost[a] = 0;
     694          _cost[_reverse[a]] = 0;
     695        }
     696      }
     697
     698      // Initialize delta value
    516699      if (scaling) {
    517700        // 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];
     701        Value max_sup = 0, max_dem = 0;
     702        for (int i = 0; i != _node_num; ++i) {
     703          if ( _excess[i] > max_sup) max_sup =  _excess[i];
     704          if (-_excess[i] > max_dem) max_dem = -_excess[i];
     705        }
     706        Value max_cap = 0;
     707        for (int j = 0; j != _res_arc_num; ++j) {
     708          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
    526709        }
    527710        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
    … …  
    534717      }
    535718
    536       return true;
    537     }
    538 
    539     bool start() {
     719      return OPTIMAL;
     720    }
     721
     722    ProblemType start() {
     723      // Execute the algorithm
     724      ProblemType pt;
    540725      if (_delta > 1)
    541         return startWithScaling();
     726        pt = startWithScaling();
    542727      else
    543         return startWithoutScaling();
    544     }
    545 
    546     /// Execute the capacity scaling algorithm.
    547     bool startWithScaling() {
    548       // Processing capacity scaling phases
    549       Node s, t;
     728        pt = startWithoutScaling();
     729
     730      // Handle non-zero lower bounds
     731      if (_have_lower) {
     732        for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
     733          if (!_forward[j]) _res_cap[j] += _lower[j];
     734        }
     735      }
     736
     737      // Shift potentials if necessary
     738      Cost pr = _pi[_root];
     739      if (_sum_supply < 0 || pr > 0) {
     740        for (int i = 0; i != _node_num; ++i) {
     741          _pi[i] -= pr;
     742        }       
     743      }
     744     
     745      return pt;
     746    }
     747
     748    // Execute the capacity scaling algorithm
     749    ProblemType startWithScaling() {
     750      // Process capacity scaling phases
     751      int s, t;
    550752      int phase_cnt = 0;
    551753      int factor = 4;
     754      ResidualDijkstra _dijkstra(*this);
    552755      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
     756        // Saturate all arcs not satisfying the optimality condition
     757        for (int u = 0; u != _node_num; ++u) {
     758          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
     759            int v = _target[a];
     760            Cost c = _cost[a] + _pi[u] - _pi[v];
     761            Value rc = _res_cap[a];
     762            if (c < 0 && rc >= _delta) {
     763              _excess[u] -= rc;
     764              _excess[v] += rc;
     765              _res_cap[a] = 0;
     766              _res_cap[_reverse[a]] += rc;
     767            }
     768          }
     769        }
     770
     771        // Find excess nodes and deficit nodes
    572772        _excess_nodes.clear();
    573773        _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);
     774        for (int u = 0; u != _node_num; ++u) {
     775          if (_excess[u] >=  _delta) _excess_nodes.push_back(u);
     776          if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
    577777        }
    578778        int next_node = 0, next_def_node = 0;
    579779
    580         // Finding augmenting shortest paths
     780        // Find augmenting shortest paths
    581781        while (next_node < int(_excess_nodes.size())) {
    582           // Checking deficit nodes
     782          // Check deficit nodes
    583783          if (_delta > 1) {
    584784            bool delta_deficit = false;
    … …  
    593793          }
    594794
    595           // Running Dijkstra
     795          // Run Dijkstra in the residual network
    596796          s = _excess_nodes[next_node];
    597           if ((t = _dijkstra->run(s, _delta)) == INVALID) {
     797          if ((t = _dijkstra.run(s, _delta)) == -1) {
    598798            if (_delta > 1) {
    599799              ++next_node;
    600800              continue;
    601801            }
    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;
     802            return INFEASIBLE;
     803          }
     804
     805          // Augment along a shortest path from s to t
     806          Value d = std::min(_excess[s], -_excess[t]);
     807          int u = t;
     808          int a;
    609809          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;
     810            while ((a = _pred[u]) != -1) {
     811              if (_res_cap[a] < d) d = _res_cap[a];
     812              u = _source[a];
    620813            }
    621814          }
    622815          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             }
     816          while ((a = _pred[u]) != -1) {
     817            _res_cap[a] -= d;
     818            _res_cap[_reverse[a]] += d;
     819            u = _source[a];
    633820          }
    634821          _excess[s] -= d;
    … …  
    639826
    640827        if (_delta == 1) break;
    641         if (++phase_cnt > _phase_num / 4) factor = 2;
     828        if (++phase_cnt == _phase_num / 4) factor = 2;
    642829        _delta = _delta <= factor ? 1 : _delta / factor;
    643830      }
    644831
    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;
     832      return OPTIMAL;
     833    }
     834
     835    // Execute the successive shortest path algorithm
     836    ProblemType startWithoutScaling() {
     837      // Find excess nodes
     838      _excess_nodes.clear();
     839      for (int i = 0; i != _node_num; ++i) {
     840        if (_excess[i] > 0) _excess_nodes.push_back(i);
     841      }
     842      if (_excess_nodes.size() == 0) return OPTIMAL;
    659843      int next_node = 0;
    660844
    661       // Finding shortest paths
    662       Node s, t;
     845      // Find shortest paths
     846      int s, t;
     847      ResidualDijkstra _dijkstra(*this);
    663848      while ( _excess[_excess_nodes[next_node]] > 0 ||
    664849              ++next_node < int(_excess_nodes.size()) )
    665850      {
    666         // Running Dijkstra
     851        // Run Dijkstra in the residual network
    667852        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;
     853        if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
     854
     855        // Augment along a shortest path from s to t
     856        Value d = std::min(_excess[s], -_excess[t]);
     857        int u = t;
     858        int a;
    674859        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;
     860          while ((a = _pred[u]) != -1) {
     861            if (_res_cap[a] < d) d = _res_cap[a];
     862            u = _source[a];
    685863          }
    686864        }
    687865        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           }
     866        while ((a = _pred[u]) != -1) {
     867          _res_cap[a] -= d;
     868          _res_cap[_reverse[a]] += d;
     869          u = _source[a];
    698870        }
    699871        _excess[s] -= d;
    … …  
    701873      }
    702874
    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;
     875      return OPTIMAL;
    709876    }
    710877
Note: See TracChangeset for help on using the changeset viewer.