lemon/capacity_scaling.h
author Peter Kovacs <kpeter@inf.elte.hu>
Thu, 12 Nov 2009 23:26:13 +0100
changeset 806 fa6f37d7a25b
parent 805 d3e32a777d0b
child 807 78071e00de00
permissions -rw-r--r--
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.
     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_algs
    23 ///
    24 /// \file
    25 /// \brief Capacity Scaling algorithm for finding a minimum cost flow.
    26 
    27 #include <vector>
    28 #include <limits>
    29 #include <lemon/core.h>
    30 #include <lemon/bin_heap.h>
    31 
    32 namespace lemon {
    33 
    34   /// \addtogroup min_cost_flow_algs
    35   /// @{
    36 
    37   /// \brief Implementation of the Capacity Scaling algorithm for
    38   /// finding a \ref min_cost_flow "minimum cost flow".
    39   ///
    40   /// \ref CapacityScaling implements the capacity scaling version
    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.
    44   ///
    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.
    49   ///
    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.
    55   ///
    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>
    61   class CapacityScaling
    62   {
    63   public:
    64 
    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   
    91   private:
    92 
    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.
   156     class ResidualDijkstra
   157     {
   158       typedef RangeMap<int> HeapCrossRef;
   159       typedef BinHeap<Cost, HeapCrossRef> Heap;
   160 
   161     private:
   162 
   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       
   175     public:
   176 
   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)
   182       {}
   183 
   184       int run(int s, Value delta = 1) {
   185         HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP);
   186         Heap heap(heap_cross_ref);
   187         heap.push(s, 0);
   188         _pred[s] = -1;
   189         _proc_nodes.clear();
   190 
   191         // Process nodes
   192         while (!heap.empty() && _excess[heap.top()] > -delta) {
   193           int u = heap.top(), v;
   194           Cost d = heap.prio() + _pi[u], dn;
   195           _dist[u] = heap.prio();
   196           _proc_nodes.push_back(u);
   197           heap.pop();
   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)) {
   204               case Heap::PRE_HEAP:
   205                 heap.push(v, d + _cost[a] - _pi[v]);
   206                 _pred[v] = a;
   207                 break;
   208               case Heap::IN_HEAP:
   209                 dn = d + _cost[a] - _pi[v];
   210                 if (dn < heap[v]) {
   211                   heap.decrease(v, dn);
   212                   _pred[v] = a;
   213                 }
   214                 break;
   215               case Heap::POST_HEAP:
   216                 break;
   217             }
   218           }
   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         }
   228 
   229         return t;
   230       }
   231 
   232     }; //class ResidualDijkstra
   233 
   234   public:
   235 
   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())
   246     {
   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) {
   399       for (NodeIt n(_graph); n != INVALID; ++n) {
   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     ///
   416     /// \param s The source node.
   417     /// \param t The target node.
   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;
   428       return *this;
   429     }
   430     
   431     /// @}
   432 
   433     /// \name Execution control
   434 
   435     /// @{
   436 
   437     /// \brief Run the algorithm.
   438     ///
   439     /// 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.
   455     ///
   456     /// \param scaling Enable or disable capacity scaling.
   457     /// If the maximum upper bound and/or the amount of total supply
   458     /// is rather small, the algorithm could be slightly faster without
   459     /// scaling.
   460     ///
   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;
   522     }
   523 
   524     /// @}
   525 
   526     /// \name Query Functions
   527     /// The results of the algorithm can be obtained using these
   528     /// functions.\n
   529     /// The \ref run() function must be called before using them.
   530 
   531     /// @{
   532 
   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.
   546     ///
   547     /// \pre \ref run() must be called before using this function.
   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.
   568     ///
   569     /// \pre \ref run() must be called before using this function.
   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.
   579     ///
   580     /// \pre \ref run() must be called before using this function.
   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.
   592     ///
   593     /// \pre \ref run() must be called before using this function.
   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.
   604     ///
   605     /// \pre \ref run() must be called before using this function.
   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       }
   611     }
   612 
   613     /// @}
   614 
   615   private:
   616 
   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
   699       if (scaling) {
   700         // With scaling
   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];
   709         }
   710         max_sup = std::min(std::min(max_sup, max_dem), max_cap);
   711         _phase_num = 0;
   712         for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
   713           ++_phase_num;
   714       } else {
   715         // Without scaling
   716         _delta = 1;
   717       }
   718 
   719       return OPTIMAL;
   720     }
   721 
   722     ProblemType start() {
   723       // Execute the algorithm
   724       ProblemType pt;
   725       if (_delta > 1)
   726         pt = startWithScaling();
   727       else
   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;
   752       int phase_cnt = 0;
   753       int factor = 4;
   754       ResidualDijkstra _dijkstra(*this);
   755       while (true) {
   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
   772         _excess_nodes.clear();
   773         _deficit_nodes.clear();
   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);
   777         }
   778         int next_node = 0, next_def_node = 0;
   779 
   780         // Find augmenting shortest paths
   781         while (next_node < int(_excess_nodes.size())) {
   782           // Check deficit nodes
   783           if (_delta > 1) {
   784             bool delta_deficit = false;
   785             for ( ; next_def_node < int(_deficit_nodes.size());
   786                     ++next_def_node ) {
   787               if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
   788                 delta_deficit = true;
   789                 break;
   790               }
   791             }
   792             if (!delta_deficit) break;
   793           }
   794 
   795           // Run Dijkstra in the residual network
   796           s = _excess_nodes[next_node];
   797           if ((t = _dijkstra.run(s, _delta)) == -1) {
   798             if (_delta > 1) {
   799               ++next_node;
   800               continue;
   801             }
   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;
   809           if (d > _delta) {
   810             while ((a = _pred[u]) != -1) {
   811               if (_res_cap[a] < d) d = _res_cap[a];
   812               u = _source[a];
   813             }
   814           }
   815           u = t;
   816           while ((a = _pred[u]) != -1) {
   817             _res_cap[a] -= d;
   818             _res_cap[_reverse[a]] += d;
   819             u = _source[a];
   820           }
   821           _excess[s] -= d;
   822           _excess[t] += d;
   823 
   824           if (_excess[s] < _delta) ++next_node;
   825         }
   826 
   827         if (_delta == 1) break;
   828         if (++phase_cnt == _phase_num / 4) factor = 2;
   829         _delta = _delta <= factor ? 1 : _delta / factor;
   830       }
   831 
   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;
   843       int next_node = 0;
   844 
   845       // Find shortest paths
   846       int s, t;
   847       ResidualDijkstra _dijkstra(*this);
   848       while ( _excess[_excess_nodes[next_node]] > 0 ||
   849               ++next_node < int(_excess_nodes.size()) )
   850       {
   851         // Run Dijkstra in the residual network
   852         s = _excess_nodes[next_node];
   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;
   859         if (d > 1) {
   860           while ((a = _pred[u]) != -1) {
   861             if (_res_cap[a] < d) d = _res_cap[a];
   862             u = _source[a];
   863           }
   864         }
   865         u = t;
   866         while ((a = _pred[u]) != -1) {
   867           _res_cap[a] -= d;
   868           _res_cap[_reverse[a]] += d;
   869           u = _source[a];
   870         }
   871         _excess[s] -= d;
   872         _excess[t] += d;
   873       }
   874 
   875       return OPTIMAL;
   876     }
   877 
   878   }; //class CapacityScaling
   879 
   880   ///@}
   881 
   882 } //namespace lemon
   883 
   884 #endif //LEMON_CAPACITY_SCALING_H