Changeset 2514:57143c09dc20 in lemon-0.x

Timestamp:

11/17/07 21:58:11 (16 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)

Files:

• benchmark/hcube.cc (modified) (1 diff)
• demo/disjoint_paths_demo.cc (modified) (4 diffs)
• demo/sub_graph_adaptor_demo.cc (modified) (2 diffs)
• doc/groups.dox (modified) (1 diff)
• doc/images/flow.eps (deleted)
• doc/images/flow.png (deleted)
• lemon/Makefile.am (modified) (2 diffs)
• lemon/dinitz_sleator_tarjan.h (added)
• lemon/dynamic_tree.h (added)
• lemon/edmonds_karp.h (modified) (10 diffs)
• lemon/goldberg_tarjan.h (added)
• lemon/preflow.h (modified) (3 diffs)
• test/preflow_test.cc (modified) (7 diffs)

benchmark/hcube.cc

@@ -113 +113 @@
    Graph::EdgeMap<int> flow(G);
    
-    Preflow<Graph,int> MF(G,nodes[0],nodes[1<<dim-1],map,flow);
+    Preflow<Graph> MF(G,map,nodes[0],nodes[1<<dim-1]);
     for(int i=0;i<10;i++)
-      MF.run(MF.NO_FLOW);
+      MF.run();
   }
   PrintTime("PREFLOW",T);

demo/disjoint_paths_demo.cc

@@ -64 +64 @@
   Graph::EdgeMap<int> flow(graph);
 
-  Preflow<Graph, int, Capacity> preflow(graph, source, target, capacity, flow); 
+  Preflow<Graph, Capacity> preflow(graph, capacity, source, target); 
   
   preflow.run();
@@ -73 +73 @@
     title("edge disjoint paths").scaleToA4().
     copyright("(C) 2003-2007 LEMON Project").drawArrows().
-    edgeColors(composeMap(functorMap(color), flow)).
+    edgeColors(composeMap(functorMap(color), preflow.flowMap())).
     coords(coords).run();
 
@@ -88 +88 @@
   SGraph::EdgeMap<int> sflow(sgraph);
 
-  Preflow<SGraph, int, SCapacity> spreflow(sgraph, SGraph::outNode(source), 
-                                           SGraph::inNode(target), 
-                                           scapacity, sflow);
+  Preflow<SGraph, SCapacity> spreflow(sgraph, scapacity, 
+                                      SGraph::outNode(source), 
+                                      SGraph::inNode(target));
 
   spreflow.run();
@@ -100 +100 @@
     title("node disjoint paths").scaleToA4().
     copyright("(C) 2003-2007 LEMON Project").drawArrows().
-    edgeColors(composeMap(functorMap(color), sflow)).
+    edgeColors(composeMap(functorMap(color), spreflow.flowMap())).
     coords(SGraph::combinedNodeMap(coords,
                                    shiftMap(coords,

demo/sub_graph_adaptor_demo.cc

@@ -110 +110 @@
 
   ConstMap<Edge, int> const_1_map(1);
-  Graph::EdgeMap<int> flow(g, 0);
   // Max flow between s and t in the graph of tight edges.
-  Preflow<SubGW, int, ConstMap<Edge, int>, Graph::EdgeMap<int> > 
-    preflow(gw, s, t, const_1_map, flow);
+  Preflow<SubGW, ConstMap<Edge, int> > preflow(gw, const_1_map, s, t);
   preflow.run();
 
@@ -121 +119 @@
        << endl;
   for(EdgeIt e(g); e!=INVALID; ++e) 
-    if (flow[e])
+    if (preflow.flow(e) != 0)
       cout << " " << g.id(e) << ", "
            << g.id(g.source(e)) << "--" 

doc/groups.dox

@@ -224 +224 @@
 feasible circulations.
 
-\image html flow.png
-\image latex flow.eps "Graph flow" width=\textwidth
+The maximum flow problem is to find a flow between a single-source and
+single-target that is maximum. Formally, there is \f$G=(V,A)\f$
+directed graph, an \f$c_a:A\rightarrow\mathbf{R}^+_0\f$ capacity
+function and given \f$s, t \in V\f$ source and target node. The
+maximum flow is the solution of the next optimization problem:
+
+\f[ 0 \le f_a \le c_a \f]
+\f[ \sum_{v\in\delta^{-}(u)}f_{vu}=\sum_{v\in\delta^{+}(u)}f_{uv} \quad u \in V \setminus \{s,t\}\f]
+\f[ \max \sum_{v\in\delta^{+}(s)}f_{uv} - \sum_{v\in\delta^{-}(s)}f_{vu}\f]
+
+The lemon contains several algorithms for solve maximum flow problems:
+- \ref lemon::EdmondsKarp "Edmonds-Karp" 
+- \ref lemon::Preflow "Goldberg's Preflow algorithm"
+- \ref lemon::DinitzSleatorTarjan "Dinitz's blocking flow algorithm with dynamic tree"
+- \ref lemon::GoldbergTarjan "Preflow algorithm with dynamic trees"
+
+In most cases the \ref lemon::Preflow "preflow" algorithm provides the
+fastest method to compute the maximum flow. All impelementations
+provides functions for query the minimum cut, which is the dual linear
+programming probelm of the maximum flow.
+
 */
 

lemon/Makefile.am

@@ -53 +53 @@
         lemon/dfs.h \
         lemon/dijkstra.h \
+        lemon/dinitz_sleator_tarjan.h \
         lemon/dist_log.h \
         lemon/dim2.h \
         lemon/dimacs.h \
+        lemon/dynamic_tree.h \
         lemon/edge_set.h \
         lemon/edmonds_karp.h \
@@ -72 +74 @@
         lemon/graph_writer.h \
         lemon/grid_ugraph.h \
+        lemon/goldberg_tarjan.h \
         lemon/hao_orlin.h \
         lemon/hypercube_graph.h \

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);
+    }    
+
+    /// @}
+
   };
 

lemon/preflow.h

@@ -20 +20 @@
 #define LEMON_PREFLOW_H
 
-#include <vector>
-#include <queue>
-
 #include <lemon/error.h>
 #include <lemon/bits/invalid.h>
@@ -28 +25 @@
 #include <lemon/maps.h>
 #include <lemon/graph_utils.h>
+#include <lemon/elevator.h>
 
 /// \file
@@ -33 +31 @@
 /// \brief Implementation of the preflow algorithm.
 
-namespace lemon {
-
-  ///\ingroup max_flow
-  ///\brief %Preflow algorithms class.
+namespace lemon { 
+  
+  /// \brief Default traits class of Preflow class.
   ///
-  ///This class provides an implementation of the \e preflow \e
-  ///algorithm producing a flow of maximum value in a directed
-  ///graph. The preflow algorithms are the fastest known max flow algorithms. 
-  ///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.
+  /// Default traits class of Preflow class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _CapacityMap>
+  struct PreflowDefaultTraits {
+
+    /// \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 eleavator type used by Preflow algorithm.
+    /// 
+    /// The elevator type used by Preflow algorithm.
+    ///
+    /// \sa Elevator
+    /// \sa LinkedElevator
+    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
+    
+    /// \brief Instantiates an Elevator.
+    ///
+    /// This function instantiates a \ref Elevator. 
+    /// \param graph The graph, to which we would like to define the elevator.
+    /// \param max_level The maximum level of the elevator.
+    static Elevator* createElevator(const Graph& graph, int max_level) {
+      return new Elevator(graph, max_level);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+  
+
+  /// \ingroup max_flow
   ///
-  ///After running \ref lemon::Preflow::phase1() "phase1()"
-  ///or \ref lemon::Preflow::run() "run()", the maximal flow
-  ///value can be obtained by calling \ref flowValue(). The minimum
-  ///value cut can be written into a <tt>bool</tt> node map by
-  ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
-  ///the inclusionwise minimum and maximum of the minimum value cuts,
-  ///resp.)
+  /// \brief %Preflow algorithms class.
   ///
-  ///\param Graph The directed graph type the algorithm runs on.
-  ///\param Num The number type of the capacities and the flow values.
-  ///\param CapacityMap The capacity map type.
-  ///\param FlowMap The flow map type.
-  ///\param Tol The tolerance type. 
+  /// This class provides an implementation of the Goldberg's \e
+  /// preflow \e algorithm producing a flow of maximum value in a
+  /// directed graph. The preflow algorithms are the fastest known max
+  /// flow algorithms. The current implementation use a mixture of the
+  /// \e "highest label" and the \e "bound decrease" heuristics.
+  /// The worst case time complexity of the algorithm is \f$O(n^3)\f$.
   ///
-  ///\author Jacint Szabo 
-  ///\todo Second template parameter is superfluous
-  template <typename Graph, typename Num,
-            typename CapacityMap=typename Graph::template EdgeMap<Num>,
-            typename FlowMap=typename Graph::template EdgeMap<Num>,
-            typename Tol=Tolerance<Num> >
+  /// The algorithm consists from two phases. After the first phase
+  /// the maximal flow value and the minimum cut can be obtained. The
+  /// second phase constructs the feasible maximum flow on each edge.
+  ///
+  /// \param _Graph The directed graph type the algorithm runs on.
+  /// \param _CapacityMap The flow map type.
+  /// \param _Traits Traits class to set various data types used by
+  /// the algorithm.  The default traits class is \ref
+  /// PreflowDefaultTraits.  See \ref PreflowDefaultTraits for the
+  /// documentation of a %Preflow traits class. 
+  ///
+  ///\author Jacint Szabo and Balazs Dezso
+#ifdef DOXYGEN
+  template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+  template <typename _Graph, 
+            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+            typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
+#endif
   class Preflow {
-  protected:
-    typedef typename Graph::Node Node;
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::EdgeIt EdgeIt;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
-    typedef typename Graph::InEdgeIt InEdgeIt;
-
-    typedef typename Graph::template NodeMap<Node> NNMap;
-    typedef typename std::vector<Node> VecNode;
-
-    const Graph* _g;
-    Node _source;
-    Node _target;
-    const CapacityMap* _capacity;
-    FlowMap* _flow;
-
-    Tol _surely;
+  public:
     
-    int _node_num;      //the number of nodes of G
-    
-    typename Graph::template NodeMap<int> level;  
-    typename Graph::template NodeMap<Num> excess;
-
-    // constants used for heuristics
-    static const int H0=20;
-    static const int H1=1;
-
-  public:
-
-    ///\ref Exception for the case when s=t.
-
-    ///\ref Exception for the case when the source equals the target.
-    class InvalidArgument : public lemon::LogicError {
+    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::Elevator Elevator;
+    typedef typename Traits::Tolerance Tolerance;
+
+    /// \brief \ref Exception for uninitialized parameters.
+    ///
+    /// This error represents problems in the initialization
+    /// of the parameters of the algorithms.
+    class UninitializedParameter : public lemon::UninitializedParameter {
     public:
       virtual const char* what() const throw() {
-        return "lemon::Preflow::InvalidArgument";
+        return "lemon::Preflow::UninitializedParameter";
       }
     };
+
+  private:
+
+    GRAPH_TYPEDEFS(typename Graph);
+
+    const Graph& _graph;
+    const CapacityMap* _capacity;
+
+    int _node_num;
+
+    Node _source, _target;
+
+    FlowMap* _flow;
+    bool _local_flow;
+
+    Elevator* _level;
+    bool _local_level;
+
+    typedef typename Graph::template NodeMap<Value> ExcessMap;
+    ExcessMap* _excess;
+
+    Tolerance _tolerance;
+
+    bool _phase;
+
+    void createStructures() {
+      _node_num = countNodes(_graph);
+
+      if (!_flow) {
+        _flow = Traits::createFlowMap(_graph);
+        _local_flow = true;
+      }
+      if (!_level) {
+        _level = Traits::createElevator(_graph, _node_num);
+        _local_level = true;
+      }
+      if (!_excess) {
+        _excess = new ExcessMap(_graph);
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+        delete _flow;
+      }
+      if (_local_level) {
+        delete _level;
+      }
+      if (_excess) {
+        delete _excess;
+      }
+    }
+
+  public:
+
+    typedef Preflow Create;
+
+    ///\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 Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
+      typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph&, int) {
+        throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type
+    template <typename _Elevator>
+    struct DefElevator 
+      : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
+      typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefStandardElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph& graph, int max_level) {
+        return new Elevator(graph, max_level);
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type. The Elevator should be standard constructor interface, ie.
+    /// the graph and the maximum level should be passed to it.
+    template <typename _Elevator>
+    struct DefStandardElevator 
+      : public Preflow<Graph, CapacityMap, 
+                       DefStandardElevatorTraits<_Elevator> > {
+      typedef Preflow<Graph, CapacityMap, 
+                      DefStandardElevatorTraits<_Elevator> > 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.
+    Preflow(const Graph& graph, const CapacityMap& capacity, 
+               Node source, Node target) 
+      : _graph(graph), _capacity(&capacity), 
+        _node_num(0), _source(source), _target(target), 
+        _flow(0), _local_flow(false),
+        _level(0), _local_level(false),
+        _excess(0), _tolerance(), _phase() {}
+
+    /// \brief Destrcutor.
+    ///
+    /// Destructor.
+    ~Preflow() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the capacity map.
+    ///
+    /// Sets the capacity map.
+    /// \return \c (*this)
+    Preflow& capacityMap(const CapacityMap& map) {
+      _capacity = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// \return \c (*this)
+    Preflow& 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 elevator.
+    ///
+    /// Sets the elevator.
+    /// \return \c (*this)
+    Preflow& elevator(Elevator& elevator) {
+      if (_local_level) {
+        delete _level;
+        _local_level = false;
+      }
+      _level = &elevator;
+      return *this;
+    }
+
+    /// \brief Returns the elevator.
+    ///
+    /// \return The elevator.
+    const Elevator& elevator() {
+      return *_level;
+    }
+
+    /// \brief Sets the source node.
+    ///
+    /// Sets the source node.
+    /// \return \c (*this)
+    Preflow& source(const Node& node) {
+      _source = node;
+      return *this;
+    }
+
+    /// \brief Sets the target node.
+    ///
+    /// Sets the target node.
+    /// \return \c (*this)
+    Preflow& target(const Node& node) {
+      _target = node;
+      return *this;
+    }
+ 
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    Preflow& 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 one of the member functions called \c
+    /// run().  
+    /// \n
+    /// If you need more control on initial solution or
+    /// execution then you have to call one \ref init() function and then
+    /// the startFirstPhase() and if you need the startSecondPhase().  
     
+    ///@{
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures.
+    ///
+    void init() {
+      createStructures();
+
+      _phase = true;
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _excess->set(n, 0);
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow->set(e, 0);
+      }
+
+      typename Graph::template NodeMap<bool> reached(_graph, false);
+
+      _level->initStart();
+      _level->initAddItem(_target);
+
+      std::vector<Node> queue;
+      reached.set(_source, true);
+
+      queue.push_back(_target);
+      reached.set(_target, true);
+      while (!queue.empty()) {
+        std::vector<Node> nqueue;
+        for (int i = 0; i < int(queue.size()); ++i) {
+          Node n = queue[i];
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node u = _graph.source(e);
+            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+              reached.set(u, true);
+              _level->initAddItem(u);
+              nqueue.push_back(u);
+            }
+          }
+        }
+        queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+        if (_tolerance.positive((*_capacity)[e])) {
+          Node u = _graph.target(e);
+          if ((*_level)[u] == _level->maxLevel()) continue;
+          _flow->set(e, (*_capacity)[e]);
+          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+          if (u != _target && !_level->active(u)) {
+            _level->activate(u);
+          }
+        }
+      }
+    }
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures and sets the initial
+    /// flow to the given \c flowMap. The \c flowMap should contain a
+    /// flow or at least a preflow, ie. in each node excluding the
+    /// target the incoming flow should greater or equal to the
+    /// outgoing flow.
+    /// \return %False when the given \c flowMap is not a preflow. 
+    template <typename FlowMap>
+    bool flowInit(const FlowMap& flowMap) {
+      createStructures();
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow->set(e, flowMap[e]);
+      }
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        Value excess = 0;
+        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+          excess += (*_flow)[e];
+        }
+        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+          excess -= (*_flow)[e];
+        }
+        if (excess < 0 && n != _source) return false;
+        _excess->set(n, excess);
+      }
+
+      typename Graph::template NodeMap<bool> reached(_graph, false);
+
+      _level->initStart();
+      _level->initAddItem(_target);
+
+      std::vector<Node> queue;
+      reached.set(_source, true);
+
+      queue.push_back(_target);
+      reached.set(_target, true);
+      while (!queue.empty()) {
+        _level->initNewLevel();
+        std::vector<Node> nqueue;
+        for (int i = 0; i < int(queue.size()); ++i) {
+          Node n = queue[i];
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node u = _graph.source(e);
+            if (!reached[u] && 
+                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+              reached.set(u, true);
+              _level->initAddItem(u);
+              nqueue.push_back(u);
+            }
+          }
+          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node v = _graph.target(e);
+            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+              reached.set(v, true);
+              _level->initAddItem(v);
+              nqueue.push_back(v);
+            }
+          }
+        }
+        queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+        Value rem = (*_capacity)[e] - (*_flow)[e]; 
+        if (_tolerance.positive(rem)) {
+          Node u = _graph.target(e);
+          if ((*_level)[u] == _level->maxLevel()) continue;
+          _flow->set(e, (*_capacity)[e]);
+          _excess->set(u, (*_excess)[u] + rem);
+          if (u != _target && !_level->active(u)) {
+            _level->activate(u);
+          }
+        }
+      }
+      for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
+        Value rem = (*_flow)[e]; 
+        if (_tolerance.positive(rem)) {
+          Node v = _graph.source(e);
+          if ((*_level)[v] == _level->maxLevel()) continue;
+          _flow->set(e, 0);
+          _excess->set(v, (*_excess)[v] + rem);
+          if (v != _target && !_level->active(v)) {
+            _level->activate(v);
+          }
+        }
+      }
+      return true;
+    }
+
+    /// \brief Starts the first phase of the preflow algorithm.
+    ///
+    /// The preflow algorithm consists of two phases, this method runs
+    /// the first phase. After the first phase the maximum flow value
+    /// and a minimum value cut can already be computed, although a
+    /// maximum flow is not yet obtained. So after calling this method
+    /// \ref flowValue() returns the value of a maximum flow and \ref
+    /// minCut() returns a minimum cut.     
+    /// \pre One of the \ref init() functions should be called. 
+    void startFirstPhase() {
+      _phase = true;
+
+      Node n = _level->highestActive();
+      int level = _level->highestActiveLevel();
+      while (n != INVALID) {
+        int num = _node_num;
+
+        while (num > 0 && n != INVALID) {
+          Value excess = (*_excess)[n];
+          int new_level = _level->maxLevel();
+
+          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.target(e);
+            if ((*_level)[v] < level) {
+              if (!_level->active(v) && v != _target) {
+                _level->activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] + excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push_1;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, (*_capacity)[e]);
+              }
+            } else if (new_level > (*_level)[v]) {
+              new_level = (*_level)[v];
+            }
+          }
+
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Value rem = (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.source(e);
+            if ((*_level)[v] < level) {
+              if (!_level->active(v) && v != _target) {
+                _level->activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] - excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push_1;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, 0);
+              }
+            } else if (new_level > (*_level)[v]) {
+              new_level = (*_level)[v];
+            }
+          }
+
+        no_more_push_1:
+
+          _excess->set(n, excess);
+
+          if (excess != 0) {
+            if (new_level + 1 < _level->maxLevel()) {
+              _level->liftHighestActive(new_level + 1);
+            } else {
+              _level->liftHighestActiveToTop();
+            }
+            if (_level->emptyLevel(level)) {
+              _level->liftToTop(level);
+            }
+          } else {
+            _level->deactivate(n);
+          }
+          
+          n = _level->highestActive();
+          level = _level->highestActiveLevel();
+          --num;
+        }
+        
+        num = _node_num * 20;
+        while (num > 0 && n != INVALID) {
+          Value excess = (*_excess)[n];
+          int new_level = _level->maxLevel();
+
+          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Value rem = (*_capacity)[e] - (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.target(e);
+            if ((*_level)[v] < level) {
+              if (!_level->active(v) && v != _target) {
+                _level->activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] + excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push_2;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, (*_capacity)[e]);
+              }
+            } else if (new_level > (*_level)[v]) {
+              new_level = (*_level)[v];
+            }
+          }
+
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Value rem = (*_flow)[e];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.source(e);
+            if ((*_level)[v] < level) {
+              if (!_level->active(v) && v != _target) {
+                _level->activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(e, (*_flow)[e] - excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push_2;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(e, 0);
+              }
+            } else if (new_level > (*_level)[v]) {
+              new_level = (*_level)[v];
+            }
+          }
+
+        no_more_push_2:
+
+          _excess->set(n, excess);
+
+          if (excess != 0) {
+            if (new_level + 1 < _level->maxLevel()) {
+              _level->liftActiveOn(level, new_level + 1);
+            } else {
+              _level->liftActiveToTop(level);
+            }
+            if (_level->emptyLevel(level)) {
+              _level->liftToTop(level);
+            }
+          } else {
+            _level->deactivate(n);
+          }
+
+          while (level >= 0 && _level->activeFree(level)) {
+            --level;
+          }
+          if (level == -1) {
+            n = _level->highestActive();
+            level = _level->highestActiveLevel();
+          } else {
+            n = _level->activeOn(level);
+          }
+          --num;
+        }
+      }
+    }
+
+    /// \brief Starts the second phase of the preflow algorithm.
+    ///
+    /// The preflow algorithm consists of two phases, this method runs
+    /// the second phase. After calling \ref init() and \ref
+    /// startFirstPhase() and then \ref startSecondPhase(), \ref
+    /// flowMap() return a maximum flow, \ref flowValue() returns the
+    /// value of a maximum flow, \ref minCut() returns a minimum cut
+    /// \pre The \ref init() and startFirstPhase() functions should be
+    /// called before.
+    void startSecondPhase() {
+      _phase = false;
+
+      typename Graph::template NodeMap<bool> reached(_graph, false);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        reached.set(n, (*_level)[n] < _level->maxLevel());
+      }
+
+      _level->initStart();
+      _level->initAddItem(_source);
+ 
+      std::vector<Node> queue;
+      queue.push_back(_source);
+      reached.set(_source, true);
+
+      while (!queue.empty()) {
+        _level->initNewLevel();
+        std::vector<Node> nqueue;
+        for (int i = 0; i < int(queue.size()); ++i) {
+          Node n = queue[i];
+          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node v = _graph.target(e);
+            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+              reached.set(v, true);
+              _level->initAddItem(v);
+              nqueue.push_back(v);
+            }
+          }
+          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            Node u = _graph.source(e);
+            if (!reached[u] && 
+                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+              reached.set(u, true);
+              _level->initAddItem(u);
+              nqueue.push_back(u);
+            }
+          }
+        }
+        queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        if ((*_excess)[n] > 0 && _target != n) {
+          _level->activate(n);
+        }
+      }
+
+      Node n;
+      while ((n = _level->highestActive()) != INVALID) {
+        Value excess = (*_excess)[n];
+        int level = _level->highestActiveLevel();
+        int new_level = _level->maxLevel();
+
+        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+          Value rem = (*_capacity)[e] - (*_flow)[e];
+          if (!_tolerance.positive(rem)) continue;
+          Node v = _graph.target(e);
+          if ((*_level)[v] < level) {
+            if (!_level->active(v) && v != _source) {
+              _level->activate(v);
+            }
+            if (!_tolerance.less(rem, excess)) {
+              _flow->set(e, (*_flow)[e] + excess);
+              _excess->set(v, (*_excess)[v] + excess);
+              excess = 0;
+              goto no_more_push;
+            } else {
+              excess -= rem;
+              _excess->set(v, (*_excess)[v] + rem);
+              _flow->set(e, (*_capacity)[e]);
+            }
+          } else if (new_level > (*_level)[v]) {
+            new_level = (*_level)[v];
+          }
+        }
+
+        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+          Value rem = (*_flow)[e];
+          if (!_tolerance.positive(rem)) continue;
+          Node v = _graph.source(e);
+          if ((*_level)[v] < level) {
+            if (!_level->active(v) && v != _source) {
+              _level->activate(v);
+            }
+            if (!_tolerance.less(rem, excess)) {
+              _flow->set(e, (*_flow)[e] - excess);
+              _excess->set(v, (*_excess)[v] + excess);
+              excess = 0;
+              goto no_more_push;
+            } else {
+              excess -= rem;
+              _excess->set(v, (*_excess)[v] + rem);
+              _flow->set(e, 0);
+            }
+          } else if (new_level > (*_level)[v]) {
+            new_level = (*_level)[v];
+          }
+        }
+
+      no_more_push:
+
+        _excess->set(n, excess);
+      
+        if (excess != 0) {
+          if (new_level + 1 < _level->maxLevel()) {
+            _level->liftHighestActive(new_level + 1);
+          } else {
+            // Calculation error 
+            _level->liftHighestActiveToTop();
+          }
+          if (_level->emptyLevel(level)) {
+            // Calculation error 
+            _level->liftToTop(level);
+          }
+        } else {
+          _level->deactivate(n);
+        }
+
+      }
+    }
+
+    /// \brief Runs the preflow algorithm.  
+    ///
+    /// Runs the preflow algorithm.
+    /// \note pf.run() is just a shortcut of the following code.
+    /// \code
+    ///   pf.init();
+    ///   pf.startFirstPhase();
+    ///   pf.startSecondPhase();
+    /// \endcode
+    void run() {
+      init();
+      startFirstPhase();
+      startSecondPhase();
+    }
+
+    /// \brief Runs the preflow algorithm to compute the minimum cut.  
+    ///
+    /// Runs the preflow algorithm to compute the minimum cut.
+    /// \note pf.runMinCut() is just a shortcut of the following code.
+    /// \code
+    ///   pf.init();
+    ///   pf.startFirstPhase();
+    /// \endcode
+    void runMinCut() {
+      init();
+      startFirstPhase();
+    }
+
+    /// @}
+
+    /// \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.
     
-    ///Indicates the property of the starting flow map.
-    
-    ///Indicates the property of the starting flow map.
-    ///
-    enum FlowEnum{
-      ///indicates an unspecified edge map. \c flow will be 
-      ///set to the constant zero flow in the beginning of
-      ///the algorithm in this case.
-      NO_FLOW,
-      ///constant zero flow
-      ZERO_FLOW,
-      ///any flow, i.e. the sum of the in-flows equals to
-      ///the sum of the out-flows in every node except the \c source and
-      ///the \c target.
-      GEN_FLOW,
-      ///any preflow, i.e. the sum of the in-flows is at 
-      ///least the sum of the out-flows in every node except the \c source.
-      PRE_FLOW
-    };
-
-    ///Indicates the state of the preflow algorithm.
-
-    ///Indicates the state of the preflow algorithm.
-    ///
-    enum StatusEnum {
-      ///before running the algorithm or
-      ///at an unspecified state.
-      AFTER_NOTHING,
-      ///right after running \ref phase1()
-      AFTER_PREFLOW_PHASE_1,      
-      ///after running \ref phase2()
-      AFTER_PREFLOW_PHASE_2
-    };
-    
-  protected: 
-    FlowEnum flow_prop;
-    StatusEnum status; // Do not needle this flag only if necessary.
-    
-  public: 
-    ///The constructor of the class.
-
-    ///The constructor of the class. 
-    ///\param _gr The directed graph the algorithm runs on. 
-    ///\param _s The source node.
-    ///\param _t The target node.
-    ///\param _cap The capacity of the edges. 
-    ///\param _f The flow of the edges. 
-    ///\param _sr Tol class.
-    ///Except the graph, all of these parameters can be reset by
-    ///calling \ref source, \ref target, \ref capacityMap and \ref
-    ///flowMap, resp.
-    Preflow(const Graph& _gr, Node _s, Node _t, 
-            const CapacityMap& _cap, FlowMap& _f,
-            const Tol &_sr=Tol()) :
-        _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
-        _flow(&_f), _surely(_sr),
-        _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), 
-        flow_prop(NO_FLOW), status(AFTER_NOTHING) { 
-        if ( _source==_target )
-          throw InvalidArgument();
-    }
-    
-    ///Give a reference to the tolerance handler class
-
-    ///Give a reference to the tolerance handler class
-    ///\sa Tolerance
-    Tol &tolerance() { return _surely; }
-
-    ///Runs the preflow algorithm.  
-
-    ///Runs the preflow algorithm.
-    ///
-    void run() {
-      phase1(flow_prop);
-      phase2();
-    }
-    
-    ///Runs the preflow algorithm.  
-    
-    ///Runs the preflow algorithm. 
-    ///\pre The starting flow map must be
-    /// - a constant zero flow if \c fp is \c ZERO_FLOW,
-    /// - an arbitrary flow if \c fp is \c GEN_FLOW,
-    /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
-    /// - any map if \c fp is NO_FLOW.
-    ///If the starting flow map is a flow or a preflow then 
-    ///the algorithm terminates faster.
-    void run(FlowEnum fp) {
-      flow_prop=fp;
-      run();
-    }
-      
-    ///Runs the first phase of the preflow algorithm.
-
-    ///The preflow algorithm consists of two phases, this method runs
-    ///the first phase. After the first phase the maximum flow value
-    ///and a minimum value cut can already be computed, although a
-    ///maximum flow is not yet obtained. So after calling this method
-    ///\ref flowValue returns the value of a maximum flow and \ref
-    ///minCut returns a minimum cut.     
-    ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
-    ///value cuts unless calling \ref phase2.  
-    ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
-    ///is needed for this phase.
-    ///\pre The starting flow must be 
-    ///- a constant zero flow if \c fp is \c ZERO_FLOW, 
-    ///- an arbitary flow if \c fp is \c GEN_FLOW, 
-    ///- an arbitary preflow if \c fp is \c PRE_FLOW, 
-    ///- any map if \c fp is NO_FLOW.
-    void phase1(FlowEnum fp)
-    {
-      flow_prop=fp;
-      phase1();
-    }
-
-    
-    ///Runs the first phase of the preflow algorithm.
-
-    ///The preflow algorithm consists of two phases, this method runs
-    ///the first phase. After the first phase the maximum flow value
-    ///and a minimum value cut can already be computed, although a
-    ///maximum flow is not yet obtained. So after calling this method
-    ///\ref flowValue returns the value of a maximum flow and \ref
-    ///minCut returns a minimum cut.
-    ///\warning \ref minMinCut() and \ref maxMinCut() do not
-    ///give minimum value cuts unless calling \ref phase2().
-    ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
-    ///is needed for this phase.
-    void phase1()
-    {
-      int heur0=int(H0*_node_num);  //time while running 'bound decrease'
-      int heur1=int(H1*_node_num);  //time while running 'highest label'
-      int heur=heur1;         //starting time interval (#of relabels)
-      int numrelabel=0;
-
-      bool what_heur=1;
-      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
-
-      bool end=false;
-      //Needed for 'bound decrease', true means no active 
-      //nodes are above bound b.
-
-      int k=_node_num-2;  //bound on the highest level under n containing a node
-      int b=k;    //bound on the highest level under n containing an active node
-
-      VecNode first(_node_num, INVALID);
-      NNMap next(*_g, INVALID);
-
-      NNMap left(*_g, INVALID);
-      NNMap right(*_g, INVALID);
-      VecNode level_list(_node_num,INVALID);
-      //List of the nodes in level i<n, set to n.
-
-      preflowPreproc(first, next, level_list, left, right);
-
-      //Push/relabel on the highest level active nodes.
-      while ( true ) {
-        if ( b == 0 ) {
-          if ( !what_heur && !end && k > 0 ) {
-            b=k;
-            end=true;
-          } else break;
-        }
-
-        if ( first[b]==INVALID ) --b;
-        else {
-          end=false;
-          Node w=first[b];
-          first[b]=next[w];
-          int newlevel=push(w, next, first);
-          if ( excess[w] != 0 ) {
-            relabel(w, newlevel, first, next, level_list, 
-                    left, right, b, k, what_heur);
-          }
-
-          ++numrelabel;
-          if ( numrelabel >= heur ) {
-            numrelabel=0;
-            if ( what_heur ) {
-              what_heur=0;
-              heur=heur0;
-              end=false;
-            } else {
-              what_heur=1;
-              heur=heur1;
-              b=k;
-            }
-          }
-        }
-      }
-      flow_prop=PRE_FLOW;
-      status=AFTER_PREFLOW_PHASE_1;
-    }
-    // Heuristics:
-    //   2 phase
-    //   gap
-    //   list 'level_list' on the nodes on level i implemented by hand
-    //   stack 'active' on the active nodes on level i      
-    //   runs heuristic 'highest label' for H1*n relabels
-    //   runs heuristic 'bound decrease' for H0*n relabels,
-    //        starts with 'highest label'
-    //   Parameters H0 and H1 are initialized to 20 and 1.
-
-
-    ///Runs the second phase of the preflow algorithm.
-
-    ///The preflow algorithm consists of two phases, this method runs
-    ///the second phase. After calling \ref phase1() and then
-    ///\ref phase2(),
-    /// \ref flowMap() return a maximum flow, \ref flowValue
-    ///returns the value of a maximum flow, \ref minCut returns a
-    ///minimum cut, while the methods \ref minMinCut and \ref
-    ///maxMinCut return the inclusionwise minimum and maximum cuts of
-    ///minimum value, resp.  \pre \ref phase1 must be called before.
-    ///
-    /// \todo The inexact computation can cause positive excess on a set of 
-    /// unpushable nodes. We may have to watch the empty level in this case 
-    /// due to avoid the terrible long running time.
-    void phase2()
-    {
-
-      int k=_node_num-2;  //bound on the highest level under n containing a node
-      int b=k;    //bound on the highest level under n of an active node
-
-    
-      VecNode first(_node_num, INVALID);
-      NNMap next(*_g, INVALID); 
-      level.set(_source,0);
-      std::queue<Node> bfs_queue;
-      bfs_queue.push(_source);
-
-      while ( !bfs_queue.empty() ) {
-
-        Node v=bfs_queue.front();
-        bfs_queue.pop();
-        int l=level[v]+1;
-
-        for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
-          if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
-          Node u=_g->source(e);
-          if ( level[u] >= _node_num ) {
-            bfs_queue.push(u);
-            level.set(u, l);
-            if ( excess[u] != 0 ) {
-              next.set(u,first[l]);
-              first[l]=u;
-            }
-          }
-        }
-
-        for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
-          if ( !_surely.positive((*_flow)[e]) ) continue;
-          Node u=_g->target(e);
-          if ( level[u] >= _node_num ) {
-            bfs_queue.push(u);
-            level.set(u, l);
-            if ( excess[u] != 0 ) {
-              next.set(u,first[l]);
-              first[l]=u;
-            }
-          }
-        }
-      }
-      b=_node_num-2;
-
-      while ( true ) {
-
-        if ( b == 0 ) break;
-        if ( first[b]==INVALID ) --b;
-        else {
-          Node w=first[b];
-          first[b]=next[w];
-          int newlevel=push(w,next, first);
-          
-          if ( newlevel == _node_num) {
-            excess.set(w, 0);
-            level.set(w,_node_num);
-          }
-          //relabel
-          if ( excess[w] != 0 ) {
-            level.set(w,++newlevel);
-            next.set(w,first[newlevel]);
-            first[newlevel]=w;
-            b=newlevel;
-          }
-        } 
-      } // while(true)
-      flow_prop=GEN_FLOW;
-      status=AFTER_PREFLOW_PHASE_2;
-    }
-
-    /// Returns the value of the maximum flow.
-
+    ///@{
+
+    /// \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 running \ref phase1.
-    Num flowValue() const {
-      return excess[_target];
-    }
-
-
-    ///Returns a minimum value cut.
-
-    ///Sets \c M to the characteristic vector of a minimum value
-    ///cut. This method can be called both after running \ref
-    ///phase1 and \ref phase2. It is much faster after
-    ///\ref phase1.  \pre M should be a bool-valued node-map. \pre
-    ///If \ref minCut() is called after \ref phase2() then M should
-    ///be initialized to false.
-    template<typename _CutMap>
-    void minCut(_CutMap& M) const {
-      switch ( status ) {
-        case AFTER_PREFLOW_PHASE_1:
-        for(NodeIt v(*_g); v!=INVALID; ++v) {
-          if (level[v] < _node_num) {
-            M.set(v, false);
-          } else {
-            M.set(v, true);
-          }
-        }
-        break;
-        case AFTER_PREFLOW_PHASE_2:
-        minMinCut(M);
-        break;
-        case AFTER_NOTHING:
-        break;
-      }
-    }
-
-    ///Returns the inclusionwise minimum of the minimum value cuts.
-
-    ///Sets \c M 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 s in the residual graph.  \pre M
-    ///should be a node map of bools initialized to false.  \pre \ref
-    ///phase2 should already be run.
-    template<typename _CutMap>
-    void minMinCut(_CutMap& M) const {
-
-      std::queue<Node> queue;
-      M.set(_source,true);
-      queue.push(_source);
-      
-      while (!queue.empty()) {
-        Node w=queue.front();
-        queue.pop();
-        
-        for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-          Node v=_g->target(e);
-          if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) {
-            queue.push(v);
-            M.set(v, true);
-          }
-        }
-        
-        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-          Node v=_g->source(e);
-          if (!M[v] && _surely.positive((*_flow)[e]) ) {
-            queue.push(v);
-            M.set(v, true);
-          }
-        }
-      }
+    /// the maximum flow already after the first phase.
+    Value flowValue() const {
+      return (*_excess)[_target];
+    }
+
+    /// \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 ((*_level)[node] == _level->maxLevel()) == _phase;
+    }
+ 
+    /// \brief Returns a minimum value cut.
+    ///
+    /// Sets the \c cutMap to the characteristic vector of a minimum value
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase(). The result after second
+    /// phase could be changed slightly if inexact computation is used.
+    /// \pre The \c cutMap should be a bool-valued node-map.
+    template <typename CutMap>
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        cutMap.set(n, minCut(n));
+      }
+    }
+
+    /// \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];
     }
     
-    ///Returns the inclusionwise maximum of the minimum value cuts.
-
-    ///Sets \c M to the characteristic vector of the minimum value cut
-    ///which is inclusionwise maximum. It is computed by processing a
-    ///backward bfs from the target node \c t in the residual graph.
-    ///\pre \ref phase2() or run() should already be run.
-    template<typename _CutMap>
-    void maxMinCut(_CutMap& M) const {
-
-      for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
-
-      std::queue<Node> queue;
-
-      M.set(_target,false);
-      queue.push(_target);
-
-      while (!queue.empty()) {
-        Node w=queue.front();
-        queue.pop();
-
-        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-          Node v=_g->source(e);
-          if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) {
-            queue.push(v);
-            M.set(v, false);
-          }
-        }
-
-        for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-          Node v=_g->target(e);
-          if (M[v] && _surely.positive((*_flow)[e]) ) {
-            queue.push(v);
-            M.set(v, false);
-          }
-        }
-      }
-    }
-
-    ///Sets the source node to \c _s.
-
-    ///Sets the source node to \c _s.
-    /// 
-    void source(Node _s) { 
-      _source=_s; 
-      if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
-      status=AFTER_NOTHING; 
-    }
-
-    ///Returns the source node.
-
-    ///Returns the source node.
-    /// 
-    Node source() const { 
-      return _source;
-    }
-
-    ///Sets the target node to \c _t.
-
-    ///Sets the target node to \c _t.
-    ///
-    void target(Node _t) { 
-      _target=_t; 
-      if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
-      status=AFTER_NOTHING; 
-    }
-
-    ///Returns the target node.
-
-    ///Returns the target node.
-    /// 
-    Node target() const { 
-      return _target;
-    }
-
-    /// Sets the edge map of the capacities to _cap.
-
-    /// Sets the edge map of the capacities to _cap.
-    /// 
-    void capacityMap(const CapacityMap& _cap) { 
-      _capacity=&_cap; 
-      status=AFTER_NOTHING; 
-    }
-    /// Returns a reference to capacity map.
-
-    /// Returns a reference to capacity map.
-    /// 
-    const CapacityMap &capacityMap() const { 
-      return *_capacity;
-    }
-
-    /// Sets the edge map of the flows to _flow.
-
-    /// Sets the edge map of the flows to _flow.
-    /// 
-    void flowMap(FlowMap& _f) { 
-      _flow=&_f; 
-      flow_prop=NO_FLOW;
-      status=AFTER_NOTHING; 
-    }
-     
-    /// Returns a reference to flow map.
-
-    /// Returns a reference to flow map.
-    /// 
-    const FlowMap &flowMap() const { 
-      return *_flow;
-    }
-
-  private:
-
-    int push(Node w, NNMap& next, VecNode& first) {
-
-      int lev=level[w];
-      Num exc=excess[w];
-      int newlevel=_node_num;       //bound on the next level of w
-
-      for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-        if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
-        Node v=_g->target(e);
-        
-        if( lev > level[v] ) { //Push is allowed now
-          
-          if ( excess[v] == 0 && v!=_target && v!=_source ) {
-            next.set(v,first[level[v]]);
-            first[level[v]]=v;
-          }
-
-          Num cap=(*_capacity)[e];
-          Num flo=(*_flow)[e];
-          Num remcap=cap-flo;
-          
-          if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
-            
-            _flow->set(e, flo+exc);
-            excess.set(v, excess[v]+exc);
-            exc=0;
-            break;
-
-          } else { //A saturating push.
-            _flow->set(e, cap);
-            excess.set(v, excess[v]+remcap);
-            exc-=remcap;
-          }
-        } else if ( newlevel > level[v] ) newlevel = level[v];
-      } //for out edges wv
-
-      if ( exc != 0 ) {
-        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-          
-          if ( !_surely.positive((*_flow)[e]) ) continue;
-          Node v=_g->source(e);
-          
-          if( lev > level[v] ) { //Push is allowed now
-
-            if ( excess[v] == 0 && v!=_target && v!=_source ) {
-              next.set(v,first[level[v]]);
-              first[level[v]]=v;
-            }
-
-            Num flo=(*_flow)[e];
-
-            if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
-
-              _flow->set(e, flo-exc);
-              excess.set(v, excess[v]+exc);
-              exc=0;
-              break;
-            } else {  //A saturating push.
-
-              excess.set(v, excess[v]+flo);
-              exc-=flo;
-              _flow->set(e,0);
-            }
-          } else if ( newlevel > level[v] ) newlevel = level[v];
-        } //for in edges vw
-
-      } // if w still has excess after the out edge for cycle
-
-      excess.set(w, exc);
-      
-      return newlevel;
-    }
-    
-    
-    
-    void preflowPreproc(VecNode& first, NNMap& next, 
-                        VecNode& level_list, NNMap& left, NNMap& right)
-    {
-      for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
-      std::queue<Node> bfs_queue;
-      
-      if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
-        //Reverse_bfs from t in the residual graph,
-        //to find the starting level.
-        level.set(_target,0);
-        bfs_queue.push(_target);
-        
-        while ( !bfs_queue.empty() ) {
-          
-          Node v=bfs_queue.front();
-          bfs_queue.pop();
-          int l=level[v]+1;
-          
-          for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
-            if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue;
-            Node w=_g->source(e);
-            if ( level[w] == _node_num && w != _source ) {
-              bfs_queue.push(w);
-              Node z=level_list[l];
-              if ( z!=INVALID ) left.set(z,w);
-              right.set(w,z);
-              level_list[l]=w;
-              level.set(w, l);
-            }
-          }
-          
-          for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
-            if ( !_surely.positive((*_flow)[e]) ) continue;
-            Node w=_g->target(e);
-            if ( level[w] == _node_num && w != _source ) {
-              bfs_queue.push(w);
-              Node z=level_list[l];
-              if ( z!=INVALID ) left.set(z,w);
-              right.set(w,z);
-              level_list[l]=w;
-              level.set(w, l);
-            }
-          }
-        } //while
-      } //if
-
-
-      switch (flow_prop) {
-        case NO_FLOW:  
-        for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
-        case ZERO_FLOW:
-        for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
-        
-        //Reverse_bfs from t, to find the starting level.
-        level.set(_target,0);
-        bfs_queue.push(_target);
-        
-        while ( !bfs_queue.empty() ) {
-          
-          Node v=bfs_queue.front();
-          bfs_queue.pop();
-          int l=level[v]+1;
-          
-          for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
-            Node w=_g->source(e);
-            if ( level[w] == _node_num && w != _source ) {
-              bfs_queue.push(w);
-              Node z=level_list[l];
-              if ( z!=INVALID ) left.set(z,w);
-              right.set(w,z);
-              level_list[l]=w;
-              level.set(w, l);
-            }
-          }
-        }
-        
-        //the starting flow
-        for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
-          Num c=(*_capacity)[e];
-          if ( !_surely.positive(c) ) continue;
-          Node w=_g->target(e);
-          if ( level[w] < _node_num ) {
-            if ( excess[w] == 0 && w!=_target ) { //putting into the stack
-              next.set(w,first[level[w]]);
-              first[level[w]]=w;
-            }
-            _flow->set(e, c);
-            excess.set(w, excess[w]+c);
-          }
-        }
-        break;
-
-        case GEN_FLOW:
-        for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
-        {
-          Num exc=0;
-          for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
-          for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
-          if (!_surely.positive(exc)) {
-            exc = 0;
-          }
-          excess.set(_target,exc);
-        }
-
-        //the starting flow
-        for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)  {
-          Num rem=(*_capacity)[e]-(*_flow)[e];
-          if ( !_surely.positive(rem) ) continue;
-          Node w=_g->target(e);
-          if ( level[w] < _node_num ) {
-            if ( excess[w] == 0 && w!=_target ) { //putting into the stack
-              next.set(w,first[level[w]]);
-              first[level[w]]=w;
-            }   
-            _flow->set(e, (*_capacity)[e]);
-            excess.set(w, excess[w]+rem);
-          }
-        }
-        
-        for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
-          if ( !_surely.positive((*_flow)[e]) ) continue;
-          Node w=_g->source(e);
-          if ( level[w] < _node_num ) {
-            if ( excess[w] == 0 && w!=_target ) {
-              next.set(w,first[level[w]]);
-              first[level[w]]=w;
-            }  
-            excess.set(w, excess[w]+(*_flow)[e]);
-            _flow->set(e, 0);
-          }
-        }
-        break;
-
-        case PRE_FLOW:  
-        //the starting flow
-        for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
-          Num rem=(*_capacity)[e]-(*_flow)[e];
-          if ( !_surely.positive(rem) ) continue;
-          Node w=_g->target(e);
-          if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
-        }
-        
-        for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
-          if ( !_surely.positive((*_flow)[e]) ) continue;
-          Node w=_g->source(e);
-          if ( level[w] < _node_num ) _flow->set(e, 0);
-        }
-        
-        //computing the excess
-        for(NodeIt w(*_g); w!=INVALID; ++w) {
-          Num exc=0;
-          for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
-          for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
-          if (!_surely.positive(exc)) {
-            exc = 0;
-          }
-          excess.set(w,exc);
-          
-          //putting the active nodes into the stack
-          int lev=level[w];
-            if ( exc != 0 && lev < _node_num && Node(w) != _target ) {
-              next.set(w,first[lev]);
-              first[lev]=w;
-            }
-        }
-        break;
-      } //switch
-    } //preflowPreproc
-
-
-    void relabel(Node w, int newlevel, VecNode& first, NNMap& next, 
-                 VecNode& level_list, NNMap& left,
-                 NNMap& right, int& b, int& k, bool what_heur )
-    {
-
-      int lev=level[w];
-
-      Node right_n=right[w];
-      Node left_n=left[w];
-
-      //unlacing starts
-      if ( right_n!=INVALID ) {
-        if ( left_n!=INVALID ) {
-          right.set(left_n, right_n);
-          left.set(right_n, left_n);
-        } else {
-          level_list[lev]=right_n;
-          left.set(right_n, INVALID);
-        }
-      } else {
-        if ( left_n!=INVALID ) {
-          right.set(left_n, INVALID);
-        } else {
-          level_list[lev]=INVALID;
-        }
-      }
-      //unlacing ends
-
-      if ( level_list[lev]==INVALID ) {
-
-        //gapping starts
-        for (int i=lev; i!=k ; ) {
-          Node v=level_list[++i];
-          while ( v!=INVALID ) {
-            level.set(v,_node_num);
-            v=right[v];
-          }
-          level_list[i]=INVALID;
-          if ( !what_heur ) first[i]=INVALID;
-        }
-
-        level.set(w,_node_num);
-        b=lev-1;
-        k=b;
-        //gapping ends
-
-      } else {
-
-        if ( newlevel == _node_num ) level.set(w,_node_num);
-        else {
-          level.set(w,++newlevel);
-          next.set(w,first[newlevel]);
-          first[newlevel]=w;
-          if ( what_heur ) b=newlevel;
-          if ( k < newlevel ) ++k;      //now k=newlevel
-          Node z=level_list[newlevel];
-          if ( z!=INVALID ) left.set(z,w);
-          right.set(w,z);
-          left.set(w,INVALID);
-          level_list[newlevel]=w;
-        }
-      }
-    } //relabel
-
-  }; 
-
-  ///\ingroup max_flow
-  ///\brief Function type interface for Preflow algorithm.
-  ///
-  ///Function type interface for Preflow algorithm.
-  ///\sa Preflow
-  template<class GR, class CM, class FM>
-  Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
-                            typename GR::Node source,
-                            typename GR::Node target,
-                            const CM &cap,
-                            FM &flow
-                            )
-  {
-    return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
-  }
-
-} //namespace lemon
-
-#endif //LEMON_PREFLOW_H
+    /// @}    
+  };
+}
+
+#endif

test/preflow_test.cc

@@ -29 +29 @@
 using namespace lemon;
 
-void check_Preflow() 
+void checkPreflow() 
 {
   typedef int VType;
@@ -38 +38 @@
   typedef concepts::ReadMap<Edge,VType> CapMap;
   typedef concepts::ReadWriteMap<Edge,VType> FlowMap;
-  typedef concepts::ReadWriteMap<Node,bool> CutMap;
+  typedef concepts::WriteMap<Node,bool> CutMap;
  
-  typedef Preflow<Graph, int, CapMap, FlowMap> PType;
-
   Graph g;
   Node n;
+  Edge e;
   CapMap cap;
   FlowMap flow;
   CutMap cut;
 
-  PType preflow_test(g,n,n,cap,flow);
-
+  Preflow<Graph, CapMap>::DefFlowMap<FlowMap>::Create preflow_test(g,cap,n,n);
+
+  preflow_test.capacityMap(cap);
+  flow = preflow_test.flowMap();
+  preflow_test.flowMap(flow);
+  preflow_test.source(n);
+  preflow_test.target(n);
+  
+  preflow_test.init();
+  preflow_test.flowInit(cap);
+  preflow_test.startFirstPhase();
+  preflow_test.startSecondPhase();
   preflow_test.run();
+  preflow_test.runMinCut();
+
   preflow_test.flowValue();
-  preflow_test.source(n);
-  preflow_test.flowMap(flow);
-
-  preflow_test.phase1(PType::NO_FLOW);
-  preflow_test.minCut(cut);
-
-  preflow_test.phase2();
-  preflow_test.target(n);
-  preflow_test.capacityMap(cap);
-  preflow_test.minMinCut(cut);
-  preflow_test.maxMinCut(cut);
-}
-
-int cut_value ( SmartGraph& g, SmartGraph::NodeMap<bool>& cut, 
-                SmartGraph::EdgeMap<int>& cap) {
+  preflow_test.minCut(n);
+  preflow_test.minCutMap(cut);
+  preflow_test.flow(e);
+
+}
+
+int cutValue (const SmartGraph& g, 
+              const SmartGraph::NodeMap<bool>& cut, 
+              const SmartGraph::EdgeMap<int>& cap) {
   
   int c=0;
@@ -73 +78 @@
   }
   return c;
+}
+
+bool checkFlow(const SmartGraph& g, 
+               const SmartGraph::EdgeMap<int>& flow, 
+               const SmartGraph::EdgeMap<int>& cap, 
+               SmartGraph::Node s, SmartGraph::Node t) {
+
+  for (SmartGraph::EdgeIt e(g); e != INVALID; ++e) {
+    if (flow[e] < 0 || flow[e] > cap[e]) return false;
+  }
+
+  for (SmartGraph::NodeIt n(g); n != INVALID; ++n) {
+    if (n == s || n == t) continue;
+    int sum = 0;
+    for (SmartGraph::OutEdgeIt e(g, n); e != INVALID; ++e) {
+      sum += flow[e];
+    }
+    for (SmartGraph::InEdgeIt e(g, n); e != INVALID; ++e) {
+      sum -= flow[e];
+    }
+    if (sum != 0) return false;
+  }
+  return true;
 }
 
@@ -86 +114 @@
   typedef Graph::NodeMap<bool> CutMap;
 
-  typedef Preflow<Graph, int> PType;
+  typedef Preflow<Graph, CapMap> PType;
 
   std::string f_name;
@@ -103 +131 @@
   readDimacs(file, g, cap, s, t);
 
-  FlowMap flow(g,0);
-
- 
-
-  PType preflow_test(g, s, t, cap, flow);
-  preflow_test.run(PType::ZERO_FLOW);
+  PType preflow_test(g, cap, s, t);
+  preflow_test.run();
     
-  CutMap min_cut(g,false);
-  preflow_test.minCut(min_cut); 
-  int min_cut_value=cut_value(g,min_cut,cap);
-   
-  CutMap min_min_cut(g,false);
-  preflow_test.minMinCut(min_min_cut); 
-  int min_min_cut_value=cut_value(g,min_min_cut,cap);
-   
-  CutMap max_min_cut(g,false);
-  preflow_test.maxMinCut(max_min_cut); 
-  int max_min_cut_value=cut_value(g,max_min_cut,cap);
-
-  check(preflow_test.flowValue() == min_cut_value &&
-        min_cut_value == min_min_cut_value &&
-        min_min_cut_value == max_min_cut_value,
+  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
+        "The flow is not feasible.");
+
+  CutMap min_cut(g);
+  preflow_test.minCutMap(min_cut); 
+  int min_cut_value=cutValue(g,min_cut,cap);
+  
+  check(preflow_test.flowValue() == min_cut_value,
         "The max flow value is not equal to the three min cut values.");
 
+  FlowMap flow(g);
+  flow = preflow_test.flowMap();
+
   int flow_value=preflow_test.flowValue();
 
-
-
   for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e]; 
-  preflow_test.capacityMap(cap);  
-
-  preflow_test.phase1(PType::PRE_FLOW);
-
-  CutMap min_cut1(g,false);
-  preflow_test.minCut(min_cut1); 
-  min_cut_value=cut_value(g,min_cut1,cap);
+  preflow_test.flowInit(flow);
+  preflow_test.startFirstPhase();
+
+  CutMap min_cut1(g);
+  preflow_test.minCutMap(min_cut1); 
+  min_cut_value=cutValue(g,min_cut1,cap);
    
   check(preflow_test.flowValue() == min_cut_value &&
@@ -144 +161 @@
         "The max flow value or the min cut value is wrong.");
 
-  preflow_test.phase2();
-
-  CutMap min_cut2(g,false);
-  preflow_test.minCut(min_cut2); 
-  min_cut_value=cut_value(g,min_cut2,cap);
+  preflow_test.startSecondPhase();
+
+  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
+        "The flow is not feasible.");
+
+  CutMap min_cut2(g);
+  preflow_test.minCutMap(min_cut2); 
+  min_cut_value=cutValue(g,min_cut2,cap);
    
-  CutMap min_min_cut2(g,false);
-  preflow_test.minMinCut(min_min_cut2); 
-  min_min_cut_value=cut_value(g,min_min_cut2,cap);
- 
-  preflow_test.maxMinCut(max_min_cut); 
-  max_min_cut_value=cut_value(g,max_min_cut,cap);
-
   check(preflow_test.flowValue() == min_cut_value &&
-        min_cut_value == min_min_cut_value &&
-        min_min_cut_value == max_min_cut_value &&
         min_cut_value == 2*flow_value,
         "The max flow value or the three min cut values were not doubled");
 
-
-
-  EdgeIt e(g);
-  for( int i=1; i==10; ++i ) {
-    flow.set(e,0);
-    ++e;
-  }
 
   preflow_test.flowMap(flow); 
@@ -186 +190 @@
   preflow_test.run();
 
-  CutMap min_cut3(g,false);
-  preflow_test.minCut(min_cut3); 
-  min_cut_value=cut_value(g,min_cut3,cap);
+  CutMap min_cut3(g);
+  preflow_test.minCutMap(min_cut3); 
+  min_cut_value=cutValue(g,min_cut3,cap);
    
-  CutMap min_min_cut3(g,false);
-  preflow_test.minMinCut(min_min_cut3); 
-  min_min_cut_value=cut_value(g,min_min_cut3,cap);
-   
-  preflow_test.maxMinCut(max_min_cut); 
-  max_min_cut_value=cut_value(g,max_min_cut,cap);
-
-  check(preflow_test.flowValue() == min_cut_value &&
-        min_cut_value == min_min_cut_value &&
-        min_min_cut_value == max_min_cut_value,
+
+  check(preflow_test.flowValue() == min_cut_value,
         "The max flow value or the three min cut values are incorrect.");
-}
+  
+  return 0;
+}