Changeset 2514:57143c09dc20 in lemon-0.x for lemon/edmonds_karp.h

Timestamp:

11/17/07 21:58:11 (17 years ago)

Author:

Balazs Dezso

Branch:

default

Phase:

public

Convert:

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

Message:

Redesign the maximum flow algorithms

Redesigned interface
Preflow changed to use elevator
Edmonds-Karp does not use the ResGraphAdaptor?
Goldberg-Tarjan algorithm (Preflow with Dynamic Trees)
Dinitz-Sleator-Tarjan (Blocking flow with Dynamic Tree)

File:

• lemon/edmonds_karp.h (modified) (10 diffs)

lemon/edmonds_karp.h

@@ -24 +24 @@
 /// \brief Implementation of the Edmonds-Karp algorithm.
 
-#include <lemon/graph_adaptor.h>
 #include <lemon/tolerance.h>
-#include <lemon/bfs.h>
+#include <vector>
 
 namespace lemon {
 
+  /// \brief Default traits class of EdmondsKarp class.
+  ///
+  /// Default traits class of EdmondsKarp class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _CapacityMap>
+  struct EdmondsKarpDefaultTraits {
+
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \brief The type of the map that stores the edge capacities.
+    ///
+    /// The type of the map that stores the edge capacities.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _CapacityMap CapacityMap;
+
+    /// \brief The type of the length of the edges.
+    typedef typename CapacityMap::Value Value;
+
+    /// \brief The map type that stores the flow values.
+    ///
+    /// The map type that stores the flow values. 
+    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param graph The graph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Graph& graph) {
+      return new FlowMap(graph);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+
   /// \ingroup max_flow
+  ///
   /// \brief Edmonds-Karp algorithms class.
   ///
   /// This class provides an implementation of the \e Edmonds-Karp \e
   /// algorithm producing a flow of maximum value in a directed
-  /// graph. The Edmonds-Karp algorithm is slower than the Preflow algorithm
-  /// but it has an advantage of the step-by-step execution control with
-  /// feasible flow solutions. The \e source node, the \e target node, the \e
-  /// capacity of the edges and the \e starting \e flow value of the
-  /// edges should be passed to the algorithm through the
-  /// constructor.
+  /// graphs. The Edmonds-Karp algorithm is slower than the Preflow
+  /// algorithm but it has an advantage of the step-by-step execution
+  /// control with feasible flow solutions. The \e source node, the \e
+  /// target node, the \e capacity of the edges and the \e starting \e
+  /// flow value of the edges should be passed to the algorithm
+  /// through the constructor.
   ///
-  /// The time complexity of the algorithm is \f$ O(n * e^2) \f$ in
+  /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
   /// worst case.  Always try the preflow algorithm instead of this if
   /// you just want to compute the optimal flow.
   ///
   /// \param _Graph The directed graph type the algorithm runs on.
-  /// \param _Number The number type of the capacities and the flow values.
   /// \param _CapacityMap The capacity map type.
-  /// \param _FlowMap The flow map type.
-  /// \param _Tolerance The tolerance class to handle computation problems.
+  /// \param _Traits Traits class to set various data types used by
+  /// the algorithm.  The default traits class is \ref
+  /// EdmondsKarpDefaultTraits.  See \ref EdmondsKarpDefaultTraits for the
+  /// documentation of a Edmonds-Karp traits class. 
   ///
   /// \author Balazs Dezso 
 #ifdef DOXYGEN
-  template <typename _Graph, typename _Number,
-            typename _CapacityMap, typename _FlowMap, typename _Tolerance>
-#else
-  template <typename _Graph, typename _Number,
-            typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,
-            typename _FlowMap = typename _Graph::template EdgeMap<_Number>,
-            typename _Tolerance = Tolerance<_Number> >
+  template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else 
+  template <typename _Graph,
+            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+            typename _Traits = EdmondsKarpDefaultTraits<_Graph, _CapacityMap> >
 #endif
   class EdmondsKarp {
   public:
+
+    typedef _Traits Traits;
+    typedef typename Traits::Graph Graph;
+    typedef typename Traits::CapacityMap CapacityMap;
+    typedef typename Traits::Value Value; 
+
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Tolerance Tolerance;
 
     /// \brief \ref Exception for the case when the source equals the target.
@@ -77 +125 @@
 
 
-    /// \brief The graph type the algorithm runs on. 
-    typedef _Graph Graph;
-    /// \brief The value type of the algorithms.
-    typedef _Number Number;
-    /// \brief The capacity map on the edges.
-    typedef _CapacityMap CapacityMap;
-    /// \brief The flow map on the edges.
-    typedef _FlowMap FlowMap;
-    /// \brief The tolerance used by the algorithm.
-    typedef _Tolerance Tolerance;
-
-    typedef ResGraphAdaptor<Graph, Number, CapacityMap, 
-                            FlowMap, Tolerance> ResGraph;
-
   private:
 
-    typedef typename Graph::Node Node;
-    typedef typename Graph::Edge Edge;
+    GRAPH_TYPEDEFS(typename Graph);
+    typedef typename Graph::template NodeMap<Edge> PredMap;
     
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::EdgeIt EdgeIt;
-    typedef typename Graph::InEdgeIt InEdgeIt;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
+    const Graph& _graph;
+    const CapacityMap* _capacity;
+
+    Node _source, _target;
+
+    FlowMap* _flow;
+    bool _local_flow;
+
+    PredMap* _pred;
+    std::vector<Node> _queue;
+    
+    Tolerance _tolerance;
+    Value _flow_value;
+
+    void createStructures() {
+      if (!_flow) {
+        _flow = Traits::createFlowMap(_graph);
+        _local_flow = true;
+      }
+      if (!_pred) {
+        _pred = new PredMap(_graph);
+      }
+      _queue.resize(countNodes(_graph));
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+        delete _flow;
+      }
+      if (_pred) {
+        delete _pred;
+      }
+    }
     
   public:
 
+    ///\name Named template parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct DefFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Graph&) {
+        throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename _FlowMap>
+    struct DefFlowMap 
+      : public EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
+      typedef EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > 
+      Create;
+    };
+
+
+    /// @}
+
     /// \brief The constructor of the class.
     ///
     /// The constructor of the class. 
     /// \param graph The directed graph the algorithm runs on. 
+    /// \param capacity The capacity of the edges. 
     /// \param source The source node.
     /// \param target The target node.
-    /// \param capacity The capacity of the edges. 
-    /// \param flow The flow of the edges. 
-    /// \param tolerance Tolerance class.
-    EdmondsKarp(const Graph& graph, Node source, Node target,
-                const CapacityMap& capacity, FlowMap& flow, 
-                const Tolerance& tolerance = Tolerance())
-      : _graph(graph), _capacity(capacity), _flow(flow), 
-        _tolerance(tolerance), _resgraph(graph, capacity, flow, tolerance), 
-        _source(source), _target(target) 
+    EdmondsKarp(const Graph& graph, const CapacityMap& capacity,
+                Node source, Node target)
+      : _graph(graph), _capacity(&capacity), _source(source), _target(target),
+        _flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
     {
       if (_source == _target) {
@@ -123 +209 @@
       }
     }
+
+    /// \brief Destrcutor.
+    ///
+    /// Destructor.
+    ~EdmondsKarp() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the capacity map.
+    ///
+    /// Sets the capacity map.
+    /// \return \c (*this)
+    EdmondsKarp& capacityMap(const CapacityMap& map) {
+      _capacity = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// \return \c (*this)
+    EdmondsKarp& flowMap(FlowMap& map) {
+      if (_local_flow) {
+        delete _flow;
+        _local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Returns the flow map.
+    ///
+    /// \return The flow map.
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /// \brief Sets the source node.
+    ///
+    /// Sets the source node.
+    /// \return \c (*this)
+    EdmondsKarp& source(const Node& node) {
+      _source = node;
+      return *this;
+    }
+
+    /// \brief Sets the target node.
+    ///
+    /// Sets the target node.
+    /// \return \c (*this)
+    EdmondsKarp& target(const Node& node) {
+      _target = node;
+      return *this;
+    }
+
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    EdmondsKarp& tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      return *this;
+    } 
+
+    /// \brief Returns the tolerance used by algorithm.
+    ///
+    /// Returns the tolerance used by algorithm.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+
+    /// \name Execution control The simplest way to execute the
+    /// algorithm is to use the \c run() member functions.
+    /// \n
+    /// If you need more control on initial solution or
+    /// execution then you have to call one \ref init() function and then
+    /// the start() or multiple times the \c augment() member function.  
+    
+    ///@{
 
     /// \brief Initializes the algorithm
@@ -128 +292 @@
     /// It sets the flow to empty flow.
     void init() {
+      createStructures();
       for (EdgeIt it(_graph); it != INVALID; ++it) {
-        _flow.set(it, 0);
-      }
-      _value = 0;
+        _flow->set(it, 0);
+      }
+      _flow_value = 0;
     }
     
     /// \brief Initializes the algorithm
     /// 
-    /// If the flow map initially flow this let the flow map
-    /// unchanged but the flow value will be set by the flow
-    /// on the outedges from the source.
-    void flowInit() {
-      _value = 0;
+    /// Initializes the flow to the \c flowMap. The \c flowMap should
+    /// contain a feasible flow, ie. in each node excluding the source
+    /// and the target the incoming flow should be equal to the
+    /// outgoing flow.
+    template <typename FlowMap>
+    void flowInit(const FlowMap& flowMap) {
+      createStructures();
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow->set(e, flowMap[e]);
+      }
+      _flow_value = 0;
       for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value += _flow[jt];
+        _flow_value += (*_flow)[jt];
       }
       for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value -= _flow[jt];
+        _flow_value -= (*_flow)[jt];
       }
     }
@@ -151 +322 @@
     /// \brief Initializes the algorithm
     /// 
-    /// If the flow map initially flow this let the flow map
-    /// unchanged but the flow value will be set by the flow
-    /// on the outedges from the source. It also checks that
-    /// the flow map really contains a flow.
-    /// \return %True when the flow map really a flow.
-    bool checkedFlowInit() {
-      _value = 0;
-      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value += _flow[jt];
-      }
-      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value -= _flow[jt];
+    /// Initializes the flow to the \c flowMap. The \c flowMap should
+    /// contain a feasible flow, ie. in each node excluding the source
+    /// and the target the incoming flow should be equal to the
+    /// outgoing flow.  
+    /// \return %False when the given flowMap does not contain
+    /// feasible flow.
+    template <typename FlowMap>
+    bool checkedFlowInit(const FlowMap& flowMap) {
+      createStructures();
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow->set(e, flowMap[e]);
       }
       for (NodeIt it(_graph); it != INVALID; ++it) {
         if (it == _source || it == _target) continue;
-        Number outFlow = 0;
+        Value outFlow = 0;
         for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
-          outFlow += _flow[jt];
+          outFlow += (*_flow)[jt];
         }
-        Number inFlow = 0;
+        Value inFlow = 0;
         for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
-          inFlow += _flow[jt];
+          inFlow += (*_flow)[jt];
         }
         if (_tolerance.different(outFlow, inFlow)) {
@@ -179 +349 @@
       }
       for (EdgeIt it(_graph); it != INVALID; ++it) {
-        if (_tolerance.less(_flow[it], 0)) return false;
-        if (_tolerance.less(_capacity[it], _flow[it])) return false;
+        if (_tolerance.less((*_flow)[it], 0)) return false;
+        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
+      }
+      _flow_value = 0;
+      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value += (*_flow)[jt];
+      }
+      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value -= (*_flow)[jt];
       }
       return true;
@@ -196 +373 @@
     /// current flow is a feasible and optimal solution.
     bool augment() {
-      typename Bfs<ResGraph>
-      ::template DefDistMap<NullMap<Node, int> >
-      ::Create bfs(_resgraph);
-
-      NullMap<Node, int> distMap;
-      bfs.distMap(distMap);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _pred->set(n, INVALID);
+      }
       
-      bfs.init();
-      bfs.addSource(_source);
-      bfs.start(_target);
-
-      if (!bfs.reached(_target)) {
-        return false;
-      }
-      Number min = _resgraph.rescap(bfs.predEdge(_target));
-      for (Node it = bfs.predNode(_target); it != _source; 
-           it = bfs.predNode(it)) {
-        if (min > _resgraph.rescap(bfs.predEdge(it))) {
-          min = _resgraph.rescap(bfs.predEdge(it));
-        }
-      } 
-      for (Node it = _target; it != _source; it = bfs.predNode(it)) {
-        _resgraph.augment(bfs.predEdge(it), min);
-      }
-      _value += min;
-      return true;
+      int first = 0, last = 1;
+      
+      _queue[0] = _source;
+      _pred->set(_source, OutEdgeIt(_graph, _source));
+
+      while (first != last && (*_pred)[_target] == INVALID) {
+        Node n = _queue[first++];
+        
+        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+          Value rem = (*_capacity)[e] - (*_flow)[e];
+          Node t = _graph.target(e);
+          if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
+            _pred->set(t, e);
+            _queue[last++] = t;
+          }
+        }
+        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+          Value rem = (*_flow)[e];
+          Node t = _graph.source(e);
+          if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
+            _pred->set(t, e);
+            _queue[last++] = t;
+          }
+        }
+      }
+
+      if ((*_pred)[_target] != INVALID) {
+        Node n = _target;
+        Edge e = (*_pred)[n];
+
+        Value prem = (*_capacity)[e] - (*_flow)[e];
+        n = _graph.source(e);
+        while (n != _source) {
+          e = (*_pred)[n];
+          if (_graph.target(e) == n) {
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (rem < prem) prem = rem;
+            n = _graph.source(e);
+          } else {
+            Value rem = (*_flow)[e];
+            if (rem < prem) prem = rem;
+            n = _graph.target(e);   
+          } 
+        }
+
+        n = _target;
+        e = (*_pred)[n];
+
+        _flow->set(e, (*_flow)[e] + prem);
+        n = _graph.source(e);
+        while (n != _source) {
+          e = (*_pred)[n];
+          if (_graph.target(e) == n) {
+            _flow->set(e, (*_flow)[e] + prem);
+            n = _graph.source(e);
+          } else {
+            _flow->set(e, (*_flow)[e] - prem);
+            n = _graph.target(e);   
+          } 
+        }
+
+        _flow_value += prem;    
+        return true;
+      } else {
+        return false;
+      }
     }
 
@@ -229 +450 @@
     void start() {
       while (augment()) {}
-    }
-
-    /// \brief Gives back the current flow value.
-    ///
-    /// Gives back the current flow _value.
-    Number flowValue() const {
-      return _value;
     }
 
@@ -251 +465 @@
     }
 
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Dijkstra algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Returns the value of the maximum flow.
+    ///
+    /// Returns the value of the maximum flow by returning the excess
+    /// of the target node \c t. This value equals to the value of
+    /// the maximum flow already after the first phase.
+    Value flowValue() const {
+      return _flow_value;
+    }
+
+
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Sets the \c flowMap to the flow on the edges. This method can
+    /// be called after the second phase of algorithm.
+    Value flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+
+    /// \brief Returns true when the node is on the source side of minimum cut.
+    ///
+
+    /// Returns true when the node is on the source side of minimum
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase().
+    bool minCut(const Node& node) const {
+      return (*_pred)[node] != INVALID;
+    }
+
     /// \brief Returns a minimum value cut.
     ///
@@ -257 +509 @@
     /// \retval cut Write node bool map. 
     template <typename CutMap>
-    void minCut(CutMap& cut) const {
-      minMinCut(cut);
-    }
-
-    /// \brief Returns the inclusionwise minimum of the minimum value cuts.
-    ///
-    /// Sets \c cut to the characteristic vector of the minimum value cut
-    /// which is inclusionwise minimum. It is computed by processing a
-    /// bfs from the source node \c source in the residual graph.  
-    /// \retval cut Write node bool map. 
-    template <typename CutMap>
-    void minMinCut(CutMap& cut) const {
-
-      typename Bfs<ResGraph>
-      ::template DefDistMap<NullMap<Node, int> >
-      ::template DefProcessedMap<CutMap>
-      ::Create bfs(_resgraph);
-
-      NullMap<Node, int> distMap;
-      bfs.distMap(distMap);
-
-      bfs.processedMap(cut);
-     
-      bfs.run(_source);
-    }
-
-    /// \brief Returns the inclusionwise minimum of the minimum value cuts.
-    ///
-    /// Sets \c cut to the characteristic vector of the minimum value cut
-    /// which is inclusionwise minimum. It is computed by processing a
-    /// bfs from the source node \c source in the residual graph.  
-    /// \retval cut Write node bool map. 
-    template <typename CutMap>
-    void maxMinCut(CutMap& cut) const {
-
-      typedef RevGraphAdaptor<const ResGraph> RevGraph;
-
-      RevGraph revgraph(_resgraph);
-
-      typename Bfs<RevGraph>
-      ::template DefDistMap<NullMap<Node, int> >
-      ::template DefPredMap<NullMap<Node, Edge> >
-      ::template DefProcessedMap<NotWriteMap<CutMap> >
-      ::Create bfs(revgraph);
-
-      NullMap<Node, int> distMap;
-      bfs.distMap(distMap);
-
-      NullMap<Node, Edge> predMap;
-      bfs.predMap(predMap);
-
-      NotWriteMap<CutMap> notcut(cut);
-      bfs.processedMap(notcut);
-     
-      bfs.run(_target);
-    }
-
-    /// \brief Returns the source node.
-    ///
-    /// Returns the source node.
-    /// 
-    Node source() const { 
-      return _source;
-    }
-
-    /// \brief Returns the target node.
-    ///
-    /// Returns the target node.
-    /// 
-    Node target() const { 
-      return _target;
-    }
-
-    /// \brief Returns a reference to capacity map.
-    ///
-    /// Returns a reference to capacity map.
-    /// 
-    const CapacityMap &capacityMap() const { 
-      return *_capacity;
-    }
-     
-    /// \brief Returns a reference to flow map.
-    ///
-    /// Returns a reference to flow map.
-    /// 
-    const FlowMap &flowMap() const { 
-      return *_flow;
-    }
-
-    /// \brief Returns the tolerance used by algorithm.
-    ///
-    /// Returns the tolerance used by algorithm.
-    const Tolerance& tolerance() const {
-      return tolerance;
-    } 
-    
-  private:
-    
-    const Graph& _graph;
-    const CapacityMap& _capacity;
-    FlowMap& _flow;
-    Tolerance _tolerance;
-    
-    ResGraph _resgraph;
-    Node _source, _target;
-    Number _value;
-    
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        cutMap.set(n, (*_pred)[n] != INVALID);
+      }
+      cutMap.set(_source, true);
+    }    
+
+    /// @}
+
   };