COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/preflow.h @ 901:30d5f950aa5f

Last change on this file since 901:30d5f950aa5f was 901:30d5f950aa5f, checked in by Alpar Juttner <alpar@…>, 9 years ago

Fix wrong initialization in Preflow (#414)

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        }
529      }
530      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
531        Value rem = (*_flow)[e];
532        if (_tolerance.positive(rem)) {
533          Node v = _graph.source(e);
534          if ((*_level)[v] == _level->maxLevel()) continue;
535          _flow->set(e, 0);
536          (*_excess)[v] += rem;
537        }
538      }
539      for (NodeIt n(_graph); n != INVALID; ++n)
540        if(n!=_source && n!=_target && _tolerance.positive((*_excess)[n]))
541          _level->activate(n);
542         
543      return true;
544    }
545
546    /// \brief Starts the first phase of the preflow algorithm.
547    ///
548    /// The preflow algorithm consists of two phases, this method runs
549    /// the first phase. After the first phase the maximum flow value
550    /// and a minimum value cut can already be computed, although a
551    /// maximum flow is not yet obtained. So after calling this method
552    /// \ref flowValue() returns the value of a maximum flow and \ref
553    /// minCut() returns a minimum cut.
554    /// \pre One of the \ref init() functions must be called before
555    /// using this function.
556    void startFirstPhase() {
557      _phase = true;
558
559      while (true) {
560        int num = _node_num;
561
562        Node n = INVALID;
563        int level = -1;
564
565        while (num > 0) {
566          n = _level->highestActive();
567          if (n == INVALID) goto first_phase_done;
568          level = _level->highestActiveLevel();
569          --num;
570         
571          Value excess = (*_excess)[n];
572          int new_level = _level->maxLevel();
573
574          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
575            Value rem = (*_capacity)[e] - (*_flow)[e];
576            if (!_tolerance.positive(rem)) continue;
577            Node v = _graph.target(e);
578            if ((*_level)[v] < level) {
579              if (!_level->active(v) && v != _target) {
580                _level->activate(v);
581              }
582              if (!_tolerance.less(rem, excess)) {
583                _flow->set(e, (*_flow)[e] + excess);
584                (*_excess)[v] += excess;
585                excess = 0;
586                goto no_more_push_1;
587              } else {
588                excess -= rem;
589                (*_excess)[v] += rem;
590                _flow->set(e, (*_capacity)[e]);
591              }
592            } else if (new_level > (*_level)[v]) {
593              new_level = (*_level)[v];
594            }
595          }
596
597          for (InArcIt e(_graph, n); e != INVALID; ++e) {
598            Value rem = (*_flow)[e];
599            if (!_tolerance.positive(rem)) continue;
600            Node v = _graph.source(e);
601            if ((*_level)[v] < level) {
602              if (!_level->active(v) && v != _target) {
603                _level->activate(v);
604              }
605              if (!_tolerance.less(rem, excess)) {
606                _flow->set(e, (*_flow)[e] - excess);
607                (*_excess)[v] += excess;
608                excess = 0;
609                goto no_more_push_1;
610              } else {
611                excess -= rem;
612                (*_excess)[v] += rem;
613                _flow->set(e, 0);
614              }
615            } else if (new_level > (*_level)[v]) {
616              new_level = (*_level)[v];
617            }
618          }
619
620        no_more_push_1:
621
622          (*_excess)[n] = excess;
623
624          if (excess != 0) {
625            if (new_level + 1 < _level->maxLevel()) {
626              _level->liftHighestActive(new_level + 1);
627            } else {
628              _level->liftHighestActiveToTop();
629            }
630            if (_level->emptyLevel(level)) {
631              _level->liftToTop(level);
632            }
633          } else {
634            _level->deactivate(n);
635          }
636        }
637
638        num = _node_num * 20;
639        while (num > 0) {
640          while (level >= 0 && _level->activeFree(level)) {
641            --level;
642          }
643          if (level == -1) {
644            n = _level->highestActive();
645            level = _level->highestActiveLevel();
646            if (n == INVALID) goto first_phase_done;
647          } else {
648            n = _level->activeOn(level);
649          }
650          --num;
651
652          Value excess = (*_excess)[n];
653          int new_level = _level->maxLevel();
654
655          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
656            Value rem = (*_capacity)[e] - (*_flow)[e];
657            if (!_tolerance.positive(rem)) continue;
658            Node v = _graph.target(e);
659            if ((*_level)[v] < level) {
660              if (!_level->active(v) && v != _target) {
661                _level->activate(v);
662              }
663              if (!_tolerance.less(rem, excess)) {
664                _flow->set(e, (*_flow)[e] + excess);
665                (*_excess)[v] += excess;
666                excess = 0;
667                goto no_more_push_2;
668              } else {
669                excess -= rem;
670                (*_excess)[v] += rem;
671                _flow->set(e, (*_capacity)[e]);
672              }
673            } else if (new_level > (*_level)[v]) {
674              new_level = (*_level)[v];
675            }
676          }
677
678          for (InArcIt e(_graph, n); e != INVALID; ++e) {
679            Value rem = (*_flow)[e];
680            if (!_tolerance.positive(rem)) continue;
681            Node v = _graph.source(e);
682            if ((*_level)[v] < level) {
683              if (!_level->active(v) && v != _target) {
684                _level->activate(v);
685              }
686              if (!_tolerance.less(rem, excess)) {
687                _flow->set(e, (*_flow)[e] - excess);
688                (*_excess)[v] += excess;
689                excess = 0;
690                goto no_more_push_2;
691              } else {
692                excess -= rem;
693                (*_excess)[v] += rem;
694                _flow->set(e, 0);
695              }
696            } else if (new_level > (*_level)[v]) {
697              new_level = (*_level)[v];
698            }
699          }
700
701        no_more_push_2:
702
703          (*_excess)[n] = excess;
704
705          if (excess != 0) {
706            if (new_level + 1 < _level->maxLevel()) {
707              _level->liftActiveOn(level, new_level + 1);
708            } else {
709              _level->liftActiveToTop(level);
710            }
711            if (_level->emptyLevel(level)) {
712              _level->liftToTop(level);
713            }
714          } else {
715            _level->deactivate(n);
716          }
717        }
718      }
719    first_phase_done:;
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.