lemon/capacity_scaling.h
author Alpar Juttner <alpar@cs.elte.hu>
Sun, 14 Feb 2010 19:23:55 +0100
changeset 842 c2ff0a365245
parent 830 75c97c3786d6
parent 825 75e6020b19b1
child 840 2914b6f0fde0
permissions -rw-r--r--
ArgParser can throw exception instead of exit(1) (#332)
     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       // Reset data structures
   323       reset();
   324     }
   325 
   326     /// \name Parameters
   327     /// The parameters of the algorithm can be specified using these
   328     /// functions.
   329 
   330     /// @{
   331 
   332     /// \brief Set the lower bounds on the arcs.
   333     ///
   334     /// This function sets the lower bounds on the arcs.
   335     /// If it is not used before calling \ref run(), the lower bounds
   336     /// will be set to zero on all arcs.
   337     ///
   338     /// \param map An arc map storing the lower bounds.
   339     /// Its \c Value type must be convertible to the \c Value type
   340     /// of the algorithm.
   341     ///
   342     /// \return <tt>(*this)</tt>
   343     template <typename LowerMap>
   344     CapacityScaling& lowerMap(const LowerMap& map) {
   345       _have_lower = true;
   346       for (ArcIt a(_graph); a != INVALID; ++a) {
   347         _lower[_arc_idf[a]] = map[a];
   348         _lower[_arc_idb[a]] = map[a];
   349       }
   350       return *this;
   351     }
   352 
   353     /// \brief Set the upper bounds (capacities) on the arcs.
   354     ///
   355     /// This function sets the upper bounds (capacities) on the arcs.
   356     /// If it is not used before calling \ref run(), the upper bounds
   357     /// will be set to \ref INF on all arcs (i.e. the flow value will be
   358     /// unbounded from above).
   359     ///
   360     /// \param map An arc map storing the upper bounds.
   361     /// Its \c Value type must be convertible to the \c Value type
   362     /// of the algorithm.
   363     ///
   364     /// \return <tt>(*this)</tt>
   365     template<typename UpperMap>
   366     CapacityScaling& upperMap(const UpperMap& map) {
   367       for (ArcIt a(_graph); a != INVALID; ++a) {
   368         _upper[_arc_idf[a]] = map[a];
   369       }
   370       return *this;
   371     }
   372 
   373     /// \brief Set the costs of the arcs.
   374     ///
   375     /// This function sets the costs of the arcs.
   376     /// If it is not used before calling \ref run(), the costs
   377     /// will be set to \c 1 on all arcs.
   378     ///
   379     /// \param map An arc map storing the costs.
   380     /// Its \c Value type must be convertible to the \c Cost type
   381     /// of the algorithm.
   382     ///
   383     /// \return <tt>(*this)</tt>
   384     template<typename CostMap>
   385     CapacityScaling& costMap(const CostMap& map) {
   386       for (ArcIt a(_graph); a != INVALID; ++a) {
   387         _cost[_arc_idf[a]] =  map[a];
   388         _cost[_arc_idb[a]] = -map[a];
   389       }
   390       return *this;
   391     }
   392 
   393     /// \brief Set the supply values of the nodes.
   394     ///
   395     /// This function sets the supply values of the nodes.
   396     /// If neither this function nor \ref stSupply() is used before
   397     /// calling \ref run(), the supply of each node will be set to zero.
   398     ///
   399     /// \param map A node map storing the supply values.
   400     /// Its \c Value type must be convertible to the \c Value type
   401     /// of the algorithm.
   402     ///
   403     /// \return <tt>(*this)</tt>
   404     template<typename SupplyMap>
   405     CapacityScaling& supplyMap(const SupplyMap& map) {
   406       for (NodeIt n(_graph); n != INVALID; ++n) {
   407         _supply[_node_id[n]] = map[n];
   408       }
   409       return *this;
   410     }
   411 
   412     /// \brief Set single source and target nodes and a supply value.
   413     ///
   414     /// This function sets a single source node and a single target node
   415     /// and the required flow value.
   416     /// If neither this function nor \ref supplyMap() is used before
   417     /// calling \ref run(), the supply of each node will be set to zero.
   418     ///
   419     /// Using this function has the same effect as using \ref supplyMap()
   420     /// with such a map in which \c k is assigned to \c s, \c -k is
   421     /// assigned to \c t and all other nodes have zero supply value.
   422     ///
   423     /// \param s The source node.
   424     /// \param t The target node.
   425     /// \param k The required amount of flow from node \c s to node \c t
   426     /// (i.e. the supply of \c s and the demand of \c t).
   427     ///
   428     /// \return <tt>(*this)</tt>
   429     CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
   430       for (int i = 0; i != _node_num; ++i) {
   431         _supply[i] = 0;
   432       }
   433       _supply[_node_id[s]] =  k;
   434       _supply[_node_id[t]] = -k;
   435       return *this;
   436     }
   437     
   438     /// @}
   439 
   440     /// \name Execution control
   441     /// The algorithm can be executed using \ref run().
   442 
   443     /// @{
   444 
   445     /// \brief Run the algorithm.
   446     ///
   447     /// This function runs the algorithm.
   448     /// The paramters can be specified using functions \ref lowerMap(),
   449     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   450     /// For example,
   451     /// \code
   452     ///   CapacityScaling<ListDigraph> cs(graph);
   453     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   454     ///     .supplyMap(sup).run();
   455     /// \endcode
   456     ///
   457     /// This function can be called more than once. All the given parameters
   458     /// are kept for the next call, unless \ref resetParams() or \ref reset()
   459     /// is used, thus only the modified parameters have to be set again.
   460     /// If the underlying digraph was also modified after the construction
   461     /// of the class (or the last \ref reset() call), then the \ref reset()
   462     /// function must be called.
   463     ///
   464     /// \param factor The capacity scaling factor. It must be larger than
   465     /// one to use scaling. If it is less or equal to one, then scaling
   466     /// will be disabled.
   467     ///
   468     /// \return \c INFEASIBLE if no feasible flow exists,
   469     /// \n \c OPTIMAL if the problem has optimal solution
   470     /// (i.e. it is feasible and bounded), and the algorithm has found
   471     /// optimal flow and node potentials (primal and dual solutions),
   472     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
   473     /// and infinite upper bound. It means that the objective function
   474     /// is unbounded on that arc, however, note that it could actually be
   475     /// bounded over the feasible flows, but this algroithm cannot handle
   476     /// these cases.
   477     ///
   478     /// \see ProblemType
   479     /// \see resetParams(), reset()
   480     ProblemType run(int factor = 4) {
   481       _factor = factor;
   482       ProblemType pt = init();
   483       if (pt != OPTIMAL) return pt;
   484       return start();
   485     }
   486 
   487     /// \brief Reset all the parameters that have been given before.
   488     ///
   489     /// This function resets all the paramaters that have been given
   490     /// before using functions \ref lowerMap(), \ref upperMap(),
   491     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   492     ///
   493     /// It is useful for multiple \ref run() calls. Basically, all the given
   494     /// parameters are kept for the next \ref run() call, unless
   495     /// \ref resetParams() or \ref reset() is used.
   496     /// If the underlying digraph was also modified after the construction
   497     /// of the class or the last \ref reset() call, then the \ref reset()
   498     /// function must be used, otherwise \ref resetParams() is sufficient.
   499     ///
   500     /// For example,
   501     /// \code
   502     ///   CapacityScaling<ListDigraph> cs(graph);
   503     ///
   504     ///   // First run
   505     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   506     ///     .supplyMap(sup).run();
   507     ///
   508     ///   // Run again with modified cost map (resetParams() is not called,
   509     ///   // so only the cost map have to be set again)
   510     ///   cost[e] += 100;
   511     ///   cs.costMap(cost).run();
   512     ///
   513     ///   // Run again from scratch using resetParams()
   514     ///   // (the lower bounds will be set to zero on all arcs)
   515     ///   cs.resetParams();
   516     ///   cs.upperMap(capacity).costMap(cost)
   517     ///     .supplyMap(sup).run();
   518     /// \endcode
   519     ///
   520     /// \return <tt>(*this)</tt>
   521     ///
   522     /// \see reset(), run()
   523     CapacityScaling& resetParams() {
   524       for (int i = 0; i != _node_num; ++i) {
   525         _supply[i] = 0;
   526       }
   527       for (int j = 0; j != _res_arc_num; ++j) {
   528         _lower[j] = 0;
   529         _upper[j] = INF;
   530         _cost[j] = _forward[j] ? 1 : -1;
   531       }
   532       _have_lower = false;
   533       return *this;
   534     }
   535 
   536     /// \brief Reset the internal data structures and all the parameters
   537     /// that have been given before.
   538     ///
   539     /// This function resets the internal data structures and all the
   540     /// paramaters that have been given before using functions \ref lowerMap(),
   541     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   542     ///
   543     /// It is useful for multiple \ref run() calls. Basically, all the given
   544     /// parameters are kept for the next \ref run() call, unless
   545     /// \ref resetParams() or \ref reset() is used.
   546     /// If the underlying digraph was also modified after the construction
   547     /// of the class or the last \ref reset() call, then the \ref reset()
   548     /// function must be used, otherwise \ref resetParams() is sufficient.
   549     ///
   550     /// See \ref resetParams() for examples.
   551     ///
   552     /// \return <tt>(*this)</tt>
   553     ///
   554     /// \see resetParams(), run()
   555     CapacityScaling& reset() {
   556       // Resize vectors
   557       _node_num = countNodes(_graph);
   558       _arc_num = countArcs(_graph);
   559       _res_arc_num = 2 * (_arc_num + _node_num);
   560       _root = _node_num;
   561       ++_node_num;
   562 
   563       _first_out.resize(_node_num + 1);
   564       _forward.resize(_res_arc_num);
   565       _source.resize(_res_arc_num);
   566       _target.resize(_res_arc_num);
   567       _reverse.resize(_res_arc_num);
   568 
   569       _lower.resize(_res_arc_num);
   570       _upper.resize(_res_arc_num);
   571       _cost.resize(_res_arc_num);
   572       _supply.resize(_node_num);
   573       
   574       _res_cap.resize(_res_arc_num);
   575       _pi.resize(_node_num);
   576       _excess.resize(_node_num);
   577       _pred.resize(_node_num);
   578 
   579       // Copy the graph
   580       int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
   581       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   582         _node_id[n] = i;
   583       }
   584       i = 0;
   585       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   586         _first_out[i] = j;
   587         for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   588           _arc_idf[a] = j;
   589           _forward[j] = true;
   590           _source[j] = i;
   591           _target[j] = _node_id[_graph.runningNode(a)];
   592         }
   593         for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   594           _arc_idb[a] = j;
   595           _forward[j] = false;
   596           _source[j] = i;
   597           _target[j] = _node_id[_graph.runningNode(a)];
   598         }
   599         _forward[j] = false;
   600         _source[j] = i;
   601         _target[j] = _root;
   602         _reverse[j] = k;
   603         _forward[k] = true;
   604         _source[k] = _root;
   605         _target[k] = i;
   606         _reverse[k] = j;
   607         ++j; ++k;
   608       }
   609       _first_out[i] = j;
   610       _first_out[_node_num] = k;
   611       for (ArcIt a(_graph); a != INVALID; ++a) {
   612         int fi = _arc_idf[a];
   613         int bi = _arc_idb[a];
   614         _reverse[fi] = bi;
   615         _reverse[bi] = fi;
   616       }
   617       
   618       // Reset parameters
   619       resetParams();
   620       return *this;
   621     }
   622 
   623     /// @}
   624 
   625     /// \name Query Functions
   626     /// The results of the algorithm can be obtained using these
   627     /// functions.\n
   628     /// The \ref run() function must be called before using them.
   629 
   630     /// @{
   631 
   632     /// \brief Return the total cost of the found flow.
   633     ///
   634     /// This function returns the total cost of the found flow.
   635     /// Its complexity is O(e).
   636     ///
   637     /// \note The return type of the function can be specified as a
   638     /// template parameter. For example,
   639     /// \code
   640     ///   cs.totalCost<double>();
   641     /// \endcode
   642     /// It is useful if the total cost cannot be stored in the \c Cost
   643     /// type of the algorithm, which is the default return type of the
   644     /// function.
   645     ///
   646     /// \pre \ref run() must be called before using this function.
   647     template <typename Number>
   648     Number totalCost() const {
   649       Number c = 0;
   650       for (ArcIt a(_graph); a != INVALID; ++a) {
   651         int i = _arc_idb[a];
   652         c += static_cast<Number>(_res_cap[i]) *
   653              (-static_cast<Number>(_cost[i]));
   654       }
   655       return c;
   656     }
   657 
   658 #ifndef DOXYGEN
   659     Cost totalCost() const {
   660       return totalCost<Cost>();
   661     }
   662 #endif
   663 
   664     /// \brief Return the flow on the given arc.
   665     ///
   666     /// This function returns the flow on the given arc.
   667     ///
   668     /// \pre \ref run() must be called before using this function.
   669     Value flow(const Arc& a) const {
   670       return _res_cap[_arc_idb[a]];
   671     }
   672 
   673     /// \brief Return the flow map (the primal solution).
   674     ///
   675     /// This function copies the flow value on each arc into the given
   676     /// map. The \c Value type of the algorithm must be convertible to
   677     /// the \c Value type of the map.
   678     ///
   679     /// \pre \ref run() must be called before using this function.
   680     template <typename FlowMap>
   681     void flowMap(FlowMap &map) const {
   682       for (ArcIt a(_graph); a != INVALID; ++a) {
   683         map.set(a, _res_cap[_arc_idb[a]]);
   684       }
   685     }
   686 
   687     /// \brief Return the potential (dual value) of the given node.
   688     ///
   689     /// This function returns the potential (dual value) of the
   690     /// given node.
   691     ///
   692     /// \pre \ref run() must be called before using this function.
   693     Cost potential(const Node& n) const {
   694       return _pi[_node_id[n]];
   695     }
   696 
   697     /// \brief Return the potential map (the dual solution).
   698     ///
   699     /// This function copies the potential (dual value) of each node
   700     /// into the given map.
   701     /// The \c Cost type of the algorithm must be convertible to the
   702     /// \c Value type of the map.
   703     ///
   704     /// \pre \ref run() must be called before using this function.
   705     template <typename PotentialMap>
   706     void potentialMap(PotentialMap &map) const {
   707       for (NodeIt n(_graph); n != INVALID; ++n) {
   708         map.set(n, _pi[_node_id[n]]);
   709       }
   710     }
   711 
   712     /// @}
   713 
   714   private:
   715 
   716     // Initialize the algorithm
   717     ProblemType init() {
   718       if (_node_num <= 1) return INFEASIBLE;
   719 
   720       // Check the sum of supply values
   721       _sum_supply = 0;
   722       for (int i = 0; i != _root; ++i) {
   723         _sum_supply += _supply[i];
   724       }
   725       if (_sum_supply > 0) return INFEASIBLE;
   726       
   727       // Initialize vectors
   728       for (int i = 0; i != _root; ++i) {
   729         _pi[i] = 0;
   730         _excess[i] = _supply[i];
   731       }
   732 
   733       // Remove non-zero lower bounds
   734       const Value MAX = std::numeric_limits<Value>::max();
   735       int last_out;
   736       if (_have_lower) {
   737         for (int i = 0; i != _root; ++i) {
   738           last_out = _first_out[i+1];
   739           for (int j = _first_out[i]; j != last_out; ++j) {
   740             if (_forward[j]) {
   741               Value c = _lower[j];
   742               if (c >= 0) {
   743                 _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
   744               } else {
   745                 _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
   746               }
   747               _excess[i] -= c;
   748               _excess[_target[j]] += c;
   749             } else {
   750               _res_cap[j] = 0;
   751             }
   752           }
   753         }
   754       } else {
   755         for (int j = 0; j != _res_arc_num; ++j) {
   756           _res_cap[j] = _forward[j] ? _upper[j] : 0;
   757         }
   758       }
   759 
   760       // Handle negative costs
   761       for (int i = 0; i != _root; ++i) {
   762         last_out = _first_out[i+1] - 1;
   763         for (int j = _first_out[i]; j != last_out; ++j) {
   764           Value rc = _res_cap[j];
   765           if (_cost[j] < 0 && rc > 0) {
   766             if (rc >= MAX) return UNBOUNDED;
   767             _excess[i] -= rc;
   768             _excess[_target[j]] += rc;
   769             _res_cap[j] = 0;
   770             _res_cap[_reverse[j]] += rc;
   771           }
   772         }
   773       }
   774       
   775       // Handle GEQ supply type
   776       if (_sum_supply < 0) {
   777         _pi[_root] = 0;
   778         _excess[_root] = -_sum_supply;
   779         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   780           int ra = _reverse[a];
   781           _res_cap[a] = -_sum_supply + 1;
   782           _res_cap[ra] = 0;
   783           _cost[a] = 0;
   784           _cost[ra] = 0;
   785         }
   786       } else {
   787         _pi[_root] = 0;
   788         _excess[_root] = 0;
   789         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   790           int ra = _reverse[a];
   791           _res_cap[a] = 1;
   792           _res_cap[ra] = 0;
   793           _cost[a] = 0;
   794           _cost[ra] = 0;
   795         }
   796       }
   797 
   798       // Initialize delta value
   799       if (_factor > 1) {
   800         // With scaling
   801         Value max_sup = 0, max_dem = 0;
   802         for (int i = 0; i != _node_num; ++i) {
   803           Value ex = _excess[i];
   804           if ( ex > max_sup) max_sup =  ex;
   805           if (-ex > max_dem) max_dem = -ex;
   806         }
   807         Value max_cap = 0;
   808         for (int j = 0; j != _res_arc_num; ++j) {
   809           if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
   810         }
   811         max_sup = std::min(std::min(max_sup, max_dem), max_cap);
   812         for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
   813       } else {
   814         // Without scaling
   815         _delta = 1;
   816       }
   817 
   818       return OPTIMAL;
   819     }
   820 
   821     ProblemType start() {
   822       // Execute the algorithm
   823       ProblemType pt;
   824       if (_delta > 1)
   825         pt = startWithScaling();
   826       else
   827         pt = startWithoutScaling();
   828 
   829       // Handle non-zero lower bounds
   830       if (_have_lower) {
   831         int limit = _first_out[_root];
   832         for (int j = 0; j != limit; ++j) {
   833           if (!_forward[j]) _res_cap[j] += _lower[j];
   834         }
   835       }
   836 
   837       // Shift potentials if necessary
   838       Cost pr = _pi[_root];
   839       if (_sum_supply < 0 || pr > 0) {
   840         for (int i = 0; i != _node_num; ++i) {
   841           _pi[i] -= pr;
   842         }        
   843       }
   844       
   845       return pt;
   846     }
   847 
   848     // Execute the capacity scaling algorithm
   849     ProblemType startWithScaling() {
   850       // Perform capacity scaling phases
   851       int s, t;
   852       ResidualDijkstra _dijkstra(*this);
   853       while (true) {
   854         // Saturate all arcs not satisfying the optimality condition
   855         int last_out;
   856         for (int u = 0; u != _node_num; ++u) {
   857           last_out = _sum_supply < 0 ?
   858             _first_out[u+1] : _first_out[u+1] - 1;
   859           for (int a = _first_out[u]; a != last_out; ++a) {
   860             int v = _target[a];
   861             Cost c = _cost[a] + _pi[u] - _pi[v];
   862             Value rc = _res_cap[a];
   863             if (c < 0 && rc >= _delta) {
   864               _excess[u] -= rc;
   865               _excess[v] += rc;
   866               _res_cap[a] = 0;
   867               _res_cap[_reverse[a]] += rc;
   868             }
   869           }
   870         }
   871 
   872         // Find excess nodes and deficit nodes
   873         _excess_nodes.clear();
   874         _deficit_nodes.clear();
   875         for (int u = 0; u != _node_num; ++u) {
   876           Value ex = _excess[u];
   877           if (ex >=  _delta) _excess_nodes.push_back(u);
   878           if (ex <= -_delta) _deficit_nodes.push_back(u);
   879         }
   880         int next_node = 0, next_def_node = 0;
   881 
   882         // Find augmenting shortest paths
   883         while (next_node < int(_excess_nodes.size())) {
   884           // Check deficit nodes
   885           if (_delta > 1) {
   886             bool delta_deficit = false;
   887             for ( ; next_def_node < int(_deficit_nodes.size());
   888                     ++next_def_node ) {
   889               if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
   890                 delta_deficit = true;
   891                 break;
   892               }
   893             }
   894             if (!delta_deficit) break;
   895           }
   896 
   897           // Run Dijkstra in the residual network
   898           s = _excess_nodes[next_node];
   899           if ((t = _dijkstra.run(s, _delta)) == -1) {
   900             if (_delta > 1) {
   901               ++next_node;
   902               continue;
   903             }
   904             return INFEASIBLE;
   905           }
   906 
   907           // Augment along a shortest path from s to t
   908           Value d = std::min(_excess[s], -_excess[t]);
   909           int u = t;
   910           int a;
   911           if (d > _delta) {
   912             while ((a = _pred[u]) != -1) {
   913               if (_res_cap[a] < d) d = _res_cap[a];
   914               u = _source[a];
   915             }
   916           }
   917           u = t;
   918           while ((a = _pred[u]) != -1) {
   919             _res_cap[a] -= d;
   920             _res_cap[_reverse[a]] += d;
   921             u = _source[a];
   922           }
   923           _excess[s] -= d;
   924           _excess[t] += d;
   925 
   926           if (_excess[s] < _delta) ++next_node;
   927         }
   928 
   929         if (_delta == 1) break;
   930         _delta = _delta <= _factor ? 1 : _delta / _factor;
   931       }
   932 
   933       return OPTIMAL;
   934     }
   935 
   936     // Execute the successive shortest path algorithm
   937     ProblemType startWithoutScaling() {
   938       // Find excess nodes
   939       _excess_nodes.clear();
   940       for (int i = 0; i != _node_num; ++i) {
   941         if (_excess[i] > 0) _excess_nodes.push_back(i);
   942       }
   943       if (_excess_nodes.size() == 0) return OPTIMAL;
   944       int next_node = 0;
   945 
   946       // Find shortest paths
   947       int s, t;
   948       ResidualDijkstra _dijkstra(*this);
   949       while ( _excess[_excess_nodes[next_node]] > 0 ||
   950               ++next_node < int(_excess_nodes.size()) )
   951       {
   952         // Run Dijkstra in the residual network
   953         s = _excess_nodes[next_node];
   954         if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
   955 
   956         // Augment along a shortest path from s to t
   957         Value d = std::min(_excess[s], -_excess[t]);
   958         int u = t;
   959         int a;
   960         if (d > 1) {
   961           while ((a = _pred[u]) != -1) {
   962             if (_res_cap[a] < d) d = _res_cap[a];
   963             u = _source[a];
   964           }
   965         }
   966         u = t;
   967         while ((a = _pred[u]) != -1) {
   968           _res_cap[a] -= d;
   969           _res_cap[_reverse[a]] += d;
   970           u = _source[a];
   971         }
   972         _excess[s] -= d;
   973         _excess[t] += d;
   974       }
   975 
   976       return OPTIMAL;
   977     }
   978 
   979   }; //class CapacityScaling
   980 
   981   ///@}
   982 
   983 } //namespace lemon
   984 
   985 #endif //LEMON_CAPACITY_SCALING_H