diff -r fa71d9612c42 -r 054566ac0934 lemon/cost_scaling.h
--- a/lemon/cost_scaling.h	Wed Feb 27 11:39:03 2008 +0000
+++ b/lemon/cost_scaling.h	Thu Feb 28 02:54:27 2008 +0000
@@ -54,11 +54,8 @@
   /// \warning
   /// - Edge capacities and costs should be \e non-negative \e integers.
   /// - Supply values should be \e signed \e integers.
-  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
-  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
-  ///   convertible to each other.
-  /// - All value types must be convertible to \c CostMap::Value, which
-  ///   must be signed type.
+  /// - The value types of the maps should be convertible to each other.
+  /// - \c CostMap::Value must be signed type.
   ///
   /// \note Edge costs are multiplied with the number of nodes during
   /// the algorithm so overflow problems may arise more easily than with
@@ -97,7 +94,7 @@
   public:
 
     /// The type of the flow map.
-    typedef CapacityEdgeMap FlowMap;
+    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
     /// The type of the potential map.
     typedef typename Graph::template NodeMap<LCost> PotentialMap;
 
@@ -107,20 +104,21 @@
     ///
     /// \ref ResidualCostMap is a map adaptor class for handling
     /// residual edge costs.
-    class ResidualCostMap : public MapBase<ResEdge, LCost>
+    template <typename Map>
+    class ResidualCostMap : public MapBase<ResEdge, typename Map::Value>
     {
     private:
 
-      const LargeCostMap &_cost_map;
+      const Map &_cost_map;
 
     public:
 
       ///\e
-      ResidualCostMap(const LargeCostMap &cost_map) :
+      ResidualCostMap(const Map &cost_map) :
         _cost_map(cost_map) {}
 
       ///\e
-      LCost operator[](const ResEdge &e) const {
+      typename Map::Value operator[](const ResEdge &e) const {
         return ResGraph::forward(e) ?  _cost_map[e] : -_cost_map[e];
       }
 
@@ -160,8 +158,8 @@
     static const int ALPHA = 4;
 
     // Paramters for heuristics
-    static const int BF_HEURISTIC_EPSILON_BOUND    = 5000;
-    static const int BF_HEURISTIC_BOUND_FACTOR = 3;
+    static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
+    static const int BF_HEURISTIC_BOUND_FACTOR  = 3;
 
   private:
 
@@ -180,16 +178,18 @@
     bool _valid_supply;
 
     // Edge map of the current flow
-    FlowMap _flow;
+    FlowMap *_flow;
+    bool _local_flow;
     // Node map of the current potentials
-    PotentialMap _potential;
+    PotentialMap *_potential;
+    bool _local_potential;
 
     // The residual graph
-    ResGraph _res_graph;
+    ResGraph *_res_graph;
     // The residual cost map
-    ResidualCostMap _res_cost;
+    ResidualCostMap<LargeCostMap> _res_cost;
     // The reduced cost map
-    ReducedCostMap _red_cost;
+    ReducedCostMap *_red_cost;
     // The excess map
     SupplyNodeMap _excess;
     // The epsilon parameter used for cost scaling
@@ -197,9 +197,9 @@
 
   public:
 
-    /// \brief General constructor of the class (with lower bounds).
+    /// \brief General constructor (with lower bounds).
     ///
-    /// General constructor of the class (with lower bounds).
+    /// General constructor (with lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
     /// \param lower The lower bounds of the edges.
@@ -212,9 +212,9 @@
                  const CostMap &cost,
                  const SupplyMap &supply ) :
       _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
-      _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
-      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
-      _red_cost(graph, _cost, _potential), _excess(graph, 0)
+      _cost(graph), _supply(graph), _flow(0), _local_flow(false),
+      _potential(0), _local_potential(false), _res_cost(_cost),
+      _excess(graph, 0)
     {
       // Removing non-zero lower bounds
       _capacity = subMap(capacity, lower);
@@ -231,9 +231,9 @@
       _valid_supply = sum == 0;
     }
 
-    /// \brief General constructor of the class (without lower bounds).
+    /// \brief General constructor (without lower bounds).
     ///
-    /// General constructor of the class (without lower bounds).
+    /// General constructor (without lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
     /// \param capacity The capacities (upper bounds) of the edges.
@@ -244,9 +244,9 @@
                  const CostMap &cost,
                  const SupplyMap &supply ) :
       _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
-      _cost(graph), _supply(supply), _flow(graph, 0), _potential(graph, 0),
-      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
-      _red_cost(graph, _cost, _potential), _excess(graph, 0)
+      _cost(graph), _supply(supply), _flow(0), _local_flow(false),
+      _potential(0), _local_potential(false), _res_cost(_cost),
+      _excess(graph, 0)
     {
       // Checking the sum of supply values
       Supply sum = 0;
@@ -254,9 +254,9 @@
       _valid_supply = sum == 0;
     }
 
-    /// \brief Simple constructor of the class (with lower bounds).
+    /// \brief Simple constructor (with lower bounds).
     ///
-    /// Simple constructor of the class (with lower bounds).
+    /// Simple constructor (with lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
     /// \param lower The lower bounds of the edges.
@@ -273,9 +273,9 @@
                  Node s, Node t,
                  Supply flow_value ) :
       _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
-      _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
-      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
-      _red_cost(graph, _cost, _potential), _excess(graph, 0)
+      _cost(graph), _supply(graph), _flow(0), _local_flow(false),
+      _potential(0), _local_potential(false), _res_cost(_cost),
+      _excess(graph, 0)
     {
       // Removing nonzero lower bounds
       _capacity = subMap(capacity, lower);
@@ -292,9 +292,9 @@
       _valid_supply = true;
     }
 
-    /// \brief Simple constructor of the class (without lower bounds).
+    /// \brief Simple constructor (without lower bounds).
     ///
-    /// Simple constructor of the class (without lower bounds).
+    /// Simple constructor (without lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
     /// \param capacity The capacities (upper bounds) of the edges.
@@ -309,24 +309,75 @@
                  Node s, Node t,
                  Supply flow_value ) :
       _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
-      _cost(graph), _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
-      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
-      _red_cost(graph, _cost, _potential), _excess(graph, 0)
+      _cost(graph), _supply(graph, 0), _flow(0), _local_flow(false),
+      _potential(0), _local_potential(false), _res_cost(_cost),
+      _excess(graph, 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 Sets the flow map.
+    ///
+    /// Sets the flow map.
+    ///
+    /// \return \c (*this)
+    CostScaling& flowMap(FlowMap &map) {
+      if (_local_flow) {
+        delete _flow;
+        _local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Sets the potential map.
+    ///
+    /// Sets 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 only way to execute the algorithm is to call the run()
+    /// function.
+
+    /// @{
+
     /// \brief Runs the algorithm.
     ///
     /// Runs the algorithm.
     ///
     /// \return \c true if a feasible flow can be found.
     bool run() {
-      init() && start();
+      return init() && start();
     }
 
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the algorithm can be obtained using these
+    /// functions.
+    /// \n run() must be called before using them.
+
+    /// @{
+
     /// \brief Returns a const reference to the edge map storing the
     /// found flow.
     ///
@@ -334,7 +385,7 @@
     ///
     /// \pre \ref run() must be called before using this function.
     const FlowMap& flowMap() const {
-      return _flow;
+      return *_flow;
     }
 
     /// \brief Returns a const reference to the node map storing the
@@ -345,7 +396,25 @@
     ///
     /// \pre \ref run() must be called before using this function.
     const PotentialMap& potentialMap() const {
-      return _potential;
+      return *_potential;
+    }
+
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Returns the flow on the edge.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Capacity flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+
+    /// \brief Returns the potential of the node.
+    ///
+    /// Returns the potential of the node.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost potential(const Node& node) const {
+      return (*_potential)[node];
     }
 
     /// \brief Returns the total cost of the found flow.
@@ -357,16 +426,31 @@
     Cost totalCost() const {
       Cost c = 0;
       for (EdgeIt e(_graph); e != INVALID; ++e)
-        c += _flow[e] * _orig_cost[e];
+        c += (*_flow)[e] * _orig_cost[e];
       return c;
     }
 
+    /// @}
+
   private:
 
     /// Initializes the algorithm.
     bool init() {
       if (!_valid_supply) return false;
 
+      // Initializing 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 ResGraph(_graph, _capacity, *_flow);
+
       // Initializing the scaled cost map and the epsilon parameter
       Cost max_cost = 0;
       int node_num = countNodes(_graph);
@@ -379,9 +463,9 @@
       // Finding a feasible flow using Circulation
       Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
                    SupplyMap >
-        circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
+        circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
                      _supply );
-      return circulation.flowMap(_flow).run();
+      return circulation.flowMap(*_flow).run();
     }
 
 
@@ -397,9 +481,9 @@
         // Performing price refinement heuristic using Bellman-Ford
         // algorithm
         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
-          typedef ShiftMap<ResidualCostMap> ShiftCostMap;
+          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
           ShiftCostMap shift_cost(_res_cost, _epsilon);
-          BellmanFord<ResGraph, ShiftCostMap> bf(_res_graph, shift_cost);
+          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
           bf.init(0);
           bool done = false;
           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
@@ -407,7 +491,7 @@
             done = bf.processNextWeakRound();
           if (done) {
             for (NodeIt n(_graph); n != INVALID; ++n)
-              _potential[n] = bf.dist(n);
+              (*_potential)[n] = bf.dist(n);
             continue;
           }
         }
@@ -415,16 +499,16 @@
         // Saturating edges not satisfying the optimality condition
         Capacity delta;
         for (EdgeIt e(_graph); e != INVALID; ++e) {
-          if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
-            delta = _capacity[e] - _flow[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];
+            (*_flow)[e] = _capacity[e];
           }
-          if (_flow[e] > 0 && -_red_cost[e] < 0) {
-            _excess[_graph.target(e)] -= _flow[e];
-            _excess[_graph.source(e)] += _flow[e];
-            _flow[e] = 0;
+          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+            _excess[_graph.target(e)] -= (*_flow)[e];
+            _excess[_graph.source(e)] += (*_flow)[e];
+            (*_flow)[e] = 0;
           }
         }
 
@@ -440,26 +524,26 @@
           // Performing push operations if there are admissible edges
           if (_excess[n] > 0) {
             for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
-              if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
-                delta = _capacity[e] - _flow[e] <= _excess[n] ?
-                        _capacity[e] - _flow[e] : _excess[n];
+              if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
+                delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
+                        _capacity[e] - (*_flow)[e] : _excess[n];
                 t = _graph.target(e);
 
                 // Push-look-ahead heuristic
                 Capacity ahead = -_excess[t];
                 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
-                  if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
-                    ahead += _capacity[oe] - _flow[oe];
+                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
+                    ahead += _capacity[oe] - (*_flow)[oe];
                 }
                 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
-                  if (_flow[ie] > 0 && -_red_cost[ie] < 0)
-                    ahead += _flow[ie];
+                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
+                    ahead += (*_flow)[ie];
                 }
                 if (ahead < 0) ahead = 0;
 
                 // Pushing flow along the edge
                 if (ahead < delta) {
-                  _flow[e] += ahead;
+                  (*_flow)[e] += ahead;
                   _excess[n] -= ahead;
                   _excess[t] += ahead;
                   active_nodes.push_front(t);
@@ -467,7 +551,7 @@
                   relabel_enabled = false;
                   break;
                 } else {
-                  _flow[e] += delta;
+                  (*_flow)[e] += delta;
                   _excess[n] -= delta;
                   _excess[t] += delta;
                   if (_excess[t] > 0 && _excess[t] <= delta)
@@ -481,25 +565,25 @@
 
           if (_excess[n] > 0) {
             for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
-              if (_flow[e] > 0 && -_red_cost[e] < 0) {
-                delta = _flow[e] <= _excess[n] ? _flow[e] : _excess[n];
+              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+                delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
                 t = _graph.source(e);
 
                 // Push-look-ahead heuristic
                 Capacity ahead = -_excess[t];
                 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
-                  if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
-                    ahead += _capacity[oe] - _flow[oe];
+                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
+                    ahead += _capacity[oe] - (*_flow)[oe];
                 }
                 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
-                  if (_flow[ie] > 0 && -_red_cost[ie] < 0)
-                    ahead += _flow[ie];
+                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
+                    ahead += (*_flow)[ie];
                 }
                 if (ahead < 0) ahead = 0;
 
                 // Pushing flow along the edge
                 if (ahead < delta) {
-                  _flow[e] -= ahead;
+                  (*_flow)[e] -= ahead;
                   _excess[n] -= ahead;
                   _excess[t] += ahead;
                   active_nodes.push_front(t);
@@ -507,7 +591,7 @@
                   relabel_enabled = false;
                   break;
                 } else {
-                  _flow[e] -= delta;
+                  (*_flow)[e] -= delta;
                   _excess[n] -= delta;
                   _excess[t] += delta;
                   if (_excess[t] > 0 && _excess[t] <= delta)
@@ -523,15 +607,15 @@
             // Performing relabel operation if the node is still active
             LCost min_red_cost = std::numeric_limits<LCost>::max();
             for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
-              if ( _capacity[oe] - _flow[oe] > 0 &&
-                   _red_cost[oe] < min_red_cost )
-                min_red_cost = _red_cost[oe];
+              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
+                   (*_red_cost)[oe] < min_red_cost )
+                min_red_cost = (*_red_cost)[oe];
             }
             for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
-              if (_flow[ie] > 0 && -_red_cost[ie] < min_red_cost)
-                min_red_cost = -_red_cost[ie];
+              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
+                min_red_cost = -(*_red_cost)[ie];
             }
-            _potential[n] -= min_red_cost + _epsilon;
+            (*_potential)[n] -= min_red_cost + _epsilon;
             hyper[n] = false;
           }
 
@@ -544,10 +628,18 @@
         }
       }
 
+      // Computing node potentials for the original costs
+      ResidualCostMap<CostMap> res_cost(_orig_cost);
+      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
+        bf(*_res_graph, res_cost);
+      bf.init(0); bf.start();
+      for (NodeIt n(_graph); n != INVALID; ++n)
+        (*_potential)[n] = bf.dist(n);
+
       // Handling non-zero lower bounds
       if (_lower) {
         for (EdgeIt e(_graph); e != INVALID; ++e)
-          _flow[e] += (*_lower)[e];
+          (*_flow)[e] += (*_lower)[e];
       }
       return true;
     }