Context Navigation

← Previous Changeset
Next Changeset →

Changeset 809:22bb98ca0101 in lemon-main

Timestamp:

11/12/09 23:30:45 (15 years ago)

Author:

Peter Kovacs <kpeter@…>

Branch:

default

Phase:

public

Rebase:

35333461653237666330353538396562623465656466303531646265343536366333386130366130

Message:

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.

File:

: 1 edited

lemon/cost_scaling.h (modified) (6 diffs)

Legend:

: Unmodified
: Added
: Removed

lemon/cost_scaling.h

-                      r808
+                      r809
 #include <lemon/maps.h>
 #include <lemon/math.h>
 #include <lemon/adaptors.h>
+#include <lemon/static_graph.h>
 #include <lemon/circulation.h>
 #include <lemon/bellman_ford.h>
 …
 namespace lemon {
+  /// \brief Default traits class of CostScaling algorithm.
+  ///
+  /// Default traits class of CostScaling algorithm.
+  /// \tparam GR Digraph type.
+  /// \tparam V The value type used for flow amounts, capacity bounds
+  /// and supply values. By default it is \c int.
+  /// \tparam C The value type used for costs and potentials.
+  /// By default it is the same as \c V.
+#ifdef DOXYGEN
+  template <typename GR, typename V = int, typename C = V>
+#else
+  template < typename GR, typename V = int, typename C = V,
+             bool integer = std::numeric_limits<C>::is_integer >
+#endif
+  struct CostScalingDefaultTraits
+  {
+    /// The type of the digraph
+    typedef GR Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef V Value;
+    /// The type of the arc costs
+    typedef C Cost;
+    /// \brief The large cost type used for internal computations
+    ///
+    /// The large cost type used for internal computations.
+    /// It is \c long \c long if the \c Cost type is integer,
+    /// otherwise it is \c double.
+    /// \c Cost must be convertible to \c LargeCost.
+    typedef double LargeCost;
+  };
+  // Default traits class for integer cost types
+  template <typename GR, typename V, typename C>
+  struct CostScalingDefaultTraits<GR, V, C, true>
+  {
+    typedef GR Digraph;
+    typedef V Value;
+    typedef C Cost;
+#ifdef LEMON_HAVE_LONG_LONG
+    typedef long long LargeCost;
+#else
+    typedef long LargeCost;
+#endif
+  };
   /// \addtogroup min_cost_flow_algs
   /// @{
   /// \brief Implementation of the cost scaling algorithm for finding a
   /// minimum cost flow.
+  /// \brief Implementation of the Cost Scaling algorithm for
+  /// finding a \ref min_cost_flow "minimum cost flow".
   ///
+  /// \ref CostScaling implements the cost scaling algorithm performing
+  /// augment/push and relabel operations for finding a minimum cost
+  /// flow.
+  /// \ref CostScaling implements a cost scaling algorithm that performs
+  /// push/augment and relabel operations for finding a minimum cost
+  /// flow. It is an efficient primal-dual solution method, which
+  /// can be viewed as the generalization of the \ref Preflow
+  /// "preflow push-relabel" algorithm for the maximum flow problem.
   ///
+  /// \tparam Digraph The digraph type the algorithm runs on.
+  /// \tparam LowerMap The type of the lower bound map.
+  /// \tparam CapacityMap The type of the capacity (upper bound) map.
+  /// \tparam CostMap The type of the cost (length) map.
+  /// \tparam SupplyMap The type of the supply map.
+  /// Most of the parameters of the problem (except for the digraph)
+  /// can be given using separate functions, and the algorithm can be
+  /// executed using the \ref run() function. If some parameters are not
+  /// specified, then default values will be used.
   ///
   /// \warning
   /// - Arc capacities and costs should be \e non-negative \e integers.
   /// - Supply values should be \e signed \e integers.
   /// - The value types of the maps should be convertible to each other.
   /// - \c CostMap::Value must be signed type.
+  /// \tparam GR The digraph type the algorithm runs on.
+  /// \tparam V The value type used for flow amounts, capacity bounds
+  /// and supply values in the algorithm. By default it is \c int.
+  /// \tparam C The value type used for costs and potentials in the
+  /// algorithm. By default it is the same as \c V.
   ///
+  /// \note Arc costs are multiplied with the number of nodes during
+  /// the algorithm so overflow problems may arise more easily than with
+  /// other minimum cost flow algorithms.
+  /// If it is available, <tt>long long int</tt> type is used instead of
+  /// <tt>long int</tt> in the inside computations.
+  ///
+  /// \author Peter Kovacs
+  template < typename Digraph,
+             typename LowerMap = typename Digraph::template ArcMap<int>,
+             typename CapacityMap = typename Digraph::template ArcMap<int>,
+             typename CostMap = typename Digraph::template ArcMap<int>,
+             typename SupplyMap = typename Digraph::template NodeMap<int> >
+  /// \warning Both value types must be signed and all input data must
+  /// be integer.
+  /// \warning This algorithm does not support negative costs for such
+  /// arcs that have infinite upper bound.
+#ifdef DOXYGEN
+  template <typename GR, typename V, typename C, typename TR>
+#else
+  template < typename GR, typename V = int, typename C = V,
+             typename TR = CostScalingDefaultTraits<GR, V, C> >
+#endif
   class CostScaling
+  {
-    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
-    typedef typename CapacityMap::Value Capacity;
-    typedef typename CostMap::Value Cost;
-    typedef typename SupplyMap::Value Supply;
-    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
-    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
-    typedef ResidualDigraph< const Digraph,
-                             CapacityArcMap, CapacityArcMap > ResDigraph;
-    typedef typename ResDigraph::Arc ResArc;
-#if defined __GNUC__ && !defined __STRICT_ANSI__
-    typedef long long int LCost;
-#else
-    typedef long int LCost;
-#endif
-    typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
   public:
+    /// The type of the flow map.
+    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
+    /// The type of the potential map.
+    typedef typename Digraph::template NodeMap<LCost> PotentialMap;
+    /// The type of the digraph
+    typedef typename TR::Digraph Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef typename TR::Value Value;
+    /// The type of the arc costs
+    typedef typename TR::Cost Cost;
+    /// \brief The large cost type
+    ///
+    /// The large cost type used for internal computations.
+    /// Using the \ref CostScalingDefaultTraits "default traits class",
+    /// it is \c long \c long if the \c Cost type is integer,
+    /// otherwise it is \c double.
+    typedef typename TR::LargeCost LargeCost;
+    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
+    typedef TR Traits;
+  public:
+    /// \brief Problem type constants for the \c run() function.
+    ///
+    /// Enum type containing the problem type constants that can be
+    /// returned by the \ref run() function of the algorithm.
+    enum ProblemType {
+      /// The problem has no feasible solution (flow).
+      INFEASIBLE,
+      /// The problem has optimal solution (i.e. it is feasible and
+      /// bounded), and the algorithm has found optimal flow and node
+      /// potentials (primal and dual solutions).
+      OPTIMAL,
+      /// The digraph contains an arc of negative cost and infinite
+      /// upper bound. It means that the objective function is unbounded
+      /// on that arc, however note that it could actually be bounded
+      /// over the feasible flows, but this algroithm cannot handle
+      /// these cases.
+      UNBOUNDED
+    };
   private:
+    /// \brief Map adaptor class for handling residual arc costs.
+    ///
+    /// Map adaptor class for handling residual arc costs.
+    template <typename Map>
+    class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
+    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
+    typedef std::vector<int> IntVector;
+    typedef std::vector<char> BoolVector;
+    typedef std::vector<Value> ValueVector;
+    typedef std::vector<Cost> CostVector;
+    typedef std::vector<LargeCost> LargeCostVector;
+  private:
+    template <typename KT, typename VT>
+    class VectorMap {
+    public:
+      typedef KT Key;
+      typedef VT Value;
+      VectorMap(std::vector<Value>& v) : _v(v) {}
+      const Value& operator[](const Key& key) const {
+        return _v[StaticDigraph::id(key)];
+      }
+      Value& operator[](const Key& key) {
+        return _v[StaticDigraph::id(key)];
+      }
+      void set(const Key& key, const Value& val) {
+        _v[StaticDigraph::id(key)] = val;
+      }
+    private:
+      std::vector<Value>& _v;
+    };
+    typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
+    typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
+  private:
+    // Data related to the underlying digraph
+    const GR &_graph;
+    int _node_num;
+    int _arc_num;
+    int _res_node_num;
+    int _res_arc_num;
+    int _root;
+    // Parameters of the problem
+    bool _have_lower;
+    Value _sum_supply;
+    // Data structures for storing the digraph
+    IntNodeMap _node_id;
+    IntArcMap _arc_idf;
+    IntArcMap _arc_idb;
+    IntVector _first_out;
+    BoolVector _forward;
+    IntVector _source;
+    IntVector _target;
+    IntVector _reverse;
+    // Node and arc data
+    ValueVector _lower;
+    ValueVector _upper;
+    CostVector _scost;
+    ValueVector _supply;
+    ValueVector _res_cap;
+    LargeCostVector _cost;
+    LargeCostVector _pi;
+    ValueVector _excess;
+    IntVector _next_out;
+    std::deque<int> _active_nodes;
+    // Data for scaling
+    LargeCost _epsilon;
+    int _alpha;
+    // Data for a StaticDigraph structure
+    typedef std::pair<int, int> IntPair;
+    StaticDigraph _sgr;
+    std::vector<IntPair> _arc_vec;
+    std::vector<LargeCost> _cost_vec;
+    LargeCostArcMap _cost_map;
+    LargeCostNodeMap _pi_map;
+  public:
+    /// \brief Constant for infinite upper bounds (capacities).
+    ///
+    /// Constant for infinite upper bounds (capacities).
+    /// It is \c std::numeric_limits<Value>::infinity() if available,
+    /// \c std::numeric_limits<Value>::max() otherwise.
+    const Value INF;
+  public:
+    /// \name Named Template Parameters
+    /// @{
+    template <typename T>
+    struct SetLargeCostTraits : public Traits {
+      typedef T LargeCost;
+    };
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// \c LargeCost type.
+    ///
+    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
+    /// type, which is used for internal computations in the algorithm.
+    /// \c Cost must be convertible to \c LargeCost.
+    template <typename T>
+    struct SetLargeCost
+      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
+      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
+    };
+    /// @}
+  public:
+    /// \brief Constructor.
+    ///
+    /// The constructor of the class.
+    ///
+    /// \param graph The digraph the algorithm runs on.
+    CostScaling(const GR& graph) :
+      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
+      _cost_map(_cost_vec), _pi_map(_pi),
+      INF(std::numeric_limits<Value>::has_infinity ?
+          std::numeric_limits<Value>::infinity() :
+          std::numeric_limits<Value>::max())
+    {
+    private:
+      const Map &_cost_map;
+    public:
+      ///\e
+      ResidualCostMap(const Map &cost_map) :
+        _cost_map(cost_map) {}
+      ///\e
+      inline typename Map::Value operator[](const ResArc &e) const {
+        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
+      }
+    }; //class ResidualCostMap
+    /// \brief Map adaptor class for handling reduced arc costs.
+    ///
+    /// Map adaptor class for handling reduced arc costs.
+    class ReducedCostMap : public MapBase<Arc, LCost>
+    {
+    private:
+      const Digraph &_gr;
+      const LargeCostMap &_cost_map;
+      const PotentialMap &_pot_map;
+    public:
+      ///\e
+      ReducedCostMap( const Digraph &gr,
+                      const LargeCostMap &cost_map,
+                      const PotentialMap &pot_map ) :
+        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
+      ///\e
+      inline LCost operator[](const Arc &e) const {
+        return _cost_map[e] + _pot_map[_gr.source(e)]
+                            - _pot_map[_gr.target(e)];
+      }
+    }; //class ReducedCostMap
+  private:
+    // The digraph the algorithm runs on
+    const Digraph &_graph;
+    // The original lower bound map
+    const LowerMap *_lower;
+    // The modified capacity map
+    CapacityArcMap _capacity;
+    // The original cost map
+    const CostMap &_orig_cost;
+    // The scaled cost map
+    LargeCostMap _cost;
+    // The modified supply map
+    SupplyNodeMap _supply;
+    bool _valid_supply;
+    // Arc map of the current flow
+    FlowMap *_flow;
+    bool _local_flow;
+    // Node map of the current potentials
+    PotentialMap *_potential;
+    bool _local_potential;
+    // The residual cost map
+    ResidualCostMap<LargeCostMap> _res_cost;
+    // The residual digraph
+    ResDigraph *_res_graph;
+    // The reduced cost map
+    ReducedCostMap *_red_cost;
+    // The excess map
+    SupplyNodeMap _excess;
+    // The epsilon parameter used for cost scaling
+    LCost _epsilon;
+    // The scaling factor
+    int _alpha;
+  public:
+    /// \brief General constructor (with lower bounds).
+    ///
+    /// General constructor (with lower bounds).
+    ///
+    /// \param digraph The digraph the algorithm runs on.
+    /// \param lower The lower bounds of the arcs.
+    /// \param capacity The capacities (upper bounds) of the arcs.
+    /// \param cost The cost (length) values of the arcs.
+    /// \param supply The supply values of the nodes (signed).
+    CostScaling( const Digraph &digraph,
+                 const LowerMap &lower,
+                 const CapacityMap &capacity,
+                 const CostMap &cost,
+                 const SupplyMap &supply ) :
+      _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
+      _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
+      _potential(NULL), _local_potential(false), _res_cost(_cost),
+      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
+    {
+      // Check the sum of supply values
+      Supply sum = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
+      _valid_supply = sum == 0;
+      // Check the value types
+      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+        "The flow type of CostScaling must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of CostScaling must be signed");
+      // Resize vectors
+      _node_num = countNodes(_graph);
+      _arc_num = countArcs(_graph);
+      _res_node_num = _node_num + 1;
+      _res_arc_num = 2 * (_arc_num + _node_num);
+      _root = _node_num;
+      _first_out.resize(_res_node_num + 1);
+      _forward.resize(_res_arc_num);
+      _source.resize(_res_arc_num);
+      _target.resize(_res_arc_num);
+      _reverse.resize(_res_arc_num);
+      _lower.resize(_res_arc_num);
+      _upper.resize(_res_arc_num);
+      _scost.resize(_res_arc_num);
+      _supply.resize(_res_node_num);
+      for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
+      for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
+      // Remove non-zero lower bounds
+      for (ArcIt e(_graph); e != INVALID; ++e) {
+        if (lower[e] != 0) {
+          _capacity[e] -= lower[e];
+          _supply[_graph.source(e)] -= lower[e];
+          _supply[_graph.target(e)] += lower[e];
+        }
+      }
+    }
+/*
+    /// \brief General constructor (without lower bounds).
+    ///
+    /// General constructor (without lower bounds).
+    ///
+    /// \param digraph The digraph the algorithm runs on.
+    /// \param capacity The capacities (upper bounds) of the arcs.
+    /// \param cost The cost (length) values of the arcs.
+    /// \param supply The supply values of the nodes (signed).
+    CostScaling( const Digraph &digraph,
+                 const CapacityMap &capacity,
+                 const CostMap &cost,
+                 const SupplyMap &supply ) :
+      _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
+      _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
+      _potential(NULL), _local_potential(false), _res_cost(_cost),
+      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
+    {
+      // Check the sum of supply values
+      Supply sum = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
+      _valid_supply = sum == 0;
+    }
+    /// \brief Simple constructor (with lower bounds).
+    ///
+    /// Simple constructor (with lower bounds).
+    ///
+    /// \param digraph The digraph the algorithm runs on.
+    /// \param lower The lower bounds of the arcs.
+    /// \param capacity The capacities (upper bounds) of the arcs.
+    /// \param cost The cost (length) values of the arcs.
+      _res_cap.resize(_res_arc_num);
+      _cost.resize(_res_arc_num);
+      _pi.resize(_res_node_num);
+      _excess.resize(_res_node_num);
+      _next_out.resize(_res_node_num);
+      _arc_vec.reserve(_res_arc_num);
+      _cost_vec.reserve(_res_arc_num);
+      // Copy the graph
+      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _node_id[n] = i;
+      }
+      i = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _first_out[i] = j;
+        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idf[a] = j;
+          _forward[j] = true;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idb[a] = j;
+          _forward[j] = false;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        _forward[j] = false;
+        _source[j] = i;
+        _target[j] = _root;
+        _reverse[j] = k;
+        _forward[k] = true;
+        _source[k] = _root;
+        _target[k] = i;
+        _reverse[k] = j;
+        ++j; ++k;
+      }
+      _first_out[i] = j;
+      _first_out[_res_node_num] = k;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int fi = _arc_idf[a];
+        int bi = _arc_idb[a];
+        _reverse[fi] = bi;
+        _reverse[bi] = fi;
+      }
+      // Reset parameters
+      reset();
+    }
+    /// \name Parameters
+    /// The parameters of the algorithm can be specified using these
+    /// functions.
+    /// @{
+    /// \brief Set the lower bounds on the arcs.
+    ///
+    /// This function sets the lower bounds on the arcs.
+    /// If it is not used before calling \ref run(), the lower bounds
+    /// will be set to zero on all arcs.
+    ///
+    /// \param map An arc map storing the lower bounds.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template <typename LowerMap>
+    CostScaling& lowerMap(const LowerMap& map) {
+      _have_lower = true;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _lower[_arc_idf[a]] = map[a];
+        _lower[_arc_idb[a]] = map[a];
+      }
+      return *this;
+    }
+    /// \brief Set the upper bounds (capacities) on the arcs.
+    ///
+    /// This function sets the upper bounds (capacities) on the arcs.
+    /// If it is not used before calling \ref run(), the upper bounds
+    /// will be set to \ref INF on all arcs (i.e. the flow value will be
+    /// unbounded from above on each arc).
+    ///
+    /// \param map An arc map storing the upper bounds.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename UpperMap>
+    CostScaling& upperMap(const UpperMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _upper[_arc_idf[a]] = map[a];
+      }
+      return *this;
+    }
+    /// \brief Set the costs of the arcs.
+    ///
+    /// This function sets the costs of the arcs.
+    /// If it is not used before calling \ref run(), the costs
+    /// will be set to \c 1 on all arcs.
+    ///
+    /// \param map An arc map storing the costs.
+    /// Its \c Value type must be convertible to the \c Cost type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename CostMap>
+    CostScaling& costMap(const CostMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _scost[_arc_idf[a]] =  map[a];
+        _scost[_arc_idb[a]] = -map[a];
+      }
+      return *this;
+    }
+    /// \brief Set the supply values of the nodes.
+    ///
+    /// This function sets the supply values of the nodes.
+    /// If neither this function nor \ref stSupply() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// \param map A node map storing the supply values.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename SupplyMap>
+    CostScaling& supplyMap(const SupplyMap& map) {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _supply[_node_id[n]] = map[n];
+      }
+      return *this;
+    }
+    /// \brief Set single source and target nodes and a supply value.
+    ///
+    /// This function sets a single source node and a single target node
+    /// and the required flow value.
+    /// If neither this function nor \ref supplyMap() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// Using this function has the same effect as using \ref supplyMap()
+    /// with such a map in which \c k is assigned to \c s, \c -k is
+    /// assigned to \c t and all other nodes have zero supply value.
+    ///
     /// \param s The source node.
     /// \param t The target node.
+    /// \param flow_value The required amount of flow from node \c s
+    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
+    CostScaling( const Digraph &digraph,
+                 const LowerMap &lower,
+                 const CapacityMap &capacity,
+                 const CostMap &cost,
+                 Node s, Node t,
+                 Supply flow_value ) :
+      _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
+      _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
+      _potential(NULL), _local_potential(false), _res_cost(_cost),
+      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
+    {
+      // Remove non-zero lower bounds
+      _supply[s] =  flow_value;
+      _supply[t] = -flow_value;
+      for (ArcIt e(_graph); e != INVALID; ++e) {
+        if (lower[e] != 0) {
+          _capacity[e] -= lower[e];
+          _supply[_graph.source(e)] -= lower[e];
+          _supply[_graph.target(e)] += lower[e];
+        }
+      }
+      _valid_supply = true;
+    }
+    /// \brief Simple constructor (without lower bounds).
+    ///
+    /// Simple constructor (without lower bounds).
+    ///
+    /// \param digraph The digraph the algorithm runs on.
+    /// \param capacity The capacities (upper bounds) of the arcs.
+    /// \param cost The cost (length) values of the arcs.
+    /// \param s The source node.
+    /// \param t The target node.
+    /// \param flow_value The required amount of flow from node \c s
+    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
+    CostScaling( const Digraph &digraph,
+                 const CapacityMap &capacity,
+                 const CostMap &cost,
+                 Node s, Node t,
+                 Supply flow_value ) :
+      _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
+      _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
+      _potential(NULL), _local_potential(false), _res_cost(_cost),
+      _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
+    {
+      _supply[s] =  flow_value;
+      _supply[t] = -flow_value;
+      _valid_supply = true;
+    }
+*/
+    /// Destructor.
+    ~CostScaling() {
+      if (_local_flow) delete _flow;
+      if (_local_potential) delete _potential;
+      delete _res_graph;
+      delete _red_cost;
+    }
+    /// \brief Set the flow map.
+    ///
+    /// Set the flow map.
+    ///
+    /// \return \c (*this)
+    CostScaling& flowMap(FlowMap &map) {
+      if (_local_flow) {
+        delete _flow;
+        _local_flow = false;
+      }
+      _flow = &map;
+    /// \param k The required amount of flow from node \c s to node \c t
+    /// (i.e. the supply of \c s and the demand of \c t).
+    ///
+    /// \return <tt>(*this)</tt>
+    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
+      for (int i = 0; i != _res_node_num; ++i) {
+        _supply[i] = 0;
+      }
+      _supply[_node_id[s]] =  k;
+      _supply[_node_id[t]] = -k;
       return *this;
+    }
+    /// \brief Set the potential map.
+    ///
+    /// Set the potential map.
+    ///
+    /// \return \c (*this)
+    CostScaling& potentialMap(PotentialMap &map) {
+      if (_local_potential) {
+        delete _potential;
+        _local_potential = false;
+      }
+      _potential = &map;
+      return *this;
+    }
+    /// @}
     /// \name Execution control
+    /// The algorithm can be executed using \ref run().
     /// @{
 …
     /// \brief Run the algorithm.
     ///
+    /// Run the algorithm.
+    /// This function runs the algorithm.
+    /// The paramters can be specified using functions \ref lowerMap(),
+    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
+    /// For example,
+    /// \code
+    ///   CostScaling<ListDigraph> cs(graph);
+    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// This function can be called more than once. All the parameters
+    /// that have been given are kept for the next call, unless
+    /// \ref reset() is called, thus only the modified parameters
+    /// have to be set again. See \ref reset() for examples.
+    /// However the underlying digraph must not be modified after this
+    /// class have been constructed, since it copies the digraph.
     ///
     /// \param partial_augment By default the algorithm performs
 …
     /// relabel operations instead.
     ///
+    /// \return \c true if a feasible flow can be found.
+    bool run(bool partial_augment = true) {
+    /// \return \c INFEASIBLE if no feasible flow exists,
+    /// \n \c OPTIMAL if the problem has optimal solution
+    /// (i.e. it is feasible and bounded), and the algorithm has found
+    /// optimal flow and node potentials (primal and dual solutions),
+    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
+    /// and infinite upper bound. It means that the objective function
+    /// is unbounded on that arc, however note that it could actually be
+    /// bounded over the feasible flows, but this algroithm cannot handle
+    /// these cases.
+    ///
+    /// \see ProblemType
+    ProblemType run(bool partial_augment = true) {
+      ProblemType pt = init();
+      if (pt != OPTIMAL) return pt;
+      start(partial_augment);
+      return OPTIMAL;
+    }
+    /// \brief Reset all the parameters that have been given before.
+    ///
+    /// This function resets all the paramaters that have been given
+    /// before using functions \ref lowerMap(), \ref upperMap(),
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
+    ///
+    /// It is useful for multiple run() calls. If this function is not
+    /// used, all the parameters given before are kept for the next
+    /// \ref run() call.
+    /// However the underlying digraph must not be modified after this
+    /// class have been constructed, since it copies and extends the graph.
+    ///
+    /// For example,
+    /// \code
+    ///   CostScaling<ListDigraph> cs(graph);
+    ///
+    ///   // First run
+    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    ///
+    ///   // Run again with modified cost map (reset() is not called,
+    ///   // so only the cost map have to be set again)
+    ///   cost[e] += 100;
+    ///   cs.costMap(cost).run();
+    ///
+    ///   // Run again from scratch using reset()
+    ///   // (the lower bounds will be set to zero on all arcs)
+    ///   cs.reset();
+    ///   cs.upperMap(capacity).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// \return <tt>(*this)</tt>
+    CostScaling& reset() {
+      for (int i = 0; i != _res_node_num; ++i) {
+        _supply[i] = 0;
+      }
+      int limit = _first_out[_root];
+      for (int j = 0; j != limit; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _scost[j] = _forward[j] ? 1 : -1;
+      }
+      for (int j = limit; j != _res_arc_num; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _scost[j] = 0;
+        _scost[_reverse[j]] = 0;
+      }
+      _have_lower = false;
+      return *this;
+    }
+    /// @}
+    /// \name Query Functions
+    /// The results of the algorithm can be obtained using these
+    /// functions.\n
+    /// The \ref run() function must be called before using them.
+    /// @{
+    /// \brief Return the total cost of the found flow.
+    ///
+    /// This function returns the total cost of the found flow.
+    /// Its complexity is O(e).
+    ///
+    /// \note The return type of the function can be specified as a
+    /// template parameter. For example,
+    /// \code
+    ///   cs.totalCost<double>();
+    /// \endcode
+    /// It is useful if the total cost cannot be stored in the \c Cost
+    /// type of the algorithm, which is the default return type of the
+    /// function.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    template <typename Number>
+    Number totalCost() const {
+      Number c = 0;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int i = _arc_idb[a];
+        c += static_cast<Number>(_res_cap[i]) *
+             (-static_cast<Number>(_scost[i]));
+      }
+      return c;
+    }
+#ifndef DOXYGEN
+    Cost totalCost() const {
+      return totalCost<Cost>();
+    }
+#endif
+    /// \brief Return the flow on the given arc.
+    ///
+    /// This function returns the flow on the given arc.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Value flow(const Arc& a) const {
+      return _res_cap[_arc_idb[a]];
+    }
+    /// \brief Return the flow map (the primal solution).
+    ///
+    /// This function copies the flow value on each arc into the given
+    /// map. The \c Value type of the algorithm must be convertible to
+    /// the \c Value type of the map.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    template <typename FlowMap>
+    void flowMap(FlowMap &map) const {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        map.set(a, _res_cap[_arc_idb[a]]);
+      }
+    }
+    /// \brief Return the potential (dual value) of the given node.
+    ///
+    /// This function returns the potential (dual value) of the
+    /// given node.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost potential(const Node& n) const {
+      return static_cast<Cost>(_pi[_node_id[n]]);
+    }
+    /// \brief Return the potential map (the dual solution).
+    ///
+    /// This function copies the potential (dual value) of each node
+    /// into the given map.
+    /// The \c Cost type of the algorithm must be convertible to the
+    /// \c Value type of the map.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    template <typename PotentialMap>
+    void potentialMap(PotentialMap &map) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
+      }
+    }
+    /// @}
+  private:
+    // Initialize the algorithm
+    ProblemType init() {
+      if (_res_node_num == 0) return INFEASIBLE;
+      // Scaling factor
+      _alpha = 8;
+      // Check the sum of supply values
+      _sum_supply = 0;
+      for (int i = 0; i != _root; ++i) {
+        _sum_supply += _supply[i];
+      }
+      if (_sum_supply > 0) return INFEASIBLE;
+      // Initialize vectors
+      for (int i = 0; i != _res_node_num; ++i) {
+        _pi[i] = 0;
+        _excess[i] = _supply[i];
+      }
+      // Remove infinite upper bounds and check negative arcs
+      const Value MAX = std::numeric_limits<Value>::max();
+      int last_out;
+      if (_have_lower) {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j]) {
+              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
+              if (c >= MAX) return UNBOUNDED;
+              _excess[i] -= c;
+              _excess[_target[j]] += c;
+            }
+          }
+        }
+      } else {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j] && _scost[j] < 0) {
+              Value c = _upper[j];
+              if (c >= MAX) return UNBOUNDED;
+              _excess[i] -= c;
+              _excess[_target[j]] += c;
+            }
+          }
+        }
+      }
+      Value ex, max_cap = 0;
+      for (int i = 0; i != _res_node_num; ++i) {
+        ex = _excess[i];
+        _excess[i] = 0;
+        if (ex < 0) max_cap -= ex;
+      }
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_upper[j] >= MAX) _upper[j] = max_cap;
+      }
+      // Initialize the large cost vector and the epsilon parameter
+      _epsilon = 0;
+      LargeCost lc;
+      for (int i = 0; i != _root; ++i) {
+        last_out = _first_out[i+1];
+        for (int j = _first_out[i]; j != last_out; ++j) {
+          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
+          _cost[j] = lc;
+          if (lc > _epsilon) _epsilon = lc;
+        }
+      }
+      _epsilon /= _alpha;
+      // Initialize maps for Circulation and remove non-zero lower bounds
+      ConstMap<Arc, Value> low(0);
+      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
+      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
+      ValueArcMap cap(_graph), flow(_graph);
+      ValueNodeMap sup(_graph);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        sup[n] = _supply[_node_id[n]];
+      }
+      if (_have_lower) {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          int j = _arc_idf[a];
+          Value c = _lower[j];
+          cap[a] = _upper[j] - c;
+          sup[_graph.source(a)] -= c;
+          sup[_graph.target(a)] += c;
+        }
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          cap[a] = _upper[_arc_idf[a]];
+        }
+      }
+      // Find a feasible flow using Circulation
+      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
+        circ(_graph, low, cap, sup);
+      if (!circ.flowMap(flow).run()) return INFEASIBLE;
+      // Set residual capacities and handle GEQ supply type
+      if (_sum_supply < 0) {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          Value fa = flow[a];
+          _res_cap[_arc_idf[a]] = cap[a] - fa;
+          _res_cap[_arc_idb[a]] = fa;
+          sup[_graph.source(a)] -= fa;
+          sup[_graph.target(a)] += fa;
+        }
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          _excess[_node_id[n]] = sup[n];
+        }
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int u = _target[a];
+          int ra = _reverse[a];
+          _res_cap[a] = -_sum_supply + 1;
+          _res_cap[ra] = -_excess[u];
+          _cost[a] = 0;
+          _cost[ra] = 0;
+          _excess[u] = 0;
+        }
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          Value fa = flow[a];
+          _res_cap[_arc_idf[a]] = cap[a] - fa;
+          _res_cap[_arc_idb[a]] = fa;
+        }
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int ra = _reverse[a];
+          _res_cap[a] = 1;
+          _res_cap[ra] = 0;
+          _cost[a] = 0;
+          _cost[ra] = 0;
+        }
+      }
+      return OPTIMAL;
+    }
+    // Execute the algorithm and transform the results
+    void start(bool partial_augment) {
+      // Execute the algorithm
       if (partial_augment) {
         return init() && startPartialAugment();
+        startPartialAugment();
       } else {
+        return init() && startPushRelabel();
+      }
+    }
+    /// @}
+    /// \name Query Functions
+    /// The result of the algorithm can be obtained using these
+    /// functions.\n
+    /// \ref lemon::CostScaling::run() "run()" must be called before
+    /// using them.
+    /// @{
+    /// \brief Return a const reference to the arc map storing the
+    /// found flow.
+    ///
+    /// Return a const reference to the arc map storing the found flow.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const FlowMap& flowMap() const {
+      return *_flow;
+    }
+    /// \brief Return a const reference to the node map storing the
+    /// found potentials (the dual solution).
+    ///
+    /// Return a const reference to the node map storing the found
+    /// potentials (the dual solution).
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const PotentialMap& potentialMap() const {
+      return *_potential;
+    }
+    /// \brief Return the flow on the given arc.
+    ///
+    /// Return the flow on the given arc.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Capacity flow(const Arc& arc) const {
+      return (*_flow)[arc];
+    }
+    /// \brief Return the potential of the given node.
+    ///
+    /// Return the potential of the given node.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost potential(const Node& node) const {
+      return (*_potential)[node];
+    }
+    /// \brief Return the total cost of the found flow.
+    ///
+    /// Return the total cost of the found flow. The complexity of the
+    /// function is \f$ O(e) \f$.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost totalCost() const {
+      Cost c = 0;
+      for (ArcIt e(_graph); e != INVALID; ++e)
+        c += (*_flow)[e] * _orig_cost[e];
+      return c;
+    }
+    /// @}
+  private:
+    /// Initialize the algorithm.
+    bool init() {
+      if (!_valid_supply) return false;
+      // The scaling factor
+      _alpha = 8;
+      // Initialize flow and potential maps
+      if (!_flow) {
+        _flow = new FlowMap(_graph);
+        _local_flow = true;
+      }
+      if (!_potential) {
+        _potential = new PotentialMap(_graph);
+        _local_potential = true;
+      }
+      _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
+      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
+      // Initialize the scaled cost map and the epsilon parameter
+      Cost max_cost = 0;
+      int node_num = countNodes(_graph);
+      for (ArcIt e(_graph); e != INVALID; ++e) {
+        _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
+        if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
+      }
+      _epsilon = max_cost * node_num;
+      // Find a feasible flow using Circulation
+      Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
+                   SupplyMap >
+        circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
+                     _supply );
+      return circulation.flowMap(*_flow).run();
+        startPushRelabel();
+      }
+      // Compute node potentials for the original costs
+      _arc_vec.clear();
+      _cost_vec.clear();
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_res_cap[j] > 0) {
+          _arc_vec.push_back(IntPair(_source[j], _target[j]));
+          _cost_vec.push_back(_scost[j]);
+        }
+      }
+      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+      typename BellmanFord<StaticDigraph, LargeCostArcMap>
+        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
+      bf.distMap(_pi_map);
+      bf.init(0);
+      bf.start();
+      // Handle non-zero lower bounds
+      if (_have_lower) {
+        int limit = _first_out[_root];
+        for (int j = 0; j != limit; ++j) {
+          if (!_forward[j]) _res_cap[j] += _lower[j];
+        }
+      }
+    }
     /// Execute the algorithm performing partial augmentation and
     /// relabel operations.
     bool startPartialAugment() {
+    /// relabel operations
+    void startPartialAugment() {
       // Paramters for heuristics
 //      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
 //      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
       // Maximum augment path length
       const int MAX_PATH_LENGTH = 4;
+      // Variables
+      typename Digraph::template NodeMap<Arc> pred_arc(_graph);
+      typename Digraph::template NodeMap<bool> forward(_graph);
+      typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
+      typename Digraph::template NodeMap<InArcIt> next_in(_graph);
+      typename Digraph::template NodeMap<bool> next_dir(_graph);
+      std::deque<Node> active_nodes;
+      std::vector<Node> path_nodes;
+//      int node_num = countNodes(_graph);
+      // Perform cost scaling phases
+      IntVector pred_arc(_res_node_num);
+      std::vector<int> path_nodes;
       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
 : _epsilon / _alpha )
+      {
-/*
         // "Early Termination" heuristic: use Bellman-Ford algorithm
         // to check if the current flow is optimal
         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
+          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
+          ShiftCostMap shift_cost(_res_cost, 1);
+          BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
+          _arc_vec.clear();
+          _cost_vec.clear();
+          for (int j = 0; j != _res_arc_num; ++j) {
+            if (_res_cap[j] > 0) {
+              _arc_vec.push_back(IntPair(_source[j], _target[j]));
+              _cost_vec.push_back(_cost[j] + 1);
+            }
+          }
+          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
           bf.init(0);
           bool done = false;
           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
+          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
           for (int i = 0; i < K && !done; ++i)
             done = bf.processNextWeakRound();
           if (done) break;
+        }
+*/
         // Saturate arcs not satisfying the optimality condition
+        Capacity delta;
+        for (ArcIt e(_graph); e != INVALID; ++e) {
+          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
+            delta = _capacity[e] - (*_flow)[e];
+            _excess[_graph.source(e)] -= delta;
+            _excess[_graph.target(e)] += delta;
+            (*_flow)[e] = _capacity[e];
+        for (int a = 0; a != _res_arc_num; ++a) {
+          if (_res_cap[a] > 0 &&
+              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+            Value delta = _res_cap[a];
+            _excess[_source[a]] -= delta;
+            _excess[_target[a]] += delta;
+            _res_cap[a] = 0;
+            _res_cap[_reverse[a]] += delta;
+          }
+          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+            _excess[_graph.target(e)] -= (*_flow)[e];
+            _excess[_graph.source(e)] += (*_flow)[e];
+            (*_flow)[e] = 0;
+        }
+        // Find active nodes (i.e. nodes with positive excess)
+        for (int u = 0; u != _res_node_num; ++u) {
+          if (_excess[u] > 0) _active_nodes.push_back(u);
+        }
+        // Initialize the next arcs
+        for (int u = 0; u != _res_node_num; ++u) {
+          _next_out[u] = _first_out[u];
+        }
+        // Perform partial augment and relabel operations
+        while (true) {
+          // Select an active node (FIFO selection)
+          while (_active_nodes.size() > 0 &&
+                 _excess[_active_nodes.front()] <= 0) {
+            _active_nodes.pop_front();
+          }
+        }
+        // Find active nodes (i.e. nodes with positive excess)
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          if (_excess[n] > 0) active_nodes.push_back(n);
+        }
+        // Initialize the next arc maps
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          next_out[n] = OutArcIt(_graph, n);
+          next_in[n] = InArcIt(_graph, n);
+          next_dir[n] = true;
+        }
+        // Perform partial augment and relabel operations
+        while (active_nodes.size() > 0) {
+          // Select an active node (FIFO selection)
+          if (_excess[active_nodes[0]] <= 0) {
+            active_nodes.pop_front();
+            continue;
+          }
+          Node start = active_nodes[0];
+          if (_active_nodes.size() == 0) break;
+          int start = _active_nodes.front();
           path_nodes.clear();
           path_nodes.push_back(start);
           // Find an augmenting path from the start node
+          Node u, tip = start;
+          LCost min_red_cost;
+          while ( _excess[tip] >= 0 &&
+                  int(path_nodes.size()) <= MAX_PATH_LENGTH )
+          {
+            if (next_dir[tip]) {
+              for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
+                if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
+                  u = _graph.target(e);
+                  pred_arc[u] = e;
+                  forward[u] = true;
+                  next_out[tip] = e;
+                  tip = u;
+                  path_nodes.push_back(tip);
+                  goto next_step;
+                }
+              }
+              next_dir[tip] = false;
+            }
+            for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
+              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+                u = _graph.source(e);
+                pred_arc[u] = e;
+                forward[u] = false;
+                next_in[tip] = e;
+          int tip = start;
+          while (_excess[tip] >= 0 &&
+                 int(path_nodes.size()) <= MAX_PATH_LENGTH) {
+            int u;
+            LargeCost min_red_cost, rc;
+            int last_out = _sum_supply < 0 ?
+              _first_out[tip+1] : _first_out[tip+1] - 1;
+            for (int a = _next_out[tip]; a != last_out; ++a) {
+              if (_res_cap[a] > 0 &&
+                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+                u = _target[a];
+                pred_arc[u] = a;
+                _next_out[tip] = a;
                 tip = u;
                 path_nodes.push_back(tip);
 …
             // Relabel tip node
+            min_red_cost = std::numeric_limits<LCost>::max() / 2;
+            for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
+              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
+                   (*_red_cost)[oe] < min_red_cost )
+                min_red_cost = (*_red_cost)[oe];
+            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
+            for (int a = _first_out[tip]; a != last_out; ++a) {
+              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
+              if (_res_cap[a] > 0 && rc < min_red_cost) {
+                min_red_cost = rc;
+              }
+            }
+            for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
+              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
+                min_red_cost = -(*_red_cost)[ie];
+            }
+            (*_potential)[tip] -= min_red_cost + _epsilon;
+            // Reset the next arc maps
+            next_out[tip] = OutArcIt(_graph, tip);
+            next_in[tip] = InArcIt(_graph, tip);
+            next_dir[tip] = true;
+            _pi[tip] -= min_red_cost + _epsilon;
+            // Reset the next arc of tip
+            _next_out[tip] = _first_out[tip];
             // Step back
             if (tip != start) {
               path_nodes.pop_back();
               tip = path_nodes[path_nodes.size()-1];
+              tip = path_nodes.back();
+            }
+          next_step:
+            continue;
+          next_step: ;
+          }
           // Augment along the found path (as much flow as possible)
+          Capacity delta;
+          Value delta;
+          int u, v = path_nodes.front(), pa;
           for (int i = 1; i < int(path_nodes.size()); ++i) {
+            u = path_nodes[i];
+            delta = forward[u] ?
+              _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
+              (*_flow)[pred_arc[u]];
+            delta = std::min(delta, _excess[path_nodes[i-1]]);
+            (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
+            _excess[path_nodes[i-1]] -= delta;
+            _excess[u] += delta;
+            if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
+            u = v;
+            v = path_nodes[i];
+            pa = pred_arc[v];
+            delta = std::min(_res_cap[pa], _excess[u]);
+            _res_cap[pa] -= delta;
+            _res_cap[_reverse[pa]] += delta;
+            _excess[u] -= delta;
+            _excess[v] += delta;
+            if (_excess[v] > 0 && _excess[v] <= delta)
+              _active_nodes.push_back(v);
+          }
+        }
+      }
+      // Compute node potentials for the original costs
+      ResidualCostMap<CostMap> res_cost(_orig_cost);
+      BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
+        bf(*_res_graph, res_cost);
+      bf.init(0); bf.start();
+      for (NodeIt n(_graph); n != INVALID; ++n)
+        (*_potential)[n] = bf.dist(n);
+      // Handle non-zero lower bounds
+      if (_lower) {
+        for (ArcIt e(_graph); e != INVALID; ++e)
+          (*_flow)[e] += (*_lower)[e];
+      }
+      return true;
+    }
+    /// Execute the algorithm performing push and relabel operations.
+    bool startPushRelabel() {
+    }
+    /// Execute the algorithm performing push and relabel operations
+    void startPushRelabel() {
       // Paramters for heuristics
+//      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+//      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+      typename Digraph::template NodeMap<bool> hyper(_graph, false);
+      typename Digraph::template NodeMap<Arc> pred_arc(_graph);
+      typename Digraph::template NodeMap<bool> forward(_graph);
+      typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
+      typename Digraph::template NodeMap<InArcIt> next_in(_graph);
+      typename Digraph::template NodeMap<bool> next_dir(_graph);
+      std::deque<Node> active_nodes;
+//      int node_num = countNodes(_graph);
+      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+      // Perform cost scaling phases
+      BoolVector hyper(_res_node_num, false);
       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
 : _epsilon / _alpha )
+      {
-/*
         // "Early Termination" heuristic: use Bellman-Ford algorithm
         // to check if the current flow is optimal
         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
+          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
+          ShiftCostMap shift_cost(_res_cost, 1);
+          BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
+          _arc_vec.clear();
+          _cost_vec.clear();
+          for (int j = 0; j != _res_arc_num; ++j) {
+            if (_res_cap[j] > 0) {
+              _arc_vec.push_back(IntPair(_source[j], _target[j]));
+              _cost_vec.push_back(_cost[j] + 1);
+            }
+          }
+          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
           bf.init(0);
           bool done = false;
           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
+          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
           for (int i = 0; i < K && !done; ++i)
             done = bf.processNextWeakRound();
           if (done) break;
+        }
-*/
         // Saturate arcs not satisfying the optimality condition
+        Capacity delta;
+        for (ArcIt e(_graph); e != INVALID; ++e) {
+          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
+            delta = _capacity[e] - (*_flow)[e];
+            _excess[_graph.source(e)] -= delta;
+            _excess[_graph.target(e)] += delta;
+            (*_flow)[e] = _capacity[e];
+        for (int a = 0; a != _res_arc_num; ++a) {
+          if (_res_cap[a] > 0 &&
+              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+            Value delta = _res_cap[a];
+            _excess[_source[a]] -= delta;
+            _excess[_target[a]] += delta;
+            _res_cap[a] = 0;
+            _res_cap[_reverse[a]] += delta;
+          }
-          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
-            _excess[_graph.target(e)] -= (*_flow)[e];
-            _excess[_graph.source(e)] += (*_flow)[e];
-            (*_flow)[e] = 0;
+          }
+        }
         // Find active nodes (i.e. nodes with positive excess)
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          if (_excess[n] > 0) active_nodes.push_back(n);
+        }
+        // Initialize the next arc maps
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          next_out[n] = OutArcIt(_graph, n);
+          next_in[n] = InArcIt(_graph, n);
+          next_dir[n] = true;
+        for (int u = 0; u != _res_node_num; ++u) {
+          if (_excess[u] > 0) _active_nodes.push_back(u);
+        }
+        // Initialize the next arcs
+        for (int u = 0; u != _res_node_num; ++u) {
+          _next_out[u] = _first_out[u];
+        }
         // Perform push and relabel operations
+        while (active_nodes.size() > 0) {
+        while (_active_nodes.size() > 0) {
+          LargeCost min_red_cost, rc;
+          Value delta;
+          int n, t, a, last_out = _res_arc_num;
           // Select an active node (FIFO selection)
+          Node n = active_nodes[0], t;
+          bool relabel_enabled = true;
+        next_node:
+          n = _active_nodes.front();
+          last_out = _sum_supply < 0 ?
+            _first_out[n+1] : _first_out[n+1] - 1;
           // Perform push operations if there are admissible arcs
           if (_excess[n] > 0 && next_dir[n]) {
             OutArcIt e = next_out[n];
             for ( ; e != INVALID; ++e) {
               if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
                 delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
                 t = _graph.target(e);
+          if (_excess[n] > 0) {
+            for (a = _next_out[n]; a != last_out; ++a) {
+              if (_res_cap[a] > 0 &&
+                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+                delta = std::min(_res_cap[a], _excess[n]);
+                t = _target[a];
                 // Push-look-ahead heuristic
                 Capacity ahead = -_excess[t];
                 for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
                   if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
                     ahead += _capacity[oe] - (*_flow)[oe];
+                }
                 for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
                   if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
                     ahead += (*_flow)[ie];
+                Value ahead = -_excess[t];
+                int last_out_t = _sum_supply < 0 ?
+                  _first_out[t+1] : _first_out[t+1] - 1;
+                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
+                  if (_res_cap[ta] > 0 &&
+                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
+                    ahead += _res_cap[ta];
+                  if (ahead >= delta) break;
+                }
                 if (ahead < 0) ahead = 0;
 …
                 // Push flow along the arc
                 if (ahead < delta) {
+                  (*_flow)[e] += ahead;
+                  _res_cap[a] -= ahead;
+                  _res_cap[_reverse[a]] += ahead;
                   _excess[n] -= ahead;
                   _excess[t] += ahead;
                   active_nodes.push_front(t);
+                  _active_nodes.push_front(t);
                   hyper[t] = true;
                   relabel_enabled = false;
                   break;
+                  _next_out[n] = a;
+                  goto next_node;
                 } else {
+                  (*_flow)[e] += delta;
+                  _res_cap[a] -= delta;
+                  _res_cap[_reverse[a]] += delta;
                   _excess[n] -= delta;
                   _excess[t] += delta;
                   if (_excess[t] > 0 && _excess[t] <= delta)
                     active_nodes.push_back(t);
+                    _active_nodes.push_back(t);
+                }
+                if (_excess[n] == 0) break;
+                if (_excess[n] == 0) {
+                  _next_out[n] = a;
+                  goto remove_nodes;
+                }
+              }
+            }
+            if (e != INVALID) {
+              next_out[n] = e;
+            } else {
+              next_dir[n] = false;
+            }
+            _next_out[n] = a;
+          }
+          if (_excess[n] > 0 && !next_dir[n]) {
+            InArcIt e = next_in[n];
+            for ( ; e != INVALID; ++e) {
+              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+                delta = std::min((*_flow)[e], _excess[n]);
+                t = _graph.source(e);
+                // Push-look-ahead heuristic
+                Capacity ahead = -_excess[t];
+                for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
+                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
+                    ahead += _capacity[oe] - (*_flow)[oe];
+                }
+                for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
+                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
+                    ahead += (*_flow)[ie];
+                }
+                if (ahead < 0) ahead = 0;
+                // Push flow along the arc
+                if (ahead < delta) {
+                  (*_flow)[e] -= ahead;
+                  _excess[n] -= ahead;
+                  _excess[t] += ahead;
+                  active_nodes.push_front(t);
+                  hyper[t] = true;
+                  relabel_enabled = false;
+                  break;
+                } else {
+                  (*_flow)[e] -= delta;
+                  _excess[n] -= delta;
+                  _excess[t] += delta;
+                  if (_excess[t] > 0 && _excess[t] <= delta)
+                    active_nodes.push_back(t);
+                }
+                if (_excess[n] == 0) break;
+          // Relabel the node if it is still active (or hyper)
+          if (_excess[n] > 0 || hyper[n]) {
+            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
+            for (int a = _first_out[n]; a != last_out; ++a) {
+              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
+              if (_res_cap[a] > 0 && rc < min_red_cost) {
+                min_red_cost = rc;
+              }
+            }
+            next_in[n] = e;
+            _pi[n] -= min_red_cost + _epsilon;
+            hyper[n] = false;
+            // Reset the next arc
+            _next_out[n] = _first_out[n];
+          }
+          // Relabel the node if it is still active (or hyper)
+          if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
+            LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
+            for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
+              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
+                   (*_red_cost)[oe] < min_red_cost )
+                min_red_cost = (*_red_cost)[oe];
+            }
+            for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
+              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
+                min_red_cost = -(*_red_cost)[ie];
+            }
+            (*_potential)[n] -= min_red_cost + _epsilon;
+            hyper[n] = false;
+            // Reset the next arc maps
+            next_out[n] = OutArcIt(_graph, n);
+            next_in[n] = InArcIt(_graph, n);
+            next_dir[n] = true;
+          // Remove nodes that are not active nor hyper
+        remove_nodes:
+          while ( _active_nodes.size() > 0 &&
+                  _excess[_active_nodes.front()] <= 0 &&
+                  !hyper[_active_nodes.front()] ) {
+            _active_nodes.pop_front();
+          }
+          // Remove nodes that are not active nor hyper
+          while ( active_nodes.size() > 0 &&
+                  _excess[active_nodes[0]] <= 0 &&
+                  !hyper[active_nodes[0]] ) {
+            active_nodes.pop_front();
+          }
+        }
+      }
+      // Compute node potentials for the original costs
+      ResidualCostMap<CostMap> res_cost(_orig_cost);
+      BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
+        bf(*_res_graph, res_cost);
+      bf.init(0); bf.start();
+      for (NodeIt n(_graph); n != INVALID; ++n)
+        (*_potential)[n] = bf.dist(n);
+      // Handle non-zero lower bounds
+      if (_lower) {
+        for (ArcIt e(_graph); e != INVALID; ++e)
+          (*_flow)[e] += (*_lower)[e];
+      }
+      return true;
+        }
+      }
+    }

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 809:22bb98ca0101 in lemon-main

Legend:

lemon/cost_scaling.h

Download in other formats: