COIN-OR::LEMON - Graph Library

source: lemon/lemon/preflow.h @ 982:bb70ad62c95f

Last change on this file since 982:bb70ad62c95f was 982:bb70ad62c95f, checked in by Balazs Dezso <deba@…>, 14 years ago

Fix critical bug in preflow (#372)

The wrong transition between the bound decrease and highest active
heuristics caused the bug. The last node chosen in bound decrease mode
is used in the first iteration in highest active mode.

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