lemon/cost_scaling.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sun, 23 Oct 2011 22:42:14 +0200
changeset 955 7f6eeffe3cd1
parent 937 1226290a9b7d
child 1003 16f55008c863
permissions -rw-r--r--
Fix windows-specific compilation error (#429)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2010
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_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   /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
   101   /// implementations available in LEMON for this problem.
   102   ///
   103   /// Most of the parameters of the problem (except for the digraph)
   104   /// can be given using separate functions, and the algorithm can be
   105   /// executed using the \ref run() function. If some parameters are not
   106   /// specified, then default values will be used.
   107   ///
   108   /// \tparam GR The digraph type the algorithm runs on.
   109   /// \tparam V The number type used for flow amounts, capacity bounds
   110   /// and supply values in the algorithm. By default, it is \c int.
   111   /// \tparam C The number type used for costs and potentials in the
   112   /// algorithm. By default, it is the same as \c V.
   113   /// \tparam TR The traits class that defines various types used by the
   114   /// algorithm. By default, it is \ref CostScalingDefaultTraits
   115   /// "CostScalingDefaultTraits<GR, V, C>".
   116   /// In most cases, this parameter should not be set directly,
   117   /// consider to use the named template parameters instead.
   118   ///
   119   /// \warning Both \c V and \c C must be signed number types.
   120   /// \warning All input data (capacities, supply values, and costs) must
   121   /// be integer.
   122   /// \warning This algorithm does not support negative costs for
   123   /// arcs having infinite upper bound.
   124   ///
   125   /// \note %CostScaling provides three different internal methods,
   126   /// from which the most efficient one is used by default.
   127   /// For more information, see \ref Method.
   128 #ifdef DOXYGEN
   129   template <typename GR, typename V, typename C, typename TR>
   130 #else
   131   template < typename GR, typename V = int, typename C = V,
   132              typename TR = CostScalingDefaultTraits<GR, V, C> >
   133 #endif
   134   class CostScaling
   135   {
   136   public:
   137 
   138     /// The type of the digraph
   139     typedef typename TR::Digraph Digraph;
   140     /// The type of the flow amounts, capacity bounds and supply values
   141     typedef typename TR::Value Value;
   142     /// The type of the arc costs
   143     typedef typename TR::Cost Cost;
   144 
   145     /// \brief The large cost type
   146     ///
   147     /// The large cost type used for internal computations.
   148     /// By default, it is \c long \c long if the \c Cost type is integer,
   149     /// otherwise it is \c double.
   150     typedef typename TR::LargeCost LargeCost;
   151 
   152     /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
   153     typedef TR Traits;
   154 
   155   public:
   156 
   157     /// \brief Problem type constants for the \c run() function.
   158     ///
   159     /// Enum type containing the problem type constants that can be
   160     /// returned by the \ref run() function of the algorithm.
   161     enum ProblemType {
   162       /// The problem has no feasible solution (flow).
   163       INFEASIBLE,
   164       /// The problem has optimal solution (i.e. it is feasible and
   165       /// bounded), and the algorithm has found optimal flow and node
   166       /// potentials (primal and dual solutions).
   167       OPTIMAL,
   168       /// The digraph contains an arc of negative cost and infinite
   169       /// upper bound. It means that the objective function is unbounded
   170       /// on that arc, however, note that it could actually be bounded
   171       /// over the feasible flows, but this algroithm cannot handle
   172       /// these cases.
   173       UNBOUNDED
   174     };
   175 
   176     /// \brief Constants for selecting the internal method.
   177     ///
   178     /// Enum type containing constants for selecting the internal method
   179     /// for the \ref run() function.
   180     ///
   181     /// \ref CostScaling provides three internal methods that differ mainly
   182     /// in their base operations, which are used in conjunction with the
   183     /// relabel operation.
   184     /// By default, the so called \ref PARTIAL_AUGMENT
   185     /// "Partial Augment-Relabel" method is used, which turned out to be
   186     /// the most efficient and the most robust on various test inputs.
   187     /// However, the other methods can be selected using the \ref run()
   188     /// function with the proper parameter.
   189     enum Method {
   190       /// Local push operations are used, i.e. flow is moved only on one
   191       /// admissible arc at once.
   192       PUSH,
   193       /// Augment operations are used, i.e. flow is moved on admissible
   194       /// paths from a node with excess to a node with deficit.
   195       AUGMENT,
   196       /// Partial augment operations are used, i.e. flow is moved on
   197       /// admissible paths started from a node with excess, but the
   198       /// lengths of these paths are limited. This method can be viewed
   199       /// as a combined version of the previous two operations.
   200       PARTIAL_AUGMENT
   201     };
   202 
   203   private:
   204 
   205     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   206 
   207     typedef std::vector<int> IntVector;
   208     typedef std::vector<Value> ValueVector;
   209     typedef std::vector<Cost> CostVector;
   210     typedef std::vector<LargeCost> LargeCostVector;
   211     typedef std::vector<char> BoolVector;
   212     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
   213 
   214   private:
   215 
   216     template <typename KT, typename VT>
   217     class StaticVectorMap {
   218     public:
   219       typedef KT Key;
   220       typedef VT Value;
   221 
   222       StaticVectorMap(std::vector<Value>& v) : _v(v) {}
   223 
   224       const Value& operator[](const Key& key) const {
   225         return _v[StaticDigraph::id(key)];
   226       }
   227 
   228       Value& operator[](const Key& key) {
   229         return _v[StaticDigraph::id(key)];
   230       }
   231 
   232       void set(const Key& key, const Value& val) {
   233         _v[StaticDigraph::id(key)] = val;
   234       }
   235 
   236     private:
   237       std::vector<Value>& _v;
   238     };
   239 
   240     typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
   241 
   242   private:
   243 
   244     // Data related to the underlying digraph
   245     const GR &_graph;
   246     int _node_num;
   247     int _arc_num;
   248     int _res_node_num;
   249     int _res_arc_num;
   250     int _root;
   251 
   252     // Parameters of the problem
   253     bool _have_lower;
   254     Value _sum_supply;
   255     int _sup_node_num;
   256 
   257     // Data structures for storing the digraph
   258     IntNodeMap _node_id;
   259     IntArcMap _arc_idf;
   260     IntArcMap _arc_idb;
   261     IntVector _first_out;
   262     BoolVector _forward;
   263     IntVector _source;
   264     IntVector _target;
   265     IntVector _reverse;
   266 
   267     // Node and arc data
   268     ValueVector _lower;
   269     ValueVector _upper;
   270     CostVector _scost;
   271     ValueVector _supply;
   272 
   273     ValueVector _res_cap;
   274     LargeCostVector _cost;
   275     LargeCostVector _pi;
   276     ValueVector _excess;
   277     IntVector _next_out;
   278     std::deque<int> _active_nodes;
   279 
   280     // Data for scaling
   281     LargeCost _epsilon;
   282     int _alpha;
   283 
   284     IntVector _buckets;
   285     IntVector _bucket_next;
   286     IntVector _bucket_prev;
   287     IntVector _rank;
   288     int _max_rank;
   289 
   290   public:
   291 
   292     /// \brief Constant for infinite upper bounds (capacities).
   293     ///
   294     /// Constant for infinite upper bounds (capacities).
   295     /// It is \c std::numeric_limits<Value>::infinity() if available,
   296     /// \c std::numeric_limits<Value>::max() otherwise.
   297     const Value INF;
   298 
   299   public:
   300 
   301     /// \name Named Template Parameters
   302     /// @{
   303 
   304     template <typename T>
   305     struct SetLargeCostTraits : public Traits {
   306       typedef T LargeCost;
   307     };
   308 
   309     /// \brief \ref named-templ-param "Named parameter" for setting
   310     /// \c LargeCost type.
   311     ///
   312     /// \ref named-templ-param "Named parameter" for setting \c LargeCost
   313     /// type, which is used for internal computations in the algorithm.
   314     /// \c Cost must be convertible to \c LargeCost.
   315     template <typename T>
   316     struct SetLargeCost
   317       : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
   318       typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
   319     };
   320 
   321     /// @}
   322 
   323   protected:
   324 
   325     CostScaling() {}
   326 
   327   public:
   328 
   329     /// \brief Constructor.
   330     ///
   331     /// The constructor of the class.
   332     ///
   333     /// \param graph The digraph the algorithm runs on.
   334     CostScaling(const GR& graph) :
   335       _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
   336       INF(std::numeric_limits<Value>::has_infinity ?
   337           std::numeric_limits<Value>::infinity() :
   338           std::numeric_limits<Value>::max())
   339     {
   340       // Check the number types
   341       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
   342         "The flow type of CostScaling must be signed");
   343       LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   344         "The cost type of CostScaling must be signed");
   345 
   346       // Reset data structures
   347       reset();
   348     }
   349 
   350     /// \name Parameters
   351     /// The parameters of the algorithm can be specified using these
   352     /// functions.
   353 
   354     /// @{
   355 
   356     /// \brief Set the lower bounds on the arcs.
   357     ///
   358     /// This function sets the lower bounds on the arcs.
   359     /// If it is not used before calling \ref run(), the lower bounds
   360     /// will be set to zero on all arcs.
   361     ///
   362     /// \param map An arc map storing the lower bounds.
   363     /// Its \c Value type must be convertible to the \c Value type
   364     /// of the algorithm.
   365     ///
   366     /// \return <tt>(*this)</tt>
   367     template <typename LowerMap>
   368     CostScaling& lowerMap(const LowerMap& map) {
   369       _have_lower = true;
   370       for (ArcIt a(_graph); a != INVALID; ++a) {
   371         _lower[_arc_idf[a]] = map[a];
   372         _lower[_arc_idb[a]] = map[a];
   373       }
   374       return *this;
   375     }
   376 
   377     /// \brief Set the upper bounds (capacities) on the arcs.
   378     ///
   379     /// This function sets the upper bounds (capacities) on the arcs.
   380     /// If it is not used before calling \ref run(), the upper bounds
   381     /// will be set to \ref INF on all arcs (i.e. the flow value will be
   382     /// unbounded from above).
   383     ///
   384     /// \param map An arc map storing the upper bounds.
   385     /// Its \c Value type must be convertible to the \c Value type
   386     /// of the algorithm.
   387     ///
   388     /// \return <tt>(*this)</tt>
   389     template<typename UpperMap>
   390     CostScaling& upperMap(const UpperMap& map) {
   391       for (ArcIt a(_graph); a != INVALID; ++a) {
   392         _upper[_arc_idf[a]] = map[a];
   393       }
   394       return *this;
   395     }
   396 
   397     /// \brief Set the costs of the arcs.
   398     ///
   399     /// This function sets the costs of the arcs.
   400     /// If it is not used before calling \ref run(), the costs
   401     /// will be set to \c 1 on all arcs.
   402     ///
   403     /// \param map An arc map storing the costs.
   404     /// Its \c Value type must be convertible to the \c Cost type
   405     /// of the algorithm.
   406     ///
   407     /// \return <tt>(*this)</tt>
   408     template<typename CostMap>
   409     CostScaling& costMap(const CostMap& map) {
   410       for (ArcIt a(_graph); a != INVALID; ++a) {
   411         _scost[_arc_idf[a]] =  map[a];
   412         _scost[_arc_idb[a]] = -map[a];
   413       }
   414       return *this;
   415     }
   416 
   417     /// \brief Set the supply values of the nodes.
   418     ///
   419     /// This function sets the supply values of the nodes.
   420     /// If neither this function nor \ref stSupply() is used before
   421     /// calling \ref run(), the supply of each node will be set to zero.
   422     ///
   423     /// \param map A node map storing the supply values.
   424     /// Its \c Value type must be convertible to the \c Value type
   425     /// of the algorithm.
   426     ///
   427     /// \return <tt>(*this)</tt>
   428     template<typename SupplyMap>
   429     CostScaling& supplyMap(const SupplyMap& map) {
   430       for (NodeIt n(_graph); n != INVALID; ++n) {
   431         _supply[_node_id[n]] = map[n];
   432       }
   433       return *this;
   434     }
   435 
   436     /// \brief Set single source and target nodes and a supply value.
   437     ///
   438     /// This function sets a single source node and a single target node
   439     /// and the required flow value.
   440     /// If neither this function nor \ref supplyMap() is used before
   441     /// calling \ref run(), the supply of each node will be set to zero.
   442     ///
   443     /// Using this function has the same effect as using \ref supplyMap()
   444     /// with a map in which \c k is assigned to \c s, \c -k is
   445     /// assigned to \c t and all other nodes have zero supply value.
   446     ///
   447     /// \param s The source node.
   448     /// \param t The target node.
   449     /// \param k The required amount of flow from node \c s to node \c t
   450     /// (i.e. the supply of \c s and the demand of \c t).
   451     ///
   452     /// \return <tt>(*this)</tt>
   453     CostScaling& stSupply(const Node& s, const Node& t, Value k) {
   454       for (int i = 0; i != _res_node_num; ++i) {
   455         _supply[i] = 0;
   456       }
   457       _supply[_node_id[s]] =  k;
   458       _supply[_node_id[t]] = -k;
   459       return *this;
   460     }
   461 
   462     /// @}
   463 
   464     /// \name Execution control
   465     /// The algorithm can be executed using \ref run().
   466 
   467     /// @{
   468 
   469     /// \brief Run the algorithm.
   470     ///
   471     /// This function runs the algorithm.
   472     /// The paramters can be specified using functions \ref lowerMap(),
   473     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   474     /// For example,
   475     /// \code
   476     ///   CostScaling<ListDigraph> cs(graph);
   477     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   478     ///     .supplyMap(sup).run();
   479     /// \endcode
   480     ///
   481     /// This function can be called more than once. All the given parameters
   482     /// are kept for the next call, unless \ref resetParams() or \ref reset()
   483     /// is used, thus only the modified parameters have to be set again.
   484     /// If the underlying digraph was also modified after the construction
   485     /// of the class (or the last \ref reset() call), then the \ref reset()
   486     /// function must be called.
   487     ///
   488     /// \param method The internal method that will be used in the
   489     /// algorithm. For more information, see \ref Method.
   490     /// \param factor The cost scaling factor. It must be at least two.
   491     ///
   492     /// \return \c INFEASIBLE if no feasible flow exists,
   493     /// \n \c OPTIMAL if the problem has optimal solution
   494     /// (i.e. it is feasible and bounded), and the algorithm has found
   495     /// optimal flow and node potentials (primal and dual solutions),
   496     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
   497     /// and infinite upper bound. It means that the objective function
   498     /// is unbounded on that arc, however, note that it could actually be
   499     /// bounded over the feasible flows, but this algroithm cannot handle
   500     /// these cases.
   501     ///
   502     /// \see ProblemType, Method
   503     /// \see resetParams(), reset()
   504     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) {
   505       LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2");
   506       _alpha = factor;
   507       ProblemType pt = init();
   508       if (pt != OPTIMAL) return pt;
   509       start(method);
   510       return OPTIMAL;
   511     }
   512 
   513     /// \brief Reset all the parameters that have been given before.
   514     ///
   515     /// This function resets all the paramaters that have been given
   516     /// before using functions \ref lowerMap(), \ref upperMap(),
   517     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   518     ///
   519     /// It is useful for multiple \ref run() calls. Basically, all the given
   520     /// parameters are kept for the next \ref run() call, unless
   521     /// \ref resetParams() or \ref reset() is used.
   522     /// If the underlying digraph was also modified after the construction
   523     /// of the class or the last \ref reset() call, then the \ref reset()
   524     /// function must be used, otherwise \ref resetParams() is sufficient.
   525     ///
   526     /// For example,
   527     /// \code
   528     ///   CostScaling<ListDigraph> cs(graph);
   529     ///
   530     ///   // First run
   531     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   532     ///     .supplyMap(sup).run();
   533     ///
   534     ///   // Run again with modified cost map (resetParams() is not called,
   535     ///   // so only the cost map have to be set again)
   536     ///   cost[e] += 100;
   537     ///   cs.costMap(cost).run();
   538     ///
   539     ///   // Run again from scratch using resetParams()
   540     ///   // (the lower bounds will be set to zero on all arcs)
   541     ///   cs.resetParams();
   542     ///   cs.upperMap(capacity).costMap(cost)
   543     ///     .supplyMap(sup).run();
   544     /// \endcode
   545     ///
   546     /// \return <tt>(*this)</tt>
   547     ///
   548     /// \see reset(), run()
   549     CostScaling& resetParams() {
   550       for (int i = 0; i != _res_node_num; ++i) {
   551         _supply[i] = 0;
   552       }
   553       int limit = _first_out[_root];
   554       for (int j = 0; j != limit; ++j) {
   555         _lower[j] = 0;
   556         _upper[j] = INF;
   557         _scost[j] = _forward[j] ? 1 : -1;
   558       }
   559       for (int j = limit; j != _res_arc_num; ++j) {
   560         _lower[j] = 0;
   561         _upper[j] = INF;
   562         _scost[j] = 0;
   563         _scost[_reverse[j]] = 0;
   564       }
   565       _have_lower = false;
   566       return *this;
   567     }
   568 
   569     /// \brief Reset the internal data structures and all the parameters
   570     /// that have been given before.
   571     ///
   572     /// This function resets the internal data structures and all the
   573     /// paramaters that have been given before using functions \ref lowerMap(),
   574     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   575     ///
   576     /// It is useful for multiple \ref run() calls. By default, all the given
   577     /// parameters are kept for the next \ref run() call, unless
   578     /// \ref resetParams() or \ref reset() is used.
   579     /// If the underlying digraph was also modified after the construction
   580     /// of the class or the last \ref reset() call, then the \ref reset()
   581     /// function must be used, otherwise \ref resetParams() is sufficient.
   582     ///
   583     /// See \ref resetParams() for examples.
   584     ///
   585     /// \return <tt>(*this)</tt>
   586     ///
   587     /// \see resetParams(), run()
   588     CostScaling& reset() {
   589       // Resize vectors
   590       _node_num = countNodes(_graph);
   591       _arc_num = countArcs(_graph);
   592       _res_node_num = _node_num + 1;
   593       _res_arc_num = 2 * (_arc_num + _node_num);
   594       _root = _node_num;
   595 
   596       _first_out.resize(_res_node_num + 1);
   597       _forward.resize(_res_arc_num);
   598       _source.resize(_res_arc_num);
   599       _target.resize(_res_arc_num);
   600       _reverse.resize(_res_arc_num);
   601 
   602       _lower.resize(_res_arc_num);
   603       _upper.resize(_res_arc_num);
   604       _scost.resize(_res_arc_num);
   605       _supply.resize(_res_node_num);
   606 
   607       _res_cap.resize(_res_arc_num);
   608       _cost.resize(_res_arc_num);
   609       _pi.resize(_res_node_num);
   610       _excess.resize(_res_node_num);
   611       _next_out.resize(_res_node_num);
   612 
   613       // Copy the graph
   614       int i = 0, j = 0, k = 2 * _arc_num + _node_num;
   615       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   616         _node_id[n] = i;
   617       }
   618       i = 0;
   619       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   620         _first_out[i] = j;
   621         for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   622           _arc_idf[a] = j;
   623           _forward[j] = true;
   624           _source[j] = i;
   625           _target[j] = _node_id[_graph.runningNode(a)];
   626         }
   627         for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   628           _arc_idb[a] = j;
   629           _forward[j] = false;
   630           _source[j] = i;
   631           _target[j] = _node_id[_graph.runningNode(a)];
   632         }
   633         _forward[j] = false;
   634         _source[j] = i;
   635         _target[j] = _root;
   636         _reverse[j] = k;
   637         _forward[k] = true;
   638         _source[k] = _root;
   639         _target[k] = i;
   640         _reverse[k] = j;
   641         ++j; ++k;
   642       }
   643       _first_out[i] = j;
   644       _first_out[_res_node_num] = k;
   645       for (ArcIt a(_graph); a != INVALID; ++a) {
   646         int fi = _arc_idf[a];
   647         int bi = _arc_idb[a];
   648         _reverse[fi] = bi;
   649         _reverse[bi] = fi;
   650       }
   651 
   652       // Reset parameters
   653       resetParams();
   654       return *this;
   655     }
   656 
   657     /// @}
   658 
   659     /// \name Query Functions
   660     /// The results of the algorithm can be obtained using these
   661     /// functions.\n
   662     /// The \ref run() function must be called before using them.
   663 
   664     /// @{
   665 
   666     /// \brief Return the total cost of the found flow.
   667     ///
   668     /// This function returns the total cost of the found flow.
   669     /// Its complexity is O(e).
   670     ///
   671     /// \note The return type of the function can be specified as a
   672     /// template parameter. For example,
   673     /// \code
   674     ///   cs.totalCost<double>();
   675     /// \endcode
   676     /// It is useful if the total cost cannot be stored in the \c Cost
   677     /// type of the algorithm, which is the default return type of the
   678     /// function.
   679     ///
   680     /// \pre \ref run() must be called before using this function.
   681     template <typename Number>
   682     Number totalCost() const {
   683       Number c = 0;
   684       for (ArcIt a(_graph); a != INVALID; ++a) {
   685         int i = _arc_idb[a];
   686         c += static_cast<Number>(_res_cap[i]) *
   687              (-static_cast<Number>(_scost[i]));
   688       }
   689       return c;
   690     }
   691 
   692 #ifndef DOXYGEN
   693     Cost totalCost() const {
   694       return totalCost<Cost>();
   695     }
   696 #endif
   697 
   698     /// \brief Return the flow on the given arc.
   699     ///
   700     /// This function returns the flow on the given arc.
   701     ///
   702     /// \pre \ref run() must be called before using this function.
   703     Value flow(const Arc& a) const {
   704       return _res_cap[_arc_idb[a]];
   705     }
   706 
   707     /// \brief Return the flow map (the primal solution).
   708     ///
   709     /// This function copies the flow value on each arc into the given
   710     /// map. The \c Value type of the algorithm must be convertible to
   711     /// the \c Value type of the map.
   712     ///
   713     /// \pre \ref run() must be called before using this function.
   714     template <typename FlowMap>
   715     void flowMap(FlowMap &map) const {
   716       for (ArcIt a(_graph); a != INVALID; ++a) {
   717         map.set(a, _res_cap[_arc_idb[a]]);
   718       }
   719     }
   720 
   721     /// \brief Return the potential (dual value) of the given node.
   722     ///
   723     /// This function returns the potential (dual value) of the
   724     /// given node.
   725     ///
   726     /// \pre \ref run() must be called before using this function.
   727     Cost potential(const Node& n) const {
   728       return static_cast<Cost>(_pi[_node_id[n]]);
   729     }
   730 
   731     /// \brief Return the potential map (the dual solution).
   732     ///
   733     /// This function copies the potential (dual value) of each node
   734     /// into the given map.
   735     /// The \c Cost type of the algorithm must be convertible to the
   736     /// \c Value type of the map.
   737     ///
   738     /// \pre \ref run() must be called before using this function.
   739     template <typename PotentialMap>
   740     void potentialMap(PotentialMap &map) const {
   741       for (NodeIt n(_graph); n != INVALID; ++n) {
   742         map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
   743       }
   744     }
   745 
   746     /// @}
   747 
   748   private:
   749 
   750     // Initialize the algorithm
   751     ProblemType init() {
   752       if (_res_node_num <= 1) return INFEASIBLE;
   753 
   754       // Check the sum of supply values
   755       _sum_supply = 0;
   756       for (int i = 0; i != _root; ++i) {
   757         _sum_supply += _supply[i];
   758       }
   759       if (_sum_supply > 0) return INFEASIBLE;
   760 
   761 
   762       // Initialize vectors
   763       for (int i = 0; i != _res_node_num; ++i) {
   764         _pi[i] = 0;
   765         _excess[i] = _supply[i];
   766       }
   767 
   768       // Remove infinite upper bounds and check negative arcs
   769       const Value MAX = std::numeric_limits<Value>::max();
   770       int last_out;
   771       if (_have_lower) {
   772         for (int i = 0; i != _root; ++i) {
   773           last_out = _first_out[i+1];
   774           for (int j = _first_out[i]; j != last_out; ++j) {
   775             if (_forward[j]) {
   776               Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
   777               if (c >= MAX) return UNBOUNDED;
   778               _excess[i] -= c;
   779               _excess[_target[j]] += c;
   780             }
   781           }
   782         }
   783       } else {
   784         for (int i = 0; i != _root; ++i) {
   785           last_out = _first_out[i+1];
   786           for (int j = _first_out[i]; j != last_out; ++j) {
   787             if (_forward[j] && _scost[j] < 0) {
   788               Value c = _upper[j];
   789               if (c >= MAX) return UNBOUNDED;
   790               _excess[i] -= c;
   791               _excess[_target[j]] += c;
   792             }
   793           }
   794         }
   795       }
   796       Value ex, max_cap = 0;
   797       for (int i = 0; i != _res_node_num; ++i) {
   798         ex = _excess[i];
   799         _excess[i] = 0;
   800         if (ex < 0) max_cap -= ex;
   801       }
   802       for (int j = 0; j != _res_arc_num; ++j) {
   803         if (_upper[j] >= MAX) _upper[j] = max_cap;
   804       }
   805 
   806       // Initialize the large cost vector and the epsilon parameter
   807       _epsilon = 0;
   808       LargeCost lc;
   809       for (int i = 0; i != _root; ++i) {
   810         last_out = _first_out[i+1];
   811         for (int j = _first_out[i]; j != last_out; ++j) {
   812           lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
   813           _cost[j] = lc;
   814           if (lc > _epsilon) _epsilon = lc;
   815         }
   816       }
   817       _epsilon /= _alpha;
   818 
   819       // Initialize maps for Circulation and remove non-zero lower bounds
   820       ConstMap<Arc, Value> low(0);
   821       typedef typename Digraph::template ArcMap<Value> ValueArcMap;
   822       typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
   823       ValueArcMap cap(_graph), flow(_graph);
   824       ValueNodeMap sup(_graph);
   825       for (NodeIt n(_graph); n != INVALID; ++n) {
   826         sup[n] = _supply[_node_id[n]];
   827       }
   828       if (_have_lower) {
   829         for (ArcIt a(_graph); a != INVALID; ++a) {
   830           int j = _arc_idf[a];
   831           Value c = _lower[j];
   832           cap[a] = _upper[j] - c;
   833           sup[_graph.source(a)] -= c;
   834           sup[_graph.target(a)] += c;
   835         }
   836       } else {
   837         for (ArcIt a(_graph); a != INVALID; ++a) {
   838           cap[a] = _upper[_arc_idf[a]];
   839         }
   840       }
   841 
   842       _sup_node_num = 0;
   843       for (NodeIt n(_graph); n != INVALID; ++n) {
   844         if (sup[n] > 0) ++_sup_node_num;
   845       }
   846 
   847       // Find a feasible flow using Circulation
   848       Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
   849         circ(_graph, low, cap, sup);
   850       if (!circ.flowMap(flow).run()) return INFEASIBLE;
   851 
   852       // Set residual capacities and handle GEQ supply type
   853       if (_sum_supply < 0) {
   854         for (ArcIt a(_graph); a != INVALID; ++a) {
   855           Value fa = flow[a];
   856           _res_cap[_arc_idf[a]] = cap[a] - fa;
   857           _res_cap[_arc_idb[a]] = fa;
   858           sup[_graph.source(a)] -= fa;
   859           sup[_graph.target(a)] += fa;
   860         }
   861         for (NodeIt n(_graph); n != INVALID; ++n) {
   862           _excess[_node_id[n]] = sup[n];
   863         }
   864         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   865           int u = _target[a];
   866           int ra = _reverse[a];
   867           _res_cap[a] = -_sum_supply + 1;
   868           _res_cap[ra] = -_excess[u];
   869           _cost[a] = 0;
   870           _cost[ra] = 0;
   871           _excess[u] = 0;
   872         }
   873       } else {
   874         for (ArcIt a(_graph); a != INVALID; ++a) {
   875           Value fa = flow[a];
   876           _res_cap[_arc_idf[a]] = cap[a] - fa;
   877           _res_cap[_arc_idb[a]] = fa;
   878         }
   879         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   880           int ra = _reverse[a];
   881           _res_cap[a] = 0;
   882           _res_cap[ra] = 0;
   883           _cost[a] = 0;
   884           _cost[ra] = 0;
   885         }
   886       }
   887 
   888       // Initialize data structures for buckets
   889       _max_rank = _alpha * _res_node_num;
   890       _buckets.resize(_max_rank);
   891       _bucket_next.resize(_res_node_num + 1);
   892       _bucket_prev.resize(_res_node_num + 1);
   893       _rank.resize(_res_node_num + 1);
   894 
   895       return OPTIMAL;
   896     }
   897 
   898     // Execute the algorithm and transform the results
   899     void start(Method method) {
   900       const int MAX_PARTIAL_PATH_LENGTH = 4;
   901 
   902       switch (method) {
   903         case PUSH:
   904           startPush();
   905           break;
   906         case AUGMENT:
   907           startAugment(_res_node_num - 1);
   908           break;
   909         case PARTIAL_AUGMENT:
   910           startAugment(MAX_PARTIAL_PATH_LENGTH);
   911           break;
   912       }
   913 
   914       // Compute node potentials (dual solution)
   915       for (int i = 0; i != _res_node_num; ++i) {
   916         _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha));
   917       }
   918       bool optimal = true;
   919       for (int i = 0; optimal && i != _res_node_num; ++i) {
   920         LargeCost pi_i = _pi[i];
   921         int last_out = _first_out[i+1];
   922         for (int j = _first_out[i]; j != last_out; ++j) {
   923           if (_res_cap[j] > 0 && _scost[j] + pi_i - _pi[_target[j]] < 0) {
   924             optimal = false;
   925             break;
   926           }
   927         }
   928       }
   929 
   930       if (!optimal) {
   931         // Compute node potentials for the original costs with BellmanFord
   932         // (if it is necessary)
   933         typedef std::pair<int, int> IntPair;
   934         StaticDigraph sgr;
   935         std::vector<IntPair> arc_vec;
   936         std::vector<LargeCost> cost_vec;
   937         LargeCostArcMap cost_map(cost_vec);
   938 
   939         arc_vec.clear();
   940         cost_vec.clear();
   941         for (int j = 0; j != _res_arc_num; ++j) {
   942           if (_res_cap[j] > 0) {
   943             int u = _source[j], v = _target[j];
   944             arc_vec.push_back(IntPair(u, v));
   945             cost_vec.push_back(_scost[j] + _pi[u] - _pi[v]);
   946           }
   947         }
   948         sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end());
   949 
   950         typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create
   951           bf(sgr, cost_map);
   952         bf.init(0);
   953         bf.start();
   954 
   955         for (int i = 0; i != _res_node_num; ++i) {
   956           _pi[i] += bf.dist(sgr.node(i));
   957         }
   958       }
   959 
   960       // Shift potentials to meet the requirements of the GEQ type
   961       // optimality conditions
   962       LargeCost max_pot = _pi[_root];
   963       for (int i = 0; i != _res_node_num; ++i) {
   964         if (_pi[i] > max_pot) max_pot = _pi[i];
   965       }
   966       if (max_pot != 0) {
   967         for (int i = 0; i != _res_node_num; ++i) {
   968           _pi[i] -= max_pot;
   969         }
   970       }
   971 
   972       // Handle non-zero lower bounds
   973       if (_have_lower) {
   974         int limit = _first_out[_root];
   975         for (int j = 0; j != limit; ++j) {
   976           if (!_forward[j]) _res_cap[j] += _lower[j];
   977         }
   978       }
   979     }
   980 
   981     // Initialize a cost scaling phase
   982     void initPhase() {
   983       // Saturate arcs not satisfying the optimality condition
   984       for (int u = 0; u != _res_node_num; ++u) {
   985         int last_out = _first_out[u+1];
   986         LargeCost pi_u = _pi[u];
   987         for (int a = _first_out[u]; a != last_out; ++a) {
   988           Value delta = _res_cap[a];
   989           if (delta > 0) {
   990             int v = _target[a];
   991             if (_cost[a] + pi_u - _pi[v] < 0) {
   992               _excess[u] -= delta;
   993               _excess[v] += delta;
   994               _res_cap[a] = 0;
   995               _res_cap[_reverse[a]] += delta;
   996             }
   997           }
   998         }
   999       }
  1000 
  1001       // Find active nodes (i.e. nodes with positive excess)
  1002       for (int u = 0; u != _res_node_num; ++u) {
  1003         if (_excess[u] > 0) _active_nodes.push_back(u);
  1004       }
  1005 
  1006       // Initialize the next arcs
  1007       for (int u = 0; u != _res_node_num; ++u) {
  1008         _next_out[u] = _first_out[u];
  1009       }
  1010     }
  1011 
  1012     // Price (potential) refinement heuristic
  1013     bool priceRefinement() {
  1014 
  1015       // Stack for stroing the topological order
  1016       IntVector stack(_res_node_num);
  1017       int stack_top;
  1018 
  1019       // Perform phases
  1020       while (topologicalSort(stack, stack_top)) {
  1021 
  1022         // Compute node ranks in the acyclic admissible network and
  1023         // store the nodes in buckets
  1024         for (int i = 0; i != _res_node_num; ++i) {
  1025           _rank[i] = 0;
  1026         }
  1027         const int bucket_end = _root + 1;
  1028         for (int r = 0; r != _max_rank; ++r) {
  1029           _buckets[r] = bucket_end;
  1030         }
  1031         int top_rank = 0;
  1032         for ( ; stack_top >= 0; --stack_top) {
  1033           int u = stack[stack_top], v;
  1034           int rank_u = _rank[u];
  1035 
  1036           LargeCost rc, pi_u = _pi[u];
  1037           int last_out = _first_out[u+1];
  1038           for (int a = _first_out[u]; a != last_out; ++a) {
  1039             if (_res_cap[a] > 0) {
  1040               v = _target[a];
  1041               rc = _cost[a] + pi_u - _pi[v];
  1042               if (rc < 0) {
  1043                 LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
  1044                 if (nrc < LargeCost(_max_rank)) {
  1045                   int new_rank_v = rank_u + static_cast<int>(nrc);
  1046                   if (new_rank_v > _rank[v]) {
  1047                     _rank[v] = new_rank_v;
  1048                   }
  1049                 }
  1050               }
  1051             }
  1052           }
  1053 
  1054           if (rank_u > 0) {
  1055             top_rank = std::max(top_rank, rank_u);
  1056             int bfirst = _buckets[rank_u];
  1057             _bucket_next[u] = bfirst;
  1058             _bucket_prev[bfirst] = u;
  1059             _buckets[rank_u] = u;
  1060           }
  1061         }
  1062 
  1063         // Check if the current flow is epsilon-optimal
  1064         if (top_rank == 0) {
  1065           return true;
  1066         }
  1067 
  1068         // Process buckets in top-down order
  1069         for (int rank = top_rank; rank > 0; --rank) {
  1070           while (_buckets[rank] != bucket_end) {
  1071             // Remove the first node from the current bucket
  1072             int u = _buckets[rank];
  1073             _buckets[rank] = _bucket_next[u];
  1074 
  1075             // Search the outgoing arcs of u
  1076             LargeCost rc, pi_u = _pi[u];
  1077             int last_out = _first_out[u+1];
  1078             int v, old_rank_v, new_rank_v;
  1079             for (int a = _first_out[u]; a != last_out; ++a) {
  1080               if (_res_cap[a] > 0) {
  1081                 v = _target[a];
  1082                 old_rank_v = _rank[v];
  1083 
  1084                 if (old_rank_v < rank) {
  1085 
  1086                   // Compute the new rank of node v
  1087                   rc = _cost[a] + pi_u - _pi[v];
  1088                   if (rc < 0) {
  1089                     new_rank_v = rank;
  1090                   } else {
  1091                     LargeCost nrc = rc / _epsilon;
  1092                     new_rank_v = 0;
  1093                     if (nrc < LargeCost(_max_rank)) {
  1094                       new_rank_v = rank - 1 - static_cast<int>(nrc);
  1095                     }
  1096                   }
  1097 
  1098                   // Change the rank of node v
  1099                   if (new_rank_v > old_rank_v) {
  1100                     _rank[v] = new_rank_v;
  1101 
  1102                     // Remove v from its old bucket
  1103                     if (old_rank_v > 0) {
  1104                       if (_buckets[old_rank_v] == v) {
  1105                         _buckets[old_rank_v] = _bucket_next[v];
  1106                       } else {
  1107                         int pv = _bucket_prev[v], nv = _bucket_next[v];
  1108                         _bucket_next[pv] = nv;
  1109                         _bucket_prev[nv] = pv;
  1110                       }
  1111                     }
  1112 
  1113                     // Insert v into its new bucket
  1114                     int nv = _buckets[new_rank_v];
  1115                     _bucket_next[v] = nv;
  1116                     _bucket_prev[nv] = v;
  1117                     _buckets[new_rank_v] = v;
  1118                   }
  1119                 }
  1120               }
  1121             }
  1122 
  1123             // Refine potential of node u
  1124             _pi[u] -= rank * _epsilon;
  1125           }
  1126         }
  1127 
  1128       }
  1129 
  1130       return false;
  1131     }
  1132 
  1133     // Find and cancel cycles in the admissible network and
  1134     // determine topological order using DFS
  1135     bool topologicalSort(IntVector &stack, int &stack_top) {
  1136       const int MAX_CYCLE_CANCEL = 1;
  1137 
  1138       BoolVector reached(_res_node_num, false);
  1139       BoolVector processed(_res_node_num, false);
  1140       IntVector pred(_res_node_num);
  1141       for (int i = 0; i != _res_node_num; ++i) {
  1142         _next_out[i] = _first_out[i];
  1143       }
  1144       stack_top = -1;
  1145 
  1146       int cycle_cnt = 0;
  1147       for (int start = 0; start != _res_node_num; ++start) {
  1148         if (reached[start]) continue;
  1149 
  1150         // Start DFS search from this start node
  1151         pred[start] = -1;
  1152         int tip = start, v;
  1153         while (true) {
  1154           // Check the outgoing arcs of the current tip node
  1155           reached[tip] = true;
  1156           LargeCost pi_tip = _pi[tip];
  1157           int a, last_out = _first_out[tip+1];
  1158           for (a = _next_out[tip]; a != last_out; ++a) {
  1159             if (_res_cap[a] > 0) {
  1160               v = _target[a];
  1161               if (_cost[a] + pi_tip - _pi[v] < 0) {
  1162                 if (!reached[v]) {
  1163                   // A new node is reached
  1164                   reached[v] = true;
  1165                   pred[v] = tip;
  1166                   _next_out[tip] = a;
  1167                   tip = v;
  1168                   a = _next_out[tip];
  1169                   last_out = _first_out[tip+1];
  1170                   break;
  1171                 }
  1172                 else if (!processed[v]) {
  1173                   // A cycle is found
  1174                   ++cycle_cnt;
  1175                   _next_out[tip] = a;
  1176 
  1177                   // Find the minimum residual capacity along the cycle
  1178                   Value d, delta = _res_cap[a];
  1179                   int u, delta_node = tip;
  1180                   for (u = tip; u != v; ) {
  1181                     u = pred[u];
  1182                     d = _res_cap[_next_out[u]];
  1183                     if (d <= delta) {
  1184                       delta = d;
  1185                       delta_node = u;
  1186                     }
  1187                   }
  1188 
  1189                   // Augment along the cycle
  1190                   _res_cap[a] -= delta;
  1191                   _res_cap[_reverse[a]] += delta;
  1192                   for (u = tip; u != v; ) {
  1193                     u = pred[u];
  1194                     int ca = _next_out[u];
  1195                     _res_cap[ca] -= delta;
  1196                     _res_cap[_reverse[ca]] += delta;
  1197                   }
  1198 
  1199                   // Check the maximum number of cycle canceling
  1200                   if (cycle_cnt >= MAX_CYCLE_CANCEL) {
  1201                     return false;
  1202                   }
  1203 
  1204                   // Roll back search to delta_node
  1205                   if (delta_node != tip) {
  1206                     for (u = tip; u != delta_node; u = pred[u]) {
  1207                       reached[u] = false;
  1208                     }
  1209                     tip = delta_node;
  1210                     a = _next_out[tip] + 1;
  1211                     last_out = _first_out[tip+1];
  1212                     break;
  1213                   }
  1214                 }
  1215               }
  1216             }
  1217           }
  1218 
  1219           // Step back to the previous node
  1220           if (a == last_out) {
  1221             processed[tip] = true;
  1222             stack[++stack_top] = tip;
  1223             tip = pred[tip];
  1224             if (tip < 0) {
  1225               // Finish DFS from the current start node
  1226               break;
  1227             }
  1228             ++_next_out[tip];
  1229           }
  1230         }
  1231 
  1232       }
  1233 
  1234       return (cycle_cnt == 0);
  1235     }
  1236 
  1237     // Global potential update heuristic
  1238     void globalUpdate() {
  1239       const int bucket_end = _root + 1;
  1240 
  1241       // Initialize buckets
  1242       for (int r = 0; r != _max_rank; ++r) {
  1243         _buckets[r] = bucket_end;
  1244       }
  1245       Value total_excess = 0;
  1246       int b0 = bucket_end;
  1247       for (int i = 0; i != _res_node_num; ++i) {
  1248         if (_excess[i] < 0) {
  1249           _rank[i] = 0;
  1250           _bucket_next[i] = b0;
  1251           _bucket_prev[b0] = i;
  1252           b0 = i;
  1253         } else {
  1254           total_excess += _excess[i];
  1255           _rank[i] = _max_rank;
  1256         }
  1257       }
  1258       if (total_excess == 0) return;
  1259       _buckets[0] = b0;
  1260 
  1261       // Search the buckets
  1262       int r = 0;
  1263       for ( ; r != _max_rank; ++r) {
  1264         while (_buckets[r] != bucket_end) {
  1265           // Remove the first node from the current bucket
  1266           int u = _buckets[r];
  1267           _buckets[r] = _bucket_next[u];
  1268 
  1269           // Search the incomming arcs of u
  1270           LargeCost pi_u = _pi[u];
  1271           int last_out = _first_out[u+1];
  1272           for (int a = _first_out[u]; a != last_out; ++a) {
  1273             int ra = _reverse[a];
  1274             if (_res_cap[ra] > 0) {
  1275               int v = _source[ra];
  1276               int old_rank_v = _rank[v];
  1277               if (r < old_rank_v) {
  1278                 // Compute the new rank of v
  1279                 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
  1280                 int new_rank_v = old_rank_v;
  1281                 if (nrc < LargeCost(_max_rank)) {
  1282                   new_rank_v = r + 1 + static_cast<int>(nrc);
  1283                 }
  1284 
  1285                 // Change the rank of v
  1286                 if (new_rank_v < old_rank_v) {
  1287                   _rank[v] = new_rank_v;
  1288                   _next_out[v] = _first_out[v];
  1289 
  1290                   // Remove v from its old bucket
  1291                   if (old_rank_v < _max_rank) {
  1292                     if (_buckets[old_rank_v] == v) {
  1293                       _buckets[old_rank_v] = _bucket_next[v];
  1294                     } else {
  1295                       int pv = _bucket_prev[v], nv = _bucket_next[v];
  1296                       _bucket_next[pv] = nv;
  1297                       _bucket_prev[nv] = pv;
  1298                     }
  1299                   }
  1300 
  1301                   // Insert v into its new bucket
  1302                   int nv = _buckets[new_rank_v];
  1303                   _bucket_next[v] = nv;
  1304                   _bucket_prev[nv] = v;
  1305                   _buckets[new_rank_v] = v;
  1306                 }
  1307               }
  1308             }
  1309           }
  1310 
  1311           // Finish search if there are no more active nodes
  1312           if (_excess[u] > 0) {
  1313             total_excess -= _excess[u];
  1314             if (total_excess <= 0) break;
  1315           }
  1316         }
  1317         if (total_excess <= 0) break;
  1318       }
  1319 
  1320       // Relabel nodes
  1321       for (int u = 0; u != _res_node_num; ++u) {
  1322         int k = std::min(_rank[u], r);
  1323         if (k > 0) {
  1324           _pi[u] -= _epsilon * k;
  1325           _next_out[u] = _first_out[u];
  1326         }
  1327       }
  1328     }
  1329 
  1330     /// Execute the algorithm performing augment and relabel operations
  1331     void startAugment(int max_length) {
  1332       // Paramters for heuristics
  1333       const int PRICE_REFINEMENT_LIMIT = 2;
  1334       const double GLOBAL_UPDATE_FACTOR = 1.0;
  1335       const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
  1336         (_res_node_num + _sup_node_num * _sup_node_num));
  1337       int next_global_update_limit = global_update_skip;
  1338 
  1339       // Perform cost scaling phases
  1340       IntVector path;
  1341       BoolVector path_arc(_res_arc_num, false);
  1342       int relabel_cnt = 0;
  1343       int eps_phase_cnt = 0;
  1344       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1345                                         1 : _epsilon / _alpha )
  1346       {
  1347         ++eps_phase_cnt;
  1348 
  1349         // Price refinement heuristic
  1350         if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
  1351           if (priceRefinement()) continue;
  1352         }
  1353 
  1354         // Initialize current phase
  1355         initPhase();
  1356 
  1357         // Perform partial augment and relabel operations
  1358         while (true) {
  1359           // Select an active node (FIFO selection)
  1360           while (_active_nodes.size() > 0 &&
  1361                  _excess[_active_nodes.front()] <= 0) {
  1362             _active_nodes.pop_front();
  1363           }
  1364           if (_active_nodes.size() == 0) break;
  1365           int start = _active_nodes.front();
  1366 
  1367           // Find an augmenting path from the start node
  1368           int tip = start;
  1369           while (int(path.size()) < max_length && _excess[tip] >= 0) {
  1370             int u;
  1371             LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
  1372             LargeCost pi_tip = _pi[tip];
  1373             int last_out = _first_out[tip+1];
  1374             for (int a = _next_out[tip]; a != last_out; ++a) {
  1375               if (_res_cap[a] > 0) {
  1376                 u = _target[a];
  1377                 rc = _cost[a] + pi_tip - _pi[u];
  1378                 if (rc < 0) {
  1379                   path.push_back(a);
  1380                   _next_out[tip] = a;
  1381                   if (path_arc[a]) {
  1382                     goto augment;   // a cycle is found, stop path search
  1383                   }
  1384                   tip = u;
  1385                   path_arc[a] = true;
  1386                   goto next_step;
  1387                 }
  1388                 else if (rc < min_red_cost) {
  1389                   min_red_cost = rc;
  1390                 }
  1391               }
  1392             }
  1393 
  1394             // Relabel tip node
  1395             if (tip != start) {
  1396               int ra = _reverse[path.back()];
  1397               min_red_cost =
  1398                 std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
  1399             }
  1400             last_out = _next_out[tip];
  1401             for (int a = _first_out[tip]; a != last_out; ++a) {
  1402               if (_res_cap[a] > 0) {
  1403                 rc = _cost[a] + pi_tip - _pi[_target[a]];
  1404                 if (rc < min_red_cost) {
  1405                   min_red_cost = rc;
  1406                 }
  1407               }
  1408             }
  1409             _pi[tip] -= min_red_cost + _epsilon;
  1410             _next_out[tip] = _first_out[tip];
  1411             ++relabel_cnt;
  1412 
  1413             // Step back
  1414             if (tip != start) {
  1415               int pa = path.back();
  1416               path_arc[pa] = false;
  1417               tip = _source[pa];
  1418               path.pop_back();
  1419             }
  1420 
  1421           next_step: ;
  1422           }
  1423 
  1424           // Augment along the found path (as much flow as possible)
  1425         augment:
  1426           Value delta;
  1427           int pa, u, v = start;
  1428           for (int i = 0; i != int(path.size()); ++i) {
  1429             pa = path[i];
  1430             u = v;
  1431             v = _target[pa];
  1432             path_arc[pa] = false;
  1433             delta = std::min(_res_cap[pa], _excess[u]);
  1434             _res_cap[pa] -= delta;
  1435             _res_cap[_reverse[pa]] += delta;
  1436             _excess[u] -= delta;
  1437             _excess[v] += delta;
  1438             if (_excess[v] > 0 && _excess[v] <= delta) {
  1439               _active_nodes.push_back(v);
  1440             }
  1441           }
  1442           path.clear();
  1443 
  1444           // Global update heuristic
  1445           if (relabel_cnt >= next_global_update_limit) {
  1446             globalUpdate();
  1447             next_global_update_limit += global_update_skip;
  1448           }
  1449         }
  1450 
  1451       }
  1452 
  1453     }
  1454 
  1455     /// Execute the algorithm performing push and relabel operations
  1456     void startPush() {
  1457       // Paramters for heuristics
  1458       const int PRICE_REFINEMENT_LIMIT = 2;
  1459       const double GLOBAL_UPDATE_FACTOR = 2.0;
  1460 
  1461       const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
  1462         (_res_node_num + _sup_node_num * _sup_node_num));
  1463       int next_global_update_limit = global_update_skip;
  1464 
  1465       // Perform cost scaling phases
  1466       BoolVector hyper(_res_node_num, false);
  1467       LargeCostVector hyper_cost(_res_node_num);
  1468       int relabel_cnt = 0;
  1469       int eps_phase_cnt = 0;
  1470       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1471                                         1 : _epsilon / _alpha )
  1472       {
  1473         ++eps_phase_cnt;
  1474 
  1475         // Price refinement heuristic
  1476         if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
  1477           if (priceRefinement()) continue;
  1478         }
  1479 
  1480         // Initialize current phase
  1481         initPhase();
  1482 
  1483         // Perform push and relabel operations
  1484         while (_active_nodes.size() > 0) {
  1485           LargeCost min_red_cost, rc, pi_n;
  1486           Value delta;
  1487           int n, t, a, last_out = _res_arc_num;
  1488 
  1489         next_node:
  1490           // Select an active node (FIFO selection)
  1491           n = _active_nodes.front();
  1492           last_out = _first_out[n+1];
  1493           pi_n = _pi[n];
  1494 
  1495           // Perform push operations if there are admissible arcs
  1496           if (_excess[n] > 0) {
  1497             for (a = _next_out[n]; a != last_out; ++a) {
  1498               if (_res_cap[a] > 0 &&
  1499                   _cost[a] + pi_n - _pi[_target[a]] < 0) {
  1500                 delta = std::min(_res_cap[a], _excess[n]);
  1501                 t = _target[a];
  1502 
  1503                 // Push-look-ahead heuristic
  1504                 Value ahead = -_excess[t];
  1505                 int last_out_t = _first_out[t+1];
  1506                 LargeCost pi_t = _pi[t];
  1507                 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
  1508                   if (_res_cap[ta] > 0 &&
  1509                       _cost[ta] + pi_t - _pi[_target[ta]] < 0)
  1510                     ahead += _res_cap[ta];
  1511                   if (ahead >= delta) break;
  1512                 }
  1513                 if (ahead < 0) ahead = 0;
  1514 
  1515                 // Push flow along the arc
  1516                 if (ahead < delta && !hyper[t]) {
  1517                   _res_cap[a] -= ahead;
  1518                   _res_cap[_reverse[a]] += ahead;
  1519                   _excess[n] -= ahead;
  1520                   _excess[t] += ahead;
  1521                   _active_nodes.push_front(t);
  1522                   hyper[t] = true;
  1523                   hyper_cost[t] = _cost[a] + pi_n - pi_t;
  1524                   _next_out[n] = a;
  1525                   goto next_node;
  1526                 } else {
  1527                   _res_cap[a] -= delta;
  1528                   _res_cap[_reverse[a]] += delta;
  1529                   _excess[n] -= delta;
  1530                   _excess[t] += delta;
  1531                   if (_excess[t] > 0 && _excess[t] <= delta)
  1532                     _active_nodes.push_back(t);
  1533                 }
  1534 
  1535                 if (_excess[n] == 0) {
  1536                   _next_out[n] = a;
  1537                   goto remove_nodes;
  1538                 }
  1539               }
  1540             }
  1541             _next_out[n] = a;
  1542           }
  1543 
  1544           // Relabel the node if it is still active (or hyper)
  1545           if (_excess[n] > 0 || hyper[n]) {
  1546              min_red_cost = hyper[n] ? -hyper_cost[n] :
  1547                std::numeric_limits<LargeCost>::max();
  1548             for (int a = _first_out[n]; a != last_out; ++a) {
  1549               if (_res_cap[a] > 0) {
  1550                 rc = _cost[a] + pi_n - _pi[_target[a]];
  1551                 if (rc < min_red_cost) {
  1552                   min_red_cost = rc;
  1553                 }
  1554               }
  1555             }
  1556             _pi[n] -= min_red_cost + _epsilon;
  1557             _next_out[n] = _first_out[n];
  1558             hyper[n] = false;
  1559             ++relabel_cnt;
  1560           }
  1561 
  1562           // Remove nodes that are not active nor hyper
  1563         remove_nodes:
  1564           while ( _active_nodes.size() > 0 &&
  1565                   _excess[_active_nodes.front()] <= 0 &&
  1566                   !hyper[_active_nodes.front()] ) {
  1567             _active_nodes.pop_front();
  1568           }
  1569 
  1570           // Global update heuristic
  1571           if (relabel_cnt >= next_global_update_limit) {
  1572             globalUpdate();
  1573             for (int u = 0; u != _res_node_num; ++u)
  1574               hyper[u] = false;
  1575             next_global_update_limit += global_update_skip;
  1576           }
  1577         }
  1578       }
  1579     }
  1580 
  1581   }; //class CostScaling
  1582 
  1583   ///@}
  1584 
  1585 } //namespace lemon
  1586 
  1587 #endif //LEMON_COST_SCALING_H