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

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

    }

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