COIN-OR::LEMON - Graph Library

source: lemon/lemon/circulation.h @ 1232:fc3854d936f7

Last change on this file since 1232:fc3854d936f7 was 1159:7fdaa05a69a1, checked in by Alpar Juttner <alpar@…>, 12 years ago

Merge #449 to branches >=1.2

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      while((act=_level->highestActive())!=INVALID) {
576        int actlevel=(*_level)[act];
577        int mlevel=_node_num;
578        Value exc=(*_excess)[act];
579
580        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
581          Node v = _g.target(e);
582          Value fc=(*_up)[e]-(*_flow)[e];
583          if(!_tol.positive(fc)) continue;
584          if((*_level)[v]<actlevel) {
585            if(!_tol.less(fc, exc)) {
586              _flow->set(e, (*_flow)[e] + exc);
587              (*_excess)[v] += exc;
588              if(!_level->active(v) && _tol.positive((*_excess)[v]))
589                _level->activate(v);
590              (*_excess)[act] = 0;
591              _level->deactivate(act);
592              goto next_l;
593            }
594            else {
595              _flow->set(e, (*_up)[e]);
596              (*_excess)[v] += fc;
597              if(!_level->active(v) && _tol.positive((*_excess)[v]))
598                _level->activate(v);
599              exc-=fc;
600            }
601          }
602          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
603        }
604        for(InArcIt e(_g,act);e!=INVALID; ++e) {
605          Node v = _g.source(e);
606          Value fc=(*_flow)[e]-(*_lo)[e];
607          if(!_tol.positive(fc)) continue;
608          if((*_level)[v]<actlevel) {
609            if(!_tol.less(fc, exc)) {
610              _flow->set(e, (*_flow)[e] - exc);
611              (*_excess)[v] += exc;
612              if(!_level->active(v) && _tol.positive((*_excess)[v]))
613                _level->activate(v);
614              (*_excess)[act] = 0;
615              _level->deactivate(act);
616              goto next_l;
617            }
618            else {
619              _flow->set(e, (*_lo)[e]);
620              (*_excess)[v] += fc;
621              if(!_level->active(v) && _tol.positive((*_excess)[v]))
622                _level->activate(v);
623              exc-=fc;
624            }
625          }
626          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
627        }
628
629        (*_excess)[act] = exc;
630        if(!_tol.positive(exc)) _level->deactivate(act);
631        else if(mlevel==_node_num) {
632          _level->liftHighestActiveToTop();
633          _el = _node_num;
634          return false;
635        }
636        else {
637          _level->liftHighestActive(mlevel+1);
638          if(_level->onLevel(actlevel)==0) {
639            _el = actlevel;
640            return false;
641          }
642        }
643      next_l:
644        ;
645      }
646      return true;
647    }
648
649    /// Runs the algorithm.
650
651    /// This function runs the algorithm.
652    ///
653    /// \return \c true if a feasible circulation is found.
654    ///
655    /// \note Apart from the return value, c.run() is just a shortcut of
656    /// the following code.
657    /// \code
658    ///   c.greedyInit();
659    ///   c.start();
660    /// \endcode
661    bool run() {
662      greedyInit();
663      return start();
664    }
665
666    /// @}
667
668    /// \name Query Functions
669    /// The results of the circulation algorithm can be obtained using
670    /// these functions.\n
671    /// Either \ref run() or \ref start() should be called before
672    /// using them.
673
674    ///@{
675
676    /// \brief Returns the flow value on the given arc.
677    ///
678    /// Returns the flow value on the given arc.
679    ///
680    /// \pre Either \ref run() or \ref init() must be called before
681    /// using this function.
682    Value flow(const Arc& arc) const {
683      return (*_flow)[arc];
684    }
685
686    /// \brief Returns a const reference to the flow map.
687    ///
688    /// Returns a const reference to the arc map storing the found flow.
689    ///
690    /// \pre Either \ref run() or \ref init() must be called before
691    /// using this function.
692    const FlowMap& flowMap() const {
693      return *_flow;
694    }
695
696    /**
697       \brief Returns \c true if the given node is in a barrier.
698
699       Barrier is a set \e B of nodes for which
700
701       \f[ \sum_{uv\in A: u\in B} upper(uv) -
702           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
703
704       holds. The existence of a set with this property prooves that a
705       feasible circualtion cannot exist.
706
707       This function returns \c true if the given node is in the found
708       barrier. If a feasible circulation is found, the function
709       gives back \c false for every node.
710
711       \pre Either \ref run() or \ref init() must be called before
712       using this function.
713
714       \sa barrierMap()
715       \sa checkBarrier()
716    */
717    bool barrier(const Node& node) const
718    {
719      return (*_level)[node] >= _el;
720    }
721
722    /// \brief Gives back a barrier.
723    ///
724    /// This function sets \c bar to the characteristic vector of the
725    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
726    /// node map with \c bool (or convertible) value type.
727    ///
728    /// If a feasible circulation is found, the function gives back an
729    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
730    ///
731    /// \note This function calls \ref barrier() for each node,
732    /// so it runs in O(n) time.
733    ///
734    /// \pre Either \ref run() or \ref init() must be called before
735    /// using this function.
736    ///
737    /// \sa barrier()
738    /// \sa checkBarrier()
739    template<class BarrierMap>
740    void barrierMap(BarrierMap &bar) const
741    {
742      for(NodeIt n(_g);n!=INVALID;++n)
743        bar.set(n, (*_level)[n] >= _el);
744    }
745
746    /// @}
747
748    /// \name Checker Functions
749    /// The feasibility of the results can be checked using
750    /// these functions.\n
751    /// Either \ref run() or \ref start() should be called before
752    /// using them.
753
754    ///@{
755
756    ///Check if the found flow is a feasible circulation
757
758    ///Check if the found flow is a feasible circulation,
759    ///
760    bool checkFlow() const {
761      for(ArcIt e(_g);e!=INVALID;++e)
762        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
763      for(NodeIt n(_g);n!=INVALID;++n)
764        {
765          Value dif=-(*_supply)[n];
766          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
767          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
768          if(_tol.negative(dif)) return false;
769        }
770      return true;
771    }
772
773    ///Check whether or not the last execution provides a barrier
774
775    ///Check whether or not the last execution provides a barrier.
776    ///\sa barrier()
777    ///\sa barrierMap()
778    bool checkBarrier() const
779    {
780      Value delta=0;
781      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
782        std::numeric_limits<Value>::infinity() :
783        std::numeric_limits<Value>::max();
784      for(NodeIt n(_g);n!=INVALID;++n)
785        if(barrier(n))
786          delta-=(*_supply)[n];
787      for(ArcIt e(_g);e!=INVALID;++e)
788        {
789          Node s=_g.source(e);
790          Node t=_g.target(e);
791          if(barrier(s)&&!barrier(t)) {
792            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
793            delta+=(*_up)[e];
794          }
795          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
796        }
797      return _tol.negative(delta);
798    }
799
800    /// @}
801
802  };
803
804}
805
806#endif
Note: See TracBrowser for help on using the repository browser.