lemon/preflow.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 17 Apr 2009 18:14:35 +0200
changeset 657 dacc2cee2b4c
parent 525 9605e051942f
child 658 85cb3aa71cce
permissions -rw-r--r--
Slightly modify the interface of Circulation and Preflow (#266)
in order to synchronize them to the interface of NetworkSimplex.

Circulation:
- The "delta" notation is replaced by "supply".
- lowerCapMap(), upperCapMap() are renamed to lowerMap() and upperMap().
- Value is renamed to Flow.

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