Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 432)
+++ lemon/Makefile.am	(revision 434)
@@ -22,4 +22,5 @@
         lemon/bfs.h \
         lemon/bin_heap.h \
+        lemon/circulation.h \
         lemon/color.h \
 	lemon/concept_check.h \
@@ -37,4 +38,5 @@
 	lemon/hypercube_graph.h \
 	lemon/kruskal.h \
+	lemon/hao_orlin.h \
 	lemon/lgf_reader.h \
 	lemon/lgf_writer.h \
Index: lemon/bfs.h
===================================================================
--- lemon/bfs.h	(revision 341)
+++ lemon/bfs.h	(revision 421)
@@ -52,5 +52,5 @@
     ///Instantiates a PredMap.
 
-    ///This function instantiates a PredMap.  
+    ///This function instantiates a PredMap.
     ///\param g is the digraph, to which we would like to define the
     ///PredMap.
@@ -81,5 +81,6 @@
     ///The type of the map that indicates which nodes are reached.
 
-    ///The type of the map that indicates which nodes are reached.///It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    ///The type of the map that indicates which nodes are reached.
+    ///It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
     typedef typename Digraph::template NodeMap<bool> ReachedMap;
     ///Instantiates a ReachedMap.
@@ -119,11 +120,5 @@
   ///
   ///\tparam GR The type of the digraph the algorithm runs on.
-  ///The default value is \ref ListDigraph. The value of GR is not used
-  ///directly by \ref Bfs, it is only passed to \ref BfsDefaultTraits.
-  ///\tparam TR Traits class to set various data types used by the algorithm.
-  ///The default traits class is
-  ///\ref BfsDefaultTraits "BfsDefaultTraits<GR>".
-  ///See \ref BfsDefaultTraits for the documentation of
-  ///a Bfs traits class.
+  ///The default type is \ref ListDigraph.
 #ifdef DOXYGEN
   template <typename GR,
@@ -151,5 +146,5 @@
     typedef PredMapPath<Digraph, PredMap> Path;
 
-    ///The traits class.
+    ///The \ref BfsDefaultTraits "traits class" of the algorithm.
     typedef TR Traits;
 
@@ -213,5 +208,5 @@
     typedef Bfs Create;
 
-    ///\name Named template parameters
+    ///\name Named Template Parameters
 
     ///@{
@@ -231,4 +226,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///PredMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetPredMap : public Bfs< Digraph, SetPredMapTraits<T> > {
@@ -250,4 +246,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///DistMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetDistMap : public Bfs< Digraph, SetDistMapTraits<T> > {
@@ -269,4 +266,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///ReachedMap type.
+    ///It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
     template <class T>
     struct SetReachedMap : public Bfs< Digraph, SetReachedMapTraits<T> > {
@@ -288,4 +286,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///ProcessedMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetProcessedMap : public Bfs< Digraph, SetProcessedMapTraits<T> > {
@@ -340,7 +339,8 @@
 
     ///Sets the map that stores the predecessor arcs.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Bfs &predMap(PredMap &m)
@@ -357,7 +357,8 @@
 
     ///Sets the map that indicates which nodes are reached.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Bfs &reachedMap(ReachedMap &m)
@@ -374,7 +375,8 @@
 
     ///Sets the map that indicates which nodes are processed.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Bfs &processedMap(ProcessedMap &m)
@@ -392,7 +394,8 @@
     ///Sets the map that stores the distances of the nodes calculated by
     ///the algorithm.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Bfs &distMap(DistMap &m)
@@ -408,20 +411,17 @@
   public:
 
-    ///\name Execution control
-    ///The simplest way to execute the algorithm is to use
-    ///one of the member functions called \ref lemon::Bfs::run() "run()".
-    ///\n
-    ///If you need more control on the execution, first you must call
-    ///\ref lemon::Bfs::init() "init()", then you can add several source
-    ///nodes with \ref lemon::Bfs::addSource() "addSource()".
-    ///Finally \ref lemon::Bfs::start() "start()" will perform the
-    ///actual path computation.
+    ///\name Execution Control
+    ///The simplest way to execute the BFS algorithm is to use one of the
+    ///member functions called \ref run(Node) "run()".\n
+    ///If you need more control on the execution, first you have to call
+    ///\ref init(), then you can add several source nodes with
+    ///\ref addSource(). Finally the actual path computation can be
+    ///performed with one of the \ref start() functions.
 
     ///@{
 
+    ///\brief Initializes the internal data structures.
+    ///
     ///Initializes the internal data structures.
-
-    ///Initializes the internal data structures.
-    ///
     void init()
     {
@@ -557,14 +557,14 @@
     }
 
-    ///\brief Returns \c false if there are nodes
-    ///to be processed.
-    ///
-    ///Returns \c false if there are nodes
-    ///to be processed in the queue.
+    ///Returns \c false if there are nodes to be processed.
+
+    ///Returns \c false if there are nodes to be processed
+    ///in the queue.
     bool emptyQueue() const { return _queue_tail==_queue_head; }
 
     ///Returns the number of the nodes to be processed.
 
-    ///Returns the number of the nodes to be processed in the queue.
+    ///Returns the number of the nodes to be processed
+    ///in the queue.
     int queueSize() const { return _queue_head-_queue_tail; }
 
@@ -731,8 +731,8 @@
 
     ///\name Query Functions
-    ///The result of the %BFS algorithm can be obtained using these
+    ///The results of the BFS algorithm can be obtained using these
     ///functions.\n
-    ///Either \ref lemon::Bfs::run() "run()" or \ref lemon::Bfs::start()
-    ///"start()" must be called before using them.
+    ///Either \ref run(Node) "run()" or \ref start() should be called
+    ///before using them.
 
     ///@{
@@ -742,8 +742,8 @@
     ///Returns the shortest path to a node.
     ///
-    ///\warning \c t should be reachable from the root(s).
-    ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\warning \c t should be reached from the root(s).
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Path path(Node t) const { return Path(*G, *_pred, t); }
 
@@ -752,9 +752,9 @@
     ///Returns the distance of a node from the root(s).
     ///
-    ///\warning If node \c v is not reachable from the root(s), then
+    ///\warning If node \c v is not reached from the root(s), then
     ///the return value of this function is undefined.
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     int dist(Node v) const { return (*_dist)[v]; }
 
@@ -763,12 +763,12 @@
     ///This function returns the 'previous arc' of the shortest path
     ///tree for the node \c v, i.e. it returns the last arc of a
-    ///shortest path from the root(s) to \c v. It is \c INVALID if \c v
-    ///is not reachable from the root(s) or if \c v is a root.
+    ///shortest path from a root to \c v. It is \c INVALID if \c v
+    ///is not reached from the root(s) or if \c v is a root.
     ///
     ///The shortest path tree used here is equal to the shortest path
     ///tree used in \ref predNode().
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Arc predArc(Node v) const { return (*_pred)[v];}
 
@@ -777,12 +777,12 @@
     ///This function returns the 'previous node' of the shortest path
     ///tree for the node \c v, i.e. it returns the last but one node
-    ///from a shortest path from the root(s) to \c v. It is \c INVALID
-    ///if \c v is not reachable from the root(s) or if \c v is a root.
+    ///from a shortest path from a root to \c v. It is \c INVALID
+    ///if \c v is not reached from the root(s) or if \c v is a root.
     ///
     ///The shortest path tree used here is equal to the shortest path
     ///tree used in \ref predArc().
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Node predNode(Node v) const { return (*_pred)[v]==INVALID ? INVALID:
                                   G->source((*_pred)[v]); }
@@ -794,5 +794,5 @@
     ///of the nodes calculated by the algorithm.
     ///
-    ///\pre Either \ref run() or \ref init()
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     const DistMap &distMap() const { return *_dist;}
@@ -804,12 +804,13 @@
     ///arcs, which form the shortest path tree.
     ///
-    ///\pre Either \ref run() or \ref init()
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     const PredMap &predMap() const { return *_pred;}
 
-    ///Checks if a node is reachable from the root(s).
-
-    ///Returns \c true if \c v is reachable from the root(s).
-    ///\pre Either \ref run() or \ref start()
+    ///Checks if a node is reached from the root(s).
+
+    ///Returns \c true if \c v is reached from the root(s).
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     bool reached(Node v) const { return (*_reached)[v]; }
@@ -957,6 +958,6 @@
   /// This auxiliary class is created to implement the
   /// \ref bfs() "function-type interface" of \ref Bfs algorithm.
-  /// It does not have own \ref run() method, it uses the functions
-  /// and features of the plain \ref Bfs.
+  /// It does not have own \ref run(Node) "run()" method, it uses the
+  /// functions and features of the plain \ref Bfs.
   ///
   /// This class should only be used through the \ref bfs() function,
@@ -1178,5 +1179,5 @@
   ///  bool reached = bfs(g).path(p).dist(d).run(s,t);
   ///\endcode
-  ///\warning Don't forget to put the \ref BfsWizard::run() "run()"
+  ///\warning Don't forget to put the \ref BfsWizard::run(Node) "run()"
   ///to the end of the parameter list.
   ///\sa BfsWizard
@@ -1364,5 +1365,5 @@
     typedef BfsVisit Create;
 
-    /// \name Named template parameters
+    /// \name Named Template Parameters
 
     ///@{
@@ -1406,7 +1407,8 @@
     ///
     /// Sets the map that indicates which nodes are reached.
-    /// If you don't use this function before calling \ref run(),
-    /// it will allocate one. The destructor deallocates this
-    /// automatically allocated map, of course.
+    /// If you don't use this function before calling \ref run(Node) "run()"
+    /// or \ref init(), an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
     /// \return <tt> (*this) </tt>
     BfsVisit &reachedMap(ReachedMap &m) {
@@ -1421,14 +1423,11 @@
   public:
 
-    /// \name Execution control
-    /// The simplest way to execute the algorithm is to use
-    /// one of the member functions called \ref lemon::BfsVisit::run()
-    /// "run()".
-    /// \n
-    /// If you need more control on the execution, first you must call
-    /// \ref lemon::BfsVisit::init() "init()", then you can add several
-    /// source nodes with \ref lemon::BfsVisit::addSource() "addSource()".
-    /// Finally \ref lemon::BfsVisit::start() "start()" will perform the
-    /// actual path computation.
+    /// \name Execution Control
+    /// The simplest way to execute the BFS algorithm is to use one of the
+    /// member functions called \ref run(Node) "run()".\n
+    /// If you need more control on the execution, first you have to call
+    /// \ref init(), then you can add several source nodes with
+    /// \ref addSource(). Finally the actual path computation can be
+    /// performed with one of the \ref start() functions.
 
     /// @{
@@ -1730,15 +1729,16 @@
 
     /// \name Query Functions
-    /// The result of the %BFS algorithm can be obtained using these
+    /// The results of the BFS algorithm can be obtained using these
     /// functions.\n
-    /// Either \ref lemon::BfsVisit::run() "run()" or
-    /// \ref lemon::BfsVisit::start() "start()" must be called before
-    /// using them.
+    /// Either \ref run(Node) "run()" or \ref start() should be called
+    /// before using them.
+
     ///@{
 
-    /// \brief Checks if a node is reachable from the root(s).
-    ///
-    /// Returns \c true if \c v is reachable from the root(s).
-    /// \pre Either \ref run() or \ref start()
+    /// \brief Checks if a node is reached from the root(s).
+    ///
+    /// Returns \c true if \c v is reached from the root(s).
+    ///
+    /// \pre Either \ref run(Node) "run()" or \ref init()
     /// must be called before using this function.
     bool reached(Node v) { return (*_reached)[v]; }
Index: lemon/circulation.h
===================================================================
--- lemon/circulation.h	(revision 417)
+++ lemon/circulation.h	(revision 417)
@@ -0,0 +1,755 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * 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_CIRCULATION_H
+#define LEMON_CIRCULATION_H
+
+#include <lemon/tolerance.h>
+#include <lemon/elevator.h>
+
+///\ingroup max_flow
+///\file
+///\brief Push-relabel algorithm for finding a feasible circulation.
+///
+namespace lemon {
+
+  /// \brief Default traits class of Circulation class.
+  ///
+  /// Default traits class of Circulation class.
+  /// \tparam _Diraph Digraph type.
+  /// \tparam _LCapMap Lower bound capacity map type.
+  /// \tparam _UCapMap Upper bound capacity map type.
+  /// \tparam _DeltaMap Delta map type.
+  template <typename _Diraph, typename _LCapMap,
+            typename _UCapMap, typename _DeltaMap>
+  struct CirculationDefaultTraits {
+
+    /// \brief The type of the digraph the algorithm runs on.
+    typedef _Diraph Digraph;
+
+    /// \brief The type of the map that stores the circulation lower
+    /// bound.
+    ///
+    /// The type of the map that stores the circulation lower bound.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _LCapMap LCapMap;
+
+    /// \brief The type of the map that stores the circulation upper
+    /// bound.
+    ///
+    /// The type of the map that stores the circulation upper bound.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _UCapMap UCapMap;
+
+    /// \brief The type of the map that stores the lower bound for
+    /// the supply of the nodes.
+    ///
+    /// The type of the map that stores the lower bound for the supply
+    /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
+    /// concept.
+    typedef _DeltaMap DeltaMap;
+
+    /// \brief The type of the flow values.
+    typedef typename DeltaMap::Value Value;
+
+    /// \brief The type of the map that stores the flow values.
+    ///
+    /// The type of the map 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 elevator type used by the algorithm.
+    ///
+    /// The elevator type used by the algorithm.
+    ///
+    /// \sa Elevator
+    /// \sa LinkedElevator
+    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
+
+    /// \brief Instantiates an Elevator.
+    ///
+    /// This function instantiates an \ref Elevator.
+    /// \param digraph The digraph, to which we would like to define
+    /// the elevator.
+    /// \param max_level The maximum level of the elevator.
+    static Elevator* createElevator(const Digraph& digraph, int max_level) {
+      return new Elevator(digraph, max_level);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef lemon::Tolerance<Value> Tolerance;
+
+  };
+
+  /**
+     \brief Push-relabel algorithm for the network circulation problem.
+
+     \ingroup max_flow
+     This class implements a push-relabel algorithm for the network
+     circulation problem.
+     It is to find a feasible circulation when lower and upper bounds
+     are given for the flow values on the arcs and lower bounds
+     are given for the supply values of the nodes.
+
+     The exact formulation of this problem is the following.
+     Let \f$G=(V,A)\f$ be a digraph,
+     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
+     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
+     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
+     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
+     \geq delta(v) \quad \forall v\in V, \f]
+     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
+     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
+     \f$v\f$. It can be either positive or negative, however note that
+     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
+     have a feasible solution.
+
+     \note A special case of this problem is when
+     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
+     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
+     Thus a feasible solution for the
+     \ref min_cost_flow "minimum cost flow" problem can be calculated
+     in this way.
+
+     \tparam _Digraph The type of the digraph the algorithm runs on.
+     \tparam _LCapMap The type of the lower bound capacity map. The default
+     map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
+     \tparam _UCapMap The type of the upper bound capacity map. The default
+     map type is \c _LCapMap.
+     \tparam _DeltaMap The type of the map that stores the lower bound
+     for the supply of the nodes. The default map type is
+     \c _Digraph::ArcMap<_UCapMap::Value>.
+  */
+#ifdef DOXYGEN
+template< typename _Digraph,
+          typename _LCapMap,
+          typename _UCapMap,
+          typename _DeltaMap,
+          typename _Traits >
+#else
+template< typename _Digraph,
+          typename _LCapMap = typename _Digraph::template ArcMap<int>,
+          typename _UCapMap = _LCapMap,
+          typename _DeltaMap = typename _Digraph::
+                               template NodeMap<typename _UCapMap::Value>,
+          typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
+                                                    _UCapMap, _DeltaMap> >
+#endif
+  class Circulation {
+  public:
+
+    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
+    typedef _Traits Traits;
+    ///The type of the digraph the algorithm runs on.
+    typedef typename Traits::Digraph Digraph;
+    ///The type of the flow values.
+    typedef typename Traits::Value Value;
+
+    /// The type of the lower bound capacity map.
+    typedef typename Traits::LCapMap LCapMap;
+    /// The type of the upper bound capacity map.
+    typedef typename Traits::UCapMap UCapMap;
+    /// \brief The type of the map that stores the lower bound for
+    /// the supply of the nodes.
+    typedef typename Traits::DeltaMap DeltaMap;
+    ///The type of the flow map.
+    typedef typename Traits::FlowMap FlowMap;
+
+    ///The type of the elevator.
+    typedef typename Traits::Elevator Elevator;
+    ///The type of the tolerance.
+    typedef typename Traits::Tolerance Tolerance;
+
+  private:
+
+    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+    const Digraph &_g;
+    int _node_num;
+
+    const LCapMap *_lo;
+    const UCapMap *_up;
+    const DeltaMap *_delta;
+
+    FlowMap *_flow;
+    bool _local_flow;
+
+    Elevator* _level;
+    bool _local_level;
+
+    typedef typename Digraph::template NodeMap<Value> ExcessMap;
+    ExcessMap* _excess;
+
+    Tolerance _tol;
+    int _el;
+
+  public:
+
+    typedef Circulation Create;
+
+    ///\name Named Template Parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct SetFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Digraph&) {
+        LEMON_ASSERT(false, "FlowMap is not initialized");
+        return 0; // ignore warnings
+      }
+    };
+
+    /// \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 SetFlowMap
+      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
+                           SetFlowMapTraits<_FlowMap> > {
+      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
+                          SetFlowMapTraits<_FlowMap> > Create;
+    };
+
+    template <typename _Elevator>
+    struct SetElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Digraph&, int) {
+        LEMON_ASSERT(false, "Elevator is not initialized");
+        return 0; // ignore warnings
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type. If this named parameter is used, then an external
+    /// elevator object must be passed to the algorithm using the
+    /// \ref elevator(Elevator&) "elevator()" function before calling
+    /// \ref run() or \ref init().
+    /// \sa SetStandardElevator
+    template <typename _Elevator>
+    struct SetElevator
+      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
+                           SetElevatorTraits<_Elevator> > {
+      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
+                          SetElevatorTraits<_Elevator> > Create;
+    };
+
+    template <typename _Elevator>
+    struct SetStandardElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Digraph& digraph, int max_level) {
+        return new Elevator(digraph, max_level);
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type with automatic allocation
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type with automatic allocation.
+    /// The Elevator should have standard constructor interface to be
+    /// able to automatically created by the algorithm (i.e. the
+    /// digraph and the maximum level should be passed to it).
+    /// However an external elevator object could also be passed to the
+    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
+    /// before calling \ref run() or \ref init().
+    /// \sa SetElevator
+    template <typename _Elevator>
+    struct SetStandardElevator
+      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
+                       SetStandardElevatorTraits<_Elevator> > {
+      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
+                      SetStandardElevatorTraits<_Elevator> > Create;
+    };
+
+    /// @}
+
+  protected:
+
+    Circulation() {}
+
+  public:
+
+    /// The constructor of the class.
+
+    /// The constructor of the class.
+    /// \param g The digraph the algorithm runs on.
+    /// \param lo The lower bound capacity of the arcs.
+    /// \param up The upper bound capacity of the arcs.
+    /// \param delta The lower bound for the supply of the nodes.
+    Circulation(const Digraph &g,const LCapMap &lo,
+                const UCapMap &up,const DeltaMap &delta)
+      : _g(g), _node_num(),
+        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
+        _level(0), _local_level(false), _excess(0), _el() {}
+
+    /// Destructor.
+    ~Circulation() {
+      destroyStructures();
+    }
+
+
+  private:
+
+    void createStructures() {
+      _node_num = _el = countNodes(_g);
+
+      if (!_flow) {
+        _flow = Traits::createFlowMap(_g);
+        _local_flow = true;
+      }
+      if (!_level) {
+        _level = Traits::createElevator(_g, _node_num);
+        _local_level = true;
+      }
+      if (!_excess) {
+        _excess = new ExcessMap(_g);
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+        delete _flow;
+      }
+      if (_local_level) {
+        delete _level;
+      }
+      if (_excess) {
+        delete _excess;
+      }
+    }
+
+  public:
+
+    /// Sets the lower bound capacity map.
+
+    /// Sets the lower bound capacity map.
+    /// \return <tt>(*this)</tt>
+    Circulation& lowerCapMap(const LCapMap& map) {
+      _lo = &map;
+      return *this;
+    }
+
+    /// Sets the upper bound capacity map.
+
+    /// Sets the upper bound capacity map.
+    /// \return <tt>(*this)</tt>
+    Circulation& upperCapMap(const LCapMap& map) {
+      _up = &map;
+      return *this;
+    }
+
+    /// Sets the lower bound map for the supply of the nodes.
+
+    /// Sets the lower bound map for the supply of the nodes.
+    /// \return <tt>(*this)</tt>
+    Circulation& deltaMap(const DeltaMap& map) {
+      _delta = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// If you don't use this function before calling \ref run() or
+    /// \ref init(), an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
+    /// \return <tt>(*this)</tt>
+    Circulation& flowMap(FlowMap& map) {
+      if (_local_flow) {
+        delete _flow;
+        _local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Sets the elevator used by algorithm.
+    ///
+    /// Sets the elevator used by algorithm.
+    /// If you don't use this function before calling \ref run() or
+    /// \ref init(), an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated elevator,
+    /// of course.
+    /// \return <tt>(*this)</tt>
+    Circulation& elevator(Elevator& elevator) {
+      if (_local_level) {
+        delete _level;
+        _local_level = false;
+      }
+      _level = &elevator;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the elevator.
+    ///
+    /// Returns a const reference to the elevator.
+    ///
+    /// \pre Either \ref run() or \ref init() must be called before
+    /// using this function.
+    const Elevator& elevator() {
+      return *_level;
+    }
+
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    Circulation& tolerance(const Tolerance& tolerance) const {
+      _tol = tolerance;
+      return *this;
+    }
+
+    /// \brief Returns a const reference to the tolerance.
+    ///
+    /// Returns a const reference to the tolerance.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    }
+
+    /// \name Execution Control
+    /// The simplest way to execute the algorithm is to call \ref run().\n
+    /// If you need more control on the initial solution or the execution,
+    /// first you have to call one of the \ref init() functions, then
+    /// the \ref start() function.
+
+    ///@{
+
+    /// Initializes the internal data structures.
+
+    /// Initializes the internal data structures and sets all flow values
+    /// to the lower bound.
+    void init()
+    {
+      createStructures();
+
+      for(NodeIt n(_g);n!=INVALID;++n) {
+        _excess->set(n, (*_delta)[n]);
+      }
+
+      for (ArcIt e(_g);e!=INVALID;++e) {
+        _flow->set(e, (*_lo)[e]);
+        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
+        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
+      }
+
+      // global relabeling tested, but in general case it provides
+      // worse performance for random digraphs
+      _level->initStart();
+      for(NodeIt n(_g);n!=INVALID;++n)
+        _level->initAddItem(n);
+      _level->initFinish();
+      for(NodeIt n(_g);n!=INVALID;++n)
+        if(_tol.positive((*_excess)[n]))
+          _level->activate(n);
+    }
+
+    /// Initializes the internal data structures using a greedy approach.
+
+    /// Initializes the internal data structures using a greedy approach
+    /// to construct the initial solution.
+    void greedyInit()
+    {
+      createStructures();
+
+      for(NodeIt n(_g);n!=INVALID;++n) {
+        _excess->set(n, (*_delta)[n]);
+      }
+
+      for (ArcIt e(_g);e!=INVALID;++e) {
+        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
+          _flow->set(e, (*_up)[e]);
+          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
+          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
+        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
+          _flow->set(e, (*_lo)[e]);
+          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
+          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
+        } else {
+          Value fc = -(*_excess)[_g.target(e)];
+          _flow->set(e, fc);
+          _excess->set(_g.target(e), 0);
+          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
+        }
+      }
+
+      _level->initStart();
+      for(NodeIt n(_g);n!=INVALID;++n)
+        _level->initAddItem(n);
+      _level->initFinish();
+      for(NodeIt n(_g);n!=INVALID;++n)
+        if(_tol.positive((*_excess)[n]))
+          _level->activate(n);
+    }
+
+    ///Executes the algorithm
+
+    ///This function executes the algorithm.
+    ///
+    ///\return \c true if a feasible circulation is found.
+    ///
+    ///\sa barrier()
+    ///\sa barrierMap()
+    bool start()
+    {
+
+      Node act;
+      Node bact=INVALID;
+      Node last_activated=INVALID;
+      while((act=_level->highestActive())!=INVALID) {
+        int actlevel=(*_level)[act];
+        int mlevel=_node_num;
+        Value exc=(*_excess)[act];
+
+        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
+          Node v = _g.target(e);
+          Value fc=(*_up)[e]-(*_flow)[e];
+          if(!_tol.positive(fc)) continue;
+          if((*_level)[v]<actlevel) {
+            if(!_tol.less(fc, exc)) {
+              _flow->set(e, (*_flow)[e] + exc);
+              _excess->set(v, (*_excess)[v] + exc);
+              if(!_level->active(v) && _tol.positive((*_excess)[v]))
+                _level->activate(v);
+              _excess->set(act,0);
+              _level->deactivate(act);
+              goto next_l;
+            }
+            else {
+              _flow->set(e, (*_up)[e]);
+              _excess->set(v, (*_excess)[v] + fc);
+              if(!_level->active(v) && _tol.positive((*_excess)[v]))
+                _level->activate(v);
+              exc-=fc;
+            }
+          }
+          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
+        }
+        for(InArcIt e(_g,act);e!=INVALID; ++e) {
+          Node v = _g.source(e);
+          Value fc=(*_flow)[e]-(*_lo)[e];
+          if(!_tol.positive(fc)) continue;
+          if((*_level)[v]<actlevel) {
+            if(!_tol.less(fc, exc)) {
+              _flow->set(e, (*_flow)[e] - exc);
+              _excess->set(v, (*_excess)[v] + exc);
+              if(!_level->active(v) && _tol.positive((*_excess)[v]))
+                _level->activate(v);
+              _excess->set(act,0);
+              _level->deactivate(act);
+              goto next_l;
+            }
+            else {
+              _flow->set(e, (*_lo)[e]);
+              _excess->set(v, (*_excess)[v] + fc);
+              if(!_level->active(v) && _tol.positive((*_excess)[v]))
+                _level->activate(v);
+              exc-=fc;
+            }
+          }
+          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
+        }
+
+        _excess->set(act, exc);
+        if(!_tol.positive(exc)) _level->deactivate(act);
+        else if(mlevel==_node_num) {
+          _level->liftHighestActiveToTop();
+          _el = _node_num;
+          return false;
+        }
+        else {
+          _level->liftHighestActive(mlevel+1);
+          if(_level->onLevel(actlevel)==0) {
+            _el = actlevel;
+            return false;
+          }
+        }
+      next_l:
+        ;
+      }
+      return true;
+    }
+
+    /// Runs the algorithm.
+
+    /// This function runs the algorithm.
+    ///
+    /// \return \c true if a feasible circulation is found.
+    ///
+    /// \note Apart from the return value, c.run() is just a shortcut of
+    /// the following code.
+    /// \code
+    ///   c.greedyInit();
+    ///   c.start();
+    /// \endcode
+    bool run() {
+      greedyInit();
+      return start();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The results of the circulation algorithm can be obtained using
+    /// these functions.\n
+    /// Either \ref run() or \ref start() should be called before
+    /// using them.
+
+    ///@{
+
+    /// \brief Returns the flow on the given arc.
+    ///
+    /// Returns the flow on the given arc.
+    ///
+    /// \pre Either \ref run() or \ref init() must be called before
+    /// using this function.
+    Value flow(const Arc& arc) const {
+      return (*_flow)[arc];
+    }
+
+    /// \brief Returns a const reference to the flow map.
+    ///
+    /// Returns a const reference to the arc map storing the found flow.
+    ///
+    /// \pre Either \ref run() or \ref init() must be called before
+    /// using this function.
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /**
+       \brief Returns \c true if the given node is in a barrier.
+
+       Barrier is a set \e B of nodes for which
+
+       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
+           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
+
+       holds. The existence of a set with this property prooves that a
+       feasible circualtion cannot exist.
+
+       This function returns \c true if the given node is in the found
+       barrier. If a feasible circulation is found, the function
+       gives back \c false for every node.
+
+       \pre Either \ref run() or \ref init() must be called before
+       using this function.
+
+       \sa barrierMap()
+       \sa checkBarrier()
+    */
+    bool barrier(const Node& node)
+    {
+      return (*_level)[node] >= _el;
+    }
+
+    /// \brief Gives back a barrier.
+    ///
+    /// This function sets \c bar to the characteristic vector of the
+    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
+    /// node map with \c bool (or convertible) value type.
+    ///
+    /// If a feasible circulation is found, the function gives back an
+    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
+    ///
+    /// \note This function calls \ref barrier() for each node,
+    /// so it runs in \f$O(n)\f$ time.
+    ///
+    /// \pre Either \ref run() or \ref init() must be called before
+    /// using this function.
+    ///
+    /// \sa barrier()
+    /// \sa checkBarrier()
+    template<class BarrierMap>
+    void barrierMap(BarrierMap &bar)
+    {
+      for(NodeIt n(_g);n!=INVALID;++n)
+        bar.set(n, (*_level)[n] >= _el);
+    }
+
+    /// @}
+
+    /// \name Checker Functions
+    /// The feasibility of the results can be checked using
+    /// these functions.\n
+    /// Either \ref run() or \ref start() should be called before
+    /// using them.
+
+    ///@{
+
+    ///Check if the found flow is a feasible circulation
+
+    ///Check if the found flow is a feasible circulation,
+    ///
+    bool checkFlow() {
+      for(ArcIt e(_g);e!=INVALID;++e)
+        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
+      for(NodeIt n(_g);n!=INVALID;++n)
+        {
+          Value dif=-(*_delta)[n];
+          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
+          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
+          if(_tol.negative(dif)) return false;
+        }
+      return true;
+    }
+
+    ///Check whether or not the last execution provides a barrier
+
+    ///Check whether or not the last execution provides a barrier.
+    ///\sa barrier()
+    ///\sa barrierMap()
+    bool checkBarrier()
+    {
+      Value delta=0;
+      for(NodeIt n(_g);n!=INVALID;++n)
+        if(barrier(n))
+          delta-=(*_delta)[n];
+      for(ArcIt e(_g);e!=INVALID;++e)
+        {
+          Node s=_g.source(e);
+          Node t=_g.target(e);
+          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
+          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
+        }
+      return _tol.negative(delta);
+    }
+
+    /// @}
+
+  };
+
+}
+
+#endif
Index: lemon/dfs.h
===================================================================
--- lemon/dfs.h	(revision 319)
+++ lemon/dfs.h	(revision 421)
@@ -120,11 +120,5 @@
   ///
   ///\tparam GR The type of the digraph the algorithm runs on.
-  ///The default value is \ref ListDigraph. The value of GR is not used
-  ///directly by \ref Dfs, it is only passed to \ref DfsDefaultTraits.
-  ///\tparam TR Traits class to set various data types used by the algorithm.
-  ///The default traits class is
-  ///\ref DfsDefaultTraits "DfsDefaultTraits<GR>".
-  ///See \ref DfsDefaultTraits for the documentation of
-  ///a Dfs traits class.
+  ///The default type is \ref ListDigraph.
 #ifdef DOXYGEN
   template <typename GR,
@@ -152,5 +146,5 @@
     typedef PredMapPath<Digraph, PredMap> Path;
 
-    ///The traits class.
+    ///The \ref DfsDefaultTraits "traits class" of the algorithm.
     typedef TR Traits;
 
@@ -231,4 +225,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///PredMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetPredMap : public Dfs<Digraph, SetPredMapTraits<T> > {
@@ -250,4 +245,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///DistMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetDistMap : public Dfs< Digraph, SetDistMapTraits<T> > {
@@ -269,4 +265,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///ReachedMap type.
+    ///It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
     template <class T>
     struct SetReachedMap : public Dfs< Digraph, SetReachedMapTraits<T> > {
@@ -288,4 +285,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///ProcessedMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetProcessedMap : public Dfs< Digraph, SetProcessedMapTraits<T> > {
@@ -339,7 +337,8 @@
 
     ///Sets the map that stores the predecessor arcs.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dfs &predMap(PredMap &m)
@@ -356,7 +355,8 @@
 
     ///Sets the map that indicates which nodes are reached.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dfs &reachedMap(ReachedMap &m)
@@ -373,7 +373,8 @@
 
     ///Sets the map that indicates which nodes are processed.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dfs &processedMap(ProcessedMap &m)
@@ -391,7 +392,8 @@
     ///Sets the map that stores the distances of the nodes calculated by
     ///the algorithm.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dfs &distMap(DistMap &m)
@@ -407,20 +409,18 @@
   public:
 
-    ///\name Execution control
-    ///The simplest way to execute the algorithm is to use
-    ///one of the member functions called \ref lemon::Dfs::run() "run()".
-    ///\n
-    ///If you need more control on the execution, first you must call
-    ///\ref lemon::Dfs::init() "init()", then you can add a source node
-    ///with \ref lemon::Dfs::addSource() "addSource()".
-    ///Finally \ref lemon::Dfs::start() "start()" will perform the
-    ///actual path computation.
+    ///\name Execution Control
+    ///The simplest way to execute the DFS algorithm is to use one of the
+    ///member functions called \ref run(Node) "run()".\n
+    ///If you need more control on the execution, first you have to call
+    ///\ref init(), then you can add a source node with \ref addSource()
+    ///and perform the actual computation with \ref start().
+    ///This procedure can be repeated if there are nodes that have not
+    ///been reached.
 
     ///@{
 
+    ///\brief Initializes the internal data structures.
+    ///
     ///Initializes the internal data structures.
-
-    ///Initializes the internal data structures.
-    ///
     void init()
     {
@@ -439,9 +439,8 @@
     ///Adds a new source node to the set of nodes to be processed.
     ///
-    ///\pre The stack must be empty. (Otherwise the algorithm gives
-    ///false results.)
-    ///
-    ///\warning Distances will be wrong (or at least strange) in case of
-    ///multiple sources.
+    ///\pre The stack must be empty. Otherwise the algorithm gives
+    ///wrong results. (One of the outgoing arcs of all the source nodes
+    ///except for the last one will not be visited and distances will
+    ///also be wrong.)
     void addSource(Node s)
     {
@@ -507,14 +506,14 @@
     }
 
-    ///\brief Returns \c false if there are nodes
-    ///to be processed.
-    ///
-    ///Returns \c false if there are nodes
-    ///to be processed in the queue (stack).
+    ///Returns \c false if there are nodes to be processed.
+
+    ///Returns \c false if there are nodes to be processed
+    ///in the queue (stack).
     bool emptyQueue() const { return _stack_head<0; }
 
     ///Returns the number of the nodes to be processed.
 
-    ///Returns the number of the nodes to be processed in the queue (stack).
+    ///Returns the number of the nodes to be processed
+    ///in the queue (stack).
     int queueSize() const { return _stack_head+1; }
 
@@ -638,6 +637,6 @@
     ///
     ///The algorithm computes
-    ///- the %DFS tree,
-    ///- the distance of each node from the root in the %DFS tree.
+    ///- the %DFS tree (forest),
+    ///- the distance of each node from the root(s) in the %DFS tree.
     ///
     ///\note <tt>d.run()</tt> is just a shortcut of the following code.
@@ -664,8 +663,8 @@
 
     ///\name Query Functions
-    ///The result of the %DFS algorithm can be obtained using these
+    ///The results of the DFS algorithm can be obtained using these
     ///functions.\n
-    ///Either \ref lemon::Dfs::run() "run()" or \ref lemon::Dfs::start()
-    ///"start()" must be called before using them.
+    ///Either \ref run(Node) "run()" or \ref start() should be called
+    ///before using them.
 
     ///@{
@@ -675,19 +674,19 @@
     ///Returns the DFS path to a node.
     ///
-    ///\warning \c t should be reachable from the root.
-    ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\warning \c t should be reached from the root(s).
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Path path(Node t) const { return Path(*G, *_pred, t); }
 
-    ///The distance of a node from the root.
-
-    ///Returns the distance of a node from the root.
-    ///
-    ///\warning If node \c v is not reachable from the root, then
+    ///The distance of a node from the root(s).
+
+    ///Returns the distance of a node from the root(s).
+    ///
+    ///\warning If node \c v is not reached from the root(s), then
     ///the return value of this function is undefined.
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     int dist(Node v) const { return (*_dist)[v]; }
 
@@ -695,13 +694,13 @@
 
     ///This function returns the 'previous arc' of the %DFS tree for the
-    ///node \c v, i.e. it returns the last arc of a %DFS path from the
-    ///root to \c v. It is \c INVALID
-    ///if \c v is not reachable from the root(s) or if \c v is a root.
+    ///node \c v, i.e. it returns the last arc of a %DFS path from a
+    ///root to \c v. It is \c INVALID if \c v is not reached from the
+    ///root(s) or if \c v is a root.
     ///
     ///The %DFS tree used here is equal to the %DFS tree used in
     ///\ref predNode().
     ///
-    ///\pre Either \ref run() or \ref start() must be called before using
-    ///this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Arc predArc(Node v) const { return (*_pred)[v];}
 
@@ -710,12 +709,12 @@
     ///This function returns the 'previous node' of the %DFS
     ///tree for the node \c v, i.e. it returns the last but one node
-    ///from a %DFS path from the root to \c v. It is \c INVALID
-    ///if \c v is not reachable from the root(s) or if \c v is a root.
+    ///from a %DFS path from a root to \c v. It is \c INVALID
+    ///if \c v is not reached from the root(s) or if \c v is a root.
     ///
     ///The %DFS tree used here is equal to the %DFS tree used in
     ///\ref predArc().
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Node predNode(Node v) const { return (*_pred)[v]==INVALID ? INVALID:
                                   G->source((*_pred)[v]); }
@@ -727,5 +726,5 @@
     ///distances of the nodes calculated by the algorithm.
     ///
-    ///\pre Either \ref run() or \ref init()
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     const DistMap &distMap() const { return *_dist;}
@@ -737,12 +736,13 @@
     ///arcs, which form the DFS tree.
     ///
-    ///\pre Either \ref run() or \ref init()
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     const PredMap &predMap() const { return *_pred;}
 
-    ///Checks if a node is reachable from the root(s).
-
-    ///Returns \c true if \c v is reachable from the root(s).
-    ///\pre Either \ref run() or \ref start()
+    ///Checks if a node is reached from the root(s).
+
+    ///Returns \c true if \c v is reached from the root(s).
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     bool reached(Node v) const { return (*_reached)[v]; }
@@ -890,6 +890,6 @@
   /// This auxiliary class is created to implement the
   /// \ref dfs() "function-type interface" of \ref Dfs algorithm.
-  /// It does not have own \ref run() method, it uses the functions
-  /// and features of the plain \ref Dfs.
+  /// It does not have own \ref run(Node) "run()" method, it uses the
+  /// functions and features of the plain \ref Dfs.
   ///
   /// This class should only be used through the \ref dfs() function,
@@ -1111,6 +1111,5 @@
   ///  bool reached = dfs(g).path(p).dist(d).run(s,t);
   ///\endcode
-
-  ///\warning Don't forget to put the \ref DfsWizard::run() "run()"
+  ///\warning Don't forget to put the \ref DfsWizard::run(Node) "run()"
   ///to the end of the parameter list.
   ///\sa DfsWizard
@@ -1310,5 +1309,5 @@
     typedef DfsVisit Create;
 
-    /// \name Named template parameters
+    /// \name Named Template Parameters
 
     ///@{
@@ -1352,7 +1351,8 @@
     ///
     /// Sets the map that indicates which nodes are reached.
-    /// If you don't use this function before calling \ref run(),
-    /// it will allocate one. The destructor deallocates this
-    /// automatically allocated map, of course.
+    /// If you don't use this function before calling \ref run(Node) "run()"
+    /// or \ref init(), an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
     /// \return <tt> (*this) </tt>
     DfsVisit &reachedMap(ReachedMap &m) {
@@ -1367,14 +1367,12 @@
   public:
 
-    /// \name Execution control
-    /// The simplest way to execute the algorithm is to use
-    /// one of the member functions called \ref lemon::DfsVisit::run()
-    /// "run()".
-    /// \n
-    /// If you need more control on the execution, first you must call
-    /// \ref lemon::DfsVisit::init() "init()", then you can add several
-    /// source nodes with \ref lemon::DfsVisit::addSource() "addSource()".
-    /// Finally \ref lemon::DfsVisit::start() "start()" will perform the
-    /// actual path computation.
+    /// \name Execution Control
+    /// The simplest way to execute the DFS algorithm is to use one of the
+    /// member functions called \ref run(Node) "run()".\n
+    /// If you need more control on the execution, first you have to call
+    /// \ref init(), then you can add a source node with \ref addSource()
+    /// and perform the actual computation with \ref start().
+    /// This procedure can be repeated if there are nodes that have not
+    /// been reached.
 
     /// @{
@@ -1392,13 +1390,12 @@
     }
 
-    ///Adds a new source node.
-
-    ///Adds a new source node to the set of nodes to be processed.
-    ///
-    ///\pre The stack must be empty. (Otherwise the algorithm gives
-    ///false results.)
-    ///
-    ///\warning Distances will be wrong (or at least strange) in case of
-    ///multiple sources.
+    /// \brief Adds a new source node.
+    ///
+    /// Adds a new source node to the set of nodes to be processed.
+    ///
+    /// \pre The stack must be empty. Otherwise the algorithm gives
+    /// wrong results. (One of the outgoing arcs of all the source nodes
+    /// except for the last one will not be visited and distances will
+    /// also be wrong.)
     void addSource(Node s)
     {
@@ -1590,6 +1587,6 @@
     ///
     /// The algorithm computes
-    /// - the %DFS tree,
-    /// - the distance of each node from the root in the %DFS tree.
+    /// - the %DFS tree (forest),
+    /// - the distance of each node from the root(s) in the %DFS tree.
     ///
     /// \note <tt>d.run()</tt> is just a shortcut of the following code.
@@ -1616,15 +1613,16 @@
 
     /// \name Query Functions
-    /// The result of the %DFS algorithm can be obtained using these
+    /// The results of the DFS algorithm can be obtained using these
     /// functions.\n
-    /// Either \ref lemon::DfsVisit::run() "run()" or
-    /// \ref lemon::DfsVisit::start() "start()" must be called before
-    /// using them.
+    /// Either \ref run(Node) "run()" or \ref start() should be called
+    /// before using them.
+
     ///@{
 
-    /// \brief Checks if a node is reachable from the root(s).
-    ///
-    /// Returns \c true if \c v is reachable from the root(s).
-    /// \pre Either \ref run() or \ref start()
+    /// \brief Checks if a node is reached from the root(s).
+    ///
+    /// Returns \c true if \c v is reached from the root(s).
+    ///
+    /// \pre Either \ref run(Node) "run()" or \ref init()
     /// must be called before using this function.
     bool reached(Node v) { return (*_reached)[v]; }
Index: lemon/dijkstra.h
===================================================================
--- lemon/dijkstra.h	(revision 313)
+++ lemon/dijkstra.h	(revision 424)
@@ -55,27 +55,4 @@
   };
 
-  /// \brief Widest path operation traits for the Dijkstra algorithm class.
-  ///
-  /// This operation traits class defines all computational operations and
-  /// constants which are used in the Dijkstra algorithm for widest path
-  /// computation.
-  ///
-  /// \see DijkstraDefaultOperationTraits
-  template <typename Value>
-  struct DijkstraWidestPathOperationTraits {
-    /// \brief Gives back the maximum value of the type.
-    static Value zero() {
-      return std::numeric_limits<Value>::max();
-    }
-    /// \brief Gives back the minimum of the given two elements.
-    static Value plus(const Value& left, const Value& right) {
-      return std::min(left, right);
-    }
-    /// \brief Gives back true only if the first value is less than the second.
-    static bool less(const Value& left, const Value& right) {
-      return left < right;
-    }
-  };
-
   ///Default traits class of Dijkstra class.
 
@@ -203,18 +180,11 @@
   ///
   ///\tparam GR The type of the digraph the algorithm runs on.
-  ///The default value is \ref ListDigraph.
-  ///The value of GR is not used directly by \ref Dijkstra, it is only
-  ///passed to \ref DijkstraDefaultTraits.
-  ///\tparam LM A readable arc map that determines the lengths of the
-  ///arcs. It is read once for each arc, so the map may involve in
+  ///The default type is \ref ListDigraph.
+  ///\tparam LM A \ref concepts::ReadMap "readable" arc map that specifies
+  ///the lengths of the arcs.
+  ///It is read once for each arc, so the map may involve in
   ///relatively time consuming process to compute the arc lengths if
   ///it is necessary. The default map type is \ref
-  ///concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
-  ///The value of LM is not used directly by \ref Dijkstra, it is only
-  ///passed to \ref DijkstraDefaultTraits.
-  ///\tparam TR Traits class to set various data types used by the algorithm.
-  ///The default traits class is \ref DijkstraDefaultTraits
-  ///"DijkstraDefaultTraits<GR,LM>". See \ref DijkstraDefaultTraits
-  ///for the documentation of a Dijkstra traits class.
+  ///concepts::Digraph::ArcMap "GR::ArcMap<int>".
 #ifdef DOXYGEN
   template <typename GR, typename LM, typename TR>
@@ -250,5 +220,5 @@
     typedef typename TR::OperationTraits OperationTraits;
 
-    ///The traits class.
+    ///The \ref DijkstraDefaultTraits "traits class" of the algorithm.
     typedef TR Traits;
 
@@ -332,4 +302,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///PredMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetPredMap
@@ -352,4 +323,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///DistMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetDistMap
@@ -372,4 +344,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///ProcessedMap type.
+    ///It must meet the \ref concepts::WriteMap "WriteMap" concept.
     template <class T>
     struct SetProcessedMap
@@ -412,8 +385,12 @@
     };
     ///\brief \ref named-templ-param "Named parameter" for setting
-    ///heap and cross reference type
+    ///heap and cross reference types
     ///
     ///\ref named-templ-param "Named parameter" for setting heap and cross
-    ///reference type.
+    ///reference types. If this named parameter is used, then external
+    ///heap and cross reference objects must be passed to the algorithm
+    ///using the \ref heap() function before calling \ref run(Node) "run()"
+    ///or \ref init().
+    ///\sa SetStandardHeap
     template <class H, class CR = typename Digraph::template NodeMap<int> >
     struct SetHeap
@@ -435,10 +412,16 @@
     };
     ///\brief \ref named-templ-param "Named parameter" for setting
-    ///heap and cross reference type with automatic allocation
+    ///heap and cross reference types with automatic allocation
     ///
     ///\ref named-templ-param "Named parameter" for setting heap and cross
-    ///reference type. It can allocate the heap and the cross reference
-    ///object if the cross reference's constructor waits for the digraph as
-    ///parameter and the heap's constructor waits for the cross reference.
+    ///reference types with automatic allocation.
+    ///They should have standard constructor interfaces to be able to
+    ///automatically created by the algorithm (i.e. the digraph should be
+    ///passed to the constructor of the cross reference and the cross
+    ///reference should be passed to the constructor of the heap).
+    ///However external heap and cross reference objects could also be
+    ///passed to the algorithm using the \ref heap() function before
+    ///calling \ref run(Node) "run()" or \ref init().
+    ///\sa SetHeap
     template <class H, class CR = typename Digraph::template NodeMap<int> >
     struct SetStandardHeap
@@ -510,7 +493,8 @@
 
     ///Sets the map that stores the predecessor arcs.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dijkstra &predMap(PredMap &m)
@@ -527,7 +511,8 @@
 
     ///Sets the map that indicates which nodes are processed.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dijkstra &processedMap(ProcessedMap &m)
@@ -545,7 +530,8 @@
     ///Sets the map that stores the distances of the nodes calculated by the
     ///algorithm.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated map, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), an instance will be allocated automatically.
+    ///The destructor deallocates this automatically allocated map,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dijkstra &distMap(DistMap &m)
@@ -562,7 +548,9 @@
 
     ///Sets the heap and the cross reference used by algorithm.
-    ///If you don't use this function before calling \ref run(),
-    ///it will allocate one. The destructor deallocates this
-    ///automatically allocated heap and cross reference, of course.
+    ///If you don't use this function before calling \ref run(Node) "run()"
+    ///or \ref init(), heap and cross reference instances will be
+    ///allocated automatically.
+    ///The destructor deallocates these automatically allocated objects,
+    ///of course.
     ///\return <tt> (*this) </tt>
     Dijkstra &heap(Heap& hp, HeapCrossRef &cr)
@@ -591,20 +579,17 @@
   public:
 
-    ///\name Execution control
-    ///The simplest way to execute the algorithm is to use one of the
-    ///member functions called \ref lemon::Dijkstra::run() "run()".
-    ///\n
-    ///If you need more control on the execution, first you must call
-    ///\ref lemon::Dijkstra::init() "init()", then you can add several
-    ///source nodes with \ref lemon::Dijkstra::addSource() "addSource()".
-    ///Finally \ref lemon::Dijkstra::start() "start()" will perform the
-    ///actual path computation.
+    ///\name Execution Control
+    ///The simplest way to execute the %Dijkstra algorithm is to use
+    ///one of the member functions called \ref run(Node) "run()".\n
+    ///If you need more control on the execution, first you have to call
+    ///\ref init(), then you can add several source nodes with
+    ///\ref addSource(). Finally the actual path computation can be
+    ///performed with one of the \ref start() functions.
 
     ///@{
 
+    ///\brief Initializes the internal data structures.
+    ///
     ///Initializes the internal data structures.
-
-    ///Initializes the internal data structures.
-    ///
     void init()
     {
@@ -682,15 +667,14 @@
     }
 
-    ///\brief Returns \c false if there are nodes
-    ///to be processed.
-    ///
-    ///Returns \c false if there are nodes
-    ///to be processed in the priority heap.
+    ///Returns \c false if there are nodes to be processed.
+
+    ///Returns \c false if there are nodes to be processed
+    ///in the priority heap.
     bool emptyQueue() const { return _heap->empty(); }
 
-    ///Returns the number of the nodes to be processed in the priority heap
-
-    ///Returns the number of the nodes to be processed in the priority heap.
-    ///
+    ///Returns the number of the nodes to be processed.
+
+    ///Returns the number of the nodes to be processed
+    ///in the priority heap.
     int queueSize() const { return _heap->size(); }
 
@@ -813,9 +797,8 @@
 
     ///\name Query Functions
-    ///The result of the %Dijkstra algorithm can be obtained using these
+    ///The results of the %Dijkstra algorithm can be obtained using these
     ///functions.\n
-    ///Either \ref lemon::Dijkstra::run() "run()" or
-    ///\ref lemon::Dijkstra::start() "start()" must be called before
-    ///using them.
+    ///Either \ref run(Node) "run()" or \ref start() should be called
+    ///before using them.
 
     ///@{
@@ -825,8 +808,8 @@
     ///Returns the shortest path to a node.
     ///
-    ///\warning \c t should be reachable from the root(s).
-    ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\warning \c t should be reached from the root(s).
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Path path(Node t) const { return Path(*G, *_pred, t); }
 
@@ -835,9 +818,9 @@
     ///Returns the distance of a node from the root(s).
     ///
-    ///\warning If node \c v is not reachable from the root(s), then
+    ///\warning If node \c v is not reached from the root(s), then
     ///the return value of this function is undefined.
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Value dist(Node v) const { return (*_dist)[v]; }
 
@@ -846,12 +829,12 @@
     ///This function returns the 'previous arc' of the shortest path
     ///tree for the node \c v, i.e. it returns the last arc of a
-    ///shortest path from the root(s) to \c v. It is \c INVALID if \c v
-    ///is not reachable from the root(s) or if \c v is a root.
+    ///shortest path from a root to \c v. It is \c INVALID if \c v
+    ///is not reached from the root(s) or if \c v is a root.
     ///
     ///The shortest path tree used here is equal to the shortest path
     ///tree used in \ref predNode().
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Arc predArc(Node v) const { return (*_pred)[v]; }
 
@@ -860,12 +843,12 @@
     ///This function returns the 'previous node' of the shortest path
     ///tree for the node \c v, i.e. it returns the last but one node
-    ///from a shortest path from the root(s) to \c v. It is \c INVALID
-    ///if \c v is not reachable from the root(s) or if \c v is a root.
+    ///from a shortest path from a root to \c v. It is \c INVALID
+    ///if \c v is not reached from the root(s) or if \c v is a root.
     ///
     ///The shortest path tree used here is equal to the shortest path
     ///tree used in \ref predArc().
     ///
-    ///\pre Either \ref run() or \ref start() must be called before
-    ///using this function.
+    ///\pre Either \ref run(Node) "run()" or \ref init()
+    ///must be called before using this function.
     Node predNode(Node v) const { return (*_pred)[v]==INVALID ? INVALID:
                                   G->source((*_pred)[v]); }
@@ -877,5 +860,5 @@
     ///of the nodes calculated by the algorithm.
     ///
-    ///\pre Either \ref run() or \ref init()
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     const DistMap &distMap() const { return *_dist;}
@@ -887,12 +870,13 @@
     ///arcs, which form the shortest path tree.
     ///
-    ///\pre Either \ref run() or \ref init()
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     const PredMap &predMap() const { return *_pred;}
 
-    ///Checks if a node is reachable from the root(s).
-
-    ///Returns \c true if \c v is reachable from the root(s).
-    ///\pre Either \ref run() or \ref start()
+    ///Checks if a node is reached from the root(s).
+
+    ///Returns \c true if \c v is reached from the root(s).
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     bool reached(Node v) const { return (*_heap_cross_ref)[v] !=
@@ -903,5 +887,6 @@
     ///Returns \c true if \c v is processed, i.e. the shortest
     ///path to \c v has already found.
-    ///\pre Either \ref run() or \ref init()
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function.
     bool processed(Node v) const { return (*_heap_cross_ref)[v] ==
@@ -912,5 +897,6 @@
     ///Returns the current distance of a node from the root(s).
     ///It may be decreased in the following processes.
-    ///\pre Either \ref run() or \ref init()
+    ///
+    ///\pre Either \ref run(Node) "run()" or \ref init()
     ///must be called before using this function and
     ///node \c v must be reached but not necessarily processed.
@@ -1095,6 +1081,6 @@
   /// This auxiliary class is created to implement the
   /// \ref dijkstra() "function-type interface" of \ref Dijkstra algorithm.
-  /// It does not have own \ref run() method, it uses the functions
-  /// and features of the plain \ref Dijkstra.
+  /// It does not have own \ref run(Node) "run()" method, it uses the
+  /// functions and features of the plain \ref Dijkstra.
   ///
   /// This class should only be used through the \ref dijkstra() function,
@@ -1291,5 +1277,5 @@
   ///  bool reached = dijkstra(g,length).path(p).dist(d).run(s,t);
   ///\endcode
-  ///\warning Don't forget to put the \ref DijkstraWizard::run() "run()"
+  ///\warning Don't forget to put the \ref DijkstraWizard::run(Node) "run()"
   ///to the end of the parameter list.
   ///\sa DijkstraWizard
Index: lemon/hao_orlin.h
===================================================================
--- lemon/hao_orlin.h	(revision 428)
+++ lemon/hao_orlin.h	(revision 428)
@@ -0,0 +1,966 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * 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_HAO_ORLIN_H
+#define LEMON_HAO_ORLIN_H
+
+#include <vector>
+#include <list>
+#include <limits>
+
+#include <lemon/maps.h>
+#include <lemon/core.h>
+#include <lemon/tolerance.h>
+
+/// \file
+/// \ingroup min_cut
+/// \brief Implementation of the Hao-Orlin algorithm.
+///
+/// Implementation of the Hao-Orlin algorithm class for testing network
+/// reliability.
+
+namespace lemon {
+
+  /// \ingroup min_cut
+  ///
+  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
+  ///
+  /// Hao-Orlin calculates a minimum cut in a directed graph
+  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
+  /// consists of two phases: in the first phase it determines a
+  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
+  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
+  /// out-degree) and in the second phase it determines a minimum cut
+  /// with \f$ source \f$ on the sink-side (i.e. a set
+  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
+  /// out-degree). Obviously, the smaller of these two cuts will be a
+  /// minimum cut of \f$ D \f$. The algorithm is a modified
+  /// push-relabel preflow algorithm and our implementation calculates
+  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
+  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
+  /// purpose of such algorithm is testing network reliability. For an
+  /// undirected graph you can run just the first phase of the
+  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
+  /// which solves the undirected problem in
+  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
+  /// NagamochiIbaraki algorithm class.
+  ///
+  /// \param _Digraph is the graph type of the algorithm.
+  /// \param _CapacityMap is an edge map of capacities which should
+  /// be any numreric type. The default type is _Digraph::ArcMap<int>.
+  /// \param _Tolerance is the handler of the inexact computation. The
+  /// default type for this is Tolerance<CapacityMap::Value>.
+#ifdef DOXYGEN
+  template <typename _Digraph, typename _CapacityMap, typename _Tolerance>
+#else
+  template <typename _Digraph,
+            typename _CapacityMap = typename _Digraph::template ArcMap<int>,
+            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
+#endif
+  class HaoOrlin {
+  private:
+
+    typedef _Digraph Digraph;
+    typedef _CapacityMap CapacityMap;
+    typedef _Tolerance Tolerance;
+
+    typedef typename CapacityMap::Value Value;
+
+    TEMPLATE_GRAPH_TYPEDEFS(Digraph);
+
+    const Digraph& _graph;
+    const CapacityMap* _capacity;
+
+    typedef typename Digraph::template ArcMap<Value> FlowMap;
+    FlowMap* _flow;
+
+    Node _source;
+
+    int _node_num;
+
+    // Bucketing structure
+    std::vector<Node> _first, _last;
+    typename Digraph::template NodeMap<Node>* _next;
+    typename Digraph::template NodeMap<Node>* _prev;
+    typename Digraph::template NodeMap<bool>* _active;
+    typename Digraph::template NodeMap<int>* _bucket;
+
+    std::vector<bool> _dormant;
+
+    std::list<std::list<int> > _sets;
+    std::list<int>::iterator _highest;
+
+    typedef typename Digraph::template NodeMap<Value> ExcessMap;
+    ExcessMap* _excess;
+
+    typedef typename Digraph::template NodeMap<bool> SourceSetMap;
+    SourceSetMap* _source_set;
+
+    Value _min_cut;
+
+    typedef typename Digraph::template NodeMap<bool> MinCutMap;
+    MinCutMap* _min_cut_map;
+
+    Tolerance _tolerance;
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor of the algorithm class.
+    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
+             const Tolerance& tolerance = Tolerance()) :
+      _graph(graph), _capacity(&capacity), _flow(0), _source(),
+      _node_num(), _first(), _last(), _next(0), _prev(0),
+      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
+      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
+      _tolerance(tolerance) {}
+
+    ~HaoOrlin() {
+      if (_min_cut_map) {
+        delete _min_cut_map;
+      }
+      if (_source_set) {
+        delete _source_set;
+      }
+      if (_excess) {
+        delete _excess;
+      }
+      if (_next) {
+        delete _next;
+      }
+      if (_prev) {
+        delete _prev;
+      }
+      if (_active) {
+        delete _active;
+      }
+      if (_bucket) {
+        delete _bucket;
+      }
+      if (_flow) {
+        delete _flow;
+      }
+    }
+
+  private:
+
+    void activate(const Node& i) {
+      _active->set(i, true);
+
+      int bucket = (*_bucket)[i];
+
+      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
+      //unlace
+      _next->set((*_prev)[i], (*_next)[i]);
+      if ((*_next)[i] != INVALID) {
+        _prev->set((*_next)[i], (*_prev)[i]);
+      } else {
+        _last[bucket] = (*_prev)[i];
+      }
+      //lace
+      _next->set(i, _first[bucket]);
+      _prev->set(_first[bucket], i);
+      _prev->set(i, INVALID);
+      _first[bucket] = i;
+    }
+
+    void deactivate(const Node& i) {
+      _active->set(i, false);
+      int bucket = (*_bucket)[i];
+
+      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
+
+      //unlace
+      _prev->set((*_next)[i], (*_prev)[i]);
+      if ((*_prev)[i] != INVALID) {
+        _next->set((*_prev)[i], (*_next)[i]);
+      } else {
+        _first[bucket] = (*_next)[i];
+      }
+      //lace
+      _prev->set(i, _last[bucket]);
+      _next->set(_last[bucket], i);
+      _next->set(i, INVALID);
+      _last[bucket] = i;
+    }
+
+    void addItem(const Node& i, int bucket) {
+      (*_bucket)[i] = bucket;
+      if (_last[bucket] != INVALID) {
+        _prev->set(i, _last[bucket]);
+        _next->set(_last[bucket], i);
+        _next->set(i, INVALID);
+        _last[bucket] = i;
+      } else {
+        _prev->set(i, INVALID);
+        _first[bucket] = i;
+        _next->set(i, INVALID);
+        _last[bucket] = i;
+      }
+    }
+
+    void findMinCutOut() {
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _excess->set(n, 0);
+      }
+
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _flow->set(a, 0);
+      }
+
+      int bucket_num = 0;
+      std::vector<Node> queue(_node_num);
+      int qfirst = 0, qlast = 0, qsep = 0;
+
+      {
+        typename Digraph::template NodeMap<bool> reached(_graph, false);
+
+        reached.set(_source, true);
+        bool first_set = true;
+
+        for (NodeIt t(_graph); t != INVALID; ++t) {
+          if (reached[t]) continue;
+          _sets.push_front(std::list<int>());
+          
+          queue[qlast++] = t;
+          reached.set(t, true);
+
+          while (qfirst != qlast) {
+            if (qsep == qfirst) {
+              ++bucket_num;
+              _sets.front().push_front(bucket_num);
+              _dormant[bucket_num] = !first_set;
+              _first[bucket_num] = _last[bucket_num] = INVALID;
+              qsep = qlast;
+            }
+
+            Node n = queue[qfirst++];
+            addItem(n, bucket_num);
+
+            for (InArcIt a(_graph, n); a != INVALID; ++a) {
+              Node u = _graph.source(a);
+              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
+                reached.set(u, true);
+                queue[qlast++] = u;
+              }
+            }
+          }
+          first_set = false;
+        }
+
+        ++bucket_num;
+        _bucket->set(_source, 0);
+        _dormant[0] = true;
+      }
+      _source_set->set(_source, true);
+
+      Node target = _last[_sets.back().back()];
+      {
+        for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
+          if (_tolerance.positive((*_capacity)[a])) {
+            Node u = _graph.target(a);
+            _flow->set(a, (*_capacity)[a]);
+            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
+            if (!(*_active)[u] && u != _source) {
+              activate(u);
+            }
+          }
+        }
+
+        if ((*_active)[target]) {
+          deactivate(target);
+        }
+
+        _highest = _sets.back().begin();
+        while (_highest != _sets.back().end() &&
+               !(*_active)[_first[*_highest]]) {
+          ++_highest;
+        }
+      }
+
+      while (true) {
+        while (_highest != _sets.back().end()) {
+          Node n = _first[*_highest];
+          Value excess = (*_excess)[n];
+          int next_bucket = _node_num;
+
+          int under_bucket;
+          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
+            under_bucket = -1;
+          } else {
+            under_bucket = *(++std::list<int>::iterator(_highest));
+          }
+
+          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
+            Node v = _graph.target(a);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_capacity)[a] - (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(a, (*_flow)[a] + excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(a, (*_capacity)[a]);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+
+          for (InArcIt a(_graph, n); a != INVALID; ++a) {
+            Node v = _graph.source(a);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(a, (*_flow)[a] - excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(a, 0);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+
+        no_more_push:
+
+          _excess->set(n, excess);
+
+          if (excess != 0) {
+            if ((*_next)[n] == INVALID) {
+              typename std::list<std::list<int> >::iterator new_set =
+                _sets.insert(--_sets.end(), std::list<int>());
+              new_set->splice(new_set->end(), _sets.back(),
+                              _sets.back().begin(), ++_highest);
+              for (std::list<int>::iterator it = new_set->begin();
+                   it != new_set->end(); ++it) {
+                _dormant[*it] = true;
+              }
+              while (_highest != _sets.back().end() &&
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else if (next_bucket == _node_num) {
+              _first[(*_bucket)[n]] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+
+              std::list<std::list<int> >::iterator new_set =
+                _sets.insert(--_sets.end(), std::list<int>());
+
+              new_set->push_front(bucket_num);
+              _bucket->set(n, bucket_num);
+              _first[bucket_num] = _last[bucket_num] = n;
+              _next->set(n, INVALID);
+              _prev->set(n, INVALID);
+              _dormant[bucket_num] = true;
+              ++bucket_num;
+
+              while (_highest != _sets.back().end() &&
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else {
+              _first[*_highest] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+
+              while (next_bucket != *_highest) {
+                --_highest;
+              }
+
+              if (_highest == _sets.back().begin()) {
+                _sets.back().push_front(bucket_num);
+                _dormant[bucket_num] = false;
+                _first[bucket_num] = _last[bucket_num] = INVALID;
+                ++bucket_num;
+              }
+              --_highest;
+
+              _bucket->set(n, *_highest);
+              _next->set(n, _first[*_highest]);
+              if (_first[*_highest] != INVALID) {
+                _prev->set(_first[*_highest], n);
+              } else {
+                _last[*_highest] = n;
+              }
+              _first[*_highest] = n;
+            }
+          } else {
+
+            deactivate(n);
+            if (!(*_active)[_first[*_highest]]) {
+              ++_highest;
+              if (_highest != _sets.back().end() &&
+                  !(*_active)[_first[*_highest]]) {
+                _highest = _sets.back().end();
+              }
+            }
+          }
+        }
+
+        if ((*_excess)[target] < _min_cut) {
+          _min_cut = (*_excess)[target];
+          for (NodeIt i(_graph); i != INVALID; ++i) {
+            _min_cut_map->set(i, true);
+          }
+          for (std::list<int>::iterator it = _sets.back().begin();
+               it != _sets.back().end(); ++it) {
+            Node n = _first[*it];
+            while (n != INVALID) {
+              _min_cut_map->set(n, false);
+              n = (*_next)[n];
+            }
+          }
+        }
+
+        {
+          Node new_target;
+          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
+            if ((*_next)[target] == INVALID) {
+              _last[(*_bucket)[target]] = (*_prev)[target];
+              new_target = (*_prev)[target];
+            } else {
+              _prev->set((*_next)[target], (*_prev)[target]);
+              new_target = (*_next)[target];
+            }
+            if ((*_prev)[target] == INVALID) {
+              _first[(*_bucket)[target]] = (*_next)[target];
+            } else {
+              _next->set((*_prev)[target], (*_next)[target]);
+            }
+          } else {
+            _sets.back().pop_back();
+            if (_sets.back().empty()) {
+              _sets.pop_back();
+              if (_sets.empty())
+                break;
+              for (std::list<int>::iterator it = _sets.back().begin();
+                   it != _sets.back().end(); ++it) {
+                _dormant[*it] = false;
+              }
+            }
+            new_target = _last[_sets.back().back()];
+          }
+
+          _bucket->set(target, 0);
+
+          _source_set->set(target, true);
+          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
+            Value rem = (*_capacity)[a] - (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.target(a);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(a, (*_capacity)[a]);
+          }
+
+          for (InArcIt a(_graph, target); a != INVALID; ++a) {
+            Value rem = (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.source(a);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(a, 0);
+          }
+
+          target = new_target;
+          if ((*_active)[target]) {
+            deactivate(target);
+          }
+
+          _highest = _sets.back().begin();
+          while (_highest != _sets.back().end() &&
+                 !(*_active)[_first[*_highest]]) {
+            ++_highest;
+          }
+        }
+      }
+    }
+
+    void findMinCutIn() {
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _excess->set(n, 0);
+      }
+
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _flow->set(a, 0);
+      }
+
+      int bucket_num = 0;
+      std::vector<Node> queue(_node_num);
+      int qfirst = 0, qlast = 0, qsep = 0;
+
+      {
+        typename Digraph::template NodeMap<bool> reached(_graph, false);
+
+        reached.set(_source, true);
+
+        bool first_set = true;
+
+        for (NodeIt t(_graph); t != INVALID; ++t) {
+          if (reached[t]) continue;
+          _sets.push_front(std::list<int>());
+          
+          queue[qlast++] = t;
+          reached.set(t, true);
+
+          while (qfirst != qlast) {
+            if (qsep == qfirst) {
+              ++bucket_num;
+              _sets.front().push_front(bucket_num);
+              _dormant[bucket_num] = !first_set;
+              _first[bucket_num] = _last[bucket_num] = INVALID;
+              qsep = qlast;
+            }
+
+            Node n = queue[qfirst++];
+            addItem(n, bucket_num);
+
+            for (OutArcIt a(_graph, n); a != INVALID; ++a) {
+              Node u = _graph.target(a);
+              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
+                reached.set(u, true);
+                queue[qlast++] = u;
+              }
+            }
+          }
+          first_set = false;
+        }
+
+        ++bucket_num;
+        _bucket->set(_source, 0);
+        _dormant[0] = true;
+      }
+      _source_set->set(_source, true);
+
+      Node target = _last[_sets.back().back()];
+      {
+        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
+          if (_tolerance.positive((*_capacity)[a])) {
+            Node u = _graph.source(a);
+            _flow->set(a, (*_capacity)[a]);
+            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
+            if (!(*_active)[u] && u != _source) {
+              activate(u);
+            }
+          }
+        }
+        if ((*_active)[target]) {
+          deactivate(target);
+        }
+
+        _highest = _sets.back().begin();
+        while (_highest != _sets.back().end() &&
+               !(*_active)[_first[*_highest]]) {
+          ++_highest;
+        }
+      }
+
+
+      while (true) {
+        while (_highest != _sets.back().end()) {
+          Node n = _first[*_highest];
+          Value excess = (*_excess)[n];
+          int next_bucket = _node_num;
+
+          int under_bucket;
+          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
+            under_bucket = -1;
+          } else {
+            under_bucket = *(++std::list<int>::iterator(_highest));
+          }
+
+          for (InArcIt a(_graph, n); a != INVALID; ++a) {
+            Node v = _graph.source(a);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_capacity)[a] - (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(a, (*_flow)[a] + excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(a, (*_capacity)[a]);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+
+          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
+            Node v = _graph.target(a);
+            if (_dormant[(*_bucket)[v]]) continue;
+            Value rem = (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            if ((*_bucket)[v] == under_bucket) {
+              if (!(*_active)[v] && v != target) {
+                activate(v);
+              }
+              if (!_tolerance.less(rem, excess)) {
+                _flow->set(a, (*_flow)[a] - excess);
+                _excess->set(v, (*_excess)[v] + excess);
+                excess = 0;
+                goto no_more_push;
+              } else {
+                excess -= rem;
+                _excess->set(v, (*_excess)[v] + rem);
+                _flow->set(a, 0);
+              }
+            } else if (next_bucket > (*_bucket)[v]) {
+              next_bucket = (*_bucket)[v];
+            }
+          }
+
+        no_more_push:
+
+          _excess->set(n, excess);
+
+          if (excess != 0) {
+            if ((*_next)[n] == INVALID) {
+              typename std::list<std::list<int> >::iterator new_set =
+                _sets.insert(--_sets.end(), std::list<int>());
+              new_set->splice(new_set->end(), _sets.back(),
+                              _sets.back().begin(), ++_highest);
+              for (std::list<int>::iterator it = new_set->begin();
+                   it != new_set->end(); ++it) {
+                _dormant[*it] = true;
+              }
+              while (_highest != _sets.back().end() &&
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else if (next_bucket == _node_num) {
+              _first[(*_bucket)[n]] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+
+              std::list<std::list<int> >::iterator new_set =
+                _sets.insert(--_sets.end(), std::list<int>());
+
+              new_set->push_front(bucket_num);
+              _bucket->set(n, bucket_num);
+              _first[bucket_num] = _last[bucket_num] = n;
+              _next->set(n, INVALID);
+              _prev->set(n, INVALID);
+              _dormant[bucket_num] = true;
+              ++bucket_num;
+
+              while (_highest != _sets.back().end() &&
+                     !(*_active)[_first[*_highest]]) {
+                ++_highest;
+              }
+            } else {
+              _first[*_highest] = (*_next)[n];
+              _prev->set((*_next)[n], INVALID);
+
+              while (next_bucket != *_highest) {
+                --_highest;
+              }
+              if (_highest == _sets.back().begin()) {
+                _sets.back().push_front(bucket_num);
+                _dormant[bucket_num] = false;
+                _first[bucket_num] = _last[bucket_num] = INVALID;
+                ++bucket_num;
+              }
+              --_highest;
+
+              _bucket->set(n, *_highest);
+              _next->set(n, _first[*_highest]);
+              if (_first[*_highest] != INVALID) {
+                _prev->set(_first[*_highest], n);
+              } else {
+                _last[*_highest] = n;
+              }
+              _first[*_highest] = n;
+            }
+          } else {
+
+            deactivate(n);
+            if (!(*_active)[_first[*_highest]]) {
+              ++_highest;
+              if (_highest != _sets.back().end() &&
+                  !(*_active)[_first[*_highest]]) {
+                _highest = _sets.back().end();
+              }
+            }
+          }
+        }
+
+        if ((*_excess)[target] < _min_cut) {
+          _min_cut = (*_excess)[target];
+          for (NodeIt i(_graph); i != INVALID; ++i) {
+            _min_cut_map->set(i, false);
+          }
+          for (std::list<int>::iterator it = _sets.back().begin();
+               it != _sets.back().end(); ++it) {
+            Node n = _first[*it];
+            while (n != INVALID) {
+              _min_cut_map->set(n, true);
+              n = (*_next)[n];
+            }
+          }
+        }
+
+        {
+          Node new_target;
+          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
+            if ((*_next)[target] == INVALID) {
+              _last[(*_bucket)[target]] = (*_prev)[target];
+              new_target = (*_prev)[target];
+            } else {
+              _prev->set((*_next)[target], (*_prev)[target]);
+              new_target = (*_next)[target];
+            }
+            if ((*_prev)[target] == INVALID) {
+              _first[(*_bucket)[target]] = (*_next)[target];
+            } else {
+              _next->set((*_prev)[target], (*_next)[target]);
+            }
+          } else {
+            _sets.back().pop_back();
+            if (_sets.back().empty()) {
+              _sets.pop_back();
+              if (_sets.empty())
+                break;
+              for (std::list<int>::iterator it = _sets.back().begin();
+                   it != _sets.back().end(); ++it) {
+                _dormant[*it] = false;
+              }
+            }
+            new_target = _last[_sets.back().back()];
+          }
+
+          _bucket->set(target, 0);
+
+          _source_set->set(target, true);
+          for (InArcIt a(_graph, target); a != INVALID; ++a) {
+            Value rem = (*_capacity)[a] - (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.source(a);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(a, (*_capacity)[a]);
+          }
+
+          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
+            Value rem = (*_flow)[a];
+            if (!_tolerance.positive(rem)) continue;
+            Node v = _graph.target(a);
+            if (!(*_active)[v] && !(*_source_set)[v]) {
+              activate(v);
+            }
+            _excess->set(v, (*_excess)[v] + rem);
+            _flow->set(a, 0);
+          }
+
+          target = new_target;
+          if ((*_active)[target]) {
+            deactivate(target);
+          }
+
+          _highest = _sets.back().begin();
+          while (_highest != _sets.back().end() &&
+                 !(*_active)[_first[*_highest]]) {
+            ++_highest;
+          }
+        }
+      }
+    }
+
+  public:
+
+    /// \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 the execution,
+    /// first you must call \ref init(), then the \ref calculateIn() or
+    /// \ref calculateOut() functions.
+
+    /// @{
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures. It creates
+    /// the maps, residual graph adaptors and some bucket structures
+    /// for the algorithm.
+    void init() {
+      init(NodeIt(_graph));
+    }
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures. It creates
+    /// the maps, residual graph adaptor and some bucket structures
+    /// for the algorithm. Node \c source  is used as the push-relabel
+    /// algorithm's source.
+    void init(const Node& source) {
+      _source = source;
+
+      _node_num = countNodes(_graph);
+
+      _first.resize(_node_num);
+      _last.resize(_node_num);
+
+      _dormant.resize(_node_num);
+
+      if (!_flow) {
+        _flow = new FlowMap(_graph);
+      }
+      if (!_next) {
+        _next = new typename Digraph::template NodeMap<Node>(_graph);
+      }
+      if (!_prev) {
+        _prev = new typename Digraph::template NodeMap<Node>(_graph);
+      }
+      if (!_active) {
+        _active = new typename Digraph::template NodeMap<bool>(_graph);
+      }
+      if (!_bucket) {
+        _bucket = new typename Digraph::template NodeMap<int>(_graph);
+      }
+      if (!_excess) {
+        _excess = new ExcessMap(_graph);
+      }
+      if (!_source_set) {
+        _source_set = new SourceSetMap(_graph);
+      }
+      if (!_min_cut_map) {
+        _min_cut_map = new MinCutMap(_graph);
+      }
+
+      _min_cut = std::numeric_limits<Value>::max();
+    }
+
+
+    /// \brief Calculates a minimum cut with \f$ source \f$ on the
+    /// source-side.
+    ///
+    /// Calculates a minimum cut with \f$ source \f$ on the
+    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with
+    /// \f$ source \in X \f$ and minimal out-degree).
+    void calculateOut() {
+      findMinCutOut();
+    }
+
+    /// \brief Calculates a minimum cut with \f$ source \f$ on the
+    /// target-side.
+    ///
+    /// Calculates a minimum cut with \f$ source \f$ on the
+    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with
+    /// \f$ source \in X \f$ and minimal out-degree).
+    void calculateIn() {
+      findMinCutIn();
+    }
+
+
+    /// \brief Runs the algorithm.
+    ///
+    /// Runs the algorithm. It finds nodes \c source and \c target
+    /// arbitrarily and then calls \ref init(), \ref calculateOut()
+    /// and \ref calculateIn().
+    void run() {
+      init();
+      calculateOut();
+      calculateIn();
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// Runs the algorithm. It uses the given \c source node, finds a
+    /// proper \c target and then calls the \ref init(), \ref
+    /// calculateOut() and \ref calculateIn().
+    void run(const Node& s) {
+      init(s);
+      calculateOut();
+      calculateIn();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %HaoOrlin algorithm
+    /// can be obtained using these functions.
+    /// \n
+    /// Before using these functions, either \ref run(), \ref
+    /// calculateOut() or \ref calculateIn() must be called.
+
+    /// @{
+
+    /// \brief Returns the value of the minimum value cut.
+    ///
+    /// Returns the value of the minimum value cut.
+    Value minCutValue() const {
+      return _min_cut;
+    }
+
+
+    /// \brief Returns a minimum cut.
+    ///
+    /// Sets \c nodeMap to the characteristic vector of a minimum
+    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
+    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
+    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
+    /// bool-valued node-map.
+    template <typename NodeMap>
+    Value minCutMap(NodeMap& nodeMap) const {
+      for (NodeIt it(_graph); it != INVALID; ++it) {
+        nodeMap.set(it, (*_min_cut_map)[it]);
+      }
+      return _min_cut;
+    }
+
+    /// @}
+
+  }; //class HaoOrlin
+
+
+} //namespace lemon
+
+#endif //LEMON_HAO_ORLIN_H
