COIN-OR::LEMON - Graph Library

source: lemon/lemon/preflow.h @ 761:98a30824fe36

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

Small doc improvements

File size: 29.7 KB
RevLine 
[404]1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
[463]5 * Copyright (C) 2003-2009
[404]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.
[525]34  /// \tparam GR Digraph type.
[606]35  /// \tparam CAP Capacity map type.
36  template <typename GR, typename CAP>
[404]37  struct PreflowDefaultTraits {
38
[408]39    /// \brief The type of the digraph the algorithm runs on.
[525]40    typedef GR Digraph;
[404]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.
[606]46    typedef CAP CapacityMap;
[404]47
[408]48    /// \brief The type of the flow values.
[688]49    typedef typename CapacityMap::Value Value;
[404]50
[408]51    /// \brief The type of the map that stores the flow values.
[404]52    ///
[408]53    /// The type of the map that stores the flow values.
[404]54    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
[760]55#ifdef DOXYGEN
56    typedef GR::ArcMap<Value> FlowMap;
57#else
[688]58    typedef typename Digraph::template ArcMap<Value> FlowMap;
[760]59#endif
[404]60
61    /// \brief Instantiates a FlowMap.
62    ///
63    /// This function instantiates a \ref FlowMap.
[657]64    /// \param digraph The digraph for which we would like to define
[404]65    /// the flow map.
66    static FlowMap* createFlowMap(const Digraph& digraph) {
67      return new FlowMap(digraph);
68    }
69
[408]70    /// \brief The elevator type used by Preflow algorithm.
[404]71    ///
72    /// The elevator type used by Preflow algorithm.
73    ///
[760]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
[404]80
81    /// \brief Instantiates an Elevator.
82    ///
[408]83    /// This function instantiates an \ref Elevator.
[657]84    /// \param digraph The digraph for which we would like to define
[404]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.
[688]94    typedef lemon::Tolerance<Value> Tolerance;
[404]95
96  };
97
98
99  /// \ingroup max_flow
100  ///
[408]101  /// \brief %Preflow algorithm class.
[404]102  ///
[408]103  /// This class provides an implementation of Goldberg-Tarjan's \e preflow
[606]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
[404]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  ///
[408]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.
[404]114  ///
[525]115  /// \tparam GR The type of the digraph the algorithm runs on.
[606]116  /// \tparam CAP The type of the capacity map. The default map
[525]117  /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
[404]118#ifdef DOXYGEN
[606]119  template <typename GR, typename CAP, typename TR>
[404]120#else
[525]121  template <typename GR,
[606]122            typename CAP = typename GR::template ArcMap<int>,
123            typename TR = PreflowDefaultTraits<GR, CAP> >
[404]124#endif
125  class Preflow {
126  public:
127
[408]128    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
[525]129    typedef TR Traits;
[408]130    ///The type of the digraph the algorithm runs on.
[404]131    typedef typename Traits::Digraph Digraph;
[408]132    ///The type of the capacity map.
[404]133    typedef typename Traits::CapacityMap CapacityMap;
[408]134    ///The type of the flow values.
[688]135    typedef typename Traits::Value Value;
[404]136
[408]137    ///The type of the flow map.
[404]138    typedef typename Traits::FlowMap FlowMap;
[408]139    ///The type of the elevator.
[404]140    typedef typename Traits::Elevator Elevator;
[408]141    ///The type of the tolerance.
[404]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
[688]161    typedef typename Digraph::template NodeMap<Value> ExcessMap;
[404]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
[408]201    ///\name Named Template Parameters
[404]202
203    ///@{
204
[606]205    template <typename T>
[406]206    struct SetFlowMapTraits : public Traits {
[606]207      typedef T FlowMap;
[404]208      static FlowMap *createFlowMap(const Digraph&) {
[405]209        LEMON_ASSERT(false, "FlowMap is not initialized");
210        return 0; // ignore warnings
[404]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
[408]218    /// type.
[606]219    template <typename T>
[406]220    struct SetFlowMap
[606]221      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
[404]222      typedef Preflow<Digraph, CapacityMap,
[606]223                      SetFlowMapTraits<T> > Create;
[404]224    };
225
[606]226    template <typename T>
[406]227    struct SetElevatorTraits : public Traits {
[606]228      typedef T Elevator;
[404]229      static Elevator *createElevator(const Digraph&, int) {
[405]230        LEMON_ASSERT(false, "Elevator is not initialized");
231        return 0; // ignore warnings
[404]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
[408]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
[606]244    template <typename T>
[406]245    struct SetElevator
[606]246      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
[404]247      typedef Preflow<Digraph, CapacityMap,
[606]248                      SetElevatorTraits<T> > Create;
[404]249    };
250
[606]251    template <typename T>
[406]252    struct SetStandardElevatorTraits : public Traits {
[606]253      typedef T Elevator;
[404]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
[408]260    /// Elevator type with automatic allocation
[404]261    ///
262    /// \ref named-templ-param "Named parameter" for setting Elevator
[408]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
[606]271    template <typename T>
[406]272    struct SetStandardElevator
[404]273      : public Preflow<Digraph, CapacityMap,
[606]274                       SetStandardElevatorTraits<T> > {
[404]275      typedef Preflow<Digraph, CapacityMap,
[606]276                      SetStandardElevatorTraits<T> > Create;
[404]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,
[408]296            Node source, Node target)
[404]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
[408]303    /// \brief Destructor.
[404]304    ///
305    /// Destructor.
306    ~Preflow() {
307      destroyStructures();
308    }
309
310    /// \brief Sets the capacity map.
311    ///
312    /// Sets the capacity map.
[408]313    /// \return <tt>(*this)</tt>
[404]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.
[408]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>
[404]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
[408]336    /// \brief Sets the source node.
[404]337    ///
[408]338    /// Sets the source node.
339    /// \return <tt>(*this)</tt>
340    Preflow& source(const Node& node) {
341      _source = node;
342      return *this;
[404]343    }
344
[408]345    /// \brief Sets the target node.
[404]346    ///
[408]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>
[404]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
[408]371    /// \brief Returns a const reference to the elevator.
[404]372    ///
[408]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.
[437]377    const Elevator& elevator() const {
[404]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
[408]389    /// \brief Returns a const reference to the tolerance.
[404]390    ///
[408]391    /// Returns a const reference to the tolerance.
[404]392    const Tolerance& tolerance() const {
393      return tolerance;
394    }
395
[408]396    /// \name Execution Control
397    /// The simplest way to execute the preflow algorithm is to use
398    /// \ref run() or \ref runMinCut().\n
[760]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
[408]401    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
[404]402
403    ///@{
404
405    /// \brief Initializes the internal data structures.
406    ///
[408]407    /// Initializes the internal data structures and sets the initial
408    /// flow to zero on each arc.
[404]409    void init() {
410      createStructures();
411
412      _phase = true;
413      for (NodeIt n(_graph); n != INVALID; ++n) {
[628]414        (*_excess)[n] = 0;
[404]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;
[628]427      reached[_source] = true;
[404]428
429      queue.push_back(_target);
[628]430      reached[_target] = true;
[404]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])) {
[628]439              reached[u] = true;
[404]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]);
[628]454          (*_excess)[u] += (*_capacity)[e];
[404]455          if (u != _target && !_level->active(u)) {
456            _level->activate(u);
457          }
458        }
459      }
460    }
461
[408]462    /// \brief Initializes the internal data structures using the
463    /// given flow map.
[404]464    ///
465    /// Initializes the internal data structures and sets the initial
466    /// flow to the given \c flowMap. The \c flowMap should contain a
[408]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
[404]469    /// outgoing flow.
[408]470    /// \return \c false if the given \c flowMap is not a preflow.
[404]471    template <typename FlowMap>
[407]472    bool init(const FlowMap& flowMap) {
[404]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) {
[688]480        Value excess = 0;
[404]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;
[628]488        (*_excess)[n] = excess;
[404]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;
[628]497      reached[_source] = true;
[404]498
499      queue.push_back(_target);
[628]500      reached[_target] = true;
[404]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])) {
[628]510              reached[u] = true;
[404]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])) {
[628]518              reached[v] = true;
[404]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) {
[688]529        Value rem = (*_capacity)[e] - (*_flow)[e];
[404]530        if (_tolerance.positive(rem)) {
531          Node u = _graph.target(e);
532          if ((*_level)[u] == _level->maxLevel()) continue;
533          _flow->set(e, (*_capacity)[e]);
[628]534          (*_excess)[u] += rem;
[404]535          if (u != _target && !_level->active(u)) {
536            _level->activate(u);
537          }
538        }
539      }
540      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
[688]541        Value rem = (*_flow)[e];
[404]542        if (_tolerance.positive(rem)) {
543          Node v = _graph.source(e);
544          if ((*_level)[v] == _level->maxLevel()) continue;
545          _flow->set(e, 0);
[628]546          (*_excess)[v] += rem;
[404]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.
[408]563    /// \pre One of the \ref init() functions must be called before
564    /// using this function.
[404]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) {
[688]574          Value excess = (*_excess)[n];
[404]575          int new_level = _level->maxLevel();
576
577          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
[688]578            Value rem = (*_capacity)[e] - (*_flow)[e];
[404]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);
[628]587                (*_excess)[v] += excess;
[404]588                excess = 0;
589                goto no_more_push_1;
590              } else {
591                excess -= rem;
[628]592                (*_excess)[v] += rem;
[404]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) {
[688]601            Value rem = (*_flow)[e];
[404]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);
[628]610                (*_excess)[v] += excess;
[404]611                excess = 0;
612                goto no_more_push_1;
613              } else {
614                excess -= rem;
[628]615                (*_excess)[v] += rem;
[404]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
[628]625          (*_excess)[n] = excess;
[404]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) {
[688]647          Value excess = (*_excess)[n];
[404]648          int new_level = _level->maxLevel();
649
650          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
[688]651            Value rem = (*_capacity)[e] - (*_flow)[e];
[404]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);
[628]660                (*_excess)[v] += excess;
[404]661                excess = 0;
662                goto no_more_push_2;
663              } else {
664                excess -= rem;
[628]665                (*_excess)[v] += rem;
[404]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) {
[688]674            Value rem = (*_flow)[e];
[404]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);
[628]683                (*_excess)[v] += excess;
[404]684                excess = 0;
685                goto no_more_push_2;
686              } else {
687                excess -= rem;
[628]688                (*_excess)[v] += rem;
[404]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
[628]698          (*_excess)[n] = excess;
[404]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
[408]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
[404]733    /// value of a maximum flow, \ref minCut() returns a minimum cut
[408]734    /// \pre One of the \ref init() functions and \ref startFirstPhase()
735    /// must be called before using this function.
[404]736    void startSecondPhase() {
737      _phase = false;
738
739      typename Digraph::template NodeMap<bool> reached(_graph);
740      for (NodeIt n(_graph); n != INVALID; ++n) {
[628]741        reached[n] = (*_level)[n] < _level->maxLevel();
[404]742      }
743
744      _level->initStart();
745      _level->initAddItem(_source);
746
747      std::vector<Node> queue;
748      queue.push_back(_source);
[628]749      reached[_source] = true;
[404]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])) {
[628]759              reached[v] = true;
[404]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])) {
[628]768              reached[u] = true;
[404]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) {
[688]788        Value excess = (*_excess)[n];
[404]789        int level = _level->highestActiveLevel();
790        int new_level = _level->maxLevel();
791
792        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
[688]793          Value rem = (*_capacity)[e] - (*_flow)[e];
[404]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);
[628]802              (*_excess)[v] += excess;
[404]803              excess = 0;
804              goto no_more_push;
805            } else {
806              excess -= rem;
[628]807              (*_excess)[v] += rem;
[404]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) {
[688]816          Value rem = (*_flow)[e];
[404]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);
[628]825              (*_excess)[v] += excess;
[404]826              excess = 0;
827              goto no_more_push;
828            } else {
829              excess -= rem;
[628]830              (*_excess)[v] += rem;
[404]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
[628]840        (*_excess)[n] = excess;
[404]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
[408]891    /// The results of the preflow algorithm can be obtained using these
[404]892    /// functions.\n
[408]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.
[404]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
[408]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.
[688]907    Value flowValue() const {
[404]908      return (*_excess)[_target];
909    }
910
[688]911    /// \brief Returns the flow value on the given arc.
[404]912    ///
[688]913    /// Returns the flow value on the given arc. This method can
[408]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.
[688]918    Value flow(const Arc& arc) const {
[408]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.
[437]929    const FlowMap& flowMap() const {
[408]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
[404]938    /// startFirstPhase() and \ref startSecondPhase().
[408]939    ///
940    /// \pre Either \ref run() or \ref init() must be called before
941    /// using this function.
[404]942    bool minCut(const Node& node) const {
943      return ((*_level)[node] == _level->maxLevel()) == _phase;
944    }
945
[408]946    /// \brief Gives back a minimum value cut.
[404]947    ///
[408]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
[606]957    /// O(n) time.
[408]958    ///
959    /// \pre Either \ref run() or \ref init() must be called before
960    /// using this function.
[404]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.