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