COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/preflow.h @ 392:db3251947eba

Last change on this file since 392:db3251947eba was 392:db3251947eba, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Rename flowInit() to init() in Preflow (#176)

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