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