lemon/circulation.h
author Akos Ladanyi <ladanyi@tmit.bme.hu>
Mon, 03 Nov 2008 11:59:54 +0000
changeset 522 22f932bbb305
parent 420 6a2a33ad261b
child 492 9605e051942f
permissions -rw-r--r--
Test for euler.h (#65)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_CIRCULATION_H
    20 #define LEMON_CIRCULATION_H
    21 
    22 #include <lemon/tolerance.h>
    23 #include <lemon/elevator.h>
    24 
    25 ///\ingroup max_flow
    26 ///\file
    27 ///\brief Push-relabel algorithm for finding a feasible circulation.
    28 ///
    29 namespace lemon {
    30 
    31   /// \brief Default traits class of Circulation class.
    32   ///
    33   /// Default traits class of Circulation class.
    34   /// \tparam _Diraph Digraph type.
    35   /// \tparam _LCapMap Lower bound capacity map type.
    36   /// \tparam _UCapMap Upper bound capacity map type.
    37   /// \tparam _DeltaMap Delta map type.
    38   template <typename _Diraph, typename _LCapMap,
    39             typename _UCapMap, typename _DeltaMap>
    40   struct CirculationDefaultTraits {
    41 
    42     /// \brief The type of the digraph the algorithm runs on.
    43     typedef _Diraph Digraph;
    44 
    45     /// \brief The type of the map that stores the circulation lower
    46     /// bound.
    47     ///
    48     /// The type of the map that stores the circulation lower bound.
    49     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    50     typedef _LCapMap LCapMap;
    51 
    52     /// \brief The type of the map that stores the circulation upper
    53     /// bound.
    54     ///
    55     /// The type of the map that stores the circulation upper bound.
    56     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    57     typedef _UCapMap UCapMap;
    58 
    59     /// \brief The type of the map that stores the lower bound for
    60     /// the supply of the nodes.
    61     ///
    62     /// The type of the map that stores the lower bound for the supply
    63     /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
    64     /// concept.
    65     typedef _DeltaMap DeltaMap;
    66 
    67     /// \brief The type of the flow values.
    68     typedef typename DeltaMap::Value Value;
    69 
    70     /// \brief The type of the map that stores the flow values.
    71     ///
    72     /// The type of the map that stores the flow values.
    73     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
    74     typedef typename Digraph::template ArcMap<Value> FlowMap;
    75 
    76     /// \brief Instantiates a FlowMap.
    77     ///
    78     /// This function instantiates a \ref FlowMap.
    79     /// \param digraph The digraph, to which we would like to define
    80     /// the flow map.
    81     static FlowMap* createFlowMap(const Digraph& digraph) {
    82       return new FlowMap(digraph);
    83     }
    84 
    85     /// \brief The elevator type used by the algorithm.
    86     ///
    87     /// The elevator type used by the algorithm.
    88     ///
    89     /// \sa Elevator
    90     /// \sa LinkedElevator
    91     typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
    92 
    93     /// \brief Instantiates an Elevator.
    94     ///
    95     /// This function instantiates an \ref Elevator.
    96     /// \param digraph The digraph, to which we would like to define
    97     /// the elevator.
    98     /// \param max_level The maximum level of the elevator.
    99     static Elevator* createElevator(const Digraph& digraph, int max_level) {
   100       return new Elevator(digraph, max_level);
   101     }
   102 
   103     /// \brief The tolerance used by the algorithm
   104     ///
   105     /// The tolerance used by the algorithm to handle inexact computation.
   106     typedef lemon::Tolerance<Value> Tolerance;
   107 
   108   };
   109 
   110   /**
   111      \brief Push-relabel algorithm for the network circulation problem.
   112 
   113      \ingroup max_flow
   114      This class implements a push-relabel algorithm for the network
   115      circulation problem.
   116      It is to find a feasible circulation when lower and upper bounds
   117      are given for the flow values on the arcs and lower bounds
   118      are given for the supply values of the nodes.
   119 
   120      The exact formulation of this problem is the following.
   121      Let \f$G=(V,A)\f$ be a digraph,
   122      \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
   123      \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
   124      \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
   125      \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
   126      \geq delta(v) \quad \forall v\in V, \f]
   127      \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
   128      \note \f$delta(v)\f$ specifies a lower bound for the supply of node
   129      \f$v\f$. It can be either positive or negative, however note that
   130      \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
   131      have a feasible solution.
   132 
   133      \note A special case of this problem is when
   134      \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
   135      will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
   136      Thus a feasible solution for the
   137      \ref min_cost_flow "minimum cost flow" problem can be calculated
   138      in this way.
   139 
   140      \tparam _Digraph The type of the digraph the algorithm runs on.
   141      \tparam _LCapMap The type of the lower bound capacity map. The default
   142      map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
   143      \tparam _UCapMap The type of the upper bound capacity map. The default
   144      map type is \c _LCapMap.
   145      \tparam _DeltaMap The type of the map that stores the lower bound
   146      for the supply of the nodes. The default map type is
   147      \c _Digraph::ArcMap<_UCapMap::Value>.
   148   */
   149 #ifdef DOXYGEN
   150 template< typename _Digraph,
   151           typename _LCapMap,
   152           typename _UCapMap,
   153           typename _DeltaMap,
   154           typename _Traits >
   155 #else
   156 template< typename _Digraph,
   157           typename _LCapMap = typename _Digraph::template ArcMap<int>,
   158           typename _UCapMap = _LCapMap,
   159           typename _DeltaMap = typename _Digraph::
   160                                template NodeMap<typename _UCapMap::Value>,
   161           typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
   162                                                     _UCapMap, _DeltaMap> >
   163 #endif
   164   class Circulation {
   165   public:
   166 
   167     ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
   168     typedef _Traits Traits;
   169     ///The type of the digraph the algorithm runs on.
   170     typedef typename Traits::Digraph Digraph;
   171     ///The type of the flow values.
   172     typedef typename Traits::Value Value;
   173 
   174     /// The type of the lower bound capacity map.
   175     typedef typename Traits::LCapMap LCapMap;
   176     /// The type of the upper bound capacity map.
   177     typedef typename Traits::UCapMap UCapMap;
   178     /// \brief The type of the map that stores the lower bound for
   179     /// the supply of the nodes.
   180     typedef typename Traits::DeltaMap DeltaMap;
   181     ///The type of the flow map.
   182     typedef typename Traits::FlowMap FlowMap;
   183 
   184     ///The type of the elevator.
   185     typedef typename Traits::Elevator Elevator;
   186     ///The type of the tolerance.
   187     typedef typename Traits::Tolerance Tolerance;
   188 
   189   private:
   190 
   191     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   192 
   193     const Digraph &_g;
   194     int _node_num;
   195 
   196     const LCapMap *_lo;
   197     const UCapMap *_up;
   198     const DeltaMap *_delta;
   199 
   200     FlowMap *_flow;
   201     bool _local_flow;
   202 
   203     Elevator* _level;
   204     bool _local_level;
   205 
   206     typedef typename Digraph::template NodeMap<Value> ExcessMap;
   207     ExcessMap* _excess;
   208 
   209     Tolerance _tol;
   210     int _el;
   211 
   212   public:
   213 
   214     typedef Circulation Create;
   215 
   216     ///\name Named Template Parameters
   217 
   218     ///@{
   219 
   220     template <typename _FlowMap>
   221     struct SetFlowMapTraits : public Traits {
   222       typedef _FlowMap FlowMap;
   223       static FlowMap *createFlowMap(const Digraph&) {
   224         LEMON_ASSERT(false, "FlowMap is not initialized");
   225         return 0; // ignore warnings
   226       }
   227     };
   228 
   229     /// \brief \ref named-templ-param "Named parameter" for setting
   230     /// FlowMap type
   231     ///
   232     /// \ref named-templ-param "Named parameter" for setting FlowMap
   233     /// type.
   234     template <typename _FlowMap>
   235     struct SetFlowMap
   236       : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
   237                            SetFlowMapTraits<_FlowMap> > {
   238       typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
   239                           SetFlowMapTraits<_FlowMap> > Create;
   240     };
   241 
   242     template <typename _Elevator>
   243     struct SetElevatorTraits : public Traits {
   244       typedef _Elevator Elevator;
   245       static Elevator *createElevator(const Digraph&, int) {
   246         LEMON_ASSERT(false, "Elevator is not initialized");
   247         return 0; // ignore warnings
   248       }
   249     };
   250 
   251     /// \brief \ref named-templ-param "Named parameter" for setting
   252     /// Elevator type
   253     ///
   254     /// \ref named-templ-param "Named parameter" for setting Elevator
   255     /// type. If this named parameter is used, then an external
   256     /// elevator object must be passed to the algorithm using the
   257     /// \ref elevator(Elevator&) "elevator()" function before calling
   258     /// \ref run() or \ref init().
   259     /// \sa SetStandardElevator
   260     template <typename _Elevator>
   261     struct SetElevator
   262       : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
   263                            SetElevatorTraits<_Elevator> > {
   264       typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
   265                           SetElevatorTraits<_Elevator> > Create;
   266     };
   267 
   268     template <typename _Elevator>
   269     struct SetStandardElevatorTraits : public Traits {
   270       typedef _Elevator Elevator;
   271       static Elevator *createElevator(const Digraph& digraph, int max_level) {
   272         return new Elevator(digraph, max_level);
   273       }
   274     };
   275 
   276     /// \brief \ref named-templ-param "Named parameter" for setting
   277     /// Elevator type with automatic allocation
   278     ///
   279     /// \ref named-templ-param "Named parameter" for setting Elevator
   280     /// type with automatic allocation.
   281     /// The Elevator should have standard constructor interface to be
   282     /// able to automatically created by the algorithm (i.e. the
   283     /// digraph and the maximum level should be passed to it).
   284     /// However an external elevator object could also be passed to the
   285     /// algorithm with the \ref elevator(Elevator&) "elevator()" function
   286     /// before calling \ref run() or \ref init().
   287     /// \sa SetElevator
   288     template <typename _Elevator>
   289     struct SetStandardElevator
   290       : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
   291                        SetStandardElevatorTraits<_Elevator> > {
   292       typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
   293                       SetStandardElevatorTraits<_Elevator> > Create;
   294     };
   295 
   296     /// @}
   297 
   298   protected:
   299 
   300     Circulation() {}
   301 
   302   public:
   303 
   304     /// The constructor of the class.
   305 
   306     /// The constructor of the class.
   307     /// \param g The digraph the algorithm runs on.
   308     /// \param lo The lower bound capacity of the arcs.
   309     /// \param up The upper bound capacity of the arcs.
   310     /// \param delta The lower bound for the supply of the nodes.
   311     Circulation(const Digraph &g,const LCapMap &lo,
   312                 const UCapMap &up,const DeltaMap &delta)
   313       : _g(g), _node_num(),
   314         _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
   315         _level(0), _local_level(false), _excess(0), _el() {}
   316 
   317     /// Destructor.
   318     ~Circulation() {
   319       destroyStructures();
   320     }
   321 
   322 
   323   private:
   324 
   325     void createStructures() {
   326       _node_num = _el = countNodes(_g);
   327 
   328       if (!_flow) {
   329         _flow = Traits::createFlowMap(_g);
   330         _local_flow = true;
   331       }
   332       if (!_level) {
   333         _level = Traits::createElevator(_g, _node_num);
   334         _local_level = true;
   335       }
   336       if (!_excess) {
   337         _excess = new ExcessMap(_g);
   338       }
   339     }
   340 
   341     void destroyStructures() {
   342       if (_local_flow) {
   343         delete _flow;
   344       }
   345       if (_local_level) {
   346         delete _level;
   347       }
   348       if (_excess) {
   349         delete _excess;
   350       }
   351     }
   352 
   353   public:
   354 
   355     /// Sets the lower bound capacity map.
   356 
   357     /// Sets the lower bound capacity map.
   358     /// \return <tt>(*this)</tt>
   359     Circulation& lowerCapMap(const LCapMap& map) {
   360       _lo = &map;
   361       return *this;
   362     }
   363 
   364     /// Sets the upper bound capacity map.
   365 
   366     /// Sets the upper bound capacity map.
   367     /// \return <tt>(*this)</tt>
   368     Circulation& upperCapMap(const LCapMap& map) {
   369       _up = &map;
   370       return *this;
   371     }
   372 
   373     /// Sets the lower bound map for the supply of the nodes.
   374 
   375     /// Sets the lower bound map for the supply of the nodes.
   376     /// \return <tt>(*this)</tt>
   377     Circulation& deltaMap(const DeltaMap& map) {
   378       _delta = &map;
   379       return *this;
   380     }
   381 
   382     /// \brief Sets the flow map.
   383     ///
   384     /// Sets the flow map.
   385     /// If you don't use this function before calling \ref run() or
   386     /// \ref init(), an instance will be allocated automatically.
   387     /// The destructor deallocates this automatically allocated map,
   388     /// of course.
   389     /// \return <tt>(*this)</tt>
   390     Circulation& flowMap(FlowMap& map) {
   391       if (_local_flow) {
   392         delete _flow;
   393         _local_flow = false;
   394       }
   395       _flow = &map;
   396       return *this;
   397     }
   398 
   399     /// \brief Sets the elevator used by algorithm.
   400     ///
   401     /// Sets the elevator used by algorithm.
   402     /// If you don't use this function before calling \ref run() or
   403     /// \ref init(), an instance will be allocated automatically.
   404     /// The destructor deallocates this automatically allocated elevator,
   405     /// of course.
   406     /// \return <tt>(*this)</tt>
   407     Circulation& elevator(Elevator& elevator) {
   408       if (_local_level) {
   409         delete _level;
   410         _local_level = false;
   411       }
   412       _level = &elevator;
   413       return *this;
   414     }
   415 
   416     /// \brief Returns a const reference to the elevator.
   417     ///
   418     /// Returns a const reference to the elevator.
   419     ///
   420     /// \pre Either \ref run() or \ref init() must be called before
   421     /// using this function.
   422     const Elevator& elevator() const {
   423       return *_level;
   424     }
   425 
   426     /// \brief Sets the tolerance used by algorithm.
   427     ///
   428     /// Sets the tolerance used by algorithm.
   429     Circulation& tolerance(const Tolerance& tolerance) const {
   430       _tol = tolerance;
   431       return *this;
   432     }
   433 
   434     /// \brief Returns a const reference to the tolerance.
   435     ///
   436     /// Returns a const reference to the tolerance.
   437     const Tolerance& tolerance() const {
   438       return tolerance;
   439     }
   440 
   441     /// \name Execution Control
   442     /// The simplest way to execute the algorithm is to call \ref run().\n
   443     /// If you need more control on the initial solution or the execution,
   444     /// first you have to call one of the \ref init() functions, then
   445     /// the \ref start() function.
   446 
   447     ///@{
   448 
   449     /// Initializes the internal data structures.
   450 
   451     /// Initializes the internal data structures and sets all flow values
   452     /// to the lower bound.
   453     void init()
   454     {
   455       createStructures();
   456 
   457       for(NodeIt n(_g);n!=INVALID;++n) {
   458         _excess->set(n, (*_delta)[n]);
   459       }
   460 
   461       for (ArcIt e(_g);e!=INVALID;++e) {
   462         _flow->set(e, (*_lo)[e]);
   463         _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
   464         _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
   465       }
   466 
   467       // global relabeling tested, but in general case it provides
   468       // worse performance for random digraphs
   469       _level->initStart();
   470       for(NodeIt n(_g);n!=INVALID;++n)
   471         _level->initAddItem(n);
   472       _level->initFinish();
   473       for(NodeIt n(_g);n!=INVALID;++n)
   474         if(_tol.positive((*_excess)[n]))
   475           _level->activate(n);
   476     }
   477 
   478     /// Initializes the internal data structures using a greedy approach.
   479 
   480     /// Initializes the internal data structures using a greedy approach
   481     /// to construct the initial solution.
   482     void greedyInit()
   483     {
   484       createStructures();
   485 
   486       for(NodeIt n(_g);n!=INVALID;++n) {
   487         _excess->set(n, (*_delta)[n]);
   488       }
   489 
   490       for (ArcIt e(_g);e!=INVALID;++e) {
   491         if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
   492           _flow->set(e, (*_up)[e]);
   493           _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
   494           _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
   495         } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
   496           _flow->set(e, (*_lo)[e]);
   497           _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
   498           _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
   499         } else {
   500           Value fc = -(*_excess)[_g.target(e)];
   501           _flow->set(e, fc);
   502           _excess->set(_g.target(e), 0);
   503           _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
   504         }
   505       }
   506 
   507       _level->initStart();
   508       for(NodeIt n(_g);n!=INVALID;++n)
   509         _level->initAddItem(n);
   510       _level->initFinish();
   511       for(NodeIt n(_g);n!=INVALID;++n)
   512         if(_tol.positive((*_excess)[n]))
   513           _level->activate(n);
   514     }
   515 
   516     ///Executes the algorithm
   517 
   518     ///This function executes the algorithm.
   519     ///
   520     ///\return \c true if a feasible circulation is found.
   521     ///
   522     ///\sa barrier()
   523     ///\sa barrierMap()
   524     bool start()
   525     {
   526 
   527       Node act;
   528       Node bact=INVALID;
   529       Node last_activated=INVALID;
   530       while((act=_level->highestActive())!=INVALID) {
   531         int actlevel=(*_level)[act];
   532         int mlevel=_node_num;
   533         Value exc=(*_excess)[act];
   534 
   535         for(OutArcIt e(_g,act);e!=INVALID; ++e) {
   536           Node v = _g.target(e);
   537           Value fc=(*_up)[e]-(*_flow)[e];
   538           if(!_tol.positive(fc)) continue;
   539           if((*_level)[v]<actlevel) {
   540             if(!_tol.less(fc, exc)) {
   541               _flow->set(e, (*_flow)[e] + exc);
   542               _excess->set(v, (*_excess)[v] + exc);
   543               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   544                 _level->activate(v);
   545               _excess->set(act,0);
   546               _level->deactivate(act);
   547               goto next_l;
   548             }
   549             else {
   550               _flow->set(e, (*_up)[e]);
   551               _excess->set(v, (*_excess)[v] + fc);
   552               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   553                 _level->activate(v);
   554               exc-=fc;
   555             }
   556           }
   557           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
   558         }
   559         for(InArcIt e(_g,act);e!=INVALID; ++e) {
   560           Node v = _g.source(e);
   561           Value fc=(*_flow)[e]-(*_lo)[e];
   562           if(!_tol.positive(fc)) continue;
   563           if((*_level)[v]<actlevel) {
   564             if(!_tol.less(fc, exc)) {
   565               _flow->set(e, (*_flow)[e] - exc);
   566               _excess->set(v, (*_excess)[v] + exc);
   567               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   568                 _level->activate(v);
   569               _excess->set(act,0);
   570               _level->deactivate(act);
   571               goto next_l;
   572             }
   573             else {
   574               _flow->set(e, (*_lo)[e]);
   575               _excess->set(v, (*_excess)[v] + fc);
   576               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   577                 _level->activate(v);
   578               exc-=fc;
   579             }
   580           }
   581           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
   582         }
   583 
   584         _excess->set(act, exc);
   585         if(!_tol.positive(exc)) _level->deactivate(act);
   586         else if(mlevel==_node_num) {
   587           _level->liftHighestActiveToTop();
   588           _el = _node_num;
   589           return false;
   590         }
   591         else {
   592           _level->liftHighestActive(mlevel+1);
   593           if(_level->onLevel(actlevel)==0) {
   594             _el = actlevel;
   595             return false;
   596           }
   597         }
   598       next_l:
   599         ;
   600       }
   601       return true;
   602     }
   603 
   604     /// Runs the algorithm.
   605 
   606     /// This function runs the algorithm.
   607     ///
   608     /// \return \c true if a feasible circulation is found.
   609     ///
   610     /// \note Apart from the return value, c.run() is just a shortcut of
   611     /// the following code.
   612     /// \code
   613     ///   c.greedyInit();
   614     ///   c.start();
   615     /// \endcode
   616     bool run() {
   617       greedyInit();
   618       return start();
   619     }
   620 
   621     /// @}
   622 
   623     /// \name Query Functions
   624     /// The results of the circulation algorithm can be obtained using
   625     /// these functions.\n
   626     /// Either \ref run() or \ref start() should be called before
   627     /// using them.
   628 
   629     ///@{
   630 
   631     /// \brief Returns the flow on the given arc.
   632     ///
   633     /// Returns the flow on the given arc.
   634     ///
   635     /// \pre Either \ref run() or \ref init() must be called before
   636     /// using this function.
   637     Value flow(const Arc& arc) const {
   638       return (*_flow)[arc];
   639     }
   640 
   641     /// \brief Returns a const reference to the flow map.
   642     ///
   643     /// Returns a const reference to the arc map storing the found flow.
   644     ///
   645     /// \pre Either \ref run() or \ref init() must be called before
   646     /// using this function.
   647     const FlowMap& flowMap() const {
   648       return *_flow;
   649     }
   650 
   651     /**
   652        \brief Returns \c true if the given node is in a barrier.
   653 
   654        Barrier is a set \e B of nodes for which
   655 
   656        \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
   657            \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
   658 
   659        holds. The existence of a set with this property prooves that a
   660        feasible circualtion cannot exist.
   661 
   662        This function returns \c true if the given node is in the found
   663        barrier. If a feasible circulation is found, the function
   664        gives back \c false for every node.
   665 
   666        \pre Either \ref run() or \ref init() must be called before
   667        using this function.
   668 
   669        \sa barrierMap()
   670        \sa checkBarrier()
   671     */
   672     bool barrier(const Node& node) const
   673     {
   674       return (*_level)[node] >= _el;
   675     }
   676 
   677     /// \brief Gives back a barrier.
   678     ///
   679     /// This function sets \c bar to the characteristic vector of the
   680     /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
   681     /// node map with \c bool (or convertible) value type.
   682     ///
   683     /// If a feasible circulation is found, the function gives back an
   684     /// empty set, so \c bar[v] will be \c false for all nodes \c v.
   685     ///
   686     /// \note This function calls \ref barrier() for each node,
   687     /// so it runs in \f$O(n)\f$ time.
   688     ///
   689     /// \pre Either \ref run() or \ref init() must be called before
   690     /// using this function.
   691     ///
   692     /// \sa barrier()
   693     /// \sa checkBarrier()
   694     template<class BarrierMap>
   695     void barrierMap(BarrierMap &bar) const
   696     {
   697       for(NodeIt n(_g);n!=INVALID;++n)
   698         bar.set(n, (*_level)[n] >= _el);
   699     }
   700 
   701     /// @}
   702 
   703     /// \name Checker Functions
   704     /// The feasibility of the results can be checked using
   705     /// these functions.\n
   706     /// Either \ref run() or \ref start() should be called before
   707     /// using them.
   708 
   709     ///@{
   710 
   711     ///Check if the found flow is a feasible circulation
   712 
   713     ///Check if the found flow is a feasible circulation,
   714     ///
   715     bool checkFlow() const {
   716       for(ArcIt e(_g);e!=INVALID;++e)
   717         if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
   718       for(NodeIt n(_g);n!=INVALID;++n)
   719         {
   720           Value dif=-(*_delta)[n];
   721           for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
   722           for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
   723           if(_tol.negative(dif)) return false;
   724         }
   725       return true;
   726     }
   727 
   728     ///Check whether or not the last execution provides a barrier
   729 
   730     ///Check whether or not the last execution provides a barrier.
   731     ///\sa barrier()
   732     ///\sa barrierMap()
   733     bool checkBarrier() const
   734     {
   735       Value delta=0;
   736       for(NodeIt n(_g);n!=INVALID;++n)
   737         if(barrier(n))
   738           delta-=(*_delta)[n];
   739       for(ArcIt e(_g);e!=INVALID;++e)
   740         {
   741           Node s=_g.source(e);
   742           Node t=_g.target(e);
   743           if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
   744           else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
   745         }
   746       return _tol.negative(delta);
   747     }
   748 
   749     /// @}
   750 
   751   };
   752 
   753 }
   754 
   755 #endif