COIN-OR::LEMON - Graph Library

source: lemon/lemon/circulation.h @ 1127:be7dd3a8d6a3

Last change on this file since 1127:be7dd3a8d6a3 was 956:141f9c0db4a3, checked in by Alpar Juttner <alpar@…>, 14 years ago

Unify the sources (#339)

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