Changeset 2581:054566ac0934 in lemon-0.x for lemon/cost_scaling.h

Timestamp:

02/28/08 03:54:27 (16 years ago)

Author:

Peter Kovacs

Branch:

default

Phase:

public

Convert:

svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3463

Message:

Query improvements in the min cost flow algorithms.

• External flow and potential maps can be used.
• New query functions: flow() and potential().
• CycleCanceling? also provides dual solution (node potentials).
• Doc improvements.

File:

• lemon/cost_scaling.h (modified) (33 diffs)

lemon/cost_scaling.h

@@ -55 +55 @@
   /// - 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
@@ -98 +95 @@
 
     /// 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;
@@ -108 +105 @@
     /// \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];
       }
@@ -161 +159 @@
 
     // 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:
@@ -181 +179 @@
 
     // 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;
@@ -198 +198 @@
   public:
 
-    /// \brief General constructor of the class (with lower bounds).
-    ///
-    /// General constructor of the class (with lower bounds).
+    /// \brief General constructor (with lower bounds).
+    ///
+    /// General constructor (with lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
@@ -213 +213 @@
                  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
@@ -232 +232 @@
     }
 
-    /// \brief General constructor of the class (without lower bounds).
-    ///
-    /// General constructor of the class (without lower bounds).
+    /// \brief General constructor (without lower bounds).
+    ///
+    /// General constructor (without lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
@@ -245 +245 @@
                  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
@@ -255 +255 @@
     }
 
-    /// \brief Simple constructor of the class (with lower bounds).
-    ///
-    /// Simple constructor of the class (with lower bounds).
+    /// \brief Simple constructor (with lower bounds).
+    ///
+    /// Simple constructor (with lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
@@ -274 +274 @@
                  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
@@ -293 +293 @@
     }
 
-    /// \brief Simple constructor of the class (without lower bounds).
-    ///
-    /// Simple constructor of the class (without lower bounds).
+    /// \brief Simple constructor (without lower bounds).
+    ///
+    /// Simple constructor (without lower bounds).
     ///
     /// \param graph The directed graph the algorithm runs on.
@@ -310 +310 @@
                  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;
@@ -319 +319 @@
     }
 
+    /// 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.
     ///
@@ -325 +367 @@
     /// \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
@@ -335 +386 @@
     /// \pre \ref run() must be called before using this function.
     const FlowMap& flowMap() const {
-      return _flow;
+      return *_flow;
     }
 
@@ -346 +397 @@
     /// \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];
     }
 
@@ -358 +427 @@
       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:
@@ -367 +438 @@
     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
@@ -380 +464 @@
       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();
     }
 
@@ -398 +482 @@
         // 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;
@@ -408 +492 @@
           if (done) {
             for (NodeIt n(_graph); n != INVALID; ++n)
-              _potential[n] = bf.dist(n);
+              (*_potential)[n] = bf.dist(n);
             continue;
           }
@@ -416 +500 @@
         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;
           }
         }
@@ -441 +525 @@
           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);
 
@@ -449 +533 @@
                 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;
@@ -460 +544 @@
                 // Pushing flow along the edge
                 if (ahead < delta) {
-                  _flow[e] += ahead;
+                  (*_flow)[e] += ahead;
                   _excess[n] -= ahead;
                   _excess[t] += ahead;
@@ -468 +552 @@
                   break;
                 } else {
-                  _flow[e] += delta;
+                  (*_flow)[e] += delta;
                   _excess[n] -= delta;
                   _excess[t] += delta;
@@ -482 +566 @@
           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);
 
@@ -489 +573 @@
                 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;
@@ -500 +584 @@
                 // Pushing flow along the edge
                 if (ahead < delta) {
-                  _flow[e] -= ahead;
+                  (*_flow)[e] -= ahead;
                   _excess[n] -= ahead;
                   _excess[t] += ahead;
@@ -508 +592 @@
                   break;
                 } else {
-                  _flow[e] -= delta;
+                  (*_flow)[e] -= delta;
                   _excess[n] -= delta;
                   _excess[t] += delta;
@@ -524 +608 @@
             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;
           }
@@ -545 +629 @@
       }
 
+      // 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;