COIN-OR::LEMON - Graph Library

source: lemon/lemon/circulation.h @ 1337:4add05447ca0

Last change on this file since 1337:4add05447ca0 was 1270:dceba191c00d, checked in by Alpar Juttner <alpar@…>, 11 years ago

Apply unify-sources.sh to the source tree

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