diff -r 0643a9c2c3ae -r aef153f430e1 lemon/cycle_canceling.h
--- a/lemon/cycle_canceling.h	Fri Nov 13 00:09:35 2009 +0100
+++ b/lemon/cycle_canceling.h	Fri Nov 13 00:10:33 2009 +0100
@@ -19,441 +19,817 @@
 #ifndef LEMON_CYCLE_CANCELING_H
 #define LEMON_CYCLE_CANCELING_H
 
-/// \ingroup min_cost_flow
-///
+/// \ingroup min_cost_flow_algs
 /// \file
-/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
+/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
 
 #include <vector>
+#include <limits>
+
+#include <lemon/core.h>
+#include <lemon/maps.h>
+#include <lemon/path.h>
+#include <lemon/math.h>
+#include <lemon/static_graph.h>
 #include <lemon/adaptors.h>
-#include <lemon/path.h>
-
 #include <lemon/circulation.h>
 #include <lemon/bellman_ford.h>
 #include <lemon/howard.h>
 
 namespace lemon {
 
-  /// \addtogroup min_cost_flow
+  /// \addtogroup min_cost_flow_algs
   /// @{
 
-  /// \brief Implementation of a cycle-canceling algorithm for
-  /// finding a minimum cost flow.
+  /// \brief Implementation of cycle-canceling algorithms for
+  /// finding a \ref min_cost_flow "minimum cost flow".
   ///
-  /// \ref CycleCanceling implements a cycle-canceling algorithm for
-  /// finding a minimum cost flow.
+  /// \ref CycleCanceling implements three different cycle-canceling
+  /// algorithms for finding a \ref min_cost_flow "minimum cost flow".
+  /// The most efficent one (both theoretically and practically)
+  /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
+  /// thus it is the default method.
+  /// It is strongly polynomial, but in practice, it is typically much
+  /// slower than the scaling algorithms and NetworkSimplex.
   ///
-  /// \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 number type used for flow amounts, capacity bounds
+  /// and supply values in the algorithm. By default, it is \c int.
+  /// \tparam C The number type used for costs and potentials in the
+  /// algorithm. By default, it is the same as \c V.
   ///
-  /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
-  /// used for negative cycle detection with limited iteration number.
-  /// However \ref CycleCanceling also provides the "Minimum Mean
-  /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
-  /// but rather slower in practice.
-  /// To use this version of the algorithm, call \ref run() with \c true
-  /// parameter.
+  /// \warning Both number 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.
   ///
-  /// \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> >
+  /// \note For more information about the three available methods,
+  /// see \ref Method.
+#ifdef DOXYGEN
+  template <typename GR, typename V, typename C>
+#else
+  template <typename GR, typename V = int, typename C = V>
+#endif
   class CycleCanceling
   {
-    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+  public:
 
-    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::Node ResNode;
-    typedef typename ResDigraph::NodeIt ResNodeIt;
-    typedef typename ResDigraph::Arc ResArc;
-    typedef typename ResDigraph::ArcIt ResArcIt;
+    /// 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;
 
   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<Cost> PotentialMap;
+    /// \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
+    };
+
+    /// \brief Constants for selecting the used method.
+    ///
+    /// Enum type containing constants for selecting the used method
+    /// for the \ref run() function.
+    ///
+    /// \ref CycleCanceling provides three different cycle-canceling
+    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
+    /// is used, which proved to be the most efficient and the most robust
+    /// on various test inputs.
+    /// However, the other methods can be selected using the \ref run()
+    /// function with the proper parameter.
+    enum Method {
+      /// A simple cycle-canceling method, which uses the
+      /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
+      /// number for detecting negative cycles in the residual network.
+      SIMPLE_CYCLE_CANCELING,
+      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
+      /// well-known strongly polynomial method. It improves along a
+      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
+      /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
+      MINIMUM_MEAN_CYCLE_CANCELING,
+      /// The "Cancel And Tighten" algorithm, which can be viewed as an
+      /// improved version of the previous method.
+      /// It is faster both in theory and in practice, its running time
+      /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
+      CANCEL_AND_TIGHTEN
+    };
 
   private:
 
-    /// \brief Map adaptor class for handling residual arc costs.
-    ///
-    /// Map adaptor class for handling residual arc costs.
-    class ResidualCostMap : public MapBase<ResArc, Cost>
-    {
-    private:
+    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
+    
+    typedef std::vector<int> IntVector;
+    typedef std::vector<char> CharVector;
+    typedef std::vector<double> DoubleVector;
+    typedef std::vector<Value> ValueVector;
+    typedef std::vector<Cost> CostVector;
 
-      const CostMap &_cost_map;
-
+  private:
+  
+    template <typename KT, typename VT>
+    class VectorMap {
     public:
-
-      ///\e
-      ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
-
-      ///\e
-      Cost operator[](const ResArc &e) const {
-        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
+      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)];
       }
 
-    }; //class ResidualCostMap
+      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, Cost> CostNodeMap;
+    typedef VectorMap<StaticDigraph::Arc, Cost> CostArcMap;
 
   private:
 
-    // The maximum number of iterations for the first execution of the
-    // Bellman-Ford algorithm. It should be at least 2.
-    static const int BF_FIRST_LIMIT  = 2;
-    // The iteration limit for the Bellman-Ford algorithm is multiplied
-    // by BF_LIMIT_FACTOR/100 in every round.
-    static const int BF_LIMIT_FACTOR = 150;
 
-  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;
 
-    // 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 &_cost;
-    // The modified supply map
-    SupplyNodeMap _supply;
-    bool _valid_supply;
+    // Parameters of the problem
+    bool _have_lower;
+    Value _sum_supply;
 
-    // Arc map of the current flow
-    FlowMap *_flow;
-    bool _local_flow;
-    // Node map of the current potentials
-    PotentialMap *_potential;
-    bool _local_potential;
+    // Data structures for storing the digraph
+    IntNodeMap _node_id;
+    IntArcMap _arc_idf;
+    IntArcMap _arc_idb;
+    IntVector _first_out;
+    CharVector _forward;
+    IntVector _source;
+    IntVector _target;
+    IntVector _reverse;
 
-    // The residual digraph
-    ResDigraph *_res_graph;
-    // The residual cost map
-    ResidualCostMap _res_cost;
+    // Node and arc data
+    ValueVector _lower;
+    ValueVector _upper;
+    CostVector _cost;
+    ValueVector _supply;
+
+    ValueVector _res_cap;
+    CostVector _pi;
+
+    // Data for a StaticDigraph structure
+    typedef std::pair<int, int> IntPair;
+    StaticDigraph _sgr;
+    std::vector<IntPair> _arc_vec;
+    std::vector<Cost> _cost_vec;
+    IntVector _id_vec;
+    CostArcMap _cost_map;
+    CostNodeMap _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:
 
-    /// \brief General constructor (with lower bounds).
+    /// \brief Constructor.
     ///
-    /// General constructor (with lower bounds).
+    /// The constructor of the class.
     ///
-    /// \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).
-    CycleCanceling( const Digraph &digraph,
-                    const LowerMap &lower,
-                    const CapacityMap &capacity,
-                    const CostMap &cost,
-                    const SupplyMap &supply ) :
-      _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
-      _supply(digraph), _flow(NULL), _local_flow(false),
-      _potential(NULL), _local_potential(false),
-      _res_graph(NULL), _res_cost(_cost)
+    /// \param graph The digraph the algorithm runs on.
+    CycleCanceling(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())
     {
-      // Check the sum of supply values
-      Supply sum = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        _supply[n] = supply[n];
-        sum += _supply[n];
+      // Check the number types
+      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+        "The flow type of CycleCanceling must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of CycleCanceling 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);
+      _cost.resize(_res_arc_num);
+      _supply.resize(_res_node_num);
+      
+      _res_cap.resize(_res_arc_num);
+      _pi.resize(_res_node_num);
+
+      _arc_vec.reserve(_res_arc_num);
+      _cost_vec.reserve(_res_arc_num);
+      _id_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;
       }
-      _valid_supply = sum == 0;
-
-      // Remove non-zero lower bounds
-      for (ArcIt e(_graph); e != INVALID; ++e) {
-        _capacity[e] = capacity[e];
-        if (lower[e] != 0) {
-          _capacity[e] -= lower[e];
-          _supply[_graph.source(e)] -= lower[e];
-          _supply[_graph.target(e)] += lower[e];
+      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;
       }
-    }
-/*
-    /// \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).
-    CycleCanceling( const Digraph &digraph,
-                    const CapacityMap &capacity,
-                    const CostMap &cost,
-                    const SupplyMap &supply ) :
-      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
-      _supply(supply), _flow(NULL), _local_flow(false),
-      _potential(NULL), _local_potential(false), _res_graph(NULL),
-      _res_cost(_cost)
-    {
-      // Check the sum of supply values
-      Supply sum = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
-      _valid_supply = sum == 0;
+      _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();
     }
 
-    /// \brief Simple constructor (with lower bounds).
+    /// \name Parameters
+    /// The parameters of the algorithm can be specified using these
+    /// functions.
+
+    /// @{
+
+    /// \brief Set the lower bounds on the arcs.
     ///
-    /// Simple constructor (with lower bounds).
+    /// 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 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 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).
-    CycleCanceling( 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), _cost(cost),
-      _supply(digraph, 0), _flow(NULL), _local_flow(false),
-      _potential(NULL), _local_potential(false), _res_graph(NULL),
-      _res_cost(_cost)
-    {
-      // 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];
-        }
+    /// \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>
+    CycleCanceling& 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];
       }
-      _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).
-    CycleCanceling( const Digraph &digraph,
-                    const CapacityMap &capacity,
-                    const CostMap &cost,
-                    Node s, Node t,
-                    Supply flow_value ) :
-      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
-      _supply(digraph, 0), _flow(NULL), _local_flow(false),
-      _potential(NULL), _local_potential(false), _res_graph(NULL),
-      _res_cost(_cost)
-    {
-      _supply[s] =  flow_value;
-      _supply[t] = -flow_value;
-      _valid_supply = true;
-    }
-*/
-    /// Destructor.
-    ~CycleCanceling() {
-      if (_local_flow) delete _flow;
-      if (_local_potential) delete _potential;
-      delete _res_graph;
-    }
-
-    /// \brief Set the flow map.
-    ///
-    /// Set the flow map.
-    ///
-    /// \return \c (*this)
-    CycleCanceling& flowMap(FlowMap &map) {
-      if (_local_flow) {
-        delete _flow;
-        _local_flow = false;
-      }
-      _flow = &map;
       return *this;
     }
 
-    /// \brief Set the potential map.
+    /// \brief Set the upper bounds (capacities) on the arcs.
     ///
-    /// Set the potential map.
+    /// 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).
     ///
-    /// \return \c (*this)
-    CycleCanceling& potentialMap(PotentialMap &map) {
-      if (_local_potential) {
-        delete _potential;
-        _local_potential = false;
+    /// \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>
+    CycleCanceling& upperMap(const UpperMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _upper[_arc_idf[a]] = map[a];
       }
-      _potential = &map;
       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>
+    CycleCanceling& costMap(const CostMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _cost[_arc_idf[a]] =  map[a];
+        _cost[_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>
+    CycleCanceling& 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 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>
+    CycleCanceling& 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;
+    }
+    
+    /// @}
+
     /// \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
+    ///   CycleCanceling<ListDigraph> cc(graph);
+    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
     ///
-    /// \param min_mean_cc Set this parameter to \c true to run the
-    /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
-    /// polynomial, but rather slower in practice.
+    /// 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 and extends the graph.
     ///
-    /// \return \c true if a feasible flow can be found.
-    bool run(bool min_mean_cc = false) {
-      return init() && start(min_mean_cc);
+    /// \param method The cycle-canceling method that will be used.
+    /// For more information, see \ref Method.
+    ///
+    /// \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, Method
+    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
+      ProblemType pt = init();
+      if (pt != OPTIMAL) return pt;
+      start(method);
+      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
+    ///   CycleCanceling<ListDigraph> cs(graph);
+    ///
+    ///   // First run
+    ///   cc.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;
+    ///   cc.costMap(cost).run();
+    ///
+    ///   // Run again from scratch using reset()
+    ///   // (the lower bounds will be set to zero on all arcs)
+    ///   cc.reset();
+    ///   cc.upperMap(capacity).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// \return <tt>(*this)</tt>
+    CycleCanceling& 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;
+        _cost[j] = _forward[j] ? 1 : -1;
+      }
+      for (int j = limit; j != _res_arc_num; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _cost[j] = 0;
+        _cost[_reverse[j]] = 0;
+      }      
+      _have_lower = false;
+      return *this;
     }
 
     /// @}
 
     /// \name Query Functions
-    /// The result of the algorithm can be obtained using these
+    /// The results of the algorithm can be obtained using these
     /// functions.\n
-    /// \ref lemon::CycleCanceling::run() "run()" must be called before
-    /// using them.
+    /// The \ref run() function must be called before using them.
 
     /// @{
 
-    /// \brief Return a const reference to the arc map storing the
-    /// found flow.
+    /// \brief Return the total cost of the found flow.
     ///
-    /// Return a const reference to the arc map storing 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
+    ///   cc.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.
-    const FlowMap& flowMap() const {
-      return *_flow;
+    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>(_cost[i]));
+      }
+      return c;
     }
 
-    /// \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;
+#ifndef DOXYGEN
+    Cost totalCost() const {
+      return totalCost<Cost>();
     }
+#endif
 
     /// \brief Return the flow on the given arc.
     ///
-    /// 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.
-    Capacity flow(const Arc& arc) const {
-      return (*_flow)[arc];
+    Value flow(const Arc& a) const {
+      return _res_cap[_arc_idb[a]];
     }
 
-    /// \brief Return the potential of the given node.
+    /// \brief Return the flow map (the primal solution).
     ///
-    /// Return the potential of the given node.
+    /// 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.
-    Cost potential(const Node& node) const {
-      return (*_potential)[node];
+    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 total cost of the found flow.
+    /// \brief Return the potential (dual value) of the given node.
     ///
-    /// Return the total cost of the found flow. The complexity of the
-    /// function is \f$ O(e) \f$.
+    /// This function returns the potential (dual value) of the
+    /// given node.
     ///
     /// \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] * _cost[e];
-      return c;
+    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.
-    bool init() {
-      if (!_valid_supply) return false;
+    // Initialize the algorithm
+    ProblemType init() {
+      if (_res_node_num <= 1) return INFEASIBLE;
 
-      // Initializing flow and potential maps
-      if (!_flow) {
-        _flow = new FlowMap(_graph);
-        _local_flow = true;
+      // Check the sum of supply values
+      _sum_supply = 0;
+      for (int i = 0; i != _root; ++i) {
+        _sum_supply += _supply[i];
       }
-      if (!_potential) {
-        _potential = new PotentialMap(_graph);
-        _local_potential = true;
+      if (_sum_supply > 0) return INFEASIBLE;
+      
+
+      // Initialize vectors
+      for (int i = 0; i != _res_node_num; ++i) {
+        _pi[i] = 0;
+      }
+      ValueVector excess(_supply);
+      
+      // 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 = _cost[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] && _cost[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];
+        if (ex < 0) max_cap -= ex;
+      }
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_upper[j] >= MAX) _upper[j] = max_cap;
       }
 
-      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
+      // 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]];
+        }
+      }
 
-      // Finding 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();
+      // 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;
+        }
+      } 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;
+    }
+    
+    // Build a StaticDigraph structure containing the current
+    // residual network
+    void buildResidualNetwork() {
+      _arc_vec.clear();
+      _cost_vec.clear();
+      _id_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]);
+          _id_vec.push_back(j);
+        }
+      }
+      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
     }
 
-    bool start(bool min_mean_cc) {
-      if (min_mean_cc)
-        startMinMean();
-      else
-        start();
+    // Execute the algorithm and transform the results
+    void start(Method method) {
+      // Execute the algorithm
+      switch (method) {
+        case SIMPLE_CYCLE_CANCELING:
+          startSimpleCycleCanceling();
+          break;
+        case MINIMUM_MEAN_CYCLE_CANCELING:
+          startMinMeanCycleCanceling();
+          break;
+        case CANCEL_AND_TIGHTEN:
+          startCancelAndTighten();
+          break;
+      }
 
-      // Handling non-zero lower bounds
-      if (_lower) {
-        for (ArcIt e(_graph); e != INVALID; ++e)
-          (*_flow)[e] += (*_lower)[e];
+      // Compute node potentials
+      if (method != SIMPLE_CYCLE_CANCELING) {
+        buildResidualNetwork();
+        typename BellmanFord<StaticDigraph, CostArcMap>
+          ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
+        bf.distMap(_pi_map);
+        bf.init(0);
+        bf.start();
       }
-      return true;
+
+      // 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];
+        }
+      }
     }
 
-    /// \brief Execute the algorithm using \ref BellmanFord.
-    ///
-    /// Execute the algorithm using the \ref BellmanFord
-    /// "Bellman-Ford" algorithm for negative cycle detection with
-    /// successively larger limit for the number of iterations.
-    void start() {
-      typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
-      typename ResDigraph::template NodeMap<int> visited(*_res_graph);
-      std::vector<ResArc> cycle;
-      int node_num = countNodes(_graph);
+    // Execute the "Simple Cycle Canceling" method
+    void startSimpleCycleCanceling() {
+      // Constants for computing the iteration limits
+      const int BF_FIRST_LIMIT  = 2;
+      const double BF_LIMIT_FACTOR = 1.5;
+      
+      typedef VectorMap<StaticDigraph::Arc, Value> FilterMap;
+      typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
+      typedef VectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
+      typedef typename BellmanFord<ResDigraph, CostArcMap>
+        ::template SetDistMap<CostNodeMap>
+        ::template SetPredMap<PredMap>::Create BF;
+      
+      // Build the residual network
+      _arc_vec.clear();
+      _cost_vec.clear();
+      for (int j = 0; j != _res_arc_num; ++j) {
+        _arc_vec.push_back(IntPair(_source[j], _target[j]));
+        _cost_vec.push_back(_cost[j]);
+      }
+      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+      FilterMap filter_map(_res_cap);
+      ResDigraph rgr(_sgr, filter_map);
+      std::vector<int> cycle;
+      std::vector<StaticDigraph::Arc> pred(_res_arc_num);
+      PredMap pred_map(pred);
+      BF bf(rgr, _cost_map);
+      bf.distMap(_pi_map).predMap(pred_map);
 
       int length_bound = BF_FIRST_LIMIT;
       bool optimal = false;
       while (!optimal) {
-        BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
-        bf.predMap(pred);
         bf.init(0);
         int iter_num = 0;
         bool cycle_found = false;
         while (!cycle_found) {
-          int curr_iter_num = iter_num + length_bound <= node_num ?
-                              length_bound : node_num - iter_num;
+          // Perform some iterations of the Bellman-Ford algorithm
+          int curr_iter_num = iter_num + length_bound <= _node_num ?
+            length_bound : _node_num - iter_num;
           iter_num += curr_iter_num;
           int real_iter_num = curr_iter_num;
           for (int i = 0; i < curr_iter_num; ++i) {
@@ -465,89 +841,290 @@
           if (real_iter_num < curr_iter_num) {
             // Optimal flow is found
             optimal = true;
-            // Setting node potentials
-            for (NodeIt n(_graph); n != INVALID; ++n)
-              (*_potential)[n] = bf.dist(n);
             break;
           } else {
-            // Searching for node disjoint negative cycles
-            for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
-              visited[n] = 0;
+            // Search for node disjoint negative cycles
+            std::vector<int> state(_res_node_num, 0);
             int id = 0;
-            for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
-              if (visited[n] > 0) continue;
-              visited[n] = ++id;
-              ResNode u = pred[n] == INVALID ?
-                          INVALID : _res_graph->source(pred[n]);
-              while (u != INVALID && visited[u] == 0) {
-                visited[u] = id;
-                u = pred[u] == INVALID ?
-                    INVALID : _res_graph->source(pred[u]);
+            for (int u = 0; u != _res_node_num; ++u) {
+              if (state[u] != 0) continue;
+              ++id;
+              int v = u;
+              for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
+                   -1 : rgr.id(rgr.source(pred[v]))) {
+                state[v] = id;
               }
-              if (u != INVALID && visited[u] == id) {
-                // Finding the negative cycle
+              if (v != -1 && state[v] == id) {
+                // A negative cycle is found
                 cycle_found = true;
                 cycle.clear();
-                ResArc e = pred[u];
-                cycle.push_back(e);
-                Capacity d = _res_graph->residualCapacity(e);
-                while (_res_graph->source(e) != u) {
-                  cycle.push_back(e = pred[_res_graph->source(e)]);
-                  if (_res_graph->residualCapacity(e) < d)
-                    d = _res_graph->residualCapacity(e);
+                StaticDigraph::Arc a = pred[v];
+                Value d, delta = _res_cap[rgr.id(a)];
+                cycle.push_back(rgr.id(a));
+                while (rgr.id(rgr.source(a)) != v) {
+                  a = pred_map[rgr.source(a)];
+                  d = _res_cap[rgr.id(a)];
+                  if (d < delta) delta = d;
+                  cycle.push_back(rgr.id(a));
                 }
 
-                // Augmenting along the cycle
-                for (int i = 0; i < int(cycle.size()); ++i)
-                  _res_graph->augment(cycle[i], d);
+                // Augment along the cycle
+                for (int i = 0; i < int(cycle.size()); ++i) {
+                  int j = cycle[i];
+                  _res_cap[j] -= delta;
+                  _res_cap[_reverse[j]] += delta;
+                }
               }
             }
           }
 
-          if (!cycle_found)
-            length_bound = length_bound * BF_LIMIT_FACTOR / 100;
+          // Increase iteration limit if no cycle is found
+          if (!cycle_found) {
+            length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
+          }
         }
       }
     }
 
-    /// \brief Execute the algorithm using \ref Howard.
-    ///
-    /// Execute the algorithm using \ref Howard for negative
-    /// cycle detection.
-    void startMinMean() {
-      typedef Path<ResDigraph> ResPath;
-      Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
-      ResPath cycle;
+    // Execute the "Minimum Mean Cycle Canceling" method
+    void startMinMeanCycleCanceling() {
+      typedef SimplePath<StaticDigraph> SPath;
+      typedef typename SPath::ArcIt SPathArcIt;
+      typedef typename Howard<StaticDigraph, CostArcMap>
+        ::template SetPath<SPath>::Create MMC;
+      
+      SPath cycle;
+      MMC mmc(_sgr, _cost_map);
+      mmc.cycle(cycle);
+      buildResidualNetwork();
+      while (mmc.findMinMean() && mmc.cycleLength() < 0) {
+        // Find the cycle
+        mmc.findCycle();
 
-      mmc.cycle(cycle);
-      if (mmc.findMinMean()) {
-        while (mmc.cycleLength() < 0) {
-          // Finding the cycle
-          mmc.findCycle();
+        // Compute delta value
+        Value delta = INF;
+        for (SPathArcIt a(cycle); a != INVALID; ++a) {
+          Value d = _res_cap[_id_vec[_sgr.id(a)]];
+          if (d < delta) delta = d;
+        }
 
-          // Finding the largest flow amount that can be augmented
-          // along the cycle
-          Capacity delta = 0;
-          for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
-            if (delta == 0 || _res_graph->residualCapacity(e) < delta)
-              delta = _res_graph->residualCapacity(e);
+        // Augment along the cycle
+        for (SPathArcIt a(cycle); a != INVALID; ++a) {
+          int j = _id_vec[_sgr.id(a)];
+          _res_cap[j] -= delta;
+          _res_cap[_reverse[j]] += delta;
+        }
+
+        // Rebuild the residual network        
+        buildResidualNetwork();
+      }
+    }
+
+    // Execute the "Cancel And Tighten" method
+    void startCancelAndTighten() {
+      // Constants for the min mean cycle computations
+      const double LIMIT_FACTOR = 1.0;
+      const int MIN_LIMIT = 5;
+
+      // Contruct auxiliary data vectors
+      DoubleVector pi(_res_node_num, 0.0);
+      IntVector level(_res_node_num);
+      CharVector reached(_res_node_num);
+      CharVector processed(_res_node_num);
+      IntVector pred_node(_res_node_num);
+      IntVector pred_arc(_res_node_num);
+      std::vector<int> stack(_res_node_num);
+      std::vector<int> proc_vector(_res_node_num);
+
+      // Initialize epsilon
+      double epsilon = 0;
+      for (int a = 0; a != _res_arc_num; ++a) {
+        if (_res_cap[a] > 0 && -_cost[a] > epsilon)
+          epsilon = -_cost[a];
+      }
+
+      // Start phases
+      Tolerance<double> tol;
+      tol.epsilon(1e-6);
+      int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
+      if (limit < MIN_LIMIT) limit = MIN_LIMIT;
+      int iter = limit;
+      while (epsilon * _res_node_num >= 1) {
+        // Find and cancel cycles in the admissible network using DFS
+        for (int u = 0; u != _res_node_num; ++u) {
+          reached[u] = false;
+          processed[u] = false;
+        }
+        int stack_head = -1;
+        int proc_head = -1;
+        for (int start = 0; start != _res_node_num; ++start) {
+          if (reached[start]) continue;
+
+          // New start node
+          reached[start] = true;
+          pred_arc[start] = -1;
+          pred_node[start] = -1;
+
+          // Find the first admissible outgoing arc
+          double p = pi[start];
+          int a = _first_out[start];
+          int last_out = _first_out[start+1];
+          for (; a != last_out && (_res_cap[a] == 0 ||
+               !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+          if (a == last_out) {
+            processed[start] = true;
+            proc_vector[++proc_head] = start;
+            continue;
+          }
+          stack[++stack_head] = a;
+
+          while (stack_head >= 0) {
+            int sa = stack[stack_head];
+            int u = _source[sa];
+            int v = _target[sa];
+
+            if (!reached[v]) {
+              // A new node is reached
+              reached[v] = true;
+              pred_node[v] = u;
+              pred_arc[v] = sa;
+              p = pi[v];
+              a = _first_out[v];
+              last_out = _first_out[v+1];
+              for (; a != last_out && (_res_cap[a] == 0 ||
+                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+              stack[++stack_head] = a == last_out ? -1 : a;
+            } else {
+              if (!processed[v]) {
+                // A cycle is found
+                int n, w = u;
+                Value d, delta = _res_cap[sa];
+                for (n = u; n != v; n = pred_node[n]) {
+                  d = _res_cap[pred_arc[n]];
+                  if (d <= delta) {
+                    delta = d;
+                    w = pred_node[n];
+                  }
+                }
+
+                // Augment along the cycle
+                _res_cap[sa] -= delta;
+                _res_cap[_reverse[sa]] += delta;
+                for (n = u; n != v; n = pred_node[n]) {
+                  int pa = pred_arc[n];
+                  _res_cap[pa] -= delta;
+                  _res_cap[_reverse[pa]] += delta;
+                }
+                for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
+                  --stack_head;
+                  reached[n] = false;
+                }
+                u = w;
+              }
+              v = u;
+
+              // Find the next admissible outgoing arc
+              p = pi[v];
+              a = stack[stack_head] + 1;
+              last_out = _first_out[v+1];
+              for (; a != last_out && (_res_cap[a] == 0 ||
+                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+              stack[stack_head] = a == last_out ? -1 : a;
+            }
+
+            while (stack_head >= 0 && stack[stack_head] == -1) {
+              processed[v] = true;
+              proc_vector[++proc_head] = v;
+              if (--stack_head >= 0) {
+                // Find the next admissible outgoing arc
+                v = _source[stack[stack_head]];
+                p = pi[v];
+                a = stack[stack_head] + 1;
+                last_out = _first_out[v+1];
+                for (; a != last_out && (_res_cap[a] == 0 ||
+                     !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+                stack[stack_head] = a == last_out ? -1 : a;
+              }
+            }
+          }
+        }
+
+        // Tighten potentials and epsilon
+        if (--iter > 0) {
+          for (int u = 0; u != _res_node_num; ++u) {
+            level[u] = 0;
+          }
+          for (int i = proc_head; i > 0; --i) {
+            int u = proc_vector[i];
+            double p = pi[u];
+            int l = level[u] + 1;
+            int last_out = _first_out[u+1];
+            for (int a = _first_out[u]; a != last_out; ++a) {
+              int v = _target[a];
+              if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
+                  l > level[v]) level[v] = l;
+            }
           }
 
-          // Augmenting along the cycle
-          for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
-            _res_graph->augment(e, delta);
+          // Modify potentials
+          double q = std::numeric_limits<double>::max();
+          for (int u = 0; u != _res_node_num; ++u) {
+            int lu = level[u];
+            double p, pu = pi[u];
+            int last_out = _first_out[u+1];
+            for (int a = _first_out[u]; a != last_out; ++a) {
+              if (_res_cap[a] == 0) continue;
+              int v = _target[a];
+              int ld = lu - level[v];
+              if (ld > 0) {
+                p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
+                if (p < q) q = p;
+              }
+            }
+          }
+          for (int u = 0; u != _res_node_num; ++u) {
+            pi[u] -= q * level[u];
+          }
 
-          // Finding the minimum cycle mean for the modified residual
-          // digraph
-          if (!mmc.findMinMean()) break;
+          // Modify epsilon
+          epsilon = 0;
+          for (int u = 0; u != _res_node_num; ++u) {
+            double curr, pu = pi[u];
+            int last_out = _first_out[u+1];
+            for (int a = _first_out[u]; a != last_out; ++a) {
+              if (_res_cap[a] == 0) continue;
+              curr = _cost[a] + pu - pi[_target[a]];
+              if (-curr > epsilon) epsilon = -curr;
+            }
+          }
+        } else {
+          typedef Howard<StaticDigraph, CostArcMap> MMC;
+          typedef typename BellmanFord<StaticDigraph, CostArcMap>
+            ::template SetDistMap<CostNodeMap>::Create BF;
+
+          // Set epsilon to the minimum cycle mean
+          buildResidualNetwork();
+          MMC mmc(_sgr, _cost_map);
+          mmc.findMinMean();
+          epsilon = -mmc.cycleMean();
+          Cost cycle_cost = mmc.cycleLength();
+          int cycle_size = mmc.cycleArcNum();
+          
+          // Compute feasible potentials for the current epsilon
+          for (int i = 0; i != int(_cost_vec.size()); ++i) {
+            _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
+          }
+          BF bf(_sgr, _cost_map);
+          bf.distMap(_pi_map);
+          bf.init(0);
+          bf.start();
+          for (int u = 0; u != _res_node_num; ++u) {
+            pi[u] = static_cast<double>(_pi[u]) / cycle_size;
+          }
+        
+          iter = limit;
         }
       }
-
-      // Computing node potentials
-      BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
-      bf.init(0); bf.start();
-      for (NodeIt n(_graph); n != INVALID; ++n)
-        (*_potential)[n] = bf.dist(n);
     }
 
   }; //class CycleCanceling