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 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. alpar@389: /// \param _Graph Digraph type. alpar@389: /// \param _CapacityMap Type of capacity map. alpar@389: template 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 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 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 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 alpar@389: #else alpar@389: template , 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 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 alpar@389: struct DefFlowMapTraits : 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 alpar@389: struct DefFlowMap alpar@389: : public Preflow > { alpar@389: typedef Preflow > Create; alpar@389: }; alpar@389: alpar@389: template alpar@389: struct DefElevatorTraits : 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 alpar@389: struct DefElevator alpar@389: : public Preflow > { alpar@389: typedef Preflow > Create; alpar@389: }; alpar@389: alpar@389: template alpar@389: struct DefStandardElevatorTraits : 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 alpar@389: struct DefStandardElevator 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, 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 = ↦ 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 = ↦ 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 reached(_graph, false); alpar@389: alpar@389: _level->initStart(); alpar@389: _level->initAddItem(_target); alpar@389: alpar@389: std::vector 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 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 alpar@389: bool flowInit(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 reached(_graph, false); alpar@389: alpar@389: _level->initStart(); alpar@389: _level->initAddItem(_target); alpar@389: alpar@389: std::vector 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 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 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 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 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 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