///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.

   ///This class implements a preflow algorithm

   ///This class implements a push-relabel algorithm

   template<class Graph,

   template<class _Graph,

 	   class Value,

 	   class _LCapMap=typename _Graph::template EdgeMap<int>,

 	   class FlowMap=typename Graph::template EdgeMap<Value>,

 	   class _UCapMap=_LCapMap,

 	   class LCapMap=typename Graph::template EdgeMap<Value>,

 	   class _DeltaMap=typename _Graph::template NodeMap<

 	   class UCapMap=LCapMap,

 	     typename _UCapMap::Value>,

 	   class DeltaMap=typename Graph::template NodeMap<Value>

 	   class _Traits=CirculationDefaultTraits<_Graph, _LCapMap, 

 	    >

 						  _UCapMap, _DeltaMap> >

     typedef typename Graph::Node Node;

     typedef typename Graph::NodeIt NodeIt;

     typedef _Traits Traits;

     typedef typename Graph::Edge Edge;

     typedef typename Traits::Graph Graph;

     typedef typename Graph::EdgeIt EdgeIt;

     GRAPH_TYPEDEFS(typename Graph);

     typedef typename Graph::InEdgeIt InEdgeIt;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     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 LCapMap &_lo;

     const LCapMap *_lo;

     const UCapMap &_up;

     const UCapMap *_up;

     const DeltaMap &_delta;

     const DeltaMap *_delta;

     FlowMap &_x;

     Tolerance<Value> _tol;

     FlowMap *_flow;

     Elevator<Graph,typename Graph::Node> _levels;

     bool _local_flow;

     typename Graph::template NodeMap<Value> _excess;

     Elevator* _level;

     bool _local_level;

     ExcessMap* _excess;

     Tolerance _tol;

     int _el;

     ///\e

     Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,

     typedef Circulation Create;

 		const DeltaMap &delta,FlowMap &x) :

       _g(g),

     ///\name Named template parameters

       _node_num(countNodes(g)),

       _lo(lo),_up(up),_delta(delta),_x(x),

     ///@{

       _levels(g,_node_num),

       _excess(g)

     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();

     }

     void addExcess(Node n,Value v)

     void createStructures() {

     {

       _node_num = _el = countNodes(_g);

       if(_tol.positive(_excess[n]+=v))

 	{

       if (!_flow) {

 	  if(!_levels.active(n)) _levels.activate(n);

 	_flow = Traits::createFlowMap(_g);

 	}

 	_local_flow = true;

       else if(_levels.active(n)) _levels.deactivate(n);

       }

     }

       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 \c (*this)

     Circulation& lowerCapMap(const LCapMap& map) {

       _lo = ↦

       return *this;

     }

     /// Sets the upper bound capacity map.

     /// Sets the upper bound capacity map.

     /// \return \c (*this)

     Circulation& upperCapMap(const LCapMap& map) {

       _up = ↦

       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 = ↦

       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 = ↦

       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.

       _x=_lo;

       for (EdgeIt e(_g);e!=INVALID;++e) {

 	_flow->set(e, (*_lo)[e]);

       for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=_delta[n];

 	_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() {

 	{

 	if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;

 	  _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);

     }

   public:

     ///Check if \c x is a feasible circulation

     template<class FT>

     bool checkX(FT &x) {

       for(EdgeIt e(_g);e!=INVALID;++e)

 	if(x[e]<_lo[e]||x[e]>_up[e]) return false;

 	  Value dif=-_delta[n];

 	  Value dif=-(*_delta)[n];

 	  for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];

 	  for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];

 	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];

 	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];

     };

     }

     ///Check if the default \c x is a feasible circulation

     bool checkX() { return checkX(_x); }

     ///Check whether or not the last execution provides a barrier

     ///Check if \c bar is a real barrier

     ///Check whether or not the last execution provides a barrier

     ///Check if \c bar is a real barrier

     template<class GT>

     bool checkBarrier() 

     bool checkBarrier(GT &bar) 

 	if(bar[n])

 	if(barrier(n))

 	  delta-=_delta[n];

 	  delta-=(*_delta)[n];

 	  if(bar[s]&&!bar[t]) delta+=_up[e];

 	  if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];

 	  else if(bar[t]&&!bar[s]) delta-=_lo[e];

 	  else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];

     }  

     }

     ///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);

     }

     ///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;

     }