alpar@2375: /* -*- C++ -*-
alpar@2375:  *
alpar@2391:  * This file is a part of LEMON, a generic C++ optimization library
alpar@2391:  *
alpar@2391:  * Copyright (C) 2003-2007
alpar@2391:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@2375:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@2375:  *
alpar@2375:  * Permission to use, modify and distribute this software is granted
alpar@2375:  * provided that this copyright notice appears in all copies. For
alpar@2375:  * precise terms see the accompanying LICENSE file.
alpar@2375:  *
alpar@2375:  * This software is provided "AS IS" with no warranty of any kind,
alpar@2375:  * express or implied, and with no claim as to its suitability for any
alpar@2375:  * purpose.
alpar@2375:  *
alpar@2375:  */
alpar@2375: 
alpar@2375: #ifndef LEMON_CIRCULATION_H
alpar@2375: #define LEMON_CIRCULATION_H
alpar@2375: 
alpar@2375: #include <lemon/graph_utils.h>
alpar@2375: #include <iostream>
alpar@2375: #include <queue>
alpar@2375: #include <lemon/tolerance.h>
alpar@2375: #include <lemon/elevator.h>
alpar@2375: 
deba@2376: ///\ingroup max_flow
alpar@2375: ///\file
alpar@2375: ///\brief Push-prelabel algorithm for finding a feasible circulation.
alpar@2375: ///
alpar@2375: namespace lemon {
alpar@2375: 
deba@2526:   /// \brief Default traits class of Circulation class.
deba@2526:   ///
deba@2526:   /// Default traits class of Circulation class.
deba@2526:   /// \param _Graph Graph type.
deba@2526:   /// \param _CapacityMap Type of capacity map.
deba@2526:   template <typename _Graph, typename _LCapMap, 
deba@2526: 	    typename _UCapMap, typename _DeltaMap>
deba@2526:   struct CirculationDefaultTraits {
deba@2526: 
deba@2526:     /// \brief The graph type the algorithm runs on. 
deba@2526:     typedef _Graph Graph;
deba@2526: 
deba@2526:     /// \brief The type of the map that stores the circulation lower
deba@2526:     /// bound.
deba@2526:     ///
deba@2526:     /// The type of the map that stores the circulation lower bound.
deba@2526:     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
deba@2526:     typedef _LCapMap LCapMap;
deba@2526: 
deba@2526:     /// \brief The type of the map that stores the circulation upper
deba@2526:     /// bound.
deba@2526:     ///
deba@2526:     /// The type of the map that stores the circulation upper bound.
deba@2526:     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
deba@2526:     typedef _UCapMap UCapMap;
deba@2526: 
deba@2526:     /// \brief The type of the map that stores the upper bound of
deba@2526:     /// node excess.
deba@2526:     ///
deba@2526:     /// The type of the map that stores the lower bound of node
deba@2526:     /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
deba@2526:     /// concept.
deba@2526:     typedef _DeltaMap DeltaMap;
deba@2526: 
deba@2526:     /// \brief The type of the length of the edges.
deba@2526:     typedef typename DeltaMap::Value Value;
deba@2526: 
deba@2526:     /// \brief The map type that stores the flow values.
deba@2526:     ///
deba@2526:     /// The map type that stores the flow values. 
deba@2526:     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
deba@2526:     typedef typename Graph::template EdgeMap<Value> FlowMap;
deba@2526: 
deba@2526:     /// \brief Instantiates a FlowMap.
deba@2526:     ///
deba@2526:     /// This function instantiates a \ref FlowMap. 
deba@2526:     /// \param graph The graph, to which we would like to define the flow map.
deba@2526:     static FlowMap* createFlowMap(const Graph& graph) {
deba@2526:       return new FlowMap(graph);
deba@2526:     }
deba@2526: 
deba@2526:     /// \brief The eleavator type used by Circulation algorithm.
deba@2526:     /// 
deba@2526:     /// The elevator type used by Circulation algorithm.
deba@2526:     ///
deba@2526:     /// \sa Elevator
deba@2526:     /// \sa LinkedElevator
deba@2526:     typedef Elevator<Graph, typename Graph::Node> Elevator;
deba@2526:     
deba@2526:     /// \brief Instantiates an Elevator.
deba@2526:     ///
deba@2526:     /// This function instantiates a \ref Elevator. 
deba@2526:     /// \param graph The graph, to which we would like to define the elevator.
deba@2526:     /// \param max_level The maximum level of the elevator.
deba@2526:     static Elevator* createElevator(const Graph& graph, int max_level) {
deba@2526:       return new Elevator(graph, max_level);
deba@2526:     }
deba@2526: 
deba@2526:     /// \brief The tolerance used by the algorithm
deba@2526:     ///
deba@2526:     /// The tolerance used by the algorithm to handle inexact computation.
deba@2526:     typedef Tolerance<Value> Tolerance;
deba@2526: 
deba@2526:   };
deba@2526:   
deba@2526:   ///Push-relabel algorithm for the Network Circulation Problem.
alpar@2375:   
deba@2378:   ///\ingroup max_flow
deba@2526:   ///This class implements a push-relabel algorithm
alpar@2375:   ///for the Network Circulation Problem.
alpar@2375:   ///The exact formulation of this problem is the following.
alpar@2450:   /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq -delta(v)\quad \forall v\in V \f]
alpar@2375:   /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
alpar@2375:   ///
deba@2526:   template<class _Graph,
deba@2526: 	   class _LCapMap=typename _Graph::template EdgeMap<int>,
deba@2526: 	   class _UCapMap=_LCapMap,
deba@2526: 	   class _DeltaMap=typename _Graph::template NodeMap<
deba@2526: 	     typename _UCapMap::Value>,
deba@2526: 	   class _Traits=CirculationDefaultTraits<_Graph, _LCapMap, 
deba@2526: 						  _UCapMap, _DeltaMap> >
alpar@2375:   class Circulation {
deba@2526: 
deba@2526:     typedef _Traits Traits;
deba@2526:     typedef typename Traits::Graph Graph;
deba@2526:     GRAPH_TYPEDEFS(typename Graph);
deba@2526: 
deba@2526:     typedef typename Traits::Value Value;
deba@2526: 
deba@2526:     typedef typename Traits::LCapMap LCapMap;
deba@2526:     typedef typename Traits::UCapMap UCapMap;
deba@2526:     typedef typename Traits::DeltaMap DeltaMap;
deba@2526:     typedef typename Traits::FlowMap FlowMap;
deba@2526:     typedef typename Traits::Elevator Elevator;
deba@2526:     typedef typename Traits::Tolerance Tolerance;
deba@2526: 
deba@2526:     typedef typename Graph::template NodeMap<Value> ExcessMap;
alpar@2375: 
alpar@2375:     const Graph &_g;
alpar@2375:     int _node_num;
alpar@2375: 
deba@2526:     const LCapMap *_lo;
deba@2526:     const UCapMap *_up;
deba@2526:     const DeltaMap *_delta;
deba@2526: 
deba@2526:     FlowMap *_flow;
deba@2526:     bool _local_flow;
deba@2526: 
deba@2526:     Elevator* _level;
deba@2526:     bool _local_level;
deba@2526: 
deba@2526:     ExcessMap* _excess;
deba@2526: 
deba@2526:     Tolerance _tol;
deba@2526:     int _el;
alpar@2375: 
alpar@2375:   public:
deba@2526: 
deba@2526:     typedef Circulation Create;
deba@2526: 
deba@2526:     ///\name Named template parameters
deba@2526: 
deba@2526:     ///@{
deba@2526: 
deba@2526:     template <typename _FlowMap>
deba@2526:     struct DefFlowMapTraits : public Traits {
deba@2526:       typedef _FlowMap FlowMap;
deba@2526:       static FlowMap *createFlowMap(const Graph&) {
deba@2526: 	throw UninitializedParameter();
deba@2526:       }
deba@2526:     };
deba@2526: 
deba@2526:     /// \brief \ref named-templ-param "Named parameter" for setting
deba@2526:     /// FlowMap type
deba@2526:     ///
deba@2526:     /// \ref named-templ-param "Named parameter" for setting FlowMap
deba@2526:     /// type
deba@2526:     template <typename _FlowMap>
deba@2526:     struct DefFlowMap 
deba@2526:       : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
deba@2526: 			   DefFlowMapTraits<_FlowMap> > {
deba@2526:       typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
deba@2526: 			  DefFlowMapTraits<_FlowMap> > Create;
deba@2526:     };
deba@2526: 
deba@2526:     template <typename _Elevator>
deba@2526:     struct DefElevatorTraits : public Traits {
deba@2526:       typedef _Elevator Elevator;
deba@2526:       static Elevator *createElevator(const Graph&, int) {
deba@2526: 	throw UninitializedParameter();
deba@2526:       }
deba@2526:     };
deba@2526: 
deba@2526:     /// \brief \ref named-templ-param "Named parameter" for setting
deba@2526:     /// Elevator type
deba@2526:     ///
deba@2526:     /// \ref named-templ-param "Named parameter" for setting Elevator
deba@2526:     /// type
deba@2526:     template <typename _Elevator>
deba@2526:     struct DefElevator 
deba@2526:       : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
deba@2526: 			   DefElevatorTraits<_Elevator> > {
deba@2526:       typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
deba@2526: 			  DefElevatorTraits<_Elevator> > Create;
deba@2526:     };
deba@2526: 
deba@2526:     template <typename _Elevator>
deba@2526:     struct DefStandardElevatorTraits : public Traits {
deba@2526:       typedef _Elevator Elevator;
deba@2526:       static Elevator *createElevator(const Graph& graph, int max_level) {
deba@2526: 	return new Elevator(graph, max_level);
deba@2526:       }
deba@2526:     };
deba@2526: 
deba@2526:     /// \brief \ref named-templ-param "Named parameter" for setting
deba@2526:     /// Elevator type
deba@2526:     ///
deba@2526:     /// \ref named-templ-param "Named parameter" for setting Elevator
deba@2526:     /// type. The Elevator should be standard constructor interface, ie.
deba@2526:     /// the graph and the maximum level should be passed to it.
deba@2526:     template <typename _Elevator>
deba@2526:     struct DefStandardElevator 
deba@2526:       : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 
deba@2526: 		       DefStandardElevatorTraits<_Elevator> > {
deba@2526:       typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
deba@2526: 		      DefStandardElevatorTraits<_Elevator> > Create;
deba@2526:     };    
deba@2526: 
deba@2526:     /// @}
deba@2526: 
deba@2527:   protected:
deba@2527: 
deba@2527:     Circulation() {}
deba@2527:   
deba@2527:   public:
deba@2527: 
deba@2526:     /// The constructor of the class.
deba@2526: 
deba@2526:     /// The constructor of the class. 
deba@2526:     /// \param g The directed graph the algorithm runs on. 
deba@2526:     /// \param lo The lower bound capacity of the edges. 
deba@2526:     /// \param up The upper bound capacity of the edges.
deba@2526:     /// \param delta The lower bound on node excess.
deba@2526:     Circulation(const Graph &g,const LCapMap &lo,
deba@2526: 		const UCapMap &up,const DeltaMap &delta) 
deba@2526:       : _g(g), _node_num(),
deba@2526: 	_lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
deba@2526: 	_level(0), _local_level(false), _excess(0), _el() {}
deba@2526: 
deba@2526:     /// Destrcutor.
deba@2526:     ~Circulation() {
deba@2526:       destroyStructures();
alpar@2375:     }
deba@2526: 
alpar@2375:   private:
alpar@2375: 
deba@2526:     void createStructures() {
deba@2526:       _node_num = _el = countNodes(_g);
deba@2526: 
deba@2526:       if (!_flow) {
deba@2526: 	_flow = Traits::createFlowMap(_g);
deba@2526: 	_local_flow = true;
deba@2526:       }
deba@2526:       if (!_level) {
deba@2526: 	_level = Traits::createElevator(_g, _node_num);
deba@2526: 	_local_level = true;
deba@2526:       }
deba@2526:       if (!_excess) {
deba@2526: 	_excess = new ExcessMap(_g);
deba@2526:       }
alpar@2375:     }
alpar@2375: 
deba@2526:     void destroyStructures() {
deba@2526:       if (_local_flow) {
deba@2526: 	delete _flow;
deba@2526:       }
deba@2526:       if (_local_level) {
deba@2526: 	delete _level;
deba@2526:       }
deba@2526:       if (_excess) {
deba@2526: 	delete _excess;
deba@2526:       }
alpar@2375:     }
alpar@2375: 
alpar@2375:   public:
deba@2526: 
deba@2526:     /// Sets the lower bound capacity map.
deba@2526: 
deba@2526:     /// Sets the lower bound capacity map.
deba@2526:     /// \return \c (*this)
deba@2526:     Circulation& lowerCapMap(const LCapMap& map) {
deba@2526:       _lo = &map;
deba@2526:       return *this;
deba@2526:     }
deba@2526: 
deba@2526:     /// Sets the upper bound capacity map.
deba@2526: 
deba@2526:     /// Sets the upper bound capacity map.
deba@2526:     /// \return \c (*this)
deba@2526:     Circulation& upperCapMap(const LCapMap& map) {
deba@2526:       _up = &map;
deba@2526:       return *this;
deba@2526:     }
deba@2526: 
deba@2526:     /// Sets the lower bound map on excess.
deba@2526: 
deba@2526:     /// Sets the lower bound map on excess.
deba@2526:     /// \return \c (*this)
deba@2526:     Circulation& deltaMap(const DeltaMap& map) {
deba@2526:       _delta = &map;
deba@2526:       return *this;
deba@2526:     }
deba@2526: 
deba@2526:     /// Sets the flow map.
deba@2526: 
deba@2526:     /// Sets the flow map.
deba@2526:     /// \return \c (*this)
deba@2526:     Circulation& flowMap(FlowMap& map) {
deba@2526:       if (_local_flow) {
deba@2526: 	delete _flow;
deba@2526: 	_local_flow = false;
deba@2526:       }
deba@2526:       _flow = &map;
deba@2526:       return *this;
deba@2526:     }
deba@2526: 
deba@2526:     /// Returns the flow map.
deba@2526: 
deba@2526:     /// \return The flow map.
deba@2526:     ///
deba@2526:     const FlowMap& flowMap() {
deba@2526:       return *_flow;
deba@2526:     }
deba@2526: 
deba@2526:     /// Sets the elevator.
deba@2526: 
deba@2526:     /// Sets the elevator.
deba@2526:     /// \return \c (*this)
deba@2526:     Circulation& elevator(Elevator& elevator) {
deba@2526:       if (_local_level) {
deba@2526: 	delete _level;
deba@2526: 	_local_level = false;
deba@2526:       }
deba@2526:       _level = &elevator;
deba@2526:       return *this;
deba@2526:     }
deba@2526: 
deba@2526:     /// Returns the elevator.
deba@2526: 
deba@2526:     /// \return The elevator.
deba@2526:     ///
deba@2526:     const Elevator& elevator() {
deba@2526:       return *_level;
deba@2526:     }
deba@2526:     
deba@2526:     /// Sets the tolerance used by algorithm.
deba@2526: 
deba@2526:     /// Sets the tolerance used by algorithm.
deba@2526:     ///
deba@2526:     Circulation& tolerance(const Tolerance& tolerance) const {
deba@2526:       _tol = tolerance;
deba@2526:       return *this;
deba@2526:     } 
deba@2526: 
deba@2526:     /// Returns the tolerance used by algorithm.
deba@2526: 
deba@2526:     /// Returns the tolerance used by algorithm.
deba@2526:     ///
deba@2526:     const Tolerance& tolerance() const {
deba@2526:       return tolerance;
deba@2526:     } 
deba@2526:     
kpeter@2533:     /// \name Execution control
kpeter@2533:     /// The simplest way to execute the algorithm is to use one of the
kpeter@2533:     /// member functions called \c run().
deba@2526:     /// \n 
deba@2526:     /// If you need more control on initial solution or execution then
deba@2526:     /// you have to call one \ref init() function and then the start()
deba@2526:     /// function.
deba@2526:     
deba@2526:     ///@{
deba@2526:   
deba@2526:     /// Initializes the internal data structures.
deba@2526: 
deba@2526:     /// Initializes the internal data structures. This function sets
deba@2526:     /// all flow values to the lower bound.  
deba@2526:     /// \return This function returns false if the initialization
deba@2526:     /// process found a barrier.
deba@2526:     void init() 
deba@2526:     {
deba@2526:       createStructures();
deba@2526: 
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n) {
deba@2526: 	_excess->set(n, (*_delta)[n]);
deba@2526:       }
deba@2526:      
deba@2526:       for (EdgeIt e(_g);e!=INVALID;++e) {
deba@2526: 	_flow->set(e, (*_lo)[e]);
deba@2526: 	_excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
deba@2526: 	_excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
deba@2526:       }
deba@2526: 
deba@2526:       typename Graph::template NodeMap<bool> reached(_g, false);
deba@2526: 
deba@2526: 
deba@2526:       // global relabeling tested, but in general case it provides
deba@2526:       // worse performance for random graphs
deba@2526:       _level->initStart();
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n)
deba@2526: 	_level->initAddItem(n);
deba@2526:       _level->initFinish();
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n)
deba@2526: 	if(_tol.positive((*_excess)[n]))
deba@2526: 	  _level->activate(n);
deba@2526:     }
deba@2526: 
deba@2526:     /// Initializes the internal data structures.
deba@2526: 
deba@2526:     /// Initializes the internal data structures. This functions uses
deba@2526:     /// greedy approach to construct the initial solution. 
deba@2526:     void greedyInit() 
deba@2526:     {
deba@2526:       createStructures();
deba@2526:      
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n) {
deba@2526: 	_excess->set(n, (*_delta)[n]);
deba@2526:       }
deba@2526:      
deba@2526:       for (EdgeIt e(_g);e!=INVALID;++e) {
deba@2526: 	if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
deba@2526: 	  _flow->set(e, (*_up)[e]);
deba@2526: 	  _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
deba@2526: 	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
deba@2526: 	} else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
deba@2526: 	  _flow->set(e, (*_lo)[e]);
deba@2526: 	  _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
deba@2526: 	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
deba@2526: 	} else {
deba@2526: 	  Value fc = -(*_excess)[_g.target(e)];
deba@2526: 	  _flow->set(e, fc);
deba@2526: 	  _excess->set(_g.target(e), 0);
deba@2526: 	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
deba@2526: 	}
deba@2526:       }
deba@2526:      
deba@2526:       _level->initStart();
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n)
deba@2526: 	_level->initAddItem(n);
deba@2526:       _level->initFinish();
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n)
deba@2526: 	if(_tol.positive((*_excess)[n]))
deba@2526: 	  _level->activate(n);
deba@2526:     }
deba@2526: 
deba@2526:     ///Starts the algorithm
deba@2526: 
deba@2526:     ///This function starts the algorithm.
deba@2526:     ///\return This function returns true if it found a feasible circulation.
deba@2526:     ///
deba@2526:     ///\sa barrier()
deba@2526:     bool start() 
deba@2526:     {
deba@2526:       
deba@2526:       Node act;
deba@2526:       Node bact=INVALID;
deba@2526:       Node last_activated=INVALID;
deba@2526:       while((act=_level->highestActive())!=INVALID) {
deba@2526: 	int actlevel=(*_level)[act];
deba@2526: 	int mlevel=_node_num;
deba@2526: 	Value exc=(*_excess)[act];
deba@2526: 	
deba@2526: 	for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
deba@2526: 	  Node v = _g.target(e);
deba@2526: 	  Value fc=(*_up)[e]-(*_flow)[e];
deba@2526: 	  if(!_tol.positive(fc)) continue;
deba@2526: 	  if((*_level)[v]<actlevel) {
deba@2526: 	    if(!_tol.less(fc, exc)) {
deba@2526: 	      _flow->set(e, (*_flow)[e] + exc);
deba@2526: 	      _excess->set(v, (*_excess)[v] + exc);
deba@2526: 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))
deba@2526: 		_level->activate(v);
deba@2526: 	      _excess->set(act,0);
deba@2526: 	      _level->deactivate(act);
deba@2526: 	      goto next_l;
deba@2526: 	    }
deba@2526: 	    else {
deba@2526: 	      _flow->set(e, (*_up)[e]);
deba@2526: 	      _excess->set(v, (*_excess)[v] + exc);
deba@2526: 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))
deba@2526: 		_level->activate(v);
deba@2526: 	      exc-=fc;
deba@2526: 	    }
deba@2526: 	  } 
deba@2526: 	  else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
deba@2526: 	}
deba@2526: 	for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
deba@2526: 	  Node v = _g.source(e);
deba@2526: 	  Value fc=(*_flow)[e]-(*_lo)[e];
deba@2526: 	  if(!_tol.positive(fc)) continue;
deba@2526: 	  if((*_level)[v]<actlevel) {
deba@2526: 	    if(!_tol.less(fc, exc)) {
deba@2526: 	      _flow->set(e, (*_flow)[e] - exc);
deba@2526: 	      _excess->set(v, (*_excess)[v] + exc);
deba@2526: 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))
deba@2526: 		_level->activate(v);
deba@2526: 	      _excess->set(act,0);
deba@2526: 	      _level->deactivate(act);
deba@2526: 	      goto next_l;
deba@2526: 	    }
deba@2526: 	    else {
deba@2526: 	      _flow->set(e, (*_lo)[e]);
deba@2526: 	      _excess->set(v, (*_excess)[v] + fc);
deba@2526: 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))
deba@2526: 		_level->activate(v);
deba@2526: 	      exc-=fc;
deba@2526: 	    }
deba@2526: 	  } 
deba@2526: 	  else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
deba@2526: 	}
deba@2526:    
deba@2526: 	_excess->set(act, exc);
deba@2526: 	if(!_tol.positive(exc)) _level->deactivate(act);
deba@2526: 	else if(mlevel==_node_num) {
deba@2526: 	  _level->liftHighestActiveToTop();
deba@2526: 	  _el = _node_num;
deba@2526: 	  return false;
deba@2526: 	}
deba@2526: 	else {
deba@2526: 	  _level->liftHighestActive(mlevel+1);
deba@2526: 	  if(_level->onLevel(actlevel)==0) {
deba@2526: 	    _el = actlevel;
deba@2526: 	    return false;
deba@2526: 	  }
deba@2526: 	}
deba@2526:       next_l:
deba@2526: 	;
deba@2526:       }
deba@2526:       return true;
deba@2526:     }
deba@2526: 
deba@2526:     /// Runs the circulation algorithm.  
deba@2526: 
deba@2526:     /// Runs the circulation algorithm.
deba@2526:     /// \note fc.run() is just a shortcut of the following code.
deba@2526:     /// \code
deba@2526:     ///   fc.greedyInit();
deba@2526:     ///   return fc.start();
deba@2526:     /// \endcode
deba@2526:     bool run() {
deba@2526:       greedyInit();
deba@2526:       return start();
deba@2526:     }
deba@2526: 
deba@2526:     /// @}
deba@2526: 
deba@2526:     /// \name Query Functions
deba@2526:     /// The result of the %Circulation algorithm can be obtained using
deba@2526:     /// these functions.
deba@2526:     /// \n
deba@2526:     /// Before the use of these functions,
deba@2526:     /// either run() or start() must be called.
deba@2526:     
deba@2526:     ///@{
deba@2526:     
deba@2526:     ///Returns a barrier
deba@2526:     
deba@2526:     ///Barrier is a set \e B of nodes for which
deba@2526:     /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
deba@2526:     ///holds. The existence of a set with this property prooves that a feasible
deba@2526:     ///flow cannot exists.
deba@2526:     ///\sa checkBarrier()
deba@2526:     ///\sa run()
deba@2526:     template<class GT>
deba@2526:     void barrierMap(GT &bar) 
deba@2526:     {
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n)
deba@2526: 	bar.set(n, (*_level)[n] >= _el);
deba@2526:     }  
deba@2526: 
deba@2526:     ///Returns true if the node is in the barrier
deba@2526: 
deba@2526:     ///Returns true if the node is in the barrier
deba@2526:     ///\sa barrierMap()
deba@2526:     bool barrier(const Node& node) 
deba@2526:     {
deba@2526:       return (*_level)[node] >= _el;
deba@2526:     }  
deba@2526: 
deba@2526:     /// \brief Returns the flow on the edge.
deba@2526:     ///
deba@2526:     /// Sets the \c flowMap to the flow on the edges. This method can
deba@2526:     /// be called after the second phase of algorithm.
deba@2526:     Value flow(const Edge& edge) const {
deba@2526:       return (*_flow)[edge];
deba@2526:     }
deba@2526: 
deba@2526:     /// @}
deba@2526: 
deba@2526:     /// \name Checker Functions
deba@2526:     /// The feasibility  of the results can be checked using
deba@2526:     /// these functions.
deba@2526:     /// \n
deba@2526:     /// Before the use of these functions,
deba@2526:     /// either run() or start() must be called.
deba@2526:     
deba@2526:     ///@{
deba@2526: 
deba@2526:     ///Check if the  \c flow is a feasible circulation
deba@2526:     bool checkFlow() {
alpar@2375:       for(EdgeIt e(_g);e!=INVALID;++e)
deba@2526: 	if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	{
deba@2526: 	  Value dif=-(*_delta)[n];
deba@2526: 	  for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
deba@2526: 	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
alpar@2375: 	  if(_tol.negative(dif)) return false;
alpar@2375: 	}
alpar@2375:       return true;
deba@2526:     }
alpar@2375: 
alpar@2408:     ///Check whether or not the last execution provides a barrier
alpar@2375: 
alpar@2408:     ///Check whether or not the last execution provides a barrier
alpar@2375:     ///\sa barrier()
alpar@2375:     bool checkBarrier() 
alpar@2375:     {
deba@2526:       Value delta=0;
deba@2526:       for(NodeIt n(_g);n!=INVALID;++n)
deba@2526: 	if(barrier(n))
deba@2526: 	  delta-=(*_delta)[n];
deba@2526:       for(EdgeIt e(_g);e!=INVALID;++e)
deba@2526: 	{
deba@2526: 	  Node s=_g.source(e);
deba@2526: 	  Node t=_g.target(e);
deba@2526: 	  if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
deba@2526: 	  else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
deba@2526: 	}
deba@2526:       return _tol.negative(delta);
alpar@2375:     }
alpar@2375: 
deba@2526:     /// @}
alpar@2375: 
alpar@2375:   };
alpar@2375:   
alpar@2375: }
alpar@2375: 
alpar@2375: #endif