/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2007

  * 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 max_flow

 ///\file

 ///\brief Push-prelabel algorithm for finding a feasible circulation.

 ///

 namespace lemon {

   ///Preflow algorithm for the Network Circulation Problem.

   ///\ingroup max_flow

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