author Peter Kovacs Thu, 12 Nov 2009 23:30:45 +0100 changeset 809 22bb98ca0101 parent 808 9c428bb2b105 child 810 3b53491bf643
Entirely rework CostScaling (#180)

- Use the new interface similarly to NetworkSimplex.
- Rework the implementation using an efficient internal structure
for handling the residual network. This improvement made the
code much faster.
- Handle GEQ supply type (LEQ is not supported).
- Handle infinite upper bounds.
- Handle negative costs (for arcs of finite upper bound).
- Traits class + named parameter for the LargeCost type used in
internal computations.
- Extend the documentation.
 lemon/cost_scaling.h file | annotate | diff | comparison | revisions
     1.1 --- a/lemon/cost_scaling.h	Thu Nov 12 23:29:42 2009 +0100
1.2 +++ b/lemon/cost_scaling.h	Thu Nov 12 23:30:45 2009 +0100
1.3 @@ -30,548 +30,912 @@
1.4  #include <lemon/core.h>
1.5  #include <lemon/maps.h>
1.6  #include <lemon/math.h>
1.8 +#include <lemon/static_graph.h>
1.9  #include <lemon/circulation.h>
1.10  #include <lemon/bellman_ford.h>
1.11
1.12  namespace lemon {
1.13
1.14 +  /// \brief Default traits class of CostScaling algorithm.
1.15 +  ///
1.16 +  /// Default traits class of CostScaling algorithm.
1.17 +  /// \tparam GR Digraph type.
1.18 +  /// \tparam V The value type used for flow amounts, capacity bounds
1.19 +  /// and supply values. By default it is \c int.
1.20 +  /// \tparam C The value type used for costs and potentials.
1.21 +  /// By default it is the same as \c V.
1.22 +#ifdef DOXYGEN
1.23 +  template <typename GR, typename V = int, typename C = V>
1.24 +#else
1.25 +  template < typename GR, typename V = int, typename C = V,
1.26 +             bool integer = std::numeric_limits<C>::is_integer >
1.27 +#endif
1.28 +  struct CostScalingDefaultTraits
1.29 +  {
1.30 +    /// The type of the digraph
1.31 +    typedef GR Digraph;
1.32 +    /// The type of the flow amounts, capacity bounds and supply values
1.33 +    typedef V Value;
1.34 +    /// The type of the arc costs
1.35 +    typedef C Cost;
1.36 +
1.37 +    /// \brief The large cost type used for internal computations
1.38 +    ///
1.39 +    /// The large cost type used for internal computations.
1.40 +    /// It is \c long \c long if the \c Cost type is integer,
1.41 +    /// otherwise it is \c double.
1.42 +    /// \c Cost must be convertible to \c LargeCost.
1.43 +    typedef double LargeCost;
1.44 +  };
1.45 +
1.46 +  // Default traits class for integer cost types
1.47 +  template <typename GR, typename V, typename C>
1.48 +  struct CostScalingDefaultTraits<GR, V, C, true>
1.49 +  {
1.50 +    typedef GR Digraph;
1.51 +    typedef V Value;
1.52 +    typedef C Cost;
1.53 +#ifdef LEMON_HAVE_LONG_LONG
1.54 +    typedef long long LargeCost;
1.55 +#else
1.56 +    typedef long LargeCost;
1.57 +#endif
1.58 +  };
1.59 +
1.60 +
1.62    /// @{
1.63
1.64 -  /// \brief Implementation of the cost scaling algorithm for finding a
1.65 -  /// minimum cost flow.
1.66 +  /// \brief Implementation of the Cost Scaling algorithm for
1.67 +  /// finding a \ref min_cost_flow "minimum cost flow".
1.68    ///
1.69 -  /// \ref CostScaling implements the cost scaling algorithm performing
1.70 -  /// augment/push and relabel operations for finding a minimum cost
1.71 -  /// flow.
1.72 +  /// \ref CostScaling implements a cost scaling algorithm that performs
1.73 +  /// push/augment and relabel operations for finding a minimum cost
1.74 +  /// flow. It is an efficient primal-dual solution method, which
1.75 +  /// can be viewed as the generalization of the \ref Preflow
1.76 +  /// "preflow push-relabel" algorithm for the maximum flow problem.
1.77    ///
1.78 -  /// \tparam Digraph The digraph type the algorithm runs on.
1.79 -  /// \tparam LowerMap The type of the lower bound map.
1.80 -  /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.81 -  /// \tparam CostMap The type of the cost (length) map.
1.82 -  /// \tparam SupplyMap The type of the supply map.
1.83 +  /// Most of the parameters of the problem (except for the digraph)
1.84 +  /// can be given using separate functions, and the algorithm can be
1.85 +  /// executed using the \ref run() function. If some parameters are not
1.86 +  /// specified, then default values will be used.
1.87    ///
1.88 -  /// \warning
1.89 -  /// - Arc capacities and costs should be \e non-negative \e integers.
1.90 -  /// - Supply values should be \e signed \e integers.
1.91 -  /// - The value types of the maps should be convertible to each other.
1.92 -  /// - \c CostMap::Value must be signed type.
1.93 +  /// \tparam GR The digraph type the algorithm runs on.
1.94 +  /// \tparam V The value type used for flow amounts, capacity bounds
1.95 +  /// and supply values in the algorithm. By default it is \c int.
1.96 +  /// \tparam C The value type used for costs and potentials in the
1.97 +  /// algorithm. By default it is the same as \c V.
1.98    ///
1.99 -  /// \note Arc costs are multiplied with the number of nodes during
1.100 -  /// the algorithm so overflow problems may arise more easily than with
1.101 -  /// other minimum cost flow algorithms.
1.102 -  /// If it is available, <tt>long long int</tt> type is used instead of
1.103 -  /// <tt>long int</tt> in the inside computations.
1.104 -  ///
1.105 -  /// \author Peter Kovacs
1.106 -  template < typename Digraph,
1.107 -             typename LowerMap = typename Digraph::template ArcMap<int>,
1.108 -             typename CapacityMap = typename Digraph::template ArcMap<int>,
1.109 -             typename CostMap = typename Digraph::template ArcMap<int>,
1.110 -             typename SupplyMap = typename Digraph::template NodeMap<int> >
1.111 +  /// \warning Both value types must be signed and all input data must
1.112 +  /// be integer.
1.113 +  /// \warning This algorithm does not support negative costs for such
1.114 +  /// arcs that have infinite upper bound.
1.115 +#ifdef DOXYGEN
1.116 +  template <typename GR, typename V, typename C, typename TR>
1.117 +#else
1.118 +  template < typename GR, typename V = int, typename C = V,
1.119 +             typename TR = CostScalingDefaultTraits<GR, V, C> >
1.120 +#endif
1.121    class CostScaling
1.122    {
1.123 -    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.124 +  public:
1.125
1.126 -    typedef typename CapacityMap::Value Capacity;
1.127 -    typedef typename CostMap::Value Cost;
1.128 -    typedef typename SupplyMap::Value Supply;
1.129 -    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.130 -    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.131 +    /// The type of the digraph
1.132 +    typedef typename TR::Digraph Digraph;
1.133 +    /// The type of the flow amounts, capacity bounds and supply values
1.134 +    typedef typename TR::Value Value;
1.135 +    /// The type of the arc costs
1.136 +    typedef typename TR::Cost Cost;
1.137
1.138 -    typedef ResidualDigraph< const Digraph,
1.139 -                             CapacityArcMap, CapacityArcMap > ResDigraph;
1.140 -    typedef typename ResDigraph::Arc ResArc;
1.141 +    /// \brief The large cost type
1.142 +    ///
1.143 +    /// The large cost type used for internal computations.
1.144 +    /// Using the \ref CostScalingDefaultTraits "default traits class",
1.145 +    /// it is \c long \c long if the \c Cost type is integer,
1.146 +    /// otherwise it is \c double.
1.147 +    typedef typename TR::LargeCost LargeCost;
1.148
1.149 -#if defined __GNUC__ && !defined __STRICT_ANSI__
1.150 -    typedef long long int LCost;
1.151 -#else
1.152 -    typedef long int LCost;
1.153 -#endif
1.154 -    typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
1.155 +    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
1.156 +    typedef TR Traits;
1.157
1.158    public:
1.159
1.160 -    /// The type of the flow map.
1.161 -    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.162 -    /// The type of the potential map.
1.163 -    typedef typename Digraph::template NodeMap<LCost> PotentialMap;
1.164 +    /// \brief Problem type constants for the \c run() function.
1.165 +    ///
1.166 +    /// Enum type containing the problem type constants that can be
1.167 +    /// returned by the \ref run() function of the algorithm.
1.168 +    enum ProblemType {
1.169 +      /// The problem has no feasible solution (flow).
1.170 +      INFEASIBLE,
1.171 +      /// The problem has optimal solution (i.e. it is feasible and
1.172 +      /// bounded), and the algorithm has found optimal flow and node
1.173 +      /// potentials (primal and dual solutions).
1.174 +      OPTIMAL,
1.175 +      /// The digraph contains an arc of negative cost and infinite
1.176 +      /// upper bound. It means that the objective function is unbounded
1.177 +      /// on that arc, however note that it could actually be bounded
1.178 +      /// over the feasible flows, but this algroithm cannot handle
1.179 +      /// these cases.
1.180 +      UNBOUNDED
1.181 +    };
1.182
1.183    private:
1.184
1.185 -    /// \brief Map adaptor class for handling residual arc costs.
1.186 -    ///
1.187 -    /// Map adaptor class for handling residual arc costs.
1.188 -    template <typename Map>
1.189 -    class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
1.190 -    {
1.191 -    private:
1.192 +    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.193
1.194 -      const Map &_cost_map;
1.195 +    typedef std::vector<int> IntVector;
1.196 +    typedef std::vector<char> BoolVector;
1.197 +    typedef std::vector<Value> ValueVector;
1.198 +    typedef std::vector<Cost> CostVector;
1.199 +    typedef std::vector<LargeCost> LargeCostVector;
1.200
1.201 +  private:
1.202 +
1.203 +    template <typename KT, typename VT>
1.204 +    class VectorMap {
1.205      public:
1.206 -
1.207 -      ///\e
1.208 -      ResidualCostMap(const Map &cost_map) :
1.209 -        _cost_map(cost_map) {}
1.210 -
1.211 -      ///\e
1.212 -      inline typename Map::Value operator[](const ResArc &e) const {
1.213 -        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.214 +      typedef KT Key;
1.215 +      typedef VT Value;
1.216 +
1.217 +      VectorMap(std::vector<Value>& v) : _v(v) {}
1.218 +
1.219 +      const Value& operator[](const Key& key) const {
1.220 +        return _v[StaticDigraph::id(key)];
1.221        }
1.222
1.223 -    }; //class ResidualCostMap
1.224 -
1.225 -    /// \brief Map adaptor class for handling reduced arc costs.
1.226 -    ///
1.227 -    /// Map adaptor class for handling reduced arc costs.
1.228 -    class ReducedCostMap : public MapBase<Arc, LCost>
1.229 -    {
1.230 -    private:
1.231 -
1.232 -      const Digraph &_gr;
1.233 -      const LargeCostMap &_cost_map;
1.234 -      const PotentialMap &_pot_map;
1.235 -
1.236 -    public:
1.237 -
1.238 -      ///\e
1.239 -      ReducedCostMap( const Digraph &gr,
1.240 -                      const LargeCostMap &cost_map,
1.241 -                      const PotentialMap &pot_map ) :
1.242 -        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
1.243 -
1.244 -      ///\e
1.245 -      inline LCost operator[](const Arc &e) const {
1.246 -        return _cost_map[e] + _pot_map[_gr.source(e)]
1.247 -                            - _pot_map[_gr.target(e)];
1.248 +      Value& operator[](const Key& key) {
1.249 +        return _v[StaticDigraph::id(key)];
1.250 +      }
1.251 +
1.252 +      void set(const Key& key, const Value& val) {
1.253 +        _v[StaticDigraph::id(key)] = val;
1.254        }
1.255
1.256 -    }; //class ReducedCostMap
1.257 +    private:
1.258 +      std::vector<Value>& _v;
1.259 +    };
1.260 +
1.261 +    typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
1.262 +    typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
1.263
1.264    private:
1.265
1.266 -    // The digraph the algorithm runs on
1.267 -    const Digraph &_graph;
1.268 -    // The original lower bound map
1.269 -    const LowerMap *_lower;
1.270 -    // The modified capacity map
1.271 -    CapacityArcMap _capacity;
1.272 -    // The original cost map
1.273 -    const CostMap &_orig_cost;
1.274 -    // The scaled cost map
1.275 -    LargeCostMap _cost;
1.276 -    // The modified supply map
1.277 -    SupplyNodeMap _supply;
1.278 -    bool _valid_supply;
1.279 +    // Data related to the underlying digraph
1.280 +    const GR &_graph;
1.281 +    int _node_num;
1.282 +    int _arc_num;
1.283 +    int _res_node_num;
1.284 +    int _res_arc_num;
1.285 +    int _root;
1.286
1.287 -    // Arc map of the current flow
1.288 -    FlowMap *_flow;
1.289 -    bool _local_flow;
1.290 -    // Node map of the current potentials
1.291 -    PotentialMap *_potential;
1.292 -    bool _local_potential;
1.293 +    // Parameters of the problem
1.294 +    bool _have_lower;
1.295 +    Value _sum_supply;
1.296
1.297 -    // The residual cost map
1.298 -    ResidualCostMap<LargeCostMap> _res_cost;
1.299 -    // The residual digraph
1.300 -    ResDigraph *_res_graph;
1.301 -    // The reduced cost map
1.302 -    ReducedCostMap *_red_cost;
1.303 -    // The excess map
1.304 -    SupplyNodeMap _excess;
1.305 -    // The epsilon parameter used for cost scaling
1.306 -    LCost _epsilon;
1.307 -    // The scaling factor
1.308 +    // Data structures for storing the digraph
1.309 +    IntNodeMap _node_id;
1.310 +    IntArcMap _arc_idf;
1.311 +    IntArcMap _arc_idb;
1.312 +    IntVector _first_out;
1.313 +    BoolVector _forward;
1.314 +    IntVector _source;
1.315 +    IntVector _target;
1.316 +    IntVector _reverse;
1.317 +
1.318 +    // Node and arc data
1.319 +    ValueVector _lower;
1.320 +    ValueVector _upper;
1.321 +    CostVector _scost;
1.322 +    ValueVector _supply;
1.323 +
1.324 +    ValueVector _res_cap;
1.325 +    LargeCostVector _cost;
1.326 +    LargeCostVector _pi;
1.327 +    ValueVector _excess;
1.328 +    IntVector _next_out;
1.329 +    std::deque<int> _active_nodes;
1.330 +
1.331 +    // Data for scaling
1.332 +    LargeCost _epsilon;
1.333      int _alpha;
1.334
1.335 +    // Data for a StaticDigraph structure
1.336 +    typedef std::pair<int, int> IntPair;
1.337 +    StaticDigraph _sgr;
1.338 +    std::vector<IntPair> _arc_vec;
1.339 +    std::vector<LargeCost> _cost_vec;
1.340 +    LargeCostArcMap _cost_map;
1.341 +    LargeCostNodeMap _pi_map;
1.342 +
1.343 +  public:
1.344 +
1.345 +    /// \brief Constant for infinite upper bounds (capacities).
1.346 +    ///
1.347 +    /// Constant for infinite upper bounds (capacities).
1.348 +    /// It is \c std::numeric_limits<Value>::infinity() if available,
1.349 +    /// \c std::numeric_limits<Value>::max() otherwise.
1.350 +    const Value INF;
1.351 +
1.352    public:
1.353
1.354 -    /// \brief General constructor (with lower bounds).
1.355 +    /// \name Named Template Parameters
1.356 +    /// @{
1.357 +
1.358 +    template <typename T>
1.359 +    struct SetLargeCostTraits : public Traits {
1.360 +      typedef T LargeCost;
1.361 +    };
1.362 +
1.363 +    /// \brief \ref named-templ-param "Named parameter" for setting
1.364 +    /// \c LargeCost type.
1.365      ///
1.366 -    /// General constructor (with lower bounds).
1.367 +    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
1.368 +    /// type, which is used for internal computations in the algorithm.
1.369 +    /// \c Cost must be convertible to \c LargeCost.
1.370 +    template <typename T>
1.371 +    struct SetLargeCost
1.372 +      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
1.373 +      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
1.374 +    };
1.375 +
1.376 +    /// @}
1.377 +
1.378 +  public:
1.379 +
1.380 +    /// \brief Constructor.
1.381      ///
1.382 -    /// \param digraph The digraph the algorithm runs on.
1.383 -    /// \param lower The lower bounds of the arcs.
1.384 -    /// \param capacity The capacities (upper bounds) of the arcs.
1.385 -    /// \param cost The cost (length) values of the arcs.
1.386 -    /// \param supply The supply values of the nodes (signed).
1.387 -    CostScaling( const Digraph &digraph,
1.388 -                 const LowerMap &lower,
1.389 -                 const CapacityMap &capacity,
1.390 -                 const CostMap &cost,
1.391 -                 const SupplyMap &supply ) :
1.392 -      _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
1.393 -      _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
1.394 -      _potential(NULL), _local_potential(false), _res_cost(_cost),
1.395 -      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.396 +    /// The constructor of the class.
1.397 +    ///
1.398 +    /// \param graph The digraph the algorithm runs on.
1.399 +    CostScaling(const GR& graph) :
1.400 +      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
1.401 +      _cost_map(_cost_vec), _pi_map(_pi),
1.402 +      INF(std::numeric_limits<Value>::has_infinity ?
1.403 +          std::numeric_limits<Value>::infinity() :
1.404 +          std::numeric_limits<Value>::max())
1.405      {
1.406 -      // Check the sum of supply values
1.407 -      Supply sum = 0;
1.408 -      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.409 -      _valid_supply = sum == 0;
1.410 +      // Check the value types
1.411 +      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.412 +        "The flow type of CostScaling must be signed");
1.413 +      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.414 +        "The cost type of CostScaling must be signed");
1.415 +
1.416 +      // Resize vectors
1.417 +      _node_num = countNodes(_graph);
1.418 +      _arc_num = countArcs(_graph);
1.419 +      _res_node_num = _node_num + 1;
1.420 +      _res_arc_num = 2 * (_arc_num + _node_num);
1.421 +      _root = _node_num;
1.422 +
1.423 +      _first_out.resize(_res_node_num + 1);
1.424 +      _forward.resize(_res_arc_num);
1.425 +      _source.resize(_res_arc_num);
1.426 +      _target.resize(_res_arc_num);
1.427 +      _reverse.resize(_res_arc_num);
1.428 +
1.429 +      _lower.resize(_res_arc_num);
1.430 +      _upper.resize(_res_arc_num);
1.431 +      _scost.resize(_res_arc_num);
1.432 +      _supply.resize(_res_node_num);
1.433
1.434 -      for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
1.435 -      for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
1.436 +      _res_cap.resize(_res_arc_num);
1.437 +      _cost.resize(_res_arc_num);
1.438 +      _pi.resize(_res_node_num);
1.439 +      _excess.resize(_res_node_num);
1.440 +      _next_out.resize(_res_node_num);
1.441
1.442 -      // Remove non-zero lower bounds
1.443 -      for (ArcIt e(_graph); e != INVALID; ++e) {
1.444 -        if (lower[e] != 0) {
1.445 -          _capacity[e] -= lower[e];
1.446 -          _supply[_graph.source(e)] -= lower[e];
1.447 -          _supply[_graph.target(e)] += lower[e];
1.448 +      _arc_vec.reserve(_res_arc_num);
1.449 +      _cost_vec.reserve(_res_arc_num);
1.450 +
1.451 +      // Copy the graph
1.452 +      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
1.453 +      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.454 +        _node_id[n] = i;
1.455 +      }
1.456 +      i = 0;
1.457 +      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.458 +        _first_out[i] = j;
1.459 +        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.460 +          _arc_idf[a] = j;
1.461 +          _forward[j] = true;
1.462 +          _source[j] = i;
1.463 +          _target[j] = _node_id[_graph.runningNode(a)];
1.464          }
1.465 +        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.466 +          _arc_idb[a] = j;
1.467 +          _forward[j] = false;
1.468 +          _source[j] = i;
1.469 +          _target[j] = _node_id[_graph.runningNode(a)];
1.470 +        }
1.471 +        _forward[j] = false;
1.472 +        _source[j] = i;
1.473 +        _target[j] = _root;
1.474 +        _reverse[j] = k;
1.475 +        _forward[k] = true;
1.476 +        _source[k] = _root;
1.477 +        _target[k] = i;
1.478 +        _reverse[k] = j;
1.479 +        ++j; ++k;
1.480        }
1.481 -    }
1.482 -/*
1.483 -    /// \brief General constructor (without lower bounds).
1.484 -    ///
1.485 -    /// General constructor (without lower bounds).
1.486 -    ///
1.487 -    /// \param digraph The digraph the algorithm runs on.
1.488 -    /// \param capacity The capacities (upper bounds) of the arcs.
1.489 -    /// \param cost The cost (length) values of the arcs.
1.490 -    /// \param supply The supply values of the nodes (signed).
1.491 -    CostScaling( const Digraph &digraph,
1.492 -                 const CapacityMap &capacity,
1.493 -                 const CostMap &cost,
1.494 -                 const SupplyMap &supply ) :
1.495 -      _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.496 -      _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
1.497 -      _potential(NULL), _local_potential(false), _res_cost(_cost),
1.498 -      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.499 -    {
1.500 -      // Check the sum of supply values
1.501 -      Supply sum = 0;
1.502 -      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.503 -      _valid_supply = sum == 0;
1.504 +      _first_out[i] = j;
1.505 +      _first_out[_res_node_num] = k;
1.506 +      for (ArcIt a(_graph); a != INVALID; ++a) {
1.507 +        int fi = _arc_idf[a];
1.508 +        int bi = _arc_idb[a];
1.509 +        _reverse[fi] = bi;
1.510 +        _reverse[bi] = fi;
1.511 +      }
1.512 +
1.513 +      // Reset parameters
1.514 +      reset();
1.515      }
1.516
1.517 -    /// \brief Simple constructor (with lower bounds).
1.518 +    /// \name Parameters
1.519 +    /// The parameters of the algorithm can be specified using these
1.520 +    /// functions.
1.521 +
1.522 +    /// @{
1.523 +
1.524 +    /// \brief Set the lower bounds on the arcs.
1.525      ///
1.526 -    /// Simple constructor (with lower bounds).
1.527 +    /// This function sets the lower bounds on the arcs.
1.528 +    /// If it is not used before calling \ref run(), the lower bounds
1.529 +    /// will be set to zero on all arcs.
1.530      ///
1.531 -    /// \param digraph The digraph the algorithm runs on.
1.532 -    /// \param lower The lower bounds of the arcs.
1.533 -    /// \param capacity The capacities (upper bounds) of the arcs.
1.534 -    /// \param cost The cost (length) values of the arcs.
1.535 -    /// \param s The source node.
1.536 -    /// \param t The target node.
1.537 -    /// \param flow_value The required amount of flow from node \c s
1.538 -    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.539 -    CostScaling( const Digraph &digraph,
1.540 -                 const LowerMap &lower,
1.541 -                 const CapacityMap &capacity,
1.542 -                 const CostMap &cost,
1.543 -                 Node s, Node t,
1.544 -                 Supply flow_value ) :
1.545 -      _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
1.546 -      _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.547 -      _potential(NULL), _local_potential(false), _res_cost(_cost),
1.548 -      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.549 -    {
1.550 -      // Remove non-zero lower bounds
1.551 -      _supply[s] =  flow_value;
1.552 -      _supply[t] = -flow_value;
1.553 -      for (ArcIt e(_graph); e != INVALID; ++e) {
1.554 -        if (lower[e] != 0) {
1.555 -          _capacity[e] -= lower[e];
1.556 -          _supply[_graph.source(e)] -= lower[e];
1.557 -          _supply[_graph.target(e)] += lower[e];
1.558 -        }
1.559 +    /// \param map An arc map storing the lower bounds.
1.560 +    /// Its \c Value type must be convertible to the \c Value type
1.561 +    /// of the algorithm.
1.562 +    ///
1.563 +    /// \return <tt>(*this)</tt>
1.564 +    template <typename LowerMap>
1.565 +    CostScaling& lowerMap(const LowerMap& map) {
1.566 +      _have_lower = true;
1.567 +      for (ArcIt a(_graph); a != INVALID; ++a) {
1.568 +        _lower[_arc_idf[a]] = map[a];
1.569 +        _lower[_arc_idb[a]] = map[a];
1.570        }
1.571 -      _valid_supply = true;
1.572 -    }
1.573 -
1.574 -    /// \brief Simple constructor (without lower bounds).
1.575 -    ///
1.576 -    /// Simple constructor (without lower bounds).
1.577 -    ///
1.578 -    /// \param digraph The digraph the algorithm runs on.
1.579 -    /// \param capacity The capacities (upper bounds) of the arcs.
1.580 -    /// \param cost The cost (length) values of the arcs.
1.581 -    /// \param s The source node.
1.582 -    /// \param t The target node.
1.583 -    /// \param flow_value The required amount of flow from node \c s
1.584 -    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.585 -    CostScaling( const Digraph &digraph,
1.586 -                 const CapacityMap &capacity,
1.587 -                 const CostMap &cost,
1.588 -                 Node s, Node t,
1.589 -                 Supply flow_value ) :
1.590 -      _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.591 -      _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.592 -      _potential(NULL), _local_potential(false), _res_cost(_cost),
1.593 -      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.594 -    {
1.595 -      _supply[s] =  flow_value;
1.596 -      _supply[t] = -flow_value;
1.597 -      _valid_supply = true;
1.598 -    }
1.599 -*/
1.600 -    /// Destructor.
1.601 -    ~CostScaling() {
1.602 -      if (_local_flow) delete _flow;
1.603 -      if (_local_potential) delete _potential;
1.604 -      delete _res_graph;
1.605 -      delete _red_cost;
1.606 -    }
1.607 -
1.608 -    /// \brief Set the flow map.
1.609 -    ///
1.610 -    /// Set the flow map.
1.611 -    ///
1.612 -    /// \return \c (*this)
1.613 -    CostScaling& flowMap(FlowMap &map) {
1.614 -      if (_local_flow) {
1.615 -        delete _flow;
1.616 -        _local_flow = false;
1.617 -      }
1.618 -      _flow = &map;
1.619        return *this;
1.620      }
1.621
1.622 -    /// \brief Set the potential map.
1.623 +    /// \brief Set the upper bounds (capacities) on the arcs.
1.624      ///
1.625 -    /// Set the potential map.
1.626 +    /// This function sets the upper bounds (capacities) on the arcs.
1.627 +    /// If it is not used before calling \ref run(), the upper bounds
1.628 +    /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.629 +    /// unbounded from above on each arc).
1.630      ///
1.631 -    /// \return \c (*this)
1.632 -    CostScaling& potentialMap(PotentialMap &map) {
1.633 -      if (_local_potential) {
1.634 -        delete _potential;
1.635 -        _local_potential = false;
1.636 +    /// \param map An arc map storing the upper bounds.
1.637 +    /// Its \c Value type must be convertible to the \c Value type
1.638 +    /// of the algorithm.
1.639 +    ///
1.640 +    /// \return <tt>(*this)</tt>
1.641 +    template<typename UpperMap>
1.642 +    CostScaling& upperMap(const UpperMap& map) {
1.643 +      for (ArcIt a(_graph); a != INVALID; ++a) {
1.644 +        _upper[_arc_idf[a]] = map[a];
1.645        }
1.646 -      _potential = &map;
1.647        return *this;
1.648      }
1.649
1.650 +    /// \brief Set the costs of the arcs.
1.651 +    ///
1.652 +    /// This function sets the costs of the arcs.
1.653 +    /// If it is not used before calling \ref run(), the costs
1.654 +    /// will be set to \c 1 on all arcs.
1.655 +    ///
1.656 +    /// \param map An arc map storing the costs.
1.657 +    /// Its \c Value type must be convertible to the \c Cost type
1.658 +    /// of the algorithm.
1.659 +    ///
1.660 +    /// \return <tt>(*this)</tt>
1.661 +    template<typename CostMap>
1.662 +    CostScaling& costMap(const CostMap& map) {
1.663 +      for (ArcIt a(_graph); a != INVALID; ++a) {
1.664 +        _scost[_arc_idf[a]] =  map[a];
1.665 +        _scost[_arc_idb[a]] = -map[a];
1.666 +      }
1.667 +      return *this;
1.668 +    }
1.669 +
1.670 +    /// \brief Set the supply values of the nodes.
1.671 +    ///
1.672 +    /// This function sets the supply values of the nodes.
1.673 +    /// If neither this function nor \ref stSupply() is used before
1.674 +    /// calling \ref run(), the supply of each node will be set to zero.
1.675 +    ///
1.676 +    /// \param map A node map storing the supply values.
1.677 +    /// Its \c Value type must be convertible to the \c Value type
1.678 +    /// of the algorithm.
1.679 +    ///
1.680 +    /// \return <tt>(*this)</tt>
1.681 +    template<typename SupplyMap>
1.682 +    CostScaling& supplyMap(const SupplyMap& map) {
1.683 +      for (NodeIt n(_graph); n != INVALID; ++n) {
1.684 +        _supply[_node_id[n]] = map[n];
1.685 +      }
1.686 +      return *this;
1.687 +    }
1.688 +
1.689 +    /// \brief Set single source and target nodes and a supply value.
1.690 +    ///
1.691 +    /// This function sets a single source node and a single target node
1.692 +    /// and the required flow value.
1.693 +    /// If neither this function nor \ref supplyMap() is used before
1.694 +    /// calling \ref run(), the supply of each node will be set to zero.
1.695 +    ///
1.696 +    /// Using this function has the same effect as using \ref supplyMap()
1.697 +    /// with such a map in which \c k is assigned to \c s, \c -k is
1.698 +    /// assigned to \c t and all other nodes have zero supply value.
1.699 +    ///
1.700 +    /// \param s The source node.
1.701 +    /// \param t The target node.
1.702 +    /// \param k The required amount of flow from node \c s to node \c t
1.703 +    /// (i.e. the supply of \c s and the demand of \c t).
1.704 +    ///
1.705 +    /// \return <tt>(*this)</tt>
1.706 +    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
1.707 +      for (int i = 0; i != _res_node_num; ++i) {
1.708 +        _supply[i] = 0;
1.709 +      }
1.710 +      _supply[_node_id[s]] =  k;
1.711 +      _supply[_node_id[t]] = -k;
1.712 +      return *this;
1.713 +    }
1.714 +
1.715 +    /// @}
1.716 +
1.717      /// \name Execution control
1.718 +    /// The algorithm can be executed using \ref run().
1.719
1.720      /// @{
1.721
1.722      /// \brief Run the algorithm.
1.723      ///
1.724 -    /// Run the algorithm.
1.725 +    /// This function runs the algorithm.
1.726 +    /// The paramters can be specified using functions \ref lowerMap(),
1.727 +    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.728 +    /// For example,
1.729 +    /// \code
1.730 +    ///   CostScaling<ListDigraph> cs(graph);
1.731 +    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.732 +    ///     .supplyMap(sup).run();
1.733 +    /// \endcode
1.734 +    ///
1.735 +    /// This function can be called more than once. All the parameters
1.736 +    /// that have been given are kept for the next call, unless
1.737 +    /// \ref reset() is called, thus only the modified parameters
1.738 +    /// have to be set again. See \ref reset() for examples.
1.739 +    /// However the underlying digraph must not be modified after this
1.740 +    /// class have been constructed, since it copies the digraph.
1.741      ///
1.742      /// \param partial_augment By default the algorithm performs
1.743      /// partial augment and relabel operations in the cost scaling
1.744      /// phases. Set this parameter to \c false for using local push and
1.746      ///
1.747 -    /// \return \c true if a feasible flow can be found.
1.748 -    bool run(bool partial_augment = true) {
1.749 -      if (partial_augment) {
1.750 -        return init() && startPartialAugment();
1.751 -      } else {
1.752 -        return init() && startPushRelabel();
1.753 +    /// \return \c INFEASIBLE if no feasible flow exists,
1.754 +    /// \n \c OPTIMAL if the problem has optimal solution
1.755 +    /// (i.e. it is feasible and bounded), and the algorithm has found
1.756 +    /// optimal flow and node potentials (primal and dual solutions),
1.757 +    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
1.758 +    /// and infinite upper bound. It means that the objective function
1.759 +    /// is unbounded on that arc, however note that it could actually be
1.760 +    /// bounded over the feasible flows, but this algroithm cannot handle
1.761 +    /// these cases.
1.762 +    ///
1.763 +    /// \see ProblemType
1.764 +    ProblemType run(bool partial_augment = true) {
1.765 +      ProblemType pt = init();
1.766 +      if (pt != OPTIMAL) return pt;
1.767 +      start(partial_augment);
1.768 +      return OPTIMAL;
1.769 +    }
1.770 +
1.771 +    /// \brief Reset all the parameters that have been given before.
1.772 +    ///
1.773 +    /// This function resets all the paramaters that have been given
1.774 +    /// before using functions \ref lowerMap(), \ref upperMap(),
1.775 +    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.776 +    ///
1.777 +    /// It is useful for multiple run() calls. If this function is not
1.778 +    /// used, all the parameters given before are kept for the next
1.779 +    /// \ref run() call.
1.780 +    /// However the underlying digraph must not be modified after this
1.781 +    /// class have been constructed, since it copies and extends the graph.
1.782 +    ///
1.783 +    /// For example,
1.784 +    /// \code
1.785 +    ///   CostScaling<ListDigraph> cs(graph);
1.786 +    ///
1.787 +    ///   // First run
1.788 +    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.789 +    ///     .supplyMap(sup).run();
1.790 +    ///
1.791 +    ///   // Run again with modified cost map (reset() is not called,
1.792 +    ///   // so only the cost map have to be set again)
1.793 +    ///   cost[e] += 100;
1.794 +    ///   cs.costMap(cost).run();
1.795 +    ///
1.796 +    ///   // Run again from scratch using reset()
1.797 +    ///   // (the lower bounds will be set to zero on all arcs)
1.798 +    ///   cs.reset();
1.799 +    ///   cs.upperMap(capacity).costMap(cost)
1.800 +    ///     .supplyMap(sup).run();
1.801 +    /// \endcode
1.802 +    ///
1.803 +    /// \return <tt>(*this)</tt>
1.804 +    CostScaling& reset() {
1.805 +      for (int i = 0; i != _res_node_num; ++i) {
1.806 +        _supply[i] = 0;
1.807        }
1.808 +      int limit = _first_out[_root];
1.809 +      for (int j = 0; j != limit; ++j) {
1.810 +        _lower[j] = 0;
1.811 +        _upper[j] = INF;
1.812 +        _scost[j] = _forward[j] ? 1 : -1;
1.813 +      }
1.814 +      for (int j = limit; j != _res_arc_num; ++j) {
1.815 +        _lower[j] = 0;
1.816 +        _upper[j] = INF;
1.817 +        _scost[j] = 0;
1.818 +        _scost[_reverse[j]] = 0;
1.819 +      }
1.820 +      _have_lower = false;
1.821 +      return *this;
1.822      }
1.823
1.824      /// @}
1.825
1.826      /// \name Query Functions
1.827 -    /// The result of the algorithm can be obtained using these
1.828 +    /// The results of the algorithm can be obtained using these
1.829      /// functions.\n
1.830 -    /// \ref lemon::CostScaling::run() "run()" must be called before
1.831 -    /// using them.
1.832 +    /// The \ref run() function must be called before using them.
1.833
1.834      /// @{
1.835
1.836 -    /// \brief Return a const reference to the arc map storing the
1.837 -    /// found flow.
1.838 +    /// \brief Return the total cost of the found flow.
1.839      ///
1.840 -    /// Return a const reference to the arc map storing the found flow.
1.841 +    /// This function returns the total cost of the found flow.
1.842 +    /// Its complexity is O(e).
1.843 +    ///
1.844 +    /// \note The return type of the function can be specified as a
1.845 +    /// template parameter. For example,
1.846 +    /// \code
1.847 +    ///   cs.totalCost<double>();
1.848 +    /// \endcode
1.849 +    /// It is useful if the total cost cannot be stored in the \c Cost
1.850 +    /// type of the algorithm, which is the default return type of the
1.851 +    /// function.
1.852      ///
1.853      /// \pre \ref run() must be called before using this function.
1.854 -    const FlowMap& flowMap() const {
1.855 -      return *_flow;
1.856 +    template <typename Number>
1.857 +    Number totalCost() const {
1.858 +      Number c = 0;
1.859 +      for (ArcIt a(_graph); a != INVALID; ++a) {
1.860 +        int i = _arc_idb[a];
1.861 +        c += static_cast<Number>(_res_cap[i]) *
1.862 +             (-static_cast<Number>(_scost[i]));
1.863 +      }
1.864 +      return c;
1.865      }
1.866
1.867 -    /// \brief Return a const reference to the node map storing the
1.868 -    /// found potentials (the dual solution).
1.869 -    ///
1.870 -    /// Return a const reference to the node map storing the found
1.871 -    /// potentials (the dual solution).
1.872 -    ///
1.873 -    /// \pre \ref run() must be called before using this function.
1.874 -    const PotentialMap& potentialMap() const {
1.875 -      return *_potential;
1.876 +#ifndef DOXYGEN
1.877 +    Cost totalCost() const {
1.879      }
1.880 +#endif
1.881
1.882      /// \brief Return the flow on the given arc.
1.883      ///
1.884 -    /// Return the flow on the given arc.
1.885 +    /// This function returns the flow on the given arc.
1.886      ///
1.887      /// \pre \ref run() must be called before using this function.
1.888 -    Capacity flow(const Arc& arc) const {
1.889 -      return (*_flow)[arc];
1.890 +    Value flow(const Arc& a) const {
1.891 +      return _res_cap[_arc_idb[a]];
1.892      }
1.893
1.894 -    /// \brief Return the potential of the given node.
1.895 +    /// \brief Return the flow map (the primal solution).
1.896      ///
1.897 -    /// Return the potential of the given node.
1.898 +    /// This function copies the flow value on each arc into the given
1.899 +    /// map. The \c Value type of the algorithm must be convertible to
1.900 +    /// the \c Value type of the map.
1.901      ///
1.902      /// \pre \ref run() must be called before using this function.
1.903 -    Cost potential(const Node& node) const {
1.904 -      return (*_potential)[node];
1.905 +    template <typename FlowMap>
1.906 +    void flowMap(FlowMap &map) const {
1.907 +      for (ArcIt a(_graph); a != INVALID; ++a) {
1.908 +        map.set(a, _res_cap[_arc_idb[a]]);
1.909 +      }
1.910      }
1.911
1.912 -    /// \brief Return the total cost of the found flow.
1.913 +    /// \brief Return the potential (dual value) of the given node.
1.914      ///
1.915 -    /// Return the total cost of the found flow. The complexity of the
1.916 -    /// function is \f$O(e) \f$.
1.917 +    /// This function returns the potential (dual value) of the
1.918 +    /// given node.
1.919      ///
1.920      /// \pre \ref run() must be called before using this function.
1.921 -    Cost totalCost() const {
1.922 -      Cost c = 0;
1.923 -      for (ArcIt e(_graph); e != INVALID; ++e)
1.924 -        c += (*_flow)[e] * _orig_cost[e];
1.925 -      return c;
1.926 +    Cost potential(const Node& n) const {
1.927 +      return static_cast<Cost>(_pi[_node_id[n]]);
1.928 +    }
1.929 +
1.930 +    /// \brief Return the potential map (the dual solution).
1.931 +    ///
1.932 +    /// This function copies the potential (dual value) of each node
1.933 +    /// into the given map.
1.934 +    /// The \c Cost type of the algorithm must be convertible to the
1.935 +    /// \c Value type of the map.
1.936 +    ///
1.937 +    /// \pre \ref run() must be called before using this function.
1.938 +    template <typename PotentialMap>
1.939 +    void potentialMap(PotentialMap &map) const {
1.940 +      for (NodeIt n(_graph); n != INVALID; ++n) {
1.941 +        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
1.942 +      }
1.943      }
1.944
1.945      /// @}
1.946
1.947    private:
1.948
1.949 -    /// Initialize the algorithm.
1.950 -    bool init() {
1.951 -      if (!_valid_supply) return false;
1.952 -      // The scaling factor
1.953 +    // Initialize the algorithm
1.954 +    ProblemType init() {
1.955 +      if (_res_node_num == 0) return INFEASIBLE;
1.956 +
1.957 +      // Scaling factor
1.958        _alpha = 8;
1.959
1.960 -      // Initialize flow and potential maps
1.961 -      if (!_flow) {
1.962 -        _flow = new FlowMap(_graph);
1.963 -        _local_flow = true;
1.964 +      // Check the sum of supply values
1.965 +      _sum_supply = 0;
1.966 +      for (int i = 0; i != _root; ++i) {
1.967 +        _sum_supply += _supply[i];
1.968        }
1.969 -      if (!_potential) {
1.970 -        _potential = new PotentialMap(_graph);
1.971 -        _local_potential = true;
1.972 +      if (_sum_supply > 0) return INFEASIBLE;
1.973 +
1.974 +
1.975 +      // Initialize vectors
1.976 +      for (int i = 0; i != _res_node_num; ++i) {
1.977 +        _pi[i] = 0;
1.978 +        _excess[i] = _supply[i];
1.979 +      }
1.980 +
1.981 +      // Remove infinite upper bounds and check negative arcs
1.982 +      const Value MAX = std::numeric_limits<Value>::max();
1.983 +      int last_out;
1.984 +      if (_have_lower) {
1.985 +        for (int i = 0; i != _root; ++i) {
1.986 +          last_out = _first_out[i+1];
1.987 +          for (int j = _first_out[i]; j != last_out; ++j) {
1.988 +            if (_forward[j]) {
1.989 +              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
1.990 +              if (c >= MAX) return UNBOUNDED;
1.991 +              _excess[i] -= c;
1.992 +              _excess[_target[j]] += c;
1.993 +            }
1.994 +          }
1.995 +        }
1.996 +      } else {
1.997 +        for (int i = 0; i != _root; ++i) {
1.998 +          last_out = _first_out[i+1];
1.999 +          for (int j = _first_out[i]; j != last_out; ++j) {
1.1000 +            if (_forward[j] && _scost[j] < 0) {
1.1001 +              Value c = _upper[j];
1.1002 +              if (c >= MAX) return UNBOUNDED;
1.1003 +              _excess[i] -= c;
1.1004 +              _excess[_target[j]] += c;
1.1005 +            }
1.1006 +          }
1.1007 +        }
1.1008 +      }
1.1009 +      Value ex, max_cap = 0;
1.1010 +      for (int i = 0; i != _res_node_num; ++i) {
1.1011 +        ex = _excess[i];
1.1012 +        _excess[i] = 0;
1.1013 +        if (ex < 0) max_cap -= ex;
1.1014 +      }
1.1015 +      for (int j = 0; j != _res_arc_num; ++j) {
1.1016 +        if (_upper[j] >= MAX) _upper[j] = max_cap;
1.1017        }
1.1018
1.1019 -      _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
1.1020 -      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
1.1021 +      // Initialize the large cost vector and the epsilon parameter
1.1022 +      _epsilon = 0;
1.1023 +      LargeCost lc;
1.1024 +      for (int i = 0; i != _root; ++i) {
1.1025 +        last_out = _first_out[i+1];
1.1026 +        for (int j = _first_out[i]; j != last_out; ++j) {
1.1027 +          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
1.1028 +          _cost[j] = lc;
1.1029 +          if (lc > _epsilon) _epsilon = lc;
1.1030 +        }
1.1031 +      }
1.1032 +      _epsilon /= _alpha;
1.1033
1.1034 -      // Initialize the scaled cost map and the epsilon parameter
1.1035 -      Cost max_cost = 0;
1.1036 -      int node_num = countNodes(_graph);
1.1037 -      for (ArcIt e(_graph); e != INVALID; ++e) {
1.1038 -        _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
1.1039 -        if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
1.1040 +      // Initialize maps for Circulation and remove non-zero lower bounds
1.1041 +      ConstMap<Arc, Value> low(0);
1.1042 +      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
1.1043 +      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
1.1044 +      ValueArcMap cap(_graph), flow(_graph);
1.1045 +      ValueNodeMap sup(_graph);
1.1046 +      for (NodeIt n(_graph); n != INVALID; ++n) {
1.1047 +        sup[n] = _supply[_node_id[n]];
1.1048        }
1.1049 -      _epsilon = max_cost * node_num;
1.1050 +      if (_have_lower) {
1.1051 +        for (ArcIt a(_graph); a != INVALID; ++a) {
1.1052 +          int j = _arc_idf[a];
1.1053 +          Value c = _lower[j];
1.1054 +          cap[a] = _upper[j] - c;
1.1055 +          sup[_graph.source(a)] -= c;
1.1056 +          sup[_graph.target(a)] += c;
1.1057 +        }
1.1058 +      } else {
1.1059 +        for (ArcIt a(_graph); a != INVALID; ++a) {
1.1060 +          cap[a] = _upper[_arc_idf[a]];
1.1061 +        }
1.1062 +      }
1.1063
1.1064        // Find a feasible flow using Circulation
1.1065 -      Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
1.1066 -                   SupplyMap >
1.1067 -        circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
1.1068 -                     _supply );
1.1069 -      return circulation.flowMap(*_flow).run();
1.1070 +      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
1.1071 +        circ(_graph, low, cap, sup);
1.1072 +      if (!circ.flowMap(flow).run()) return INFEASIBLE;
1.1073 +
1.1074 +      // Set residual capacities and handle GEQ supply type
1.1075 +      if (_sum_supply < 0) {
1.1076 +        for (ArcIt a(_graph); a != INVALID; ++a) {
1.1077 +          Value fa = flow[a];
1.1078 +          _res_cap[_arc_idf[a]] = cap[a] - fa;
1.1079 +          _res_cap[_arc_idb[a]] = fa;
1.1080 +          sup[_graph.source(a)] -= fa;
1.1081 +          sup[_graph.target(a)] += fa;
1.1082 +        }
1.1083 +        for (NodeIt n(_graph); n != INVALID; ++n) {
1.1084 +          _excess[_node_id[n]] = sup[n];
1.1085 +        }
1.1086 +        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1087 +          int u = _target[a];
1.1088 +          int ra = _reverse[a];
1.1089 +          _res_cap[a] = -_sum_supply + 1;
1.1090 +          _res_cap[ra] = -_excess[u];
1.1091 +          _cost[a] = 0;
1.1092 +          _cost[ra] = 0;
1.1093 +          _excess[u] = 0;
1.1094 +        }
1.1095 +      } else {
1.1096 +        for (ArcIt a(_graph); a != INVALID; ++a) {
1.1097 +          Value fa = flow[a];
1.1098 +          _res_cap[_arc_idf[a]] = cap[a] - fa;
1.1099 +          _res_cap[_arc_idb[a]] = fa;
1.1100 +        }
1.1101 +        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1102 +          int ra = _reverse[a];
1.1103 +          _res_cap[a] = 1;
1.1104 +          _res_cap[ra] = 0;
1.1105 +          _cost[a] = 0;
1.1106 +          _cost[ra] = 0;
1.1107 +        }
1.1108 +      }
1.1109 +
1.1110 +      return OPTIMAL;
1.1111 +    }
1.1112 +
1.1113 +    // Execute the algorithm and transform the results
1.1114 +    void start(bool partial_augment) {
1.1115 +      // Execute the algorithm
1.1116 +      if (partial_augment) {
1.1117 +        startPartialAugment();
1.1118 +      } else {
1.1119 +        startPushRelabel();
1.1120 +      }
1.1121 +
1.1122 +      // Compute node potentials for the original costs
1.1123 +      _arc_vec.clear();
1.1124 +      _cost_vec.clear();
1.1125 +      for (int j = 0; j != _res_arc_num; ++j) {
1.1126 +        if (_res_cap[j] > 0) {
1.1127 +          _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1128 +          _cost_vec.push_back(_scost[j]);
1.1129 +        }
1.1130 +      }
1.1131 +      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1132 +
1.1133 +      typename BellmanFord<StaticDigraph, LargeCostArcMap>
1.1134 +        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
1.1135 +      bf.distMap(_pi_map);
1.1136 +      bf.init(0);
1.1137 +      bf.start();
1.1138 +
1.1139 +      // Handle non-zero lower bounds
1.1140 +      if (_have_lower) {
1.1141 +        int limit = _first_out[_root];
1.1142 +        for (int j = 0; j != limit; ++j) {
1.1143 +          if (!_forward[j]) _res_cap[j] += _lower[j];
1.1144 +        }
1.1145 +      }
1.1146      }
1.1147
1.1148      /// Execute the algorithm performing partial augmentation and
1.1149 -    /// relabel operations.
1.1150 -    bool startPartialAugment() {
1.1151 +    /// relabel operations
1.1152 +    void startPartialAugment() {
1.1153        // Paramters for heuristics
1.1154 -//      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1155 -//      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
1.1156 +      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1157 +      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
1.1158        // Maximum augment path length
1.1159        const int MAX_PATH_LENGTH = 4;
1.1160
1.1161 -      // Variables
1.1162 -      typename Digraph::template NodeMap<Arc> pred_arc(_graph);
1.1163 -      typename Digraph::template NodeMap<bool> forward(_graph);
1.1164 -      typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
1.1165 -      typename Digraph::template NodeMap<InArcIt> next_in(_graph);
1.1166 -      typename Digraph::template NodeMap<bool> next_dir(_graph);
1.1167 -      std::deque<Node> active_nodes;
1.1168 -      std::vector<Node> path_nodes;
1.1169 -
1.1170 -//      int node_num = countNodes(_graph);
1.1171 +      // Perform cost scaling phases
1.1172 +      IntVector pred_arc(_res_node_num);
1.1173 +      std::vector<int> path_nodes;
1.1174        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.1175                                          1 : _epsilon / _alpha )
1.1176        {
1.1177 -/*
1.1178          // "Early Termination" heuristic: use Bellman-Ford algorithm
1.1179          // to check if the current flow is optimal
1.1180          if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.1181 -          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.1182 -          ShiftCostMap shift_cost(_res_cost, 1);
1.1183 -          BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.1184 +          _arc_vec.clear();
1.1185 +          _cost_vec.clear();
1.1186 +          for (int j = 0; j != _res_arc_num; ++j) {
1.1187 +            if (_res_cap[j] > 0) {
1.1188 +              _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1189 +              _cost_vec.push_back(_cost[j] + 1);
1.1190 +            }
1.1191 +          }
1.1192 +          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1193 +
1.1194 +          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.1195            bf.init(0);
1.1196            bool done = false;
1.1197 -          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.1198 +          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1.1199            for (int i = 0; i < K && !done; ++i)
1.1200              done = bf.processNextWeakRound();
1.1201            if (done) break;
1.1202          }
1.1203 -*/
1.1204 +
1.1205          // Saturate arcs not satisfying the optimality condition
1.1206 -        Capacity delta;
1.1207 -        for (ArcIt e(_graph); e != INVALID; ++e) {
1.1208 -          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1209 -            delta = _capacity[e] - (*_flow)[e];
1.1210 -            _excess[_graph.source(e)] -= delta;
1.1211 -            _excess[_graph.target(e)] += delta;
1.1212 -            (*_flow)[e] = _capacity[e];
1.1213 -          }
1.1214 -          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1215 -            _excess[_graph.target(e)] -= (*_flow)[e];
1.1216 -            _excess[_graph.source(e)] += (*_flow)[e];
1.1217 -            (*_flow)[e] = 0;
1.1218 +        for (int a = 0; a != _res_arc_num; ++a) {
1.1219 +          if (_res_cap[a] > 0 &&
1.1220 +              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1221 +            Value delta = _res_cap[a];
1.1222 +            _excess[_source[a]] -= delta;
1.1223 +            _excess[_target[a]] += delta;
1.1224 +            _res_cap[a] = 0;
1.1225 +            _res_cap[_reverse[a]] += delta;
1.1226            }
1.1227          }
1.1228 -
1.1229 +
1.1230          // Find active nodes (i.e. nodes with positive excess)
1.1231 -        for (NodeIt n(_graph); n != INVALID; ++n) {
1.1232 -          if (_excess[n] > 0) active_nodes.push_back(n);
1.1233 +        for (int u = 0; u != _res_node_num; ++u) {
1.1234 +          if (_excess[u] > 0) _active_nodes.push_back(u);
1.1235          }
1.1236
1.1237 -        // Initialize the next arc maps
1.1238 -        for (NodeIt n(_graph); n != INVALID; ++n) {
1.1239 -          next_out[n] = OutArcIt(_graph, n);
1.1240 -          next_in[n] = InArcIt(_graph, n);
1.1241 -          next_dir[n] = true;
1.1242 +        // Initialize the next arcs
1.1243 +        for (int u = 0; u != _res_node_num; ++u) {
1.1244 +          _next_out[u] = _first_out[u];
1.1245          }
1.1246
1.1247          // Perform partial augment and relabel operations
1.1248 -        while (active_nodes.size() > 0) {
1.1249 +        while (true) {
1.1250            // Select an active node (FIFO selection)
1.1251 -          if (_excess[active_nodes[0]] <= 0) {
1.1252 -            active_nodes.pop_front();
1.1253 -            continue;
1.1254 +          while (_active_nodes.size() > 0 &&
1.1255 +                 _excess[_active_nodes.front()] <= 0) {
1.1256 +            _active_nodes.pop_front();
1.1257            }
1.1258 -          Node start = active_nodes[0];
1.1259 +          if (_active_nodes.size() == 0) break;
1.1260 +          int start = _active_nodes.front();
1.1261            path_nodes.clear();
1.1262            path_nodes.push_back(start);
1.1263
1.1264            // Find an augmenting path from the start node
1.1265 -          Node u, tip = start;
1.1266 -          LCost min_red_cost;
1.1267 -          while ( _excess[tip] >= 0 &&
1.1268 -                  int(path_nodes.size()) <= MAX_PATH_LENGTH )
1.1269 -          {
1.1270 -            if (next_dir[tip]) {
1.1271 -              for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
1.1272 -                if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1273 -                  u = _graph.target(e);
1.1274 -                  pred_arc[u] = e;
1.1275 -                  forward[u] = true;
1.1276 -                  next_out[tip] = e;
1.1277 -                  tip = u;
1.1278 -                  path_nodes.push_back(tip);
1.1279 -                  goto next_step;
1.1280 -                }
1.1281 -              }
1.1282 -              next_dir[tip] = false;
1.1283 -            }
1.1284 -            for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
1.1285 -              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1286 -                u = _graph.source(e);
1.1287 -                pred_arc[u] = e;
1.1288 -                forward[u] = false;
1.1289 -                next_in[tip] = e;
1.1290 +          int tip = start;
1.1291 +          while (_excess[tip] >= 0 &&
1.1292 +                 int(path_nodes.size()) <= MAX_PATH_LENGTH) {
1.1293 +            int u;
1.1294 +            LargeCost min_red_cost, rc;
1.1295 +            int last_out = _sum_supply < 0 ?
1.1296 +              _first_out[tip+1] : _first_out[tip+1] - 1;
1.1297 +            for (int a = _next_out[tip]; a != last_out; ++a) {
1.1298 +              if (_res_cap[a] > 0 &&
1.1299 +                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1300 +                u = _target[a];
1.1301 +                pred_arc[u] = a;
1.1302 +                _next_out[tip] = a;
1.1303                  tip = u;
1.1304                  path_nodes.push_back(tip);
1.1305                  goto next_step;
1.1306 @@ -579,266 +943,186 @@
1.1307              }
1.1308
1.1309              // Relabel tip node
1.1310 -            min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.1311 -            for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
1.1312 -              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.1313 -                   (*_red_cost)[oe] < min_red_cost )
1.1314 -                min_red_cost = (*_red_cost)[oe];
1.1315 +            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1.1316 +            for (int a = _first_out[tip]; a != last_out; ++a) {
1.1317 +              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1.1318 +              if (_res_cap[a] > 0 && rc < min_red_cost) {
1.1319 +                min_red_cost = rc;
1.1320 +              }
1.1321              }
1.1322 -            for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
1.1323 -              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.1324 -                min_red_cost = -(*_red_cost)[ie];
1.1325 -            }
1.1326 -            (*_potential)[tip] -= min_red_cost + _epsilon;
1.1327 +            _pi[tip] -= min_red_cost + _epsilon;
1.1328
1.1329 -            // Reset the next arc maps
1.1330 -            next_out[tip] = OutArcIt(_graph, tip);
1.1331 -            next_in[tip] = InArcIt(_graph, tip);
1.1332 -            next_dir[tip] = true;
1.1333 +            // Reset the next arc of tip
1.1334 +            _next_out[tip] = _first_out[tip];
1.1335
1.1336              // Step back
1.1337              if (tip != start) {
1.1338                path_nodes.pop_back();
1.1339 -              tip = path_nodes[path_nodes.size()-1];
1.1340 +              tip = path_nodes.back();
1.1341              }
1.1342
1.1343 -          next_step:
1.1344 -            continue;
1.1345 +          next_step: ;
1.1346            }
1.1347
1.1348            // Augment along the found path (as much flow as possible)
1.1349 -          Capacity delta;
1.1350 +          Value delta;
1.1351 +          int u, v = path_nodes.front(), pa;
1.1352            for (int i = 1; i < int(path_nodes.size()); ++i) {
1.1353 -            u = path_nodes[i];
1.1354 -            delta = forward[u] ?
1.1355 -              _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
1.1356 -              (*_flow)[pred_arc[u]];
1.1357 -            delta = std::min(delta, _excess[path_nodes[i-1]]);
1.1358 -            (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
1.1359 -            _excess[path_nodes[i-1]] -= delta;
1.1360 -            _excess[u] += delta;
1.1361 -            if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
1.1362 +            u = v;
1.1363 +            v = path_nodes[i];
1.1364 +            pa = pred_arc[v];
1.1365 +            delta = std::min(_res_cap[pa], _excess[u]);
1.1366 +            _res_cap[pa] -= delta;
1.1367 +            _res_cap[_reverse[pa]] += delta;
1.1368 +            _excess[u] -= delta;
1.1369 +            _excess[v] += delta;
1.1370 +            if (_excess[v] > 0 && _excess[v] <= delta)
1.1371 +              _active_nodes.push_back(v);
1.1372            }
1.1373          }
1.1374        }
1.1375 -
1.1376 -      // Compute node potentials for the original costs
1.1377 -      ResidualCostMap<CostMap> res_cost(_orig_cost);
1.1378 -      BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
1.1379 -        bf(*_res_graph, res_cost);
1.1380 -      bf.init(0); bf.start();
1.1381 -      for (NodeIt n(_graph); n != INVALID; ++n)
1.1382 -        (*_potential)[n] = bf.dist(n);
1.1383 -
1.1384 -      // Handle non-zero lower bounds
1.1385 -      if (_lower) {
1.1386 -        for (ArcIt e(_graph); e != INVALID; ++e)
1.1387 -          (*_flow)[e] += (*_lower)[e];
1.1388 -      }
1.1389 -      return true;
1.1390      }
1.1391
1.1392 -    /// Execute the algorithm performing push and relabel operations.
1.1393 -    bool startPushRelabel() {
1.1394 +    /// Execute the algorithm performing push and relabel operations
1.1395 +    void startPushRelabel() {
1.1396        // Paramters for heuristics
1.1397 -//      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1398 -//      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
1.1399 +      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1400 +      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
1.1401
1.1402 -      typename Digraph::template NodeMap<bool> hyper(_graph, false);
1.1403 -      typename Digraph::template NodeMap<Arc> pred_arc(_graph);
1.1404 -      typename Digraph::template NodeMap<bool> forward(_graph);
1.1405 -      typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
1.1406 -      typename Digraph::template NodeMap<InArcIt> next_in(_graph);
1.1407 -      typename Digraph::template NodeMap<bool> next_dir(_graph);
1.1408 -      std::deque<Node> active_nodes;
1.1409 -
1.1410 -//      int node_num = countNodes(_graph);
1.1411 +      // Perform cost scaling phases
1.1412 +      BoolVector hyper(_res_node_num, false);
1.1413        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.1414                                          1 : _epsilon / _alpha )
1.1415        {
1.1416 -/*
1.1417          // "Early Termination" heuristic: use Bellman-Ford algorithm
1.1418          // to check if the current flow is optimal
1.1419          if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.1420 -          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.1421 -          ShiftCostMap shift_cost(_res_cost, 1);
1.1422 -          BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.1423 +          _arc_vec.clear();
1.1424 +          _cost_vec.clear();
1.1425 +          for (int j = 0; j != _res_arc_num; ++j) {
1.1426 +            if (_res_cap[j] > 0) {
1.1427 +              _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1428 +              _cost_vec.push_back(_cost[j] + 1);
1.1429 +            }
1.1430 +          }
1.1431 +          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1432 +
1.1433 +          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.1434            bf.init(0);
1.1435            bool done = false;
1.1436 -          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.1437 +          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1.1438            for (int i = 0; i < K && !done; ++i)
1.1439              done = bf.processNextWeakRound();
1.1440            if (done) break;
1.1441          }
1.1442 -*/
1.1443
1.1444          // Saturate arcs not satisfying the optimality condition
1.1445 -        Capacity delta;
1.1446 -        for (ArcIt e(_graph); e != INVALID; ++e) {
1.1447 -          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1448 -            delta = _capacity[e] - (*_flow)[e];
1.1449 -            _excess[_graph.source(e)] -= delta;
1.1450 -            _excess[_graph.target(e)] += delta;
1.1451 -            (*_flow)[e] = _capacity[e];
1.1452 -          }
1.1453 -          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1454 -            _excess[_graph.target(e)] -= (*_flow)[e];
1.1455 -            _excess[_graph.source(e)] += (*_flow)[e];
1.1456 -            (*_flow)[e] = 0;
1.1457 +        for (int a = 0; a != _res_arc_num; ++a) {
1.1458 +          if (_res_cap[a] > 0 &&
1.1459 +              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1460 +            Value delta = _res_cap[a];
1.1461 +            _excess[_source[a]] -= delta;
1.1462 +            _excess[_target[a]] += delta;
1.1463 +            _res_cap[a] = 0;
1.1464 +            _res_cap[_reverse[a]] += delta;
1.1465            }
1.1466          }
1.1467
1.1468          // Find active nodes (i.e. nodes with positive excess)
1.1469 -        for (NodeIt n(_graph); n != INVALID; ++n) {
1.1470 -          if (_excess[n] > 0) active_nodes.push_back(n);
1.1471 +        for (int u = 0; u != _res_node_num; ++u) {
1.1472 +          if (_excess[u] > 0) _active_nodes.push_back(u);
1.1473          }
1.1474
1.1475 -        // Initialize the next arc maps
1.1476 -        for (NodeIt n(_graph); n != INVALID; ++n) {
1.1477 -          next_out[n] = OutArcIt(_graph, n);
1.1478 -          next_in[n] = InArcIt(_graph, n);
1.1479 -          next_dir[n] = true;
1.1480 +        // Initialize the next arcs
1.1481 +        for (int u = 0; u != _res_node_num; ++u) {
1.1482 +          _next_out[u] = _first_out[u];
1.1483          }
1.1484
1.1485          // Perform push and relabel operations
1.1486 -        while (active_nodes.size() > 0) {
1.1487 +        while (_active_nodes.size() > 0) {
1.1488 +          LargeCost min_red_cost, rc;
1.1489 +          Value delta;
1.1490 +          int n, t, a, last_out = _res_arc_num;
1.1491 +
1.1492            // Select an active node (FIFO selection)
1.1493 -          Node n = active_nodes[0], t;
1.1494 -          bool relabel_enabled = true;
1.1495 +        next_node:
1.1496 +          n = _active_nodes.front();
1.1497 +          last_out = _sum_supply < 0 ?
1.1498 +            _first_out[n+1] : _first_out[n+1] - 1;
1.1499
1.1500            // Perform push operations if there are admissible arcs
1.1501 -          if (_excess[n] > 0 && next_dir[n]) {
1.1502 -            OutArcIt e = next_out[n];
1.1503 -            for ( ; e != INVALID; ++e) {
1.1504 -              if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1505 -                delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
1.1506 -                t = _graph.target(e);
1.1507 +          if (_excess[n] > 0) {
1.1508 +            for (a = _next_out[n]; a != last_out; ++a) {
1.1509 +              if (_res_cap[a] > 0 &&
1.1510 +                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1511 +                delta = std::min(_res_cap[a], _excess[n]);
1.1512 +                t = _target[a];
1.1513
1.1515 -                Capacity ahead = -_excess[t];
1.1516 -                for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
1.1517 -                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.1518 -                    ahead += _capacity[oe] - (*_flow)[oe];
1.1519 -                }
1.1520 -                for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
1.1521 -                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.1523 +                Value ahead = -_excess[t];
1.1524 +                int last_out_t = _sum_supply < 0 ?
1.1525 +                  _first_out[t+1] : _first_out[t+1] - 1;
1.1526 +                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1.1527 +                  if (_res_cap[ta] > 0 &&
1.1528 +                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
1.1530 +                  if (ahead >= delta) break;
1.1531                  }
1.1533
1.1534                  // Push flow along the arc
1.1535                  if (ahead < delta) {
1.1541 -                  active_nodes.push_front(t);
1.1542 +                  _active_nodes.push_front(t);
1.1543                    hyper[t] = true;
1.1544 -                  relabel_enabled = false;
1.1545 -                  break;
1.1546 +                  _next_out[n] = a;
1.1547 +                  goto next_node;
1.1548                  } else {
1.1549 -                  (*_flow)[e] += delta;
1.1550 +                  _res_cap[a] -= delta;
1.1551 +                  _res_cap[_reverse[a]] += delta;
1.1552                    _excess[n] -= delta;
1.1553                    _excess[t] += delta;
1.1554                    if (_excess[t] > 0 && _excess[t] <= delta)
1.1555 -                    active_nodes.push_back(t);
1.1556 +                    _active_nodes.push_back(t);
1.1557                  }
1.1558
1.1559 -                if (_excess[n] == 0) break;
1.1560 +                if (_excess[n] == 0) {
1.1561 +                  _next_out[n] = a;
1.1562 +                  goto remove_nodes;
1.1563 +                }
1.1564                }
1.1565              }
1.1566 -            if (e != INVALID) {
1.1567 -              next_out[n] = e;
1.1568 -            } else {
1.1569 -              next_dir[n] = false;
1.1570 -            }
1.1571 -          }
1.1572 -
1.1573 -          if (_excess[n] > 0 && !next_dir[n]) {
1.1574 -            InArcIt e = next_in[n];
1.1575 -            for ( ; e != INVALID; ++e) {
1.1576 -              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1577 -                delta = std::min((*_flow)[e], _excess[n]);
1.1578 -                t = _graph.source(e);
1.1579 -
1.1581 -                Capacity ahead = -_excess[t];
1.1582 -                for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
1.1583 -                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.1584 -                    ahead += _capacity[oe] - (*_flow)[oe];
1.1585 -                }
1.1586 -                for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
1.1587 -                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.1589 -                }
1.1591 -
1.1592 -                // Push flow along the arc
1.1593 -                if (ahead < delta) {
1.1597 -                  active_nodes.push_front(t);
1.1598 -                  hyper[t] = true;
1.1599 -                  relabel_enabled = false;
1.1600 -                  break;
1.1601 -                } else {
1.1602 -                  (*_flow)[e] -= delta;
1.1603 -                  _excess[n] -= delta;
1.1604 -                  _excess[t] += delta;
1.1605 -                  if (_excess[t] > 0 && _excess[t] <= delta)
1.1606 -                    active_nodes.push_back(t);
1.1607 -                }
1.1608 -
1.1609 -                if (_excess[n] == 0) break;
1.1610 -              }
1.1611 -            }
1.1612 -            next_in[n] = e;
1.1613 +            _next_out[n] = a;
1.1614            }
1.1615
1.1616            // Relabel the node if it is still active (or hyper)
1.1617 -          if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
1.1618 -            LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.1619 -            for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
1.1620 -              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.1621 -                   (*_red_cost)[oe] < min_red_cost )
1.1622 -                min_red_cost = (*_red_cost)[oe];
1.1623 +          if (_excess[n] > 0 || hyper[n]) {
1.1624 +            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1.1625 +            for (int a = _first_out[n]; a != last_out; ++a) {
1.1626 +              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1.1627 +              if (_res_cap[a] > 0 && rc < min_red_cost) {
1.1628 +                min_red_cost = rc;
1.1629 +              }
1.1630              }
1.1631 -            for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
1.1632 -              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.1633 -                min_red_cost = -(*_red_cost)[ie];
1.1634 -            }
1.1635 -            (*_potential)[n] -= min_red_cost + _epsilon;
1.1636 +            _pi[n] -= min_red_cost + _epsilon;
1.1637              hyper[n] = false;
1.1638
1.1639 -            // Reset the next arc maps
1.1640 -            next_out[n] = OutArcIt(_graph, n);
1.1641 -            next_in[n] = InArcIt(_graph, n);
1.1642 -            next_dir[n] = true;
1.1643 +            // Reset the next arc
1.1644 +            _next_out[n] = _first_out[n];
1.1645            }
1.1646 -
1.1647 +
1.1648            // Remove nodes that are not active nor hyper
1.1649 -          while ( active_nodes.size() > 0 &&
1.1650 -                  _excess[active_nodes[0]] <= 0 &&
1.1651 -                  !hyper[active_nodes[0]] ) {
1.1652 -            active_nodes.pop_front();
1.1653 +        remove_nodes:
1.1654 +          while ( _active_nodes.size() > 0 &&
1.1655 +                  _excess[_active_nodes.front()] <= 0 &&
1.1656 +                  !hyper[_active_nodes.front()] ) {
1.1657 +            _active_nodes.pop_front();
1.1658            }
1.1659          }
1.1660        }
1.1661 -
1.1662 -      // Compute node potentials for the original costs
1.1663 -      ResidualCostMap<CostMap> res_cost(_orig_cost);
1.1664 -      BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
1.1665 -        bf(*_res_graph, res_cost);
1.1666 -      bf.init(0); bf.start();
1.1667 -      for (NodeIt n(_graph); n != INVALID; ++n)
1.1668 -        (*_potential)[n] = bf.dist(n);
1.1669 -
1.1670 -      // Handle non-zero lower bounds
1.1671 -      if (_lower) {
1.1672 -        for (ArcIt e(_graph); e != INVALID; ++e)
1.1673 -          (*_flow)[e] += (*_lower)[e];
1.1674 -      }
1.1675 -      return true;
1.1676      }
1.1677
1.1678    }; //class CostScaling