lemon/cost_scaling.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sun, 09 Jan 2011 00:56:52 +0100
changeset 1034 ef200e268af2
parent 863 a93f1a27d831
child 919 e0cef67fe565
child 921 140c953ad5d1
child 931 f112c18bc304
permissions -rw-r--r--
Unifications and improvements in TSP algorithms (#386)
     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   /// 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<Value> ValueVector;
   205     typedef std::vector<Cost> CostVector;
   206     typedef std::vector<LargeCost> LargeCostVector;
   207     typedef std::vector<char> BoolVector;
   208     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
   209 
   210   private:
   211 
   212     template <typename KT, typename VT>
   213     class StaticVectorMap {
   214     public:
   215       typedef KT Key;
   216       typedef VT Value;
   217 
   218       StaticVectorMap(std::vector<Value>& v) : _v(v) {}
   219 
   220       const Value& operator[](const Key& key) const {
   221         return _v[StaticDigraph::id(key)];
   222       }
   223 
   224       Value& operator[](const Key& key) {
   225         return _v[StaticDigraph::id(key)];
   226       }
   227 
   228       void set(const Key& key, const Value& val) {
   229         _v[StaticDigraph::id(key)] = val;
   230       }
   231 
   232     private:
   233       std::vector<Value>& _v;
   234     };
   235 
   236     typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
   237     typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
   238 
   239   private:
   240 
   241     // Data related to the underlying digraph
   242     const GR &_graph;
   243     int _node_num;
   244     int _arc_num;
   245     int _res_node_num;
   246     int _res_arc_num;
   247     int _root;
   248 
   249     // Parameters of the problem
   250     bool _have_lower;
   251     Value _sum_supply;
   252     int _sup_node_num;
   253 
   254     // Data structures for storing the digraph
   255     IntNodeMap _node_id;
   256     IntArcMap _arc_idf;
   257     IntArcMap _arc_idb;
   258     IntVector _first_out;
   259     BoolVector _forward;
   260     IntVector _source;
   261     IntVector _target;
   262     IntVector _reverse;
   263 
   264     // Node and arc data
   265     ValueVector _lower;
   266     ValueVector _upper;
   267     CostVector _scost;
   268     ValueVector _supply;
   269 
   270     ValueVector _res_cap;
   271     LargeCostVector _cost;
   272     LargeCostVector _pi;
   273     ValueVector _excess;
   274     IntVector _next_out;
   275     std::deque<int> _active_nodes;
   276 
   277     // Data for scaling
   278     LargeCost _epsilon;
   279     int _alpha;
   280 
   281     IntVector _buckets;
   282     IntVector _bucket_next;
   283     IntVector _bucket_prev;
   284     IntVector _rank;
   285     int _max_rank;
   286 
   287     // Data for a StaticDigraph structure
   288     typedef std::pair<int, int> IntPair;
   289     StaticDigraph _sgr;
   290     std::vector<IntPair> _arc_vec;
   291     std::vector<LargeCost> _cost_vec;
   292     LargeCostArcMap _cost_map;
   293     LargeCostNodeMap _pi_map;
   294 
   295   public:
   296 
   297     /// \brief Constant for infinite upper bounds (capacities).
   298     ///
   299     /// Constant for infinite upper bounds (capacities).
   300     /// It is \c std::numeric_limits<Value>::infinity() if available,
   301     /// \c std::numeric_limits<Value>::max() otherwise.
   302     const Value INF;
   303 
   304   public:
   305 
   306     /// \name Named Template Parameters
   307     /// @{
   308 
   309     template <typename T>
   310     struct SetLargeCostTraits : public Traits {
   311       typedef T LargeCost;
   312     };
   313 
   314     /// \brief \ref named-templ-param "Named parameter" for setting
   315     /// \c LargeCost type.
   316     ///
   317     /// \ref named-templ-param "Named parameter" for setting \c LargeCost
   318     /// type, which is used for internal computations in the algorithm.
   319     /// \c Cost must be convertible to \c LargeCost.
   320     template <typename T>
   321     struct SetLargeCost
   322       : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
   323       typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
   324     };
   325 
   326     /// @}
   327 
   328   protected:
   329 
   330     CostScaling() {}
   331 
   332   public:
   333 
   334     /// \brief Constructor.
   335     ///
   336     /// The constructor of the class.
   337     ///
   338     /// \param graph The digraph the algorithm runs on.
   339     CostScaling(const GR& graph) :
   340       _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
   341       _cost_map(_cost_vec), _pi_map(_pi),
   342       INF(std::numeric_limits<Value>::has_infinity ?
   343           std::numeric_limits<Value>::infinity() :
   344           std::numeric_limits<Value>::max())
   345     {
   346       // Check the number types
   347       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
   348         "The flow type of CostScaling must be signed");
   349       LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   350         "The cost type of CostScaling must be signed");
   351 
   352       // Reset data structures
   353       reset();
   354     }
   355 
   356     /// \name Parameters
   357     /// The parameters of the algorithm can be specified using these
   358     /// functions.
   359 
   360     /// @{
   361 
   362     /// \brief Set the lower bounds on the arcs.
   363     ///
   364     /// This function sets the lower bounds on the arcs.
   365     /// If it is not used before calling \ref run(), the lower bounds
   366     /// will be set to zero on all arcs.
   367     ///
   368     /// \param map An arc map storing the lower bounds.
   369     /// Its \c Value type must be convertible to the \c Value type
   370     /// of the algorithm.
   371     ///
   372     /// \return <tt>(*this)</tt>
   373     template <typename LowerMap>
   374     CostScaling& lowerMap(const LowerMap& map) {
   375       _have_lower = true;
   376       for (ArcIt a(_graph); a != INVALID; ++a) {
   377         _lower[_arc_idf[a]] = map[a];
   378         _lower[_arc_idb[a]] = map[a];
   379       }
   380       return *this;
   381     }
   382 
   383     /// \brief Set the upper bounds (capacities) on the arcs.
   384     ///
   385     /// This function sets the upper bounds (capacities) on the arcs.
   386     /// If it is not used before calling \ref run(), the upper bounds
   387     /// will be set to \ref INF on all arcs (i.e. the flow value will be
   388     /// unbounded from above).
   389     ///
   390     /// \param map An arc map storing the upper bounds.
   391     /// Its \c Value type must be convertible to the \c Value type
   392     /// of the algorithm.
   393     ///
   394     /// \return <tt>(*this)</tt>
   395     template<typename UpperMap>
   396     CostScaling& upperMap(const UpperMap& map) {
   397       for (ArcIt a(_graph); a != INVALID; ++a) {
   398         _upper[_arc_idf[a]] = map[a];
   399       }
   400       return *this;
   401     }
   402 
   403     /// \brief Set the costs of the arcs.
   404     ///
   405     /// This function sets the costs of the arcs.
   406     /// If it is not used before calling \ref run(), the costs
   407     /// will be set to \c 1 on all arcs.
   408     ///
   409     /// \param map An arc map storing the costs.
   410     /// Its \c Value type must be convertible to the \c Cost type
   411     /// of the algorithm.
   412     ///
   413     /// \return <tt>(*this)</tt>
   414     template<typename CostMap>
   415     CostScaling& costMap(const CostMap& map) {
   416       for (ArcIt a(_graph); a != INVALID; ++a) {
   417         _scost[_arc_idf[a]] =  map[a];
   418         _scost[_arc_idb[a]] = -map[a];
   419       }
   420       return *this;
   421     }
   422 
   423     /// \brief Set the supply values of the nodes.
   424     ///
   425     /// This function sets the supply values of the nodes.
   426     /// If neither this function nor \ref stSupply() is used before
   427     /// calling \ref run(), the supply of each node will be set to zero.
   428     ///
   429     /// \param map A node map storing the supply values.
   430     /// Its \c Value type must be convertible to the \c Value type
   431     /// of the algorithm.
   432     ///
   433     /// \return <tt>(*this)</tt>
   434     template<typename SupplyMap>
   435     CostScaling& supplyMap(const SupplyMap& map) {
   436       for (NodeIt n(_graph); n != INVALID; ++n) {
   437         _supply[_node_id[n]] = map[n];
   438       }
   439       return *this;
   440     }
   441 
   442     /// \brief Set single source and target nodes and a supply value.
   443     ///
   444     /// This function sets a single source node and a single target node
   445     /// and the required flow value.
   446     /// If neither this function nor \ref supplyMap() is used before
   447     /// calling \ref run(), the supply of each node will be set to zero.
   448     ///
   449     /// Using this function has the same effect as using \ref supplyMap()
   450     /// with such a map in which \c k is assigned to \c s, \c -k is
   451     /// assigned to \c t and all other nodes have zero supply value.
   452     ///
   453     /// \param s The source node.
   454     /// \param t The target node.
   455     /// \param k The required amount of flow from node \c s to node \c t
   456     /// (i.e. the supply of \c s and the demand of \c t).
   457     ///
   458     /// \return <tt>(*this)</tt>
   459     CostScaling& stSupply(const Node& s, const Node& t, Value k) {
   460       for (int i = 0; i != _res_node_num; ++i) {
   461         _supply[i] = 0;
   462       }
   463       _supply[_node_id[s]] =  k;
   464       _supply[_node_id[t]] = -k;
   465       return *this;
   466     }
   467 
   468     /// @}
   469 
   470     /// \name Execution control
   471     /// The algorithm can be executed using \ref run().
   472 
   473     /// @{
   474 
   475     /// \brief Run the algorithm.
   476     ///
   477     /// This function runs the algorithm.
   478     /// The paramters can be specified using functions \ref lowerMap(),
   479     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
   480     /// For example,
   481     /// \code
   482     ///   CostScaling<ListDigraph> cs(graph);
   483     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   484     ///     .supplyMap(sup).run();
   485     /// \endcode
   486     ///
   487     /// This function can be called more than once. All the given parameters
   488     /// are kept for the next call, unless \ref resetParams() or \ref reset()
   489     /// is used, thus only the modified parameters have to be set again.
   490     /// If the underlying digraph was also modified after the construction
   491     /// of the class (or the last \ref reset() call), then the \ref reset()
   492     /// function must be called.
   493     ///
   494     /// \param method The internal method that will be used in the
   495     /// algorithm. For more information, see \ref Method.
   496     /// \param factor The cost scaling factor. It must be larger than one.
   497     ///
   498     /// \return \c INFEASIBLE if no feasible flow exists,
   499     /// \n \c OPTIMAL if the problem has optimal solution
   500     /// (i.e. it is feasible and bounded), and the algorithm has found
   501     /// optimal flow and node potentials (primal and dual solutions),
   502     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
   503     /// and infinite upper bound. It means that the objective function
   504     /// is unbounded on that arc, however, note that it could actually be
   505     /// bounded over the feasible flows, but this algroithm cannot handle
   506     /// these cases.
   507     ///
   508     /// \see ProblemType, Method
   509     /// \see resetParams(), reset()
   510     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
   511       _alpha = factor;
   512       ProblemType pt = init();
   513       if (pt != OPTIMAL) return pt;
   514       start(method);
   515       return OPTIMAL;
   516     }
   517 
   518     /// \brief Reset all the parameters that have been given before.
   519     ///
   520     /// This function resets all the paramaters that have been given
   521     /// before using functions \ref lowerMap(), \ref upperMap(),
   522     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   523     ///
   524     /// It is useful for multiple \ref run() calls. Basically, all the given
   525     /// parameters are kept for the next \ref run() call, unless
   526     /// \ref resetParams() or \ref reset() is used.
   527     /// If the underlying digraph was also modified after the construction
   528     /// of the class or the last \ref reset() call, then the \ref reset()
   529     /// function must be used, otherwise \ref resetParams() is sufficient.
   530     ///
   531     /// For example,
   532     /// \code
   533     ///   CostScaling<ListDigraph> cs(graph);
   534     ///
   535     ///   // First run
   536     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
   537     ///     .supplyMap(sup).run();
   538     ///
   539     ///   // Run again with modified cost map (resetParams() is not called,
   540     ///   // so only the cost map have to be set again)
   541     ///   cost[e] += 100;
   542     ///   cs.costMap(cost).run();
   543     ///
   544     ///   // Run again from scratch using resetParams()
   545     ///   // (the lower bounds will be set to zero on all arcs)
   546     ///   cs.resetParams();
   547     ///   cs.upperMap(capacity).costMap(cost)
   548     ///     .supplyMap(sup).run();
   549     /// \endcode
   550     ///
   551     /// \return <tt>(*this)</tt>
   552     ///
   553     /// \see reset(), run()
   554     CostScaling& resetParams() {
   555       for (int i = 0; i != _res_node_num; ++i) {
   556         _supply[i] = 0;
   557       }
   558       int limit = _first_out[_root];
   559       for (int j = 0; j != limit; ++j) {
   560         _lower[j] = 0;
   561         _upper[j] = INF;
   562         _scost[j] = _forward[j] ? 1 : -1;
   563       }
   564       for (int j = limit; j != _res_arc_num; ++j) {
   565         _lower[j] = 0;
   566         _upper[j] = INF;
   567         _scost[j] = 0;
   568         _scost[_reverse[j]] = 0;
   569       }
   570       _have_lower = false;
   571       return *this;
   572     }
   573 
   574     /// \brief Reset all the parameters that have been given before.
   575     ///
   576     /// This function resets all the paramaters that have been given
   577     /// before using functions \ref lowerMap(), \ref upperMap(),
   578     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
   579     ///
   580     /// It is useful for multiple run() calls. If this function is not
   581     /// used, all the parameters given before are kept for the next
   582     /// \ref run() call.
   583     /// However, the underlying digraph must not be modified after this
   584     /// class have been constructed, since it copies and extends the graph.
   585     /// \return <tt>(*this)</tt>
   586     CostScaling& reset() {
   587       // Resize vectors
   588       _node_num = countNodes(_graph);
   589       _arc_num = countArcs(_graph);
   590       _res_node_num = _node_num + 1;
   591       _res_arc_num = 2 * (_arc_num + _node_num);
   592       _root = _node_num;
   593 
   594       _first_out.resize(_res_node_num + 1);
   595       _forward.resize(_res_arc_num);
   596       _source.resize(_res_arc_num);
   597       _target.resize(_res_arc_num);
   598       _reverse.resize(_res_arc_num);
   599 
   600       _lower.resize(_res_arc_num);
   601       _upper.resize(_res_arc_num);
   602       _scost.resize(_res_arc_num);
   603       _supply.resize(_res_node_num);
   604 
   605       _res_cap.resize(_res_arc_num);
   606       _cost.resize(_res_arc_num);
   607       _pi.resize(_res_node_num);
   608       _excess.resize(_res_node_num);
   609       _next_out.resize(_res_node_num);
   610 
   611       _arc_vec.reserve(_res_arc_num);
   612       _cost_vec.reserve(_res_arc_num);
   613 
   614       // Copy the graph
   615       int i = 0, j = 0, k = 2 * _arc_num + _node_num;
   616       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   617         _node_id[n] = i;
   618       }
   619       i = 0;
   620       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   621         _first_out[i] = j;
   622         for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   623           _arc_idf[a] = j;
   624           _forward[j] = true;
   625           _source[j] = i;
   626           _target[j] = _node_id[_graph.runningNode(a)];
   627         }
   628         for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
   629           _arc_idb[a] = j;
   630           _forward[j] = false;
   631           _source[j] = i;
   632           _target[j] = _node_id[_graph.runningNode(a)];
   633         }
   634         _forward[j] = false;
   635         _source[j] = i;
   636         _target[j] = _root;
   637         _reverse[j] = k;
   638         _forward[k] = true;
   639         _source[k] = _root;
   640         _target[k] = i;
   641         _reverse[k] = j;
   642         ++j; ++k;
   643       }
   644       _first_out[i] = j;
   645       _first_out[_res_node_num] = k;
   646       for (ArcIt a(_graph); a != INVALID; ++a) {
   647         int fi = _arc_idf[a];
   648         int bi = _arc_idb[a];
   649         _reverse[fi] = bi;
   650         _reverse[bi] = fi;
   651       }
   652 
   653       // Reset parameters
   654       resetParams();
   655       return *this;
   656     }
   657 
   658     /// @}
   659 
   660     /// \name Query Functions
   661     /// The results of the algorithm can be obtained using these
   662     /// functions.\n
   663     /// The \ref run() function must be called before using them.
   664 
   665     /// @{
   666 
   667     /// \brief Return the total cost of the found flow.
   668     ///
   669     /// This function returns the total cost of the found flow.
   670     /// Its complexity is O(e).
   671     ///
   672     /// \note The return type of the function can be specified as a
   673     /// template parameter. For example,
   674     /// \code
   675     ///   cs.totalCost<double>();
   676     /// \endcode
   677     /// It is useful if the total cost cannot be stored in the \c Cost
   678     /// type of the algorithm, which is the default return type of the
   679     /// function.
   680     ///
   681     /// \pre \ref run() must be called before using this function.
   682     template <typename Number>
   683     Number totalCost() const {
   684       Number c = 0;
   685       for (ArcIt a(_graph); a != INVALID; ++a) {
   686         int i = _arc_idb[a];
   687         c += static_cast<Number>(_res_cap[i]) *
   688              (-static_cast<Number>(_scost[i]));
   689       }
   690       return c;
   691     }
   692 
   693 #ifndef DOXYGEN
   694     Cost totalCost() const {
   695       return totalCost<Cost>();
   696     }
   697 #endif
   698 
   699     /// \brief Return the flow on the given arc.
   700     ///
   701     /// This function returns the flow on the given arc.
   702     ///
   703     /// \pre \ref run() must be called before using this function.
   704     Value flow(const Arc& a) const {
   705       return _res_cap[_arc_idb[a]];
   706     }
   707 
   708     /// \brief Return the flow map (the primal solution).
   709     ///
   710     /// This function copies the flow value on each arc into the given
   711     /// map. The \c Value type of the algorithm must be convertible to
   712     /// the \c Value type of the map.
   713     ///
   714     /// \pre \ref run() must be called before using this function.
   715     template <typename FlowMap>
   716     void flowMap(FlowMap &map) const {
   717       for (ArcIt a(_graph); a != INVALID; ++a) {
   718         map.set(a, _res_cap[_arc_idb[a]]);
   719       }
   720     }
   721 
   722     /// \brief Return the potential (dual value) of the given node.
   723     ///
   724     /// This function returns the potential (dual value) of the
   725     /// given node.
   726     ///
   727     /// \pre \ref run() must be called before using this function.
   728     Cost potential(const Node& n) const {
   729       return static_cast<Cost>(_pi[_node_id[n]]);
   730     }
   731 
   732     /// \brief Return the potential map (the dual solution).
   733     ///
   734     /// This function copies the potential (dual value) of each node
   735     /// into the given map.
   736     /// The \c Cost type of the algorithm must be convertible to the
   737     /// \c Value type of the map.
   738     ///
   739     /// \pre \ref run() must be called before using this function.
   740     template <typename PotentialMap>
   741     void potentialMap(PotentialMap &map) const {
   742       for (NodeIt n(_graph); n != INVALID; ++n) {
   743         map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
   744       }
   745     }
   746 
   747     /// @}
   748 
   749   private:
   750 
   751     // Initialize the algorithm
   752     ProblemType init() {
   753       if (_res_node_num <= 1) return INFEASIBLE;
   754 
   755       // Check the sum of supply values
   756       _sum_supply = 0;
   757       for (int i = 0; i != _root; ++i) {
   758         _sum_supply += _supply[i];
   759       }
   760       if (_sum_supply > 0) return INFEASIBLE;
   761 
   762 
   763       // Initialize vectors
   764       for (int i = 0; i != _res_node_num; ++i) {
   765         _pi[i] = 0;
   766         _excess[i] = _supply[i];
   767       }
   768 
   769       // Remove infinite upper bounds and check negative arcs
   770       const Value MAX = std::numeric_limits<Value>::max();
   771       int last_out;
   772       if (_have_lower) {
   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]) {
   777               Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
   778               if (c >= MAX) return UNBOUNDED;
   779               _excess[i] -= c;
   780               _excess[_target[j]] += c;
   781             }
   782           }
   783         }
   784       } else {
   785         for (int i = 0; i != _root; ++i) {
   786           last_out = _first_out[i+1];
   787           for (int j = _first_out[i]; j != last_out; ++j) {
   788             if (_forward[j] && _scost[j] < 0) {
   789               Value c = _upper[j];
   790               if (c >= MAX) return UNBOUNDED;
   791               _excess[i] -= c;
   792               _excess[_target[j]] += c;
   793             }
   794           }
   795         }
   796       }
   797       Value ex, max_cap = 0;
   798       for (int i = 0; i != _res_node_num; ++i) {
   799         ex = _excess[i];
   800         _excess[i] = 0;
   801         if (ex < 0) max_cap -= ex;
   802       }
   803       for (int j = 0; j != _res_arc_num; ++j) {
   804         if (_upper[j] >= MAX) _upper[j] = max_cap;
   805       }
   806 
   807       // Initialize the large cost vector and the epsilon parameter
   808       _epsilon = 0;
   809       LargeCost lc;
   810       for (int i = 0; i != _root; ++i) {
   811         last_out = _first_out[i+1];
   812         for (int j = _first_out[i]; j != last_out; ++j) {
   813           lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
   814           _cost[j] = lc;
   815           if (lc > _epsilon) _epsilon = lc;
   816         }
   817       }
   818       _epsilon /= _alpha;
   819 
   820       // Initialize maps for Circulation and remove non-zero lower bounds
   821       ConstMap<Arc, Value> low(0);
   822       typedef typename Digraph::template ArcMap<Value> ValueArcMap;
   823       typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
   824       ValueArcMap cap(_graph), flow(_graph);
   825       ValueNodeMap sup(_graph);
   826       for (NodeIt n(_graph); n != INVALID; ++n) {
   827         sup[n] = _supply[_node_id[n]];
   828       }
   829       if (_have_lower) {
   830         for (ArcIt a(_graph); a != INVALID; ++a) {
   831           int j = _arc_idf[a];
   832           Value c = _lower[j];
   833           cap[a] = _upper[j] - c;
   834           sup[_graph.source(a)] -= c;
   835           sup[_graph.target(a)] += c;
   836         }
   837       } else {
   838         for (ArcIt a(_graph); a != INVALID; ++a) {
   839           cap[a] = _upper[_arc_idf[a]];
   840         }
   841       }
   842 
   843       _sup_node_num = 0;
   844       for (NodeIt n(_graph); n != INVALID; ++n) {
   845         if (sup[n] > 0) ++_sup_node_num;
   846       }
   847 
   848       // Find a feasible flow using Circulation
   849       Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
   850         circ(_graph, low, cap, sup);
   851       if (!circ.flowMap(flow).run()) return INFEASIBLE;
   852 
   853       // Set residual capacities and handle GEQ supply type
   854       if (_sum_supply < 0) {
   855         for (ArcIt a(_graph); a != INVALID; ++a) {
   856           Value fa = flow[a];
   857           _res_cap[_arc_idf[a]] = cap[a] - fa;
   858           _res_cap[_arc_idb[a]] = fa;
   859           sup[_graph.source(a)] -= fa;
   860           sup[_graph.target(a)] += fa;
   861         }
   862         for (NodeIt n(_graph); n != INVALID; ++n) {
   863           _excess[_node_id[n]] = sup[n];
   864         }
   865         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   866           int u = _target[a];
   867           int ra = _reverse[a];
   868           _res_cap[a] = -_sum_supply + 1;
   869           _res_cap[ra] = -_excess[u];
   870           _cost[a] = 0;
   871           _cost[ra] = 0;
   872           _excess[u] = 0;
   873         }
   874       } else {
   875         for (ArcIt a(_graph); a != INVALID; ++a) {
   876           Value fa = flow[a];
   877           _res_cap[_arc_idf[a]] = cap[a] - fa;
   878           _res_cap[_arc_idb[a]] = fa;
   879         }
   880         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   881           int ra = _reverse[a];
   882           _res_cap[a] = 0;
   883           _res_cap[ra] = 0;
   884           _cost[a] = 0;
   885           _cost[ra] = 0;
   886         }
   887       }
   888 
   889       return OPTIMAL;
   890     }
   891 
   892     // Execute the algorithm and transform the results
   893     void start(Method method) {
   894       // Maximum path length for partial augment
   895       const int MAX_PATH_LENGTH = 4;
   896 
   897       // Initialize data structures for buckets
   898       _max_rank = _alpha * _res_node_num;
   899       _buckets.resize(_max_rank);
   900       _bucket_next.resize(_res_node_num + 1);
   901       _bucket_prev.resize(_res_node_num + 1);
   902       _rank.resize(_res_node_num + 1);
   903 
   904       // Execute the algorithm
   905       switch (method) {
   906         case PUSH:
   907           startPush();
   908           break;
   909         case AUGMENT:
   910           startAugment();
   911           break;
   912         case PARTIAL_AUGMENT:
   913           startAugment(MAX_PATH_LENGTH);
   914           break;
   915       }
   916 
   917       // Compute node potentials for the original costs
   918       _arc_vec.clear();
   919       _cost_vec.clear();
   920       for (int j = 0; j != _res_arc_num; ++j) {
   921         if (_res_cap[j] > 0) {
   922           _arc_vec.push_back(IntPair(_source[j], _target[j]));
   923           _cost_vec.push_back(_scost[j]);
   924         }
   925       }
   926       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   927 
   928       typename BellmanFord<StaticDigraph, LargeCostArcMap>
   929         ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
   930       bf.distMap(_pi_map);
   931       bf.init(0);
   932       bf.start();
   933 
   934       // Handle non-zero lower bounds
   935       if (_have_lower) {
   936         int limit = _first_out[_root];
   937         for (int j = 0; j != limit; ++j) {
   938           if (!_forward[j]) _res_cap[j] += _lower[j];
   939         }
   940       }
   941     }
   942 
   943     // Initialize a cost scaling phase
   944     void initPhase() {
   945       // Saturate arcs not satisfying the optimality condition
   946       for (int u = 0; u != _res_node_num; ++u) {
   947         int last_out = _first_out[u+1];
   948         LargeCost pi_u = _pi[u];
   949         for (int a = _first_out[u]; a != last_out; ++a) {
   950           int v = _target[a];
   951           if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
   952             Value delta = _res_cap[a];
   953             _excess[u] -= delta;
   954             _excess[v] += delta;
   955             _res_cap[a] = 0;
   956             _res_cap[_reverse[a]] += delta;
   957           }
   958         }
   959       }
   960 
   961       // Find active nodes (i.e. nodes with positive excess)
   962       for (int u = 0; u != _res_node_num; ++u) {
   963         if (_excess[u] > 0) _active_nodes.push_back(u);
   964       }
   965 
   966       // Initialize the next arcs
   967       for (int u = 0; u != _res_node_num; ++u) {
   968         _next_out[u] = _first_out[u];
   969       }
   970     }
   971 
   972     // Early termination heuristic
   973     bool earlyTermination() {
   974       const double EARLY_TERM_FACTOR = 3.0;
   975 
   976       // Build a static residual graph
   977       _arc_vec.clear();
   978       _cost_vec.clear();
   979       for (int j = 0; j != _res_arc_num; ++j) {
   980         if (_res_cap[j] > 0) {
   981           _arc_vec.push_back(IntPair(_source[j], _target[j]));
   982           _cost_vec.push_back(_cost[j] + 1);
   983         }
   984       }
   985       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   986 
   987       // Run Bellman-Ford algorithm to check if the current flow is optimal
   988       BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   989       bf.init(0);
   990       bool done = false;
   991       int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
   992       for (int i = 0; i < K && !done; ++i) {
   993         done = bf.processNextWeakRound();
   994       }
   995       return done;
   996     }
   997 
   998     // Global potential update heuristic
   999     void globalUpdate() {
  1000       int bucket_end = _root + 1;
  1001 
  1002       // Initialize buckets
  1003       for (int r = 0; r != _max_rank; ++r) {
  1004         _buckets[r] = bucket_end;
  1005       }
  1006       Value total_excess = 0;
  1007       for (int i = 0; i != _res_node_num; ++i) {
  1008         if (_excess[i] < 0) {
  1009           _rank[i] = 0;
  1010           _bucket_next[i] = _buckets[0];
  1011           _bucket_prev[_buckets[0]] = i;
  1012           _buckets[0] = i;
  1013         } else {
  1014           total_excess += _excess[i];
  1015           _rank[i] = _max_rank;
  1016         }
  1017       }
  1018       if (total_excess == 0) return;
  1019 
  1020       // Search the buckets
  1021       int r = 0;
  1022       for ( ; r != _max_rank; ++r) {
  1023         while (_buckets[r] != bucket_end) {
  1024           // Remove the first node from the current bucket
  1025           int u = _buckets[r];
  1026           _buckets[r] = _bucket_next[u];
  1027 
  1028           // Search the incomming arcs of u
  1029           LargeCost pi_u = _pi[u];
  1030           int last_out = _first_out[u+1];
  1031           for (int a = _first_out[u]; a != last_out; ++a) {
  1032             int ra = _reverse[a];
  1033             if (_res_cap[ra] > 0) {
  1034               int v = _source[ra];
  1035               int old_rank_v = _rank[v];
  1036               if (r < old_rank_v) {
  1037                 // Compute the new rank of v
  1038                 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
  1039                 int new_rank_v = old_rank_v;
  1040                 if (nrc < LargeCost(_max_rank))
  1041                   new_rank_v = r + 1 + int(nrc);
  1042 
  1043                 // Change the rank of v
  1044                 if (new_rank_v < old_rank_v) {
  1045                   _rank[v] = new_rank_v;
  1046                   _next_out[v] = _first_out[v];
  1047 
  1048                   // Remove v from its old bucket
  1049                   if (old_rank_v < _max_rank) {
  1050                     if (_buckets[old_rank_v] == v) {
  1051                       _buckets[old_rank_v] = _bucket_next[v];
  1052                     } else {
  1053                       _bucket_next[_bucket_prev[v]] = _bucket_next[v];
  1054                       _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
  1055                     }
  1056                   }
  1057 
  1058                   // Insert v to its new bucket
  1059                   _bucket_next[v] = _buckets[new_rank_v];
  1060                   _bucket_prev[_buckets[new_rank_v]] = v;
  1061                   _buckets[new_rank_v] = v;
  1062                 }
  1063               }
  1064             }
  1065           }
  1066 
  1067           // Finish search if there are no more active nodes
  1068           if (_excess[u] > 0) {
  1069             total_excess -= _excess[u];
  1070             if (total_excess <= 0) break;
  1071           }
  1072         }
  1073         if (total_excess <= 0) break;
  1074       }
  1075 
  1076       // Relabel nodes
  1077       for (int u = 0; u != _res_node_num; ++u) {
  1078         int k = std::min(_rank[u], r);
  1079         if (k > 0) {
  1080           _pi[u] -= _epsilon * k;
  1081           _next_out[u] = _first_out[u];
  1082         }
  1083       }
  1084     }
  1085 
  1086     /// Execute the algorithm performing augment and relabel operations
  1087     void startAugment(int max_length = std::numeric_limits<int>::max()) {
  1088       // Paramters for heuristics
  1089       const int EARLY_TERM_EPSILON_LIMIT = 1000;
  1090       const double GLOBAL_UPDATE_FACTOR = 3.0;
  1091 
  1092       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
  1093         (_res_node_num + _sup_node_num * _sup_node_num));
  1094       int next_update_limit = global_update_freq;
  1095 
  1096       int relabel_cnt = 0;
  1097 
  1098       // Perform cost scaling phases
  1099       std::vector<int> path;
  1100       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1101                                         1 : _epsilon / _alpha )
  1102       {
  1103         // Early termination heuristic
  1104         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
  1105           if (earlyTermination()) break;
  1106         }
  1107 
  1108         // Initialize current phase
  1109         initPhase();
  1110 
  1111         // Perform partial augment and relabel operations
  1112         while (true) {
  1113           // Select an active node (FIFO selection)
  1114           while (_active_nodes.size() > 0 &&
  1115                  _excess[_active_nodes.front()] <= 0) {
  1116             _active_nodes.pop_front();
  1117           }
  1118           if (_active_nodes.size() == 0) break;
  1119           int start = _active_nodes.front();
  1120 
  1121           // Find an augmenting path from the start node
  1122           path.clear();
  1123           int tip = start;
  1124           while (_excess[tip] >= 0 && int(path.size()) < max_length) {
  1125             int u;
  1126             LargeCost min_red_cost, rc, pi_tip = _pi[tip];
  1127             int last_out = _first_out[tip+1];
  1128             for (int a = _next_out[tip]; a != last_out; ++a) {
  1129               u = _target[a];
  1130               if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
  1131                 path.push_back(a);
  1132                 _next_out[tip] = a;
  1133                 tip = u;
  1134                 goto next_step;
  1135               }
  1136             }
  1137 
  1138             // Relabel tip node
  1139             min_red_cost = std::numeric_limits<LargeCost>::max();
  1140             if (tip != start) {
  1141               int ra = _reverse[path.back()];
  1142               min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
  1143             }
  1144             for (int a = _first_out[tip]; a != last_out; ++a) {
  1145               rc = _cost[a] + pi_tip - _pi[_target[a]];
  1146               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1147                 min_red_cost = rc;
  1148               }
  1149             }
  1150             _pi[tip] -= min_red_cost + _epsilon;
  1151             _next_out[tip] = _first_out[tip];
  1152             ++relabel_cnt;
  1153 
  1154             // Step back
  1155             if (tip != start) {
  1156               tip = _source[path.back()];
  1157               path.pop_back();
  1158             }
  1159 
  1160           next_step: ;
  1161           }
  1162 
  1163           // Augment along the found path (as much flow as possible)
  1164           Value delta;
  1165           int pa, u, v = start;
  1166           for (int i = 0; i != int(path.size()); ++i) {
  1167             pa = path[i];
  1168             u = v;
  1169             v = _target[pa];
  1170             delta = std::min(_res_cap[pa], _excess[u]);
  1171             _res_cap[pa] -= delta;
  1172             _res_cap[_reverse[pa]] += delta;
  1173             _excess[u] -= delta;
  1174             _excess[v] += delta;
  1175             if (_excess[v] > 0 && _excess[v] <= delta)
  1176               _active_nodes.push_back(v);
  1177           }
  1178 
  1179           // Global update heuristic
  1180           if (relabel_cnt >= next_update_limit) {
  1181             globalUpdate();
  1182             next_update_limit += global_update_freq;
  1183           }
  1184         }
  1185       }
  1186     }
  1187 
  1188     /// Execute the algorithm performing push and relabel operations
  1189     void startPush() {
  1190       // Paramters for heuristics
  1191       const int EARLY_TERM_EPSILON_LIMIT = 1000;
  1192       const double GLOBAL_UPDATE_FACTOR = 2.0;
  1193 
  1194       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
  1195         (_res_node_num + _sup_node_num * _sup_node_num));
  1196       int next_update_limit = global_update_freq;
  1197 
  1198       int relabel_cnt = 0;
  1199 
  1200       // Perform cost scaling phases
  1201       BoolVector hyper(_res_node_num, false);
  1202       LargeCostVector hyper_cost(_res_node_num);
  1203       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1204                                         1 : _epsilon / _alpha )
  1205       {
  1206         // Early termination heuristic
  1207         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
  1208           if (earlyTermination()) break;
  1209         }
  1210 
  1211         // Initialize current phase
  1212         initPhase();
  1213 
  1214         // Perform push and relabel operations
  1215         while (_active_nodes.size() > 0) {
  1216           LargeCost min_red_cost, rc, pi_n;
  1217           Value delta;
  1218           int n, t, a, last_out = _res_arc_num;
  1219 
  1220         next_node:
  1221           // Select an active node (FIFO selection)
  1222           n = _active_nodes.front();
  1223           last_out = _first_out[n+1];
  1224           pi_n = _pi[n];
  1225 
  1226           // Perform push operations if there are admissible arcs
  1227           if (_excess[n] > 0) {
  1228             for (a = _next_out[n]; a != last_out; ++a) {
  1229               if (_res_cap[a] > 0 &&
  1230                   _cost[a] + pi_n - _pi[_target[a]] < 0) {
  1231                 delta = std::min(_res_cap[a], _excess[n]);
  1232                 t = _target[a];
  1233 
  1234                 // Push-look-ahead heuristic
  1235                 Value ahead = -_excess[t];
  1236                 int last_out_t = _first_out[t+1];
  1237                 LargeCost pi_t = _pi[t];
  1238                 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
  1239                   if (_res_cap[ta] > 0 &&
  1240                       _cost[ta] + pi_t - _pi[_target[ta]] < 0)
  1241                     ahead += _res_cap[ta];
  1242                   if (ahead >= delta) break;
  1243                 }
  1244                 if (ahead < 0) ahead = 0;
  1245 
  1246                 // Push flow along the arc
  1247                 if (ahead < delta && !hyper[t]) {
  1248                   _res_cap[a] -= ahead;
  1249                   _res_cap[_reverse[a]] += ahead;
  1250                   _excess[n] -= ahead;
  1251                   _excess[t] += ahead;
  1252                   _active_nodes.push_front(t);
  1253                   hyper[t] = true;
  1254                   hyper_cost[t] = _cost[a] + pi_n - pi_t;
  1255                   _next_out[n] = a;
  1256                   goto next_node;
  1257                 } else {
  1258                   _res_cap[a] -= delta;
  1259                   _res_cap[_reverse[a]] += delta;
  1260                   _excess[n] -= delta;
  1261                   _excess[t] += delta;
  1262                   if (_excess[t] > 0 && _excess[t] <= delta)
  1263                     _active_nodes.push_back(t);
  1264                 }
  1265 
  1266                 if (_excess[n] == 0) {
  1267                   _next_out[n] = a;
  1268                   goto remove_nodes;
  1269                 }
  1270               }
  1271             }
  1272             _next_out[n] = a;
  1273           }
  1274 
  1275           // Relabel the node if it is still active (or hyper)
  1276           if (_excess[n] > 0 || hyper[n]) {
  1277              min_red_cost = hyper[n] ? -hyper_cost[n] :
  1278                std::numeric_limits<LargeCost>::max();
  1279             for (int a = _first_out[n]; a != last_out; ++a) {
  1280               rc = _cost[a] + pi_n - _pi[_target[a]];
  1281               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1282                 min_red_cost = rc;
  1283               }
  1284             }
  1285             _pi[n] -= min_red_cost + _epsilon;
  1286             _next_out[n] = _first_out[n];
  1287             hyper[n] = false;
  1288             ++relabel_cnt;
  1289           }
  1290 
  1291           // Remove nodes that are not active nor hyper
  1292         remove_nodes:
  1293           while ( _active_nodes.size() > 0 &&
  1294                   _excess[_active_nodes.front()] <= 0 &&
  1295                   !hyper[_active_nodes.front()] ) {
  1296             _active_nodes.pop_front();
  1297           }
  1298 
  1299           // Global update heuristic
  1300           if (relabel_cnt >= next_update_limit) {
  1301             globalUpdate();
  1302             for (int u = 0; u != _res_node_num; ++u)
  1303               hyper[u] = false;
  1304             next_update_limit += global_update_freq;
  1305           }
  1306         }
  1307       }
  1308     }
  1309 
  1310   }; //class CostScaling
  1311 
  1312   ///@}
  1313 
  1314 } //namespace lemon
  1315 
  1316 #endif //LEMON_COST_SCALING_H