COIN-OR::LEMON - Graph Library

Changeset 2526:b7727edd44f2 in lemon-0.x for lemon/circulation.h


Ignore:
Timestamp:
11/28/07 18:51:02 (12 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3402
Message:

Redesign Circulation interface according to new flow interface
New greedy approach initialization

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/circulation.h

    r2512 r2526  
    3232namespace lemon {
    3333
    34   ///Preflow algorithm for the Network Circulation Problem.
     34  /// \brief Default traits class of Circulation class.
     35  ///
     36  /// Default traits class of Circulation class.
     37  /// \param _Graph Graph type.
     38  /// \param _CapacityMap Type of capacity map.
     39  template <typename _Graph, typename _LCapMap,
     40            typename _UCapMap, typename _DeltaMap>
     41  struct CirculationDefaultTraits {
     42
     43    /// \brief The graph type the algorithm runs on.
     44    typedef _Graph Graph;
     45
     46    /// \brief The type of the map that stores the circulation lower
     47    /// bound.
     48    ///
     49    /// The type of the map that stores the circulation lower bound.
     50    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
     51    typedef _LCapMap LCapMap;
     52
     53    /// \brief The type of the map that stores the circulation upper
     54    /// bound.
     55    ///
     56    /// The type of the map that stores the circulation upper bound.
     57    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
     58    typedef _UCapMap UCapMap;
     59
     60    /// \brief The type of the map that stores the upper bound of
     61    /// node excess.
     62    ///
     63    /// The type of the map that stores the lower bound of node
     64    /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
     65    /// concept.
     66    typedef _DeltaMap DeltaMap;
     67
     68    /// \brief The type of the length of the edges.
     69    typedef typename DeltaMap::Value Value;
     70
     71    /// \brief The map type that stores the flow values.
     72    ///
     73    /// The map type that stores the flow values.
     74    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
     75    typedef typename Graph::template EdgeMap<Value> FlowMap;
     76
     77    /// \brief Instantiates a FlowMap.
     78    ///
     79    /// This function instantiates a \ref FlowMap.
     80    /// \param graph The graph, to which we would like to define the flow map.
     81    static FlowMap* createFlowMap(const Graph& graph) {
     82      return new FlowMap(graph);
     83    }
     84
     85    /// \brief The eleavator type used by Circulation algorithm.
     86    ///
     87    /// The elevator type used by Circulation algorithm.
     88    ///
     89    /// \sa Elevator
     90    /// \sa LinkedElevator
     91    typedef Elevator<Graph, typename Graph::Node> Elevator;
     92   
     93    /// \brief Instantiates an Elevator.
     94    ///
     95    /// This function instantiates a \ref Elevator.
     96    /// \param graph The graph, to which we would like to define the elevator.
     97    /// \param max_level The maximum level of the elevator.
     98    static Elevator* createElevator(const Graph& graph, int max_level) {
     99      return new Elevator(graph, max_level);
     100    }
     101
     102    /// \brief The tolerance used by the algorithm
     103    ///
     104    /// The tolerance used by the algorithm to handle inexact computation.
     105    typedef Tolerance<Value> Tolerance;
     106
     107  };
     108 
     109  ///Push-relabel algorithm for the Network Circulation Problem.
    35110 
    36111  ///\ingroup max_flow
    37   ///This class implements a preflow algorithm
     112  ///This class implements a push-relabel algorithm
    38113  ///for the Network Circulation Problem.
    39114  ///The exact formulation of this problem is the following.
     
    41116  /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
    42117  ///
    43   template<class Graph,
    44            class Value,
    45            class FlowMap=typename Graph::template EdgeMap<Value>,
    46            class LCapMap=typename Graph::template EdgeMap<Value>,
    47            class UCapMap=LCapMap,
    48            class DeltaMap=typename Graph::template NodeMap<Value>
    49             >
     118  template<class _Graph,
     119           class _LCapMap=typename _Graph::template EdgeMap<int>,
     120           class _UCapMap=_LCapMap,
     121           class _DeltaMap=typename _Graph::template NodeMap<
     122             typename _UCapMap::Value>,
     123           class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
     124                                                  _UCapMap, _DeltaMap> >
    50125  class Circulation {
    51     typedef typename Graph::Node Node;
    52     typedef typename Graph::NodeIt NodeIt;
    53     typedef typename Graph::Edge Edge;
    54     typedef typename Graph::EdgeIt EdgeIt;
    55     typedef typename Graph::InEdgeIt InEdgeIt;
    56     typedef typename Graph::OutEdgeIt OutEdgeIt;
    57    
     126
     127    typedef _Traits Traits;
     128    typedef typename Traits::Graph Graph;
     129    GRAPH_TYPEDEFS(typename Graph);
     130
     131    typedef typename Traits::Value Value;
     132
     133    typedef typename Traits::LCapMap LCapMap;
     134    typedef typename Traits::UCapMap UCapMap;
     135    typedef typename Traits::DeltaMap DeltaMap;
     136    typedef typename Traits::FlowMap FlowMap;
     137    typedef typename Traits::Elevator Elevator;
     138    typedef typename Traits::Tolerance Tolerance;
     139
     140    typedef typename Graph::template NodeMap<Value> ExcessMap;
    58141
    59142    const Graph &_g;
    60143    int _node_num;
    61144
    62     const LCapMap &_lo;
    63     const UCapMap &_up;
    64     const DeltaMap &_delta;
    65     FlowMap &_x;
    66     Tolerance<Value> _tol;
    67     Elevator<Graph,typename Graph::Node> _levels;
    68     typename Graph::template NodeMap<Value> _excess;
     145    const LCapMap *_lo;
     146    const UCapMap *_up;
     147    const DeltaMap *_delta;
     148
     149    FlowMap *_flow;
     150    bool _local_flow;
     151
     152    Elevator* _level;
     153    bool _local_level;
     154
     155    ExcessMap* _excess;
     156
     157    Tolerance _tol;
     158    int _el;
    69159
    70160  public:
    71     ///\e
    72     Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
    73                 const DeltaMap &delta,FlowMap &x) :
    74       _g(g),
    75       _node_num(countNodes(g)),
    76       _lo(lo),_up(up),_delta(delta),_x(x),
    77       _levels(g,_node_num),
    78       _excess(g)
    79     {
    80     }
    81    
     161
     162    typedef Circulation Create;
     163
     164    ///\name Named template parameters
     165
     166    ///@{
     167
     168    template <typename _FlowMap>
     169    struct DefFlowMapTraits : public Traits {
     170      typedef _FlowMap FlowMap;
     171      static FlowMap *createFlowMap(const Graph&) {
     172        throw UninitializedParameter();
     173      }
     174    };
     175
     176    /// \brief \ref named-templ-param "Named parameter" for setting
     177    /// FlowMap type
     178    ///
     179    /// \ref named-templ-param "Named parameter" for setting FlowMap
     180    /// type
     181    template <typename _FlowMap>
     182    struct DefFlowMap
     183      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
     184                           DefFlowMapTraits<_FlowMap> > {
     185      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
     186                          DefFlowMapTraits<_FlowMap> > Create;
     187    };
     188
     189    template <typename _Elevator>
     190    struct DefElevatorTraits : public Traits {
     191      typedef _Elevator Elevator;
     192      static Elevator *createElevator(const Graph&, int) {
     193        throw UninitializedParameter();
     194      }
     195    };
     196
     197    /// \brief \ref named-templ-param "Named parameter" for setting
     198    /// Elevator type
     199    ///
     200    /// \ref named-templ-param "Named parameter" for setting Elevator
     201    /// type
     202    template <typename _Elevator>
     203    struct DefElevator
     204      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
     205                           DefElevatorTraits<_Elevator> > {
     206      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
     207                          DefElevatorTraits<_Elevator> > Create;
     208    };
     209
     210    template <typename _Elevator>
     211    struct DefStandardElevatorTraits : public Traits {
     212      typedef _Elevator Elevator;
     213      static Elevator *createElevator(const Graph& graph, int max_level) {
     214        return new Elevator(graph, max_level);
     215      }
     216    };
     217
     218    /// \brief \ref named-templ-param "Named parameter" for setting
     219    /// Elevator type
     220    ///
     221    /// \ref named-templ-param "Named parameter" for setting Elevator
     222    /// type. The Elevator should be standard constructor interface, ie.
     223    /// the graph and the maximum level should be passed to it.
     224    template <typename _Elevator>
     225    struct DefStandardElevator
     226      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
     227                       DefStandardElevatorTraits<_Elevator> > {
     228      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
     229                      DefStandardElevatorTraits<_Elevator> > Create;
     230    };   
     231
     232    /// @}
     233
     234    /// The constructor of the class.
     235
     236    /// The constructor of the class.
     237    /// \param g The directed graph the algorithm runs on.
     238    /// \param lo The lower bound capacity of the edges.
     239    /// \param up The upper bound capacity of the edges.
     240    /// \param delta The lower bound on node excess.
     241    Circulation(const Graph &g,const LCapMap &lo,
     242                const UCapMap &up,const DeltaMap &delta)
     243      : _g(g), _node_num(),
     244        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
     245        _level(0), _local_level(false), _excess(0), _el() {}
     246
     247    /// Destrcutor.
     248    ~Circulation() {
     249      destroyStructures();
     250    }
     251
    82252  private:
    83253
    84     void addExcess(Node n,Value v)
    85     {
    86       if(_tol.positive(_excess[n]+=v))
    87         {
    88           if(!_levels.active(n)) _levels.activate(n);
    89         }
    90       else if(_levels.active(n)) _levels.deactivate(n);
    91     }
    92    
     254    void createStructures() {
     255      _node_num = _el = countNodes(_g);
     256
     257      if (!_flow) {
     258        _flow = Traits::createFlowMap(_g);
     259        _local_flow = true;
     260      }
     261      if (!_level) {
     262        _level = Traits::createElevator(_g, _node_num);
     263        _local_level = true;
     264      }
     265      if (!_excess) {
     266        _excess = new ExcessMap(_g);
     267      }
     268    }
     269
     270    void destroyStructures() {
     271      if (_local_flow) {
     272        delete _flow;
     273      }
     274      if (_local_level) {
     275        delete _level;
     276      }
     277      if (_excess) {
     278        delete _excess;
     279      }
     280    }
     281
     282  public:
     283
     284    /// Sets the lower bound capacity map.
     285
     286    /// Sets the lower bound capacity map.
     287    /// \return \c (*this)
     288    Circulation& lowerCapMap(const LCapMap& map) {
     289      _lo = &map;
     290      return *this;
     291    }
     292
     293    /// Sets the upper bound capacity map.
     294
     295    /// Sets the upper bound capacity map.
     296    /// \return \c (*this)
     297    Circulation& upperCapMap(const LCapMap& map) {
     298      _up = &map;
     299      return *this;
     300    }
     301
     302    /// Sets the lower bound map on excess.
     303
     304    /// Sets the lower bound map on excess.
     305    /// \return \c (*this)
     306    Circulation& deltaMap(const DeltaMap& map) {
     307      _delta = &map;
     308      return *this;
     309    }
     310
     311    /// Sets the flow map.
     312
     313    /// Sets the flow map.
     314    /// \return \c (*this)
     315    Circulation& flowMap(FlowMap& map) {
     316      if (_local_flow) {
     317        delete _flow;
     318        _local_flow = false;
     319      }
     320      _flow = &map;
     321      return *this;
     322    }
     323
     324    /// Returns the flow map.
     325
     326    /// \return The flow map.
     327    ///
     328    const FlowMap& flowMap() {
     329      return *_flow;
     330    }
     331
     332    /// Sets the elevator.
     333
     334    /// Sets the elevator.
     335    /// \return \c (*this)
     336    Circulation& elevator(Elevator& elevator) {
     337      if (_local_level) {
     338        delete _level;
     339        _local_level = false;
     340      }
     341      _level = &elevator;
     342      return *this;
     343    }
     344
     345    /// Returns the elevator.
     346
     347    /// \return The elevator.
     348    ///
     349    const Elevator& elevator() {
     350      return *_level;
     351    }
     352   
     353    /// Sets the tolerance used by algorithm.
     354
     355    /// Sets the tolerance used by algorithm.
     356    ///
     357    Circulation& tolerance(const Tolerance& tolerance) const {
     358      _tol = tolerance;
     359      return *this;
     360    }
     361
     362    /// Returns the tolerance used by algorithm.
     363
     364    /// Returns the tolerance used by algorithm.
     365    ///
     366    const Tolerance& tolerance() const {
     367      return tolerance;
     368    }
     369   
     370    /// \name Execution control The simplest way to execute the
     371    /// algorithm is to use one of the member functions called \c
     372    /// run(). 
     373    /// \n
     374    /// If you need more control on initial solution or execution then
     375    /// you have to call one \ref init() function and then the start()
     376    /// function.
     377   
     378    ///@{
     379 
     380    /// Initializes the internal data structures.
     381
     382    /// Initializes the internal data structures. This function sets
     383    /// all flow values to the lower bound. 
     384    /// \return This function returns false if the initialization
     385    /// process found a barrier.
    93386    void init()
    94387    {
     388      createStructures();
     389
     390      for(NodeIt n(_g);n!=INVALID;++n) {
     391        _excess->set(n, (*_delta)[n]);
     392      }
    95393     
    96       _x=_lo;
    97 
    98       for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=_delta[n];
    99 
     394      for (EdgeIt e(_g);e!=INVALID;++e) {
     395        _flow->set(e, (*_lo)[e]);
     396        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
     397        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
     398      }
     399
     400      typename Graph::template NodeMap<bool> reached(_g, false);
     401
     402
     403      // global relabeling tested, but in general case it provides
     404      // worse performance for random graphs
     405      _level->initStart();
     406      for(NodeIt n(_g);n!=INVALID;++n)
     407        _level->initAddItem(n);
     408      _level->initFinish();
     409      for(NodeIt n(_g);n!=INVALID;++n)
     410        if(_tol.positive((*_excess)[n]))
     411          _level->activate(n);
     412    }
     413
     414    /// Initializes the internal data structures.
     415
     416    /// Initializes the internal data structures. This functions uses
     417    /// greedy approach to construct the initial solution.
     418    void greedyInit()
     419    {
     420      createStructures();
     421     
     422      for(NodeIt n(_g);n!=INVALID;++n) {
     423        _excess->set(n, (*_delta)[n]);
     424      }
     425     
     426      for (EdgeIt e(_g);e!=INVALID;++e) {
     427        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
     428          _flow->set(e, (*_up)[e]);
     429          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
     430          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
     431        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
     432          _flow->set(e, (*_lo)[e]);
     433          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
     434          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
     435        } else {
     436          Value fc = -(*_excess)[_g.target(e)];
     437          _flow->set(e, fc);
     438          _excess->set(_g.target(e), 0);
     439          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
     440        }
     441      }
     442     
     443      _level->initStart();
     444      for(NodeIt n(_g);n!=INVALID;++n)
     445        _level->initAddItem(n);
     446      _level->initFinish();
     447      for(NodeIt n(_g);n!=INVALID;++n)
     448        if(_tol.positive((*_excess)[n]))
     449          _level->activate(n);
     450    }
     451
     452    ///Starts the algorithm
     453
     454    ///This function starts the algorithm.
     455    ///\return This function returns true if it found a feasible circulation.
     456    ///
     457    ///\sa barrier()
     458    bool start()
     459    {
     460     
     461      Node act;
     462      Node bact=INVALID;
     463      Node last_activated=INVALID;
     464      while((act=_level->highestActive())!=INVALID) {
     465        int actlevel=(*_level)[act];
     466        int mlevel=_node_num;
     467        Value exc=(*_excess)[act];
     468       
     469        for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
     470          Node v = _g.target(e);
     471          Value fc=(*_up)[e]-(*_flow)[e];
     472          if(!_tol.positive(fc)) continue;
     473          if((*_level)[v]<actlevel) {
     474            if(!_tol.less(fc, exc)) {
     475              _flow->set(e, (*_flow)[e] + exc);
     476              _excess->set(v, (*_excess)[v] + exc);
     477              if(!_level->active(v) && _tol.positive((*_excess)[v]))
     478                _level->activate(v);
     479              _excess->set(act,0);
     480              _level->deactivate(act);
     481              goto next_l;
     482            }
     483            else {
     484              _flow->set(e, (*_up)[e]);
     485              _excess->set(v, (*_excess)[v] + exc);
     486              if(!_level->active(v) && _tol.positive((*_excess)[v]))
     487                _level->activate(v);
     488              exc-=fc;
     489            }
     490          }
     491          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
     492        }
     493        for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
     494          Node v = _g.source(e);
     495          Value fc=(*_flow)[e]-(*_lo)[e];
     496          if(!_tol.positive(fc)) continue;
     497          if((*_level)[v]<actlevel) {
     498            if(!_tol.less(fc, exc)) {
     499              _flow->set(e, (*_flow)[e] - exc);
     500              _excess->set(v, (*_excess)[v] + exc);
     501              if(!_level->active(v) && _tol.positive((*_excess)[v]))
     502                _level->activate(v);
     503              _excess->set(act,0);
     504              _level->deactivate(act);
     505              goto next_l;
     506            }
     507            else {
     508              _flow->set(e, (*_lo)[e]);
     509              _excess->set(v, (*_excess)[v] + fc);
     510              if(!_level->active(v) && _tol.positive((*_excess)[v]))
     511                _level->activate(v);
     512              exc-=fc;
     513            }
     514          }
     515          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
     516        }
     517   
     518        _excess->set(act, exc);
     519        if(!_tol.positive(exc)) _level->deactivate(act);
     520        else if(mlevel==_node_num) {
     521          _level->liftHighestActiveToTop();
     522          _el = _node_num;
     523          return false;
     524        }
     525        else {
     526          _level->liftHighestActive(mlevel+1);
     527          if(_level->onLevel(actlevel)==0) {
     528            _el = actlevel;
     529            return false;
     530          }
     531        }
     532      next_l:
     533        ;
     534      }
     535      return true;
     536    }
     537
     538    /// Runs the circulation algorithm. 
     539
     540    /// Runs the circulation algorithm.
     541    /// \note fc.run() is just a shortcut of the following code.
     542    /// \code
     543    ///   fc.greedyInit();
     544    ///   return fc.start();
     545    /// \endcode
     546    bool run() {
     547      greedyInit();
     548      return start();
     549    }
     550
     551    /// @}
     552
     553    /// \name Query Functions
     554    /// The result of the %Circulation algorithm can be obtained using
     555    /// these functions.
     556    /// \n
     557    /// Before the use of these functions,
     558    /// either run() or start() must be called.
     559   
     560    ///@{
     561   
     562    ///Returns a barrier
     563   
     564    ///Barrier is a set \e B of nodes for which
     565    /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
     566    ///holds. The existence of a set with this property prooves that a feasible
     567    ///flow cannot exists.
     568    ///\sa checkBarrier()
     569    ///\sa run()
     570    template<class GT>
     571    void barrierMap(GT &bar)
     572    {
     573      for(NodeIt n(_g);n!=INVALID;++n)
     574        bar.set(n, (*_level)[n] >= _el);
     575    } 
     576
     577    ///Returns true if the node is in the barrier
     578
     579    ///Returns true if the node is in the barrier
     580    ///\sa barrierMap()
     581    bool barrier(const Node& node)
     582    {
     583      return (*_level)[node] >= _el;
     584    } 
     585
     586    /// \brief Returns the flow on the edge.
     587    ///
     588    /// Sets the \c flowMap to the flow on the edges. This method can
     589    /// be called after the second phase of algorithm.
     590    Value flow(const Edge& edge) const {
     591      return (*_flow)[edge];
     592    }
     593
     594    /// @}
     595
     596    /// \name Checker Functions
     597    /// The feasibility  of the results can be checked using
     598    /// these functions.
     599    /// \n
     600    /// Before the use of these functions,
     601    /// either run() or start() must be called.
     602   
     603    ///@{
     604
     605    ///Check if the  \c flow is a feasible circulation
     606    bool checkFlow() {
    100607      for(EdgeIt e(_g);e!=INVALID;++e)
    101         {
    102           _excess[_g.target(e)]+=_x[e];
    103           _excess[_g.source(e)]-=_x[e];
    104         }
    105      
    106       _levels.initStart();
    107       for(NodeIt n(_g);n!=INVALID;++n)
    108         _levels.initAddItem(n);
    109       _levels.initFinish();
    110       for(NodeIt n(_g);n!=INVALID;++n)
    111         if(_tol.positive(_excess[n]))
    112           _levels.activate(n);
    113     }
    114 
    115   public:
    116     ///Check if \c x is a feasible circulation
    117     template<class FT>
    118     bool checkX(FT &x) {
    119       for(EdgeIt e(_g);e!=INVALID;++e)
    120         if(x[e]<_lo[e]||x[e]>_up[e]) return false;
     608        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
    121609      for(NodeIt n(_g);n!=INVALID;++n)
    122610        {
    123           Value dif=-_delta[n];
    124           for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];
    125           for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];
     611          Value dif=-(*_delta)[n];
     612          for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
     613          for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
    126614          if(_tol.negative(dif)) return false;
    127615        }
    128616      return true;
    129     };
    130     ///Check if the default \c x is a feasible circulation
    131     bool checkX() { return checkX(_x); }
    132 
    133     ///Check if \c bar is a real barrier
    134 
    135     ///Check if \c bar is a real barrier
     617    }
     618
     619    ///Check whether or not the last execution provides a barrier
     620
     621    ///Check whether or not the last execution provides a barrier
    136622    ///\sa barrier()
    137     template<class GT>
    138     bool checkBarrier(GT &bar)
     623    bool checkBarrier()
    139624    {
    140625      Value delta=0;
    141626      for(NodeIt n(_g);n!=INVALID;++n)
    142         if(bar[n])
    143           delta-=_delta[n];
     627        if(barrier(n))
     628          delta-=(*_delta)[n];
    144629      for(EdgeIt e(_g);e!=INVALID;++e)
    145630        {
    146631          Node s=_g.source(e);
    147632          Node t=_g.target(e);
    148           if(bar[s]&&!bar[t]) delta+=_up[e];
    149           else if(bar[t]&&!bar[s]) delta-=_lo[e];
     633          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
     634          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
    150635        }
    151636      return _tol.negative(delta);
    152     } 
    153     ///Check whether or not the last execution provides a barrier
    154 
    155     ///Check whether or not the last execution provides a barrier
    156     ///\sa barrier()
    157     bool checkBarrier()
    158     {
    159       typename Graph:: template NodeMap<bool> bar(_g);
    160       barrier(bar);
    161       return checkBarrier(bar);
    162     }
    163     ///Run the algorithm
    164 
    165     ///This function runs the algorithm.
    166     ///\return This function returns -1 if it found a feasible circulation.
    167     /// nonnegative values (including 0) mean that no feasible solution is
    168     /// found. In this case the return value means an "empty level".
    169     ///
    170     ///\sa barrier()
    171     int run()
    172     {
    173       init();
    174      
    175 #ifdef LEMON_CIRCULATION_DEBUG
    176       for(NodeIt n(_g);n!=INVALID;++n)
    177         std::cerr<< _levels[n] << ' ';
    178       std::cerr << std::endl;
    179 #endif     
    180       Node act;
    181       Node bact=INVALID;
    182       Node last_activated=INVALID;
    183       while((act=_levels.highestActive())!=INVALID) {
    184         int actlevel=_levels[act];
    185         int tlevel;
    186         int mlevel=_node_num;
    187         Value exc=_excess[act];
    188         Value fc;
    189        
    190 #ifdef LEMON_CIRCULATION_DEBUG
    191         for(NodeIt n(_g);n!=INVALID;++n)
    192           std::cerr<< _levels[n] << ' ';
    193         std::cerr << std::endl;
    194         std::cerr << "Process node " << _g.id(act)
    195                   << " on level " << actlevel
    196                   << " with excess " << exc
    197                   << std::endl;
    198 #endif
    199         for(OutEdgeIt e(_g,act);e!=INVALID; ++e)
    200           if((fc=_up[e]-_x[e])>0)
    201             if((tlevel=_levels[_g.target(e)])<actlevel)
    202               if(fc<=exc) {
    203                 _x[e]=_up[e];
    204                 addExcess(_g.target(e),fc);
    205                 exc-=fc;
    206 #ifdef LEMON_CIRCULATION_DEBUG
    207                 std::cerr << "  Push " << fc
    208                           << " toward " << _g.id(_g.target(e)) << std::endl;
    209 #endif
    210               }
    211               else {
    212                 _x[e]+=exc;
    213                 addExcess(_g.target(e),exc);
    214                 //exc=0;
    215                 _excess[act]=0;
    216                 _levels.deactivate(act);
    217 #ifdef LEMON_CIRCULATION_DEBUG
    218                 std::cerr << "  Push " << exc
    219                           << " toward " << _g.id(_g.target(e)) << std::endl;
    220                 std::cerr << "  Deactivate\n";
    221 #endif
    222                 goto next_l;
    223               }
    224             else if(tlevel<mlevel) mlevel=tlevel;
    225        
    226         for(InEdgeIt e(_g,act);e!=INVALID; ++e)
    227           if((fc=_x[e]-_lo[e])>0)
    228             if((tlevel=_levels[_g.source(e)])<actlevel)
    229               if(fc<=exc) {
    230                 _x[e]=_lo[e];
    231                 addExcess(_g.source(e),fc);
    232                 exc-=fc;
    233 #ifdef LEMON_CIRCULATION_DEBUG
    234                 std::cerr << "  Push " << fc
    235                           << " toward " << _g.id(_g.source(e)) << std::endl;
    236 #endif
    237               }
    238               else {
    239                 _x[e]-=exc;
    240                 addExcess(_g.source(e),exc);
    241                 //exc=0;
    242                 _excess[act]=0;
    243                 _levels.deactivate(act);
    244 #ifdef LEMON_CIRCULATION_DEBUG
    245                 std::cerr << "  Push " << exc
    246                           << " toward " << _g.id(_g.source(e)) << std::endl;
    247                 std::cerr << "  Deactivate\n";
    248 #endif
    249                 goto next_l;
    250               }
    251             else if(tlevel<mlevel) mlevel=tlevel;
    252    
    253         _excess[act]=exc;
    254         if(!_tol.positive(exc)) _levels.deactivate(act);
    255         else if(mlevel==_node_num) {
    256           _levels.liftHighestActiveToTop();
    257 #ifdef LEMON_CIRCULATION_DEBUG
    258           std::cerr << "  Lift to level " << _node_num << std::endl;
    259 #endif
    260           return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
    261         }
    262         else {
    263           _levels.liftHighestActive(mlevel+1);
    264 #ifdef LEMON_CIRCULATION_DEBUG
    265           std::cerr << "  Lift to level " << mlevel+1 << std::endl;
    266 #endif
    267           if(_levels.onLevel(actlevel)==0)
    268             return actlevel;
    269         }
    270       next_l:
    271         ;
    272       }
    273 #ifdef LEMON_CIRCULATION_DEBUG
    274       std::cerr << "Feasible flow found.\n";
    275 #endif
    276       return -1;
    277     }
    278    
    279     ///Return a barrier
    280    
    281     ///Barrier is a set \e B of nodes for which
    282     /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
    283     ///holds. The existence of a set with this property prooves that a feasible
    284     ///flow cannot exists.
    285     ///\pre The run() must have been executed, and its return value was -1.
    286     ///\sa checkBarrier()
    287     ///\sa run()
    288     template<class GT>
    289     void barrier(GT &bar,int empty_level=-1)
    290     {
    291       if(empty_level==-1)
    292         for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
    293       for(NodeIt n(_g);n!=INVALID;++n)
    294         bar[n] = _levels[n]>empty_level;
    295     } 
     637    }
     638
     639    /// @}
    296640
    297641  };
Note: See TracChangeset for help on using the changeset viewer.