lemon/circulation.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sat, 20 Feb 2010 18:39:03 +0100
changeset 839 f3bc4e9b5f3a
parent 715 ece80147fb08
child 825 75e6020b19b1
permissions -rw-r--r--
New heuristics for MCF algorithms (#340)
and some implementation improvements.

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