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@440:  * Copyright (C) 2003-2009
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.
kpeter@393:   /// \tparam _Digraph Digraph type.
kpeter@393:   /// \tparam _CapacityMap Capacity map type.
kpeter@393:   template <typename _Digraph, typename _CapacityMap>
alpar@389:   struct PreflowDefaultTraits {
alpar@389: 
kpeter@393:     /// \brief The type of the digraph the algorithm runs on.
kpeter@393:     typedef _Digraph 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: 
kpeter@393:     /// \brief The type of the flow values.
alpar@389:     typedef typename CapacityMap::Value Value;
alpar@389: 
kpeter@393:     /// \brief The type of the map that stores the flow values.
alpar@389:     ///
kpeter@393:     /// The type of the map 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: 
kpeter@393:     /// \brief The elevator 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:     ///
kpeter@393:     /// This function instantiates an \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:   ///
kpeter@393:   /// \brief %Preflow algorithm class.
alpar@389:   ///
kpeter@393:   /// This class provides an implementation of Goldberg-Tarjan's \e preflow
kpeter@393:   /// \e push-relabel algorithm producing a flow of maximum value in a
kpeter@393:   /// digraph. The preflow algorithms are the fastest known maximum
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:   ///
kpeter@393:   /// The algorithm consists of two phases. After the first phase
kpeter@393:   /// the maximum flow value and the minimum cut is obtained. The
kpeter@393:   /// second phase constructs a feasible maximum flow on each arc.
alpar@389:   ///
kpeter@393:   /// \tparam _Digraph The type of the digraph the algorithm runs on.
kpeter@393:   /// \tparam _CapacityMap The type of the capacity map. The default map
kpeter@393:   /// type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
alpar@389: #ifdef DOXYGEN
kpeter@393:   template <typename _Digraph, typename _CapacityMap, typename _Traits>
alpar@389: #else
kpeter@393:   template <typename _Digraph,
kpeter@393:             typename _CapacityMap = typename _Digraph::template ArcMap<int>,
kpeter@393:             typename _Traits = PreflowDefaultTraits<_Digraph, _CapacityMap> >
alpar@389: #endif
alpar@389:   class Preflow {
alpar@389:   public:
alpar@389: 
kpeter@393:     ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
alpar@389:     typedef _Traits Traits;
kpeter@393:     ///The type of the digraph the algorithm runs on.
alpar@389:     typedef typename Traits::Digraph Digraph;
kpeter@393:     ///The type of the capacity map.
alpar@389:     typedef typename Traits::CapacityMap CapacityMap;
kpeter@393:     ///The type of the flow values.
alpar@389:     typedef typename Traits::Value Value;
alpar@389: 
kpeter@393:     ///The type of the flow map.
alpar@389:     typedef typename Traits::FlowMap FlowMap;
kpeter@393:     ///The type of the elevator.
alpar@389:     typedef typename Traits::Elevator Elevator;
kpeter@393:     ///The type of the tolerance.
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: 
kpeter@393:     ///\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
kpeter@393:     /// 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
kpeter@393:     /// type. If this named parameter is used, then an external
kpeter@393:     /// elevator object must be passed to the algorithm using the
kpeter@393:     /// \ref elevator(Elevator&) "elevator()" function before calling
kpeter@393:     /// \ref run() or \ref init().
kpeter@393:     /// \sa SetStandardElevator
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
kpeter@393:     /// Elevator type with automatic allocation
alpar@389:     ///
alpar@389:     /// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@393:     /// type with automatic allocation.
kpeter@393:     /// The Elevator should have standard constructor interface to be
kpeter@393:     /// able to automatically created by the algorithm (i.e. the
kpeter@393:     /// digraph and the maximum level should be passed to it).
kpeter@393:     /// However an external elevator object could also be passed to the
kpeter@393:     /// algorithm with the \ref elevator(Elevator&) "elevator()" function
kpeter@393:     /// before calling \ref run() or \ref init().
kpeter@393:     /// \sa SetElevator
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,
kpeter@393:             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: 
kpeter@393:     /// \brief Destructor.
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.
kpeter@393:     /// \return <tt>(*this)</tt>
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.
kpeter@393:     /// If you don't use this function before calling \ref run() or
kpeter@393:     /// \ref init(), an instance will be allocated automatically.
kpeter@393:     /// The destructor deallocates this automatically allocated map,
kpeter@393:     /// of course.
kpeter@393:     /// \return <tt>(*this)</tt>
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: 
kpeter@393:     /// \brief Sets the source node.
alpar@389:     ///
kpeter@393:     /// Sets the source node.
kpeter@393:     /// \return <tt>(*this)</tt>
kpeter@393:     Preflow& source(const Node& node) {
kpeter@393:       _source = node;
kpeter@393:       return *this;
alpar@389:     }
alpar@389: 
kpeter@393:     /// \brief Sets the target node.
alpar@389:     ///
kpeter@393:     /// Sets the target node.
kpeter@393:     /// \return <tt>(*this)</tt>
kpeter@393:     Preflow& target(const Node& node) {
kpeter@393:       _target = node;
kpeter@393:       return *this;
kpeter@393:     }
kpeter@393: 
kpeter@393:     /// \brief Sets the elevator used by algorithm.
kpeter@393:     ///
kpeter@393:     /// Sets the elevator used by algorithm.
kpeter@393:     /// If you don't use this function before calling \ref run() or
kpeter@393:     /// \ref init(), an instance will be allocated automatically.
kpeter@393:     /// The destructor deallocates this automatically allocated elevator,
kpeter@393:     /// of course.
kpeter@393:     /// \return <tt>(*this)</tt>
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: 
kpeter@393:     /// \brief Returns a const reference to the elevator.
alpar@389:     ///
kpeter@393:     /// Returns a const reference to the elevator.
kpeter@393:     ///
kpeter@393:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@393:     /// using this function.
kpeter@420:     const Elevator& elevator() const {
alpar@389:       return *_level;
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: 
kpeter@393:     /// \brief Returns a const reference to the tolerance.
alpar@389:     ///
kpeter@393:     /// Returns a const reference to the tolerance.
alpar@389:     const Tolerance& tolerance() const {
alpar@389:       return tolerance;
alpar@389:     }
alpar@389: 
kpeter@393:     /// \name Execution Control
kpeter@393:     /// The simplest way to execute the preflow algorithm is to use
kpeter@393:     /// \ref run() or \ref runMinCut().\n
kpeter@393:     /// If you need more control on the initial solution or the execution,
kpeter@393:     /// first you have to call one of the \ref init() functions, then
kpeter@393:     /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
alpar@389: 
alpar@389:     ///@{
alpar@389: 
alpar@389:     /// \brief Initializes the internal data structures.
alpar@389:     ///
kpeter@393:     /// Initializes the internal data structures and sets the initial
kpeter@393:     /// flow to zero on each arc.
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: 
kpeter@393:     /// \brief Initializes the internal data structures using the
kpeter@393:     /// given flow map.
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
kpeter@393:     /// flow or at least a preflow, i.e. at each node excluding the
kpeter@393:     /// source node the incoming flow should greater or equal to the
alpar@389:     /// outgoing flow.
kpeter@393:     /// \return \c false if 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.
kpeter@393:     /// \pre One of the \ref init() functions must be called before
kpeter@393:     /// using this function.
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
kpeter@393:     /// the second phase. After calling one of the \ref init() functions
kpeter@393:     /// and \ref startFirstPhase() and then \ref startSecondPhase(),
kpeter@393:     /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
alpar@389:     /// value of a maximum flow, \ref minCut() returns a minimum cut
kpeter@393:     /// \pre One of the \ref init() functions and \ref startFirstPhase()
kpeter@393:     /// must be called before using this function.
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
kpeter@393:     /// The results of the preflow algorithm can be obtained using these
alpar@389:     /// functions.\n
kpeter@393:     /// Either one of the \ref run() "run*()" functions or one of the
kpeter@393:     /// \ref startFirstPhase() "start*()" functions should be called
kpeter@393:     /// before using them.
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
kpeter@393:     /// of the target node. This value equals to the value of
kpeter@393:     /// the maximum flow already after the first phase of the algorithm.
kpeter@393:     ///
kpeter@393:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@393:     /// using this function.
alpar@389:     Value flowValue() const {
alpar@389:       return (*_excess)[_target];
alpar@389:     }
alpar@389: 
kpeter@393:     /// \brief Returns the flow on the given arc.
alpar@389:     ///
kpeter@393:     /// Returns the flow on the given arc. This method can
kpeter@393:     /// be called after the second phase of the algorithm.
kpeter@393:     ///
kpeter@393:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@393:     /// using this function.
kpeter@393:     Value flow(const Arc& arc) const {
kpeter@393:       return (*_flow)[arc];
kpeter@393:     }
kpeter@393: 
kpeter@393:     /// \brief Returns a const reference to the flow map.
kpeter@393:     ///
kpeter@393:     /// Returns a const reference to the arc map storing the found flow.
kpeter@393:     /// This method can be called after the second phase of the algorithm.
kpeter@393:     ///
kpeter@393:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@393:     /// using this function.
kpeter@420:     const FlowMap& flowMap() const {
kpeter@393:       return *_flow;
kpeter@393:     }
kpeter@393: 
kpeter@393:     /// \brief Returns \c true when the node is on the source side of the
kpeter@393:     /// minimum cut.
kpeter@393:     ///
kpeter@393:     /// Returns true when the node is on the source side of the found
kpeter@393:     /// minimum cut. This method can be called both after running \ref
alpar@389:     /// startFirstPhase() and \ref startSecondPhase().
kpeter@393:     ///
kpeter@393:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@393:     /// using this function.
alpar@389:     bool minCut(const Node& node) const {
alpar@389:       return ((*_level)[node] == _level->maxLevel()) == _phase;
alpar@389:     }
alpar@389: 
kpeter@393:     /// \brief Gives back a minimum value cut.
alpar@389:     ///
kpeter@393:     /// Sets \c cutMap to the characteristic vector of a minimum value
kpeter@393:     /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
kpeter@393:     /// node map with \c bool (or convertible) value type.
kpeter@393:     ///
kpeter@393:     /// This method can be called both after running \ref startFirstPhase()
kpeter@393:     /// and \ref startSecondPhase(). The result after the second phase
kpeter@393:     /// could be slightly different if inexact computation is used.
kpeter@393:     ///
kpeter@393:     /// \note This function calls \ref minCut() for each node, so it runs in
kpeter@393:     /// \f$O(n)\f$ time.
kpeter@393:     ///
kpeter@393:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@393:     /// using this function.
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:     /// @}
alpar@389:   };
alpar@389: }
alpar@389: 
alpar@389: #endif