alpar@399: /* -*- mode: C++; indent-tabs-mode: nil; -*- alpar@399: * alpar@399: * This file is a part of LEMON, a generic C++ optimization library. alpar@399: * alpar@440: * Copyright (C) 2003-2009 alpar@399: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport alpar@399: * (Egervary Research Group on Combinatorial Optimization, EGRES). alpar@399: * alpar@399: * Permission to use, modify and distribute this software is granted alpar@399: * provided that this copyright notice appears in all copies. For alpar@399: * precise terms see the accompanying LICENSE file. alpar@399: * alpar@399: * This software is provided "AS IS" with no warranty of any kind, alpar@399: * express or implied, and with no claim as to its suitability for any alpar@399: * purpose. alpar@399: * alpar@399: */ alpar@399: alpar@399: #ifndef LEMON_CIRCULATION_H alpar@399: #define LEMON_CIRCULATION_H alpar@399: alpar@399: #include alpar@399: #include kpeter@622: #include alpar@399: alpar@399: ///\ingroup max_flow alpar@399: ///\file kpeter@402: ///\brief Push-relabel algorithm for finding a feasible circulation. alpar@399: /// alpar@399: namespace lemon { alpar@399: alpar@399: /// \brief Default traits class of Circulation class. alpar@399: /// alpar@399: /// Default traits class of Circulation class. kpeter@610: /// kpeter@610: /// \tparam GR Type of the digraph the algorithm runs on. kpeter@610: /// \tparam LM The type of the lower bound map. kpeter@610: /// \tparam UM The type of the upper bound (capacity) map. kpeter@610: /// \tparam SM The type of the supply map. kpeter@492: template alpar@399: struct CirculationDefaultTraits { alpar@399: kpeter@402: /// \brief The type of the digraph the algorithm runs on. kpeter@492: typedef GR Digraph; alpar@399: kpeter@610: /// \brief The type of the lower bound map. alpar@399: /// kpeter@610: /// The type of the map that stores the lower bounds on the arcs. kpeter@610: /// It must conform to the \ref concepts::ReadMap "ReadMap" concept. kpeter@610: typedef LM LowerMap; alpar@399: kpeter@610: /// \brief The type of the upper bound (capacity) map. alpar@399: /// kpeter@610: /// The type of the map that stores the upper bounds (capacities) kpeter@610: /// on the arcs. kpeter@610: /// It must conform to the \ref concepts::ReadMap "ReadMap" concept. kpeter@610: typedef UM UpperMap; alpar@399: kpeter@610: /// \brief The type of supply map. alpar@399: /// kpeter@610: /// The type of the map that stores the signed supply values of the kpeter@610: /// nodes. kpeter@610: /// It must conform to the \ref concepts::ReadMap "ReadMap" concept. kpeter@610: typedef SM SupplyMap; alpar@399: kpeter@641: /// \brief The type of the flow and supply values. kpeter@641: typedef typename SupplyMap::Value Value; alpar@399: kpeter@402: /// \brief The type of the map that stores the flow values. alpar@399: /// kpeter@402: /// The type of the map that stores the flow values. kpeter@610: /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" kpeter@610: /// 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@399: alpar@399: /// \brief Instantiates a FlowMap. alpar@399: /// alpar@399: /// This function instantiates a \ref FlowMap. kpeter@610: /// \param digraph The digraph for which we would like to define alpar@399: /// the flow map. alpar@399: static FlowMap* createFlowMap(const Digraph& digraph) { alpar@399: return new FlowMap(digraph); alpar@399: } alpar@399: kpeter@402: /// \brief The elevator type used by the algorithm. alpar@399: /// kpeter@402: /// The elevator type used by the algorithm. alpar@399: /// kpeter@713: /// \sa Elevator, LinkedElevator kpeter@713: #ifdef DOXYGEN kpeter@713: typedef lemon::Elevator Elevator; kpeter@713: #else alpar@399: typedef lemon::Elevator Elevator; kpeter@713: #endif alpar@399: alpar@399: /// \brief Instantiates an Elevator. alpar@399: /// kpeter@402: /// This function instantiates an \ref Elevator. kpeter@610: /// \param digraph The digraph for which we would like to define alpar@399: /// the elevator. alpar@399: /// \param max_level The maximum level of the elevator. alpar@399: static Elevator* createElevator(const Digraph& digraph, int max_level) { alpar@399: return new Elevator(digraph, max_level); alpar@399: } alpar@399: alpar@399: /// \brief The tolerance used by the algorithm alpar@399: /// alpar@399: /// The tolerance used by the algorithm to handle inexact computation. kpeter@641: typedef lemon::Tolerance Tolerance; alpar@399: alpar@399: }; alpar@399: kpeter@402: /** kpeter@402: \brief Push-relabel algorithm for the network circulation problem. alpar@399: alpar@399: \ingroup max_flow kpeter@610: This class implements a push-relabel algorithm for the \e network kpeter@610: \e circulation problem. kpeter@402: It is to find a feasible circulation when lower and upper bounds kpeter@610: are given for the flow values on the arcs and lower bounds are kpeter@610: given for the difference between the outgoing and incoming flow kpeter@610: at the nodes. kpeter@402: alpar@399: The exact formulation of this problem is the following. kpeter@622: Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$ kpeter@622: \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and kpeter@622: upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$ kpeter@610: holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$ kpeter@610: denotes the signed supply values of the nodes. kpeter@610: If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ kpeter@610: supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with kpeter@610: \f$-sup(u)\f$ demand. kpeter@622: A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$ kpeter@610: solution of the following problem. kpeter@402: kpeter@610: \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) kpeter@610: \geq sup(u) \quad \forall u\in V, \f] kpeter@610: \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f] kpeter@610: kpeter@610: The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be kpeter@610: zero or negative in order to have a feasible solution (since the sum kpeter@610: of the expressions on the left-hand side of the inequalities is zero). kpeter@610: It means that the total demand must be greater or equal to the total kpeter@610: supply and all the supplies have to be carried out from the supply nodes, kpeter@610: but there could be demands that are not satisfied. kpeter@610: If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand kpeter@610: constraints have to be satisfied with equality, i.e. all demands kpeter@610: have to be satisfied and all supplies have to be used. kpeter@610: kpeter@610: If you need the opposite inequalities in the supply/demand constraints kpeter@610: (i.e. the total demand is less than the total supply and all the demands kpeter@610: have to be satisfied while there could be supplies that are not used), kpeter@610: then you could easily transform the problem to the above form by reversing kpeter@610: the direction of the arcs and taking the negative of the supply values kpeter@610: (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). kpeter@610: kpeter@622: This algorithm either calculates a feasible circulation, or provides kpeter@622: a \ref barrier() "barrier", which prooves that a feasible soultion kpeter@622: cannot exist. kpeter@622: kpeter@610: Note that this algorithm also provides a feasible solution for the kpeter@610: \ref min_cost_flow "minimum cost flow problem". kpeter@402: kpeter@492: \tparam GR The type of the digraph the algorithm runs on. kpeter@610: \tparam LM The type of the lower bound map. The default kpeter@492: map type is \ref concepts::Digraph::ArcMap "GR::ArcMap". kpeter@610: \tparam UM The type of the upper bound (capacity) map. kpeter@610: The default map type is \c LM. kpeter@610: \tparam SM The type of the supply map. The default map type is kpeter@492: \ref concepts::Digraph::NodeMap "GR::NodeMap". kpeter@825: \tparam TR The traits class that defines various types used by the kpeter@825: algorithm. By default, it is \ref CirculationDefaultTraits kpeter@825: "CirculationDefaultTraits". kpeter@825: In most cases, this parameter should not be set directly, kpeter@825: consider to use the named template parameters instead. alpar@399: */ kpeter@402: #ifdef DOXYGEN kpeter@492: template< typename GR, kpeter@492: typename LM, kpeter@492: typename UM, kpeter@610: typename SM, kpeter@492: typename TR > kpeter@402: #else kpeter@492: template< typename GR, kpeter@492: typename LM = typename GR::template ArcMap, kpeter@492: typename UM = LM, kpeter@610: typename SM = typename GR::template NodeMap, kpeter@610: typename TR = CirculationDefaultTraits > kpeter@402: #endif alpar@399: class Circulation { kpeter@402: public: alpar@399: kpeter@402: ///The \ref CirculationDefaultTraits "traits class" of the algorithm. kpeter@492: typedef TR Traits; kpeter@402: ///The type of the digraph the algorithm runs on. alpar@399: typedef typename Traits::Digraph Digraph; kpeter@641: ///The type of the flow and supply values. kpeter@641: typedef typename Traits::Value Value; alpar@399: kpeter@610: ///The type of the lower bound map. kpeter@610: typedef typename Traits::LowerMap LowerMap; kpeter@610: ///The type of the upper bound (capacity) map. kpeter@610: typedef typename Traits::UpperMap UpperMap; kpeter@610: ///The type of the supply map. kpeter@610: typedef typename Traits::SupplyMap SupplyMap; kpeter@402: ///The type of the flow map. alpar@399: typedef typename Traits::FlowMap FlowMap; kpeter@402: kpeter@402: ///The type of the elevator. alpar@399: typedef typename Traits::Elevator Elevator; kpeter@402: ///The type of the tolerance. alpar@399: typedef typename Traits::Tolerance Tolerance; alpar@399: kpeter@402: private: kpeter@402: kpeter@402: TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); alpar@399: alpar@399: const Digraph &_g; alpar@399: int _node_num; alpar@399: kpeter@610: const LowerMap *_lo; kpeter@610: const UpperMap *_up; kpeter@610: const SupplyMap *_supply; alpar@399: alpar@399: FlowMap *_flow; alpar@399: bool _local_flow; alpar@399: alpar@399: Elevator* _level; alpar@399: bool _local_level; alpar@399: kpeter@641: typedef typename Digraph::template NodeMap ExcessMap; alpar@399: ExcessMap* _excess; alpar@399: alpar@399: Tolerance _tol; alpar@399: int _el; alpar@399: alpar@399: public: alpar@399: alpar@399: typedef Circulation Create; alpar@399: kpeter@402: ///\name Named Template Parameters alpar@399: alpar@399: ///@{ alpar@399: kpeter@559: template alpar@401: struct SetFlowMapTraits : public Traits { kpeter@559: typedef T FlowMap; alpar@399: static FlowMap *createFlowMap(const Digraph&) { alpar@399: LEMON_ASSERT(false, "FlowMap is not initialized"); alpar@399: return 0; // ignore warnings alpar@399: } alpar@399: }; alpar@399: alpar@399: /// \brief \ref named-templ-param "Named parameter" for setting alpar@399: /// FlowMap type alpar@399: /// alpar@399: /// \ref named-templ-param "Named parameter" for setting FlowMap kpeter@402: /// type. kpeter@559: template alpar@401: struct SetFlowMap kpeter@610: : public Circulation > { kpeter@610: typedef Circulation > Create; alpar@399: }; alpar@399: kpeter@559: template alpar@401: struct SetElevatorTraits : public Traits { kpeter@559: typedef T Elevator; alpar@399: static Elevator *createElevator(const Digraph&, int) { alpar@399: LEMON_ASSERT(false, "Elevator is not initialized"); alpar@399: return 0; // ignore warnings alpar@399: } alpar@399: }; alpar@399: alpar@399: /// \brief \ref named-templ-param "Named parameter" for setting alpar@399: /// Elevator type alpar@399: /// alpar@399: /// \ref named-templ-param "Named parameter" for setting Elevator kpeter@402: /// type. If this named parameter is used, then an external kpeter@402: /// elevator object must be passed to the algorithm using the kpeter@402: /// \ref elevator(Elevator&) "elevator()" function before calling kpeter@402: /// \ref run() or \ref init(). kpeter@402: /// \sa SetStandardElevator kpeter@559: template alpar@401: struct SetElevator kpeter@610: : public Circulation > { kpeter@610: typedef Circulation > Create; alpar@399: }; alpar@399: kpeter@559: template alpar@401: struct SetStandardElevatorTraits : public Traits { kpeter@559: typedef T Elevator; alpar@399: static Elevator *createElevator(const Digraph& digraph, int max_level) { alpar@399: return new Elevator(digraph, max_level); alpar@399: } alpar@399: }; alpar@399: alpar@399: /// \brief \ref named-templ-param "Named parameter" for setting kpeter@402: /// Elevator type with automatic allocation alpar@399: /// alpar@399: /// \ref named-templ-param "Named parameter" for setting Elevator kpeter@402: /// type with automatic allocation. kpeter@402: /// The Elevator should have standard constructor interface to be kpeter@402: /// able to automatically created by the algorithm (i.e. the kpeter@402: /// digraph and the maximum level should be passed to it). kpeter@786: /// However, an external elevator object could also be passed to the kpeter@402: /// algorithm with the \ref elevator(Elevator&) "elevator()" function kpeter@402: /// before calling \ref run() or \ref init(). kpeter@402: /// \sa SetElevator kpeter@559: template alpar@401: struct SetStandardElevator kpeter@610: : public Circulation > { kpeter@610: typedef Circulation > Create; alpar@399: }; alpar@399: alpar@399: /// @} alpar@399: alpar@399: protected: alpar@399: alpar@399: Circulation() {} alpar@399: alpar@399: public: alpar@399: kpeter@610: /// Constructor. alpar@399: alpar@399: /// The constructor of the class. kpeter@610: /// kpeter@610: /// \param graph The digraph the algorithm runs on. kpeter@610: /// \param lower The lower bounds for the flow values on the arcs. kpeter@610: /// \param upper The upper bounds (capacities) for the flow values kpeter@610: /// on the arcs. kpeter@610: /// \param supply The signed supply values of the nodes. kpeter@610: Circulation(const Digraph &graph, const LowerMap &lower, kpeter@610: const UpperMap &upper, const SupplyMap &supply) kpeter@610: : _g(graph), _lo(&lower), _up(&upper), _supply(&supply), kpeter@610: _flow(NULL), _local_flow(false), _level(NULL), _local_level(false), kpeter@610: _excess(NULL) {} alpar@399: kpeter@402: /// Destructor. alpar@399: ~Circulation() { alpar@399: destroyStructures(); alpar@399: } alpar@399: kpeter@402: alpar@399: private: alpar@399: kpeter@622: bool checkBoundMaps() { kpeter@622: for (ArcIt e(_g);e!=INVALID;++e) { kpeter@622: if (_tol.less((*_up)[e], (*_lo)[e])) return false; kpeter@622: } kpeter@622: return true; kpeter@622: } kpeter@622: alpar@399: void createStructures() { alpar@399: _node_num = _el = countNodes(_g); alpar@399: alpar@399: if (!_flow) { alpar@399: _flow = Traits::createFlowMap(_g); alpar@399: _local_flow = true; alpar@399: } alpar@399: if (!_level) { alpar@399: _level = Traits::createElevator(_g, _node_num); alpar@399: _local_level = true; alpar@399: } alpar@399: if (!_excess) { alpar@399: _excess = new ExcessMap(_g); alpar@399: } alpar@399: } alpar@399: alpar@399: void destroyStructures() { alpar@399: if (_local_flow) { alpar@399: delete _flow; alpar@399: } alpar@399: if (_local_level) { alpar@399: delete _level; alpar@399: } alpar@399: if (_excess) { alpar@399: delete _excess; alpar@399: } alpar@399: } alpar@399: alpar@399: public: alpar@399: kpeter@610: /// Sets the lower bound map. alpar@399: kpeter@610: /// Sets the lower bound map. kpeter@402: /// \return (*this) kpeter@610: Circulation& lowerMap(const LowerMap& map) { alpar@399: _lo = ↦ alpar@399: return *this; alpar@399: } alpar@399: kpeter@610: /// Sets the upper bound (capacity) map. alpar@399: kpeter@610: /// Sets the upper bound (capacity) map. kpeter@402: /// \return (*this) kpeter@622: Circulation& upperMap(const UpperMap& map) { alpar@399: _up = ↦ alpar@399: return *this; alpar@399: } alpar@399: kpeter@610: /// Sets the supply map. alpar@399: kpeter@610: /// Sets the supply map. kpeter@402: /// \return (*this) kpeter@610: Circulation& supplyMap(const SupplyMap& map) { kpeter@610: _supply = ↦ alpar@399: return *this; alpar@399: } alpar@399: kpeter@402: /// \brief Sets the flow map. kpeter@402: /// alpar@399: /// Sets the flow map. kpeter@402: /// If you don't use this function before calling \ref run() or kpeter@402: /// \ref init(), an instance will be allocated automatically. kpeter@402: /// The destructor deallocates this automatically allocated map, kpeter@402: /// of course. kpeter@402: /// \return (*this) alpar@399: Circulation& flowMap(FlowMap& map) { alpar@399: if (_local_flow) { alpar@399: delete _flow; alpar@399: _local_flow = false; alpar@399: } alpar@399: _flow = ↦ alpar@399: return *this; alpar@399: } alpar@399: kpeter@402: /// \brief Sets the elevator used by algorithm. alpar@399: /// kpeter@402: /// Sets the elevator used by algorithm. kpeter@402: /// If you don't use this function before calling \ref run() or kpeter@402: /// \ref init(), an instance will be allocated automatically. kpeter@402: /// The destructor deallocates this automatically allocated elevator, kpeter@402: /// of course. kpeter@402: /// \return (*this) alpar@399: Circulation& elevator(Elevator& elevator) { alpar@399: if (_local_level) { alpar@399: delete _level; alpar@399: _local_level = false; alpar@399: } alpar@399: _level = &elevator; alpar@399: return *this; alpar@399: } alpar@399: kpeter@402: /// \brief Returns a const reference to the elevator. alpar@399: /// kpeter@402: /// Returns a const reference to the elevator. kpeter@402: /// kpeter@402: /// \pre Either \ref run() or \ref init() must be called before kpeter@402: /// using this function. kpeter@420: const Elevator& elevator() const { alpar@399: return *_level; alpar@399: } alpar@399: kpeter@689: /// \brief Sets the tolerance used by the algorithm. kpeter@402: /// kpeter@689: /// Sets the tolerance object used by the algorithm. kpeter@689: /// \return (*this) kpeter@688: Circulation& tolerance(const Tolerance& tolerance) { alpar@399: _tol = tolerance; alpar@399: return *this; alpar@399: } alpar@399: kpeter@402: /// \brief Returns a const reference to the tolerance. alpar@399: /// kpeter@689: /// Returns a const reference to the tolerance object used by kpeter@689: /// the algorithm. alpar@399: const Tolerance& tolerance() const { kpeter@688: return _tol; alpar@399: } alpar@399: kpeter@402: /// \name Execution Control kpeter@402: /// The simplest way to execute the algorithm is to call \ref run().\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@402: /// the \ref start() function. alpar@399: alpar@399: ///@{ alpar@399: alpar@399: /// Initializes the internal data structures. alpar@399: kpeter@402: /// Initializes the internal data structures and sets all flow values kpeter@402: /// to the lower bound. alpar@399: void init() alpar@399: { kpeter@622: LEMON_DEBUG(checkBoundMaps(), kpeter@622: "Upper bounds must be greater or equal to the lower bounds"); kpeter@622: alpar@399: createStructures(); alpar@399: alpar@399: for(NodeIt n(_g);n!=INVALID;++n) { alpar@611: (*_excess)[n] = (*_supply)[n]; alpar@399: } alpar@399: alpar@399: for (ArcIt e(_g);e!=INVALID;++e) { alpar@399: _flow->set(e, (*_lo)[e]); kpeter@581: (*_excess)[_g.target(e)] += (*_flow)[e]; kpeter@581: (*_excess)[_g.source(e)] -= (*_flow)[e]; alpar@399: } alpar@399: alpar@399: // global relabeling tested, but in general case it provides alpar@399: // worse performance for random digraphs alpar@399: _level->initStart(); alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: _level->initAddItem(n); alpar@399: _level->initFinish(); alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: if(_tol.positive((*_excess)[n])) alpar@399: _level->activate(n); alpar@399: } alpar@399: kpeter@402: /// Initializes the internal data structures using a greedy approach. alpar@399: kpeter@402: /// Initializes the internal data structures using a greedy approach kpeter@402: /// to construct the initial solution. alpar@399: void greedyInit() alpar@399: { kpeter@622: LEMON_DEBUG(checkBoundMaps(), kpeter@622: "Upper bounds must be greater or equal to the lower bounds"); kpeter@622: alpar@399: createStructures(); alpar@399: alpar@399: for(NodeIt n(_g);n!=INVALID;++n) { alpar@611: (*_excess)[n] = (*_supply)[n]; alpar@399: } alpar@399: alpar@399: for (ArcIt e(_g);e!=INVALID;++e) { kpeter@622: if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) { alpar@399: _flow->set(e, (*_up)[e]); kpeter@581: (*_excess)[_g.target(e)] += (*_up)[e]; kpeter@581: (*_excess)[_g.source(e)] -= (*_up)[e]; kpeter@622: } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) { alpar@399: _flow->set(e, (*_lo)[e]); kpeter@581: (*_excess)[_g.target(e)] += (*_lo)[e]; kpeter@581: (*_excess)[_g.source(e)] -= (*_lo)[e]; alpar@399: } else { kpeter@641: Value fc = -(*_excess)[_g.target(e)]; alpar@399: _flow->set(e, fc); kpeter@581: (*_excess)[_g.target(e)] = 0; kpeter@581: (*_excess)[_g.source(e)] -= fc; alpar@399: } alpar@399: } alpar@399: alpar@399: _level->initStart(); alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: _level->initAddItem(n); alpar@399: _level->initFinish(); alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: if(_tol.positive((*_excess)[n])) alpar@399: _level->activate(n); alpar@399: } alpar@399: kpeter@402: ///Executes the algorithm alpar@399: kpeter@402: ///This function executes the algorithm. kpeter@402: /// kpeter@402: ///\return \c true if a feasible circulation is found. alpar@399: /// alpar@399: ///\sa barrier() kpeter@402: ///\sa barrierMap() alpar@399: bool start() alpar@399: { alpar@399: alpar@399: Node act; alpar@399: Node bact=INVALID; alpar@399: Node last_activated=INVALID; alpar@399: while((act=_level->highestActive())!=INVALID) { alpar@399: int actlevel=(*_level)[act]; alpar@399: int mlevel=_node_num; kpeter@641: Value exc=(*_excess)[act]; alpar@399: alpar@399: for(OutArcIt e(_g,act);e!=INVALID; ++e) { alpar@399: Node v = _g.target(e); kpeter@641: Value fc=(*_up)[e]-(*_flow)[e]; alpar@399: if(!_tol.positive(fc)) continue; alpar@399: if((*_level)[v]set(e, (*_flow)[e] + exc); kpeter@581: (*_excess)[v] += exc; alpar@399: if(!_level->active(v) && _tol.positive((*_excess)[v])) alpar@399: _level->activate(v); kpeter@581: (*_excess)[act] = 0; alpar@399: _level->deactivate(act); alpar@399: goto next_l; alpar@399: } alpar@399: else { alpar@399: _flow->set(e, (*_up)[e]); kpeter@581: (*_excess)[v] += fc; alpar@399: if(!_level->active(v) && _tol.positive((*_excess)[v])) alpar@399: _level->activate(v); alpar@399: exc-=fc; alpar@399: } alpar@399: } alpar@399: else if((*_level)[v]set(e, (*_flow)[e] - exc); kpeter@581: (*_excess)[v] += exc; alpar@399: if(!_level->active(v) && _tol.positive((*_excess)[v])) alpar@399: _level->activate(v); kpeter@581: (*_excess)[act] = 0; alpar@399: _level->deactivate(act); alpar@399: goto next_l; alpar@399: } alpar@399: else { alpar@399: _flow->set(e, (*_lo)[e]); kpeter@581: (*_excess)[v] += fc; alpar@399: if(!_level->active(v) && _tol.positive((*_excess)[v])) alpar@399: _level->activate(v); alpar@399: exc-=fc; alpar@399: } alpar@399: } alpar@399: else if((*_level)[v]deactivate(act); alpar@399: else if(mlevel==_node_num) { alpar@399: _level->liftHighestActiveToTop(); alpar@399: _el = _node_num; alpar@399: return false; alpar@399: } alpar@399: else { alpar@399: _level->liftHighestActive(mlevel+1); alpar@399: if(_level->onLevel(actlevel)==0) { alpar@399: _el = actlevel; alpar@399: return false; alpar@399: } alpar@399: } alpar@399: next_l: alpar@399: ; alpar@399: } alpar@399: return true; alpar@399: } alpar@399: kpeter@402: /// Runs the algorithm. alpar@399: kpeter@402: /// This function runs the algorithm. kpeter@402: /// kpeter@402: /// \return \c true if a feasible circulation is found. kpeter@402: /// kpeter@402: /// \note Apart from the return value, c.run() is just a shortcut of kpeter@402: /// the following code. alpar@399: /// \code kpeter@402: /// c.greedyInit(); kpeter@402: /// c.start(); alpar@399: /// \endcode alpar@399: bool run() { alpar@399: greedyInit(); alpar@399: return start(); alpar@399: } alpar@399: alpar@399: /// @} alpar@399: alpar@399: /// \name Query Functions kpeter@402: /// The results of the circulation algorithm can be obtained using kpeter@402: /// these functions.\n kpeter@402: /// Either \ref run() or \ref start() should be called before kpeter@402: /// using them. alpar@399: alpar@399: ///@{ alpar@399: kpeter@641: /// \brief Returns the flow value on the given arc. kpeter@402: /// kpeter@641: /// Returns the flow value on the given arc. kpeter@402: /// kpeter@402: /// \pre Either \ref run() or \ref init() must be called before kpeter@402: /// using this function. kpeter@641: Value flow(const Arc& arc) const { kpeter@402: return (*_flow)[arc]; kpeter@402: } kpeter@402: kpeter@402: /// \brief Returns a const reference to the flow map. kpeter@402: /// kpeter@402: /// Returns a const reference to the arc map storing the found flow. kpeter@402: /// kpeter@402: /// \pre Either \ref run() or \ref init() must be called before kpeter@402: /// using this function. kpeter@420: const FlowMap& flowMap() const { kpeter@402: return *_flow; kpeter@402: } kpeter@402: alpar@399: /** kpeter@402: \brief Returns \c true if the given node is in a barrier. kpeter@402: alpar@399: Barrier is a set \e B of nodes for which kpeter@402: kpeter@610: \f[ \sum_{uv\in A: u\in B} upper(uv) - kpeter@610: \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f] kpeter@402: kpeter@402: holds. The existence of a set with this property prooves that a kpeter@402: feasible circualtion cannot exist. kpeter@402: kpeter@402: This function returns \c true if the given node is in the found kpeter@402: barrier. If a feasible circulation is found, the function kpeter@402: gives back \c false for every node. kpeter@402: kpeter@402: \pre Either \ref run() or \ref init() must be called before kpeter@402: using this function. kpeter@402: kpeter@402: \sa barrierMap() alpar@399: \sa checkBarrier() alpar@399: */ kpeter@420: bool barrier(const Node& node) const kpeter@402: { kpeter@402: return (*_level)[node] >= _el; kpeter@402: } kpeter@402: kpeter@402: /// \brief Gives back a barrier. kpeter@402: /// kpeter@402: /// This function sets \c bar to the characteristic vector of the kpeter@402: /// found barrier. \c bar should be a \ref concepts::WriteMap "writable" kpeter@402: /// node map with \c bool (or convertible) value type. kpeter@402: /// kpeter@402: /// If a feasible circulation is found, the function gives back an kpeter@402: /// empty set, so \c bar[v] will be \c false for all nodes \c v. kpeter@402: /// kpeter@402: /// \note This function calls \ref barrier() for each node, kpeter@559: /// so it runs in O(n) time. kpeter@402: /// kpeter@402: /// \pre Either \ref run() or \ref init() must be called before kpeter@402: /// using this function. kpeter@402: /// kpeter@402: /// \sa barrier() kpeter@402: /// \sa checkBarrier() kpeter@402: template kpeter@420: void barrierMap(BarrierMap &bar) const alpar@399: { alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: bar.set(n, (*_level)[n] >= _el); alpar@399: } alpar@399: alpar@399: /// @} alpar@399: alpar@399: /// \name Checker Functions kpeter@402: /// The feasibility of the results can be checked using kpeter@402: /// these functions.\n kpeter@402: /// Either \ref run() or \ref start() should be called before kpeter@402: /// using them. alpar@399: alpar@399: ///@{ alpar@399: kpeter@402: ///Check if the found flow is a feasible circulation kpeter@402: kpeter@402: ///Check if the found flow is a feasible circulation, kpeter@402: /// kpeter@420: bool checkFlow() const { alpar@399: for(ArcIt e(_g);e!=INVALID;++e) alpar@399: if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false; alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: { kpeter@641: Value dif=-(*_supply)[n]; alpar@399: for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e]; alpar@399: for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e]; alpar@399: if(_tol.negative(dif)) return false; alpar@399: } alpar@399: return true; alpar@399: } alpar@399: alpar@399: ///Check whether or not the last execution provides a barrier alpar@399: kpeter@402: ///Check whether or not the last execution provides a barrier. alpar@399: ///\sa barrier() kpeter@402: ///\sa barrierMap() kpeter@420: bool checkBarrier() const alpar@399: { kpeter@641: Value delta=0; kpeter@641: Value inf_cap = std::numeric_limits::has_infinity ? kpeter@641: std::numeric_limits::infinity() : kpeter@641: std::numeric_limits::max(); alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: if(barrier(n)) kpeter@610: delta-=(*_supply)[n]; alpar@399: for(ArcIt e(_g);e!=INVALID;++e) alpar@399: { alpar@399: Node s=_g.source(e); alpar@399: Node t=_g.target(e); kpeter@622: if(barrier(s)&&!barrier(t)) { kpeter@622: if (_tol.less(inf_cap - (*_up)[e], delta)) return false; kpeter@622: delta+=(*_up)[e]; kpeter@622: } alpar@399: else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e]; alpar@399: } alpar@399: return _tol.negative(delta); alpar@399: } alpar@399: alpar@399: /// @} alpar@399: alpar@399: }; alpar@399: alpar@399: } alpar@399: alpar@399: #endif