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