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