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