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