alpar@2375: /* -*- C++ -*-
alpar@2375:  * lemon/preflow_matching.h - Part of LEMON, a generic C++ optimization library
alpar@2375:  *
alpar@2375:  * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@2375:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@2375:  *
alpar@2375:  * Permission to use, modify and distribute this software is granted
alpar@2375:  * provided that this copyright notice appears in all copies. For
alpar@2375:  * precise terms see the accompanying LICENSE file.
alpar@2375:  *
alpar@2375:  * This software is provided "AS IS" with no warranty of any kind,
alpar@2375:  * express or implied, and with no claim as to its suitability for any
alpar@2375:  * purpose.
alpar@2375:  *
alpar@2375:  */
alpar@2375: 
alpar@2375: #ifndef LEMON_CIRCULATION_H
alpar@2375: #define LEMON_CIRCULATION_H
alpar@2375: 
alpar@2375: #include <lemon/graph_utils.h>
alpar@2375: #include <iostream>
alpar@2375: #include <queue>
alpar@2375: #include <lemon/tolerance.h>
alpar@2375: #include <lemon/elevator.h>
alpar@2375: 
alpar@2375: ///\ingroup flowalgs
alpar@2375: ///\file
alpar@2375: ///\brief Push-prelabel algorithm for finding a feasible circulation.
alpar@2375: ///
alpar@2375: namespace lemon {
alpar@2375: 
alpar@2375:   ///Preflow algorithm for the Network Circulation Problem.
alpar@2375:   
alpar@2375:   ///\ingroup flowalgs
alpar@2375:   ///This class implements a preflow algorithm
alpar@2375:   ///for the Network Circulation Problem.
alpar@2375:   ///The exact formulation of this problem is the following.
alpar@2375:   /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq delta(v)\quad \forall v\in V \f]
alpar@2375:   /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
alpar@2375:   ///
alpar@2375:   template<class Graph,
alpar@2375: 	   class Value,
alpar@2375: 	   class FlowMap=typename Graph::template EdgeMap<Value>,
alpar@2375: 	   class LCapMap=typename Graph::template EdgeMap<Value>,
alpar@2375: 	   class UCapMap=LCapMap,
alpar@2375: 	   class DeltaMap=typename Graph::template NodeMap<Value>
alpar@2375: 	    >
alpar@2375:   class Circulation {
alpar@2375:     typedef typename Graph::Node Node;
alpar@2375:     typedef typename Graph::NodeIt NodeIt;
alpar@2375:     typedef typename Graph::Edge Edge;
alpar@2375:     typedef typename Graph::EdgeIt EdgeIt;
alpar@2375:     typedef typename Graph::InEdgeIt InEdgeIt;
alpar@2375:     typedef typename Graph::OutEdgeIt OutEdgeIt;
alpar@2375:     
alpar@2375: 
alpar@2375:     const Graph &_g;
alpar@2375:     int _node_num;
alpar@2375: 
alpar@2375:     const LCapMap &_lo;
alpar@2375:     const UCapMap &_up;
alpar@2375:     const DeltaMap &_delta;
alpar@2375:     FlowMap &_x;
alpar@2375:     Tolerance<Value> _tol;
alpar@2375:     Elevator<Graph,typename Graph::Node> _levels;
alpar@2375:     typename Graph::template NodeMap<Value> _excess;
alpar@2375: 
alpar@2375:   public:
alpar@2375:     ///\e
alpar@2375:     Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
alpar@2375: 		const DeltaMap &delta,FlowMap &x) :
alpar@2375:       _g(g),
alpar@2375:       _node_num(countNodes(g)),
alpar@2375:       _lo(lo),_up(up),_delta(delta),_x(x),
alpar@2375:       _levels(g,_node_num),
alpar@2375:       _excess(g)
alpar@2375:     {
alpar@2375:     }
alpar@2375:     
alpar@2375:   private:
alpar@2375: 
alpar@2375:     void addExcess(Node n,Value v)
alpar@2375:     {
alpar@2375:       if(_tol.positive(_excess[n]+=v))
alpar@2375: 	{
alpar@2375: 	  if(!_levels.active(n)) _levels.activate(n);
alpar@2375: 	}
alpar@2375:       else if(_levels.active(n)) _levels.deactivate(n);
alpar@2375:     }
alpar@2375:     
alpar@2375:     void init() 
alpar@2375:     {
alpar@2375:      
alpar@2375:       _x=_lo;
alpar@2375: 
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=-_delta[n];
alpar@2375: 
alpar@2375:       for(EdgeIt e(_g);e!=INVALID;++e)
alpar@2375: 	{
alpar@2375: 	  _excess[_g.target(e)]+=_x[e];
alpar@2375: 	  _excess[_g.source(e)]-=_x[e];
alpar@2375: 	}
alpar@2375:      
alpar@2375:       _levels.initStart();
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	_levels.initAddItem(n);
alpar@2375:       _levels.initFinish();
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	if(_tol.positive(_excess[n]))
alpar@2375: 	  _levels.activate(n);
alpar@2375:     }
alpar@2375: 
alpar@2375:   public:
alpar@2375:     ///Check if \c x is a feasible circulation
alpar@2375:     template<class FT>
alpar@2375:     bool checkX(FT &x) {
alpar@2375:       for(EdgeIt e(_g);e!=INVALID;++e)
alpar@2375: 	if(x[e]<_lo[e]||x[e]>_up[e]) return false;
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	{
alpar@2375: 	  Value dif=_delta[n];
alpar@2375: 	  for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];
alpar@2375: 	  for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];
alpar@2375: 	  if(_tol.negative(dif)) return false;
alpar@2375: 	}
alpar@2375:       return true;
alpar@2375:     };
alpar@2375:     ///Check if the default \c x is a feasible circulation
alpar@2375:     bool checkX() { return checkX(_x); }
alpar@2375: 
alpar@2375:     ///Chech if \c bar is a real barrier
alpar@2375: 
alpar@2375:     ///Chech if \c bar is a real barrier
alpar@2375:     ///\sa barrier()
alpar@2375:     template<class GT>
alpar@2375:     bool checkBarrier(GT &bar) 
alpar@2375:     {
alpar@2375:       Value delta=0;
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	if(bar[n])
alpar@2375: 	  delta+=_delta[n];
alpar@2375:       for(EdgeIt e(_g);e!=INVALID;++e)
alpar@2375: 	{
alpar@2375: 	  Node s=_g.source(e);
alpar@2375: 	  Node t=_g.target(e);
alpar@2375: 	  if(bar[s]&&!bar[t]) delta+=_up[e];
alpar@2375: 	  else if(bar[t]&&!bar[s]) delta-=_lo[e];
alpar@2375: 	}
alpar@2375:       return _tol.negative(delta);
alpar@2375:     }  
alpar@2375:     ///Chech whether or not the last execution provides a barrier
alpar@2375: 
alpar@2375:     ///Chech whether or not the last execution provides a barrier
alpar@2375:     ///\sa barrier()
alpar@2375:     bool checkBarrier() 
alpar@2375:     {
alpar@2375:       typename Graph:: template NodeMap<bool> bar(_g);
alpar@2375:       barrier(bar);
alpar@2375:       return checkBarrier(bar);
alpar@2375:     }
alpar@2375:     ///Run the algorithm
alpar@2375: 
alpar@2375:     ///This function runs the algorithm.
alpar@2375:     ///\return This function returns -1 if it found a feasible circulation.
alpar@2375:     /// nonnegative values (including 0) mean that no feasible solution is
alpar@2375:     /// found. In this case the return value means an "empty level".
alpar@2375:     ///
alpar@2375:     ///\sa barrier()
alpar@2375:     int run() 
alpar@2375:     {
alpar@2375:       init();
alpar@2375:       
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	std::cerr<< _levels[n] << ' ';
alpar@2375:       std::cerr << std::endl;
alpar@2375: #endif      
alpar@2375:       Node act;
alpar@2375:       Node bact=INVALID;
alpar@2375:       Node last_activated=INVALID;
alpar@2375:       while((act=_levels.highestActive())!=INVALID) {
alpar@2375: 	int actlevel=_levels[act];
alpar@2375: 	int tlevel;
alpar@2375: 	int mlevel=_node_num;
alpar@2375: 	Value exc=_excess[act];
alpar@2375: 	Value fc;
alpar@2375: 	
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 	for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	  std::cerr<< _levels[n] << ' ';
alpar@2375: 	std::cerr << std::endl;
alpar@2375: 	std::cerr << "Process node " << _g.id(act)
alpar@2375: 		  << " on level " << actlevel
alpar@2375: 		  << " with excess " << exc
alpar@2375: 		  << std::endl;
alpar@2375: #endif
alpar@2375: 	for(OutEdgeIt e(_g,act);e!=INVALID; ++e)
alpar@2375: 	  if((fc=_up[e]-_x[e])>0)
alpar@2375: 	    if((tlevel=_levels[_g.target(e)])<actlevel)
alpar@2375: 	      if(fc<=exc) {
alpar@2375: 		_x[e]=_up[e];
alpar@2375: 		addExcess(_g.target(e),fc);
alpar@2375: 		exc-=fc;
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 		std::cerr << "  Push " << fc
alpar@2375: 			  << " toward " << _g.id(_g.target(e)) << std::endl;
alpar@2375: #endif
alpar@2375: 	      }
alpar@2375: 	      else {
alpar@2375: 		_x[e]+=exc;
alpar@2375: 		addExcess(_g.target(e),exc);
alpar@2375: 		//exc=0;
alpar@2375: 		_excess[act]=0;
alpar@2375: 		_levels.deactivate(act);
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 		std::cerr << "  Push " << exc
alpar@2375: 			  << " toward " << _g.id(_g.target(e)) << std::endl;
alpar@2375: 		std::cerr << "  Deactivate\n";
alpar@2375: #endif
alpar@2375: 		goto next_l;
alpar@2375: 	      }
alpar@2375: 	    else if(tlevel<mlevel) mlevel=tlevel;
alpar@2375: 	
alpar@2375: 	for(InEdgeIt e(_g,act);e!=INVALID; ++e)
alpar@2375: 	  if((fc=_x[e]-_lo[e])>0)
alpar@2375: 	    if((tlevel=_levels[_g.source(e)])<actlevel)
alpar@2375: 	      if(fc<=exc) {
alpar@2375: 		_x[e]=_lo[e];
alpar@2375: 		addExcess(_g.source(e),fc);
alpar@2375: 		exc-=fc;
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 		std::cerr << "  Push " << fc
alpar@2375: 			  << " toward " << _g.id(_g.source(e)) << std::endl;
alpar@2375: #endif
alpar@2375: 	      }
alpar@2375: 	      else {
alpar@2375: 		_x[e]-=exc;
alpar@2375: 		addExcess(_g.source(e),exc);
alpar@2375: 		//exc=0;
alpar@2375: 		_excess[act]=0;
alpar@2375: 		_levels.deactivate(act);
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 		std::cerr << "  Push " << exc
alpar@2375: 			  << " toward " << _g.id(_g.source(e)) << std::endl;
alpar@2375: 		std::cerr << "  Deactivate\n";
alpar@2375: #endif
alpar@2375: 		goto next_l;
alpar@2375:               }
alpar@2375: 	    else if(tlevel<mlevel) mlevel=tlevel;
alpar@2375:    
alpar@2375: 	_excess[act]=exc;
alpar@2375: 	if(!_tol.positive(exc)) _levels.deactivate(act);
alpar@2375: 	else if(mlevel==_node_num) {
alpar@2375: 	  _levels.liftHighestActiveTo(_node_num);
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 	  std::cerr << "  Lift to level " << _node_num << std::endl;
alpar@2375: #endif
alpar@2375: 	  return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
alpar@2375: 	}
alpar@2375: 	else {
alpar@2375: 	  _levels.liftHighestActiveTo(mlevel+1);
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375: 	  std::cerr << "  Lift to level " << mlevel+1 << std::endl;
alpar@2375: #endif
alpar@2375: 	  if(_levels.onLevel(actlevel)==0)
alpar@2375: 	    return actlevel;
alpar@2375: 	}
alpar@2375:       next_l:
alpar@2375: 	;
alpar@2375:       }
alpar@2375: #ifdef LEMON_CIRCULATION_DEBUG
alpar@2375:       std::cerr << "Feasible flow found.\n";
alpar@2375: #endif
alpar@2375:       return -1;
alpar@2375:     }
alpar@2375:     
alpar@2375:     ///Return a barrier
alpar@2375:     
alpar@2375:     ///Barrier is a set \e B of nodes for which
alpar@2375:     /// \f[ \sum_{v\in B}delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
alpar@2375:     ///holds. The existence of a set with this property prooves that a feasible
alpar@2375:     ///flow cannot exists.
alpar@2375:     ///\pre The run() must have been executed, and its return value was -1.
alpar@2375:     ///\sa checkBarrier()
alpar@2375:     ///\sa run()
alpar@2375:     template<class GT>
alpar@2375:     void barrier(GT &bar,int empty_level=-1) 
alpar@2375:     {
alpar@2375:       if(empty_level==-1)
alpar@2375: 	for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
alpar@2375:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@2375: 	bar[n] = _levels[n]>empty_level;
alpar@2375:     }  
alpar@2375: 
alpar@2375:   };
alpar@2375:   
alpar@2375: }
alpar@2375: 
alpar@2375: #endif