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