COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/preflow.h @ 871:86613aa28a0c

Last change on this file since 871:86613aa28a0c was 689:86c49553fea5, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Test file + doc improvements (#307)

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