COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 2521:05c0ba99cc27

Last change on this file since 2521:05c0ba99cc27 was 2518:4c0a23bd70b5, checked in by Balazs Dezso, 12 years ago

Bugfix in min cut computation

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 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^3)\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        std::vector<Node> nqueue;
411        for (int i = 0; i < int(queue.size()); ++i) {
412          Node n = queue[i];
413          for (InEdgeIt 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 (OutEdgeIt 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 flowInit(const FlowMap& flowMap) {
449      createStructures();
450
451      for (EdgeIt 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 (InEdgeIt e(_graph, n); e != INVALID; ++e) {
458          excess += (*_flow)[e];
459        }
460        for (OutEdgeIt 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 Graph::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 (InEdgeIt 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 (OutEdgeIt 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 (OutEdgeIt 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 (InEdgeIt 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 (OutEdgeIt 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 (InEdgeIt 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 (OutEdgeIt 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 (InEdgeIt 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 Graph::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 (OutEdgeIt 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 (InEdgeIt 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->markToBottom(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 (OutEdgeIt 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 (InEdgeIt 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 %Dijkstra 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 edge.
906    ///
907    /// Sets the \c flowMap to the flow on the edges. This method can
908    /// be called after the second phase of algorithm.
909    Value flow(const Edge& edge) const {
910      return (*_flow)[edge];
911    }
912   
913    /// @}   
914  };
915}
916
917#endif
Note: See TracBrowser for help on using the repository browser.