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