lemon/cost_scaling.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 19 Feb 2010 14:08:32 +0100
changeset 844 a6eb9698c321
parent 830 75c97c3786d6
parent 825 75e6020b19b1
child 840 2914b6f0fde0
permissions -rw-r--r--
Support tolerance technique for BellmanFord (#51)

A new operation traits class BellmanFordToleranceOperationTraits
is introduced, which uses the tolerance technique in its less()
function. This class can be used with the SetOperationTraits
named template parameter.
     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_COST_SCALING_H
    20 #define LEMON_COST_SCALING_H
    21 
    22 /// \ingroup min_cost_flow_algs
    23 /// \file
    24 /// \brief Cost scaling algorithm for finding a minimum cost flow.
    25 
    26 #include <vector>
    27 #include <deque>
    28 #include <limits>
    29 
    30 #include <lemon/core.h>
    31 #include <lemon/maps.h>
    32 #include <lemon/math.h>
    33 #include <lemon/static_graph.h>
    34 #include <lemon/circulation.h>
    35 #include <lemon/bellman_ford.h>
    36 
    37 namespace lemon {
    38 
    39   /// \brief Default traits class of CostScaling algorithm.
    40   ///
    41   /// Default traits class of CostScaling algorithm.
    42   /// \tparam GR Digraph type.
    43   /// \tparam V The number type used for flow amounts, capacity bounds
    44   /// and supply values. By default it is \c int.
    45   /// \tparam C The number type used for costs and potentials.
    46   /// By default it is the same as \c V.
    47 #ifdef DOXYGEN
    48   template <typename GR, typename V = int, typename C = V>
    49 #else
    50   template < typename GR, typename V = int, typename C = V,
    51              bool integer = std::numeric_limits<C>::is_integer >
    52 #endif
    53   struct CostScalingDefaultTraits
    54   {
    55     /// The type of the digraph
    56     typedef GR Digraph;
    57     /// The type of the flow amounts, capacity bounds and supply values
    58     typedef V Value;
    59     /// The type of the arc costs
    60     typedef C Cost;
    61 
    62     /// \brief The large cost type used for internal computations
    63     ///
    64     /// The large cost type used for internal computations.
    65     /// It is \c long \c long if the \c Cost type is integer,
    66     /// otherwise it is \c double.
    67     /// \c Cost must be convertible to \c LargeCost.
    68     typedef double LargeCost;
    69   };
    70 
    71   // Default traits class for integer cost types
    72   template <typename GR, typename V, typename C>
    73   struct CostScalingDefaultTraits<GR, V, C, true>
    74   {
    75     typedef GR Digraph;
    76     typedef V Value;
    77     typedef C Cost;
    78 #ifdef LEMON_HAVE_LONG_LONG
    79     typedef long long LargeCost;
    80 #else
    81     typedef long LargeCost;
    82 #endif
    83   };
    84 
    85 
    86   /// \addtogroup min_cost_flow_algs
    87   /// @{
    88 
    89   /// \brief Implementation of the Cost Scaling algorithm for
    90   /// finding a \ref min_cost_flow "minimum cost flow".
    91   ///
    92   /// \ref CostScaling implements a cost scaling algorithm that performs
    93   /// push/augment and relabel operations for finding a \ref min_cost_flow
    94   /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
    95   /// \ref goldberg97efficient, \ref bunnagel98efficient. 
    96   /// It is a highly efficient primal-dual solution method, which
    97   /// can be viewed as the generalization of the \ref Preflow
    98   /// "preflow push-relabel" algorithm for the maximum flow problem.
    99   ///
   100   /// Most of the parameters of the problem (except for the digraph)
   101   /// can be given using separate functions, and the algorithm can be
   102   /// executed using the \ref run() function. If some parameters are not
   103   /// specified, then default values will be used.
   104   ///
   105   /// \tparam GR The digraph type the algorithm runs on.
   106   /// \tparam V The number type used for flow amounts, capacity bounds
   107   /// and supply values in the algorithm. By default, it is \c int.
   108   /// \tparam C The number type used for costs and potentials in the
   109   /// algorithm. By default, it is the same as \c V.
   110   /// \tparam TR The traits class that defines various types used by the
   111   /// algorithm. By default, it is \ref CostScalingDefaultTraits
   112   /// "CostScalingDefaultTraits<GR, V, C>".
   113   /// In most cases, this parameter should not be set directly,
   114   /// consider to use the named template parameters instead.
   115   ///
   116   /// \warning Both number types must be signed and all input data must
   117   /// be integer.
   118   /// \warning This algorithm does not support negative costs for such
   119   /// arcs that have infinite upper bound.
   120   ///
   121   /// \note %CostScaling provides three different internal methods,
   122   /// from which the most efficient one is used by default.
   123   /// For more information, see \ref Method.
   124 #ifdef DOXYGEN
   125   template <typename GR, typename V, typename C, typename TR>
   126 #else
   127   template < typename GR, typename V = int, typename C = V,
   128              typename TR = CostScalingDefaultTraits<GR, V, C> >
   129 #endif
   130   class CostScaling
   131   {
   132   public:
   133 
   134     /// The type of the digraph
   135     typedef typename TR::Digraph Digraph;
   136     /// The type of the flow amounts, capacity bounds and supply values
   137     typedef typename TR::Value Value;
   138     /// The type of the arc costs
   139     typedef typename TR::Cost Cost;
   140 
   141     /// \brief The large cost type
   142     ///
   143     /// The large cost type used for internal computations.
   144     /// By default, it is \c long \c long if the \c Cost type is integer,
   145     /// otherwise it is \c double.
   146     typedef typename TR::LargeCost LargeCost;
   147 
   148     /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
   149     typedef TR Traits;
   150 
   151   public:
   152 
   153     /// \brief Problem type constants for the \c run() function.
   154     ///
   155     /// Enum type containing the problem type constants that can be
   156     /// returned by the \ref run() function of the algorithm.
   157     enum ProblemType {
   158       /// The problem has no feasible solution (flow).
   159       INFEASIBLE,
   160       /// The problem has optimal solution (i.e. it is feasible and
   161       /// bounded), and the algorithm has found optimal flow and node
   162       /// potentials (primal and dual solutions).
   163       OPTIMAL,
   164       /// The digraph contains an arc of negative cost and infinite
   165       /// upper bound. It means that the objective function is unbounded
   166       /// on that arc, however, note that it could actually be bounded
   167       /// over the feasible flows, but this algroithm cannot handle
   168       /// these cases.
   169       UNBOUNDED
   170     };
   171 
   172     /// \brief Constants for selecting the internal method.
   173     ///
   174     /// Enum type containing constants for selecting the internal method
   175     /// for the \ref run() function.
   176     ///
   177     /// \ref CostScaling provides three internal methods that differ mainly
   178     /// in their base operations, which are used in conjunction with the
   179     /// relabel operation.
   180     /// By default, the so called \ref PARTIAL_AUGMENT
   181     /// "Partial Augment-Relabel" method is used, which proved to be
   182     /// the most efficient and the most robust on various test inputs.
   183     /// However, the other methods can be selected using the \ref run()
   184     /// function with the proper parameter.
   185     enum Method {
   186       /// Local push operations are used, i.e. flow is moved only on one
   187       /// admissible arc at once.
   188       PUSH,
   189       /// Augment operations are used, i.e. flow is moved on admissible
   190       /// paths from a node with excess to a node with deficit.
   191       AUGMENT,
   192       /// Partial augment operations are used, i.e. flow is moved on 
   193       /// admissible paths started from a node with excess, but the
   194       /// lengths of these paths are limited. This method can be viewed
   195       /// as a combined version of the previous two operations.
   196       PARTIAL_AUGMENT
   197     };
   198 
   199   private:
   200 
   201     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   202 
   203     typedef std::vector<int> IntVector;
   204     typedef std::vector<char> BoolVector;
   205     typedef std::vector<Value> ValueVector;
   206     typedef std::vector<Cost> CostVector;
   207     typedef std::vector<LargeCost> LargeCostVector;
   208 
   209   private:
   210   
   211     template <typename KT, typename VT>
   212     class StaticVectorMap {
   213     public:
   214       typedef KT Key;
   215       typedef VT Value;
   216       
   217       StaticVectorMap(std::vector<Value>& v) : _v(v) {}
   218       
   219       const Value& operator[](const Key& key) const {
   220         return _v[StaticDigraph::id(key)];
   221       }
   222 
   223       Value& operator[](const Key& key) {
   224         return _v[StaticDigraph::id(key)];
   225       }
   226       
   227       void set(const Key& key, const Value& val) {
   228         _v[StaticDigraph::id(key)] = val;
   229       }
   230 
   231     private:
   232       std::vector<Value>& _v;
   233     };
   234 
   235     typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
   236     typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
   237 
   238   private:
   239 
   240     // Data related to the underlying digraph
   241     const GR &_graph;
   242     int _node_num;
   243     int _arc_num;
   244     int _res_node_num;
   245     int _res_arc_num;
   246     int _root;
   247 
   248     // Parameters of the problem
   249     bool _have_lower;
   250     Value _sum_supply;
   251 
   252     // Data structures for storing the digraph
   253     IntNodeMap _node_id;
   254     IntArcMap _arc_idf;
   255     IntArcMap _arc_idb;
   256     IntVector _first_out;
   257     BoolVector _forward;
   258     IntVector _source;
   259     IntVector _target;
   260     IntVector _reverse;
   261 
   262     // Node and arc data
   263     ValueVector _lower;
   264     ValueVector _upper;
   265     CostVector _scost;
   266     ValueVector _supply;
   267 
   268     ValueVector _res_cap;
   269     LargeCostVector _cost;
   270     LargeCostVector _pi;
   271     ValueVector _excess;
   272     IntVector _next_out;
   273     std::deque<int> _active_nodes;
   274 
   275     // Data for scaling
   276     LargeCost _epsilon;
   277     int _alpha;
   278 
   279     // Data for a StaticDigraph structure
   280     typedef std::pair<int, int> IntPair;
   281     StaticDigraph _sgr;
   282     std::vector<IntPair> _arc_vec;
   283     std::vector<LargeCost> _cost_vec;
   284     LargeCostArcMap _cost_map;
   285     LargeCostNodeMap _pi_map;
   286   
   287   public:
   288   
   289     /// \brief Constant for infinite upper bounds (capacities).
   290     ///
   291     /// Constant for infinite upper bounds (capacities).
   292     /// It is \c std::numeric_limits<Value>::infinity() if available,
   293     /// \c std::numeric_limits<Value>::max() otherwise.
   294     const Value INF;
   295 
   296   public:
   297 
   298     /// \name Named Template Parameters
   299     /// @{
   300 
   301     template <typename T>
   302     struct SetLargeCostTraits : public Traits {
   303       typedef T LargeCost;
   304     };
   305 
   306     /// \brief \ref named-templ-param "Named parameter" for setting
   307     /// \c LargeCost type.
   308     ///
   309     /// \ref named-templ-param "Named parameter" for setting \c LargeCost
   310     /// type, which is used for internal computations in the algorithm.
   311     /// \c Cost must be convertible to \c LargeCost.
   312     template <typename T>
   313     struct SetLargeCost
   314       : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
   315       typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
   316     };
   317 
   318     /// @}
   319 
   320   public:
   321 
   322     /// \brief Constructor.
   323     ///
   324     /// The constructor of the class.
   325     ///
   326     /// \param graph The digraph the algorithm runs on.
   327     CostScaling(const GR& graph) :
   328       _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
   329       _cost_map(_cost_vec), _pi_map(_pi),
   330       INF(std::numeric_limits<Value>::has_infinity ?
   331           std::numeric_limits<Value>::infinity() :
   332           std::numeric_limits<Value>::max())
   333     {
   334       // Check the number types
   335       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
   336         "The flow type of CostScaling must be signed");
   337       LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   338         "The cost type of CostScaling must be signed");
   339       
   340       // Reset data structures
   341       reset();
   342     }
   343 
   344     /// \name Parameters
   345     /// The parameters of the algorithm can be specified using these
   346     /// functions.
   347 
   348     /// @{
   349 
   350     /// \brief Set the lower bounds on the arcs.
   351     ///
   352     /// This function sets the lower bounds on the arcs.
   353     /// If it is not used before calling \ref run(), the lower bounds
   354     /// will be set to zero on all arcs.
   355     ///
   356     /// \param map An arc map storing the lower bounds.
   357     /// Its \c Value type must be convertible to the \c Value type
   358     /// of the algorithm.
   359     ///
   360     /// \return <tt>(*this)</tt>
   361     template <typename LowerMap>
   362     CostScaling& lowerMap(const LowerMap& map) {
   363       _have_lower = true;
   364       for (ArcIt a(_graph); a != INVALID; ++a) {
   365         _lower[_arc_idf[a]] = map[a];
   366         _lower[_arc_idb[a]] = map[a];
   367       }
   368       return *this;
   369     }
   370 
   371     /// \brief Set the upper bounds (capacities) on the arcs.
   372     ///
   373     /// This function sets the upper bounds (capacities) on the arcs.
   374     /// If it is not used before calling \ref run(), the upper bounds
   375     /// will be set to \ref INF on all arcs (i.e. the flow value will be
   376     /// unbounded from above).
   377     ///
   378     /// \param map An arc map storing the upper bounds.
   379     /// Its \c Value type must be convertible to the \c Value type
   380     /// of the algorithm.
   381     ///
   382     /// \return <tt>(*this)</tt>
   383     template<typename UpperMap>
   384     CostScaling& upperMap(const UpperMap& map) {
   385       for (ArcIt a(_graph); a != INVALID; ++a) {
   386         _upper[_arc_idf[a]] = map[a];
   387       }
   388       return *this;
   389     }
   390 
   391     /// \brief Set the costs of the arcs.
   392     ///
   393     /// This function sets the costs of the arcs.
   394     /// If it is not used before calling \ref run(), the costs
   395     /// will be set to \c 1 on all arcs.
   396     ///
   397     /// \param map An arc map storing the costs.
   398     /// Its \c Value type must be convertible to the \c Cost type
   399     /// of the algorithm.
   400     ///
   401     /// \return <tt>(*this)</tt>
   402     template<typename CostMap>
   403     CostScaling& costMap(const CostMap& map) {
   404       for (ArcIt a(_graph); a != INVALID; ++a) {
   405         _scost[_arc_idf[a]] =  map[a];
   406         _scost[_arc_idb[a]] = -map[a];
   407       }
   408       return *this;
   409     }
   410 
   411     /// \brief Set the supply values of the nodes.
   412     ///
   413     /// This function sets the supply values of the nodes.
   414     /// If neither this function nor \ref stSupply() is used before
   415     /// calling \ref run(), the supply of each node will be set to zero.
   416     ///
   417     /// \param map A node map storing the supply values.
   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 SupplyMap>
   423     CostScaling& supplyMap(const SupplyMap& map) {
   424       for (NodeIt n(_graph); n != INVALID; ++n) {
   425         _supply[_node_id[n]] = map[n];
   426       }
   427       return *this;
   428     }
   429 
   430     /// \brief Set single source and target nodes and a supply value.
   431     ///
   432     /// This function sets a single source node and a single target node
   433     /// and the required flow value.
   434     /// If neither this function nor \ref supplyMap() is used before
   435     /// calling \ref run(), the supply of each node will be set to zero.
   436     ///
   437     /// Using this function has the same effect as using \ref supplyMap()
   438     /// with such a map in which \c k is assigned to \c s, \c -k is
   439     /// assigned to \c t and all other nodes have zero supply value.
   440     ///
   441     /// \param s The source node.
   442     /// \param t The target node.
   443     /// \param k The required amount of flow from node \c s to node \c t
   444     /// (i.e. the supply of \c s and the demand of \c t).
   445     ///
   446     /// \return <tt>(*this)</tt>
   447     CostScaling& stSupply(const Node& s, const Node& t, Value k) {
   448       for (int i = 0; i != _res_node_num; ++i) {
   449         _supply[i] = 0;
   450       }
   451       _supply[_node_id[s]] =  k;
   452       _supply[_node_id[t]] = -k;
   453       return *this;
   454     }
   455     
   456     /// @}
   457 
   458     /// \name Execution control
   459     /// The algorithm can be executed using \ref run().
   460 
   461     /// @{
   462 
   463     /// \brief Run the algorithm.
   464     ///
   465     /// This function runs the algorithm.
   466     /// The paramters can be specified using functions \ref lowerMap(),
   467     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   468     /// For example,
   469     /// \code
   470     ///   CostScaling<ListDigraph> cs(graph);
   471     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   472     ///     .supplyMap(sup).run();
   473     /// \endcode
   474     ///
   475     /// This function can be called more than once. All the given parameters
   476     /// are kept for the next call, unless \ref resetParams() or \ref reset()
   477     /// is used, thus only the modified parameters have to be set again.
   478     /// If the underlying digraph was also modified after the construction
   479     /// of the class (or the last \ref reset() call), then the \ref reset()
   480     /// function must be called.
   481     ///
   482     /// \param method The internal method that will be used in the
   483     /// algorithm. For more information, see \ref Method.
   484     /// \param factor The cost scaling factor. It must be larger than one.
   485     ///
   486     /// \return \c INFEASIBLE if no feasible flow exists,
   487     /// \n \c OPTIMAL if the problem has optimal solution
   488     /// (i.e. it is feasible and bounded), and the algorithm has found
   489     /// optimal flow and node potentials (primal and dual solutions),
   490     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
   491     /// and infinite upper bound. It means that the objective function
   492     /// is unbounded on that arc, however, note that it could actually be
   493     /// bounded over the feasible flows, but this algroithm cannot handle
   494     /// these cases.
   495     ///
   496     /// \see ProblemType, Method
   497     /// \see resetParams(), reset()
   498     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
   499       _alpha = factor;
   500       ProblemType pt = init();
   501       if (pt != OPTIMAL) return pt;
   502       start(method);
   503       return OPTIMAL;
   504     }
   505 
   506     /// \brief Reset all the parameters that have been given before.
   507     ///
   508     /// This function resets all the paramaters that have been given
   509     /// before using functions \ref lowerMap(), \ref upperMap(),
   510     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   511     ///
   512     /// It is useful for multiple \ref run() calls. Basically, all the given
   513     /// parameters are kept for the next \ref run() call, unless
   514     /// \ref resetParams() or \ref reset() is used.
   515     /// If the underlying digraph was also modified after the construction
   516     /// of the class or the last \ref reset() call, then the \ref reset()
   517     /// function must be used, otherwise \ref resetParams() is sufficient.
   518     ///
   519     /// For example,
   520     /// \code
   521     ///   CostScaling<ListDigraph> cs(graph);
   522     ///
   523     ///   // First run
   524     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   525     ///     .supplyMap(sup).run();
   526     ///
   527     ///   // Run again with modified cost map (resetParams() is not called,
   528     ///   // so only the cost map have to be set again)
   529     ///   cost[e] += 100;
   530     ///   cs.costMap(cost).run();
   531     ///
   532     ///   // Run again from scratch using resetParams()
   533     ///   // (the lower bounds will be set to zero on all arcs)
   534     ///   cs.resetParams();
   535     ///   cs.upperMap(capacity).costMap(cost)
   536     ///     .supplyMap(sup).run();
   537     /// \endcode
   538     ///
   539     /// \return <tt>(*this)</tt>
   540     ///
   541     /// \see reset(), run()
   542     CostScaling& resetParams() {
   543       for (int i = 0; i != _res_node_num; ++i) {
   544         _supply[i] = 0;
   545       }
   546       int limit = _first_out[_root];
   547       for (int j = 0; j != limit; ++j) {
   548         _lower[j] = 0;
   549         _upper[j] = INF;
   550         _scost[j] = _forward[j] ? 1 : -1;
   551       }
   552       for (int j = limit; j != _res_arc_num; ++j) {
   553         _lower[j] = 0;
   554         _upper[j] = INF;
   555         _scost[j] = 0;
   556         _scost[_reverse[j]] = 0;
   557       }      
   558       _have_lower = false;
   559       return *this;
   560     }
   561 
   562     /// \brief Reset all the parameters that have been given before.
   563     ///
   564     /// This function resets all the paramaters that have been given
   565     /// before using functions \ref lowerMap(), \ref upperMap(),
   566     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   567     ///
   568     /// It is useful for multiple run() calls. If this function is not
   569     /// used, all the parameters given before are kept for the next
   570     /// \ref run() call.
   571     /// However, the underlying digraph must not be modified after this
   572     /// class have been constructed, since it copies and extends the graph.
   573     /// \return <tt>(*this)</tt>
   574     CostScaling& reset() {
   575       // Resize vectors
   576       _node_num = countNodes(_graph);
   577       _arc_num = countArcs(_graph);
   578       _res_node_num = _node_num + 1;
   579       _res_arc_num = 2 * (_arc_num + _node_num);
   580       _root = _node_num;
   581 
   582       _first_out.resize(_res_node_num + 1);
   583       _forward.resize(_res_arc_num);
   584       _source.resize(_res_arc_num);
   585       _target.resize(_res_arc_num);
   586       _reverse.resize(_res_arc_num);
   587 
   588       _lower.resize(_res_arc_num);
   589       _upper.resize(_res_arc_num);
   590       _scost.resize(_res_arc_num);
   591       _supply.resize(_res_node_num);
   592       
   593       _res_cap.resize(_res_arc_num);
   594       _cost.resize(_res_arc_num);
   595       _pi.resize(_res_node_num);
   596       _excess.resize(_res_node_num);
   597       _next_out.resize(_res_node_num);
   598 
   599       _arc_vec.reserve(_res_arc_num);
   600       _cost_vec.reserve(_res_arc_num);
   601 
   602       // Copy the graph
   603       int i = 0, j = 0, k = 2 * _arc_num + _node_num;
   604       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   605         _node_id[n] = i;
   606       }
   607       i = 0;
   608       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   609         _first_out[i] = j;
   610         for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   611           _arc_idf[a] = j;
   612           _forward[j] = true;
   613           _source[j] = i;
   614           _target[j] = _node_id[_graph.runningNode(a)];
   615         }
   616         for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   617           _arc_idb[a] = j;
   618           _forward[j] = false;
   619           _source[j] = i;
   620           _target[j] = _node_id[_graph.runningNode(a)];
   621         }
   622         _forward[j] = false;
   623         _source[j] = i;
   624         _target[j] = _root;
   625         _reverse[j] = k;
   626         _forward[k] = true;
   627         _source[k] = _root;
   628         _target[k] = i;
   629         _reverse[k] = j;
   630         ++j; ++k;
   631       }
   632       _first_out[i] = j;
   633       _first_out[_res_node_num] = k;
   634       for (ArcIt a(_graph); a != INVALID; ++a) {
   635         int fi = _arc_idf[a];
   636         int bi = _arc_idb[a];
   637         _reverse[fi] = bi;
   638         _reverse[bi] = fi;
   639       }
   640       
   641       // Reset parameters
   642       resetParams();
   643       return *this;
   644     }
   645 
   646     /// @}
   647 
   648     /// \name Query Functions
   649     /// The results of the algorithm can be obtained using these
   650     /// functions.\n
   651     /// The \ref run() function must be called before using them.
   652 
   653     /// @{
   654 
   655     /// \brief Return the total cost of the found flow.
   656     ///
   657     /// This function returns the total cost of the found flow.
   658     /// Its complexity is O(e).
   659     ///
   660     /// \note The return type of the function can be specified as a
   661     /// template parameter. For example,
   662     /// \code
   663     ///   cs.totalCost<double>();
   664     /// \endcode
   665     /// It is useful if the total cost cannot be stored in the \c Cost
   666     /// type of the algorithm, which is the default return type of the
   667     /// function.
   668     ///
   669     /// \pre \ref run() must be called before using this function.
   670     template <typename Number>
   671     Number totalCost() const {
   672       Number c = 0;
   673       for (ArcIt a(_graph); a != INVALID; ++a) {
   674         int i = _arc_idb[a];
   675         c += static_cast<Number>(_res_cap[i]) *
   676              (-static_cast<Number>(_scost[i]));
   677       }
   678       return c;
   679     }
   680 
   681 #ifndef DOXYGEN
   682     Cost totalCost() const {
   683       return totalCost<Cost>();
   684     }
   685 #endif
   686 
   687     /// \brief Return the flow on the given arc.
   688     ///
   689     /// This function returns the flow on the given arc.
   690     ///
   691     /// \pre \ref run() must be called before using this function.
   692     Value flow(const Arc& a) const {
   693       return _res_cap[_arc_idb[a]];
   694     }
   695 
   696     /// \brief Return the flow map (the primal solution).
   697     ///
   698     /// This function copies the flow value on each arc into the given
   699     /// map. The \c Value type of the algorithm must be convertible to
   700     /// the \c Value type of the map.
   701     ///
   702     /// \pre \ref run() must be called before using this function.
   703     template <typename FlowMap>
   704     void flowMap(FlowMap &map) const {
   705       for (ArcIt a(_graph); a != INVALID; ++a) {
   706         map.set(a, _res_cap[_arc_idb[a]]);
   707       }
   708     }
   709 
   710     /// \brief Return the potential (dual value) of the given node.
   711     ///
   712     /// This function returns the potential (dual value) of the
   713     /// given node.
   714     ///
   715     /// \pre \ref run() must be called before using this function.
   716     Cost potential(const Node& n) const {
   717       return static_cast<Cost>(_pi[_node_id[n]]);
   718     }
   719 
   720     /// \brief Return the potential map (the dual solution).
   721     ///
   722     /// This function copies the potential (dual value) of each node
   723     /// into the given map.
   724     /// The \c Cost type of the algorithm must be convertible to the
   725     /// \c Value type of the map.
   726     ///
   727     /// \pre \ref run() must be called before using this function.
   728     template <typename PotentialMap>
   729     void potentialMap(PotentialMap &map) const {
   730       for (NodeIt n(_graph); n != INVALID; ++n) {
   731         map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
   732       }
   733     }
   734 
   735     /// @}
   736 
   737   private:
   738 
   739     // Initialize the algorithm
   740     ProblemType init() {
   741       if (_res_node_num <= 1) return INFEASIBLE;
   742 
   743       // Check the sum of supply values
   744       _sum_supply = 0;
   745       for (int i = 0; i != _root; ++i) {
   746         _sum_supply += _supply[i];
   747       }
   748       if (_sum_supply > 0) return INFEASIBLE;
   749       
   750 
   751       // Initialize vectors
   752       for (int i = 0; i != _res_node_num; ++i) {
   753         _pi[i] = 0;
   754         _excess[i] = _supply[i];
   755       }
   756       
   757       // Remove infinite upper bounds and check negative arcs
   758       const Value MAX = std::numeric_limits<Value>::max();
   759       int last_out;
   760       if (_have_lower) {
   761         for (int i = 0; i != _root; ++i) {
   762           last_out = _first_out[i+1];
   763           for (int j = _first_out[i]; j != last_out; ++j) {
   764             if (_forward[j]) {
   765               Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
   766               if (c >= MAX) return UNBOUNDED;
   767               _excess[i] -= c;
   768               _excess[_target[j]] += c;
   769             }
   770           }
   771         }
   772       } else {
   773         for (int i = 0; i != _root; ++i) {
   774           last_out = _first_out[i+1];
   775           for (int j = _first_out[i]; j != last_out; ++j) {
   776             if (_forward[j] && _scost[j] < 0) {
   777               Value c = _upper[j];
   778               if (c >= MAX) return UNBOUNDED;
   779               _excess[i] -= c;
   780               _excess[_target[j]] += c;
   781             }
   782           }
   783         }
   784       }
   785       Value ex, max_cap = 0;
   786       for (int i = 0; i != _res_node_num; ++i) {
   787         ex = _excess[i];
   788         _excess[i] = 0;
   789         if (ex < 0) max_cap -= ex;
   790       }
   791       for (int j = 0; j != _res_arc_num; ++j) {
   792         if (_upper[j] >= MAX) _upper[j] = max_cap;
   793       }
   794 
   795       // Initialize the large cost vector and the epsilon parameter
   796       _epsilon = 0;
   797       LargeCost lc;
   798       for (int i = 0; i != _root; ++i) {
   799         last_out = _first_out[i+1];
   800         for (int j = _first_out[i]; j != last_out; ++j) {
   801           lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
   802           _cost[j] = lc;
   803           if (lc > _epsilon) _epsilon = lc;
   804         }
   805       }
   806       _epsilon /= _alpha;
   807 
   808       // Initialize maps for Circulation and remove non-zero lower bounds
   809       ConstMap<Arc, Value> low(0);
   810       typedef typename Digraph::template ArcMap<Value> ValueArcMap;
   811       typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
   812       ValueArcMap cap(_graph), flow(_graph);
   813       ValueNodeMap sup(_graph);
   814       for (NodeIt n(_graph); n != INVALID; ++n) {
   815         sup[n] = _supply[_node_id[n]];
   816       }
   817       if (_have_lower) {
   818         for (ArcIt a(_graph); a != INVALID; ++a) {
   819           int j = _arc_idf[a];
   820           Value c = _lower[j];
   821           cap[a] = _upper[j] - c;
   822           sup[_graph.source(a)] -= c;
   823           sup[_graph.target(a)] += c;
   824         }
   825       } else {
   826         for (ArcIt a(_graph); a != INVALID; ++a) {
   827           cap[a] = _upper[_arc_idf[a]];
   828         }
   829       }
   830 
   831       // Find a feasible flow using Circulation
   832       Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
   833         circ(_graph, low, cap, sup);
   834       if (!circ.flowMap(flow).run()) return INFEASIBLE;
   835 
   836       // Set residual capacities and handle GEQ supply type
   837       if (_sum_supply < 0) {
   838         for (ArcIt a(_graph); a != INVALID; ++a) {
   839           Value fa = flow[a];
   840           _res_cap[_arc_idf[a]] = cap[a] - fa;
   841           _res_cap[_arc_idb[a]] = fa;
   842           sup[_graph.source(a)] -= fa;
   843           sup[_graph.target(a)] += fa;
   844         }
   845         for (NodeIt n(_graph); n != INVALID; ++n) {
   846           _excess[_node_id[n]] = sup[n];
   847         }
   848         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   849           int u = _target[a];
   850           int ra = _reverse[a];
   851           _res_cap[a] = -_sum_supply + 1;
   852           _res_cap[ra] = -_excess[u];
   853           _cost[a] = 0;
   854           _cost[ra] = 0;
   855           _excess[u] = 0;
   856         }
   857       } else {
   858         for (ArcIt a(_graph); a != INVALID; ++a) {
   859           Value fa = flow[a];
   860           _res_cap[_arc_idf[a]] = cap[a] - fa;
   861           _res_cap[_arc_idb[a]] = fa;
   862         }
   863         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   864           int ra = _reverse[a];
   865           _res_cap[a] = 1;
   866           _res_cap[ra] = 0;
   867           _cost[a] = 0;
   868           _cost[ra] = 0;
   869         }
   870       }
   871       
   872       return OPTIMAL;
   873     }
   874 
   875     // Execute the algorithm and transform the results
   876     void start(Method method) {
   877       // Maximum path length for partial augment
   878       const int MAX_PATH_LENGTH = 4;
   879       
   880       // Execute the algorithm
   881       switch (method) {
   882         case PUSH:
   883           startPush();
   884           break;
   885         case AUGMENT:
   886           startAugment();
   887           break;
   888         case PARTIAL_AUGMENT:
   889           startAugment(MAX_PATH_LENGTH);
   890           break;
   891       }
   892 
   893       // Compute node potentials for the original costs
   894       _arc_vec.clear();
   895       _cost_vec.clear();
   896       for (int j = 0; j != _res_arc_num; ++j) {
   897         if (_res_cap[j] > 0) {
   898           _arc_vec.push_back(IntPair(_source[j], _target[j]));
   899           _cost_vec.push_back(_scost[j]);
   900         }
   901       }
   902       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   903 
   904       typename BellmanFord<StaticDigraph, LargeCostArcMap>
   905         ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
   906       bf.distMap(_pi_map);
   907       bf.init(0);
   908       bf.start();
   909 
   910       // Handle non-zero lower bounds
   911       if (_have_lower) {
   912         int limit = _first_out[_root];
   913         for (int j = 0; j != limit; ++j) {
   914           if (!_forward[j]) _res_cap[j] += _lower[j];
   915         }
   916       }
   917     }
   918 
   919     /// Execute the algorithm performing augment and relabel operations
   920     void startAugment(int max_length = std::numeric_limits<int>::max()) {
   921       // Paramters for heuristics
   922       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
   923       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
   924 
   925       // Perform cost scaling phases
   926       IntVector pred_arc(_res_node_num);
   927       std::vector<int> path_nodes;
   928       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   929                                         1 : _epsilon / _alpha )
   930       {
   931         // "Early Termination" heuristic: use Bellman-Ford algorithm
   932         // to check if the current flow is optimal
   933         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
   934           _arc_vec.clear();
   935           _cost_vec.clear();
   936           for (int j = 0; j != _res_arc_num; ++j) {
   937             if (_res_cap[j] > 0) {
   938               _arc_vec.push_back(IntPair(_source[j], _target[j]));
   939               _cost_vec.push_back(_cost[j] + 1);
   940             }
   941           }
   942           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   943 
   944           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   945           bf.init(0);
   946           bool done = false;
   947           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
   948           for (int i = 0; i < K && !done; ++i)
   949             done = bf.processNextWeakRound();
   950           if (done) break;
   951         }
   952 
   953         // Saturate arcs not satisfying the optimality condition
   954         for (int a = 0; a != _res_arc_num; ++a) {
   955           if (_res_cap[a] > 0 &&
   956               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   957             Value delta = _res_cap[a];
   958             _excess[_source[a]] -= delta;
   959             _excess[_target[a]] += delta;
   960             _res_cap[a] = 0;
   961             _res_cap[_reverse[a]] += delta;
   962           }
   963         }
   964         
   965         // Find active nodes (i.e. nodes with positive excess)
   966         for (int u = 0; u != _res_node_num; ++u) {
   967           if (_excess[u] > 0) _active_nodes.push_back(u);
   968         }
   969 
   970         // Initialize the next arcs
   971         for (int u = 0; u != _res_node_num; ++u) {
   972           _next_out[u] = _first_out[u];
   973         }
   974 
   975         // Perform partial augment and relabel operations
   976         while (true) {
   977           // Select an active node (FIFO selection)
   978           while (_active_nodes.size() > 0 &&
   979                  _excess[_active_nodes.front()] <= 0) {
   980             _active_nodes.pop_front();
   981           }
   982           if (_active_nodes.size() == 0) break;
   983           int start = _active_nodes.front();
   984           path_nodes.clear();
   985           path_nodes.push_back(start);
   986 
   987           // Find an augmenting path from the start node
   988           int tip = start;
   989           while (_excess[tip] >= 0 &&
   990                  int(path_nodes.size()) <= max_length) {
   991             int u;
   992             LargeCost min_red_cost, rc;
   993             int last_out = _sum_supply < 0 ?
   994               _first_out[tip+1] : _first_out[tip+1] - 1;
   995             for (int a = _next_out[tip]; a != last_out; ++a) {
   996               if (_res_cap[a] > 0 &&
   997                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   998                 u = _target[a];
   999                 pred_arc[u] = a;
  1000                 _next_out[tip] = a;
  1001                 tip = u;
  1002                 path_nodes.push_back(tip);
  1003                 goto next_step;
  1004               }
  1005             }
  1006 
  1007             // Relabel tip node
  1008             min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
  1009             for (int a = _first_out[tip]; a != last_out; ++a) {
  1010               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
  1011               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1012                 min_red_cost = rc;
  1013               }
  1014             }
  1015             _pi[tip] -= min_red_cost + _epsilon;
  1016 
  1017             // Reset the next arc of tip
  1018             _next_out[tip] = _first_out[tip];
  1019 
  1020             // Step back
  1021             if (tip != start) {
  1022               path_nodes.pop_back();
  1023               tip = path_nodes.back();
  1024             }
  1025 
  1026           next_step: ;
  1027           }
  1028 
  1029           // Augment along the found path (as much flow as possible)
  1030           Value delta;
  1031           int u, v = path_nodes.front(), pa;
  1032           for (int i = 1; i < int(path_nodes.size()); ++i) {
  1033             u = v;
  1034             v = path_nodes[i];
  1035             pa = pred_arc[v];
  1036             delta = std::min(_res_cap[pa], _excess[u]);
  1037             _res_cap[pa] -= delta;
  1038             _res_cap[_reverse[pa]] += delta;
  1039             _excess[u] -= delta;
  1040             _excess[v] += delta;
  1041             if (_excess[v] > 0 && _excess[v] <= delta)
  1042               _active_nodes.push_back(v);
  1043           }
  1044         }
  1045       }
  1046     }
  1047 
  1048     /// Execute the algorithm performing push and relabel operations
  1049     void startPush() {
  1050       // Paramters for heuristics
  1051       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
  1052       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
  1053 
  1054       // Perform cost scaling phases
  1055       BoolVector hyper(_res_node_num, false);
  1056       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1057                                         1 : _epsilon / _alpha )
  1058       {
  1059         // "Early Termination" heuristic: use Bellman-Ford algorithm
  1060         // to check if the current flow is optimal
  1061         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
  1062           _arc_vec.clear();
  1063           _cost_vec.clear();
  1064           for (int j = 0; j != _res_arc_num; ++j) {
  1065             if (_res_cap[j] > 0) {
  1066               _arc_vec.push_back(IntPair(_source[j], _target[j]));
  1067               _cost_vec.push_back(_cost[j] + 1);
  1068             }
  1069           }
  1070           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
  1071 
  1072           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
  1073           bf.init(0);
  1074           bool done = false;
  1075           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
  1076           for (int i = 0; i < K && !done; ++i)
  1077             done = bf.processNextWeakRound();
  1078           if (done) break;
  1079         }
  1080 
  1081         // Saturate arcs not satisfying the optimality condition
  1082         for (int a = 0; a != _res_arc_num; ++a) {
  1083           if (_res_cap[a] > 0 &&
  1084               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
  1085             Value delta = _res_cap[a];
  1086             _excess[_source[a]] -= delta;
  1087             _excess[_target[a]] += delta;
  1088             _res_cap[a] = 0;
  1089             _res_cap[_reverse[a]] += delta;
  1090           }
  1091         }
  1092 
  1093         // Find active nodes (i.e. nodes with positive excess)
  1094         for (int u = 0; u != _res_node_num; ++u) {
  1095           if (_excess[u] > 0) _active_nodes.push_back(u);
  1096         }
  1097 
  1098         // Initialize the next arcs
  1099         for (int u = 0; u != _res_node_num; ++u) {
  1100           _next_out[u] = _first_out[u];
  1101         }
  1102 
  1103         // Perform push and relabel operations
  1104         while (_active_nodes.size() > 0) {
  1105           LargeCost min_red_cost, rc;
  1106           Value delta;
  1107           int n, t, a, last_out = _res_arc_num;
  1108 
  1109           // Select an active node (FIFO selection)
  1110         next_node:
  1111           n = _active_nodes.front();
  1112           last_out = _sum_supply < 0 ?
  1113             _first_out[n+1] : _first_out[n+1] - 1;
  1114 
  1115           // Perform push operations if there are admissible arcs
  1116           if (_excess[n] > 0) {
  1117             for (a = _next_out[n]; a != last_out; ++a) {
  1118               if (_res_cap[a] > 0 &&
  1119                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
  1120                 delta = std::min(_res_cap[a], _excess[n]);
  1121                 t = _target[a];
  1122 
  1123                 // Push-look-ahead heuristic
  1124                 Value ahead = -_excess[t];
  1125                 int last_out_t = _sum_supply < 0 ?
  1126                   _first_out[t+1] : _first_out[t+1] - 1;
  1127                 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
  1128                   if (_res_cap[ta] > 0 && 
  1129                       _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
  1130                     ahead += _res_cap[ta];
  1131                   if (ahead >= delta) break;
  1132                 }
  1133                 if (ahead < 0) ahead = 0;
  1134 
  1135                 // Push flow along the arc
  1136                 if (ahead < delta) {
  1137                   _res_cap[a] -= ahead;
  1138                   _res_cap[_reverse[a]] += ahead;
  1139                   _excess[n] -= ahead;
  1140                   _excess[t] += ahead;
  1141                   _active_nodes.push_front(t);
  1142                   hyper[t] = true;
  1143                   _next_out[n] = a;
  1144                   goto next_node;
  1145                 } else {
  1146                   _res_cap[a] -= delta;
  1147                   _res_cap[_reverse[a]] += delta;
  1148                   _excess[n] -= delta;
  1149                   _excess[t] += delta;
  1150                   if (_excess[t] > 0 && _excess[t] <= delta)
  1151                     _active_nodes.push_back(t);
  1152                 }
  1153 
  1154                 if (_excess[n] == 0) {
  1155                   _next_out[n] = a;
  1156                   goto remove_nodes;
  1157                 }
  1158               }
  1159             }
  1160             _next_out[n] = a;
  1161           }
  1162 
  1163           // Relabel the node if it is still active (or hyper)
  1164           if (_excess[n] > 0 || hyper[n]) {
  1165             min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
  1166             for (int a = _first_out[n]; a != last_out; ++a) {
  1167               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
  1168               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1169                 min_red_cost = rc;
  1170               }
  1171             }
  1172             _pi[n] -= min_red_cost + _epsilon;
  1173             hyper[n] = false;
  1174 
  1175             // Reset the next arc
  1176             _next_out[n] = _first_out[n];
  1177           }
  1178         
  1179           // Remove nodes that are not active nor hyper
  1180         remove_nodes:
  1181           while ( _active_nodes.size() > 0 &&
  1182                   _excess[_active_nodes.front()] <= 0 &&
  1183                   !hyper[_active_nodes.front()] ) {
  1184             _active_nodes.pop_front();
  1185           }
  1186         }
  1187       }
  1188     }
  1189 
  1190   }; //class CostScaling
  1191 
  1192   ///@}
  1193 
  1194 } //namespace lemon
  1195 
  1196 #endif //LEMON_COST_SCALING_H