lemon/capacity_scaling.h
author Alpar Juttner <alpar@cs.elte.hu>
Fri, 12 Feb 2010 21:53:15 +0100
changeset 826 02109e17027f
parent 821 072ec8120958
child 831 cc9e0c15d747
permissions -rw-r--r--
Merge
     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   /// \brief Default traits class of CapacityScaling algorithm.
    35   ///
    36   /// Default traits class of CapacityScaling algorithm.
    37   /// \tparam GR Digraph type.
    38   /// \tparam V The number type used for flow amounts, capacity bounds
    39   /// and supply values. By default it is \c int.
    40   /// \tparam C The number type used for costs and potentials.
    41   /// By default it is the same as \c V.
    42   template <typename GR, typename V = int, typename C = V>
    43   struct CapacityScalingDefaultTraits
    44   {
    45     /// The type of the digraph
    46     typedef GR Digraph;
    47     /// The type of the flow amounts, capacity bounds and supply values
    48     typedef V Value;
    49     /// The type of the arc costs
    50     typedef C Cost;
    51 
    52     /// \brief The type of the heap used for internal Dijkstra computations.
    53     ///
    54     /// The type of the heap used for internal Dijkstra computations.
    55     /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
    56     /// its priority type must be \c Cost and its cross reference type
    57     /// must be \ref RangeMap "RangeMap<int>".
    58     typedef BinHeap<Cost, RangeMap<int> > Heap;
    59   };
    60 
    61   /// \addtogroup min_cost_flow_algs
    62   /// @{
    63 
    64   /// \brief Implementation of the Capacity Scaling algorithm for
    65   /// finding a \ref min_cost_flow "minimum cost flow".
    66   ///
    67   /// \ref CapacityScaling implements the capacity scaling version
    68   /// of the successive shortest path algorithm for finding a
    69   /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
    70   /// \ref edmondskarp72theoretical. It is an efficient dual
    71   /// solution method.
    72   ///
    73   /// Most of the parameters of the problem (except for the digraph)
    74   /// can be given using separate functions, and the algorithm can be
    75   /// executed using the \ref run() function. If some parameters are not
    76   /// specified, then default values will be used.
    77   ///
    78   /// \tparam GR The digraph type the algorithm runs on.
    79   /// \tparam V The number type used for flow amounts, capacity bounds
    80   /// and supply values in the algorithm. By default, it is \c int.
    81   /// \tparam C The number type used for costs and potentials in the
    82   /// algorithm. By default, it is the same as \c V.
    83   /// \tparam TR The traits class that defines various types used by the
    84   /// algorithm. By default, it is \ref CapacityScalingDefaultTraits
    85   /// "CapacityScalingDefaultTraits<GR, V, C>".
    86   /// In most cases, this parameter should not be set directly,
    87   /// consider to use the named template parameters instead.
    88   ///
    89   /// \warning Both number types must be signed and all input data must
    90   /// be integer.
    91   /// \warning This algorithm does not support negative costs for such
    92   /// arcs that have infinite upper bound.
    93 #ifdef DOXYGEN
    94   template <typename GR, typename V, typename C, typename TR>
    95 #else
    96   template < typename GR, typename V = int, typename C = V,
    97              typename TR = CapacityScalingDefaultTraits<GR, V, C> >
    98 #endif
    99   class CapacityScaling
   100   {
   101   public:
   102 
   103     /// The type of the digraph
   104     typedef typename TR::Digraph Digraph;
   105     /// The type of the flow amounts, capacity bounds and supply values
   106     typedef typename TR::Value Value;
   107     /// The type of the arc costs
   108     typedef typename TR::Cost Cost;
   109 
   110     /// The type of the heap used for internal Dijkstra computations
   111     typedef typename TR::Heap Heap;
   112 
   113     /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
   114     typedef TR Traits;
   115 
   116   public:
   117 
   118     /// \brief Problem type constants for the \c run() function.
   119     ///
   120     /// Enum type containing the problem type constants that can be
   121     /// returned by the \ref run() function of the algorithm.
   122     enum ProblemType {
   123       /// The problem has no feasible solution (flow).
   124       INFEASIBLE,
   125       /// The problem has optimal solution (i.e. it is feasible and
   126       /// bounded), and the algorithm has found optimal flow and node
   127       /// potentials (primal and dual solutions).
   128       OPTIMAL,
   129       /// The digraph contains an arc of negative cost and infinite
   130       /// upper bound. It means that the objective function is unbounded
   131       /// on that arc, however, note that it could actually be bounded
   132       /// over the feasible flows, but this algroithm cannot handle
   133       /// these cases.
   134       UNBOUNDED
   135     };
   136   
   137   private:
   138 
   139     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   140 
   141     typedef std::vector<int> IntVector;
   142     typedef std::vector<char> BoolVector;
   143     typedef std::vector<Value> ValueVector;
   144     typedef std::vector<Cost> CostVector;
   145 
   146   private:
   147 
   148     // Data related to the underlying digraph
   149     const GR &_graph;
   150     int _node_num;
   151     int _arc_num;
   152     int _res_arc_num;
   153     int _root;
   154 
   155     // Parameters of the problem
   156     bool _have_lower;
   157     Value _sum_supply;
   158 
   159     // Data structures for storing the digraph
   160     IntNodeMap _node_id;
   161     IntArcMap _arc_idf;
   162     IntArcMap _arc_idb;
   163     IntVector _first_out;
   164     BoolVector _forward;
   165     IntVector _source;
   166     IntVector _target;
   167     IntVector _reverse;
   168 
   169     // Node and arc data
   170     ValueVector _lower;
   171     ValueVector _upper;
   172     CostVector _cost;
   173     ValueVector _supply;
   174 
   175     ValueVector _res_cap;
   176     CostVector _pi;
   177     ValueVector _excess;
   178     IntVector _excess_nodes;
   179     IntVector _deficit_nodes;
   180 
   181     Value _delta;
   182     int _factor;
   183     IntVector _pred;
   184 
   185   public:
   186   
   187     /// \brief Constant for infinite upper bounds (capacities).
   188     ///
   189     /// Constant for infinite upper bounds (capacities).
   190     /// It is \c std::numeric_limits<Value>::infinity() if available,
   191     /// \c std::numeric_limits<Value>::max() otherwise.
   192     const Value INF;
   193 
   194   private:
   195 
   196     // Special implementation of the Dijkstra algorithm for finding
   197     // shortest paths in the residual network of the digraph with
   198     // respect to the reduced arc costs and modifying the node
   199     // potentials according to the found distance labels.
   200     class ResidualDijkstra
   201     {
   202     private:
   203 
   204       int _node_num;
   205       bool _geq;
   206       const IntVector &_first_out;
   207       const IntVector &_target;
   208       const CostVector &_cost;
   209       const ValueVector &_res_cap;
   210       const ValueVector &_excess;
   211       CostVector &_pi;
   212       IntVector &_pred;
   213       
   214       IntVector _proc_nodes;
   215       CostVector _dist;
   216       
   217     public:
   218 
   219       ResidualDijkstra(CapacityScaling& cs) :
   220         _node_num(cs._node_num), _geq(cs._sum_supply < 0),
   221         _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
   222         _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
   223         _pred(cs._pred), _dist(cs._node_num)
   224       {}
   225 
   226       int run(int s, Value delta = 1) {
   227         RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
   228         Heap heap(heap_cross_ref);
   229         heap.push(s, 0);
   230         _pred[s] = -1;
   231         _proc_nodes.clear();
   232 
   233         // Process nodes
   234         while (!heap.empty() && _excess[heap.top()] > -delta) {
   235           int u = heap.top(), v;
   236           Cost d = heap.prio() + _pi[u], dn;
   237           _dist[u] = heap.prio();
   238           _proc_nodes.push_back(u);
   239           heap.pop();
   240 
   241           // Traverse outgoing residual arcs
   242           int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
   243           for (int a = _first_out[u]; a != last_out; ++a) {
   244             if (_res_cap[a] < delta) continue;
   245             v = _target[a];
   246             switch (heap.state(v)) {
   247               case Heap::PRE_HEAP:
   248                 heap.push(v, d + _cost[a] - _pi[v]);
   249                 _pred[v] = a;
   250                 break;
   251               case Heap::IN_HEAP:
   252                 dn = d + _cost[a] - _pi[v];
   253                 if (dn < heap[v]) {
   254                   heap.decrease(v, dn);
   255                   _pred[v] = a;
   256                 }
   257                 break;
   258               case Heap::POST_HEAP:
   259                 break;
   260             }
   261           }
   262         }
   263         if (heap.empty()) return -1;
   264 
   265         // Update potentials of processed nodes
   266         int t = heap.top();
   267         Cost dt = heap.prio();
   268         for (int i = 0; i < int(_proc_nodes.size()); ++i) {
   269           _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
   270         }
   271 
   272         return t;
   273       }
   274 
   275     }; //class ResidualDijkstra
   276 
   277   public:
   278 
   279     /// \name Named Template Parameters
   280     /// @{
   281 
   282     template <typename T>
   283     struct SetHeapTraits : public Traits {
   284       typedef T Heap;
   285     };
   286 
   287     /// \brief \ref named-templ-param "Named parameter" for setting
   288     /// \c Heap type.
   289     ///
   290     /// \ref named-templ-param "Named parameter" for setting \c Heap
   291     /// type, which is used for internal Dijkstra computations.
   292     /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
   293     /// its priority type must be \c Cost and its cross reference type
   294     /// must be \ref RangeMap "RangeMap<int>".
   295     template <typename T>
   296     struct SetHeap
   297       : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
   298       typedef  CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
   299     };
   300 
   301     /// @}
   302 
   303   public:
   304 
   305     /// \brief Constructor.
   306     ///
   307     /// The constructor of the class.
   308     ///
   309     /// \param graph The digraph the algorithm runs on.
   310     CapacityScaling(const GR& graph) :
   311       _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
   312       INF(std::numeric_limits<Value>::has_infinity ?
   313           std::numeric_limits<Value>::infinity() :
   314           std::numeric_limits<Value>::max())
   315     {
   316       // Check the number types
   317       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
   318         "The flow type of CapacityScaling must be signed");
   319       LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   320         "The cost type of CapacityScaling must be signed");
   321 
   322       // Resize vectors
   323       _node_num = countNodes(_graph);
   324       _arc_num = countArcs(_graph);
   325       _res_arc_num = 2 * (_arc_num + _node_num);
   326       _root = _node_num;
   327       ++_node_num;
   328 
   329       _first_out.resize(_node_num + 1);
   330       _forward.resize(_res_arc_num);
   331       _source.resize(_res_arc_num);
   332       _target.resize(_res_arc_num);
   333       _reverse.resize(_res_arc_num);
   334 
   335       _lower.resize(_res_arc_num);
   336       _upper.resize(_res_arc_num);
   337       _cost.resize(_res_arc_num);
   338       _supply.resize(_node_num);
   339       
   340       _res_cap.resize(_res_arc_num);
   341       _pi.resize(_node_num);
   342       _excess.resize(_node_num);
   343       _pred.resize(_node_num);
   344 
   345       // Copy the graph
   346       int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
   347       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   348         _node_id[n] = i;
   349       }
   350       i = 0;
   351       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   352         _first_out[i] = j;
   353         for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   354           _arc_idf[a] = j;
   355           _forward[j] = true;
   356           _source[j] = i;
   357           _target[j] = _node_id[_graph.runningNode(a)];
   358         }
   359         for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   360           _arc_idb[a] = j;
   361           _forward[j] = false;
   362           _source[j] = i;
   363           _target[j] = _node_id[_graph.runningNode(a)];
   364         }
   365         _forward[j] = false;
   366         _source[j] = i;
   367         _target[j] = _root;
   368         _reverse[j] = k;
   369         _forward[k] = true;
   370         _source[k] = _root;
   371         _target[k] = i;
   372         _reverse[k] = j;
   373         ++j; ++k;
   374       }
   375       _first_out[i] = j;
   376       _first_out[_node_num] = k;
   377       for (ArcIt a(_graph); a != INVALID; ++a) {
   378         int fi = _arc_idf[a];
   379         int bi = _arc_idb[a];
   380         _reverse[fi] = bi;
   381         _reverse[bi] = fi;
   382       }
   383       
   384       // Reset parameters
   385       reset();
   386     }
   387 
   388     /// \name Parameters
   389     /// The parameters of the algorithm can be specified using these
   390     /// functions.
   391 
   392     /// @{
   393 
   394     /// \brief Set the lower bounds on the arcs.
   395     ///
   396     /// This function sets the lower bounds on the arcs.
   397     /// If it is not used before calling \ref run(), the lower bounds
   398     /// will be set to zero on all arcs.
   399     ///
   400     /// \param map An arc map storing the lower bounds.
   401     /// Its \c Value type must be convertible to the \c Value type
   402     /// of the algorithm.
   403     ///
   404     /// \return <tt>(*this)</tt>
   405     template <typename LowerMap>
   406     CapacityScaling& lowerMap(const LowerMap& map) {
   407       _have_lower = true;
   408       for (ArcIt a(_graph); a != INVALID; ++a) {
   409         _lower[_arc_idf[a]] = map[a];
   410         _lower[_arc_idb[a]] = map[a];
   411       }
   412       return *this;
   413     }
   414 
   415     /// \brief Set the upper bounds (capacities) on the arcs.
   416     ///
   417     /// This function sets the upper bounds (capacities) on the arcs.
   418     /// If it is not used before calling \ref run(), the upper bounds
   419     /// will be set to \ref INF on all arcs (i.e. the flow value will be
   420     /// unbounded from above).
   421     ///
   422     /// \param map An arc map storing the upper bounds.
   423     /// Its \c Value type must be convertible to the \c Value type
   424     /// of the algorithm.
   425     ///
   426     /// \return <tt>(*this)</tt>
   427     template<typename UpperMap>
   428     CapacityScaling& upperMap(const UpperMap& map) {
   429       for (ArcIt a(_graph); a != INVALID; ++a) {
   430         _upper[_arc_idf[a]] = map[a];
   431       }
   432       return *this;
   433     }
   434 
   435     /// \brief Set the costs of the arcs.
   436     ///
   437     /// This function sets the costs of the arcs.
   438     /// If it is not used before calling \ref run(), the costs
   439     /// will be set to \c 1 on all arcs.
   440     ///
   441     /// \param map An arc map storing the costs.
   442     /// Its \c Value type must be convertible to the \c Cost type
   443     /// of the algorithm.
   444     ///
   445     /// \return <tt>(*this)</tt>
   446     template<typename CostMap>
   447     CapacityScaling& costMap(const CostMap& map) {
   448       for (ArcIt a(_graph); a != INVALID; ++a) {
   449         _cost[_arc_idf[a]] =  map[a];
   450         _cost[_arc_idb[a]] = -map[a];
   451       }
   452       return *this;
   453     }
   454 
   455     /// \brief Set the supply values of the nodes.
   456     ///
   457     /// This function sets the supply values of the nodes.
   458     /// If neither this function nor \ref stSupply() is used before
   459     /// calling \ref run(), the supply of each node will be set to zero.
   460     ///
   461     /// \param map A node map storing the supply values.
   462     /// Its \c Value type must be convertible to the \c Value type
   463     /// of the algorithm.
   464     ///
   465     /// \return <tt>(*this)</tt>
   466     template<typename SupplyMap>
   467     CapacityScaling& supplyMap(const SupplyMap& map) {
   468       for (NodeIt n(_graph); n != INVALID; ++n) {
   469         _supply[_node_id[n]] = map[n];
   470       }
   471       return *this;
   472     }
   473 
   474     /// \brief Set single source and target nodes and a supply value.
   475     ///
   476     /// This function sets a single source node and a single target node
   477     /// and the required flow value.
   478     /// If neither this function nor \ref supplyMap() is used before
   479     /// calling \ref run(), the supply of each node will be set to zero.
   480     ///
   481     /// Using this function has the same effect as using \ref supplyMap()
   482     /// with such a map in which \c k is assigned to \c s, \c -k is
   483     /// assigned to \c t and all other nodes have zero supply value.
   484     ///
   485     /// \param s The source node.
   486     /// \param t The target node.
   487     /// \param k The required amount of flow from node \c s to node \c t
   488     /// (i.e. the supply of \c s and the demand of \c t).
   489     ///
   490     /// \return <tt>(*this)</tt>
   491     CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
   492       for (int i = 0; i != _node_num; ++i) {
   493         _supply[i] = 0;
   494       }
   495       _supply[_node_id[s]] =  k;
   496       _supply[_node_id[t]] = -k;
   497       return *this;
   498     }
   499     
   500     /// @}
   501 
   502     /// \name Execution control
   503     /// The algorithm can be executed using \ref run().
   504 
   505     /// @{
   506 
   507     /// \brief Run the algorithm.
   508     ///
   509     /// This function runs the algorithm.
   510     /// The paramters can be specified using functions \ref lowerMap(),
   511     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   512     /// For example,
   513     /// \code
   514     ///   CapacityScaling<ListDigraph> cs(graph);
   515     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   516     ///     .supplyMap(sup).run();
   517     /// \endcode
   518     ///
   519     /// This function can be called more than once. All the parameters
   520     /// that have been given are kept for the next call, unless
   521     /// \ref reset() is called, thus only the modified parameters
   522     /// have to be set again. See \ref reset() for examples.
   523     /// However, the underlying digraph must not be modified after this
   524     /// class have been constructed, since it copies and extends the graph.
   525     ///
   526     /// \param factor The capacity scaling factor. It must be larger than
   527     /// one to use scaling. If it is less or equal to one, then scaling
   528     /// will be disabled.
   529     ///
   530     /// \return \c INFEASIBLE if no feasible flow exists,
   531     /// \n \c OPTIMAL if the problem has optimal solution
   532     /// (i.e. it is feasible and bounded), and the algorithm has found
   533     /// optimal flow and node potentials (primal and dual solutions),
   534     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
   535     /// and infinite upper bound. It means that the objective function
   536     /// is unbounded on that arc, however, note that it could actually be
   537     /// bounded over the feasible flows, but this algroithm cannot handle
   538     /// these cases.
   539     ///
   540     /// \see ProblemType
   541     ProblemType run(int factor = 4) {
   542       _factor = factor;
   543       ProblemType pt = init();
   544       if (pt != OPTIMAL) return pt;
   545       return start();
   546     }
   547 
   548     /// \brief Reset all the parameters that have been given before.
   549     ///
   550     /// This function resets all the paramaters that have been given
   551     /// before using functions \ref lowerMap(), \ref upperMap(),
   552     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   553     ///
   554     /// It is useful for multiple run() calls. If this function is not
   555     /// used, all the parameters given before are kept for the next
   556     /// \ref run() call.
   557     /// However, the underlying digraph must not be modified after this
   558     /// class have been constructed, since it copies and extends the graph.
   559     ///
   560     /// For example,
   561     /// \code
   562     ///   CapacityScaling<ListDigraph> cs(graph);
   563     ///
   564     ///   // First run
   565     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   566     ///     .supplyMap(sup).run();
   567     ///
   568     ///   // Run again with modified cost map (reset() is not called,
   569     ///   // so only the cost map have to be set again)
   570     ///   cost[e] += 100;
   571     ///   cs.costMap(cost).run();
   572     ///
   573     ///   // Run again from scratch using reset()
   574     ///   // (the lower bounds will be set to zero on all arcs)
   575     ///   cs.reset();
   576     ///   cs.upperMap(capacity).costMap(cost)
   577     ///     .supplyMap(sup).run();
   578     /// \endcode
   579     ///
   580     /// \return <tt>(*this)</tt>
   581     CapacityScaling& reset() {
   582       for (int i = 0; i != _node_num; ++i) {
   583         _supply[i] = 0;
   584       }
   585       for (int j = 0; j != _res_arc_num; ++j) {
   586         _lower[j] = 0;
   587         _upper[j] = INF;
   588         _cost[j] = _forward[j] ? 1 : -1;
   589       }
   590       _have_lower = false;
   591       return *this;
   592     }
   593 
   594     /// @}
   595 
   596     /// \name Query Functions
   597     /// The results of the algorithm can be obtained using these
   598     /// functions.\n
   599     /// The \ref run() function must be called before using them.
   600 
   601     /// @{
   602 
   603     /// \brief Return the total cost of the found flow.
   604     ///
   605     /// This function returns the total cost of the found flow.
   606     /// Its complexity is O(e).
   607     ///
   608     /// \note The return type of the function can be specified as a
   609     /// template parameter. For example,
   610     /// \code
   611     ///   cs.totalCost<double>();
   612     /// \endcode
   613     /// It is useful if the total cost cannot be stored in the \c Cost
   614     /// type of the algorithm, which is the default return type of the
   615     /// function.
   616     ///
   617     /// \pre \ref run() must be called before using this function.
   618     template <typename Number>
   619     Number totalCost() const {
   620       Number c = 0;
   621       for (ArcIt a(_graph); a != INVALID; ++a) {
   622         int i = _arc_idb[a];
   623         c += static_cast<Number>(_res_cap[i]) *
   624              (-static_cast<Number>(_cost[i]));
   625       }
   626       return c;
   627     }
   628 
   629 #ifndef DOXYGEN
   630     Cost totalCost() const {
   631       return totalCost<Cost>();
   632     }
   633 #endif
   634 
   635     /// \brief Return the flow on the given arc.
   636     ///
   637     /// This function returns the flow on the given arc.
   638     ///
   639     /// \pre \ref run() must be called before using this function.
   640     Value flow(const Arc& a) const {
   641       return _res_cap[_arc_idb[a]];
   642     }
   643 
   644     /// \brief Return the flow map (the primal solution).
   645     ///
   646     /// This function copies the flow value on each arc into the given
   647     /// map. The \c Value type of the algorithm must be convertible to
   648     /// the \c Value type of the map.
   649     ///
   650     /// \pre \ref run() must be called before using this function.
   651     template <typename FlowMap>
   652     void flowMap(FlowMap &map) const {
   653       for (ArcIt a(_graph); a != INVALID; ++a) {
   654         map.set(a, _res_cap[_arc_idb[a]]);
   655       }
   656     }
   657 
   658     /// \brief Return the potential (dual value) of the given node.
   659     ///
   660     /// This function returns the potential (dual value) of the
   661     /// given node.
   662     ///
   663     /// \pre \ref run() must be called before using this function.
   664     Cost potential(const Node& n) const {
   665       return _pi[_node_id[n]];
   666     }
   667 
   668     /// \brief Return the potential map (the dual solution).
   669     ///
   670     /// This function copies the potential (dual value) of each node
   671     /// into the given map.
   672     /// The \c Cost type of the algorithm must be convertible to the
   673     /// \c Value type of the map.
   674     ///
   675     /// \pre \ref run() must be called before using this function.
   676     template <typename PotentialMap>
   677     void potentialMap(PotentialMap &map) const {
   678       for (NodeIt n(_graph); n != INVALID; ++n) {
   679         map.set(n, _pi[_node_id[n]]);
   680       }
   681     }
   682 
   683     /// @}
   684 
   685   private:
   686 
   687     // Initialize the algorithm
   688     ProblemType init() {
   689       if (_node_num <= 1) return INFEASIBLE;
   690 
   691       // Check the sum of supply values
   692       _sum_supply = 0;
   693       for (int i = 0; i != _root; ++i) {
   694         _sum_supply += _supply[i];
   695       }
   696       if (_sum_supply > 0) return INFEASIBLE;
   697       
   698       // Initialize vectors
   699       for (int i = 0; i != _root; ++i) {
   700         _pi[i] = 0;
   701         _excess[i] = _supply[i];
   702       }
   703 
   704       // Remove non-zero lower bounds
   705       const Value MAX = std::numeric_limits<Value>::max();
   706       int last_out;
   707       if (_have_lower) {
   708         for (int i = 0; i != _root; ++i) {
   709           last_out = _first_out[i+1];
   710           for (int j = _first_out[i]; j != last_out; ++j) {
   711             if (_forward[j]) {
   712               Value c = _lower[j];
   713               if (c >= 0) {
   714                 _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
   715               } else {
   716                 _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
   717               }
   718               _excess[i] -= c;
   719               _excess[_target[j]] += c;
   720             } else {
   721               _res_cap[j] = 0;
   722             }
   723           }
   724         }
   725       } else {
   726         for (int j = 0; j != _res_arc_num; ++j) {
   727           _res_cap[j] = _forward[j] ? _upper[j] : 0;
   728         }
   729       }
   730 
   731       // Handle negative costs
   732       for (int i = 0; i != _root; ++i) {
   733         last_out = _first_out[i+1] - 1;
   734         for (int j = _first_out[i]; j != last_out; ++j) {
   735           Value rc = _res_cap[j];
   736           if (_cost[j] < 0 && rc > 0) {
   737             if (rc >= MAX) return UNBOUNDED;
   738             _excess[i] -= rc;
   739             _excess[_target[j]] += rc;
   740             _res_cap[j] = 0;
   741             _res_cap[_reverse[j]] += rc;
   742           }
   743         }
   744       }
   745       
   746       // Handle GEQ supply type
   747       if (_sum_supply < 0) {
   748         _pi[_root] = 0;
   749         _excess[_root] = -_sum_supply;
   750         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   751           int ra = _reverse[a];
   752           _res_cap[a] = -_sum_supply + 1;
   753           _res_cap[ra] = 0;
   754           _cost[a] = 0;
   755           _cost[ra] = 0;
   756         }
   757       } else {
   758         _pi[_root] = 0;
   759         _excess[_root] = 0;
   760         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   761           int ra = _reverse[a];
   762           _res_cap[a] = 1;
   763           _res_cap[ra] = 0;
   764           _cost[a] = 0;
   765           _cost[ra] = 0;
   766         }
   767       }
   768 
   769       // Initialize delta value
   770       if (_factor > 1) {
   771         // With scaling
   772         Value max_sup = 0, max_dem = 0;
   773         for (int i = 0; i != _node_num; ++i) {
   774           Value ex = _excess[i];
   775           if ( ex > max_sup) max_sup =  ex;
   776           if (-ex > max_dem) max_dem = -ex;
   777         }
   778         Value max_cap = 0;
   779         for (int j = 0; j != _res_arc_num; ++j) {
   780           if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
   781         }
   782         max_sup = std::min(std::min(max_sup, max_dem), max_cap);
   783         for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
   784       } else {
   785         // Without scaling
   786         _delta = 1;
   787       }
   788 
   789       return OPTIMAL;
   790     }
   791 
   792     ProblemType start() {
   793       // Execute the algorithm
   794       ProblemType pt;
   795       if (_delta > 1)
   796         pt = startWithScaling();
   797       else
   798         pt = startWithoutScaling();
   799 
   800       // Handle non-zero lower bounds
   801       if (_have_lower) {
   802         int limit = _first_out[_root];
   803         for (int j = 0; j != limit; ++j) {
   804           if (!_forward[j]) _res_cap[j] += _lower[j];
   805         }
   806       }
   807 
   808       // Shift potentials if necessary
   809       Cost pr = _pi[_root];
   810       if (_sum_supply < 0 || pr > 0) {
   811         for (int i = 0; i != _node_num; ++i) {
   812           _pi[i] -= pr;
   813         }        
   814       }
   815       
   816       return pt;
   817     }
   818 
   819     // Execute the capacity scaling algorithm
   820     ProblemType startWithScaling() {
   821       // Perform capacity scaling phases
   822       int s, t;
   823       ResidualDijkstra _dijkstra(*this);
   824       while (true) {
   825         // Saturate all arcs not satisfying the optimality condition
   826         int last_out;
   827         for (int u = 0; u != _node_num; ++u) {
   828           last_out = _sum_supply < 0 ?
   829             _first_out[u+1] : _first_out[u+1] - 1;
   830           for (int a = _first_out[u]; a != last_out; ++a) {
   831             int v = _target[a];
   832             Cost c = _cost[a] + _pi[u] - _pi[v];
   833             Value rc = _res_cap[a];
   834             if (c < 0 && rc >= _delta) {
   835               _excess[u] -= rc;
   836               _excess[v] += rc;
   837               _res_cap[a] = 0;
   838               _res_cap[_reverse[a]] += rc;
   839             }
   840           }
   841         }
   842 
   843         // Find excess nodes and deficit nodes
   844         _excess_nodes.clear();
   845         _deficit_nodes.clear();
   846         for (int u = 0; u != _node_num; ++u) {
   847           Value ex = _excess[u];
   848           if (ex >=  _delta) _excess_nodes.push_back(u);
   849           if (ex <= -_delta) _deficit_nodes.push_back(u);
   850         }
   851         int next_node = 0, next_def_node = 0;
   852 
   853         // Find augmenting shortest paths
   854         while (next_node < int(_excess_nodes.size())) {
   855           // Check deficit nodes
   856           if (_delta > 1) {
   857             bool delta_deficit = false;
   858             for ( ; next_def_node < int(_deficit_nodes.size());
   859                     ++next_def_node ) {
   860               if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
   861                 delta_deficit = true;
   862                 break;
   863               }
   864             }
   865             if (!delta_deficit) break;
   866           }
   867 
   868           // Run Dijkstra in the residual network
   869           s = _excess_nodes[next_node];
   870           if ((t = _dijkstra.run(s, _delta)) == -1) {
   871             if (_delta > 1) {
   872               ++next_node;
   873               continue;
   874             }
   875             return INFEASIBLE;
   876           }
   877 
   878           // Augment along a shortest path from s to t
   879           Value d = std::min(_excess[s], -_excess[t]);
   880           int u = t;
   881           int a;
   882           if (d > _delta) {
   883             while ((a = _pred[u]) != -1) {
   884               if (_res_cap[a] < d) d = _res_cap[a];
   885               u = _source[a];
   886             }
   887           }
   888           u = t;
   889           while ((a = _pred[u]) != -1) {
   890             _res_cap[a] -= d;
   891             _res_cap[_reverse[a]] += d;
   892             u = _source[a];
   893           }
   894           _excess[s] -= d;
   895           _excess[t] += d;
   896 
   897           if (_excess[s] < _delta) ++next_node;
   898         }
   899 
   900         if (_delta == 1) break;
   901         _delta = _delta <= _factor ? 1 : _delta / _factor;
   902       }
   903 
   904       return OPTIMAL;
   905     }
   906 
   907     // Execute the successive shortest path algorithm
   908     ProblemType startWithoutScaling() {
   909       // Find excess nodes
   910       _excess_nodes.clear();
   911       for (int i = 0; i != _node_num; ++i) {
   912         if (_excess[i] > 0) _excess_nodes.push_back(i);
   913       }
   914       if (_excess_nodes.size() == 0) return OPTIMAL;
   915       int next_node = 0;
   916 
   917       // Find shortest paths
   918       int s, t;
   919       ResidualDijkstra _dijkstra(*this);
   920       while ( _excess[_excess_nodes[next_node]] > 0 ||
   921               ++next_node < int(_excess_nodes.size()) )
   922       {
   923         // Run Dijkstra in the residual network
   924         s = _excess_nodes[next_node];
   925         if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
   926 
   927         // Augment along a shortest path from s to t
   928         Value d = std::min(_excess[s], -_excess[t]);
   929         int u = t;
   930         int a;
   931         if (d > 1) {
   932           while ((a = _pred[u]) != -1) {
   933             if (_res_cap[a] < d) d = _res_cap[a];
   934             u = _source[a];
   935           }
   936         }
   937         u = t;
   938         while ((a = _pred[u]) != -1) {
   939           _res_cap[a] -= d;
   940           _res_cap[_reverse[a]] += d;
   941           u = _source[a];
   942         }
   943         _excess[s] -= d;
   944         _excess[t] += d;
   945       }
   946 
   947       return OPTIMAL;
   948     }
   949 
   950   }; //class CapacityScaling
   951 
   952   ///@}
   953 
   954 } //namespace lemon
   955 
   956 #endif //LEMON_CAPACITY_SCALING_H