# HG changeset patch
# User deba
# Date 1144082063 0
# Node ID b71f8ff62046628163ebc8e62c96a0b66c62db2c
# Parent  7bf1f64962c26a7441290c019c07487e329c6397
Edmonds-Karp MaxFlow
ResGraphAdaptor with Tolerance

diff -r 7bf1f64962c2 -r b71f8ff62046 lemon/Makefile.am
--- a/lemon/Makefile.am	Mon Apr 03 16:05:26 2006 +0000
+++ b/lemon/Makefile.am	Mon Apr 03 16:34:23 2006 +0000
@@ -35,6 +35,7 @@
 	dimacs.h \
 	dag_shortest_path.h \
 	edge_set.h \
+	edmonds_karp.h \
 	error.h \
 	eps.h \
 	fib_heap.h \
diff -r 7bf1f64962c2 -r b71f8ff62046 lemon/edmonds_karp.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/edmonds_karp.h	Mon Apr 03 16:34:23 2006 +0000
@@ -0,0 +1,390 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2006
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_EDMONDS_KARP_H
+#define LEMON_EDMONDS_KARP_H
+
+/// \file
+/// \ingroup flowalgs
+/// \brief Implementation of the Edmonds-Karp algorithm.
+
+#include <lemon/graph_adaptor.h>
+#include <lemon/tolerance.h>
+#include <lemon/bfs.h>
+
+namespace lemon {
+
+  /// \ingroup flowalgs
+  /// \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. It is possible to change these quantities using the
+  /// functions \ref source, \ref target, \ref capacityMap and \ref
+  /// flowMap.
+  ///
+  /// The time complexity of the algorithm is O(n * e^2) in worst case. 
+  /// Always try the preflow algorithm instead of this if you does not
+  /// have some additional reason than to compute the optimal flow which
+  /// has O(n^3) time complexity.
+  ///
+  /// \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.
+  ///
+  /// \author Balazs Dezso 
+  template <typename _Graph, typename _Number,
+	    typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,
+            typename _FlowMap = typename _Graph::template EdgeMap<_Number>,
+	    typename _Tolerance = Tolerance<_Number> >
+  class EdmondsKarp {
+  public:
+
+    /// \brief \ref Exception for the case when the source equals the target.
+    ///
+    /// \ref Exception for the case when the source equals the target.
+    ///
+    class InvalidArgument : public lemon::LogicError {
+    public:
+      virtual const char* exceptionName() const {
+	return "lemon::EdmondsKarp::InvalidArgument";
+      }
+    };
+
+
+    /// \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;
+    
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::EdgeIt EdgeIt;
+    typedef typename Graph::InEdgeIt InEdgeIt;
+    typedef typename Graph::OutEdgeIt OutEdgeIt;
+    
+  public:
+
+    /// \brief The constructor of the class.
+    ///
+    /// The constructor of the class. 
+    /// \param _graph The directed graph the algorithm runs on. 
+    /// \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.
+    /// Except the graph, all of these parameters can be reset by
+    /// calling \ref source, \ref target, \ref capacityMap and \ref
+    /// flowMap, resp.
+    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) 
+    {
+      if (_source == _target) {
+        throw InvalidArgument();
+      }
+    }
+
+    /// \brief Initializes the algorithm
+    /// 
+    /// It sets the flow to empty flow.
+    void init() {
+      for (EdgeIt it(_graph); it != INVALID; ++it) {
+        _flow.set(it, 0);
+      }
+      _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;
+      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _value += _flow[jt];
+      }
+      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _value -= _flow[jt];
+      }
+    }
+
+    /// \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];
+      }
+      for (NodeIt it(_graph); it != INVALID; ++it) {
+        if (it == _source || it == _target) continue;
+        Number outFlow = 0;
+        for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
+          outFlow += _flow[jt];
+        }
+        Number inFlow = 0;
+        for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
+          inFlow += _flow[jt];
+        }
+        if (_tolerance.different(outFlow, inFlow)) {
+          return false;
+        }
+      }
+      for (EdgeIt it(_graph); it != INVALID; ++it) {
+        if (_tolerance.less(_flow[it], 0)) return false;
+        if (_tolerance.less(_capacity[it], _flow[it])) return false;
+      }
+      return true;
+    }
+
+    /// \brief Augment the solution on an edge shortest path.
+    /// 
+    /// Augment the solution on an edge shortest path. It search an
+    /// edge shortest path between the source and the target
+    /// in the residual graph with the bfs algoritm.
+    /// Then it increase the flow on this path with the minimal residual
+    /// capacity on the path. If there is not such path it gives back
+    /// false.
+    /// \return %False when the augmenting is not success so the
+    /// 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);
+      
+      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;
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// It runs augmenting phases until the optimal solution is reached. 
+    void start() {
+      while (augment()) {}
+    }
+
+    /// \brief Gives back the current flow value.
+    ///
+    /// Gives back the current flow _value.
+    Number flowValue() const {
+      return _value;
+    }
+
+    /// \brief runs the algorithm.
+    /// 
+    /// It is just a shorthand for:
+    /// \code 
+    /// ek.init();
+    /// ek.start();
+    /// \endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// \brief Returns a minimum value cut.
+    ///
+    /// Sets \c cut to the characteristic vector of a minimum value cut
+    /// It simply calls the minMinCut member.
+    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.  
+    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.  
+    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 Sets the source node to \c _source.
+    ///
+    /// Sets the source node to \c _source.
+    void source(Node source) { 
+      _source = source; 
+    }
+
+    /// \brief Returns the source node.
+    ///
+    /// Returns the source node.
+    /// 
+    Node source() const { 
+      return _source;
+    }
+
+    /// \brief Sets the target node to \c target.
+    ///
+    /// Sets the target node to \c target.
+    void target(Node target) { 
+      _target = target; 
+    }
+
+    /// \brief Returns the target node.
+    ///
+    /// Returns the target node.
+    /// 
+    Node target() const { 
+      return _target;
+    }
+
+    /// \brief Sets the edge map of the capacities to _cap.
+    ///
+    /// Sets the edge map of the capacities to _cap.
+    /// 
+    void capacityMap(const CapacityMap& capacity) { 
+      _capacity = &capacity; 
+    }
+
+    /// \brief Returns a reference to capacity map.
+    ///
+    /// Returns a reference to capacity map.
+    /// 
+    const CapacityMap &capacityMap() const { 
+      return *_capacity;
+    }
+
+    /// \brief Sets the edge map of the flows to \c flow.
+    ///
+    /// Sets the edge map of the flows to \c flow.
+    /// 
+    void flowMap(FlowMap& flow) { 
+      _flow = &flow; 
+    }
+     
+    /// \brief Returns a reference to flow map.
+    ///
+    /// Returns a reference to flow map.
+    /// 
+    const FlowMap &flowMap() const { 
+      return *_flow;
+    }
+
+  private:
+    
+    const Graph& _graph;
+    const CapacityMap& _capacity;
+    FlowMap& _flow;
+    Tolerance _tolerance;
+    
+    ResGraph _resgraph;
+    Node _source, _target;
+    Number _value;
+    
+  };
+
+}
+
+#endif
diff -r 7bf1f64962c2 -r b71f8ff62046 lemon/graph_adaptor.h
--- a/lemon/graph_adaptor.h	Mon Apr 03 16:05:26 2006 +0000
+++ b/lemon/graph_adaptor.h	Mon Apr 03 16:34:23 2006 +0000
@@ -34,6 +34,8 @@
 #include <lemon/bits/graph_adaptor_extender.h>
 #include <lemon/bits/graph_extender.h>
 
+#include <lemon/tolerance.h>
+
 #include <iostream>
 
 namespace lemon {
@@ -71,6 +73,9 @@
 
   public:
     GraphAdaptorBase(Graph& _graph) : graph(&_graph) { }
+
+    Graph& getGraph() { return *graph; }
+    const Graph& getGraph() const { return *graph; }
  
     typedef typename Graph::Node Node;
     typedef typename Graph::Edge Edge;
@@ -1422,41 +1427,52 @@
     return UndirGraphAdaptor<const Graph>(graph);
   }
 
-  template<typename Graph, typename Number,
-	   typename CapacityMap, typename FlowMap>
+  template<typename Graph, typename Number,  
+           typename CapacityMap, typename FlowMap, 
+           typename Tolerance = Tolerance<Number> >
   class ResForwardFilter {
     const CapacityMap* capacity;
     const FlowMap* flow;
+    Tolerance tolerance;
   public:
     typedef typename Graph::Edge Key;
     typedef bool Value;
 
-    ResForwardFilter(const CapacityMap& _capacity, const FlowMap& _flow) 
-      : capacity(&_capacity), flow(&_flow) { }
-    ResForwardFilter() : capacity(0), flow(0) { }
+    ResForwardFilter(const CapacityMap& _capacity, const FlowMap& _flow,
+                     const Tolerance& _tolerance = Tolerance()) 
+      : capacity(&_capacity), flow(&_flow), tolerance(_tolerance) { }
+
+    ResForwardFilter(const Tolerance& _tolerance) 
+      : capacity(0), flow(0), tolerance(_tolerance)  { }
+
     void setCapacity(const CapacityMap& _capacity) { capacity = &_capacity; }
     void setFlow(const FlowMap& _flow) { flow = &_flow; }
+
     bool operator[](const typename Graph::Edge& e) const {
-      return (Number((*flow)[e]) < Number((*capacity)[e]));
+      return tolerance.less((*flow)[e], (*capacity)[e]);
     }
   };
 
   template<typename Graph, typename Number,
-	   typename CapacityMap, typename FlowMap>
+	   typename CapacityMap, typename FlowMap,
+           typename Tolerance = Tolerance<Number> >
   class ResBackwardFilter {
     const CapacityMap* capacity;
     const FlowMap* flow;
+    Tolerance tolerance;
   public:
     typedef typename Graph::Edge Key;
     typedef bool Value;
 
-    ResBackwardFilter(const CapacityMap& _capacity, const FlowMap& _flow) 
-      : capacity(&_capacity), flow(&_flow) { }
-    ResBackwardFilter() : capacity(0), flow(0) { }
+    ResBackwardFilter(const CapacityMap& _capacity, const FlowMap& _flow,
+                      const Tolerance& _tolerance = Tolerance())
+      : capacity(&_capacity), flow(&_flow), tolerance(_tolerance) { }
+    ResBackwardFilter(const Tolerance& _tolerance = Tolerance())
+      : capacity(0), flow(0), tolerance(_tolerance) { }
     void setCapacity(const CapacityMap& _capacity) { capacity = &_capacity; }
     void setFlow(const FlowMap& _flow) { flow = &_flow; }
     bool operator[](const typename Graph::Edge& e) const {
-      return (Number(0) < Number((*flow)[e]));
+      return tolerance.less(0, Number((*flow)[e]));
     }
   };
 
@@ -1497,21 +1513,22 @@
   ///\author Marton Makai
   ///
   template<typename Graph, typename Number, 
-	   typename CapacityMap, typename FlowMap>
+	   typename CapacityMap, typename FlowMap,
+           typename Tolerance = Tolerance<Number> >
   class ResGraphAdaptor : 
     public EdgeSubGraphAdaptor< 
-    UndirGraphAdaptor<Graph>, 
-    typename UndirGraphAdaptor<Graph>::template CombinedEdgeMap<
-    ResForwardFilter<Graph, Number, CapacityMap, FlowMap>,  
-    ResBackwardFilter<Graph, Number, CapacityMap, FlowMap> > > {
+    UndirGraphAdaptor<const Graph>, 
+    typename UndirGraphAdaptor<const Graph>::template CombinedEdgeMap<
+    ResForwardFilter<const Graph, Number, CapacityMap, FlowMap>,  
+    ResBackwardFilter<const Graph, Number, CapacityMap, FlowMap> > > {
   public:
 
-    typedef UndirGraphAdaptor<Graph> UGraph;
+    typedef UndirGraphAdaptor<const Graph> UGraph;
 
-    typedef ResForwardFilter<Graph, Number, CapacityMap, FlowMap> 
+    typedef ResForwardFilter<const Graph, Number, CapacityMap, FlowMap> 
     ForwardFilter;
 
-    typedef ResBackwardFilter<Graph, Number, CapacityMap, FlowMap> 
+    typedef ResBackwardFilter<const Graph, Number, CapacityMap, FlowMap> 
     BackwardFilter;
 
     typedef typename UGraph::
@@ -1544,44 +1561,77 @@
 
   public:
 
-    ResGraphAdaptor(Graph& _graph, const CapacityMap& _capacity, 
-                    FlowMap& _flow) 
-      : Parent(), capacity(&_capacity), flow(&_flow),
-        ugraph(_graph),
-        forward_filter(_capacity, _flow), 
-        backward_filter(_capacity, _flow),
-        edge_filter(forward_filter, backward_filter) {
+    const Graph& getGraph() const { return ugraph.getGraph(); }
+
+    /// \brief Constructor of the residual graph.
+    ///
+    /// Constructor of the residual graph. The parameters are the graph type,
+    /// the flow map, the capacity map and a tolerance object.
+    ResGraphAdaptor(const Graph& _graph, const CapacityMap& _capacity, 
+                    FlowMap& _flow, const Tolerance& _tolerance = Tolerance()) 
+      : Parent(), capacity(&_capacity), flow(&_flow), ugraph(_graph),
+        forward_filter(_capacity, _flow, _tolerance), 
+        backward_filter(_capacity, _flow, _tolerance),
+        edge_filter(forward_filter, backward_filter)
+    {
       Parent::setGraph(ugraph);
       Parent::setEdgeFilterMap(edge_filter);
     }
 
     typedef typename Parent::Edge Edge;
 
+    /// \brief Gives back the residual capacity of the edge.
+    ///
+    /// Gives back the residual capacity of the edge.
+    Number rescap(const Edge& edge) const {
+      if (UGraph::direction(edge)) {
+        return (*capacity)[edge]-(*flow)[edge]; 
+      } else {
+        return (*flow)[edge];
+      }
+    } 
+
+    /// \brief Augment on the given edge in the residual graph.
+    ///
+    /// Augment on the given edge in the residual graph. It increase
+    /// or decrease the flow on the original edge depend on the direction
+    /// of the residual edge.
     void augment(const Edge& e, Number a) const {
       if (UGraph::direction(e)) {
-	flow->set(e, (*flow)[e]+a);
+        flow->set(e, (*flow)[e] + a);
       } else {  
-	flow->set(e, (*flow)[e]-a);
+        flow->set(e, (*flow)[e] - a);
       }
     }
 
+    /// \brief Returns the direction of the edge.
+    ///
+    /// Returns true when the edge is same oriented as the original edge.
     static bool forward(const Edge& e) {
       return UGraph::direction(e);
     }
 
+    /// \brief Returns the direction of the edge.
+    ///
+    /// Returns true when the edge is opposite oriented as the original edge.
     static bool backward(const Edge& e) {
       return !UGraph::direction(e);
     }
 
+    /// \brief Gives back the forward oriented residual edge.
+    ///
+    /// Gives back the forward oriented residual edge.
     static Edge forward(const typename Graph::Edge& e) {
       return UGraph::direct(e, true);
     }
 
+    /// \brief Gives back the backward oriented residual edge.
+    ///
+    /// Gives back the backward oriented residual edge.
     static Edge backward(const typename Graph::Edge& e) {
       return UGraph::direct(e, false);
     }
 
-
     /// \brief Residual capacity map.
     ///
     /// In generic residual graphs the residual capacity can be obtained 
@@ -1595,11 +1645,8 @@
       ResCap(const ResGraphAdaptor& _res_graph) 
         : res_graph(&_res_graph) {}
       
-      Number operator[](const Edge& e) const { 
-	if (ResGraphAdaptor::direction(e)) 
-	  return (*(res_graph->capacity))[e]-(*(res_graph->flow))[e]; 
-	else 
-	  return (*(res_graph->flow))[e]; 
+      Number operator[](const Edge& e) const {
+        return res_graph->rescap(e);
       }
       
     };