lemon/circulation.h
author Balazs Dezso <deba@inf.elte.hu>
Wed, 09 Sep 2009 15:32:03 +0200
changeset 797 30cb42e3e43a
parent 688 1f08e846df29
child 715 ece80147fb08
permissions -rw-r--r--
Port planarity related algorithms from SVN 3509 (#62)
     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 the algorithm.
   454     ///
   455     /// Sets the tolerance object used by the algorithm.
   456     /// \return <tt>(*this)</tt>
   457     Circulation& tolerance(const Tolerance& tolerance) {
   458       _tol = tolerance;
   459       return *this;
   460     }
   461 
   462     /// \brief Returns a const reference to the tolerance.
   463     ///
   464     /// Returns a const reference to the tolerance object used by
   465     /// the algorithm.
   466     const Tolerance& tolerance() const {
   467       return _tol;
   468     }
   469 
   470     /// \name Execution Control
   471     /// The simplest way to execute the algorithm is to call \ref run().\n
   472     /// If you need more control on the initial solution or the execution,
   473     /// first you have to call one of the \ref init() functions, then
   474     /// the \ref start() function.
   475 
   476     ///@{
   477 
   478     /// Initializes the internal data structures.
   479 
   480     /// Initializes the internal data structures and sets all flow values
   481     /// to the lower bound.
   482     void init()
   483     {
   484       LEMON_DEBUG(checkBoundMaps(),
   485         "Upper bounds must be greater or equal to the lower bounds");
   486 
   487       createStructures();
   488 
   489       for(NodeIt n(_g);n!=INVALID;++n) {
   490         (*_excess)[n] = (*_supply)[n];
   491       }
   492 
   493       for (ArcIt e(_g);e!=INVALID;++e) {
   494         _flow->set(e, (*_lo)[e]);
   495         (*_excess)[_g.target(e)] += (*_flow)[e];
   496         (*_excess)[_g.source(e)] -= (*_flow)[e];
   497       }
   498 
   499       // global relabeling tested, but in general case it provides
   500       // worse performance for random digraphs
   501       _level->initStart();
   502       for(NodeIt n(_g);n!=INVALID;++n)
   503         _level->initAddItem(n);
   504       _level->initFinish();
   505       for(NodeIt n(_g);n!=INVALID;++n)
   506         if(_tol.positive((*_excess)[n]))
   507           _level->activate(n);
   508     }
   509 
   510     /// Initializes the internal data structures using a greedy approach.
   511 
   512     /// Initializes the internal data structures using a greedy approach
   513     /// to construct the initial solution.
   514     void greedyInit()
   515     {
   516       LEMON_DEBUG(checkBoundMaps(),
   517         "Upper bounds must be greater or equal to the lower bounds");
   518 
   519       createStructures();
   520 
   521       for(NodeIt n(_g);n!=INVALID;++n) {
   522         (*_excess)[n] = (*_supply)[n];
   523       }
   524 
   525       for (ArcIt e(_g);e!=INVALID;++e) {
   526         if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
   527           _flow->set(e, (*_up)[e]);
   528           (*_excess)[_g.target(e)] += (*_up)[e];
   529           (*_excess)[_g.source(e)] -= (*_up)[e];
   530         } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
   531           _flow->set(e, (*_lo)[e]);
   532           (*_excess)[_g.target(e)] += (*_lo)[e];
   533           (*_excess)[_g.source(e)] -= (*_lo)[e];
   534         } else {
   535           Value fc = -(*_excess)[_g.target(e)];
   536           _flow->set(e, fc);
   537           (*_excess)[_g.target(e)] = 0;
   538           (*_excess)[_g.source(e)] -= fc;
   539         }
   540       }
   541 
   542       _level->initStart();
   543       for(NodeIt n(_g);n!=INVALID;++n)
   544         _level->initAddItem(n);
   545       _level->initFinish();
   546       for(NodeIt n(_g);n!=INVALID;++n)
   547         if(_tol.positive((*_excess)[n]))
   548           _level->activate(n);
   549     }
   550 
   551     ///Executes the algorithm
   552 
   553     ///This function executes the algorithm.
   554     ///
   555     ///\return \c true if a feasible circulation is found.
   556     ///
   557     ///\sa barrier()
   558     ///\sa barrierMap()
   559     bool start()
   560     {
   561 
   562       Node act;
   563       Node bact=INVALID;
   564       Node last_activated=INVALID;
   565       while((act=_level->highestActive())!=INVALID) {
   566         int actlevel=(*_level)[act];
   567         int mlevel=_node_num;
   568         Value exc=(*_excess)[act];
   569 
   570         for(OutArcIt e(_g,act);e!=INVALID; ++e) {
   571           Node v = _g.target(e);
   572           Value fc=(*_up)[e]-(*_flow)[e];
   573           if(!_tol.positive(fc)) continue;
   574           if((*_level)[v]<actlevel) {
   575             if(!_tol.less(fc, exc)) {
   576               _flow->set(e, (*_flow)[e] + exc);
   577               (*_excess)[v] += exc;
   578               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   579                 _level->activate(v);
   580               (*_excess)[act] = 0;
   581               _level->deactivate(act);
   582               goto next_l;
   583             }
   584             else {
   585               _flow->set(e, (*_up)[e]);
   586               (*_excess)[v] += fc;
   587               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   588                 _level->activate(v);
   589               exc-=fc;
   590             }
   591           }
   592           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
   593         }
   594         for(InArcIt e(_g,act);e!=INVALID; ++e) {
   595           Node v = _g.source(e);
   596           Value fc=(*_flow)[e]-(*_lo)[e];
   597           if(!_tol.positive(fc)) continue;
   598           if((*_level)[v]<actlevel) {
   599             if(!_tol.less(fc, exc)) {
   600               _flow->set(e, (*_flow)[e] - exc);
   601               (*_excess)[v] += exc;
   602               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   603                 _level->activate(v);
   604               (*_excess)[act] = 0;
   605               _level->deactivate(act);
   606               goto next_l;
   607             }
   608             else {
   609               _flow->set(e, (*_lo)[e]);
   610               (*_excess)[v] += fc;
   611               if(!_level->active(v) && _tol.positive((*_excess)[v]))
   612                 _level->activate(v);
   613               exc-=fc;
   614             }
   615           }
   616           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
   617         }
   618 
   619         (*_excess)[act] = exc;
   620         if(!_tol.positive(exc)) _level->deactivate(act);
   621         else if(mlevel==_node_num) {
   622           _level->liftHighestActiveToTop();
   623           _el = _node_num;
   624           return false;
   625         }
   626         else {
   627           _level->liftHighestActive(mlevel+1);
   628           if(_level->onLevel(actlevel)==0) {
   629             _el = actlevel;
   630             return false;
   631           }
   632         }
   633       next_l:
   634         ;
   635       }
   636       return true;
   637     }
   638 
   639     /// Runs the algorithm.
   640 
   641     /// This function runs the algorithm.
   642     ///
   643     /// \return \c true if a feasible circulation is found.
   644     ///
   645     /// \note Apart from the return value, c.run() is just a shortcut of
   646     /// the following code.
   647     /// \code
   648     ///   c.greedyInit();
   649     ///   c.start();
   650     /// \endcode
   651     bool run() {
   652       greedyInit();
   653       return start();
   654     }
   655 
   656     /// @}
   657 
   658     /// \name Query Functions
   659     /// The results of the circulation algorithm can be obtained using
   660     /// these functions.\n
   661     /// Either \ref run() or \ref start() should be called before
   662     /// using them.
   663 
   664     ///@{
   665 
   666     /// \brief Returns the flow value on the given arc.
   667     ///
   668     /// Returns the flow value on the given arc.
   669     ///
   670     /// \pre Either \ref run() or \ref init() must be called before
   671     /// using this function.
   672     Value flow(const Arc& arc) const {
   673       return (*_flow)[arc];
   674     }
   675 
   676     /// \brief Returns a const reference to the flow map.
   677     ///
   678     /// Returns a const reference to the arc map storing the found flow.
   679     ///
   680     /// \pre Either \ref run() or \ref init() must be called before
   681     /// using this function.
   682     const FlowMap& flowMap() const {
   683       return *_flow;
   684     }
   685 
   686     /**
   687        \brief Returns \c true if the given node is in a barrier.
   688 
   689        Barrier is a set \e B of nodes for which
   690 
   691        \f[ \sum_{uv\in A: u\in B} upper(uv) -
   692            \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
   693 
   694        holds. The existence of a set with this property prooves that a
   695        feasible circualtion cannot exist.
   696 
   697        This function returns \c true if the given node is in the found
   698        barrier. If a feasible circulation is found, the function
   699        gives back \c false for every node.
   700 
   701        \pre Either \ref run() or \ref init() must be called before
   702        using this function.
   703 
   704        \sa barrierMap()
   705        \sa checkBarrier()
   706     */
   707     bool barrier(const Node& node) const
   708     {
   709       return (*_level)[node] >= _el;
   710     }
   711 
   712     /// \brief Gives back a barrier.
   713     ///
   714     /// This function sets \c bar to the characteristic vector of the
   715     /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
   716     /// node map with \c bool (or convertible) value type.
   717     ///
   718     /// If a feasible circulation is found, the function gives back an
   719     /// empty set, so \c bar[v] will be \c false for all nodes \c v.
   720     ///
   721     /// \note This function calls \ref barrier() for each node,
   722     /// so it runs in O(n) time.
   723     ///
   724     /// \pre Either \ref run() or \ref init() must be called before
   725     /// using this function.
   726     ///
   727     /// \sa barrier()
   728     /// \sa checkBarrier()
   729     template<class BarrierMap>
   730     void barrierMap(BarrierMap &bar) const
   731     {
   732       for(NodeIt n(_g);n!=INVALID;++n)
   733         bar.set(n, (*_level)[n] >= _el);
   734     }
   735 
   736     /// @}
   737 
   738     /// \name Checker Functions
   739     /// The feasibility of the results can be checked using
   740     /// these functions.\n
   741     /// Either \ref run() or \ref start() should be called before
   742     /// using them.
   743 
   744     ///@{
   745 
   746     ///Check if the found flow is a feasible circulation
   747 
   748     ///Check if the found flow is a feasible circulation,
   749     ///
   750     bool checkFlow() const {
   751       for(ArcIt e(_g);e!=INVALID;++e)
   752         if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
   753       for(NodeIt n(_g);n!=INVALID;++n)
   754         {
   755           Value dif=-(*_supply)[n];
   756           for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
   757           for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
   758           if(_tol.negative(dif)) return false;
   759         }
   760       return true;
   761     }
   762 
   763     ///Check whether or not the last execution provides a barrier
   764 
   765     ///Check whether or not the last execution provides a barrier.
   766     ///\sa barrier()
   767     ///\sa barrierMap()
   768     bool checkBarrier() const
   769     {
   770       Value delta=0;
   771       Value inf_cap = std::numeric_limits<Value>::has_infinity ?
   772         std::numeric_limits<Value>::infinity() :
   773         std::numeric_limits<Value>::max();
   774       for(NodeIt n(_g);n!=INVALID;++n)
   775         if(barrier(n))
   776           delta-=(*_supply)[n];
   777       for(ArcIt e(_g);e!=INVALID;++e)
   778         {
   779           Node s=_g.source(e);
   780           Node t=_g.target(e);
   781           if(barrier(s)&&!barrier(t)) {
   782             if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
   783             delta+=(*_up)[e];
   784           }
   785           else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
   786         }
   787       return _tol.negative(delta);
   788     }
   789 
   790     /// @}
   791 
   792   };
   793 
   794 }
   795 
   796 #endif