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