COIN-OR::LEMON - Graph Library

source: lemon/lemon/preflow.h @ 481:7afc121e0689

Last change on this file since 481:7afc121e0689 was 463:88ed40ad0d4f, checked in by Alpar Juttner <alpar@…>, 11 years ago

Happy New Year again

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