lemon/preflow.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sat, 08 Jan 2011 22:49:09 +0100
changeset 1032 62ba43576f85
parent 877 141f9c0db4a3
parent 891 bb70ad62c95f
child 924 a80381c43760
permissions -rw-r--r--
Doc group for TSP algorithms (#386)
     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-2010
     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       while (true) {
   580         int num = _node_num;
   581 
   582         Node n = INVALID;
   583         int level = -1;
   584 
   585         while (num > 0) {
   586           n = _level->highestActive();
   587           if (n == INVALID) goto first_phase_done;
   588           level = _level->highestActiveLevel();
   589           --num;
   590           
   591           Value excess = (*_excess)[n];
   592           int new_level = _level->maxLevel();
   593 
   594           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   595             Value rem = (*_capacity)[e] - (*_flow)[e];
   596             if (!_tolerance.positive(rem)) continue;
   597             Node v = _graph.target(e);
   598             if ((*_level)[v] < level) {
   599               if (!_level->active(v) && v != _target) {
   600                 _level->activate(v);
   601               }
   602               if (!_tolerance.less(rem, excess)) {
   603                 _flow->set(e, (*_flow)[e] + excess);
   604                 (*_excess)[v] += excess;
   605                 excess = 0;
   606                 goto no_more_push_1;
   607               } else {
   608                 excess -= rem;
   609                 (*_excess)[v] += rem;
   610                 _flow->set(e, (*_capacity)[e]);
   611               }
   612             } else if (new_level > (*_level)[v]) {
   613               new_level = (*_level)[v];
   614             }
   615           }
   616 
   617           for (InArcIt e(_graph, n); e != INVALID; ++e) {
   618             Value rem = (*_flow)[e];
   619             if (!_tolerance.positive(rem)) continue;
   620             Node v = _graph.source(e);
   621             if ((*_level)[v] < level) {
   622               if (!_level->active(v) && v != _target) {
   623                 _level->activate(v);
   624               }
   625               if (!_tolerance.less(rem, excess)) {
   626                 _flow->set(e, (*_flow)[e] - excess);
   627                 (*_excess)[v] += excess;
   628                 excess = 0;
   629                 goto no_more_push_1;
   630               } else {
   631                 excess -= rem;
   632                 (*_excess)[v] += rem;
   633                 _flow->set(e, 0);
   634               }
   635             } else if (new_level > (*_level)[v]) {
   636               new_level = (*_level)[v];
   637             }
   638           }
   639 
   640         no_more_push_1:
   641 
   642           (*_excess)[n] = excess;
   643 
   644           if (excess != 0) {
   645             if (new_level + 1 < _level->maxLevel()) {
   646               _level->liftHighestActive(new_level + 1);
   647             } else {
   648               _level->liftHighestActiveToTop();
   649             }
   650             if (_level->emptyLevel(level)) {
   651               _level->liftToTop(level);
   652             }
   653           } else {
   654             _level->deactivate(n);
   655           }
   656         }
   657 
   658         num = _node_num * 20;
   659         while (num > 0) {
   660           while (level >= 0 && _level->activeFree(level)) {
   661             --level;
   662           }
   663           if (level == -1) {
   664             n = _level->highestActive();
   665             level = _level->highestActiveLevel();
   666             if (n == INVALID) goto first_phase_done;
   667           } else {
   668             n = _level->activeOn(level);
   669           }
   670           --num;
   671 
   672           Value excess = (*_excess)[n];
   673           int new_level = _level->maxLevel();
   674 
   675           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   676             Value rem = (*_capacity)[e] - (*_flow)[e];
   677             if (!_tolerance.positive(rem)) continue;
   678             Node v = _graph.target(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, (*_capacity)[e]);
   692               }
   693             } else if (new_level > (*_level)[v]) {
   694               new_level = (*_level)[v];
   695             }
   696           }
   697 
   698           for (InArcIt e(_graph, n); e != INVALID; ++e) {
   699             Value rem = (*_flow)[e];
   700             if (!_tolerance.positive(rem)) continue;
   701             Node v = _graph.source(e);
   702             if ((*_level)[v] < level) {
   703               if (!_level->active(v) && v != _target) {
   704                 _level->activate(v);
   705               }
   706               if (!_tolerance.less(rem, excess)) {
   707                 _flow->set(e, (*_flow)[e] - excess);
   708                 (*_excess)[v] += excess;
   709                 excess = 0;
   710                 goto no_more_push_2;
   711               } else {
   712                 excess -= rem;
   713                 (*_excess)[v] += rem;
   714                 _flow->set(e, 0);
   715               }
   716             } else if (new_level > (*_level)[v]) {
   717               new_level = (*_level)[v];
   718             }
   719           }
   720 
   721         no_more_push_2:
   722 
   723           (*_excess)[n] = excess;
   724 
   725           if (excess != 0) {
   726             if (new_level + 1 < _level->maxLevel()) {
   727               _level->liftActiveOn(level, new_level + 1);
   728             } else {
   729               _level->liftActiveToTop(level);
   730             }
   731             if (_level->emptyLevel(level)) {
   732               _level->liftToTop(level);
   733             }
   734           } else {
   735             _level->deactivate(n);
   736           }
   737         }
   738       }
   739     first_phase_done:;
   740     }
   741 
   742     /// \brief Starts the second phase of the preflow algorithm.
   743     ///
   744     /// The preflow algorithm consists of two phases, this method runs
   745     /// the second phase. After calling one of the \ref init() functions
   746     /// and \ref startFirstPhase() and then \ref startSecondPhase(),
   747     /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
   748     /// value of a maximum flow, \ref minCut() returns a minimum cut
   749     /// \pre One of the \ref init() functions and \ref startFirstPhase()
   750     /// must be called before using this function.
   751     void startSecondPhase() {
   752       _phase = false;
   753 
   754       typename Digraph::template NodeMap<bool> reached(_graph);
   755       for (NodeIt n(_graph); n != INVALID; ++n) {
   756         reached[n] = (*_level)[n] < _level->maxLevel();
   757       }
   758 
   759       _level->initStart();
   760       _level->initAddItem(_source);
   761 
   762       std::vector<Node> queue;
   763       queue.push_back(_source);
   764       reached[_source] = true;
   765 
   766       while (!queue.empty()) {
   767         _level->initNewLevel();
   768         std::vector<Node> nqueue;
   769         for (int i = 0; i < int(queue.size()); ++i) {
   770           Node n = queue[i];
   771           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   772             Node v = _graph.target(e);
   773             if (!reached[v] && _tolerance.positive((*_flow)[e])) {
   774               reached[v] = true;
   775               _level->initAddItem(v);
   776               nqueue.push_back(v);
   777             }
   778           }
   779           for (InArcIt e(_graph, n); e != INVALID; ++e) {
   780             Node u = _graph.source(e);
   781             if (!reached[u] &&
   782                 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
   783               reached[u] = true;
   784               _level->initAddItem(u);
   785               nqueue.push_back(u);
   786             }
   787           }
   788         }
   789         queue.swap(nqueue);
   790       }
   791       _level->initFinish();
   792 
   793       for (NodeIt n(_graph); n != INVALID; ++n) {
   794         if (!reached[n]) {
   795           _level->dirtyTopButOne(n);
   796         } else if ((*_excess)[n] > 0 && _target != n) {
   797           _level->activate(n);
   798         }
   799       }
   800 
   801       Node n;
   802       while ((n = _level->highestActive()) != INVALID) {
   803         Value excess = (*_excess)[n];
   804         int level = _level->highestActiveLevel();
   805         int new_level = _level->maxLevel();
   806 
   807         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   808           Value rem = (*_capacity)[e] - (*_flow)[e];
   809           if (!_tolerance.positive(rem)) continue;
   810           Node v = _graph.target(e);
   811           if ((*_level)[v] < level) {
   812             if (!_level->active(v) && v != _source) {
   813               _level->activate(v);
   814             }
   815             if (!_tolerance.less(rem, excess)) {
   816               _flow->set(e, (*_flow)[e] + excess);
   817               (*_excess)[v] += excess;
   818               excess = 0;
   819               goto no_more_push;
   820             } else {
   821               excess -= rem;
   822               (*_excess)[v] += rem;
   823               _flow->set(e, (*_capacity)[e]);
   824             }
   825           } else if (new_level > (*_level)[v]) {
   826             new_level = (*_level)[v];
   827           }
   828         }
   829 
   830         for (InArcIt e(_graph, n); e != INVALID; ++e) {
   831           Value rem = (*_flow)[e];
   832           if (!_tolerance.positive(rem)) continue;
   833           Node v = _graph.source(e);
   834           if ((*_level)[v] < level) {
   835             if (!_level->active(v) && v != _source) {
   836               _level->activate(v);
   837             }
   838             if (!_tolerance.less(rem, excess)) {
   839               _flow->set(e, (*_flow)[e] - excess);
   840               (*_excess)[v] += excess;
   841               excess = 0;
   842               goto no_more_push;
   843             } else {
   844               excess -= rem;
   845               (*_excess)[v] += rem;
   846               _flow->set(e, 0);
   847             }
   848           } else if (new_level > (*_level)[v]) {
   849             new_level = (*_level)[v];
   850           }
   851         }
   852 
   853       no_more_push:
   854 
   855         (*_excess)[n] = excess;
   856 
   857         if (excess != 0) {
   858           if (new_level + 1 < _level->maxLevel()) {
   859             _level->liftHighestActive(new_level + 1);
   860           } else {
   861             // Calculation error
   862             _level->liftHighestActiveToTop();
   863           }
   864           if (_level->emptyLevel(level)) {
   865             // Calculation error
   866             _level->liftToTop(level);
   867           }
   868         } else {
   869           _level->deactivate(n);
   870         }
   871 
   872       }
   873     }
   874 
   875     /// \brief Runs the preflow algorithm.
   876     ///
   877     /// Runs the preflow algorithm.
   878     /// \note pf.run() is just a shortcut of the following code.
   879     /// \code
   880     ///   pf.init();
   881     ///   pf.startFirstPhase();
   882     ///   pf.startSecondPhase();
   883     /// \endcode
   884     void run() {
   885       init();
   886       startFirstPhase();
   887       startSecondPhase();
   888     }
   889 
   890     /// \brief Runs the preflow algorithm to compute the minimum cut.
   891     ///
   892     /// Runs the preflow algorithm to compute the minimum cut.
   893     /// \note pf.runMinCut() is just a shortcut of the following code.
   894     /// \code
   895     ///   pf.init();
   896     ///   pf.startFirstPhase();
   897     /// \endcode
   898     void runMinCut() {
   899       init();
   900       startFirstPhase();
   901     }
   902 
   903     /// @}
   904 
   905     /// \name Query Functions
   906     /// The results of the preflow algorithm can be obtained using these
   907     /// functions.\n
   908     /// Either one of the \ref run() "run*()" functions or one of the
   909     /// \ref startFirstPhase() "start*()" functions should be called
   910     /// before using them.
   911 
   912     ///@{
   913 
   914     /// \brief Returns the value of the maximum flow.
   915     ///
   916     /// Returns the value of the maximum flow by returning the excess
   917     /// of the target node. This value equals to the value of
   918     /// the maximum flow already after the first phase of the algorithm.
   919     ///
   920     /// \pre Either \ref run() or \ref init() must be called before
   921     /// using this function.
   922     Value flowValue() const {
   923       return (*_excess)[_target];
   924     }
   925 
   926     /// \brief Returns the flow value on the given arc.
   927     ///
   928     /// Returns the flow value on the given arc. This method can
   929     /// be called after the second phase of the algorithm.
   930     ///
   931     /// \pre Either \ref run() or \ref init() must be called before
   932     /// using this function.
   933     Value flow(const Arc& arc) const {
   934       return (*_flow)[arc];
   935     }
   936 
   937     /// \brief Returns a const reference to the flow map.
   938     ///
   939     /// Returns a const reference to the arc map storing the found flow.
   940     /// This method can be called after the second phase of the algorithm.
   941     ///
   942     /// \pre Either \ref run() or \ref init() must be called before
   943     /// using this function.
   944     const FlowMap& flowMap() const {
   945       return *_flow;
   946     }
   947 
   948     /// \brief Returns \c true when the node is on the source side of the
   949     /// minimum cut.
   950     ///
   951     /// Returns true when the node is on the source side of the found
   952     /// minimum cut. This method can be called both after running \ref
   953     /// startFirstPhase() and \ref startSecondPhase().
   954     ///
   955     /// \pre Either \ref run() or \ref init() must be called before
   956     /// using this function.
   957     bool minCut(const Node& node) const {
   958       return ((*_level)[node] == _level->maxLevel()) == _phase;
   959     }
   960 
   961     /// \brief Gives back a minimum value cut.
   962     ///
   963     /// Sets \c cutMap to the characteristic vector of a minimum value
   964     /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
   965     /// node map with \c bool (or convertible) value type.
   966     ///
   967     /// This method can be called both after running \ref startFirstPhase()
   968     /// and \ref startSecondPhase(). The result after the second phase
   969     /// could be slightly different if inexact computation is used.
   970     ///
   971     /// \note This function calls \ref minCut() for each node, so it runs in
   972     /// O(n) time.
   973     ///
   974     /// \pre Either \ref run() or \ref init() must be called before
   975     /// using this function.
   976     template <typename CutMap>
   977     void minCutMap(CutMap& cutMap) const {
   978       for (NodeIt n(_graph); n != INVALID; ++n) {
   979         cutMap.set(n, minCut(n));
   980       }
   981     }
   982 
   983     /// @}
   984   };
   985 }
   986 
   987 #endif