# HG changeset patch
# User kpeter
# Date 1196859799 0
# Node ID 716024e7c08087cc9dcbdf6a146c575d721de2fd
# Parent  edad4c3e926d280402a62e77607f75f4c036ffb3
Redesigned CapacityScaling algorithm with almost the same interface.
The new version does not use the ResidualGraphAdaptor for performance reasons.
Scaling can be enabled and disabled with a parameter of the run() function.

diff -r edad4c3e926d -r 716024e7c080 lemon/capacity_scaling.h
--- a/lemon/capacity_scaling.h	Wed Dec 05 12:57:24 2007 +0000
+++ b/lemon/capacity_scaling.h	Wed Dec 05 13:03:19 2007 +0000
@@ -25,17 +25,9 @@
 /// \brief The capacity scaling algorithm for finding a minimum cost
 /// flow.
 
+#include <lemon/graph_adaptor.h>
+#include <lemon/bin_heap.h>
 #include <vector>
-#include <lemon/graph_adaptor.h>
-#include <lemon/dijkstra.h>
-
-#define WITH_SCALING
-
-#ifdef WITH_SCALING
-#define SCALING_FACTOR	  2
-#endif
-
-//#define _DEBUG_ITER_
 
 namespace lemon {
 
@@ -46,9 +38,9 @@
   /// successive shortest path algorithm for finding a minimum cost
   /// flow.
   ///
-  /// \ref lemon::CapacityScaling "CapacityScaling" implements the
-  /// capacity scaling version of the successive shortest path
-  /// algorithm for finding a minimum cost flow.
+  /// \ref CapacityScaling implements the capacity scaling version
+  /// of the successive shortest path algorithm for finding a minimum
+  /// cost flow.
   ///
   /// \param Graph The directed graph type the algorithm runs on.
   /// \param LowerMap The type of the lower bound map.
@@ -58,20 +50,20 @@
   ///
   /// \warning
   /// - Edge capacities and costs should be nonnegative integers.
-  ///	However \c CostMap::Value should be signed type.
+  ///   However \c CostMap::Value should be signed type.
   /// - Supply values should be signed integers.
   /// - \c LowerMap::Value must be convertible to
-  ///	\c CapacityMap::Value and \c CapacityMap::Value must be
-  ///	convertible to \c SupplyMap::Value.
+  ///   \c CapacityMap::Value and \c CapacityMap::Value must be
+  ///   convertible to \c SupplyMap::Value.
   ///
   /// \author Peter Kovacs
 
   template < typename Graph,
-	     typename LowerMap = typename Graph::template EdgeMap<int>,
-	     typename CapacityMap = LowerMap,
-	     typename CostMap = typename Graph::template EdgeMap<int>,
-	     typename SupplyMap = typename Graph::template NodeMap
-				  <typename CapacityMap::Value> >
+             typename LowerMap = typename Graph::template EdgeMap<int>,
+             typename CapacityMap = LowerMap,
+             typename CostMap = typename Graph::template EdgeMap<int>,
+             typename SupplyMap = typename Graph::template NodeMap
+                                  <typename CapacityMap::Value> >
   class CapacityScaling
   {
     typedef typename Graph::Node Node;
@@ -87,16 +79,16 @@
     typedef typename SupplyMap::Value Supply;
     typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
     typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
-
-    typedef ResGraphAdaptor< const Graph, Capacity,
-			     CapacityRefMap, CapacityRefMap > ResGraph;
-    typedef typename ResGraph::Node ResNode;
-    typedef typename ResGraph::NodeIt ResNodeIt;
-    typedef typename ResGraph::Edge ResEdge;
-    typedef typename ResGraph::EdgeIt ResEdgeIt;
+    typedef typename Graph::template NodeMap<Edge> PredMap;
 
   public:
 
+    /// \brief Type to enable or disable capacity scaling.
+    enum ScalingEnum {
+      WITH_SCALING = 0,
+      WITHOUT_SCALING = -1
+    };
+
     /// \brief The type of the flow map.
     typedef CapacityRefMap FlowMap;
     /// \brief The type of the potential map.
@@ -104,151 +96,129 @@
 
   protected:
 
-    /// \brief Map adaptor class for handling reduced edge costs.
-    class ReducedCostMap : public MapBase<ResEdge, Cost>
+    /// \brief Special implementation of the \ref Dijkstra algorithm
+    /// for finding shortest paths in the residual network of the graph
+    /// with respect to the reduced edge costs and modifying the
+    /// node potentials according to the distance of the nodes.
+    class ResidualDijkstra
     {
-    private:
+      typedef typename Graph::template NodeMap<Cost> CostNodeMap;
+      typedef typename Graph::template NodeMap<Edge> PredMap;
 
-      const ResGraph &gr;
-      const CostMap &cost_map;
-      const PotentialMap &pot_map;
+      typedef typename Graph::template NodeMap<int> HeapCrossRef;
+      typedef BinHeap<Cost, HeapCrossRef> Heap;
+
+    protected:
+
+      /// \brief The directed graph the algorithm runs on.
+      const Graph &graph;
+
+      /// \brief The flow map.
+      const FlowMap &flow;
+      /// \brief The residual capacity map.
+      const CapacityRefMap &res_cap;
+      /// \brief The cost map.
+      const CostMap &cost;
+      /// \brief The excess map.
+      const SupplyRefMap &excess;
+
+      /// \brief The potential map.
+      PotentialMap &potential;
+
+      /// \brief The distance map.
+      CostNodeMap dist;
+      /// \brief The map of predecessors edges.
+      PredMap &pred;
+      /// \brief The processed (i.e. permanently labeled) nodes.
+      std::vector<Node> proc_nodes;
 
     public:
 
-      ReducedCostMap( const ResGraph &_gr,
-		      const CostMap &_cost,
-		      const PotentialMap &_pot ) :
-	gr(_gr), cost_map(_cost), pot_map(_pot) {}
+      /// \brief The constructor of the class.
+      ResidualDijkstra( const Graph &_graph,
+                        const FlowMap &_flow,
+                        const CapacityRefMap &_res_cap,
+                        const CostMap &_cost,
+                        const SupplyMap &_excess,
+                        PotentialMap &_potential,
+                        PredMap &_pred ) :
+        graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),
+        excess(_excess), potential(_potential), dist(_graph),
+        pred(_pred)
+      {}
 
-      Cost operator[](const ResEdge &e) const {
-	return ResGraph::forward(e) ?
-	   cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
-	  -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+      /// \brief Runs the algorithm from the given source node.
+      Node run(Node s, Capacity delta) {
+        HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);
+        Heap heap(heap_cross_ref);
+        heap.push(s, 0);
+        pred[s] = INVALID;
+        proc_nodes.clear();
+
+        // Processing nodes
+        while (!heap.empty() && excess[heap.top()] > -delta) {
+          Node u = heap.top(), v;
+          Cost d = heap.prio() - potential[u], nd;
+          dist[u] = heap.prio();
+          heap.pop();
+          proc_nodes.push_back(u);
+
+          // Traversing outgoing edges
+          for (OutEdgeIt e(graph, u); e != INVALID; ++e) {
+            if (res_cap[e] >= delta) {
+              v = graph.target(e);
+              switch(heap.state(v)) {
+              case Heap::PRE_HEAP:
+                heap.push(v, d + cost[e] + potential[v]);
+                pred[v] = e;
+                break;
+              case Heap::IN_HEAP:
+                nd = d + cost[e] + potential[v];
+                if (nd < heap[v]) {
+                  heap.decrease(v, nd);
+                  pred[v] = e;
+                }
+                break;
+              case Heap::POST_HEAP:
+                break;
+              }
+            }
+          }
+
+          // Traversing incoming edges
+          for (InEdgeIt e(graph, u); e != INVALID; ++e) {
+            if (flow[e] >= delta) {
+              v = graph.source(e);
+              switch(heap.state(v)) {
+              case Heap::PRE_HEAP:
+                heap.push(v, d - cost[e] + potential[v]);
+                pred[v] = e;
+                break;
+              case Heap::IN_HEAP:
+                nd = d - cost[e] + potential[v];
+                if (nd < heap[v]) {
+                  heap.decrease(v, nd);
+                  pred[v] = e;
+                }
+                break;
+              case Heap::POST_HEAP:
+                break;
+              }
+            }
+          }
+        }
+        if (heap.empty()) return INVALID;
+
+        // Updating potentials of processed nodes
+        Node t = heap.top();
+        Cost dt = heap.prio();
+        for (int i = 0; i < proc_nodes.size(); ++i)
+          potential[proc_nodes[i]] -= dist[proc_nodes[i]] - dt;
+
+        return t;
       }
 
-    }; //class ReducedCostMap
-
-    /// \brief Map class for the \ref lemon::Dijkstra "Dijkstra"
-    /// algorithm to update node potentials.
-    class PotentialUpdateMap : public MapBase<ResNode, Cost>
-    {
-    private:
-
-      PotentialMap *pot;
-      typedef std::pair<ResNode, Cost> Pair;
-      std::vector<Pair> data;
-
-    public:
-
-      void potentialMap(PotentialMap &_pot) {
-	pot = &_pot;
-      }
-
-      void init() {
-	data.clear();
-      }
-
-      void set(const ResNode &n, const Cost &v) {
-	data.push_back(Pair(n, v));
-      }
-
-      void update() {
-	Cost val = data[data.size()-1].second;
-	for (int i = 0; i < data.size()-1; ++i)
-	  (*pot)[data[i].first] += val - data[i].second;
-      }
-
-    }; //class PotentialUpdateMap
-
-#ifdef WITH_SCALING
-    /// \brief Map adaptor class for identifing deficit nodes.
-    class DeficitBoolMap : public MapBase<ResNode, bool>
-    {
-    private:
-
-      const SupplyRefMap &imb;
-      const Capacity &delta;
-
-    public:
-
-      DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
-	imb(_imb), delta(_delta) {}
-
-      bool operator[](const ResNode &n) const {
-	return imb[n] <= -delta;
-      }
-
-    }; //class DeficitBoolMap
-
-    /// \brief Map adaptor class for filtering edges with at least
-    /// \c delta residual capacity.
-    class DeltaFilterMap : public MapBase<ResEdge, bool>
-    {
-    private:
-
-      const ResGraph &gr;
-      const Capacity &delta;
-
-    public:
-
-      DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
-	gr(_gr), delta(_delta) {}
-
-      bool operator[](const ResEdge &e) const {
-	return gr.rescap(e) >= delta;
-      }
-
-    }; //class DeltaFilterMap
-
-    typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
-      DeltaResGraph;
-
-    /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
-    class ResDijkstraTraits :
-      public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
-    {
-    public:
-
-      typedef PotentialUpdateMap DistMap;
-
-      static DistMap *createDistMap(const DeltaResGraph&) {
-	return new DistMap();
-      }
-
-    }; //class ResDijkstraTraits
-
-#else //WITHOUT_CAPACITY_SCALING
-    /// \brief Map adaptor class for identifing deficit nodes.
-    class DeficitBoolMap : public MapBase<ResNode, bool>
-    {
-    private:
-
-      const SupplyRefMap &imb;
-
-    public:
-
-      DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
-
-      bool operator[](const ResNode &n) const {
-	return imb[n] < 0;
-      }
-
-    }; //class DeficitBoolMap
-
-    /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
-    class ResDijkstraTraits :
-      public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
-    {
-    public:
-
-      typedef PotentialUpdateMap DistMap;
-
-      static DistMap *createDistMap(const ResGraph&) {
-	return new DistMap();
-      }
-
-    }; //class ResDijkstraTraits
-#endif
+    }; //class ResidualDijkstra
 
   protected:
 
@@ -269,57 +239,30 @@
     FlowMap flow;
     /// \brief The potential node map.
     PotentialMap potential;
-    /// \brief The residual graph.
-    ResGraph res_graph;
-    /// \brief The reduced cost map.
-    ReducedCostMap red_cost;
 
-    /// \brief The imbalance map.
-    SupplyRefMap imbalance;
-    /// \brief The excess nodes.
-    std::vector<ResNode> excess_nodes;
+    /// \brief The residual capacity map.
+    CapacityRefMap res_cap;
+    /// \brief The excess map.
+    SupplyRefMap excess;
+    /// \brief The excess nodes (i.e. nodes with positive excess).
+    std::vector<Node> excess_nodes;
     /// \brief The index of the next excess node.
     int next_node;
 
-#ifdef WITH_SCALING
-    typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
-      ResDijkstra;
-    /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
-    /// shortest paths in the residual graph with respect to the
-    /// reduced edge costs.
-    ResDijkstra dijkstra;
-
+    /// \brief The scaling status (enabled or disabled).
+    ScalingEnum scaling;
     /// \brief The delta parameter used for capacity scaling.
     Capacity delta;
-    /// \brief Edge filter map for filtering edges with at least
-    /// \c delta residual capacity.
-    DeltaFilterMap delta_filter;
-    /// \brief The delta residual graph (i.e. the subgraph of the
-    /// residual graph consisting of edges with at least \c delta
-    /// residual capacity).
-    DeltaResGraph dres_graph;
-    /// \brief Map for identifing deficit nodes.
-    DeficitBoolMap delta_deficit;
+    /// \brief The maximum number of phases.
+    Capacity phase_num;
+    /// \brief The deficit nodes.
+    std::vector<Node> deficit_nodes;
 
-    /// \brief The deficit nodes.
-    std::vector<ResNode> deficit_nodes;
-
-#else //WITHOUT_CAPACITY_SCALING
-    typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
-      ResDijkstra;
-    /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
-    /// shortest paths in the residual graph with respect to the
-    /// reduced edge costs.
-    ResDijkstra dijkstra;
-    /// \brief Map for identifing deficit nodes.
-    DeficitBoolMap has_deficit;
-#endif
-
-    /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
-    typename ResDijkstra::PredMap pred;
-    /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
-    /// update node potentials.
-    PotentialUpdateMap updater;
+    /// \brief Implementation of the \ref Dijkstra algorithm for
+    /// finding augmenting shortest paths in the residual network.
+    ResidualDijkstra dijkstra;
+    /// \brief The map of predecessors edges.
+    PredMap pred;
 
   public :
 
@@ -333,35 +276,27 @@
     /// \param _cost The cost (length) values of the edges.
     /// \param _supply The supply values of the nodes (signed).
     CapacityScaling( const Graph &_graph,
-		     const LowerMap &_lower,
-		     const CapacityMap &_capacity,
-		     const CostMap &_cost,
-		     const SupplyMap &_supply ) :
+                     const LowerMap &_lower,
+                     const CapacityMap &_capacity,
+                     const CostMap &_cost,
+                     const SupplyMap &_supply ) :
       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
       supply(_graph), flow(_graph, 0), potential(_graph, 0),
-      res_graph(_graph, capacity, flow),
-      red_cost(res_graph, cost, potential), imbalance(_graph),
-#ifdef WITH_SCALING
-      delta(0), delta_filter(res_graph, delta),
-      dres_graph(res_graph, delta_filter),
-      dijkstra(dres_graph, red_cost), pred(dres_graph),
-      delta_deficit(imbalance, delta)
-#else //WITHOUT_CAPACITY_SCALING
-      dijkstra(res_graph, red_cost), pred(res_graph),
-      has_deficit(imbalance)
-#endif
+      res_cap(_graph), excess(_graph), pred(_graph),
+      dijkstra(graph, flow, res_cap, cost, excess, potential, pred)
     {
       // Removing nonzero lower bounds
       capacity = subMap(_capacity, _lower);
+      res_cap = capacity;
       Supply sum = 0;
       for (NodeIt n(graph); n != INVALID; ++n) {
-	Supply s = _supply[n];
-	for (InEdgeIt e(graph, n); e != INVALID; ++e)
-	  s += _lower[e];
-	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
-	  s -= _lower[e];
-	supply[n] = s;
-	sum += s;
+        Supply s = _supply[n];
+        for (InEdgeIt e(graph, n); e != INVALID; ++e)
+          s += _lower[e];
+        for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+          s -= _lower[e];
+        supply[n] = s;
+        sum += s;
       }
       valid_supply = sum == 0;
     }
@@ -375,22 +310,13 @@
     /// \param _cost The cost (length) values of the edges.
     /// \param _supply The supply values of the nodes (signed).
     CapacityScaling( const Graph &_graph,
-		     const CapacityMap &_capacity,
-		     const CostMap &_cost,
-		     const SupplyMap &_supply ) :
+                     const CapacityMap &_capacity,
+                     const CostMap &_cost,
+                     const SupplyMap &_supply ) :
       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
       supply(_supply), flow(_graph, 0), potential(_graph, 0),
-      res_graph(_graph, capacity, flow),
-      red_cost(res_graph, cost, potential), imbalance(_graph),
-#ifdef WITH_SCALING
-      delta(0), delta_filter(res_graph, delta),
-      dres_graph(res_graph, delta_filter),
-      dijkstra(dres_graph, red_cost), pred(dres_graph),
-      delta_deficit(imbalance, delta)
-#else //WITHOUT_CAPACITY_SCALING
-      dijkstra(res_graph, red_cost), pred(res_graph),
-      has_deficit(imbalance)
-#endif
+      res_cap(_capacity), excess(_graph), pred(_graph),
+      dijkstra(graph, flow, res_cap, cost, excess, potential)
     {
       // Checking the sum of supply values
       Supply sum = 0;
@@ -412,36 +338,28 @@
     /// to node \c _t (i.e. the supply of \c _s and the demand of
     /// \c _t).
     CapacityScaling( const Graph &_graph,
-		     const LowerMap &_lower,
-		     const CapacityMap &_capacity,
-		     const CostMap &_cost,
-		     Node _s, Node _t,
-		     Supply _flow_value ) :
+                     const LowerMap &_lower,
+                     const CapacityMap &_capacity,
+                     const CostMap &_cost,
+                     Node _s, Node _t,
+                     Supply _flow_value ) :
       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
       supply(_graph), flow(_graph, 0), potential(_graph, 0),
-      res_graph(_graph, capacity, flow),
-      red_cost(res_graph, cost, potential), imbalance(_graph),
-#ifdef WITH_SCALING
-      delta(0), delta_filter(res_graph, delta),
-      dres_graph(res_graph, delta_filter),
-      dijkstra(dres_graph, red_cost), pred(dres_graph),
-      delta_deficit(imbalance, delta)
-#else //WITHOUT_CAPACITY_SCALING
-      dijkstra(res_graph, red_cost), pred(res_graph),
-      has_deficit(imbalance)
-#endif
+      res_cap(_graph), excess(_graph), pred(_graph),
+      dijkstra(graph, flow, res_cap, cost, excess, potential)
     {
       // Removing nonzero lower bounds
       capacity = subMap(_capacity, _lower);
+      res_cap = capacity;
       for (NodeIt n(graph); n != INVALID; ++n) {
-	Supply s = 0;
-	if (n == _s) s =  _flow_value;
-	if (n == _t) s = -_flow_value;
-	for (InEdgeIt e(graph, n); e != INVALID; ++e)
-	  s += _lower[e];
-	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
-	  s -= _lower[e];
-	supply[n] = s;
+        Supply s = 0;
+        if (n == _s) s =  _flow_value;
+        if (n == _t) s = -_flow_value;
+        for (InEdgeIt e(graph, n); e != INVALID; ++e)
+          s += _lower[e];
+        for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+          s -= _lower[e];
+        supply[n] = s;
       }
       valid_supply = true;
     }
@@ -459,23 +377,14 @@
     /// to node \c _t (i.e. the supply of \c _s and the demand of
     /// \c _t).
     CapacityScaling( const Graph &_graph,
-		     const CapacityMap &_capacity,
-		     const CostMap &_cost,
-		     Node _s, Node _t,
-		     Supply _flow_value ) :
+                     const CapacityMap &_capacity,
+                     const CostMap &_cost,
+                     Node _s, Node _t,
+                     Supply _flow_value ) :
       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
       supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
-      res_graph(_graph, capacity, flow),
-      red_cost(res_graph, cost, potential), imbalance(_graph),
-#ifdef WITH_SCALING
-      delta(0), delta_filter(res_graph, delta),
-      dres_graph(res_graph, delta_filter),
-      dijkstra(dres_graph, red_cost), pred(dres_graph),
-      delta_deficit(imbalance, delta)
-#else //WITHOUT_CAPACITY_SCALING
-      dijkstra(res_graph, red_cost), pred(res_graph),
-      has_deficit(imbalance)
-#endif
+      res_cap(_capacity), excess(_graph), pred(_graph),
+      dijkstra(graph, flow, res_cap, cost, excess, potential)
     {
       supply[_s] =  _flow_value;
       supply[_t] = -_flow_value;
@@ -511,7 +420,7 @@
     Cost totalCost() const {
       Cost c = 0;
       for (EdgeIt e(graph); e != INVALID; ++e)
-	c += flow[e] * cost[e];
+        c += flow[e] * cost[e];
       return c;
     }
 
@@ -519,190 +428,216 @@
     ///
     /// Runs the algorithm.
     ///
+    /// \param scaling_mode The scaling mode. In case of WITH_SCALING
+    /// capacity scaling is enabled in the algorithm (this is the
+    /// default value) otherwise it is disabled.
+    /// If the maximum edge capacity and/or the amount of total supply
+    /// is small, the algorithm could be faster without scaling.
+    ///
     /// \return \c true if a feasible flow can be found.
-    bool run() {
-      return init() && start();
+    bool run(int scaling_mode = WITH_SCALING) {
+      return init(scaling_mode) && start();
     }
 
   protected:
 
     /// \brief Initializes the algorithm.
-    bool init() {
+    bool init(int scaling_mode) {
       if (!valid_supply) return false;
-      imbalance = supply;
+      excess = supply;
 
-      // Initalizing Dijkstra class
-      updater.potentialMap(potential);
-      dijkstra.distMap(updater).predMap(pred);
-
-#ifdef WITH_SCALING
       // Initilaizing delta value
-      Supply max_sup = 0, max_dem = 0;
-      for (NodeIt n(graph); n != INVALID; ++n) {
-	if (supply[n] >  max_sup) max_sup =  supply[n];
-	if (supply[n] < -max_dem) max_dem = -supply[n];
-      }
-      if (max_dem < max_sup) max_sup = max_dem;
-      for ( delta = 1; SCALING_FACTOR * delta < max_sup;
-	    delta *= SCALING_FACTOR ) ;
-#endif
-      return true;
-    }
-
-#ifdef WITH_SCALING
-    /// \brief Executes the capacity scaling version of the successive
-    /// shortest path algorithm.
-    bool start() {
-      typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
-      typedef typename DeltaResGraph::Edge DeltaResEdge;
-#ifdef _DEBUG_ITER_
-      int dijk_num = 0;
-#endif
-
-      // Processing capacity scaling phases
-      ResNode s, t;
-      for ( ; delta >= 1; delta = delta < SCALING_FACTOR && delta > 1 ?
-				  1 : delta / SCALING_FACTOR )
-      {
-	// Saturating edges not satisfying the optimality condition
-	Capacity r;
-	for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
-	  if (red_cost[e] < 0) {
-	    res_graph.augment(e, r = res_graph.rescap(e));
-	    imbalance[dres_graph.source(e)] -= r;
-	    imbalance[dres_graph.target(e)] += r;
-	  }
-	}
-
-	// Finding excess nodes
-	excess_nodes.clear();
-	deficit_nodes.clear();
-	for (ResNodeIt n(res_graph); n != INVALID; ++n) {
-	  if (imbalance[n] >=  delta) excess_nodes.push_back(n);
-	  if (imbalance[n] <= -delta) deficit_nodes.push_back(n);
-	}
-	next_node = 0;
-
-	// Finding shortest paths
-	while (next_node < excess_nodes.size()) {
-	  // Checking deficit nodes
-	  if (delta > 1) {
-	    bool delta_def = false;
-	    for (int i = 0; i < deficit_nodes.size(); ++i) {
-	      if (imbalance[deficit_nodes[i]] <= -delta) {
-		delta_def = true;
-		break;
-	      }
-	    }
-	    if (!delta_def) break;
-	  }
-
-	  // Running Dijkstra
-	  s = excess_nodes[next_node];
-	  updater.init();
-	  dijkstra.init();
-	  dijkstra.addSource(s);
-#ifdef _DEBUG_ITER_
-	  ++dijk_num;
-#endif
-	  if ((t = dijkstra.start(delta_deficit)) == INVALID) {
-	    if (delta > 1) {
-	      ++next_node;
-	      continue;
-	    }
-	    return false;
-	  }
-
-	  // Updating node potentials
-	  updater.update();
-
-	  // Augment along a shortest path from s to t
-	  Capacity d = imbalance[s] < -imbalance[t] ?
-	    imbalance[s] : -imbalance[t];
-	  ResNode u = t;
-	  ResEdge e;
-	  if (d > delta) {
-	    while ((e = pred[u]) != INVALID) {
-	      if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
-	      u = dres_graph.source(e);
-	    }
-	  }
-	  u = t;
-	  while ((e = pred[u]) != INVALID) {
-	    res_graph.augment(e, d);
-	    u = dres_graph.source(e);
-	  }
-	  imbalance[s] -= d;
-	  imbalance[t] += d;
-	  if (imbalance[s] < delta) ++next_node;
-	}
-      }
-#ifdef _DEBUG_ITER_
-      std::cout << "Capacity Scaling algorithm finished with running Dijkstra algorithm "
-	<< dijk_num << " times." << std::endl;
-#endif
-
-      // Handling nonzero lower bounds
-      if (lower) {
-	for (EdgeIt e(graph); e != INVALID; ++e)
-	  flow[e] += (*lower)[e];
+      if (scaling_mode == WITH_SCALING) {
+        // With scaling
+        Supply max_sup = 0, max_dem = 0;
+        for (NodeIt n(graph); n != INVALID; ++n) {
+          if ( supply[n] > max_sup) max_sup =  supply[n];
+          if (-supply[n] > max_dem) max_dem = -supply[n];
+        }
+        if (max_dem < max_sup) max_sup = max_dem;
+        phase_num = 0;
+        for (delta = 1; 2 * delta <= max_sup; delta *= 2)
+          ++phase_num;
+      } else {
+        // Without scaling
+        delta = 1;
       }
       return true;
     }
 
-#else //WITHOUT_CAPACITY_SCALING
+    /// \brief Executes the algorithm.
+    bool start() {
+      if (delta > 1)
+        return startWithScaling();
+      else
+        return startWithoutScaling();
+    }
+
+    /// \brief Executes the capacity scaling version of the successive
+    /// shortest path algorithm.
+    bool startWithScaling() {
+      // Processing capacity scaling phases
+      Node s, t;
+      int phase_cnt = 0;
+      int factor = 4;
+      while (true) {
+        // Saturating all edges not satisfying the optimality condition
+        for (EdgeIt e(graph); e != INVALID; ++e) {
+          Node u = graph.source(e), v = graph.target(e);
+          Cost c = cost[e] - potential[u] + potential[v];
+          if (c < 0 && res_cap[e] >= delta) {
+            excess[u] -= res_cap[e];
+            excess[v] += res_cap[e];
+            flow[e] = capacity[e];
+            res_cap[e] = 0;
+          }
+          else if (c > 0 && flow[e] >= delta) {
+            excess[u] += flow[e];
+            excess[v] -= flow[e];
+            flow[e] = 0;
+            res_cap[e] = capacity[e];
+          }
+        }
+
+        // Finding excess nodes and deficit nodes
+        excess_nodes.clear();
+        deficit_nodes.clear();
+        for (NodeIt n(graph); n != INVALID; ++n) {
+          if (excess[n] >=  delta) excess_nodes.push_back(n);
+          if (excess[n] <= -delta) deficit_nodes.push_back(n);
+        }
+        next_node = 0;
+
+        // Finding augmenting shortest paths
+        while (next_node < excess_nodes.size()) {
+          // Checking deficit nodes
+          if (delta > 1) {
+            bool delta_deficit = false;
+            for (int i = 0; i < deficit_nodes.size(); ++i) {
+              if (excess[deficit_nodes[i]] <= -delta) {
+                delta_deficit = true;
+                break;
+              }
+            }
+            if (!delta_deficit) break;
+          }
+
+          // Running Dijkstra
+          s = excess_nodes[next_node];
+          if ((t = dijkstra.run(s, delta)) == INVALID) {
+            if (delta > 1) {
+              ++next_node;
+              continue;
+            }
+            return false;
+          }
+
+          // Augmenting along a shortest path from s to t.
+          Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
+          Node u = t;
+          Edge e;
+          if (d > delta) {
+            while ((e = pred[u]) != INVALID) {
+              Capacity rc;
+              if (u == graph.target(e)) {
+                rc = res_cap[e];
+                u = graph.source(e);
+              } else {
+                rc = flow[e];
+                u = graph.target(e);
+              }
+              if (rc < d) d = rc;
+            }
+          }
+          u = t;
+          while ((e = pred[u]) != INVALID) {
+            if (u == graph.target(e)) {
+              flow[e] += d;
+              res_cap[e] -= d;
+              u = graph.source(e);
+            } else {
+              flow[e] -= d;
+              res_cap[e] += d;
+              u = graph.target(e);
+            }
+          }
+          excess[s] -= d;
+          excess[t] += d;
+
+          if (excess[s] < delta) ++next_node;
+        }
+
+        if (delta == 1) break;
+        if (++phase_cnt > phase_num / 4) factor = 2;
+        delta = delta <= factor ? 1 : delta / factor;
+      }
+
+      // Handling nonzero lower bounds
+      if (lower) {
+        for (EdgeIt e(graph); e != INVALID; ++e)
+          flow[e] += (*lower)[e];
+      }
+      return true;
+    }
+
     /// \brief Executes the successive shortest path algorithm without
     /// capacity scaling.
-    bool start() {
+    bool startWithoutScaling() {
       // Finding excess nodes
-      for (ResNodeIt n(res_graph); n != INVALID; ++n) {
-	if (imbalance[n] > 0) excess_nodes.push_back(n);
+      for (NodeIt n(graph); n != INVALID; ++n) {
+        if (excess[n] > 0) excess_nodes.push_back(n);
       }
       if (excess_nodes.size() == 0) return true;
       next_node = 0;
 
       // Finding shortest paths
-      ResNode s, t;
-      while ( imbalance[excess_nodes[next_node]] > 0 ||
-	      ++next_node < excess_nodes.size() )
+      Node s, t;
+      while ( excess[excess_nodes[next_node]] > 0 ||
+              ++next_node < excess_nodes.size() )
       {
-	// Running Dijkstra
-	s = excess_nodes[next_node];
-	updater.init();
-	dijkstra.init();
-	dijkstra.addSource(s);
-	if ((t = dijkstra.start(has_deficit)) == INVALID)
-	  return false;
+        // Running Dijkstra
+        s = excess_nodes[next_node];
+        if ((t = dijkstra.run(s, 1)) == INVALID)
+          return false;
 
-	// Updating node potentials
-	updater.update();
-
-	// Augmenting along a shortest path from s to t
-	Capacity delta = imbalance[s] < -imbalance[t] ?
-	  imbalance[s] : -imbalance[t];
-	ResNode u = t;
-	ResEdge e;
-	while ((e = pred[u]) != INVALID) {
-	  if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
-	  u = res_graph.source(e);
-	}
-	u = t;
-	while ((e = pred[u]) != INVALID) {
-	  res_graph.augment(e, delta);
-	  u = res_graph.source(e);
-	}
-	imbalance[s] -= delta;
-	imbalance[t] += delta;
+        // Augmenting along a shortest path from s to t
+        Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
+        Node u = t;
+        Edge e;
+        while ((e = pred[u]) != INVALID) {
+          Capacity rc;
+          if (u == graph.target(e)) {
+            rc = res_cap[e];
+            u = graph.source(e);
+          } else {
+            rc = flow[e];
+            u = graph.target(e);
+          }
+          if (rc < d) d = rc;
+        }
+        u = t;
+        while ((e = pred[u]) != INVALID) {
+          if (u == graph.target(e)) {
+            flow[e] += d;
+            res_cap[e] -= d;
+            u = graph.source(e);
+          } else {
+            flow[e] -= d;
+            res_cap[e] += d;
+            u = graph.target(e);
+          }
+        }
+        excess[s] -= d;
+        excess[t] += d;
       }
 
       // Handling nonzero lower bounds
       if (lower) {
-	for (EdgeIt e(graph); e != INVALID; ++e)
-	  flow[e] += (*lower)[e];
+        for (EdgeIt e(graph); e != INVALID; ++e)
+          flow[e] += (*lower)[e];
       }
       return true;
     }
-#endif
 
   }; //class CapacityScaling