COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/preflow.h @ 611:85cb3aa71cce

Last change on this file since 611:85cb3aa71cce was 611:85cb3aa71cce, checked in by Alpar Juttner <alpar@…>, 16 years ago

Merge and fix

File size: 29.5 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 Flow;
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<Flow> 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<Flow> 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::Flow Flow;
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<Flow> 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) const {
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        Flow 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        Flow rem = (*_capacity)[e] - (*_flow)[e];
523        if (_tolerance.positive(rem)) {
524          Node u = _graph.target(e);
525          if ((*_level)[u] == _level->maxLevel()) continue;
526          _flow->set(e, (*_capacity)[e]);
527          (*_excess)[u] += rem;
528          if (u != _target && !_level->active(u)) {
529            _level->activate(u);
530          }
531        }
532      }
533      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
534        Flow rem = (*_flow)[e];
535        if (_tolerance.positive(rem)) {
536          Node v = _graph.source(e);
537          if ((*_level)[v] == _level->maxLevel()) continue;
538          _flow->set(e, 0);
539          (*_excess)[v] += rem;
540          if (v != _target && !_level->active(v)) {
541            _level->activate(v);
542          }
543        }
544      }
545      return true;
546    }
547
548    /// \brief Starts the first phase of the preflow algorithm.
549    ///
550    /// The preflow algorithm consists of two phases, this method runs
551    /// the first phase. After the first phase the maximum flow value
552    /// and a minimum value cut can already be computed, although a
553    /// maximum flow is not yet obtained. So after calling this method
554    /// \ref flowValue() returns the value of a maximum flow and \ref
555    /// minCut() returns a minimum cut.
556    /// \pre One of the \ref init() functions must be called before
557    /// using this function.
558    void startFirstPhase() {
559      _phase = true;
560
561      Node n = _level->highestActive();
562      int level = _level->highestActiveLevel();
563      while (n != INVALID) {
564        int num = _node_num;
565
566        while (num > 0 && n != INVALID) {
567          Flow excess = (*_excess)[n];
568          int new_level = _level->maxLevel();
569
570          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
571            Flow rem = (*_capacity)[e] - (*_flow)[e];
572            if (!_tolerance.positive(rem)) continue;
573            Node v = _graph.target(e);
574            if ((*_level)[v] < level) {
575              if (!_level->active(v) && v != _target) {
576                _level->activate(v);
577              }
578              if (!_tolerance.less(rem, excess)) {
579                _flow->set(e, (*_flow)[e] + excess);
580                (*_excess)[v] += excess;
581                excess = 0;
582                goto no_more_push_1;
583              } else {
584                excess -= rem;
585                (*_excess)[v] += rem;
586                _flow->set(e, (*_capacity)[e]);
587              }
588            } else if (new_level > (*_level)[v]) {
589              new_level = (*_level)[v];
590            }
591          }
592
593          for (InArcIt e(_graph, n); e != INVALID; ++e) {
594            Flow rem = (*_flow)[e];
595            if (!_tolerance.positive(rem)) continue;
596            Node v = _graph.source(e);
597            if ((*_level)[v] < level) {
598              if (!_level->active(v) && v != _target) {
599                _level->activate(v);
600              }
601              if (!_tolerance.less(rem, excess)) {
602                _flow->set(e, (*_flow)[e] - excess);
603                (*_excess)[v] += excess;
604                excess = 0;
605                goto no_more_push_1;
606              } else {
607                excess -= rem;
608                (*_excess)[v] += rem;
609                _flow->set(e, 0);
610              }
611            } else if (new_level > (*_level)[v]) {
612              new_level = (*_level)[v];
613            }
614          }
615
616        no_more_push_1:
617
618          (*_excess)[n] = excess;
619
620          if (excess != 0) {
621            if (new_level + 1 < _level->maxLevel()) {
622              _level->liftHighestActive(new_level + 1);
623            } else {
624              _level->liftHighestActiveToTop();
625            }
626            if (_level->emptyLevel(level)) {
627              _level->liftToTop(level);
628            }
629          } else {
630            _level->deactivate(n);
631          }
632
633          n = _level->highestActive();
634          level = _level->highestActiveLevel();
635          --num;
636        }
637
638        num = _node_num * 20;
639        while (num > 0 && n != INVALID) {
640          Flow excess = (*_excess)[n];
641          int new_level = _level->maxLevel();
642
643          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
644            Flow rem = (*_capacity)[e] - (*_flow)[e];
645            if (!_tolerance.positive(rem)) continue;
646            Node v = _graph.target(e);
647            if ((*_level)[v] < level) {
648              if (!_level->active(v) && v != _target) {
649                _level->activate(v);
650              }
651              if (!_tolerance.less(rem, excess)) {
652                _flow->set(e, (*_flow)[e] + excess);
653                (*_excess)[v] += excess;
654                excess = 0;
655                goto no_more_push_2;
656              } else {
657                excess -= rem;
658                (*_excess)[v] += rem;
659                _flow->set(e, (*_capacity)[e]);
660              }
661            } else if (new_level > (*_level)[v]) {
662              new_level = (*_level)[v];
663            }
664          }
665
666          for (InArcIt e(_graph, n); e != INVALID; ++e) {
667            Flow rem = (*_flow)[e];
668            if (!_tolerance.positive(rem)) continue;
669            Node v = _graph.source(e);
670            if ((*_level)[v] < level) {
671              if (!_level->active(v) && v != _target) {
672                _level->activate(v);
673              }
674              if (!_tolerance.less(rem, excess)) {
675                _flow->set(e, (*_flow)[e] - excess);
676                (*_excess)[v] += excess;
677                excess = 0;
678                goto no_more_push_2;
679              } else {
680                excess -= rem;
681                (*_excess)[v] += rem;
682                _flow->set(e, 0);
683              }
684            } else if (new_level > (*_level)[v]) {
685              new_level = (*_level)[v];
686            }
687          }
688
689        no_more_push_2:
690
691          (*_excess)[n] = excess;
692
693          if (excess != 0) {
694            if (new_level + 1 < _level->maxLevel()) {
695              _level->liftActiveOn(level, new_level + 1);
696            } else {
697              _level->liftActiveToTop(level);
698            }
699            if (_level->emptyLevel(level)) {
700              _level->liftToTop(level);
701            }
702          } else {
703            _level->deactivate(n);
704          }
705
706          while (level >= 0 && _level->activeFree(level)) {
707            --level;
708          }
709          if (level == -1) {
710            n = _level->highestActive();
711            level = _level->highestActiveLevel();
712          } else {
713            n = _level->activeOn(level);
714          }
715          --num;
716        }
717      }
718    }
719
720    /// \brief Starts the second phase of the preflow algorithm.
721    ///
722    /// The preflow algorithm consists of two phases, this method runs
723    /// the second phase. After calling one of the \ref init() functions
724    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
725    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
726    /// value of a maximum flow, \ref minCut() returns a minimum cut
727    /// \pre One of the \ref init() functions and \ref startFirstPhase()
728    /// must be called before using this function.
729    void startSecondPhase() {
730      _phase = false;
731
732      typename Digraph::template NodeMap<bool> reached(_graph);
733      for (NodeIt n(_graph); n != INVALID; ++n) {
734        reached[n] = (*_level)[n] < _level->maxLevel();
735      }
736
737      _level->initStart();
738      _level->initAddItem(_source);
739
740      std::vector<Node> queue;
741      queue.push_back(_source);
742      reached[_source] = true;
743
744      while (!queue.empty()) {
745        _level->initNewLevel();
746        std::vector<Node> nqueue;
747        for (int i = 0; i < int(queue.size()); ++i) {
748          Node n = queue[i];
749          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
750            Node v = _graph.target(e);
751            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
752              reached[v] = true;
753              _level->initAddItem(v);
754              nqueue.push_back(v);
755            }
756          }
757          for (InArcIt e(_graph, n); e != INVALID; ++e) {
758            Node u = _graph.source(e);
759            if (!reached[u] &&
760                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
761              reached[u] = true;
762              _level->initAddItem(u);
763              nqueue.push_back(u);
764            }
765          }
766        }
767        queue.swap(nqueue);
768      }
769      _level->initFinish();
770
771      for (NodeIt n(_graph); n != INVALID; ++n) {
772        if (!reached[n]) {
773          _level->dirtyTopButOne(n);
774        } else if ((*_excess)[n] > 0 && _target != n) {
775          _level->activate(n);
776        }
777      }
778
779      Node n;
780      while ((n = _level->highestActive()) != INVALID) {
781        Flow excess = (*_excess)[n];
782        int level = _level->highestActiveLevel();
783        int new_level = _level->maxLevel();
784
785        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
786          Flow rem = (*_capacity)[e] - (*_flow)[e];
787          if (!_tolerance.positive(rem)) continue;
788          Node v = _graph.target(e);
789          if ((*_level)[v] < level) {
790            if (!_level->active(v) && v != _source) {
791              _level->activate(v);
792            }
793            if (!_tolerance.less(rem, excess)) {
794              _flow->set(e, (*_flow)[e] + excess);
795              (*_excess)[v] += excess;
796              excess = 0;
797              goto no_more_push;
798            } else {
799              excess -= rem;
800              (*_excess)[v] += rem;
801              _flow->set(e, (*_capacity)[e]);
802            }
803          } else if (new_level > (*_level)[v]) {
804            new_level = (*_level)[v];
805          }
806        }
807
808        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809          Flow rem = (*_flow)[e];
810          if (!_tolerance.positive(rem)) continue;
811          Node v = _graph.source(e);
812          if ((*_level)[v] < level) {
813            if (!_level->active(v) && v != _source) {
814              _level->activate(v);
815            }
816            if (!_tolerance.less(rem, excess)) {
817              _flow->set(e, (*_flow)[e] - excess);
818              (*_excess)[v] += excess;
819              excess = 0;
820              goto no_more_push;
821            } else {
822              excess -= rem;
823              (*_excess)[v] += rem;
824              _flow->set(e, 0);
825            }
826          } else if (new_level > (*_level)[v]) {
827            new_level = (*_level)[v];
828          }
829        }
830
831      no_more_push:
832
833        (*_excess)[n] = excess;
834
835        if (excess != 0) {
836          if (new_level + 1 < _level->maxLevel()) {
837            _level->liftHighestActive(new_level + 1);
838          } else {
839            // Calculation error
840            _level->liftHighestActiveToTop();
841          }
842          if (_level->emptyLevel(level)) {
843            // Calculation error
844            _level->liftToTop(level);
845          }
846        } else {
847          _level->deactivate(n);
848        }
849
850      }
851    }
852
853    /// \brief Runs the preflow algorithm.
854    ///
855    /// Runs the preflow algorithm.
856    /// \note pf.run() is just a shortcut of the following code.
857    /// \code
858    ///   pf.init();
859    ///   pf.startFirstPhase();
860    ///   pf.startSecondPhase();
861    /// \endcode
862    void run() {
863      init();
864      startFirstPhase();
865      startSecondPhase();
866    }
867
868    /// \brief Runs the preflow algorithm to compute the minimum cut.
869    ///
870    /// Runs the preflow algorithm to compute the minimum cut.
871    /// \note pf.runMinCut() is just a shortcut of the following code.
872    /// \code
873    ///   pf.init();
874    ///   pf.startFirstPhase();
875    /// \endcode
876    void runMinCut() {
877      init();
878      startFirstPhase();
879    }
880
881    /// @}
882
883    /// \name Query Functions
884    /// The results of the preflow algorithm can be obtained using these
885    /// functions.\n
886    /// Either one of the \ref run() "run*()" functions or one of the
887    /// \ref startFirstPhase() "start*()" functions should be called
888    /// before using them.
889
890    ///@{
891
892    /// \brief Returns the value of the maximum flow.
893    ///
894    /// Returns the value of the maximum flow by returning the excess
895    /// of the target node. This value equals to the value of
896    /// the maximum flow already after the first phase of the algorithm.
897    ///
898    /// \pre Either \ref run() or \ref init() must be called before
899    /// using this function.
900    Flow flowValue() const {
901      return (*_excess)[_target];
902    }
903
904    /// \brief Returns the flow on the given arc.
905    ///
906    /// Returns the flow on the given arc. This method can
907    /// be called after the second phase of the algorithm.
908    ///
909    /// \pre Either \ref run() or \ref init() must be called before
910    /// using this function.
911    Flow flow(const Arc& arc) const {
912      return (*_flow)[arc];
913    }
914
915    /// \brief Returns a const reference to the flow map.
916    ///
917    /// Returns a const reference to the arc map storing the found flow.
918    /// This method can be called after the second phase of the algorithm.
919    ///
920    /// \pre Either \ref run() or \ref init() must be called before
921    /// using this function.
922    const FlowMap& flowMap() const {
923      return *_flow;
924    }
925
926    /// \brief Returns \c true when the node is on the source side of the
927    /// minimum cut.
928    ///
929    /// Returns true when the node is on the source side of the found
930    /// minimum cut. This method can be called both after running \ref
931    /// startFirstPhase() and \ref startSecondPhase().
932    ///
933    /// \pre Either \ref run() or \ref init() must be called before
934    /// using this function.
935    bool minCut(const Node& node) const {
936      return ((*_level)[node] == _level->maxLevel()) == _phase;
937    }
938
939    /// \brief Gives back a minimum value cut.
940    ///
941    /// Sets \c cutMap to the characteristic vector of a minimum value
942    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
943    /// node map with \c bool (or convertible) value type.
944    ///
945    /// This method can be called both after running \ref startFirstPhase()
946    /// and \ref startSecondPhase(). The result after the second phase
947    /// could be slightly different if inexact computation is used.
948    ///
949    /// \note This function calls \ref minCut() for each node, so it runs in
950    /// O(n) time.
951    ///
952    /// \pre Either \ref run() or \ref init() must be called before
953    /// using this function.
954    template <typename CutMap>
955    void minCutMap(CutMap& cutMap) const {
956      for (NodeIt n(_graph); n != INVALID; ++n) {
957        cutMap.set(n, minCut(n));
958      }
959    }
960
961    /// @}
962  };
963}
964
965#endif
Note: See TracBrowser for help on using the repository browser.