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