COIN-OR::LEMON - Graph Library

source: lemon/lemon/circulation.h @ 1064:40bbb450143e

Last change on this file since 1064:40bbb450143e was 735:1f08e846df29, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Bug fix in Preflow and Circulation (#307)

File size: 24.4 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 algorithm.
454    ///
455    /// Sets the tolerance used by algorithm.
456    Circulation& tolerance(const Tolerance& tolerance) {
457      _tol = tolerance;
458      return *this;
459    }
460
461    /// \brief Returns a const reference to the tolerance.
462    ///
463    /// Returns a const reference to the tolerance.
464    const Tolerance& tolerance() const {
465      return _tol;
466    }
467
468    /// \name Execution Control
469    /// The simplest way to execute the algorithm is to call \ref run().\n
470    /// If you need more control on the initial solution or the execution,
471    /// first you have to call one of the \ref init() functions, then
472    /// the \ref start() function.
473
474    ///@{
475
476    /// Initializes the internal data structures.
477
478    /// Initializes the internal data structures and sets all flow values
479    /// to the lower bound.
480    void init()
481    {
482      LEMON_DEBUG(checkBoundMaps(),
483        "Upper bounds must be greater or equal to the lower bounds");
484
485      createStructures();
486
487      for(NodeIt n(_g);n!=INVALID;++n) {
488        (*_excess)[n] = (*_supply)[n];
489      }
490
491      for (ArcIt e(_g);e!=INVALID;++e) {
492        _flow->set(e, (*_lo)[e]);
493        (*_excess)[_g.target(e)] += (*_flow)[e];
494        (*_excess)[_g.source(e)] -= (*_flow)[e];
495      }
496
497      // global relabeling tested, but in general case it provides
498      // worse performance for random digraphs
499      _level->initStart();
500      for(NodeIt n(_g);n!=INVALID;++n)
501        _level->initAddItem(n);
502      _level->initFinish();
503      for(NodeIt n(_g);n!=INVALID;++n)
504        if(_tol.positive((*_excess)[n]))
505          _level->activate(n);
506    }
507
508    /// Initializes the internal data structures using a greedy approach.
509
510    /// Initializes the internal data structures using a greedy approach
511    /// to construct the initial solution.
512    void greedyInit()
513    {
514      LEMON_DEBUG(checkBoundMaps(),
515        "Upper bounds must be greater or equal to the lower bounds");
516
517      createStructures();
518
519      for(NodeIt n(_g);n!=INVALID;++n) {
520        (*_excess)[n] = (*_supply)[n];
521      }
522
523      for (ArcIt e(_g);e!=INVALID;++e) {
524        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
525          _flow->set(e, (*_up)[e]);
526          (*_excess)[_g.target(e)] += (*_up)[e];
527          (*_excess)[_g.source(e)] -= (*_up)[e];
528        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
529          _flow->set(e, (*_lo)[e]);
530          (*_excess)[_g.target(e)] += (*_lo)[e];
531          (*_excess)[_g.source(e)] -= (*_lo)[e];
532        } else {
533          Value fc = -(*_excess)[_g.target(e)];
534          _flow->set(e, fc);
535          (*_excess)[_g.target(e)] = 0;
536          (*_excess)[_g.source(e)] -= fc;
537        }
538      }
539
540      _level->initStart();
541      for(NodeIt n(_g);n!=INVALID;++n)
542        _level->initAddItem(n);
543      _level->initFinish();
544      for(NodeIt n(_g);n!=INVALID;++n)
545        if(_tol.positive((*_excess)[n]))
546          _level->activate(n);
547    }
548
549    ///Executes the algorithm
550
551    ///This function executes the algorithm.
552    ///
553    ///\return \c true if a feasible circulation is found.
554    ///
555    ///\sa barrier()
556    ///\sa barrierMap()
557    bool start()
558    {
559
560      Node act;
561      Node bact=INVALID;
562      Node last_activated=INVALID;
563      while((act=_level->highestActive())!=INVALID) {
564        int actlevel=(*_level)[act];
565        int mlevel=_node_num;
566        Value exc=(*_excess)[act];
567
568        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
569          Node v = _g.target(e);
570          Value fc=(*_up)[e]-(*_flow)[e];
571          if(!_tol.positive(fc)) continue;
572          if((*_level)[v]<actlevel) {
573            if(!_tol.less(fc, exc)) {
574              _flow->set(e, (*_flow)[e] + exc);
575              (*_excess)[v] += exc;
576              if(!_level->active(v) && _tol.positive((*_excess)[v]))
577                _level->activate(v);
578              (*_excess)[act] = 0;
579              _level->deactivate(act);
580              goto next_l;
581            }
582            else {
583              _flow->set(e, (*_up)[e]);
584              (*_excess)[v] += fc;
585              if(!_level->active(v) && _tol.positive((*_excess)[v]))
586                _level->activate(v);
587              exc-=fc;
588            }
589          }
590          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
591        }
592        for(InArcIt e(_g,act);e!=INVALID; ++e) {
593          Node v = _g.source(e);
594          Value fc=(*_flow)[e]-(*_lo)[e];
595          if(!_tol.positive(fc)) continue;
596          if((*_level)[v]<actlevel) {
597            if(!_tol.less(fc, exc)) {
598              _flow->set(e, (*_flow)[e] - exc);
599              (*_excess)[v] += exc;
600              if(!_level->active(v) && _tol.positive((*_excess)[v]))
601                _level->activate(v);
602              (*_excess)[act] = 0;
603              _level->deactivate(act);
604              goto next_l;
605            }
606            else {
607              _flow->set(e, (*_lo)[e]);
608              (*_excess)[v] += fc;
609              if(!_level->active(v) && _tol.positive((*_excess)[v]))
610                _level->activate(v);
611              exc-=fc;
612            }
613          }
614          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
615        }
616
617        (*_excess)[act] = exc;
618        if(!_tol.positive(exc)) _level->deactivate(act);
619        else if(mlevel==_node_num) {
620          _level->liftHighestActiveToTop();
621          _el = _node_num;
622          return false;
623        }
624        else {
625          _level->liftHighestActive(mlevel+1);
626          if(_level->onLevel(actlevel)==0) {
627            _el = actlevel;
628            return false;
629          }
630        }
631      next_l:
632        ;
633      }
634      return true;
635    }
636
637    /// Runs the algorithm.
638
639    /// This function runs the algorithm.
640    ///
641    /// \return \c true if a feasible circulation is found.
642    ///
643    /// \note Apart from the return value, c.run() is just a shortcut of
644    /// the following code.
645    /// \code
646    ///   c.greedyInit();
647    ///   c.start();
648    /// \endcode
649    bool run() {
650      greedyInit();
651      return start();
652    }
653
654    /// @}
655
656    /// \name Query Functions
657    /// The results of the circulation algorithm can be obtained using
658    /// these functions.\n
659    /// Either \ref run() or \ref start() should be called before
660    /// using them.
661
662    ///@{
663
664    /// \brief Returns the flow value on the given arc.
665    ///
666    /// Returns the flow value on the given arc.
667    ///
668    /// \pre Either \ref run() or \ref init() must be called before
669    /// using this function.
670    Value flow(const Arc& arc) const {
671      return (*_flow)[arc];
672    }
673
674    /// \brief Returns a const reference to the flow map.
675    ///
676    /// Returns a const reference to the arc map storing the found flow.
677    ///
678    /// \pre Either \ref run() or \ref init() must be called before
679    /// using this function.
680    const FlowMap& flowMap() const {
681      return *_flow;
682    }
683
684    /**
685       \brief Returns \c true if the given node is in a barrier.
686
687       Barrier is a set \e B of nodes for which
688
689       \f[ \sum_{uv\in A: u\in B} upper(uv) -
690           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
691
692       holds. The existence of a set with this property prooves that a
693       feasible circualtion cannot exist.
694
695       This function returns \c true if the given node is in the found
696       barrier. If a feasible circulation is found, the function
697       gives back \c false for every node.
698
699       \pre Either \ref run() or \ref init() must be called before
700       using this function.
701
702       \sa barrierMap()
703       \sa checkBarrier()
704    */
705    bool barrier(const Node& node) const
706    {
707      return (*_level)[node] >= _el;
708    }
709
710    /// \brief Gives back a barrier.
711    ///
712    /// This function sets \c bar to the characteristic vector of the
713    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
714    /// node map with \c bool (or convertible) value type.
715    ///
716    /// If a feasible circulation is found, the function gives back an
717    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
718    ///
719    /// \note This function calls \ref barrier() for each node,
720    /// so it runs in O(n) time.
721    ///
722    /// \pre Either \ref run() or \ref init() must be called before
723    /// using this function.
724    ///
725    /// \sa barrier()
726    /// \sa checkBarrier()
727    template<class BarrierMap>
728    void barrierMap(BarrierMap &bar) const
729    {
730      for(NodeIt n(_g);n!=INVALID;++n)
731        bar.set(n, (*_level)[n] >= _el);
732    }
733
734    /// @}
735
736    /// \name Checker Functions
737    /// The feasibility of the results can be checked using
738    /// these functions.\n
739    /// Either \ref run() or \ref start() should be called before
740    /// using them.
741
742    ///@{
743
744    ///Check if the found flow is a feasible circulation
745
746    ///Check if the found flow is a feasible circulation,
747    ///
748    bool checkFlow() const {
749      for(ArcIt e(_g);e!=INVALID;++e)
750        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
751      for(NodeIt n(_g);n!=INVALID;++n)
752        {
753          Value dif=-(*_supply)[n];
754          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
755          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
756          if(_tol.negative(dif)) return false;
757        }
758      return true;
759    }
760
761    ///Check whether or not the last execution provides a barrier
762
763    ///Check whether or not the last execution provides a barrier.
764    ///\sa barrier()
765    ///\sa barrierMap()
766    bool checkBarrier() const
767    {
768      Value delta=0;
769      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
770        std::numeric_limits<Value>::infinity() :
771        std::numeric_limits<Value>::max();
772      for(NodeIt n(_g);n!=INVALID;++n)
773        if(barrier(n))
774          delta-=(*_supply)[n];
775      for(ArcIt e(_g);e!=INVALID;++e)
776        {
777          Node s=_g.source(e);
778          Node t=_g.target(e);
779          if(barrier(s)&&!barrier(t)) {
780            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
781            delta+=(*_up)[e];
782          }
783          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
784        }
785      return _tol.negative(delta);
786    }
787
788    /// @}
789
790  };
791
792}
793
794#endif
Note: See TracBrowser for help on using the repository browser.