COIN-OR::LEMON - Graph Library

source: lemon/lemon/preflow.h @ 833:e20173729589

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

Small doc fixes in several files (#331)

File size: 29.8 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
[736]107  /// flow algorithms. The current implementation uses a mixture of the
[404]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).
[833]267    /// However, an external elevator object could also be passed to the
[408]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
[736]381    /// \brief Sets the tolerance used by the algorithm.
[404]382    ///
[736]383    /// Sets the tolerance object used by the algorithm.
384    /// \return <tt>(*this)</tt>
[735]385    Preflow& tolerance(const Tolerance& tolerance) {
[404]386      _tolerance = tolerance;
387      return *this;
388    }
389
[408]390    /// \brief Returns a const reference to the tolerance.
[404]391    ///
[736]392    /// Returns a const reference to the tolerance object used by
393    /// the algorithm.
[404]394    const Tolerance& tolerance() const {
[735]395      return _tolerance;
[404]396    }
397
[408]398    /// \name Execution Control
399    /// The simplest way to execute the preflow algorithm is to use
400    /// \ref run() or \ref runMinCut().\n
[760]401    /// If you need better control on the initial solution or the execution,
402    /// you have to call one of the \ref init() functions first, then
[408]403    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
[404]404
405    ///@{
406
407    /// \brief Initializes the internal data structures.
408    ///
[408]409    /// Initializes the internal data structures and sets the initial
410    /// flow to zero on each arc.
[404]411    void init() {
412      createStructures();
413
414      _phase = true;
415      for (NodeIt n(_graph); n != INVALID; ++n) {
[628]416        (*_excess)[n] = 0;
[404]417      }
418
419      for (ArcIt e(_graph); e != INVALID; ++e) {
420        _flow->set(e, 0);
421      }
422
423      typename Digraph::template NodeMap<bool> reached(_graph, false);
424
425      _level->initStart();
426      _level->initAddItem(_target);
427
428      std::vector<Node> queue;
[628]429      reached[_source] = true;
[404]430
431      queue.push_back(_target);
[628]432      reached[_target] = true;
[404]433      while (!queue.empty()) {
434        _level->initNewLevel();
435        std::vector<Node> nqueue;
436        for (int i = 0; i < int(queue.size()); ++i) {
437          Node n = queue[i];
438          for (InArcIt e(_graph, n); e != INVALID; ++e) {
439            Node u = _graph.source(e);
440            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
[628]441              reached[u] = true;
[404]442              _level->initAddItem(u);
443              nqueue.push_back(u);
444            }
445          }
446        }
447        queue.swap(nqueue);
448      }
449      _level->initFinish();
450
451      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
452        if (_tolerance.positive((*_capacity)[e])) {
453          Node u = _graph.target(e);
454          if ((*_level)[u] == _level->maxLevel()) continue;
455          _flow->set(e, (*_capacity)[e]);
[628]456          (*_excess)[u] += (*_capacity)[e];
[404]457          if (u != _target && !_level->active(u)) {
458            _level->activate(u);
459          }
460        }
461      }
462    }
463
[408]464    /// \brief Initializes the internal data structures using the
465    /// given flow map.
[404]466    ///
467    /// Initializes the internal data structures and sets the initial
468    /// flow to the given \c flowMap. The \c flowMap should contain a
[408]469    /// flow or at least a preflow, i.e. at each node excluding the
470    /// source node the incoming flow should greater or equal to the
[404]471    /// outgoing flow.
[408]472    /// \return \c false if the given \c flowMap is not a preflow.
[404]473    template <typename FlowMap>
[407]474    bool init(const FlowMap& flowMap) {
[404]475      createStructures();
476
477      for (ArcIt e(_graph); e != INVALID; ++e) {
478        _flow->set(e, flowMap[e]);
479      }
480
481      for (NodeIt n(_graph); n != INVALID; ++n) {
[688]482        Value excess = 0;
[404]483        for (InArcIt e(_graph, n); e != INVALID; ++e) {
484          excess += (*_flow)[e];
485        }
486        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
487          excess -= (*_flow)[e];
488        }
489        if (excess < 0 && n != _source) return false;
[628]490        (*_excess)[n] = excess;
[404]491      }
492
493      typename Digraph::template NodeMap<bool> reached(_graph, false);
494
495      _level->initStart();
496      _level->initAddItem(_target);
497
498      std::vector<Node> queue;
[628]499      reached[_source] = true;
[404]500
501      queue.push_back(_target);
[628]502      reached[_target] = true;
[404]503      while (!queue.empty()) {
504        _level->initNewLevel();
505        std::vector<Node> nqueue;
506        for (int i = 0; i < int(queue.size()); ++i) {
507          Node n = queue[i];
508          for (InArcIt e(_graph, n); e != INVALID; ++e) {
509            Node u = _graph.source(e);
510            if (!reached[u] &&
511                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
[628]512              reached[u] = true;
[404]513              _level->initAddItem(u);
514              nqueue.push_back(u);
515            }
516          }
517          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
518            Node v = _graph.target(e);
519            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
[628]520              reached[v] = true;
[404]521              _level->initAddItem(v);
522              nqueue.push_back(v);
523            }
524          }
525        }
526        queue.swap(nqueue);
527      }
528      _level->initFinish();
529
530      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
[688]531        Value rem = (*_capacity)[e] - (*_flow)[e];
[404]532        if (_tolerance.positive(rem)) {
533          Node u = _graph.target(e);
534          if ((*_level)[u] == _level->maxLevel()) continue;
535          _flow->set(e, (*_capacity)[e]);
[628]536          (*_excess)[u] += rem;
[404]537          if (u != _target && !_level->active(u)) {
538            _level->activate(u);
539          }
540        }
541      }
542      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
[688]543        Value rem = (*_flow)[e];
[404]544        if (_tolerance.positive(rem)) {
545          Node v = _graph.source(e);
546          if ((*_level)[v] == _level->maxLevel()) continue;
547          _flow->set(e, 0);
[628]548          (*_excess)[v] += rem;
[404]549          if (v != _target && !_level->active(v)) {
550            _level->activate(v);
551          }
552        }
553      }
554      return true;
555    }
556
557    /// \brief Starts the first phase of the preflow algorithm.
558    ///
559    /// The preflow algorithm consists of two phases, this method runs
560    /// the first phase. After the first phase the maximum flow value
561    /// and a minimum value cut can already be computed, although a
562    /// maximum flow is not yet obtained. So after calling this method
563    /// \ref flowValue() returns the value of a maximum flow and \ref
564    /// minCut() returns a minimum cut.
[408]565    /// \pre One of the \ref init() functions must be called before
566    /// using this function.
[404]567    void startFirstPhase() {
568      _phase = true;
569
570      Node n = _level->highestActive();
571      int level = _level->highestActiveLevel();
572      while (n != INVALID) {
573        int num = _node_num;
574
575        while (num > 0 && n != INVALID) {
[688]576          Value excess = (*_excess)[n];
[404]577          int new_level = _level->maxLevel();
578
579          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
[688]580            Value rem = (*_capacity)[e] - (*_flow)[e];
[404]581            if (!_tolerance.positive(rem)) continue;
582            Node v = _graph.target(e);
583            if ((*_level)[v] < level) {
584              if (!_level->active(v) && v != _target) {
585                _level->activate(v);
586              }
587              if (!_tolerance.less(rem, excess)) {
588                _flow->set(e, (*_flow)[e] + excess);
[628]589                (*_excess)[v] += excess;
[404]590                excess = 0;
591                goto no_more_push_1;
592              } else {
593                excess -= rem;
[628]594                (*_excess)[v] += rem;
[404]595                _flow->set(e, (*_capacity)[e]);
596              }
597            } else if (new_level > (*_level)[v]) {
598              new_level = (*_level)[v];
599            }
600          }
601
602          for (InArcIt e(_graph, n); e != INVALID; ++e) {
[688]603            Value rem = (*_flow)[e];
[404]604            if (!_tolerance.positive(rem)) continue;
605            Node v = _graph.source(e);
606            if ((*_level)[v] < level) {
607              if (!_level->active(v) && v != _target) {
608                _level->activate(v);
609              }
610              if (!_tolerance.less(rem, excess)) {
611                _flow->set(e, (*_flow)[e] - excess);
[628]612                (*_excess)[v] += excess;
[404]613                excess = 0;
614                goto no_more_push_1;
615              } else {
616                excess -= rem;
[628]617                (*_excess)[v] += rem;
[404]618                _flow->set(e, 0);
619              }
620            } else if (new_level > (*_level)[v]) {
621              new_level = (*_level)[v];
622            }
623          }
624
625        no_more_push_1:
626
[628]627          (*_excess)[n] = excess;
[404]628
629          if (excess != 0) {
630            if (new_level + 1 < _level->maxLevel()) {
631              _level->liftHighestActive(new_level + 1);
632            } else {
633              _level->liftHighestActiveToTop();
634            }
635            if (_level->emptyLevel(level)) {
636              _level->liftToTop(level);
637            }
638          } else {
639            _level->deactivate(n);
640          }
641
642          n = _level->highestActive();
643          level = _level->highestActiveLevel();
644          --num;
645        }
646
647        num = _node_num * 20;
648        while (num > 0 && n != INVALID) {
[688]649          Value excess = (*_excess)[n];
[404]650          int new_level = _level->maxLevel();
651
652          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
[688]653            Value rem = (*_capacity)[e] - (*_flow)[e];
[404]654            if (!_tolerance.positive(rem)) continue;
655            Node v = _graph.target(e);
656            if ((*_level)[v] < level) {
657              if (!_level->active(v) && v != _target) {
658                _level->activate(v);
659              }
660              if (!_tolerance.less(rem, excess)) {
661                _flow->set(e, (*_flow)[e] + excess);
[628]662                (*_excess)[v] += excess;
[404]663                excess = 0;
664                goto no_more_push_2;
665              } else {
666                excess -= rem;
[628]667                (*_excess)[v] += rem;
[404]668                _flow->set(e, (*_capacity)[e]);
669              }
670            } else if (new_level > (*_level)[v]) {
671              new_level = (*_level)[v];
672            }
673          }
674
675          for (InArcIt e(_graph, n); e != INVALID; ++e) {
[688]676            Value rem = (*_flow)[e];
[404]677            if (!_tolerance.positive(rem)) continue;
678            Node v = _graph.source(e);
679            if ((*_level)[v] < level) {
680              if (!_level->active(v) && v != _target) {
681                _level->activate(v);
682              }
683              if (!_tolerance.less(rem, excess)) {
684                _flow->set(e, (*_flow)[e] - excess);
[628]685                (*_excess)[v] += excess;
[404]686                excess = 0;
687                goto no_more_push_2;
688              } else {
689                excess -= rem;
[628]690                (*_excess)[v] += rem;
[404]691                _flow->set(e, 0);
692              }
693            } else if (new_level > (*_level)[v]) {
694              new_level = (*_level)[v];
695            }
696          }
697
698        no_more_push_2:
699
[628]700          (*_excess)[n] = excess;
[404]701
702          if (excess != 0) {
703            if (new_level + 1 < _level->maxLevel()) {
704              _level->liftActiveOn(level, new_level + 1);
705            } else {
706              _level->liftActiveToTop(level);
707            }
708            if (_level->emptyLevel(level)) {
709              _level->liftToTop(level);
710            }
711          } else {
712            _level->deactivate(n);
713          }
714
715          while (level >= 0 && _level->activeFree(level)) {
716            --level;
717          }
718          if (level == -1) {
719            n = _level->highestActive();
720            level = _level->highestActiveLevel();
721          } else {
722            n = _level->activeOn(level);
723          }
724          --num;
725        }
726      }
727    }
728
729    /// \brief Starts the second phase of the preflow algorithm.
730    ///
731    /// The preflow algorithm consists of two phases, this method runs
[408]732    /// the second phase. After calling one of the \ref init() functions
733    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
734    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
[404]735    /// value of a maximum flow, \ref minCut() returns a minimum cut
[408]736    /// \pre One of the \ref init() functions and \ref startFirstPhase()
737    /// must be called before using this function.
[404]738    void startSecondPhase() {
739      _phase = false;
740
741      typename Digraph::template NodeMap<bool> reached(_graph);
742      for (NodeIt n(_graph); n != INVALID; ++n) {
[628]743        reached[n] = (*_level)[n] < _level->maxLevel();
[404]744      }
745
746      _level->initStart();
747      _level->initAddItem(_source);
748
749      std::vector<Node> queue;
750      queue.push_back(_source);
[628]751      reached[_source] = true;
[404]752
753      while (!queue.empty()) {
754        _level->initNewLevel();
755        std::vector<Node> nqueue;
756        for (int i = 0; i < int(queue.size()); ++i) {
757          Node n = queue[i];
758          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
759            Node v = _graph.target(e);
760            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
[628]761              reached[v] = true;
[404]762              _level->initAddItem(v);
763              nqueue.push_back(v);
764            }
765          }
766          for (InArcIt e(_graph, n); e != INVALID; ++e) {
767            Node u = _graph.source(e);
768            if (!reached[u] &&
769                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
[628]770              reached[u] = true;
[404]771              _level->initAddItem(u);
772              nqueue.push_back(u);
773            }
774          }
775        }
776        queue.swap(nqueue);
777      }
778      _level->initFinish();
779
780      for (NodeIt n(_graph); n != INVALID; ++n) {
781        if (!reached[n]) {
782          _level->dirtyTopButOne(n);
783        } else if ((*_excess)[n] > 0 && _target != n) {
784          _level->activate(n);
785        }
786      }
787
788      Node n;
789      while ((n = _level->highestActive()) != INVALID) {
[688]790        Value excess = (*_excess)[n];
[404]791        int level = _level->highestActiveLevel();
792        int new_level = _level->maxLevel();
793
794        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
[688]795          Value rem = (*_capacity)[e] - (*_flow)[e];
[404]796          if (!_tolerance.positive(rem)) continue;
797          Node v = _graph.target(e);
798          if ((*_level)[v] < level) {
799            if (!_level->active(v) && v != _source) {
800              _level->activate(v);
801            }
802            if (!_tolerance.less(rem, excess)) {
803              _flow->set(e, (*_flow)[e] + excess);
[628]804              (*_excess)[v] += excess;
[404]805              excess = 0;
806              goto no_more_push;
807            } else {
808              excess -= rem;
[628]809              (*_excess)[v] += rem;
[404]810              _flow->set(e, (*_capacity)[e]);
811            }
812          } else if (new_level > (*_level)[v]) {
813            new_level = (*_level)[v];
814          }
815        }
816
817        for (InArcIt e(_graph, n); e != INVALID; ++e) {
[688]818          Value rem = (*_flow)[e];
[404]819          if (!_tolerance.positive(rem)) continue;
820          Node v = _graph.source(e);
821          if ((*_level)[v] < level) {
822            if (!_level->active(v) && v != _source) {
823              _level->activate(v);
824            }
825            if (!_tolerance.less(rem, excess)) {
826              _flow->set(e, (*_flow)[e] - excess);
[628]827              (*_excess)[v] += excess;
[404]828              excess = 0;
829              goto no_more_push;
830            } else {
831              excess -= rem;
[628]832              (*_excess)[v] += rem;
[404]833              _flow->set(e, 0);
834            }
835          } else if (new_level > (*_level)[v]) {
836            new_level = (*_level)[v];
837          }
838        }
839
840      no_more_push:
841
[628]842        (*_excess)[n] = excess;
[404]843
844        if (excess != 0) {
845          if (new_level + 1 < _level->maxLevel()) {
846            _level->liftHighestActive(new_level + 1);
847          } else {
848            // Calculation error
849            _level->liftHighestActiveToTop();
850          }
851          if (_level->emptyLevel(level)) {
852            // Calculation error
853            _level->liftToTop(level);
854          }
855        } else {
856          _level->deactivate(n);
857        }
858
859      }
860    }
861
862    /// \brief Runs the preflow algorithm.
863    ///
864    /// Runs the preflow algorithm.
865    /// \note pf.run() is just a shortcut of the following code.
866    /// \code
867    ///   pf.init();
868    ///   pf.startFirstPhase();
869    ///   pf.startSecondPhase();
870    /// \endcode
871    void run() {
872      init();
873      startFirstPhase();
874      startSecondPhase();
875    }
876
877    /// \brief Runs the preflow algorithm to compute the minimum cut.
878    ///
879    /// Runs the preflow algorithm to compute the minimum cut.
880    /// \note pf.runMinCut() is just a shortcut of the following code.
881    /// \code
882    ///   pf.init();
883    ///   pf.startFirstPhase();
884    /// \endcode
885    void runMinCut() {
886      init();
887      startFirstPhase();
888    }
889
890    /// @}
891
892    /// \name Query Functions
[408]893    /// The results of the preflow algorithm can be obtained using these
[404]894    /// functions.\n
[408]895    /// Either one of the \ref run() "run*()" functions or one of the
896    /// \ref startFirstPhase() "start*()" functions should be called
897    /// before using them.
[404]898
899    ///@{
900
901    /// \brief Returns the value of the maximum flow.
902    ///
903    /// Returns the value of the maximum flow by returning the excess
[408]904    /// of the target node. This value equals to the value of
905    /// the maximum flow already after the first phase of the algorithm.
906    ///
907    /// \pre Either \ref run() or \ref init() must be called before
908    /// using this function.
[688]909    Value flowValue() const {
[404]910      return (*_excess)[_target];
911    }
912
[688]913    /// \brief Returns the flow value on the given arc.
[404]914    ///
[688]915    /// Returns the flow value on the given arc. This method can
[408]916    /// be called after the second phase of the algorithm.
917    ///
918    /// \pre Either \ref run() or \ref init() must be called before
919    /// using this function.
[688]920    Value flow(const Arc& arc) const {
[408]921      return (*_flow)[arc];
922    }
923
924    /// \brief Returns a const reference to the flow map.
925    ///
926    /// Returns a const reference to the arc map storing the found flow.
927    /// This method can be called after the second phase of the algorithm.
928    ///
929    /// \pre Either \ref run() or \ref init() must be called before
930    /// using this function.
[437]931    const FlowMap& flowMap() const {
[408]932      return *_flow;
933    }
934
935    /// \brief Returns \c true when the node is on the source side of the
936    /// minimum cut.
937    ///
938    /// Returns true when the node is on the source side of the found
939    /// minimum cut. This method can be called both after running \ref
[404]940    /// startFirstPhase() and \ref startSecondPhase().
[408]941    ///
942    /// \pre Either \ref run() or \ref init() must be called before
943    /// using this function.
[404]944    bool minCut(const Node& node) const {
945      return ((*_level)[node] == _level->maxLevel()) == _phase;
946    }
947
[408]948    /// \brief Gives back a minimum value cut.
[404]949    ///
[408]950    /// Sets \c cutMap to the characteristic vector of a minimum value
951    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
952    /// node map with \c bool (or convertible) value type.
953    ///
954    /// This method can be called both after running \ref startFirstPhase()
955    /// and \ref startSecondPhase(). The result after the second phase
956    /// could be slightly different if inexact computation is used.
957    ///
958    /// \note This function calls \ref minCut() for each node, so it runs in
[606]959    /// O(n) time.
[408]960    ///
961    /// \pre Either \ref run() or \ref init() must be called before
962    /// using this function.
[404]963    template <typename CutMap>
964    void minCutMap(CutMap& cutMap) const {
965      for (NodeIt n(_graph); n != INVALID; ++n) {
966        cutMap.set(n, minCut(n));
967      }
968    }
969
970    /// @}
971  };
972}
973
974#endif
Note: See TracBrowser for help on using the repository browser.