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