alpar@389: /* -*- mode: C++; indent-tabs-mode: nil; -*-
alpar@389:  *
alpar@389:  * This file is a part of LEMON, a generic C++ optimization library.
alpar@389:  *
alpar@389:  * Copyright (C) 2003-2008
alpar@389:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@389:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@389:  *
alpar@389:  * Permission to use, modify and distribute this software is granted
alpar@389:  * provided that this copyright notice appears in all copies. For
alpar@389:  * precise terms see the accompanying LICENSE file.
alpar@389:  *
alpar@389:  * This software is provided "AS IS" with no warranty of any kind,
alpar@389:  * express or implied, and with no claim as to its suitability for any
alpar@389:  * purpose.
alpar@389:  *
alpar@389:  */
alpar@389: 
alpar@389: #ifndef LEMON_PREFLOW_H
alpar@389: #define LEMON_PREFLOW_H
alpar@389: 
alpar@389: #include <lemon/tolerance.h>
alpar@389: #include <lemon/elevator.h>
alpar@389: 
alpar@389: /// \file
alpar@389: /// \ingroup max_flow
alpar@389: /// \brief Implementation of the preflow algorithm.
alpar@389: 
alpar@389: namespace lemon {
alpar@389: 
alpar@389:   /// \brief Default traits class of Preflow class.
alpar@389:   ///
alpar@389:   /// Default traits class of Preflow class.
alpar@389:   /// \param _Graph Digraph type.
alpar@389:   /// \param _CapacityMap Type of capacity map.
alpar@389:   template <typename _Graph, typename _CapacityMap>
alpar@389:   struct PreflowDefaultTraits {
alpar@389: 
alpar@389:     /// \brief The digraph type the algorithm runs on.
alpar@389:     typedef _Graph Digraph;
alpar@389: 
alpar@389:     /// \brief The type of the map that stores the arc capacities.
alpar@389:     ///
alpar@389:     /// The type of the map that stores the arc capacities.
alpar@389:     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
alpar@389:     typedef _CapacityMap CapacityMap;
alpar@389: 
alpar@389:     /// \brief The type of the length of the arcs.
alpar@389:     typedef typename CapacityMap::Value Value;
alpar@389: 
alpar@389:     /// \brief The map type that stores the flow values.
alpar@389:     ///
alpar@389:     /// The map type that stores the flow values.
alpar@389:     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
alpar@389:     typedef typename Digraph::template ArcMap<Value> FlowMap;
alpar@389: 
alpar@389:     /// \brief Instantiates a FlowMap.
alpar@389:     ///
alpar@389:     /// This function instantiates a \ref FlowMap.
alpar@389:     /// \param digraph The digraph, to which we would like to define
alpar@389:     /// the flow map.
alpar@389:     static FlowMap* createFlowMap(const Digraph& digraph) {
alpar@389:       return new FlowMap(digraph);
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief The eleavator type used by Preflow algorithm.
alpar@389:     ///
alpar@389:     /// The elevator type used by Preflow algorithm.
alpar@389:     ///
alpar@389:     /// \sa Elevator
alpar@389:     /// \sa LinkedElevator
alpar@389:     typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
alpar@389: 
alpar@389:     /// \brief Instantiates an Elevator.
alpar@389:     ///
alpar@389:     /// This function instantiates a \ref Elevator.
alpar@389:     /// \param digraph The digraph, to which we would like to define
alpar@389:     /// the elevator.
alpar@389:     /// \param max_level The maximum level of the elevator.
alpar@389:     static Elevator* createElevator(const Digraph& digraph, int max_level) {
alpar@389:       return new Elevator(digraph, max_level);
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief The tolerance used by the algorithm
alpar@389:     ///
alpar@389:     /// The tolerance used by the algorithm to handle inexact computation.
alpar@389:     typedef lemon::Tolerance<Value> Tolerance;
alpar@389: 
alpar@389:   };
alpar@389: 
alpar@389: 
alpar@389:   /// \ingroup max_flow
alpar@389:   ///
alpar@389:   /// \brief %Preflow algorithms class.
alpar@389:   ///
alpar@389:   /// This class provides an implementation of the Goldberg's \e
alpar@389:   /// preflow \e algorithm producing a flow of maximum value in a
alpar@389:   /// digraph. The preflow algorithms are the fastest known max
alpar@389:   /// flow algorithms. The current implementation use a mixture of the
alpar@389:   /// \e "highest label" and the \e "bound decrease" heuristics.
alpar@389:   /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
alpar@389:   ///
alpar@389:   /// The algorithm consists from two phases. After the first phase
alpar@389:   /// the maximal flow value and the minimum cut can be obtained. The
alpar@389:   /// second phase constructs the feasible maximum flow on each arc.
alpar@389:   ///
alpar@389:   /// \param _Graph The digraph type the algorithm runs on.
alpar@389:   /// \param _CapacityMap The flow map type.
alpar@389:   /// \param _Traits Traits class to set various data types used by
alpar@389:   /// the algorithm.  The default traits class is \ref
alpar@389:   /// PreflowDefaultTraits.  See \ref PreflowDefaultTraits for the
alpar@389:   /// documentation of a %Preflow traits class.
alpar@389:   ///
alpar@389:   ///\author Jacint Szabo and Balazs Dezso
alpar@389: #ifdef DOXYGEN
alpar@389:   template <typename _Graph, typename _CapacityMap, typename _Traits>
alpar@389: #else
alpar@389:   template <typename _Graph,
alpar@389:             typename _CapacityMap = typename _Graph::template ArcMap<int>,
alpar@389:             typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
alpar@389: #endif
alpar@389:   class Preflow {
alpar@389:   public:
alpar@389: 
alpar@389:     typedef _Traits Traits;
alpar@389:     typedef typename Traits::Digraph Digraph;
alpar@389:     typedef typename Traits::CapacityMap CapacityMap;
alpar@389:     typedef typename Traits::Value Value;
alpar@389: 
alpar@389:     typedef typename Traits::FlowMap FlowMap;
alpar@389:     typedef typename Traits::Elevator Elevator;
alpar@389:     typedef typename Traits::Tolerance Tolerance;
alpar@389: 
alpar@389:   private:
alpar@389: 
alpar@389:     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
alpar@389: 
alpar@389:     const Digraph& _graph;
alpar@389:     const CapacityMap* _capacity;
alpar@389: 
alpar@389:     int _node_num;
alpar@389: 
alpar@389:     Node _source, _target;
alpar@389: 
alpar@389:     FlowMap* _flow;
alpar@389:     bool _local_flow;
alpar@389: 
alpar@389:     Elevator* _level;
alpar@389:     bool _local_level;
alpar@389: 
alpar@389:     typedef typename Digraph::template NodeMap<Value> ExcessMap;
alpar@389:     ExcessMap* _excess;
alpar@389: 
alpar@389:     Tolerance _tolerance;
alpar@389: 
alpar@389:     bool _phase;
alpar@389: 
alpar@389: 
alpar@389:     void createStructures() {
alpar@389:       _node_num = countNodes(_graph);
alpar@389: 
alpar@389:       if (!_flow) {
alpar@389:         _flow = Traits::createFlowMap(_graph);
alpar@389:         _local_flow = true;
alpar@389:       }
alpar@389:       if (!_level) {
alpar@389:         _level = Traits::createElevator(_graph, _node_num);
alpar@389:         _local_level = true;
alpar@389:       }
alpar@389:       if (!_excess) {
alpar@389:         _excess = new ExcessMap(_graph);
alpar@389:       }
alpar@389:     }
alpar@389: 
alpar@389:     void destroyStructures() {
alpar@389:       if (_local_flow) {
alpar@389:         delete _flow;
alpar@389:       }
alpar@389:       if (_local_level) {
alpar@389:         delete _level;
alpar@389:       }
alpar@389:       if (_excess) {
alpar@389:         delete _excess;
alpar@389:       }
alpar@389:     }
alpar@389: 
alpar@389:   public:
alpar@389: 
alpar@389:     typedef Preflow Create;
alpar@389: 
alpar@389:     ///\name Named template parameters
alpar@389: 
alpar@389:     ///@{
alpar@389: 
alpar@389:     template <typename _FlowMap>
alpar@391:     struct SetFlowMapTraits : public Traits {
alpar@389:       typedef _FlowMap FlowMap;
alpar@389:       static FlowMap *createFlowMap(const Digraph&) {
alpar@390:         LEMON_ASSERT(false, "FlowMap is not initialized");
alpar@390:         return 0; // ignore warnings
alpar@389:       }
alpar@389:     };
alpar@389: 
alpar@389:     /// \brief \ref named-templ-param "Named parameter" for setting
alpar@389:     /// FlowMap type
alpar@389:     ///
alpar@389:     /// \ref named-templ-param "Named parameter" for setting FlowMap
alpar@389:     /// type
alpar@389:     template <typename _FlowMap>
alpar@391:     struct SetFlowMap
alpar@391:       : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<_FlowMap> > {
alpar@389:       typedef Preflow<Digraph, CapacityMap,
alpar@391:                       SetFlowMapTraits<_FlowMap> > Create;
alpar@389:     };
alpar@389: 
alpar@389:     template <typename _Elevator>
alpar@391:     struct SetElevatorTraits : public Traits {
alpar@389:       typedef _Elevator Elevator;
alpar@389:       static Elevator *createElevator(const Digraph&, int) {
alpar@390:         LEMON_ASSERT(false, "Elevator is not initialized");
alpar@390:         return 0; // ignore warnings
alpar@389:       }
alpar@389:     };
alpar@389: 
alpar@389:     /// \brief \ref named-templ-param "Named parameter" for setting
alpar@389:     /// Elevator type
alpar@389:     ///
alpar@389:     /// \ref named-templ-param "Named parameter" for setting Elevator
alpar@389:     /// type
alpar@389:     template <typename _Elevator>
alpar@391:     struct SetElevator
alpar@391:       : public Preflow<Digraph, CapacityMap, SetElevatorTraits<_Elevator> > {
alpar@389:       typedef Preflow<Digraph, CapacityMap,
alpar@391:                       SetElevatorTraits<_Elevator> > Create;
alpar@389:     };
alpar@389: 
alpar@389:     template <typename _Elevator>
alpar@391:     struct SetStandardElevatorTraits : public Traits {
alpar@389:       typedef _Elevator Elevator;
alpar@389:       static Elevator *createElevator(const Digraph& digraph, int max_level) {
alpar@389:         return new Elevator(digraph, max_level);
alpar@389:       }
alpar@389:     };
alpar@389: 
alpar@389:     /// \brief \ref named-templ-param "Named parameter" for setting
alpar@389:     /// Elevator type
alpar@389:     ///
alpar@389:     /// \ref named-templ-param "Named parameter" for setting Elevator
alpar@389:     /// type. The Elevator should be standard constructor interface, ie.
alpar@389:     /// the digraph and the maximum level should be passed to it.
alpar@389:     template <typename _Elevator>
alpar@391:     struct SetStandardElevator
alpar@389:       : public Preflow<Digraph, CapacityMap,
alpar@391:                        SetStandardElevatorTraits<_Elevator> > {
alpar@389:       typedef Preflow<Digraph, CapacityMap,
alpar@391:                       SetStandardElevatorTraits<_Elevator> > Create;
alpar@389:     };
alpar@389: 
alpar@389:     /// @}
alpar@389: 
alpar@389:   protected:
alpar@389: 
alpar@389:     Preflow() {}
alpar@389: 
alpar@389:   public:
alpar@389: 
alpar@389: 
alpar@389:     /// \brief The constructor of the class.
alpar@389:     ///
alpar@389:     /// The constructor of the class.
alpar@389:     /// \param digraph The digraph the algorithm runs on.
alpar@389:     /// \param capacity The capacity of the arcs.
alpar@389:     /// \param source The source node.
alpar@389:     /// \param target The target node.
alpar@389:     Preflow(const Digraph& digraph, const CapacityMap& capacity,
alpar@389:                Node source, Node target)
alpar@389:       : _graph(digraph), _capacity(&capacity),
alpar@389:         _node_num(0), _source(source), _target(target),
alpar@389:         _flow(0), _local_flow(false),
alpar@389:         _level(0), _local_level(false),
alpar@389:         _excess(0), _tolerance(), _phase() {}
alpar@389: 
alpar@389:     /// \brief Destrcutor.
alpar@389:     ///
alpar@389:     /// Destructor.
alpar@389:     ~Preflow() {
alpar@389:       destroyStructures();
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Sets the capacity map.
alpar@389:     ///
alpar@389:     /// Sets the capacity map.
alpar@389:     /// \return \c (*this)
alpar@389:     Preflow& capacityMap(const CapacityMap& map) {
alpar@389:       _capacity = &map;
alpar@389:       return *this;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Sets the flow map.
alpar@389:     ///
alpar@389:     /// Sets the flow map.
alpar@389:     /// \return \c (*this)
alpar@389:     Preflow& flowMap(FlowMap& map) {
alpar@389:       if (_local_flow) {
alpar@389:         delete _flow;
alpar@389:         _local_flow = false;
alpar@389:       }
alpar@389:       _flow = &map;
alpar@389:       return *this;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Returns the flow map.
alpar@389:     ///
alpar@389:     /// \return The flow map.
alpar@389:     const FlowMap& flowMap() {
alpar@389:       return *_flow;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Sets the elevator.
alpar@389:     ///
alpar@389:     /// Sets the elevator.
alpar@389:     /// \return \c (*this)
alpar@389:     Preflow& elevator(Elevator& elevator) {
alpar@389:       if (_local_level) {
alpar@389:         delete _level;
alpar@389:         _local_level = false;
alpar@389:       }
alpar@389:       _level = &elevator;
alpar@389:       return *this;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Returns the elevator.
alpar@389:     ///
alpar@389:     /// \return The elevator.
alpar@389:     const Elevator& elevator() {
alpar@389:       return *_level;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Sets the source node.
alpar@389:     ///
alpar@389:     /// Sets the source node.
alpar@389:     /// \return \c (*this)
alpar@389:     Preflow& source(const Node& node) {
alpar@389:       _source = node;
alpar@389:       return *this;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Sets the target node.
alpar@389:     ///
alpar@389:     /// Sets the target node.
alpar@389:     /// \return \c (*this)
alpar@389:     Preflow& target(const Node& node) {
alpar@389:       _target = node;
alpar@389:       return *this;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Sets the tolerance used by algorithm.
alpar@389:     ///
alpar@389:     /// Sets the tolerance used by algorithm.
alpar@389:     Preflow& tolerance(const Tolerance& tolerance) const {
alpar@389:       _tolerance = tolerance;
alpar@389:       return *this;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Returns the tolerance used by algorithm.
alpar@389:     ///
alpar@389:     /// Returns the tolerance used by algorithm.
alpar@389:     const Tolerance& tolerance() const {
alpar@389:       return tolerance;
alpar@389:     }
alpar@389: 
alpar@389:     /// \name Execution control The simplest way to execute the
alpar@389:     /// algorithm is to use one of the member functions called \c
alpar@389:     /// run().
alpar@389:     /// \n
alpar@389:     /// If you need more control on initial solution or
alpar@389:     /// execution then you have to call one \ref init() function and then
alpar@389:     /// the startFirstPhase() and if you need the startSecondPhase().
alpar@389: 
alpar@389:     ///@{
alpar@389: 
alpar@389:     /// \brief Initializes the internal data structures.
alpar@389:     ///
alpar@389:     /// Initializes the internal data structures.
alpar@389:     ///
alpar@389:     void init() {
alpar@389:       createStructures();
alpar@389: 
alpar@389:       _phase = true;
alpar@389:       for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@389:         _excess->set(n, 0);
alpar@389:       }
alpar@389: 
alpar@389:       for (ArcIt e(_graph); e != INVALID; ++e) {
alpar@389:         _flow->set(e, 0);
alpar@389:       }
alpar@389: 
alpar@389:       typename Digraph::template NodeMap<bool> reached(_graph, false);
alpar@389: 
alpar@389:       _level->initStart();
alpar@389:       _level->initAddItem(_target);
alpar@389: 
alpar@389:       std::vector<Node> queue;
alpar@389:       reached.set(_source, true);
alpar@389: 
alpar@389:       queue.push_back(_target);
alpar@389:       reached.set(_target, true);
alpar@389:       while (!queue.empty()) {
alpar@389:         _level->initNewLevel();
alpar@389:         std::vector<Node> nqueue;
alpar@389:         for (int i = 0; i < int(queue.size()); ++i) {
alpar@389:           Node n = queue[i];
alpar@389:           for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Node u = _graph.source(e);
alpar@389:             if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
alpar@389:               reached.set(u, true);
alpar@389:               _level->initAddItem(u);
alpar@389:               nqueue.push_back(u);
alpar@389:             }
alpar@389:           }
alpar@389:         }
alpar@389:         queue.swap(nqueue);
alpar@389:       }
alpar@389:       _level->initFinish();
alpar@389: 
alpar@389:       for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
alpar@389:         if (_tolerance.positive((*_capacity)[e])) {
alpar@389:           Node u = _graph.target(e);
alpar@389:           if ((*_level)[u] == _level->maxLevel()) continue;
alpar@389:           _flow->set(e, (*_capacity)[e]);
alpar@389:           _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
alpar@389:           if (u != _target && !_level->active(u)) {
alpar@389:             _level->activate(u);
alpar@389:           }
alpar@389:         }
alpar@389:       }
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Initializes the internal data structures.
alpar@389:     ///
alpar@389:     /// Initializes the internal data structures and sets the initial
alpar@389:     /// flow to the given \c flowMap. The \c flowMap should contain a
alpar@389:     /// flow or at least a preflow, ie. in each node excluding the
alpar@389:     /// target the incoming flow should greater or equal to the
alpar@389:     /// outgoing flow.
alpar@389:     /// \return %False when the given \c flowMap is not a preflow.
alpar@389:     template <typename FlowMap>
kpeter@392:     bool init(const FlowMap& flowMap) {
alpar@389:       createStructures();
alpar@389: 
alpar@389:       for (ArcIt e(_graph); e != INVALID; ++e) {
alpar@389:         _flow->set(e, flowMap[e]);
alpar@389:       }
alpar@389: 
alpar@389:       for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@389:         Value excess = 0;
alpar@389:         for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:           excess += (*_flow)[e];
alpar@389:         }
alpar@389:         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:           excess -= (*_flow)[e];
alpar@389:         }
alpar@389:         if (excess < 0 && n != _source) return false;
alpar@389:         _excess->set(n, excess);
alpar@389:       }
alpar@389: 
alpar@389:       typename Digraph::template NodeMap<bool> reached(_graph, false);
alpar@389: 
alpar@389:       _level->initStart();
alpar@389:       _level->initAddItem(_target);
alpar@389: 
alpar@389:       std::vector<Node> queue;
alpar@389:       reached.set(_source, true);
alpar@389: 
alpar@389:       queue.push_back(_target);
alpar@389:       reached.set(_target, true);
alpar@389:       while (!queue.empty()) {
alpar@389:         _level->initNewLevel();
alpar@389:         std::vector<Node> nqueue;
alpar@389:         for (int i = 0; i < int(queue.size()); ++i) {
alpar@389:           Node n = queue[i];
alpar@389:           for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Node u = _graph.source(e);
alpar@389:             if (!reached[u] &&
alpar@389:                 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
alpar@389:               reached.set(u, true);
alpar@389:               _level->initAddItem(u);
alpar@389:               nqueue.push_back(u);
alpar@389:             }
alpar@389:           }
alpar@389:           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Node v = _graph.target(e);
alpar@389:             if (!reached[v] && _tolerance.positive((*_flow)[e])) {
alpar@389:               reached.set(v, true);
alpar@389:               _level->initAddItem(v);
alpar@389:               nqueue.push_back(v);
alpar@389:             }
alpar@389:           }
alpar@389:         }
alpar@389:         queue.swap(nqueue);
alpar@389:       }
alpar@389:       _level->initFinish();
alpar@389: 
alpar@389:       for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
alpar@389:         Value rem = (*_capacity)[e] - (*_flow)[e];
alpar@389:         if (_tolerance.positive(rem)) {
alpar@389:           Node u = _graph.target(e);
alpar@389:           if ((*_level)[u] == _level->maxLevel()) continue;
alpar@389:           _flow->set(e, (*_capacity)[e]);
alpar@389:           _excess->set(u, (*_excess)[u] + rem);
alpar@389:           if (u != _target && !_level->active(u)) {
alpar@389:             _level->activate(u);
alpar@389:           }
alpar@389:         }
alpar@389:       }
alpar@389:       for (InArcIt e(_graph, _source); e != INVALID; ++e) {
alpar@389:         Value rem = (*_flow)[e];
alpar@389:         if (_tolerance.positive(rem)) {
alpar@389:           Node v = _graph.source(e);
alpar@389:           if ((*_level)[v] == _level->maxLevel()) continue;
alpar@389:           _flow->set(e, 0);
alpar@389:           _excess->set(v, (*_excess)[v] + rem);
alpar@389:           if (v != _target && !_level->active(v)) {
alpar@389:             _level->activate(v);
alpar@389:           }
alpar@389:         }
alpar@389:       }
alpar@389:       return true;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Starts the first phase of the preflow algorithm.
alpar@389:     ///
alpar@389:     /// The preflow algorithm consists of two phases, this method runs
alpar@389:     /// the first phase. After the first phase the maximum flow value
alpar@389:     /// and a minimum value cut can already be computed, although a
alpar@389:     /// maximum flow is not yet obtained. So after calling this method
alpar@389:     /// \ref flowValue() returns the value of a maximum flow and \ref
alpar@389:     /// minCut() returns a minimum cut.
alpar@389:     /// \pre One of the \ref init() functions should be called.
alpar@389:     void startFirstPhase() {
alpar@389:       _phase = true;
alpar@389: 
alpar@389:       Node n = _level->highestActive();
alpar@389:       int level = _level->highestActiveLevel();
alpar@389:       while (n != INVALID) {
alpar@389:         int num = _node_num;
alpar@389: 
alpar@389:         while (num > 0 && n != INVALID) {
alpar@389:           Value excess = (*_excess)[n];
alpar@389:           int new_level = _level->maxLevel();
alpar@389: 
alpar@389:           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Value rem = (*_capacity)[e] - (*_flow)[e];
alpar@389:             if (!_tolerance.positive(rem)) continue;
alpar@389:             Node v = _graph.target(e);
alpar@389:             if ((*_level)[v] < level) {
alpar@389:               if (!_level->active(v) && v != _target) {
alpar@389:                 _level->activate(v);
alpar@389:               }
alpar@389:               if (!_tolerance.less(rem, excess)) {
alpar@389:                 _flow->set(e, (*_flow)[e] + excess);
alpar@389:                 _excess->set(v, (*_excess)[v] + excess);
alpar@389:                 excess = 0;
alpar@389:                 goto no_more_push_1;
alpar@389:               } else {
alpar@389:                 excess -= rem;
alpar@389:                 _excess->set(v, (*_excess)[v] + rem);
alpar@389:                 _flow->set(e, (*_capacity)[e]);
alpar@389:               }
alpar@389:             } else if (new_level > (*_level)[v]) {
alpar@389:               new_level = (*_level)[v];
alpar@389:             }
alpar@389:           }
alpar@389: 
alpar@389:           for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Value rem = (*_flow)[e];
alpar@389:             if (!_tolerance.positive(rem)) continue;
alpar@389:             Node v = _graph.source(e);
alpar@389:             if ((*_level)[v] < level) {
alpar@389:               if (!_level->active(v) && v != _target) {
alpar@389:                 _level->activate(v);
alpar@389:               }
alpar@389:               if (!_tolerance.less(rem, excess)) {
alpar@389:                 _flow->set(e, (*_flow)[e] - excess);
alpar@389:                 _excess->set(v, (*_excess)[v] + excess);
alpar@389:                 excess = 0;
alpar@389:                 goto no_more_push_1;
alpar@389:               } else {
alpar@389:                 excess -= rem;
alpar@389:                 _excess->set(v, (*_excess)[v] + rem);
alpar@389:                 _flow->set(e, 0);
alpar@389:               }
alpar@389:             } else if (new_level > (*_level)[v]) {
alpar@389:               new_level = (*_level)[v];
alpar@389:             }
alpar@389:           }
alpar@389: 
alpar@389:         no_more_push_1:
alpar@389: 
alpar@389:           _excess->set(n, excess);
alpar@389: 
alpar@389:           if (excess != 0) {
alpar@389:             if (new_level + 1 < _level->maxLevel()) {
alpar@389:               _level->liftHighestActive(new_level + 1);
alpar@389:             } else {
alpar@389:               _level->liftHighestActiveToTop();
alpar@389:             }
alpar@389:             if (_level->emptyLevel(level)) {
alpar@389:               _level->liftToTop(level);
alpar@389:             }
alpar@389:           } else {
alpar@389:             _level->deactivate(n);
alpar@389:           }
alpar@389: 
alpar@389:           n = _level->highestActive();
alpar@389:           level = _level->highestActiveLevel();
alpar@389:           --num;
alpar@389:         }
alpar@389: 
alpar@389:         num = _node_num * 20;
alpar@389:         while (num > 0 && n != INVALID) {
alpar@389:           Value excess = (*_excess)[n];
alpar@389:           int new_level = _level->maxLevel();
alpar@389: 
alpar@389:           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Value rem = (*_capacity)[e] - (*_flow)[e];
alpar@389:             if (!_tolerance.positive(rem)) continue;
alpar@389:             Node v = _graph.target(e);
alpar@389:             if ((*_level)[v] < level) {
alpar@389:               if (!_level->active(v) && v != _target) {
alpar@389:                 _level->activate(v);
alpar@389:               }
alpar@389:               if (!_tolerance.less(rem, excess)) {
alpar@389:                 _flow->set(e, (*_flow)[e] + excess);
alpar@389:                 _excess->set(v, (*_excess)[v] + excess);
alpar@389:                 excess = 0;
alpar@389:                 goto no_more_push_2;
alpar@389:               } else {
alpar@389:                 excess -= rem;
alpar@389:                 _excess->set(v, (*_excess)[v] + rem);
alpar@389:                 _flow->set(e, (*_capacity)[e]);
alpar@389:               }
alpar@389:             } else if (new_level > (*_level)[v]) {
alpar@389:               new_level = (*_level)[v];
alpar@389:             }
alpar@389:           }
alpar@389: 
alpar@389:           for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Value rem = (*_flow)[e];
alpar@389:             if (!_tolerance.positive(rem)) continue;
alpar@389:             Node v = _graph.source(e);
alpar@389:             if ((*_level)[v] < level) {
alpar@389:               if (!_level->active(v) && v != _target) {
alpar@389:                 _level->activate(v);
alpar@389:               }
alpar@389:               if (!_tolerance.less(rem, excess)) {
alpar@389:                 _flow->set(e, (*_flow)[e] - excess);
alpar@389:                 _excess->set(v, (*_excess)[v] + excess);
alpar@389:                 excess = 0;
alpar@389:                 goto no_more_push_2;
alpar@389:               } else {
alpar@389:                 excess -= rem;
alpar@389:                 _excess->set(v, (*_excess)[v] + rem);
alpar@389:                 _flow->set(e, 0);
alpar@389:               }
alpar@389:             } else if (new_level > (*_level)[v]) {
alpar@389:               new_level = (*_level)[v];
alpar@389:             }
alpar@389:           }
alpar@389: 
alpar@389:         no_more_push_2:
alpar@389: 
alpar@389:           _excess->set(n, excess);
alpar@389: 
alpar@389:           if (excess != 0) {
alpar@389:             if (new_level + 1 < _level->maxLevel()) {
alpar@389:               _level->liftActiveOn(level, new_level + 1);
alpar@389:             } else {
alpar@389:               _level->liftActiveToTop(level);
alpar@389:             }
alpar@389:             if (_level->emptyLevel(level)) {
alpar@389:               _level->liftToTop(level);
alpar@389:             }
alpar@389:           } else {
alpar@389:             _level->deactivate(n);
alpar@389:           }
alpar@389: 
alpar@389:           while (level >= 0 && _level->activeFree(level)) {
alpar@389:             --level;
alpar@389:           }
alpar@389:           if (level == -1) {
alpar@389:             n = _level->highestActive();
alpar@389:             level = _level->highestActiveLevel();
alpar@389:           } else {
alpar@389:             n = _level->activeOn(level);
alpar@389:           }
alpar@389:           --num;
alpar@389:         }
alpar@389:       }
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Starts the second phase of the preflow algorithm.
alpar@389:     ///
alpar@389:     /// The preflow algorithm consists of two phases, this method runs
alpar@389:     /// the second phase. After calling \ref init() and \ref
alpar@389:     /// startFirstPhase() and then \ref startSecondPhase(), \ref
alpar@389:     /// flowMap() return a maximum flow, \ref flowValue() returns the
alpar@389:     /// value of a maximum flow, \ref minCut() returns a minimum cut
alpar@389:     /// \pre The \ref init() and startFirstPhase() functions should be
alpar@389:     /// called before.
alpar@389:     void startSecondPhase() {
alpar@389:       _phase = false;
alpar@389: 
alpar@389:       typename Digraph::template NodeMap<bool> reached(_graph);
alpar@389:       for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@389:         reached.set(n, (*_level)[n] < _level->maxLevel());
alpar@389:       }
alpar@389: 
alpar@389:       _level->initStart();
alpar@389:       _level->initAddItem(_source);
alpar@389: 
alpar@389:       std::vector<Node> queue;
alpar@389:       queue.push_back(_source);
alpar@389:       reached.set(_source, true);
alpar@389: 
alpar@389:       while (!queue.empty()) {
alpar@389:         _level->initNewLevel();
alpar@389:         std::vector<Node> nqueue;
alpar@389:         for (int i = 0; i < int(queue.size()); ++i) {
alpar@389:           Node n = queue[i];
alpar@389:           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Node v = _graph.target(e);
alpar@389:             if (!reached[v] && _tolerance.positive((*_flow)[e])) {
alpar@389:               reached.set(v, true);
alpar@389:               _level->initAddItem(v);
alpar@389:               nqueue.push_back(v);
alpar@389:             }
alpar@389:           }
alpar@389:           for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:             Node u = _graph.source(e);
alpar@389:             if (!reached[u] &&
alpar@389:                 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
alpar@389:               reached.set(u, true);
alpar@389:               _level->initAddItem(u);
alpar@389:               nqueue.push_back(u);
alpar@389:             }
alpar@389:           }
alpar@389:         }
alpar@389:         queue.swap(nqueue);
alpar@389:       }
alpar@389:       _level->initFinish();
alpar@389: 
alpar@389:       for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@389:         if (!reached[n]) {
alpar@389:           _level->dirtyTopButOne(n);
alpar@389:         } else if ((*_excess)[n] > 0 && _target != n) {
alpar@389:           _level->activate(n);
alpar@389:         }
alpar@389:       }
alpar@389: 
alpar@389:       Node n;
alpar@389:       while ((n = _level->highestActive()) != INVALID) {
alpar@389:         Value excess = (*_excess)[n];
alpar@389:         int level = _level->highestActiveLevel();
alpar@389:         int new_level = _level->maxLevel();
alpar@389: 
alpar@389:         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:           Value rem = (*_capacity)[e] - (*_flow)[e];
alpar@389:           if (!_tolerance.positive(rem)) continue;
alpar@389:           Node v = _graph.target(e);
alpar@389:           if ((*_level)[v] < level) {
alpar@389:             if (!_level->active(v) && v != _source) {
alpar@389:               _level->activate(v);
alpar@389:             }
alpar@389:             if (!_tolerance.less(rem, excess)) {
alpar@389:               _flow->set(e, (*_flow)[e] + excess);
alpar@389:               _excess->set(v, (*_excess)[v] + excess);
alpar@389:               excess = 0;
alpar@389:               goto no_more_push;
alpar@389:             } else {
alpar@389:               excess -= rem;
alpar@389:               _excess->set(v, (*_excess)[v] + rem);
alpar@389:               _flow->set(e, (*_capacity)[e]);
alpar@389:             }
alpar@389:           } else if (new_level > (*_level)[v]) {
alpar@389:             new_level = (*_level)[v];
alpar@389:           }
alpar@389:         }
alpar@389: 
alpar@389:         for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@389:           Value rem = (*_flow)[e];
alpar@389:           if (!_tolerance.positive(rem)) continue;
alpar@389:           Node v = _graph.source(e);
alpar@389:           if ((*_level)[v] < level) {
alpar@389:             if (!_level->active(v) && v != _source) {
alpar@389:               _level->activate(v);
alpar@389:             }
alpar@389:             if (!_tolerance.less(rem, excess)) {
alpar@389:               _flow->set(e, (*_flow)[e] - excess);
alpar@389:               _excess->set(v, (*_excess)[v] + excess);
alpar@389:               excess = 0;
alpar@389:               goto no_more_push;
alpar@389:             } else {
alpar@389:               excess -= rem;
alpar@389:               _excess->set(v, (*_excess)[v] + rem);
alpar@389:               _flow->set(e, 0);
alpar@389:             }
alpar@389:           } else if (new_level > (*_level)[v]) {
alpar@389:             new_level = (*_level)[v];
alpar@389:           }
alpar@389:         }
alpar@389: 
alpar@389:       no_more_push:
alpar@389: 
alpar@389:         _excess->set(n, excess);
alpar@389: 
alpar@389:         if (excess != 0) {
alpar@389:           if (new_level + 1 < _level->maxLevel()) {
alpar@389:             _level->liftHighestActive(new_level + 1);
alpar@389:           } else {
alpar@389:             // Calculation error
alpar@389:             _level->liftHighestActiveToTop();
alpar@389:           }
alpar@389:           if (_level->emptyLevel(level)) {
alpar@389:             // Calculation error
alpar@389:             _level->liftToTop(level);
alpar@389:           }
alpar@389:         } else {
alpar@389:           _level->deactivate(n);
alpar@389:         }
alpar@389: 
alpar@389:       }
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Runs the preflow algorithm.
alpar@389:     ///
alpar@389:     /// Runs the preflow algorithm.
alpar@389:     /// \note pf.run() is just a shortcut of the following code.
alpar@389:     /// \code
alpar@389:     ///   pf.init();
alpar@389:     ///   pf.startFirstPhase();
alpar@389:     ///   pf.startSecondPhase();
alpar@389:     /// \endcode
alpar@389:     void run() {
alpar@389:       init();
alpar@389:       startFirstPhase();
alpar@389:       startSecondPhase();
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Runs the preflow algorithm to compute the minimum cut.
alpar@389:     ///
alpar@389:     /// Runs the preflow algorithm to compute the minimum cut.
alpar@389:     /// \note pf.runMinCut() is just a shortcut of the following code.
alpar@389:     /// \code
alpar@389:     ///   pf.init();
alpar@389:     ///   pf.startFirstPhase();
alpar@389:     /// \endcode
alpar@389:     void runMinCut() {
alpar@389:       init();
alpar@389:       startFirstPhase();
alpar@389:     }
alpar@389: 
alpar@389:     /// @}
alpar@389: 
alpar@389:     /// \name Query Functions
alpar@389:     /// The result of the %Preflow algorithm can be obtained using these
alpar@389:     /// functions.\n
alpar@389:     /// Before the use of these functions,
alpar@389:     /// either run() or start() must be called.
alpar@389: 
alpar@389:     ///@{
alpar@389: 
alpar@389:     /// \brief Returns the value of the maximum flow.
alpar@389:     ///
alpar@389:     /// Returns the value of the maximum flow by returning the excess
alpar@389:     /// of the target node \c t. This value equals to the value of
alpar@389:     /// the maximum flow already after the first phase.
alpar@389:     Value flowValue() const {
alpar@389:       return (*_excess)[_target];
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Returns true when the node is on the source side of minimum cut.
alpar@389:     ///
alpar@389:     /// Returns true when the node is on the source side of minimum
alpar@389:     /// cut. This method can be called both after running \ref
alpar@389:     /// startFirstPhase() and \ref startSecondPhase().
alpar@389:     bool minCut(const Node& node) const {
alpar@389:       return ((*_level)[node] == _level->maxLevel()) == _phase;
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Returns a minimum value cut.
alpar@389:     ///
alpar@389:     /// Sets the \c cutMap to the characteristic vector of a minimum value
alpar@389:     /// cut. This method can be called both after running \ref
alpar@389:     /// startFirstPhase() and \ref startSecondPhase(). The result after second
alpar@389:     /// phase could be changed slightly if inexact computation is used.
alpar@389:     /// \pre The \c cutMap should be a bool-valued node-map.
alpar@389:     template <typename CutMap>
alpar@389:     void minCutMap(CutMap& cutMap) const {
alpar@389:       for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@389:         cutMap.set(n, minCut(n));
alpar@389:       }
alpar@389:     }
alpar@389: 
alpar@389:     /// \brief Returns the flow on the arc.
alpar@389:     ///
alpar@389:     /// Sets the \c flowMap to the flow on the arcs. This method can
alpar@389:     /// be called after the second phase of algorithm.
alpar@389:     Value flow(const Arc& arc) const {
alpar@389:       return (*_flow)[arc];
alpar@389:     }
alpar@389: 
alpar@389:     /// @}
alpar@389:   };
alpar@389: }
alpar@389: 
alpar@389: #endif