COIN-OR::LEMON - Graph Library

source: lemon/lemon/preflow.h @ 760:4ac30454f1c1

Last change on this file since 760:4ac30454f1c1 was 760:4ac30454f1c1, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Small doc improvements

File size: 29.7 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_PREFLOW_H
20#define LEMON_PREFLOW_H
21
22#include <lemon/tolerance.h>
23#include <lemon/elevator.h>
24
25/// \file
26/// \ingroup max_flow
27/// \brief Implementation of the preflow algorithm.
28
29namespace lemon {
30
31  /// \brief Default traits class of Preflow class.
32  ///
33  /// Default traits class of Preflow class.
34  /// \tparam GR Digraph type.
35  /// \tparam CAP Capacity map type.
36  template <typename GR, typename CAP>
37  struct PreflowDefaultTraits {
38
39    /// \brief The type of the digraph the algorithm runs on.
40    typedef GR Digraph;
41
42    /// \brief The type of the map that stores the arc capacities.
43    ///
44    /// The type of the map that stores the arc capacities.
45    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
46    typedef CAP CapacityMap;
47
48    /// \brief The type of the flow values.
49    typedef typename CapacityMap::Value Value;
50
51    /// \brief The type of the map that stores the flow values.
52    ///
53    /// The type of the map that stores the flow values.
54    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55#ifdef DOXYGEN
56    typedef GR::ArcMap<Value> FlowMap;
57#else
58    typedef typename Digraph::template ArcMap<Value> FlowMap;
59#endif
60
61    /// \brief Instantiates a FlowMap.
62    ///
63    /// This function instantiates a \ref FlowMap.
64    /// \param digraph The digraph for which we would like to define
65    /// the flow map.
66    static FlowMap* createFlowMap(const Digraph& digraph) {
67      return new FlowMap(digraph);
68    }
69
70    /// \brief The elevator type used by Preflow algorithm.
71    ///
72    /// The elevator type used by Preflow algorithm.
73    ///
74    /// \sa Elevator, LinkedElevator
75#ifdef DOXYGEN
76    typedef lemon::Elevator<GR, GR::Node> Elevator;
77#else
78    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
79#endif
80
81    /// \brief Instantiates an Elevator.
82    ///
83    /// This function instantiates an \ref Elevator.
84    /// \param digraph The digraph for which we would like to define
85    /// the elevator.
86    /// \param max_level The maximum level of the elevator.
87    static Elevator* createElevator(const Digraph& digraph, int max_level) {
88      return new Elevator(digraph, max_level);
89    }
90
91    /// \brief The tolerance used by the algorithm
92    ///
93    /// The tolerance used by the algorithm to handle inexact computation.
94    typedef lemon::Tolerance<Value> Tolerance;
95
96  };
97
98
99  /// \ingroup max_flow
100  ///
101  /// \brief %Preflow algorithm class.
102  ///
103  /// This class provides an implementation of Goldberg-Tarjan's \e preflow
104  /// \e push-relabel algorithm producing a \ref max_flow
105  /// "flow of maximum value" in a digraph.
106  /// The preflow algorithms are the fastest known maximum
107  /// flow algorithms. The current implementation use a mixture of the
108  /// \e "highest label" and the \e "bound decrease" heuristics.
109  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
110  ///
111  /// The algorithm consists of two phases. After the first phase
112  /// the maximum flow value and the minimum cut is obtained. The
113  /// second phase constructs a feasible maximum flow on each arc.
114  ///
115  /// \tparam GR The type of the digraph the algorithm runs on.
116  /// \tparam CAP The type of the capacity map. The default map
117  /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
118#ifdef DOXYGEN
119  template <typename GR, typename CAP, typename TR>
120#else
121  template <typename GR,
122            typename CAP = typename GR::template ArcMap<int>,
123            typename TR = PreflowDefaultTraits<GR, CAP> >
124#endif
125  class Preflow {
126  public:
127
128    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
129    typedef TR Traits;
130    ///The type of the digraph the algorithm runs on.
131    typedef typename Traits::Digraph Digraph;
132    ///The type of the capacity map.
133    typedef typename Traits::CapacityMap CapacityMap;
134    ///The type of the flow values.
135    typedef typename Traits::Value Value;
136
137    ///The type of the flow map.
138    typedef typename Traits::FlowMap FlowMap;
139    ///The type of the elevator.
140    typedef typename Traits::Elevator Elevator;
141    ///The type of the tolerance.
142    typedef typename Traits::Tolerance Tolerance;
143
144  private:
145
146    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
147
148    const Digraph& _graph;
149    const CapacityMap* _capacity;
150
151    int _node_num;
152
153    Node _source, _target;
154
155    FlowMap* _flow;
156    bool _local_flow;
157
158    Elevator* _level;
159    bool _local_level;
160
161    typedef typename Digraph::template NodeMap<Value> ExcessMap;
162    ExcessMap* _excess;
163
164    Tolerance _tolerance;
165
166    bool _phase;
167
168
169    void createStructures() {
170      _node_num = countNodes(_graph);
171
172      if (!_flow) {
173        _flow = Traits::createFlowMap(_graph);
174        _local_flow = true;
175      }
176      if (!_level) {
177        _level = Traits::createElevator(_graph, _node_num);
178        _local_level = true;
179      }
180      if (!_excess) {
181        _excess = new ExcessMap(_graph);
182      }
183    }
184
185    void destroyStructures() {
186      if (_local_flow) {
187        delete _flow;
188      }
189      if (_local_level) {
190        delete _level;
191      }
192      if (_excess) {
193        delete _excess;
194      }
195    }
196
197  public:
198
199    typedef Preflow Create;
200
201    ///\name Named Template Parameters
202
203    ///@{
204
205    template <typename T>
206    struct SetFlowMapTraits : public Traits {
207      typedef T FlowMap;
208      static FlowMap *createFlowMap(const Digraph&) {
209        LEMON_ASSERT(false, "FlowMap is not initialized");
210        return 0; // ignore warnings
211      }
212    };
213
214    /// \brief \ref named-templ-param "Named parameter" for setting
215    /// FlowMap type
216    ///
217    /// \ref named-templ-param "Named parameter" for setting FlowMap
218    /// type.
219    template <typename T>
220    struct SetFlowMap
221      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
222      typedef Preflow<Digraph, CapacityMap,
223                      SetFlowMapTraits<T> > Create;
224    };
225
226    template <typename T>
227    struct SetElevatorTraits : public Traits {
228      typedef T Elevator;
229      static Elevator *createElevator(const Digraph&, int) {
230        LEMON_ASSERT(false, "Elevator is not initialized");
231        return 0; // ignore warnings
232      }
233    };
234
235    /// \brief \ref named-templ-param "Named parameter" for setting
236    /// Elevator type
237    ///
238    /// \ref named-templ-param "Named parameter" for setting Elevator
239    /// type. If this named parameter is used, then an external
240    /// elevator object must be passed to the algorithm using the
241    /// \ref elevator(Elevator&) "elevator()" function before calling
242    /// \ref run() or \ref init().
243    /// \sa SetStandardElevator
244    template <typename T>
245    struct SetElevator
246      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
247      typedef Preflow<Digraph, CapacityMap,
248                      SetElevatorTraits<T> > Create;
249    };
250
251    template <typename T>
252    struct SetStandardElevatorTraits : public Traits {
253      typedef T Elevator;
254      static Elevator *createElevator(const Digraph& digraph, int max_level) {
255        return new Elevator(digraph, max_level);
256      }
257    };
258
259    /// \brief \ref named-templ-param "Named parameter" for setting
260    /// Elevator type with automatic allocation
261    ///
262    /// \ref named-templ-param "Named parameter" for setting Elevator
263    /// type with automatic allocation.
264    /// The Elevator should have standard constructor interface to be
265    /// able to automatically created by the algorithm (i.e. the
266    /// digraph and the maximum level should be passed to it).
267    /// However an external elevator object could also be passed to the
268    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
269    /// before calling \ref run() or \ref init().
270    /// \sa SetElevator
271    template <typename T>
272    struct SetStandardElevator
273      : public Preflow<Digraph, CapacityMap,
274                       SetStandardElevatorTraits<T> > {
275      typedef Preflow<Digraph, CapacityMap,
276                      SetStandardElevatorTraits<T> > Create;
277    };
278
279    /// @}
280
281  protected:
282
283    Preflow() {}
284
285  public:
286
287
288    /// \brief The constructor of the class.
289    ///
290    /// The constructor of the class.
291    /// \param digraph The digraph the algorithm runs on.
292    /// \param capacity The capacity of the arcs.
293    /// \param source The source node.
294    /// \param target The target node.
295    Preflow(const Digraph& digraph, const CapacityMap& capacity,
296            Node source, Node target)
297      : _graph(digraph), _capacity(&capacity),
298        _node_num(0), _source(source), _target(target),
299        _flow(0), _local_flow(false),
300        _level(0), _local_level(false),
301        _excess(0), _tolerance(), _phase() {}
302
303    /// \brief Destructor.
304    ///
305    /// Destructor.
306    ~Preflow() {
307      destroyStructures();
308    }
309
310    /// \brief Sets the capacity map.
311    ///
312    /// Sets the capacity map.
313    /// \return <tt>(*this)</tt>
314    Preflow& capacityMap(const CapacityMap& map) {
315      _capacity = &map;
316      return *this;
317    }
318
319    /// \brief Sets the flow map.
320    ///
321    /// Sets the flow map.
322    /// If you don't use this function before calling \ref run() or
323    /// \ref init(), an instance will be allocated automatically.
324    /// The destructor deallocates this automatically allocated map,
325    /// of course.
326    /// \return <tt>(*this)</tt>
327    Preflow& flowMap(FlowMap& map) {
328      if (_local_flow) {
329        delete _flow;
330        _local_flow = false;
331      }
332      _flow = &map;
333      return *this;
334    }
335
336    /// \brief Sets the source node.
337    ///
338    /// Sets the source node.
339    /// \return <tt>(*this)</tt>
340    Preflow& source(const Node& node) {
341      _source = node;
342      return *this;
343    }
344
345    /// \brief Sets the target node.
346    ///
347    /// Sets the target node.
348    /// \return <tt>(*this)</tt>
349    Preflow& target(const Node& node) {
350      _target = node;
351      return *this;
352    }
353
354    /// \brief Sets the elevator used by algorithm.
355    ///
356    /// Sets the elevator used by algorithm.
357    /// If you don't use this function before calling \ref run() or
358    /// \ref init(), an instance will be allocated automatically.
359    /// The destructor deallocates this automatically allocated elevator,
360    /// of course.
361    /// \return <tt>(*this)</tt>
362    Preflow& elevator(Elevator& elevator) {
363      if (_local_level) {
364        delete _level;
365        _local_level = false;
366      }
367      _level = &elevator;
368      return *this;
369    }
370
371    /// \brief Returns a const reference to the elevator.
372    ///
373    /// Returns a const reference to the elevator.
374    ///
375    /// \pre Either \ref run() or \ref init() must be called before
376    /// using this function.
377    const Elevator& elevator() const {
378      return *_level;
379    }
380
381    /// \brief Sets the tolerance used by algorithm.
382    ///
383    /// Sets the tolerance used by algorithm.
384    Preflow& tolerance(const Tolerance& tolerance) const {
385      _tolerance = tolerance;
386      return *this;
387    }
388
389    /// \brief Returns a const reference to the tolerance.
390    ///
391    /// Returns a const reference to the tolerance.
392    const Tolerance& tolerance() const {
393      return tolerance;
394    }
395
396    /// \name Execution Control
397    /// The simplest way to execute the preflow algorithm is to use
398    /// \ref run() or \ref runMinCut().\n
399    /// If you need better control on the initial solution or the execution,
400    /// you have to call one of the \ref init() functions first, then
401    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
402
403    ///@{
404
405    /// \brief Initializes the internal data structures.
406    ///
407    /// Initializes the internal data structures and sets the initial
408    /// flow to zero on each arc.
409    void init() {
410      createStructures();
411
412      _phase = true;
413      for (NodeIt n(_graph); n != INVALID; ++n) {
414        (*_excess)[n] = 0;
415      }
416
417      for (ArcIt e(_graph); e != INVALID; ++e) {
418        _flow->set(e, 0);
419      }
420
421      typename Digraph::template NodeMap<bool> reached(_graph, false);
422
423      _level->initStart();
424      _level->initAddItem(_target);
425
426      std::vector<Node> queue;
427      reached[_source] = true;
428
429      queue.push_back(_target);
430      reached[_target] = true;
431      while (!queue.empty()) {
432        _level->initNewLevel();
433        std::vector<Node> nqueue;
434        for (int i = 0; i < int(queue.size()); ++i) {
435          Node n = queue[i];
436          for (InArcIt e(_graph, n); e != INVALID; ++e) {
437            Node u = _graph.source(e);
438            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
439              reached[u] = true;
440              _level->initAddItem(u);
441              nqueue.push_back(u);
442            }
443          }
444        }
445        queue.swap(nqueue);
446      }
447      _level->initFinish();
448
449      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
450        if (_tolerance.positive((*_capacity)[e])) {
451          Node u = _graph.target(e);
452          if ((*_level)[u] == _level->maxLevel()) continue;
453          _flow->set(e, (*_capacity)[e]);
454          (*_excess)[u] += (*_capacity)[e];
455          if (u != _target && !_level->active(u)) {
456            _level->activate(u);
457          }
458        }
459      }
460    }
461
462    /// \brief Initializes the internal data structures using the
463    /// given flow map.
464    ///
465    /// Initializes the internal data structures and sets the initial
466    /// flow to the given \c flowMap. The \c flowMap should contain a
467    /// flow or at least a preflow, i.e. at each node excluding the
468    /// source node the incoming flow should greater or equal to the
469    /// outgoing flow.
470    /// \return \c false if the given \c flowMap is not a preflow.
471    template <typename FlowMap>
472    bool init(const FlowMap& flowMap) {
473      createStructures();
474
475      for (ArcIt e(_graph); e != INVALID; ++e) {
476        _flow->set(e, flowMap[e]);
477      }
478
479      for (NodeIt n(_graph); n != INVALID; ++n) {
480        Value excess = 0;
481        for (InArcIt e(_graph, n); e != INVALID; ++e) {
482          excess += (*_flow)[e];
483        }
484        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
485          excess -= (*_flow)[e];
486        }
487        if (excess < 0 && n != _source) return false;
488        (*_excess)[n] = excess;
489      }
490
491      typename Digraph::template NodeMap<bool> reached(_graph, false);
492
493      _level->initStart();
494      _level->initAddItem(_target);
495
496      std::vector<Node> queue;
497      reached[_source] = true;
498
499      queue.push_back(_target);
500      reached[_target] = true;
501      while (!queue.empty()) {
502        _level->initNewLevel();
503        std::vector<Node> nqueue;
504        for (int i = 0; i < int(queue.size()); ++i) {
505          Node n = queue[i];
506          for (InArcIt e(_graph, n); e != INVALID; ++e) {
507            Node u = _graph.source(e);
508            if (!reached[u] &&
509                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
510              reached[u] = true;
511              _level->initAddItem(u);
512              nqueue.push_back(u);
513            }
514          }
515          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
516            Node v = _graph.target(e);
517            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
518              reached[v] = true;
519              _level->initAddItem(v);
520              nqueue.push_back(v);
521            }
522          }
523        }
524        queue.swap(nqueue);
525      }
526      _level->initFinish();
527
528      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
529        Value rem = (*_capacity)[e] - (*_flow)[e];
530        if (_tolerance.positive(rem)) {
531          Node u = _graph.target(e);
532          if ((*_level)[u] == _level->maxLevel()) continue;
533          _flow->set(e, (*_capacity)[e]);
534          (*_excess)[u] += rem;
535          if (u != _target && !_level->active(u)) {
536            _level->activate(u);
537          }
538        }
539      }
540      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
541        Value rem = (*_flow)[e];
542        if (_tolerance.positive(rem)) {
543          Node v = _graph.source(e);
544          if ((*_level)[v] == _level->maxLevel()) continue;
545          _flow->set(e, 0);
546          (*_excess)[v] += rem;
547          if (v != _target && !_level->active(v)) {
548            _level->activate(v);
549          }
550        }
551      }
552      return true;
553    }
554
555    /// \brief Starts the first phase of the preflow algorithm.
556    ///
557    /// The preflow algorithm consists of two phases, this method runs
558    /// the first phase. After the first phase the maximum flow value
559    /// and a minimum value cut can already be computed, although a
560    /// maximum flow is not yet obtained. So after calling this method
561    /// \ref flowValue() returns the value of a maximum flow and \ref
562    /// minCut() returns a minimum cut.
563    /// \pre One of the \ref init() functions must be called before
564    /// using this function.
565    void startFirstPhase() {
566      _phase = true;
567
568      Node n = _level->highestActive();
569      int level = _level->highestActiveLevel();
570      while (n != INVALID) {
571        int num = _node_num;
572
573        while (num > 0 && n != INVALID) {
574          Value excess = (*_excess)[n];
575          int new_level = _level->maxLevel();
576
577          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
578            Value rem = (*_capacity)[e] - (*_flow)[e];
579            if (!_tolerance.positive(rem)) continue;
580            Node v = _graph.target(e);
581            if ((*_level)[v] < level) {
582              if (!_level->active(v) && v != _target) {
583                _level->activate(v);
584              }
585              if (!_tolerance.less(rem, excess)) {
586                _flow->set(e, (*_flow)[e] + excess);
587                (*_excess)[v] += excess;
588                excess = 0;
589                goto no_more_push_1;
590              } else {
591                excess -= rem;
592                (*_excess)[v] += rem;
593                _flow->set(e, (*_capacity)[e]);
594              }
595            } else if (new_level > (*_level)[v]) {
596              new_level = (*_level)[v];
597            }
598          }
599
600          for (InArcIt e(_graph, n); e != INVALID; ++e) {
601            Value rem = (*_flow)[e];
602            if (!_tolerance.positive(rem)) continue;
603            Node v = _graph.source(e);
604            if ((*_level)[v] < level) {
605              if (!_level->active(v) && v != _target) {
606                _level->activate(v);
607              }
608              if (!_tolerance.less(rem, excess)) {
609                _flow->set(e, (*_flow)[e] - excess);
610                (*_excess)[v] += excess;
611                excess = 0;
612                goto no_more_push_1;
613              } else {
614                excess -= rem;
615                (*_excess)[v] += rem;
616                _flow->set(e, 0);
617              }
618            } else if (new_level > (*_level)[v]) {
619              new_level = (*_level)[v];
620            }
621          }
622
623        no_more_push_1:
624
625          (*_excess)[n] = excess;
626
627          if (excess != 0) {
628            if (new_level + 1 < _level->maxLevel()) {
629              _level->liftHighestActive(new_level + 1);
630            } else {
631              _level->liftHighestActiveToTop();
632            }
633            if (_level->emptyLevel(level)) {
634              _level->liftToTop(level);
635            }
636          } else {
637            _level->deactivate(n);
638          }
639
640          n = _level->highestActive();
641          level = _level->highestActiveLevel();
642          --num;
643        }
644
645        num = _node_num * 20;
646        while (num > 0 && n != INVALID) {
647          Value excess = (*_excess)[n];
648          int new_level = _level->maxLevel();
649
650          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
651            Value rem = (*_capacity)[e] - (*_flow)[e];
652            if (!_tolerance.positive(rem)) continue;
653            Node v = _graph.target(e);
654            if ((*_level)[v] < level) {
655              if (!_level->active(v) && v != _target) {
656                _level->activate(v);
657              }
658              if (!_tolerance.less(rem, excess)) {
659                _flow->set(e, (*_flow)[e] + excess);
660                (*_excess)[v] += excess;
661                excess = 0;
662                goto no_more_push_2;
663              } else {
664                excess -= rem;
665                (*_excess)[v] += rem;
666                _flow->set(e, (*_capacity)[e]);
667              }
668            } else if (new_level > (*_level)[v]) {
669              new_level = (*_level)[v];
670            }
671          }
672
673          for (InArcIt e(_graph, n); e != INVALID; ++e) {
674            Value rem = (*_flow)[e];
675            if (!_tolerance.positive(rem)) continue;
676            Node v = _graph.source(e);
677            if ((*_level)[v] < level) {
678              if (!_level->active(v) && v != _target) {
679                _level->activate(v);
680              }
681              if (!_tolerance.less(rem, excess)) {
682                _flow->set(e, (*_flow)[e] - excess);
683                (*_excess)[v] += excess;
684                excess = 0;
685                goto no_more_push_2;
686              } else {
687                excess -= rem;
688                (*_excess)[v] += rem;
689                _flow->set(e, 0);
690              }
691            } else if (new_level > (*_level)[v]) {
692              new_level = (*_level)[v];
693            }
694          }
695
696        no_more_push_2:
697
698          (*_excess)[n] = excess;
699
700          if (excess != 0) {
701            if (new_level + 1 < _level->maxLevel()) {
702              _level->liftActiveOn(level, new_level + 1);
703            } else {
704              _level->liftActiveToTop(level);
705            }
706            if (_level->emptyLevel(level)) {
707              _level->liftToTop(level);
708            }
709          } else {
710            _level->deactivate(n);
711          }
712
713          while (level >= 0 && _level->activeFree(level)) {
714            --level;
715          }
716          if (level == -1) {
717            n = _level->highestActive();
718            level = _level->highestActiveLevel();
719          } else {
720            n = _level->activeOn(level);
721          }
722          --num;
723        }
724      }
725    }
726
727    /// \brief Starts the second phase of the preflow algorithm.
728    ///
729    /// The preflow algorithm consists of two phases, this method runs
730    /// the second phase. After calling one of the \ref init() functions
731    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
732    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
733    /// value of a maximum flow, \ref minCut() returns a minimum cut
734    /// \pre One of the \ref init() functions and \ref startFirstPhase()
735    /// must be called before using this function.
736    void startSecondPhase() {
737      _phase = false;
738
739      typename Digraph::template NodeMap<bool> reached(_graph);
740      for (NodeIt n(_graph); n != INVALID; ++n) {
741        reached[n] = (*_level)[n] < _level->maxLevel();
742      }
743
744      _level->initStart();
745      _level->initAddItem(_source);
746
747      std::vector<Node> queue;
748      queue.push_back(_source);
749      reached[_source] = true;
750
751      while (!queue.empty()) {
752        _level->initNewLevel();
753        std::vector<Node> nqueue;
754        for (int i = 0; i < int(queue.size()); ++i) {
755          Node n = queue[i];
756          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
757            Node v = _graph.target(e);
758            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
759              reached[v] = true;
760              _level->initAddItem(v);
761              nqueue.push_back(v);
762            }
763          }
764          for (InArcIt e(_graph, n); e != INVALID; ++e) {
765            Node u = _graph.source(e);
766            if (!reached[u] &&
767                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
768              reached[u] = true;
769              _level->initAddItem(u);
770              nqueue.push_back(u);
771            }
772          }
773        }
774        queue.swap(nqueue);
775      }
776      _level->initFinish();
777
778      for (NodeIt n(_graph); n != INVALID; ++n) {
779        if (!reached[n]) {
780          _level->dirtyTopButOne(n);
781        } else if ((*_excess)[n] > 0 && _target != n) {
782          _level->activate(n);
783        }
784      }
785
786      Node n;
787      while ((n = _level->highestActive()) != INVALID) {
788        Value excess = (*_excess)[n];
789        int level = _level->highestActiveLevel();
790        int new_level = _level->maxLevel();
791
792        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
793          Value rem = (*_capacity)[e] - (*_flow)[e];
794          if (!_tolerance.positive(rem)) continue;
795          Node v = _graph.target(e);
796          if ((*_level)[v] < level) {
797            if (!_level->active(v) && v != _source) {
798              _level->activate(v);
799            }
800            if (!_tolerance.less(rem, excess)) {
801              _flow->set(e, (*_flow)[e] + excess);
802              (*_excess)[v] += excess;
803              excess = 0;
804              goto no_more_push;
805            } else {
806              excess -= rem;
807              (*_excess)[v] += rem;
808              _flow->set(e, (*_capacity)[e]);
809            }
810          } else if (new_level > (*_level)[v]) {
811            new_level = (*_level)[v];
812          }
813        }
814
815        for (InArcIt e(_graph, n); e != INVALID; ++e) {
816          Value rem = (*_flow)[e];
817          if (!_tolerance.positive(rem)) continue;
818          Node v = _graph.source(e);
819          if ((*_level)[v] < level) {
820            if (!_level->active(v) && v != _source) {
821              _level->activate(v);
822            }
823            if (!_tolerance.less(rem, excess)) {
824              _flow->set(e, (*_flow)[e] - excess);
825              (*_excess)[v] += excess;
826              excess = 0;
827              goto no_more_push;
828            } else {
829              excess -= rem;
830              (*_excess)[v] += rem;
831              _flow->set(e, 0);
832            }
833          } else if (new_level > (*_level)[v]) {
834            new_level = (*_level)[v];
835          }
836        }
837
838      no_more_push:
839
840        (*_excess)[n] = excess;
841
842        if (excess != 0) {
843          if (new_level + 1 < _level->maxLevel()) {
844            _level->liftHighestActive(new_level + 1);
845          } else {
846            // Calculation error
847            _level->liftHighestActiveToTop();
848          }
849          if (_level->emptyLevel(level)) {
850            // Calculation error
851            _level->liftToTop(level);
852          }
853        } else {
854          _level->deactivate(n);
855        }
856
857      }
858    }
859
860    /// \brief Runs the preflow algorithm.
861    ///
862    /// Runs the preflow algorithm.
863    /// \note pf.run() is just a shortcut of the following code.
864    /// \code
865    ///   pf.init();
866    ///   pf.startFirstPhase();
867    ///   pf.startSecondPhase();
868    /// \endcode
869    void run() {
870      init();
871      startFirstPhase();
872      startSecondPhase();
873    }
874
875    /// \brief Runs the preflow algorithm to compute the minimum cut.
876    ///
877    /// Runs the preflow algorithm to compute the minimum cut.
878    /// \note pf.runMinCut() is just a shortcut of the following code.
879    /// \code
880    ///   pf.init();
881    ///   pf.startFirstPhase();
882    /// \endcode
883    void runMinCut() {
884      init();
885      startFirstPhase();
886    }
887
888    /// @}
889
890    /// \name Query Functions
891    /// The results of the preflow algorithm can be obtained using these
892    /// functions.\n
893    /// Either one of the \ref run() "run*()" functions or one of the
894    /// \ref startFirstPhase() "start*()" functions should be called
895    /// before using them.
896
897    ///@{
898
899    /// \brief Returns the value of the maximum flow.
900    ///
901    /// Returns the value of the maximum flow by returning the excess
902    /// of the target node. This value equals to the value of
903    /// the maximum flow already after the first phase of the algorithm.
904    ///
905    /// \pre Either \ref run() or \ref init() must be called before
906    /// using this function.
907    Value flowValue() const {
908      return (*_excess)[_target];
909    }
910
911    /// \brief Returns the flow value on the given arc.
912    ///
913    /// Returns the flow value on the given arc. This method can
914    /// be called after the second phase of the algorithm.
915    ///
916    /// \pre Either \ref run() or \ref init() must be called before
917    /// using this function.
918    Value flow(const Arc& arc) const {
919      return (*_flow)[arc];
920    }
921
922    /// \brief Returns a const reference to the flow map.
923    ///
924    /// Returns a const reference to the arc map storing the found flow.
925    /// This method can be called after the second phase of the algorithm.
926    ///
927    /// \pre Either \ref run() or \ref init() must be called before
928    /// using this function.
929    const FlowMap& flowMap() const {
930      return *_flow;
931    }
932
933    /// \brief Returns \c true when the node is on the source side of the
934    /// minimum cut.
935    ///
936    /// Returns true when the node is on the source side of the found
937    /// minimum cut. This method can be called both after running \ref
938    /// startFirstPhase() and \ref startSecondPhase().
939    ///
940    /// \pre Either \ref run() or \ref init() must be called before
941    /// using this function.
942    bool minCut(const Node& node) const {
943      return ((*_level)[node] == _level->maxLevel()) == _phase;
944    }
945
946    /// \brief Gives back a minimum value cut.
947    ///
948    /// Sets \c cutMap to the characteristic vector of a minimum value
949    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
950    /// node map with \c bool (or convertible) value type.
951    ///
952    /// This method can be called both after running \ref startFirstPhase()
953    /// and \ref startSecondPhase(). The result after the second phase
954    /// could be slightly different if inexact computation is used.
955    ///
956    /// \note This function calls \ref minCut() for each node, so it runs in
957    /// O(n) time.
958    ///
959    /// \pre Either \ref run() or \ref init() must be called before
960    /// using this function.
961    template <typename CutMap>
962    void minCutMap(CutMap& cutMap) const {
963      for (NodeIt n(_graph); n != INVALID; ++n) {
964        cutMap.set(n, minCut(n));
965      }
966    }
967
968    /// @}
969  };
970}
971
972#endif
Note: See TracBrowser for help on using the repository browser.