COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 2522:616c019215c4

Last change on this file since 2522:616c019215c4 was 2522:616c019215c4, checked in by Balazs Dezso, 12 years ago

Performance bug in Preflow
The initial relabeling moved each node to the lowest level
Doc bug fix

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