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