Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 2363)
+++ lemon/Makefile.am	(revision 2375)
@@ -34,4 +34,5 @@
 
 lemon_HEADERS += \
+	lemon/circulation.h \
 	lemon/bellman_ford.h \
 	lemon/bfs.h \
Index: lemon/circulation.h
===================================================================
--- lemon/circulation.h	(revision 2375)
+++ lemon/circulation.h	(revision 2375)
@@ -0,0 +1,299 @@
+/* -*- C++ -*-
+ * lemon/preflow_matching.h - Part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2005 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/graph_utils.h>
+#include <iostream>
+#include <queue>
+#include <lemon/tolerance.h>
+#include <lemon/elevator.h>
+
+///\ingroup flowalgs
+///\file
+///\brief Push-prelabel algorithm for finding a feasible circulation.
+///
+namespace lemon {
+
+  ///Preflow algorithm for the Network Circulation Problem.
+  
+  ///\ingroup flowalgs
+  ///This class implements a preflow 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>
+	    >
+  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;
+    
+
+    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;
+
+  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)
+    {
+    }
+    
+  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 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);
+    }
+
+  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;
+      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];
+	  if(_tol.negative(dif)) return false;
+	}
+      return true;
+    };
+    ///Check if the default \c x is a feasible circulation
+    bool checkX() { return checkX(_x); }
+
+    ///Chech if \c bar is a real barrier
+
+    ///Chech 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);
+    }  
+    ///Chech whether or not the last execution provides a barrier
+
+    ///Chech 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.liftHighestActiveTo(_node_num);
+#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.liftHighestActiveTo(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;
+    }  
+
+  };
+  
+}
+
+#endif
