COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/preflow.h @ 717:684964884a2e

Last change on this file since 717:684964884a2e was 715:ece80147fb08, checked in by Alpar Juttner <alpar@…>, 15 years ago

Merge

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