COIN-OR::LEMON - Graph Library

source: lemon/lemon/preflow.h @ 404:660db48f324f

Last change on this file since 404:660db48f324f was 404:660db48f324f, checked in by Alpar Juttner <alpar@…>, 11 years ago

Port preflow push max flow alg. from svn -r3516 (#176)
Namely,

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