lemon/circulation.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 24 Apr 2009 12:22:06 +0200
changeset 613 b1811c363299
parent 581 aa1804409f29
parent 610 dacc2cee2b4c
child 622 28f58740b6f8
permissions -rw-r--r--
Bug fix in NetworkSimplex (#234)
     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   ///
    35   /// \tparam GR Type of the digraph the algorithm runs on.
    36   /// \tparam LM The type of the lower bound map.
    37   /// \tparam UM The type of the upper bound (capacity) map.
    38   /// \tparam SM The type of the supply map.
    39   template <typename GR, typename LM,
    40             typename UM, typename SM>
    41   struct CirculationDefaultTraits {
    42 
    43     /// \brief The type of the digraph the algorithm runs on.
    44     typedef GR Digraph;
    45 
    46     /// \brief The type of the lower bound map.
    47     ///
    48     /// The type of the map that stores the lower bounds on the arcs.
    49     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
    50     typedef LM LowerMap;
    51 
    52     /// \brief The type of the upper bound (capacity) map.
    53     ///
    54     /// The type of the map that stores the upper bounds (capacities)
    55     /// on the arcs.
    56     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
    57     typedef UM UpperMap;
    58 
    59     /// \brief The type of supply map.
    60     ///
    61     /// The type of the map that stores the signed supply values of the 
    62     /// nodes. 
    63     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
    64     typedef SM SupplyMap;
    65 
    66     /// \brief The type of the flow values.
    67     typedef typename SupplyMap::Value Flow;
    68 
    69     /// \brief The type of the map that stores the flow values.
    70     ///
    71     /// The type of the map that stores the flow values.
    72     /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
    73     /// concept.
    74     typedef typename Digraph::template ArcMap<Flow> FlowMap;
    75 
    76     /// \brief Instantiates a FlowMap.
    77     ///
    78     /// This function instantiates a \ref FlowMap.
    79     /// \param digraph The digraph for 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 for 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<Flow> 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 \e network
   115      \e 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 are
   118      given for the difference between the outgoing and incoming flow
   119      at the nodes.
   120 
   121      The exact formulation of this problem is the following.
   122      Let \f$G=(V,A)\f$ be a digraph,
   123      \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$ denote the lower and
   124      upper bounds on the arcs, for which \f$0 \leq lower(uv) \leq upper(uv)\f$
   125      holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
   126      denotes the signed supply values of the nodes.
   127      If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
   128      supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
   129      \f$-sup(u)\f$ demand.
   130      A feasible circulation is an \f$f: A\rightarrow\mathbf{R}^+_0\f$
   131      solution of the following problem.
   132 
   133      \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
   134      \geq sup(u) \quad \forall u\in V, \f]
   135      \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
   136      
   137      The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
   138      zero or negative in order to have a feasible solution (since the sum
   139      of the expressions on the left-hand side of the inequalities is zero).
   140      It means that the total demand must be greater or equal to the total
   141      supply and all the supplies have to be carried out from the supply nodes,
   142      but there could be demands that are not satisfied.
   143      If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
   144      constraints have to be satisfied with equality, i.e. all demands
   145      have to be satisfied and all supplies have to be used.
   146      
   147      If you need the opposite inequalities in the supply/demand constraints
   148      (i.e. the total demand is less than the total supply and all the demands
   149      have to be satisfied while there could be supplies that are not used),
   150      then you could easily transform the problem to the above form by reversing
   151      the direction of the arcs and taking the negative of the supply values
   152      (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
   153 
   154      Note that this algorithm also provides a feasible solution for the
   155      \ref min_cost_flow "minimum cost flow problem".
   156 
   157      \tparam GR The type of the digraph the algorithm runs on.
   158      \tparam LM The type of the lower bound map. The default
   159      map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
   160      \tparam UM The type of the upper bound (capacity) map.
   161      The default map type is \c LM.
   162      \tparam SM The type of the supply map. The default map type is
   163      \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
   164   */
   165 #ifdef DOXYGEN
   166 template< typename GR,
   167           typename LM,
   168           typename UM,
   169           typename SM,
   170           typename TR >
   171 #else
   172 template< typename GR,
   173           typename LM = typename GR::template ArcMap<int>,
   174           typename UM = LM,
   175           typename SM = typename GR::template NodeMap<typename UM::Value>,
   176           typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
   177 #endif
   178   class Circulation {
   179   public:
   180 
   181     ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
   182     typedef TR Traits;
   183     ///The type of the digraph the algorithm runs on.
   184     typedef typename Traits::Digraph Digraph;
   185     ///The type of the flow values.
   186     typedef typename Traits::Flow Flow;
   187 
   188     ///The type of the lower bound map.
   189     typedef typename Traits::LowerMap LowerMap;
   190     ///The type of the upper bound (capacity) map.
   191     typedef typename Traits::UpperMap UpperMap;
   192     ///The type of the supply map.
   193     typedef typename Traits::SupplyMap SupplyMap;
   194     ///The type of the flow map.
   195     typedef typename Traits::FlowMap FlowMap;
   196 
   197     ///The type of the elevator.
   198     typedef typename Traits::Elevator Elevator;
   199     ///The type of the tolerance.
   200     typedef typename Traits::Tolerance Tolerance;
   201 
   202   private:
   203 
   204     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   205 
   206     const Digraph &_g;
   207     int _node_num;
   208 
   209     const LowerMap *_lo;
   210     const UpperMap *_up;
   211     const SupplyMap *_supply;
   212 
   213     FlowMap *_flow;
   214     bool _local_flow;
   215 
   216     Elevator* _level;
   217     bool _local_level;
   218 
   219     typedef typename Digraph::template NodeMap<Flow> ExcessMap;
   220     ExcessMap* _excess;
   221 
   222     Tolerance _tol;
   223     int _el;
   224 
   225   public:
   226 
   227     typedef Circulation Create;
   228 
   229     ///\name Named Template Parameters
   230 
   231     ///@{
   232 
   233     template <typename T>
   234     struct SetFlowMapTraits : public Traits {
   235       typedef T FlowMap;
   236       static FlowMap *createFlowMap(const Digraph&) {
   237         LEMON_ASSERT(false, "FlowMap is not initialized");
   238         return 0; // ignore warnings
   239       }
   240     };
   241 
   242     /// \brief \ref named-templ-param "Named parameter" for setting
   243     /// FlowMap type
   244     ///
   245     /// \ref named-templ-param "Named parameter" for setting FlowMap
   246     /// type.
   247     template <typename T>
   248     struct SetFlowMap
   249       : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
   250                            SetFlowMapTraits<T> > {
   251       typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
   252                           SetFlowMapTraits<T> > Create;
   253     };
   254 
   255     template <typename T>
   256     struct SetElevatorTraits : public Traits {
   257       typedef T Elevator;
   258       static Elevator *createElevator(const Digraph&, int) {
   259         LEMON_ASSERT(false, "Elevator is not initialized");
   260         return 0; // ignore warnings
   261       }
   262     };
   263 
   264     /// \brief \ref named-templ-param "Named parameter" for setting
   265     /// Elevator type
   266     ///
   267     /// \ref named-templ-param "Named parameter" for setting Elevator
   268     /// type. If this named parameter is used, then an external
   269     /// elevator object must be passed to the algorithm using the
   270     /// \ref elevator(Elevator&) "elevator()" function before calling
   271     /// \ref run() or \ref init().
   272     /// \sa SetStandardElevator
   273     template <typename T>
   274     struct SetElevator
   275       : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
   276                            SetElevatorTraits<T> > {
   277       typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
   278                           SetElevatorTraits<T> > Create;
   279     };
   280 
   281     template <typename T>
   282     struct SetStandardElevatorTraits : public Traits {
   283       typedef T Elevator;
   284       static Elevator *createElevator(const Digraph& digraph, int max_level) {
   285         return new Elevator(digraph, max_level);
   286       }
   287     };
   288 
   289     /// \brief \ref named-templ-param "Named parameter" for setting
   290     /// Elevator type with automatic allocation
   291     ///
   292     /// \ref named-templ-param "Named parameter" for setting Elevator
   293     /// type with automatic allocation.
   294     /// The Elevator should have standard constructor interface to be
   295     /// able to automatically created by the algorithm (i.e. the
   296     /// digraph and the maximum level should be passed to it).
   297     /// However an external elevator object could also be passed to the
   298     /// algorithm with the \ref elevator(Elevator&) "elevator()" function
   299     /// before calling \ref run() or \ref init().
   300     /// \sa SetElevator
   301     template <typename T>
   302     struct SetStandardElevator
   303       : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
   304                        SetStandardElevatorTraits<T> > {
   305       typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
   306                       SetStandardElevatorTraits<T> > Create;
   307     };
   308 
   309     /// @}
   310 
   311   protected:
   312 
   313     Circulation() {}
   314 
   315   public:
   316 
   317     /// Constructor.
   318 
   319     /// The constructor of the class.
   320     ///
   321     /// \param graph The digraph the algorithm runs on.
   322     /// \param lower The lower bounds for the flow values on the arcs.
   323     /// \param upper The upper bounds (capacities) for the flow values 
   324     /// on the arcs.
   325     /// \param supply The signed supply values of the nodes.
   326     Circulation(const Digraph &graph, const LowerMap &lower,
   327                 const UpperMap &upper, const SupplyMap &supply)
   328       : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
   329         _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
   330         _excess(NULL) {}
   331 
   332     /// Destructor.
   333     ~Circulation() {
   334       destroyStructures();
   335     }
   336 
   337 
   338   private:
   339 
   340     void createStructures() {
   341       _node_num = _el = countNodes(_g);
   342 
   343       if (!_flow) {
   344         _flow = Traits::createFlowMap(_g);
   345         _local_flow = true;
   346       }
   347       if (!_level) {
   348         _level = Traits::createElevator(_g, _node_num);
   349         _local_level = true;
   350       }
   351       if (!_excess) {
   352         _excess = new ExcessMap(_g);
   353       }
   354     }
   355 
   356     void destroyStructures() {
   357       if (_local_flow) {
   358         delete _flow;
   359       }
   360       if (_local_level) {
   361         delete _level;
   362       }
   363       if (_excess) {
   364         delete _excess;
   365       }
   366     }
   367 
   368   public:
   369 
   370     /// Sets the lower bound map.
   371 
   372     /// Sets the lower bound map.
   373     /// \return <tt>(*this)</tt>
   374     Circulation& lowerMap(const LowerMap& map) {
   375       _lo = &map;
   376       return *this;
   377     }
   378 
   379     /// Sets the upper bound (capacity) map.
   380 
   381     /// Sets the upper bound (capacity) map.
   382     /// \return <tt>(*this)</tt>
   383     Circulation& upperMap(const LowerMap& map) {
   384       _up = &map;
   385       return *this;
   386     }
   387 
   388     /// Sets the supply map.
   389 
   390     /// Sets the supply map.
   391     /// \return <tt>(*this)</tt>
   392     Circulation& supplyMap(const SupplyMap& map) {
   393       _supply = &map;
   394       return *this;
   395     }
   396 
   397     /// \brief Sets the flow map.
   398     ///
   399     /// Sets the flow map.
   400     /// If you don't use this function before calling \ref run() or
   401     /// \ref init(), an instance will be allocated automatically.
   402     /// The destructor deallocates this automatically allocated map,
   403     /// of course.
   404     /// \return <tt>(*this)</tt>
   405     Circulation& flowMap(FlowMap& map) {
   406       if (_local_flow) {
   407         delete _flow;
   408         _local_flow = false;
   409       }
   410       _flow = &map;
   411       return *this;
   412     }
   413 
   414     /// \brief Sets the elevator used by algorithm.
   415     ///
   416     /// Sets the elevator used by algorithm.
   417     /// If you don't use this function before calling \ref run() or
   418     /// \ref init(), an instance will be allocated automatically.
   419     /// The destructor deallocates this automatically allocated elevator,
   420     /// of course.
   421     /// \return <tt>(*this)</tt>
   422     Circulation& elevator(Elevator& elevator) {
   423       if (_local_level) {
   424         delete _level;
   425         _local_level = false;
   426       }
   427       _level = &elevator;
   428       return *this;
   429     }
   430 
   431     /// \brief Returns a const reference to the elevator.
   432     ///
   433     /// Returns a const reference to the elevator.
   434     ///
   435     /// \pre Either \ref run() or \ref init() must be called before
   436     /// using this function.
   437     const Elevator& elevator() const {
   438       return *_level;
   439     }
   440 
   441     /// \brief Sets the tolerance used by algorithm.
   442     ///
   443     /// Sets the tolerance used by algorithm.
   444     Circulation& tolerance(const Tolerance& tolerance) const {
   445       _tol = tolerance;
   446       return *this;
   447     }
   448 
   449     /// \brief Returns a const reference to the tolerance.
   450     ///
   451     /// Returns a const reference to the tolerance.
   452     const Tolerance& tolerance() const {
   453       return tolerance;
   454     }
   455 
   456     /// \name Execution Control
   457     /// The simplest way to execute the algorithm is to call \ref run().\n
   458     /// If you need more control on the initial solution or the execution,
   459     /// first you have to call one of the \ref init() functions, then
   460     /// the \ref start() function.
   461 
   462     ///@{
   463 
   464     /// Initializes the internal data structures.
   465 
   466     /// Initializes the internal data structures and sets all flow values
   467     /// to the lower bound.
   468     void init()
   469     {
   470       createStructures();
   471 
   472       for(NodeIt n(_g);n!=INVALID;++n) {
   473         (*_excess)[n] = (*_supply)[n];
   474       }
   475 
   476       for (ArcIt e(_g);e!=INVALID;++e) {
   477         _flow->set(e, (*_lo)[e]);
   478         (*_excess)[_g.target(e)] += (*_flow)[e];
   479         (*_excess)[_g.source(e)] -= (*_flow)[e];
   480       }
   481 
   482       // global relabeling tested, but in general case it provides
   483       // worse performance for random digraphs
   484       _level->initStart();
   485       for(NodeIt n(_g);n!=INVALID;++n)
   486         _level->initAddItem(n);
   487       _level->initFinish();
   488       for(NodeIt n(_g);n!=INVALID;++n)
   489         if(_tol.positive((*_excess)[n]))
   490           _level->activate(n);
   491     }
   492 
   493     /// Initializes the internal data structures using a greedy approach.
   494 
   495     /// Initializes the internal data structures using a greedy approach
   496     /// to construct the initial solution.
   497     void greedyInit()
   498     {
   499       createStructures();
   500 
   501       for(NodeIt n(_g);n!=INVALID;++n) {
   502         (*_excess)[n] = (*_supply)[n];
   503       }
   504 
   505       for (ArcIt e(_g);e!=INVALID;++e) {
   506         if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
   507           _flow->set(e, (*_up)[e]);
   508           (*_excess)[_g.target(e)] += (*_up)[e];
   509           (*_excess)[_g.source(e)] -= (*_up)[e];
   510         } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
   511           _flow->set(e, (*_lo)[e]);
   512           (*_excess)[_g.target(e)] += (*_lo)[e];
   513           (*_excess)[_g.source(e)] -= (*_lo)[e];
   514         } else {
   515           Flow fc = -(*_excess)[_g.target(e)];
   516           _flow->set(e, fc);
   517           (*_excess)[_g.target(e)] = 0;
   518           (*_excess)[_g.source(e)] -= fc;
   519         }
   520       }
   521 
   522       _level->initStart();
   523       for(NodeIt n(_g);n!=INVALID;++n)
   524         _level->initAddItem(n);
   525       _level->initFinish();
   526       for(NodeIt n(_g);n!=INVALID;++n)
   527         if(_tol.positive((*_excess)[n]))
   528           _level->activate(n);
   529     }
   530 
   531     ///Executes the algorithm
   532 
   533     ///This function executes the algorithm.
   534     ///
   535     ///\return \c true if a feasible circulation is found.
   536     ///
   537     ///\sa barrier()
   538     ///\sa barrierMap()
   539     bool start()
   540     {
   541 
   542       Node act;
   543       Node bact=INVALID;
   544       Node last_activated=INVALID;
   545       while((act=_level->highestActive())!=INVALID) {
   546         int actlevel=(*_level)[act];
   547         int mlevel=_node_num;
   548         Flow exc=(*_excess)[act];
   549 
   550         for(OutArcIt e(_g,act);e!=INVALID; ++e) {
   551           Node v = _g.target(e);
   552           Flow fc=(*_up)[e]-(*_flow)[e];
   553           if(!_tol.positive(fc)) continue;
   554           if((*_level)[v]<actlevel) {
   555             if(!_tol.less(fc, exc)) {
   556               _flow->set(e, (*_flow)[e] + exc);
   557               (*_excess)[v] += exc;
   558               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   559                 _level->activate(v);
   560               (*_excess)[act] = 0;
   561               _level->deactivate(act);
   562               goto next_l;
   563             }
   564             else {
   565               _flow->set(e, (*_up)[e]);
   566               (*_excess)[v] += fc;
   567               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   568                 _level->activate(v);
   569               exc-=fc;
   570             }
   571           }
   572           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
   573         }
   574         for(InArcIt e(_g,act);e!=INVALID; ++e) {
   575           Node v = _g.source(e);
   576           Flow fc=(*_flow)[e]-(*_lo)[e];
   577           if(!_tol.positive(fc)) continue;
   578           if((*_level)[v]<actlevel) {
   579             if(!_tol.less(fc, exc)) {
   580               _flow->set(e, (*_flow)[e] - exc);
   581               (*_excess)[v] += exc;
   582               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   583                 _level->activate(v);
   584               (*_excess)[act] = 0;
   585               _level->deactivate(act);
   586               goto next_l;
   587             }
   588             else {
   589               _flow->set(e, (*_lo)[e]);
   590               (*_excess)[v] += fc;
   591               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   592                 _level->activate(v);
   593               exc-=fc;
   594             }
   595           }
   596           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
   597         }
   598 
   599         (*_excess)[act] = exc;
   600         if(!_tol.positive(exc)) _level->deactivate(act);
   601         else if(mlevel==_node_num) {
   602           _level->liftHighestActiveToTop();
   603           _el = _node_num;
   604           return false;
   605         }
   606         else {
   607           _level->liftHighestActive(mlevel+1);
   608           if(_level->onLevel(actlevel)==0) {
   609             _el = actlevel;
   610             return false;
   611           }
   612         }
   613       next_l:
   614         ;
   615       }
   616       return true;
   617     }
   618 
   619     /// Runs the algorithm.
   620 
   621     /// This function runs the algorithm.
   622     ///
   623     /// \return \c true if a feasible circulation is found.
   624     ///
   625     /// \note Apart from the return value, c.run() is just a shortcut of
   626     /// the following code.
   627     /// \code
   628     ///   c.greedyInit();
   629     ///   c.start();
   630     /// \endcode
   631     bool run() {
   632       greedyInit();
   633       return start();
   634     }
   635 
   636     /// @}
   637 
   638     /// \name Query Functions
   639     /// The results of the circulation algorithm can be obtained using
   640     /// these functions.\n
   641     /// Either \ref run() or \ref start() should be called before
   642     /// using them.
   643 
   644     ///@{
   645 
   646     /// \brief Returns the flow on the given arc.
   647     ///
   648     /// Returns the flow on the given arc.
   649     ///
   650     /// \pre Either \ref run() or \ref init() must be called before
   651     /// using this function.
   652     Flow flow(const Arc& arc) const {
   653       return (*_flow)[arc];
   654     }
   655 
   656     /// \brief Returns a const reference to the flow map.
   657     ///
   658     /// Returns a const reference to the arc map storing the found flow.
   659     ///
   660     /// \pre Either \ref run() or \ref init() must be called before
   661     /// using this function.
   662     const FlowMap& flowMap() const {
   663       return *_flow;
   664     }
   665 
   666     /**
   667        \brief Returns \c true if the given node is in a barrier.
   668 
   669        Barrier is a set \e B of nodes for which
   670 
   671        \f[ \sum_{uv\in A: u\in B} upper(uv) -
   672            \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
   673 
   674        holds. The existence of a set with this property prooves that a
   675        feasible circualtion cannot exist.
   676 
   677        This function returns \c true if the given node is in the found
   678        barrier. If a feasible circulation is found, the function
   679        gives back \c false for every node.
   680 
   681        \pre Either \ref run() or \ref init() must be called before
   682        using this function.
   683 
   684        \sa barrierMap()
   685        \sa checkBarrier()
   686     */
   687     bool barrier(const Node& node) const
   688     {
   689       return (*_level)[node] >= _el;
   690     }
   691 
   692     /// \brief Gives back a barrier.
   693     ///
   694     /// This function sets \c bar to the characteristic vector of the
   695     /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
   696     /// node map with \c bool (or convertible) value type.
   697     ///
   698     /// If a feasible circulation is found, the function gives back an
   699     /// empty set, so \c bar[v] will be \c false for all nodes \c v.
   700     ///
   701     /// \note This function calls \ref barrier() for each node,
   702     /// so it runs in O(n) time.
   703     ///
   704     /// \pre Either \ref run() or \ref init() must be called before
   705     /// using this function.
   706     ///
   707     /// \sa barrier()
   708     /// \sa checkBarrier()
   709     template<class BarrierMap>
   710     void barrierMap(BarrierMap &bar) const
   711     {
   712       for(NodeIt n(_g);n!=INVALID;++n)
   713         bar.set(n, (*_level)[n] >= _el);
   714     }
   715 
   716     /// @}
   717 
   718     /// \name Checker Functions
   719     /// The feasibility of the results can be checked using
   720     /// these functions.\n
   721     /// Either \ref run() or \ref start() should be called before
   722     /// using them.
   723 
   724     ///@{
   725 
   726     ///Check if the found flow is a feasible circulation
   727 
   728     ///Check if the found flow is a feasible circulation,
   729     ///
   730     bool checkFlow() const {
   731       for(ArcIt e(_g);e!=INVALID;++e)
   732         if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
   733       for(NodeIt n(_g);n!=INVALID;++n)
   734         {
   735           Flow dif=-(*_supply)[n];
   736           for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
   737           for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
   738           if(_tol.negative(dif)) return false;
   739         }
   740       return true;
   741     }
   742 
   743     ///Check whether or not the last execution provides a barrier
   744 
   745     ///Check whether or not the last execution provides a barrier.
   746     ///\sa barrier()
   747     ///\sa barrierMap()
   748     bool checkBarrier() const
   749     {
   750       Flow delta=0;
   751       for(NodeIt n(_g);n!=INVALID;++n)
   752         if(barrier(n))
   753           delta-=(*_supply)[n];
   754       for(ArcIt e(_g);e!=INVALID;++e)
   755         {
   756           Node s=_g.source(e);
   757           Node t=_g.target(e);
   758           if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
   759           else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
   760         }
   761       return _tol.negative(delta);
   762     }
   763 
   764     /// @}
   765 
   766   };
   767 
   768 }
   769 
   770 #endif