diff -r 10715b6bcd86 -r b7727edd44f2 lemon/circulation.h
--- a/lemon/circulation.h	Wed Nov 28 17:40:41 2007 +0000
+++ b/lemon/circulation.h	Wed Nov 28 17:51:02 2007 +0000
@@ -31,268 +31,612 @@
 ///
 namespace lemon {
 
-  ///Preflow algorithm for the Network Circulation Problem.
+  /// \brief Default traits class of Circulation class.
+  ///
+  /// Default traits class of Circulation class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _LCapMap, 
+	    typename _UCapMap, typename _DeltaMap>
+  struct CirculationDefaultTraits {
+
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \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 upper bound of
+    /// node excess.
+    ///
+    /// The type of the map that stores the lower bound of node
+    /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
+    /// concept.
+    typedef _DeltaMap DeltaMap;
+
+    /// \brief The type of the length of the edges.
+    typedef typename DeltaMap::Value Value;
+
+    /// \brief The map type that stores the flow values.
+    ///
+    /// The map type that stores the flow values. 
+    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param graph The graph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Graph& graph) {
+      return new FlowMap(graph);
+    }
+
+    /// \brief The eleavator type used by Circulation algorithm.
+    /// 
+    /// The elevator type used by Circulation algorithm.
+    ///
+    /// \sa Elevator
+    /// \sa LinkedElevator
+    typedef Elevator<Graph, typename Graph::Node> Elevator;
+    
+    /// \brief Instantiates an Elevator.
+    ///
+    /// This function instantiates a \ref Elevator. 
+    /// \param graph The graph, to which we would like to define the elevator.
+    /// \param max_level The maximum level of the elevator.
+    static Elevator* createElevator(const Graph& graph, int max_level) {
+      return new Elevator(graph, max_level);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+  
+  ///Push-relabel algorithm for the Network Circulation Problem.
   
   ///\ingroup max_flow
-  ///This class implements a preflow algorithm
+  ///This class implements a push-relabel algorithm
   ///for the Network Circulation Problem.
   ///The exact formulation of this problem is the following.
   /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq -delta(v)\quad \forall v\in V \f]
   /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
   ///
-  template<class Graph,
-	   class Value,
-	   class FlowMap=typename Graph::template EdgeMap<Value>,
-	   class LCapMap=typename Graph::template EdgeMap<Value>,
-	   class UCapMap=LCapMap,
-	   class DeltaMap=typename Graph::template NodeMap<Value>
-	    >
+  template<class _Graph,
+	   class _LCapMap=typename _Graph::template EdgeMap<int>,
+	   class _UCapMap=_LCapMap,
+	   class _DeltaMap=typename _Graph::template NodeMap<
+	     typename _UCapMap::Value>,
+	   class _Traits=CirculationDefaultTraits<_Graph, _LCapMap, 
+						  _UCapMap, _DeltaMap> >
   class Circulation {
-    typedef typename Graph::Node Node;
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::Edge Edge;
-    typedef typename Graph::EdgeIt EdgeIt;
-    typedef typename Graph::InEdgeIt InEdgeIt;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
-    
+
+    typedef _Traits Traits;
+    typedef typename Traits::Graph Graph;
+    GRAPH_TYPEDEFS(typename Graph);
+
+    typedef typename Traits::Value Value;
+
+    typedef typename Traits::LCapMap LCapMap;
+    typedef typename Traits::UCapMap UCapMap;
+    typedef typename Traits::DeltaMap DeltaMap;
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Elevator Elevator;
+    typedef typename Traits::Tolerance Tolerance;
+
+    typedef typename Graph::template NodeMap<Value> ExcessMap;
 
     const Graph &_g;
     int _node_num;
 
-    const LCapMap &_lo;
-    const UCapMap &_up;
-    const DeltaMap &_delta;
-    FlowMap &_x;
-    Tolerance<Value> _tol;
-    Elevator<Graph,typename Graph::Node> _levels;
-    typename Graph::template NodeMap<Value> _excess;
+    const LCapMap *_lo;
+    const UCapMap *_up;
+    const DeltaMap *_delta;
+
+    FlowMap *_flow;
+    bool _local_flow;
+
+    Elevator* _level;
+    bool _local_level;
+
+    ExcessMap* _excess;
+
+    Tolerance _tol;
+    int _el;
 
   public:
-    ///\e
-    Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
-		const DeltaMap &delta,FlowMap &x) :
-      _g(g),
-      _node_num(countNodes(g)),
-      _lo(lo),_up(up),_delta(delta),_x(x),
-      _levels(g,_node_num),
-      _excess(g)
-    {
+
+    typedef Circulation Create;
+
+    ///\name Named template parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct DefFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Graph&) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename _FlowMap>
+    struct DefFlowMap 
+      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
+			   DefFlowMapTraits<_FlowMap> > {
+      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
+			  DefFlowMapTraits<_FlowMap> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph&, int) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type
+    template <typename _Elevator>
+    struct DefElevator 
+      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
+			   DefElevatorTraits<_Elevator> > {
+      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
+			  DefElevatorTraits<_Elevator> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefStandardElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph& graph, int max_level) {
+	return new Elevator(graph, max_level);
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type. The Elevator should be standard constructor interface, ie.
+    /// the graph and the maximum level should be passed to it.
+    template <typename _Elevator>
+    struct DefStandardElevator 
+      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
+		       DefStandardElevatorTraits<_Elevator> > {
+      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
+		      DefStandardElevatorTraits<_Elevator> > Create;
+    };    
+
+    /// @}
+
+    /// The constructor of the class.
+
+    /// The constructor of the class. 
+    /// \param g The directed graph the algorithm runs on. 
+    /// \param lo The lower bound capacity of the edges. 
+    /// \param up The upper bound capacity of the edges.
+    /// \param delta The lower bound on node excess.
+    Circulation(const Graph &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() {}
+
+    /// Destrcutor.
+    ~Circulation() {
+      destroyStructures();
     }
-    
+
   private:
 
-    void addExcess(Node n,Value v)
-    {
-      if(_tol.positive(_excess[n]+=v))
-	{
-	  if(!_levels.active(n)) _levels.activate(n);
-	}
-      else if(_levels.active(n)) _levels.deactivate(n);
+    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 init() 
-    {
-     
-      _x=_lo;
 
-      for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=_delta[n];
-
-      for(EdgeIt e(_g);e!=INVALID;++e)
-	{
-	  _excess[_g.target(e)]+=_x[e];
-	  _excess[_g.source(e)]-=_x[e];
-	}
-     
-      _levels.initStart();
-      for(NodeIt n(_g);n!=INVALID;++n)
-	_levels.initAddItem(n);
-      _levels.initFinish();
-      for(NodeIt n(_g);n!=INVALID;++n)
-	if(_tol.positive(_excess[n]))
-	  _levels.activate(n);
+    void destroyStructures() {
+      if (_local_flow) {
+	delete _flow;
+      }
+      if (_local_level) {
+	delete _level;
+      }
+      if (_excess) {
+	delete _excess;
+      }
     }
 
   public:
-    ///Check if \c x is a feasible circulation
-    template<class FT>
-    bool checkX(FT &x) {
+
+    /// Sets the lower bound capacity map.
+
+    /// Sets the lower bound capacity map.
+    /// \return \c (*this)
+    Circulation& lowerCapMap(const LCapMap& map) {
+      _lo = &map;
+      return *this;
+    }
+
+    /// Sets the upper bound capacity map.
+
+    /// Sets the upper bound capacity map.
+    /// \return \c (*this)
+    Circulation& upperCapMap(const LCapMap& map) {
+      _up = &map;
+      return *this;
+    }
+
+    /// Sets the lower bound map on excess.
+
+    /// Sets the lower bound map on excess.
+    /// \return \c (*this)
+    Circulation& deltaMap(const DeltaMap& map) {
+      _delta = &map;
+      return *this;
+    }
+
+    /// Sets the flow map.
+
+    /// Sets the flow map.
+    /// \return \c (*this)
+    Circulation& flowMap(FlowMap& map) {
+      if (_local_flow) {
+	delete _flow;
+	_local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// Returns the flow map.
+
+    /// \return The flow map.
+    ///
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /// Sets the elevator.
+
+    /// Sets the elevator.
+    /// \return \c (*this)
+    Circulation& elevator(Elevator& elevator) {
+      if (_local_level) {
+	delete _level;
+	_local_level = false;
+      }
+      _level = &elevator;
+      return *this;
+    }
+
+    /// Returns the elevator.
+
+    /// \return The elevator.
+    ///
+    const Elevator& elevator() {
+      return *_level;
+    }
+    
+    /// Sets the tolerance used by algorithm.
+
+    /// Sets the tolerance used by algorithm.
+    ///
+    Circulation& tolerance(const Tolerance& tolerance) const {
+      _tol = tolerance;
+      return *this;
+    } 
+
+    /// Returns the tolerance used by algorithm.
+
+    /// Returns the tolerance used by algorithm.
+    ///
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+    
+    /// \name Execution control The simplest way to execute the
+    /// algorithm is to use one of the member functions called \c
+    /// run().  
+    /// \n 
+    /// If you need more control on initial solution or execution then
+    /// you have to call one \ref init() function and then the start()
+    /// function.
+    
+    ///@{
+  
+    /// Initializes the internal data structures.
+
+    /// Initializes the internal data structures. This function sets
+    /// all flow values to the lower bound.  
+    /// \return This function returns false if the initialization
+    /// process found a barrier.
+    void init() 
+    {
+      createStructures();
+
+      for(NodeIt n(_g);n!=INVALID;++n) {
+	_excess->set(n, (*_delta)[n]);
+      }
+     
+      for (EdgeIt 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]);
+      }
+
+      typename Graph::template NodeMap<bool> reached(_g, false);
+
+
+      // global relabeling tested, but in general case it provides
+      // worse performance for random graphs
+      _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.
+
+    /// Initializes the internal data structures. This functions uses
+    /// greedy approach to construct the initial solution. 
+    void greedyInit() 
+    {
+      createStructures();
+     
+      for(NodeIt n(_g);n!=INVALID;++n) {
+	_excess->set(n, (*_delta)[n]);
+      }
+     
+      for (EdgeIt 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)] + (*_up)[e]);
+	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[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);
+    }
+
+    ///Starts the algorithm
+
+    ///This function starts the algorithm.
+    ///\return This function returns true if it found a feasible circulation.
+    ///
+    ///\sa barrier()
+    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(OutEdgeIt 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] + exc);
+	      if(!_level->active(v) && _tol.positive((*_excess)[v]))
+		_level->activate(v);
+	      exc-=fc;
+	    }
+	  } 
+	  else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
+	}
+	for(InEdgeIt 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 circulation algorithm.  
+
+    /// Runs the circulation algorithm.
+    /// \note fc.run() is just a shortcut of the following code.
+    /// \code
+    ///   fc.greedyInit();
+    ///   return fc.start();
+    /// \endcode
+    bool run() {
+      greedyInit();
+      return start();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Circulation algorithm can be obtained using
+    /// these functions.
+    /// \n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+    
+    ///Returns a barrier
+    
+    ///Barrier is a set \e B of nodes for which
+    /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
+    ///holds. The existence of a set with this property prooves that a feasible
+    ///flow cannot exists.
+    ///\sa checkBarrier()
+    ///\sa run()
+    template<class GT>
+    void barrierMap(GT &bar) 
+    {
+      for(NodeIt n(_g);n!=INVALID;++n)
+	bar.set(n, (*_level)[n] >= _el);
+    }  
+
+    ///Returns true if the node is in the barrier
+
+    ///Returns true if the node is in the barrier
+    ///\sa barrierMap()
+    bool barrier(const Node& node) 
+    {
+      return (*_level)[node] >= _el;
+    }  
+
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Sets the \c flowMap to the flow on the edges. This method can
+    /// be called after the second phase of algorithm.
+    Value flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+
+    /// @}
+
+    /// \name Checker Functions
+    /// The feasibility  of the results can be checked using
+    /// these functions.
+    /// \n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    ///Check if the  \c flow is a feasible circulation
+    bool checkFlow() {
       for(EdgeIt e(_g);e!=INVALID;++e)
-	if(x[e]<_lo[e]||x[e]>_up[e]) return false;
+	if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
       for(NodeIt n(_g);n!=INVALID;++n)
 	{
-	  Value dif=-_delta[n];
-	  for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];
-	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];
+	  Value dif=-(*_delta)[n];
+	  for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
+	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
 	  if(_tol.negative(dif)) return false;
 	}
       return true;
-    };
-    ///Check if the default \c x is a feasible circulation
-    bool checkX() { return checkX(_x); }
+    }
 
-    ///Check if \c bar is a real barrier
-
-    ///Check if \c bar is a real barrier
-    ///\sa barrier()
-    template<class GT>
-    bool checkBarrier(GT &bar) 
-    {
-      Value delta=0;
-      for(NodeIt n(_g);n!=INVALID;++n)
-	if(bar[n])
-	  delta-=_delta[n];
-      for(EdgeIt e(_g);e!=INVALID;++e)
-	{
-	  Node s=_g.source(e);
-	  Node t=_g.target(e);
-	  if(bar[s]&&!bar[t]) delta+=_up[e];
-	  else if(bar[t]&&!bar[s]) delta-=_lo[e];
-	}
-      return _tol.negative(delta);
-    }  
     ///Check whether or not the last execution provides a barrier
 
     ///Check whether or not the last execution provides a barrier
     ///\sa barrier()
     bool checkBarrier() 
     {
-      typename Graph:: template NodeMap<bool> bar(_g);
-      barrier(bar);
-      return checkBarrier(bar);
+      Value delta=0;
+      for(NodeIt n(_g);n!=INVALID;++n)
+	if(barrier(n))
+	  delta-=(*_delta)[n];
+      for(EdgeIt 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);
     }
-    ///Run the algorithm
 
-    ///This function runs the algorithm.
-    ///\return This function returns -1 if it found a feasible circulation.
-    /// nonnegative values (including 0) mean that no feasible solution is
-    /// found. In this case the return value means an "empty level".
-    ///
-    ///\sa barrier()
-    int run() 
-    {
-      init();
-      
-#ifdef LEMON_CIRCULATION_DEBUG
-      for(NodeIt n(_g);n!=INVALID;++n)
-	std::cerr<< _levels[n] << ' ';
-      std::cerr << std::endl;
-#endif      
-      Node act;
-      Node bact=INVALID;
-      Node last_activated=INVALID;
-      while((act=_levels.highestActive())!=INVALID) {
-	int actlevel=_levels[act];
-	int tlevel;
-	int mlevel=_node_num;
-	Value exc=_excess[act];
-	Value fc;
-	
-#ifdef LEMON_CIRCULATION_DEBUG
-	for(NodeIt n(_g);n!=INVALID;++n)
-	  std::cerr<< _levels[n] << ' ';
-	std::cerr << std::endl;
-	std::cerr << "Process node " << _g.id(act)
-		  << " on level " << actlevel
-		  << " with excess " << exc
-		  << std::endl;
-#endif
-	for(OutEdgeIt e(_g,act);e!=INVALID; ++e)
-	  if((fc=_up[e]-_x[e])>0)
-	    if((tlevel=_levels[_g.target(e)])<actlevel)
-	      if(fc<=exc) {
-		_x[e]=_up[e];
-		addExcess(_g.target(e),fc);
-		exc-=fc;
-#ifdef LEMON_CIRCULATION_DEBUG
-		std::cerr << "  Push " << fc
-			  << " toward " << _g.id(_g.target(e)) << std::endl;
-#endif
-	      }
-	      else {
-		_x[e]+=exc;
-		addExcess(_g.target(e),exc);
-		//exc=0;
-		_excess[act]=0;
-		_levels.deactivate(act);
-#ifdef LEMON_CIRCULATION_DEBUG
-		std::cerr << "  Push " << exc
-			  << " toward " << _g.id(_g.target(e)) << std::endl;
-		std::cerr << "  Deactivate\n";
-#endif
-		goto next_l;
-	      }
-	    else if(tlevel<mlevel) mlevel=tlevel;
-	
-	for(InEdgeIt e(_g,act);e!=INVALID; ++e)
-	  if((fc=_x[e]-_lo[e])>0)
-	    if((tlevel=_levels[_g.source(e)])<actlevel)
-	      if(fc<=exc) {
-		_x[e]=_lo[e];
-		addExcess(_g.source(e),fc);
-		exc-=fc;
-#ifdef LEMON_CIRCULATION_DEBUG
-		std::cerr << "  Push " << fc
-			  << " toward " << _g.id(_g.source(e)) << std::endl;
-#endif
-	      }
-	      else {
-		_x[e]-=exc;
-		addExcess(_g.source(e),exc);
-		//exc=0;
-		_excess[act]=0;
-		_levels.deactivate(act);
-#ifdef LEMON_CIRCULATION_DEBUG
-		std::cerr << "  Push " << exc
-			  << " toward " << _g.id(_g.source(e)) << std::endl;
-		std::cerr << "  Deactivate\n";
-#endif
-		goto next_l;
-              }
-	    else if(tlevel<mlevel) mlevel=tlevel;
-   
-	_excess[act]=exc;
-	if(!_tol.positive(exc)) _levels.deactivate(act);
-	else if(mlevel==_node_num) {
-	  _levels.liftHighestActiveToTop();
-#ifdef LEMON_CIRCULATION_DEBUG
-	  std::cerr << "  Lift to level " << _node_num << std::endl;
-#endif
-	  return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
-	}
-	else {
-	  _levels.liftHighestActive(mlevel+1);
-#ifdef LEMON_CIRCULATION_DEBUG
-	  std::cerr << "  Lift to level " << mlevel+1 << std::endl;
-#endif
-	  if(_levels.onLevel(actlevel)==0)
-	    return actlevel;
-	}
-      next_l:
-	;
-      }
-#ifdef LEMON_CIRCULATION_DEBUG
-      std::cerr << "Feasible flow found.\n";
-#endif
-      return -1;
-    }
-    
-    ///Return a barrier
-    
-    ///Barrier is a set \e B of nodes for which
-    /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
-    ///holds. The existence of a set with this property prooves that a feasible
-    ///flow cannot exists.
-    ///\pre The run() must have been executed, and its return value was -1.
-    ///\sa checkBarrier()
-    ///\sa run()
-    template<class GT>
-    void barrier(GT &bar,int empty_level=-1) 
-    {
-      if(empty_level==-1)
-	for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
-      for(NodeIt n(_g);n!=INVALID;++n)
-	bar[n] = _levels[n]>empty_level;
-    }  
+    /// @}
 
   };