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