COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/preflow.h @ 823:a7e93de12cbd

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

Add a warning about huge capacities in Preflow (#319)

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