lemon/capacity_scaling.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 13 Nov 2009 00:10:33 +0100
changeset 815 aef153f430e1
parent 812 4b1b378823dc
child 821 072ec8120958
permissions -rw-r--r--
Entirely rework cycle canceling algorithms (#180)

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