Index: lemon/edmonds_karp.h
===================================================================
--- lemon/edmonds_karp.h	(revision 1224)
+++ lemon/edmonds_karp.h	(revision 1224)
@@ -0,0 +1,515 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2010
+ * 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 max_flow
+/// \brief Implementation of the Edmonds-Karp algorithm.
+
+#include <lemon/tolerance.h>
+#include <vector>
+
+namespace lemon {
+
+  /// \brief Default traits class of EdmondsKarp class.
+  ///
+  /// Default traits class of EdmondsKarp class.
+  /// \param GR Digraph type.
+  /// \param CAP Type of capacity map.
+  template <typename GR, typename CAP>
+  struct EdmondsKarpDefaultTraits {
+
+    /// \brief The digraph type the algorithm runs on. 
+    typedef GR Digraph;
+
+    /// \brief The type of the map that stores the arc capacities.
+    ///
+    /// The type of the map that stores the arc capacities.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef CAP CapacityMap;
+
+    /// \brief The type of the length of the arcs.
+    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 Digraph::template ArcMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param digraph The digraph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Digraph& digraph) {
+      return new FlowMap(digraph);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef lemon::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 directed
+  /// digraphs. 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 arcs and the \e starting \e
+  /// flow value of the arcs should be passed to the algorithm
+  /// through the constructor.
+  ///
+  /// 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 GR The digraph type the algorithm runs on.
+  /// \param CAP The capacity map type.
+  /// \param TR 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. 
+
+#ifdef DOXYGEN
+  template <typename GR, typename CAP, typename TR>
+#else 
+  template <typename GR,
+	    typename CAP = typename GR::template ArcMap<int>,
+            typename TR = EdmondsKarpDefaultTraits<GR, CAP> >
+#endif
+  class EdmondsKarp {
+  public:
+
+    typedef TR Traits;
+    typedef typename Traits::Digraph Digraph;
+    typedef typename Traits::CapacityMap CapacityMap;
+    typedef typename Traits::Value Value; 
+
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Tolerance Tolerance;
+
+  private:
+
+    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+    typedef typename Digraph::template NodeMap<Arc> PredMap;
+    
+    const Digraph& _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 T>
+    struct DefFlowMapTraits : public Traits {
+      typedef T FlowMap;
+      static FlowMap *createFlowMap(const Digraph&) {
+	LEMON_ASSERT(false,"Uninitialized parameter.");
+        return 0;
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename T>
+    struct DefFlowMap 
+      : public EdmondsKarp<Digraph, CapacityMap, DefFlowMapTraits<T> > {
+      typedef EdmondsKarp<Digraph, CapacityMap, DefFlowMapTraits<T> > 
+      Create;
+    };
+
+
+    /// @}
+
+  protected:
+    
+    EdmondsKarp() {}
+
+  public:
+
+    /// \brief The constructor of the class.
+    ///
+    /// The constructor of the class. 
+    /// \param digraph The digraph the algorithm runs on. 
+    /// \param capacity The capacity of the arcs. 
+    /// \param source The source node.
+    /// \param target The target node.
+    EdmondsKarp(const Digraph& digraph, const CapacityMap& capacity,
+		Node source, Node target)
+      : _graph(digraph), _capacity(&capacity), _source(source), _target(target),
+	_flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
+    {
+      LEMON_ASSERT(_source != _target,"Flow source and target are the same nodes.");
+    }
+
+    /// \brief Destructor.
+    ///
+    /// 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() const {
+      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) {
+      _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
+    /// 
+    /// Sets the flow to empty flow.
+    void init() {
+      createStructures();
+      for (ArcIt it(_graph); it != INVALID; ++it) {
+        _flow->set(it, 0);
+      }
+      _flow_value = 0;
+    }
+    
+    /// \brief Initializes the algorithm
+    /// 
+    /// 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 (ArcIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
+      }
+      _flow_value = 0;
+      for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value += (*_flow)[jt];
+      }
+      for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value -= (*_flow)[jt];
+      }
+    }
+
+    /// \brief Initializes the algorithm
+    /// 
+    /// 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 (ArcIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
+      }
+      for (NodeIt it(_graph); it != INVALID; ++it) {
+        if (it == _source || it == _target) continue;
+        Value outFlow = 0;
+        for (OutArcIt jt(_graph, it); jt != INVALID; ++jt) {
+          outFlow += (*_flow)[jt];
+        }
+        Value inFlow = 0;
+        for (InArcIt jt(_graph, it); jt != INVALID; ++jt) {
+          inFlow += (*_flow)[jt];
+        }
+        if (_tolerance.different(outFlow, inFlow)) {
+          return false;
+        }
+      }
+      for (ArcIt it(_graph); it != INVALID; ++it) {
+        if (_tolerance.less((*_flow)[it], 0)) return false;
+        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
+      }
+      _flow_value = 0;
+      for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value += (*_flow)[jt];
+      }
+      for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value -= (*_flow)[jt];
+      }
+      return true;
+    }
+
+    /// \brief Augment the solution on an arc shortest path.
+    /// 
+    /// Augment the solution on an arc shortest path. It searches an
+    /// arc shortest path between the source and the target
+    /// in the residual digraph by the bfs algoritm.
+    /// Then it increases the flow on this path with the minimal residual
+    /// capacity on the path. If there is no such path it gives back
+    /// false.
+    /// \return %False when the augmenting didn't success so the
+    /// current flow is a feasible and optimal solution.
+    bool augment() {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	_pred->set(n, INVALID);
+      }
+      
+      int first = 0, last = 1;
+      
+      _queue[0] = _source;
+      _pred->set(_source, OutArcIt(_graph, _source));
+
+      while (first != last && (*_pred)[_target] == INVALID) {
+	Node n = _queue[first++];
+	
+	for (OutArcIt 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 (InArcIt 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;
+	Arc 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;
+      }
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// It runs augmenting phases until the optimal solution is reached. 
+    void start() {
+      while (augment()) {}
+    }
+
+    /// \brief Runs the algorithm.
+    /// 
+    /// It is just a shorthand for:
+    ///
+    ///\code 
+    /// ek.init();
+    /// ek.start();
+    ///\endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the Edmonds-Karp 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.
+
+    Value flowValue() const {
+      return _flow_value;
+    }
+
+
+    /// \brief Returns the flow on the arc.
+    ///
+    /// Sets the \c flowMap to the flow on the arcs.
+    Value flow(const Arc& arc) const {
+      return (*_flow)[arc];
+    }
+
+    /// \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.
+
+    bool minCut(const Node& node) const {
+      return ((*_pred)[node] != INVALID) or node == _source;
+    }
+
+    /// \brief Returns a minimum value cut.
+    ///
+    /// Sets \c cutMap to the characteristic vector of a minimum value cut.
+
+    template <typename CutMap>
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	cutMap.set(n, (*_pred)[n] != INVALID);
+      }
+      cutMap.set(_source, true);
+    }    
+
+    /// @}
+
+  };
+
+}
+
+#endif
Index: test/CMakeLists.txt
===================================================================
--- test/CMakeLists.txt	(revision 1206)
+++ test/CMakeLists.txt	(revision 1224)
@@ -26,4 +26,5 @@
   dim_test
   edge_set_test
+  edmonds_karp_test
   error_test
   euler_test
Index: test/edmonds_karp_test.cc
===================================================================
--- test/edmonds_karp_test.cc	(revision 1224)
+++ test/edmonds_karp_test.cc	(revision 1224)
@@ -0,0 +1,236 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2010
+ * 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.
+ *
+ */
+
+#include<iostream>
+
+#include "test_tools.h"
+#include<lemon/smart_graph.h>
+#include<lemon/edmonds_karp.h>
+#include <lemon/concepts/digraph.h>
+#include <lemon/concepts/maps.h>
+#include <lemon/lgf_reader.h>
+
+using namespace lemon;
+
+char test_lgf[] =
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "8\n"
+  "9\n"
+  "@arcs\n"
+  "    label capacity\n"
+  "0 1 0     20\n"
+  "0 2 1     0\n"
+  "1 1 2     3\n"
+  "1 2 3     8\n"
+  "1 3 4     8\n"
+  "2 5 5     5\n"
+  "3 2 6     5\n"
+  "3 5 7     5\n"
+  "3 6 8     5\n"
+  "4 3 9     3\n"
+  "5 7 10    3\n"
+  "5 6 11    10\n"
+  "5 8 12    10\n"
+  "6 8 13    8\n"
+  "8 9 14    20\n"
+  "8 1 15    5\n"
+  "9 5 16    5\n"
+  "@attributes\n"
+  "source 1\n"
+  "target 8\n";
+
+void checkEdmondKarpCompile() {
+  typedef int VType;
+  typedef concepts::Digraph Digraph;
+
+  typedef Digraph::Node Node;
+  typedef Digraph::Arc Arc;
+  typedef concepts::ReadMap<Arc,VType> CapMap;
+  typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
+  typedef concepts::WriteMap<Node,bool> CutMap;
+
+  Digraph g;
+  Node n;
+  Arc e;
+  CapMap cap;
+  FlowMap flow;
+  CutMap cut;
+  VType v;
+  bool b;
+  ignore_unused_variable_warning(v,b);
+  typedef EdmondsKarp<Digraph, CapMap>
+              ::DefFlowMap<FlowMap>
+              ::Create EKType;
+  EKType ek_test(g, cap, n, n);
+  const EKType& const_ek_test = ek_test;
+
+  EKType::Tolerance tol = const_ek_test.tolerance();
+  ek_test.tolerance(tol);
+
+  ek_test
+    .capacityMap(cap)
+    .flowMap(flow)
+    .source(n)
+    .target(n);
+
+  ek_test.init();
+  ek_test.start();
+
+  v = const_ek_test.flowValue();
+  v = const_ek_test.flow(e);
+
+  const FlowMap& fm = const_ek_test.flowMap();
+  b = const_ek_test.minCut(n);
+  const_ek_test.minCutMap(cut);
+
+  ignore_unused_variable_warning(fm);
+}
+
+int cutValue (const SmartDigraph& g,
+              const SmartDigraph::NodeMap<bool>& cut,
+              const SmartDigraph::ArcMap<int>& cap) {
+
+  int c=0;
+  for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {
+    if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];
+  }
+  return c;
+}
+
+bool checkFlow(const SmartDigraph& g,
+               const SmartDigraph::ArcMap<int>& flow,
+               const SmartDigraph::ArcMap<int>& cap,
+               SmartDigraph::Node s, SmartDigraph::Node t) {
+
+  for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
+    if (flow[e] < 0 || flow[e] > cap[e]) return false;
+  }
+
+  for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) {
+    if (n == s || n == t) continue;
+    int sum = 0;
+    for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) {
+      sum += flow[e];
+    }
+    for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) {
+      sum -= flow[e];
+    }
+    if (sum != 0) return false;
+  }
+  return true;
+}
+
+int main() {
+
+  typedef SmartDigraph Digraph;
+
+  typedef Digraph::Node Node;
+  typedef Digraph::NodeIt NodeIt;
+  typedef Digraph::ArcIt ArcIt;
+  typedef Digraph::ArcMap<int> CapMap;
+  typedef Digraph::ArcMap<int> FlowMap;
+  typedef Digraph::NodeMap<bool> CutMap;
+
+  typedef EdmondsKarp<Digraph, CapMap> EKType;
+
+  Digraph g;
+  Node s, t;
+  CapMap cap(g);
+  std::istringstream input(test_lgf);
+  DigraphReader<Digraph>(g,input).
+    arcMap("capacity", cap).
+    node("source",s).
+    node("target",t).
+    run();
+
+  EKType ek_test(g, cap, s, t);
+  ek_test.run();
+
+  check(checkFlow(g, ek_test.flowMap(), cap, s, t),
+        "The flow is not feasible.");
+
+  CutMap min_cut(g);
+  ek_test.minCutMap(min_cut);
+  int min_cut_value=cutValue(g,min_cut,cap);
+
+  check(ek_test.flowValue() == min_cut_value,
+        "The max flow value is not equal to the three min cut values.");
+
+  FlowMap flow(g);
+  for(ArcIt e(g); e!=INVALID; ++e) flow[e] = ek_test.flowMap()[e];
+
+  int flow_value=ek_test.flowValue();
+
+  for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
+  ek_test.flowInit(flow);
+  ek_test.start();
+
+  CutMap min_cut1(g);
+  ek_test.minCutMap(min_cut1);
+  min_cut_value=cutValue(g,min_cut1,cap);
+
+  check(ek_test.flowValue() == min_cut_value &&
+        min_cut_value == 2*flow_value,
+        "The max flow value or the min cut value is wrong.");
+
+  check(checkFlow(g, ek_test.flowMap(), cap, s, t),
+        "The flow is not feasible.");
+
+  CutMap min_cut2(g);
+  ek_test.minCutMap(min_cut2);
+  min_cut_value=cutValue(g,min_cut2,cap);
+
+  check(ek_test.flowValue() == min_cut_value &&
+        min_cut_value == 2*flow_value,
+        "The max flow value or the three min cut values were not doubled.");
+
+
+  ek_test.flowMap(flow);
+
+  NodeIt tmp1(g,s);
+  ++tmp1;
+  if ( tmp1 != INVALID ) s=tmp1;
+
+  NodeIt tmp2(g,t);
+  ++tmp2;
+  if ( tmp2 != INVALID ) t=tmp2;
+
+  ek_test.source(s);
+  ek_test.target(t);
+
+  ek_test.run();
+
+  CutMap min_cut3(g);
+  ek_test.minCutMap(min_cut3);
+  min_cut_value=cutValue(g,min_cut3,cap);
+
+
+  check(ek_test.flowValue() == min_cut_value,
+        "The max flow value or the three min cut values are incorrect.");
+
+  return 0;
+}
