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