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

   /// \brief Default traits class of Circulation class.

   ///

   /// Default traits class of Circulation class.

   /// \param _Graph Graph type.

   /// \param _CapacityMap Type of capacity map.

   template <typename _Graph, typename _LCapMap, 

 	    typename _UCapMap, typename _DeltaMap>

   struct CirculationDefaultTraits {

     /// \brief The graph type the algorithm runs on. 

     typedef _Graph Graph;

     /// \brief The type of the map that stores the circulation lower

     /// bound.

     ///

     /// The type of the map that stores the circulation lower bound.

     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.

     typedef _LCapMap LCapMap;

     /// \brief The type of the map that stores the circulation upper

     /// bound.

     ///

     /// The type of the map that stores the circulation upper bound.

     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.

     typedef _UCapMap UCapMap;

     /// \brief The type of the map that stores the upper bound of

     /// node excess.

     ///

     /// The type of the map that stores the lower bound of node

     /// excess. It must meet the \ref concepts::ReadMap "ReadMap"

     /// concept.

     typedef _DeltaMap DeltaMap;

     /// \brief The type of the length of the edges.

     typedef typename DeltaMap::Value Value;

     /// \brief The map type that stores the flow values.

     ///

     /// The map type that stores the flow values. 

     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.

     typedef typename Graph::template EdgeMap<Value> FlowMap;

     /// \brief Instantiates a FlowMap.

     ///

     /// This function instantiates a \ref FlowMap. 

     /// \param graph The graph, to which we would like to define the flow map.

     static FlowMap* createFlowMap(const Graph& graph) {

       return new FlowMap(graph);

     }

     /// \brief The eleavator type used by Circulation algorithm.

     /// 

     /// The elevator type used by Circulation algorithm.

     ///

     /// \sa Elevator

     /// \sa LinkedElevator

     typedef Elevator<Graph, typename Graph::Node> Elevator;

     /// \brief Instantiates an Elevator.

     ///

     /// This function instantiates a \ref Elevator. 

     /// \param graph The graph, to which we would like to define the elevator.

     /// \param max_level The maximum level of the elevator.

     static Elevator* createElevator(const Graph& graph, int max_level) {

       return new Elevator(graph, max_level);

     }

     /// \brief The tolerance used by the algorithm

     ///

     /// The tolerance used by the algorithm to handle inexact computation.

     typedef Tolerance<Value> Tolerance;

   };

   ///Push-relabel algorithm for the Network Circulation Problem.

   ///\ingroup max_flow

   ///This class implements a push-relabel 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 _LCapMap=typename _Graph::template EdgeMap<int>,

 	   class _UCapMap=_LCapMap,

 	   class _DeltaMap=typename _Graph::template NodeMap<

 	     typename _UCapMap::Value>,

 	   class _Traits=CirculationDefaultTraits<_Graph, _LCapMap, 

 						  _UCapMap, _DeltaMap> >

   class Circulation {

     typedef _Traits Traits;

     typedef typename Traits::Graph Graph;

     GRAPH_TYPEDEFS(typename Graph);

     typedef typename Traits::Value Value;

     typedef typename Traits::LCapMap LCapMap;

     typedef typename Traits::UCapMap UCapMap;

     typedef typename Traits::DeltaMap DeltaMap;

     typedef typename Traits::FlowMap FlowMap;

     typedef typename Traits::Elevator Elevator;

     typedef typename Traits::Tolerance Tolerance;

     typedef typename Graph::template NodeMap<Value> ExcessMap;

     const Graph &_g;

     int _node_num;

     const LCapMap *_lo;

     const UCapMap *_up;

     const DeltaMap *_delta;

     FlowMap *_flow;

     bool _local_flow;

     Elevator* _level;

     bool _local_level;

     ExcessMap* _excess;

     Tolerance _tol;

     int _el;

   public:

     typedef Circulation Create;

     ///\name Named template parameters

     ///@{

     template <typename _FlowMap>

     struct DefFlowMapTraits : public Traits {

       typedef _FlowMap FlowMap;

       static FlowMap *createFlowMap(const Graph&) {

 	throw UninitializedParameter();

       }

     };

     /// \brief \ref named-templ-param "Named parameter" for setting

     /// FlowMap type

     ///

     /// \ref named-templ-param "Named parameter" for setting FlowMap

     /// type

     template <typename _FlowMap>

     struct DefFlowMap 

       : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 

 			   DefFlowMapTraits<_FlowMap> > {

       typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap, 

 			  DefFlowMapTraits<_FlowMap> > Create;

     };

     template <typename _Elevator>

     struct DefElevatorTraits : public Traits {

       typedef _Elevator Elevator;

       static Elevator *createElevator(const Graph&, int) {

 	throw UninitializedParameter();

       }

     };

     /// \brief \ref named-templ-param "Named parameter" for setting

     /// Elevator type

     ///

     /// \ref named-templ-param "Named parameter" for setting Elevator

     /// type

     template <typename _Elevator>

     struct DefElevator 

       : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 

 			   DefElevatorTraits<_Elevator> > {

       typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,

 			  DefElevatorTraits<_Elevator> > Create;

     };

     template <typename _Elevator>

     struct DefStandardElevatorTraits : public Traits {

       typedef _Elevator Elevator;

       static Elevator *createElevator(const Graph& graph, int max_level) {

 	return new Elevator(graph, max_level);

       }

     };

     /// \brief \ref named-templ-param "Named parameter" for setting

     /// Elevator type

     ///

     /// \ref named-templ-param "Named parameter" for setting Elevator

     /// type. The Elevator should be standard constructor interface, ie.

     /// the graph and the maximum level should be passed to it.

     template <typename _Elevator>

     struct DefStandardElevator 

       : public Circulation<Graph, LCapMap, UCapMap, DeltaMap, 

 		       DefStandardElevatorTraits<_Elevator> > {

       typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,

 		      DefStandardElevatorTraits<_Elevator> > Create;

     };    

     /// @}

   protected:

     Circulation() {}

   public:

     /// The constructor of the class.

     /// The constructor of the class. 

     /// \param g The directed graph the algorithm runs on. 

     /// \param lo The lower bound capacity of the edges. 

     /// \param up The upper bound capacity of the edges.

     /// \param delta The lower bound on node excess.

     Circulation(const Graph &g,const LCapMap &lo,

 		const UCapMap &up,const DeltaMap &delta) 

       : _g(g), _node_num(),

 	_lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),

 	_level(0), _local_level(false), _excess(0), _el() {}

     /// Destrcutor.

     ~Circulation() {

       destroyStructures();

     }

   private:

     void createStructures() {

       _node_num = _el = countNodes(_g);

       if (!_flow) {

 	_flow = Traits::createFlowMap(_g);

 	_local_flow = true;

       }

       if (!_level) {

 	_level = Traits::createElevator(_g, _node_num);

 	_local_level = true;

       }

       if (!_excess) {

 	_excess = new ExcessMap(_g);

       }

     }

     void destroyStructures() {

       if (_local_flow) {

 	delete _flow;

       }

       if (_local_level) {

 	delete _level;

       }

       if (_excess) {

 	delete _excess;

       }

     }

   public:

     /// Sets the lower bound capacity map.

     /// Sets the lower bound capacity map.

     /// \return \c (*this)

     Circulation& lowerCapMap(const LCapMap& map) {

       _lo = ↦

       return *this;

     }

     /// Sets the upper bound capacity map.

     /// Sets the upper bound capacity map.

     /// \return \c (*this)

     Circulation& upperCapMap(const LCapMap& map) {

       _up = ↦

       return *this;

     }

     /// Sets the lower bound map on excess.

     /// Sets the lower bound map on excess.

     /// \return \c (*this)

     Circulation& deltaMap(const DeltaMap& map) {

       _delta = ↦

       return *this;

     }

     /// Sets the flow map.

     /// Sets the flow map.

     /// \return \c (*this)

     Circulation& flowMap(FlowMap& map) {

       if (_local_flow) {

 	delete _flow;

 	_local_flow = false;

       }

       _flow = ↦

       return *this;

     }

     /// Returns the flow map.

     /// \return The flow map.

     ///

     const FlowMap& flowMap() {

       return *_flow;

     }

     /// Sets the elevator.

     /// Sets the elevator.

     /// \return \c (*this)

     Circulation& elevator(Elevator& elevator) {

       if (_local_level) {

 	delete _level;

 	_local_level = false;

       }

       _level = &elevator;

       return *this;

     }

     /// Returns the elevator.

     /// \return The elevator.

     ///

     const Elevator& elevator() {

       return *_level;

     }

     /// Sets the tolerance used by algorithm.

     /// Sets the tolerance used by algorithm.

     ///

     Circulation& tolerance(const Tolerance& tolerance) const {

       _tol = tolerance;

       return *this;

     } 

     /// Returns the tolerance used by algorithm.

     /// Returns the tolerance used by algorithm.

     ///

     const Tolerance& tolerance() const {

       return tolerance;

     } 

     /// \name Execution control

     /// The simplest way to execute the algorithm is to use one of the

     /// member functions called \c run().

     /// \n 

     /// If you need more control on initial solution or execution then

     /// you have to call one \ref init() function and then the start()

     /// function.

     ///@{

     /// Initializes the internal data structures.

     /// Initializes the internal data structures. This function sets

     /// all flow values to the lower bound.  

     /// \return This function returns false if the initialization

     /// process found a barrier.

     void init() 

     {

       createStructures();

       for(NodeIt n(_g);n!=INVALID;++n) {

 	_excess->set(n, (*_delta)[n]);

       }

       for (EdgeIt e(_g);e!=INVALID;++e) {

 	_flow->set(e, (*_lo)[e]);

 	_excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);

 	_excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);

       }

       // global relabeling tested, but in general case it provides

       // worse performance for random graphs

       _level->initStart();

       for(NodeIt n(_g);n!=INVALID;++n)

 	_level->initAddItem(n);

       _level->initFinish();

       for(NodeIt n(_g);n!=INVALID;++n)

 	if(_tol.positive((*_excess)[n]))

 	  _level->activate(n);

     }

     /// Initializes the internal data structures.

     /// Initializes the internal data structures. This functions uses

     /// greedy approach to construct the initial solution. 

     void greedyInit() 

     {

       createStructures();

       for(NodeIt n(_g);n!=INVALID;++n) {

 	_excess->set(n, (*_delta)[n]);

       }

       for (EdgeIt e(_g);e!=INVALID;++e) {

 	if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {

 	  _flow->set(e, (*_up)[e]);

 	  _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);

 	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);

 	} else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {

 	  _flow->set(e, (*_lo)[e]);

 	  _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);

 	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);

 	} else {

 	  Value fc = -(*_excess)[_g.target(e)];

 	  _flow->set(e, fc);

 	  _excess->set(_g.target(e), 0);

 	  _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);

 	}

       }

       _level->initStart();

       for(NodeIt n(_g);n!=INVALID;++n)

 	_level->initAddItem(n);

       _level->initFinish();

       for(NodeIt n(_g);n!=INVALID;++n)

 	if(_tol.positive((*_excess)[n]))

 	  _level->activate(n);

     }

     ///Starts the algorithm

     ///This function starts the algorithm.

     ///\return This function returns true if it found a feasible circulation.

     ///

     ///\sa barrier()

     bool start() 

     {

       Node act;

       Node bact=INVALID;

       Node last_activated=INVALID;

       while((act=_level->highestActive())!=INVALID) {

 	int actlevel=(*_level)[act];

 	int mlevel=_node_num;

 	Value exc=(*_excess)[act];

 	for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {

 	  Node v = _g.target(e);

 	  Value fc=(*_up)[e]-(*_flow)[e];

 	  if(!_tol.positive(fc)) continue;

 	  if((*_level)[v]<actlevel) {

 	    if(!_tol.less(fc, exc)) {

 	      _flow->set(e, (*_flow)[e] + exc);

 	      _excess->set(v, (*_excess)[v] + exc);

 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))

 		_level->activate(v);

 	      _excess->set(act,0);

 	      _level->deactivate(act);

 	      goto next_l;

 	    }

 	    else {

 	      _flow->set(e, (*_up)[e]);

 	      _excess->set(v, (*_excess)[v] + fc);

 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))

 		_level->activate(v);

 	      exc-=fc;

 	    }

 	  } 

 	  else if((*_level)[v]<mlevel) mlevel=(*_level)[v];

 	}

 	for(InEdgeIt e(_g,act);e!=INVALID; ++e) {

 	  Node v = _g.source(e);

 	  Value fc=(*_flow)[e]-(*_lo)[e];

 	  if(!_tol.positive(fc)) continue;

 	  if((*_level)[v]<actlevel) {

 	    if(!_tol.less(fc, exc)) {

 	      _flow->set(e, (*_flow)[e] - exc);

 	      _excess->set(v, (*_excess)[v] + exc);

 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))

 		_level->activate(v);

 	      _excess->set(act,0);

 	      _level->deactivate(act);

 	      goto next_l;

 	    }

 	    else {

 	      _flow->set(e, (*_lo)[e]);

 	      _excess->set(v, (*_excess)[v] + fc);

 	      if(!_level->active(v) && _tol.positive((*_excess)[v]))

 		_level->activate(v);

 	      exc-=fc;

 	    }

 	  } 

 	  else if((*_level)[v]<mlevel) mlevel=(*_level)[v];

 	}

 	_excess->set(act, exc);

 	if(!_tol.positive(exc)) _level->deactivate(act);

 	else if(mlevel==_node_num) {

 	  _level->liftHighestActiveToTop();

 	  _el = _node_num;

 	  return false;

 	}

 	else {

 	  _level->liftHighestActive(mlevel+1);

 	  if(_level->onLevel(actlevel)==0) {

 	    _el = actlevel;

 	    return false;

 	  }

 	}

       next_l:

 	;

       }

       return true;

     }

     /// Runs the circulation algorithm.  

     /// Runs the circulation algorithm.

     /// \note fc.run() is just a shortcut of the following code.

     /// \code

     ///   fc.greedyInit();

     ///   return fc.start();

     /// \endcode

     bool run() {

       greedyInit();

       return start();

     }

     /// @}

     /// \name Query Functions

     /// The result of the %Circulation algorithm can be obtained using

     /// these functions.

     /// \n

     /// Before the use of these functions,

     /// either run() or start() must be called.

     ///@{

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

     ///\sa checkBarrier()

     ///\sa run()

     template<class GT>

     void barrierMap(GT &bar) 

     {

       for(NodeIt n(_g);n!=INVALID;++n)

 	bar.set(n, (*_level)[n] >= _el);

     }  

     ///Returns true if the node is in the barrier

     ///Returns true if the node is in the barrier

     ///\sa barrierMap()

     bool barrier(const Node& node) 

     {

       return (*_level)[node] >= _el;

     }  

     /// \brief Returns the flow on the edge.

     ///

     /// Sets the \c flowMap to the flow on the edges. This method can

     /// be called after the second phase of algorithm.

     Value flow(const Edge& edge) const {

       return (*_flow)[edge];

     }

     /// @}

     /// \name Checker Functions

     /// The feasibility  of the results can be checked using

     /// these functions.

     /// \n

     /// Before the use of these functions,

     /// either run() or start() must be called.

     ///@{

     ///Check if the  \c flow is a feasible circulation

     bool checkFlow() {

       for(EdgeIt e(_g);e!=INVALID;++e)

 	if((*_flow)[e]<(*_lo)[e]||(*_flow)[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-=(*_flow)[e];

 	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];

 	  if(_tol.negative(dif)) return false;

 	}

       return true;

     }

     ///Check whether or not the last execution provides a barrier

     ///Check whether or not the last execution provides a barrier

     ///\sa barrier()

     bool checkBarrier() 

     {

       Value delta=0;

       for(NodeIt n(_g);n!=INVALID;++n)

 	if(barrier(n))

 	  delta-=(*_delta)[n];

       for(EdgeIt e(_g);e!=INVALID;++e)

 	{

 	  Node s=_g.source(e);

 	  Node t=_g.target(e);

 	  if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];

 	  else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];

 	}

       return _tol.negative(delta);

     }

     /// @}

   };

 }

 #endif