source: lemon-0.x/lemon/circulation.h @ 2376:0ed45a6c74b1

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