lemon/preflow.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sat, 20 Feb 2010 18:39:03 +0100
changeset 839 f3bc4e9b5f3a
parent 788 c92296660262
child 825 75e6020b19b1
permissions -rw-r--r--
New heuristics for MCF algorithms (#340)
and some implementation improvements.

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