COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/circulation.h @ 622:28f58740b6f8

Last change on this file since 622:28f58740b6f8 was 622:28f58740b6f8, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Support infinite bounds in Circulation + fixes (#270, #266)

  • Support infinite capacities.
  • Bug fix in upperMap().
  • Fixes and improvements in the documentation.
File size: 24.4 KB
RevLine 
[399]1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
[440]5 * Copyright (C) 2003-2009
[399]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>
[622]24#include <limits>
[399]25
26///\ingroup max_flow
27///\file
[402]28///\brief Push-relabel algorithm for finding a feasible circulation.
[399]29///
30namespace lemon {
31
32  /// \brief Default traits class of Circulation class.
33  ///
34  /// Default traits class of Circulation class.
[610]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.
[492]40  template <typename GR, typename LM,
[610]41            typename UM, typename SM>
[399]42  struct CirculationDefaultTraits {
43
[402]44    /// \brief The type of the digraph the algorithm runs on.
[492]45    typedef GR Digraph;
[399]46
[610]47    /// \brief The type of the lower bound map.
[399]48    ///
[610]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;
[399]52
[610]53    /// \brief The type of the upper bound (capacity) map.
[399]54    ///
[610]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;
[399]59
[610]60    /// \brief The type of supply map.
[399]61    ///
[610]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;
[399]66
[402]67    /// \brief The type of the flow values.
[610]68    typedef typename SupplyMap::Value Flow;
[399]69
[402]70    /// \brief The type of the map that stores the flow values.
[399]71    ///
[402]72    /// The type of the map that stores the flow values.
[610]73    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
74    /// concept.
75    typedef typename Digraph::template ArcMap<Flow> FlowMap;
[399]76
77    /// \brief Instantiates a FlowMap.
78    ///
79    /// This function instantiates a \ref FlowMap.
[610]80    /// \param digraph The digraph for which we would like to define
[399]81    /// the flow map.
82    static FlowMap* createFlowMap(const Digraph& digraph) {
83      return new FlowMap(digraph);
84    }
85
[402]86    /// \brief The elevator type used by the algorithm.
[399]87    ///
[402]88    /// The elevator type used by the algorithm.
[399]89    ///
90    /// \sa Elevator
91    /// \sa LinkedElevator
92    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
93
94    /// \brief Instantiates an Elevator.
95    ///
[402]96    /// This function instantiates an \ref Elevator.
[610]97    /// \param digraph The digraph for which we would like to define
[399]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.
[610]107    typedef lemon::Tolerance<Flow> Tolerance;
[399]108
109  };
110
[402]111  /**
112     \brief Push-relabel algorithm for the network circulation problem.
[399]113
114     \ingroup max_flow
[610]115     This class implements a push-relabel algorithm for the \e network
116     \e circulation problem.
[402]117     It is to find a feasible circulation when lower and upper bounds
[610]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.
[402]121
[399]122     The exact formulation of this problem is the following.
[622]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$
[610]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.
[622]131     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
[610]132     solution of the following problem.
[402]133
[610]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
[622]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
[610]159     Note that this algorithm also provides a feasible solution for the
160     \ref min_cost_flow "minimum cost flow problem".
[402]161
[492]162     \tparam GR The type of the digraph the algorithm runs on.
[610]163     \tparam LM The type of the lower bound map. The default
[492]164     map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
[610]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
[492]168     \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
[399]169  */
[402]170#ifdef DOXYGEN
[492]171template< typename GR,
172          typename LM,
173          typename UM,
[610]174          typename SM,
[492]175          typename TR >
[402]176#else
[492]177template< typename GR,
178          typename LM = typename GR::template ArcMap<int>,
179          typename UM = LM,
[610]180          typename SM = typename GR::template NodeMap<typename UM::Value>,
181          typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
[402]182#endif
[399]183  class Circulation {
[402]184  public:
[399]185
[402]186    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
[492]187    typedef TR Traits;
[402]188    ///The type of the digraph the algorithm runs on.
[399]189    typedef typename Traits::Digraph Digraph;
[402]190    ///The type of the flow values.
[610]191    typedef typename Traits::Flow Flow;
[399]192
[610]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;
[402]199    ///The type of the flow map.
[399]200    typedef typename Traits::FlowMap FlowMap;
[402]201
202    ///The type of the elevator.
[399]203    typedef typename Traits::Elevator Elevator;
[402]204    ///The type of the tolerance.
[399]205    typedef typename Traits::Tolerance Tolerance;
206
[402]207  private:
208
209    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
[399]210
211    const Digraph &_g;
212    int _node_num;
213
[610]214    const LowerMap *_lo;
215    const UpperMap *_up;
216    const SupplyMap *_supply;
[399]217
218    FlowMap *_flow;
219    bool _local_flow;
220
221    Elevator* _level;
222    bool _local_level;
223
[610]224    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
[399]225    ExcessMap* _excess;
226
227    Tolerance _tol;
228    int _el;
229
230  public:
231
232    typedef Circulation Create;
233
[402]234    ///\name Named Template Parameters
[399]235
236    ///@{
237
[559]238    template <typename T>
[401]239    struct SetFlowMapTraits : public Traits {
[559]240      typedef T FlowMap;
[399]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
[402]251    /// type.
[559]252    template <typename T>
[401]253    struct SetFlowMap
[610]254      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
[559]255                           SetFlowMapTraits<T> > {
[610]256      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
[559]257                          SetFlowMapTraits<T> > Create;
[399]258    };
259
[559]260    template <typename T>
[401]261    struct SetElevatorTraits : public Traits {
[559]262      typedef T Elevator;
[399]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
[402]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
[559]278    template <typename T>
[401]279    struct SetElevator
[610]280      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
[559]281                           SetElevatorTraits<T> > {
[610]282      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
[559]283                          SetElevatorTraits<T> > Create;
[399]284    };
285
[559]286    template <typename T>
[401]287    struct SetStandardElevatorTraits : public Traits {
[559]288      typedef T Elevator;
[399]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
[402]295    /// Elevator type with automatic allocation
[399]296    ///
297    /// \ref named-templ-param "Named parameter" for setting Elevator
[402]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
[559]306    template <typename T>
[401]307    struct SetStandardElevator
[610]308      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
[559]309                       SetStandardElevatorTraits<T> > {
[610]310      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
[559]311                      SetStandardElevatorTraits<T> > Create;
[399]312    };
313
314    /// @}
315
316  protected:
317
318    Circulation() {}
319
320  public:
321
[610]322    /// Constructor.
[399]323
324    /// The constructor of the class.
[610]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) {}
[399]336
[402]337    /// Destructor.
[399]338    ~Circulation() {
339      destroyStructures();
340    }
341
[402]342
[399]343  private:
344
[622]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
[399]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
[610]382    /// Sets the lower bound map.
[399]383
[610]384    /// Sets the lower bound map.
[402]385    /// \return <tt>(*this)</tt>
[610]386    Circulation& lowerMap(const LowerMap& map) {
[399]387      _lo = &map;
388      return *this;
389    }
390
[610]391    /// Sets the upper bound (capacity) map.
[399]392
[610]393    /// Sets the upper bound (capacity) map.
[402]394    /// \return <tt>(*this)</tt>
[622]395    Circulation& upperMap(const UpperMap& map) {
[399]396      _up = &map;
397      return *this;
398    }
399
[610]400    /// Sets the supply map.
[399]401
[610]402    /// Sets the supply map.
[402]403    /// \return <tt>(*this)</tt>
[610]404    Circulation& supplyMap(const SupplyMap& map) {
405      _supply = &map;
[399]406      return *this;
407    }
408
[402]409    /// \brief Sets the flow map.
410    ///
[399]411    /// Sets the flow map.
[402]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>
[399]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
[402]426    /// \brief Sets the elevator used by algorithm.
[399]427    ///
[402]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>
[399]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
[402]443    /// \brief Returns a const reference to the elevator.
[399]444    ///
[402]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.
[420]449    const Elevator& elevator() const {
[399]450      return *_level;
451    }
452
[402]453    /// \brief Sets the tolerance used by algorithm.
454    ///
[399]455    /// Sets the tolerance used by algorithm.
456    Circulation& tolerance(const Tolerance& tolerance) const {
457      _tol = tolerance;
458      return *this;
459    }
460
[402]461    /// \brief Returns a const reference to the tolerance.
[399]462    ///
[402]463    /// Returns a const reference to the tolerance.
[399]464    const Tolerance& tolerance() const {
465      return tolerance;
466    }
467
[402]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.
[399]473
474    ///@{
475
476    /// Initializes the internal data structures.
477
[402]478    /// Initializes the internal data structures and sets all flow values
479    /// to the lower bound.
[399]480    void init()
481    {
[622]482      LEMON_DEBUG(checkBoundMaps(),
483        "Upper bounds must be greater or equal to the lower bounds");
484
[399]485      createStructures();
486
487      for(NodeIt n(_g);n!=INVALID;++n) {
[611]488        (*_excess)[n] = (*_supply)[n];
[399]489      }
490
491      for (ArcIt e(_g);e!=INVALID;++e) {
492        _flow->set(e, (*_lo)[e]);
[581]493        (*_excess)[_g.target(e)] += (*_flow)[e];
494        (*_excess)[_g.source(e)] -= (*_flow)[e];
[399]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
[402]508    /// Initializes the internal data structures using a greedy approach.
[399]509
[402]510    /// Initializes the internal data structures using a greedy approach
511    /// to construct the initial solution.
[399]512    void greedyInit()
513    {
[622]514      LEMON_DEBUG(checkBoundMaps(),
515        "Upper bounds must be greater or equal to the lower bounds");
516
[399]517      createStructures();
518
519      for(NodeIt n(_g);n!=INVALID;++n) {
[611]520        (*_excess)[n] = (*_supply)[n];
[399]521      }
522
523      for (ArcIt e(_g);e!=INVALID;++e) {
[622]524        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
[399]525          _flow->set(e, (*_up)[e]);
[581]526          (*_excess)[_g.target(e)] += (*_up)[e];
527          (*_excess)[_g.source(e)] -= (*_up)[e];
[622]528        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
[399]529          _flow->set(e, (*_lo)[e]);
[581]530          (*_excess)[_g.target(e)] += (*_lo)[e];
531          (*_excess)[_g.source(e)] -= (*_lo)[e];
[399]532        } else {
[610]533          Flow fc = -(*_excess)[_g.target(e)];
[399]534          _flow->set(e, fc);
[581]535          (*_excess)[_g.target(e)] = 0;
536          (*_excess)[_g.source(e)] -= fc;
[399]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
[402]549    ///Executes the algorithm
[399]550
[402]551    ///This function executes the algorithm.
552    ///
553    ///\return \c true if a feasible circulation is found.
[399]554    ///
555    ///\sa barrier()
[402]556    ///\sa barrierMap()
[399]557    bool start()
558    {
559
560      Node act;
561      Node bact=INVALID;
562      Node last_activated=INVALID;
563      while((act=_level->highestActive())!=INVALID) {
564        int actlevel=(*_level)[act];
565        int mlevel=_node_num;
[610]566        Flow exc=(*_excess)[act];
[399]567
568        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
569          Node v = _g.target(e);
[610]570          Flow fc=(*_up)[e]-(*_flow)[e];
[399]571          if(!_tol.positive(fc)) continue;
572          if((*_level)[v]<actlevel) {
573            if(!_tol.less(fc, exc)) {
574              _flow->set(e, (*_flow)[e] + exc);
[581]575              (*_excess)[v] += exc;
[399]576              if(!_level->active(v) && _tol.positive((*_excess)[v]))
577                _level->activate(v);
[581]578              (*_excess)[act] = 0;
[399]579              _level->deactivate(act);
580              goto next_l;
581            }
582            else {
583              _flow->set(e, (*_up)[e]);
[581]584              (*_excess)[v] += fc;
[399]585              if(!_level->active(v) && _tol.positive((*_excess)[v]))
586                _level->activate(v);
587              exc-=fc;
588            }
589          }
590          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
591        }
592        for(InArcIt e(_g,act);e!=INVALID; ++e) {
593          Node v = _g.source(e);
[610]594          Flow fc=(*_flow)[e]-(*_lo)[e];
[399]595          if(!_tol.positive(fc)) continue;
596          if((*_level)[v]<actlevel) {
597            if(!_tol.less(fc, exc)) {
598              _flow->set(e, (*_flow)[e] - exc);
[581]599              (*_excess)[v] += exc;
[399]600              if(!_level->active(v) && _tol.positive((*_excess)[v]))
601                _level->activate(v);
[581]602              (*_excess)[act] = 0;
[399]603              _level->deactivate(act);
604              goto next_l;
605            }
606            else {
607              _flow->set(e, (*_lo)[e]);
[581]608              (*_excess)[v] += fc;
[399]609              if(!_level->active(v) && _tol.positive((*_excess)[v]))
610                _level->activate(v);
611              exc-=fc;
612            }
613          }
614          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
615        }
616
[581]617        (*_excess)[act] = exc;
[399]618        if(!_tol.positive(exc)) _level->deactivate(act);
619        else if(mlevel==_node_num) {
620          _level->liftHighestActiveToTop();
621          _el = _node_num;
622          return false;
623        }
624        else {
625          _level->liftHighestActive(mlevel+1);
626          if(_level->onLevel(actlevel)==0) {
627            _el = actlevel;
628            return false;
629          }
630        }
631      next_l:
632        ;
633      }
634      return true;
635    }
636
[402]637    /// Runs the algorithm.
[399]638
[402]639    /// This function runs the algorithm.
640    ///
641    /// \return \c true if a feasible circulation is found.
642    ///
643    /// \note Apart from the return value, c.run() is just a shortcut of
644    /// the following code.
[399]645    /// \code
[402]646    ///   c.greedyInit();
647    ///   c.start();
[399]648    /// \endcode
649    bool run() {
650      greedyInit();
651      return start();
652    }
653
654    /// @}
655
656    /// \name Query Functions
[402]657    /// The results of the circulation algorithm can be obtained using
658    /// these functions.\n
659    /// Either \ref run() or \ref start() should be called before
660    /// using them.
[399]661
662    ///@{
663
[402]664    /// \brief Returns the flow on the given arc.
665    ///
666    /// Returns the flow on the given arc.
667    ///
668    /// \pre Either \ref run() or \ref init() must be called before
669    /// using this function.
[610]670    Flow flow(const Arc& arc) const {
[402]671      return (*_flow)[arc];
672    }
673
674    /// \brief Returns a const reference to the flow map.
675    ///
676    /// Returns a const reference to the arc map storing the found flow.
677    ///
678    /// \pre Either \ref run() or \ref init() must be called before
679    /// using this function.
[420]680    const FlowMap& flowMap() const {
[402]681      return *_flow;
682    }
683
[399]684    /**
[402]685       \brief Returns \c true if the given node is in a barrier.
686
[399]687       Barrier is a set \e B of nodes for which
[402]688
[610]689       \f[ \sum_{uv\in A: u\in B} upper(uv) -
690           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
[402]691
692       holds. The existence of a set with this property prooves that a
693       feasible circualtion cannot exist.
694
695       This function returns \c true if the given node is in the found
696       barrier. If a feasible circulation is found, the function
697       gives back \c false for every node.
698
699       \pre Either \ref run() or \ref init() must be called before
700       using this function.
701
702       \sa barrierMap()
[399]703       \sa checkBarrier()
704    */
[420]705    bool barrier(const Node& node) const
[402]706    {
707      return (*_level)[node] >= _el;
708    }
709
710    /// \brief Gives back a barrier.
711    ///
712    /// This function sets \c bar to the characteristic vector of the
713    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
714    /// node map with \c bool (or convertible) value type.
715    ///
716    /// If a feasible circulation is found, the function gives back an
717    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
718    ///
719    /// \note This function calls \ref barrier() for each node,
[559]720    /// so it runs in O(n) time.
[402]721    ///
722    /// \pre Either \ref run() or \ref init() must be called before
723    /// using this function.
724    ///
725    /// \sa barrier()
726    /// \sa checkBarrier()
727    template<class BarrierMap>
[420]728    void barrierMap(BarrierMap &bar) const
[399]729    {
730      for(NodeIt n(_g);n!=INVALID;++n)
731        bar.set(n, (*_level)[n] >= _el);
732    }
733
734    /// @}
735
736    /// \name Checker Functions
[402]737    /// The feasibility of the results can be checked using
738    /// these functions.\n
739    /// Either \ref run() or \ref start() should be called before
740    /// using them.
[399]741
742    ///@{
743
[402]744    ///Check if the found flow is a feasible circulation
745
746    ///Check if the found flow is a feasible circulation,
747    ///
[420]748    bool checkFlow() const {
[399]749      for(ArcIt e(_g);e!=INVALID;++e)
750        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
751      for(NodeIt n(_g);n!=INVALID;++n)
752        {
[610]753          Flow dif=-(*_supply)[n];
[399]754          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
755          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
756          if(_tol.negative(dif)) return false;
757        }
758      return true;
759    }
760
761    ///Check whether or not the last execution provides a barrier
762
[402]763    ///Check whether or not the last execution provides a barrier.
[399]764    ///\sa barrier()
[402]765    ///\sa barrierMap()
[420]766    bool checkBarrier() const
[399]767    {
[610]768      Flow delta=0;
[622]769      Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
770        std::numeric_limits<Flow>::infinity() :
771        std::numeric_limits<Flow>::max();
[399]772      for(NodeIt n(_g);n!=INVALID;++n)
773        if(barrier(n))
[610]774          delta-=(*_supply)[n];
[399]775      for(ArcIt e(_g);e!=INVALID;++e)
776        {
777          Node s=_g.source(e);
778          Node t=_g.target(e);
[622]779          if(barrier(s)&&!barrier(t)) {
780            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
781            delta+=(*_up)[e];
782          }
[399]783          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
784        }
785      return _tol.negative(delta);
786    }
787
788    /// @}
789
790  };
791
792}
793
794#endif
Note: See TracBrowser for help on using the repository browser.