COIN-OR::LEMON - Graph Library

source: lemon/lemon/circulation.h @ 736:86c49553fea5

Last change on this file since 736:86c49553fea5 was 736:86c49553fea5, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Test file + doc improvements (#307)

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