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@927: * Copyright (C) 2003-2011 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 alpar@389: #include 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@492: /// \tparam GR Digraph type. kpeter@559: /// \tparam CAP Capacity map type. kpeter@559: template alpar@389: struct PreflowDefaultTraits { alpar@389: kpeter@393: /// \brief The type of the digraph the algorithm runs on. kpeter@492: typedef GR 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. kpeter@559: typedef CAP CapacityMap; alpar@389: kpeter@393: /// \brief The type of the flow values. kpeter@641: 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. kpeter@713: #ifdef DOXYGEN kpeter@713: typedef GR::ArcMap FlowMap; kpeter@713: #else kpeter@641: typedef typename Digraph::template ArcMap FlowMap; kpeter@713: #endif alpar@389: alpar@389: /// \brief Instantiates a FlowMap. alpar@389: /// alpar@389: /// This function instantiates a \ref FlowMap. kpeter@610: /// \param digraph The digraph for 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: /// kpeter@713: /// \sa Elevator, LinkedElevator kpeter@713: #ifdef DOXYGEN kpeter@713: typedef lemon::Elevator Elevator; kpeter@713: #else kpeter@713: typedef lemon::Elevator Elevator; kpeter@713: #endif alpar@389: alpar@389: /// \brief Instantiates an Elevator. alpar@389: /// kpeter@393: /// This function instantiates an \ref Elevator. kpeter@610: /// \param digraph The digraph for 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. kpeter@641: typedef lemon::Tolerance 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@559: /// \e push-relabel algorithm producing a \ref max_flow kpeter@755: /// "flow of maximum value" in a digraph \ref clrs01algorithms, kpeter@755: /// \ref amo93networkflows, \ref goldberg88newapproach. kpeter@559: /// The preflow algorithms are the fastest known maximum kpeter@689: /// flow algorithms. The current implementation uses 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@823: /// \warning This implementation cannot handle infinite or very large kpeter@823: /// capacities (e.g. the maximum value of \c CAP::Value). kpeter@823: /// kpeter@492: /// \tparam GR The type of the digraph the algorithm runs on. kpeter@559: /// \tparam CAP The type of the capacity map. The default map kpeter@492: /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap". kpeter@825: /// \tparam TR The traits class that defines various types used by the kpeter@825: /// algorithm. By default, it is \ref PreflowDefaultTraits kpeter@825: /// "PreflowDefaultTraits". kpeter@825: /// In most cases, this parameter should not be set directly, kpeter@825: /// consider to use the named template parameters instead. alpar@389: #ifdef DOXYGEN kpeter@559: template alpar@389: #else kpeter@492: template , kpeter@559: typename TR = PreflowDefaultTraits > alpar@389: #endif alpar@389: class Preflow { alpar@389: public: alpar@389: kpeter@393: ///The \ref PreflowDefaultTraits "traits class" of the algorithm. kpeter@492: typedef TR 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. kpeter@641: 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: kpeter@641: typedef typename Digraph::template NodeMap 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: kpeter@559: template alpar@391: struct SetFlowMapTraits : public Traits { kpeter@559: typedef T 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. kpeter@559: template alpar@391: struct SetFlowMap kpeter@559: : public Preflow > { alpar@389: typedef Preflow > Create; alpar@389: }; alpar@389: kpeter@559: template alpar@391: struct SetElevatorTraits : public Traits { kpeter@559: typedef T 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 kpeter@559: template alpar@391: struct SetElevator kpeter@559: : public Preflow > { alpar@389: typedef Preflow > Create; alpar@389: }; alpar@389: kpeter@559: template alpar@391: struct SetStandardElevatorTraits : public Traits { kpeter@559: typedef T 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@786: /// 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 kpeter@559: template alpar@391: struct SetStandardElevator alpar@389: : public Preflow > { alpar@389: typedef Preflow > 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 (*this) alpar@389: Preflow& capacityMap(const CapacityMap& map) { alpar@389: _capacity = ↦ 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 (*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 = ↦ 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 (*this) 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 (*this) 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 (*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: 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: kpeter@689: /// \brief Sets the tolerance used by the algorithm. alpar@389: /// kpeter@689: /// Sets the tolerance object used by the algorithm. kpeter@689: /// \return (*this) kpeter@688: Preflow& tolerance(const Tolerance& tolerance) { 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@689: /// Returns a const reference to the tolerance object used by kpeter@689: /// the algorithm. alpar@389: const Tolerance& tolerance() const { kpeter@688: 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@713: /// If you need better control on the initial solution or the execution, kpeter@713: /// you have to call one of the \ref init() functions first, 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) { kpeter@581: (*_excess)[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 reached(_graph, false); alpar@389: alpar@389: _level->initStart(); alpar@389: _level->initAddItem(_target); alpar@389: alpar@389: std::vector queue; kpeter@581: reached[_source] = true; alpar@389: alpar@389: queue.push_back(_target); kpeter@581: reached[_target] = true; alpar@389: while (!queue.empty()) { alpar@389: _level->initNewLevel(); alpar@389: std::vector 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])) { kpeter@581: reached[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]); kpeter@581: (*_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 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) { kpeter@641: 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; kpeter@581: (*_excess)[n] = excess; alpar@389: } alpar@389: alpar@389: typename Digraph::template NodeMap reached(_graph, false); alpar@389: alpar@389: _level->initStart(); alpar@389: _level->initAddItem(_target); alpar@389: alpar@389: std::vector queue; kpeter@581: reached[_source] = true; alpar@389: alpar@389: queue.push_back(_target); kpeter@581: reached[_target] = true; alpar@389: while (!queue.empty()) { alpar@389: _level->initNewLevel(); alpar@389: std::vector 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])) { kpeter@581: reached[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])) { kpeter@581: reached[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) { kpeter@641: 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]); kpeter@581: (*_excess)[u] += rem; alpar@389: } alpar@389: } alpar@389: for (InArcIt e(_graph, _source); e != INVALID; ++e) { kpeter@641: 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); kpeter@581: (*_excess)[v] += rem; alpar@389: } alpar@389: } alpar@927: for (NodeIt n(_graph); n != INVALID; ++n) alpar@901: if(n!=_source && n!=_target && _tolerance.positive((*_excess)[n])) alpar@901: _level->activate(n); alpar@927: 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: deba@894: while (true) { alpar@389: int num = _node_num; alpar@389: deba@894: Node n = INVALID; deba@894: int level = -1; deba@894: deba@894: while (num > 0) { deba@894: n = _level->highestActive(); deba@894: if (n == INVALID) goto first_phase_done; deba@894: level = _level->highestActiveLevel(); deba@894: --num; alpar@927: kpeter@641: Value excess = (*_excess)[n]; alpar@389: int new_level = _level->maxLevel(); alpar@389: alpar@389: for (OutArcIt e(_graph, n); e != INVALID; ++e) { kpeter@641: 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); kpeter@581: (*_excess)[v] += excess; alpar@389: excess = 0; alpar@389: goto no_more_push_1; alpar@389: } else { alpar@389: excess -= rem; kpeter@581: (*_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) { kpeter@641: 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); kpeter@581: (*_excess)[v] += excess; alpar@389: excess = 0; alpar@389: goto no_more_push_1; alpar@389: } else { alpar@389: excess -= rem; kpeter@581: (*_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: kpeter@581: (*_excess)[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: alpar@389: num = _node_num * 20; deba@894: while (num > 0) { deba@894: while (level >= 0 && _level->activeFree(level)) { deba@894: --level; deba@894: } deba@894: if (level == -1) { deba@894: n = _level->highestActive(); deba@894: level = _level->highestActiveLevel(); deba@894: if (n == INVALID) goto first_phase_done; deba@894: } else { deba@894: n = _level->activeOn(level); deba@894: } deba@894: --num; deba@894: kpeter@641: Value excess = (*_excess)[n]; alpar@389: int new_level = _level->maxLevel(); alpar@389: alpar@389: for (OutArcIt e(_graph, n); e != INVALID; ++e) { kpeter@641: 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); kpeter@581: (*_excess)[v] += excess; alpar@389: excess = 0; alpar@389: goto no_more_push_2; alpar@389: } else { alpar@389: excess -= rem; kpeter@581: (*_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) { kpeter@641: 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); kpeter@581: (*_excess)[v] += excess; alpar@389: excess = 0; alpar@389: goto no_more_push_2; alpar@389: } else { alpar@389: excess -= rem; kpeter@581: (*_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: kpeter@581: (*_excess)[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: } deba@894: first_phase_done:; 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 reached(_graph); alpar@389: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@581: reached[n] = (*_level)[n] < _level->maxLevel(); alpar@389: } alpar@389: alpar@389: _level->initStart(); alpar@389: _level->initAddItem(_source); alpar@389: alpar@389: std::vector queue; alpar@389: queue.push_back(_source); kpeter@581: reached[_source] = true; alpar@389: alpar@389: while (!queue.empty()) { alpar@389: _level->initNewLevel(); alpar@389: std::vector 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])) { kpeter@581: reached[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])) { kpeter@581: reached[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) { kpeter@641: 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) { kpeter@641: 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); kpeter@581: (*_excess)[v] += excess; alpar@389: excess = 0; alpar@389: goto no_more_push; alpar@389: } else { alpar@389: excess -= rem; kpeter@581: (*_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) { kpeter@641: 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); kpeter@581: (*_excess)[v] += excess; alpar@389: excess = 0; alpar@389: goto no_more_push; alpar@389: } else { alpar@389: excess -= rem; kpeter@581: (*_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: kpeter@581: (*_excess)[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. kpeter@641: Value flowValue() const { alpar@389: return (*_excess)[_target]; alpar@389: } alpar@389: kpeter@641: /// \brief Returns the flow value on the given arc. alpar@389: /// kpeter@641: /// Returns the flow value 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@641: 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@559: /// O(n) time. kpeter@393: /// kpeter@393: /// \pre Either \ref run() or \ref init() must be called before kpeter@393: /// using this function. alpar@389: template 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