COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 2553:bfced05fa852

Last change on this file since 2553:bfced05fa852 was 2553:bfced05fa852, checked in by Alpar Juttner, 12 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

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