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