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@399: * Copyright (C) 2003-2008 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 alpar@399: #include alpar@399: #include alpar@399: alpar@399: ///\ingroup max_flow alpar@399: ///\file alpar@399: ///\brief Push-prelabel 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. alpar@399: /// \param _Graph Digraph type. alpar@399: /// \param _CapacityMap Type of capacity map. alpar@399: template alpar@399: struct CirculationDefaultTraits { alpar@399: alpar@399: /// \brief The digraph type the algorithm runs on. alpar@399: typedef _Graph Digraph; alpar@399: alpar@399: /// \brief The type of the map that stores the circulation lower alpar@399: /// bound. alpar@399: /// alpar@399: /// The type of the map that stores the circulation lower bound. alpar@399: /// It must meet the \ref concepts::ReadMap "ReadMap" concept. alpar@399: typedef _LCapMap LCapMap; alpar@399: alpar@399: /// \brief The type of the map that stores the circulation upper alpar@399: /// bound. alpar@399: /// alpar@399: /// The type of the map that stores the circulation upper bound. alpar@399: /// It must meet the \ref concepts::ReadMap "ReadMap" concept. alpar@399: typedef _UCapMap UCapMap; alpar@399: alpar@399: /// \brief The type of the map that stores the upper bound of alpar@399: /// node excess. alpar@399: /// alpar@399: /// The type of the map that stores the lower bound of node alpar@399: /// excess. It must meet the \ref concepts::ReadMap "ReadMap" alpar@399: /// concept. alpar@399: typedef _DeltaMap DeltaMap; alpar@399: alpar@399: /// \brief The type of the length of the arcs. alpar@399: typedef typename DeltaMap::Value Value; alpar@399: alpar@399: /// \brief The map type that stores the flow values. alpar@399: /// alpar@399: /// The map type that stores the flow values. alpar@399: /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. alpar@399: typedef typename Digraph::template ArcMap FlowMap; alpar@399: alpar@399: /// \brief Instantiates a FlowMap. alpar@399: /// alpar@399: /// This function instantiates a \ref FlowMap. alpar@399: /// \param digraph The digraph, to 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: alpar@399: /// \brief The eleavator type used by Circulation algorithm. alpar@399: /// alpar@399: /// The elevator type used by Circulation algorithm. alpar@399: /// alpar@399: /// \sa Elevator alpar@399: /// \sa LinkedElevator alpar@399: typedef lemon::Elevator Elevator; alpar@399: alpar@399: /// \brief Instantiates an Elevator. alpar@399: /// alpar@399: /// This function instantiates a \ref Elevator. alpar@399: /// \param digraph The digraph, to 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. alpar@399: typedef lemon::Tolerance Tolerance; alpar@399: alpar@399: }; alpar@399: alpar@399: ///Push-relabel algorithm for the Network Circulation Problem. alpar@399: alpar@399: /** alpar@399: \ingroup max_flow alpar@399: This class implements a push-relabel algorithm alpar@399: or the Network Circulation Problem. alpar@399: The exact formulation of this problem is the following. alpar@399: \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq alpar@399: -delta(v)\quad \forall v\in V \f] alpar@399: \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f] alpar@399: */ alpar@399: template, alpar@399: class _UCapMap=_LCapMap, alpar@399: class _DeltaMap=typename _Graph::template NodeMap< alpar@399: typename _UCapMap::Value>, alpar@399: class _Traits=CirculationDefaultTraits<_Graph, _LCapMap, alpar@399: _UCapMap, _DeltaMap> > alpar@399: class Circulation { alpar@399: alpar@399: typedef _Traits Traits; alpar@399: typedef typename Traits::Digraph Digraph; alpar@399: TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); alpar@399: alpar@399: typedef typename Traits::Value Value; alpar@399: alpar@399: typedef typename Traits::LCapMap LCapMap; alpar@399: typedef typename Traits::UCapMap UCapMap; alpar@399: typedef typename Traits::DeltaMap DeltaMap; alpar@399: typedef typename Traits::FlowMap FlowMap; alpar@399: typedef typename Traits::Elevator Elevator; alpar@399: typedef typename Traits::Tolerance Tolerance; alpar@399: alpar@399: typedef typename Digraph::template NodeMap ExcessMap; alpar@399: alpar@399: const Digraph &_g; alpar@399: int _node_num; alpar@399: alpar@399: const LCapMap *_lo; alpar@399: const UCapMap *_up; alpar@399: const DeltaMap *_delta; alpar@399: alpar@399: FlowMap *_flow; alpar@399: bool _local_flow; alpar@399: alpar@399: Elevator* _level; alpar@399: bool _local_level; alpar@399: 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: alpar@399: ///\name Named template parameters alpar@399: alpar@399: ///@{ alpar@399: alpar@399: template alpar@399: struct DefFlowMapTraits : public Traits { alpar@399: typedef _FlowMap 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 alpar@399: /// type alpar@399: template alpar@399: struct DefFlowMap alpar@399: : public Circulation > { alpar@399: typedef Circulation > Create; alpar@399: }; alpar@399: alpar@399: template alpar@399: struct DefElevatorTraits : public Traits { alpar@399: typedef _Elevator 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 alpar@399: /// type alpar@399: template alpar@399: struct DefElevator alpar@399: : public Circulation > { alpar@399: typedef Circulation > Create; alpar@399: }; alpar@399: alpar@399: template alpar@399: struct DefStandardElevatorTraits : public Traits { alpar@399: typedef _Elevator 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 alpar@399: /// Elevator type alpar@399: /// alpar@399: /// \ref named-templ-param "Named parameter" for setting Elevator alpar@399: /// type. The Elevator should be standard constructor interface, ie. alpar@399: /// the digraph and the maximum level should be passed to it. alpar@399: template alpar@399: struct DefStandardElevator alpar@399: : public Circulation > { alpar@399: 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: alpar@399: /// The constructor of the class. alpar@399: alpar@399: /// The constructor of the class. alpar@399: /// \param g The digraph the algorithm runs on. alpar@399: /// \param lo The lower bound capacity of the arcs. alpar@399: /// \param up The upper bound capacity of the arcs. alpar@399: /// \param delta The lower bound on node excess. alpar@399: Circulation(const Digraph &g,const LCapMap &lo, alpar@399: const UCapMap &up,const DeltaMap &delta) alpar@399: : _g(g), _node_num(), alpar@399: _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false), alpar@399: _level(0), _local_level(false), _excess(0), _el() {} alpar@399: alpar@399: /// Destrcutor. alpar@399: ~Circulation() { alpar@399: destroyStructures(); alpar@399: } alpar@399: alpar@399: private: alpar@399: 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: alpar@399: /// Sets the lower bound capacity map. alpar@399: alpar@399: /// Sets the lower bound capacity map. alpar@399: /// \return \c (*this) alpar@399: Circulation& lowerCapMap(const LCapMap& map) { alpar@399: _lo = ↦ alpar@399: return *this; alpar@399: } alpar@399: alpar@399: /// Sets the upper bound capacity map. alpar@399: alpar@399: /// Sets the upper bound capacity map. alpar@399: /// \return \c (*this) alpar@399: Circulation& upperCapMap(const LCapMap& map) { alpar@399: _up = ↦ alpar@399: return *this; alpar@399: } alpar@399: alpar@399: /// Sets the lower bound map on excess. alpar@399: alpar@399: /// Sets the lower bound map on excess. alpar@399: /// \return \c (*this) alpar@399: Circulation& deltaMap(const DeltaMap& map) { alpar@399: _delta = ↦ alpar@399: return *this; alpar@399: } alpar@399: alpar@399: /// Sets the flow map. alpar@399: alpar@399: /// Sets the flow map. alpar@399: /// \return \c (*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: alpar@399: /// Returns the flow map. alpar@399: alpar@399: /// \return The flow map. alpar@399: /// alpar@399: const FlowMap& flowMap() { alpar@399: return *_flow; alpar@399: } alpar@399: alpar@399: /// Sets the elevator. alpar@399: alpar@399: /// Sets the elevator. alpar@399: /// \return \c (*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: alpar@399: /// Returns the elevator. alpar@399: alpar@399: /// \return The elevator. alpar@399: /// alpar@399: const Elevator& elevator() { alpar@399: return *_level; alpar@399: } alpar@399: alpar@399: /// Sets the tolerance used by algorithm. alpar@399: alpar@399: /// Sets the tolerance used by algorithm. alpar@399: /// alpar@399: Circulation& tolerance(const Tolerance& tolerance) const { alpar@399: _tol = tolerance; alpar@399: return *this; alpar@399: } alpar@399: alpar@399: /// Returns the tolerance used by algorithm. alpar@399: alpar@399: /// Returns the tolerance used by algorithm. alpar@399: /// alpar@399: const Tolerance& tolerance() const { alpar@399: return tolerance; alpar@399: } alpar@399: alpar@399: /// \name Execution control alpar@399: /// The simplest way to execute the algorithm is to use one of the alpar@399: /// member functions called \c run(). alpar@399: /// \n alpar@399: /// If you need more control on initial solution or execution then alpar@399: /// you have to call one \ref init() function and then the start() alpar@399: /// function. alpar@399: alpar@399: ///@{ alpar@399: alpar@399: /// Initializes the internal data structures. alpar@399: alpar@399: /// Initializes the internal data structures. This function sets alpar@399: /// all flow values to the lower bound. alpar@399: /// \return This function returns false if the initialization alpar@399: /// process found a barrier. alpar@399: void init() alpar@399: { alpar@399: createStructures(); alpar@399: alpar@399: for(NodeIt n(_g);n!=INVALID;++n) { alpar@399: _excess->set(n, (*_delta)[n]); alpar@399: } alpar@399: alpar@399: for (ArcIt e(_g);e!=INVALID;++e) { alpar@399: _flow->set(e, (*_lo)[e]); alpar@399: _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]); alpar@399: _excess->set(_g.source(e), (*_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: alpar@399: /// Initializes the internal data structures. alpar@399: alpar@399: /// Initializes the internal data structures. This functions uses alpar@399: /// greedy approach to construct the initial solution. alpar@399: void greedyInit() alpar@399: { alpar@399: createStructures(); alpar@399: alpar@399: for(NodeIt n(_g);n!=INVALID;++n) { alpar@399: _excess->set(n, (*_delta)[n]); alpar@399: } alpar@399: alpar@399: for (ArcIt e(_g);e!=INVALID;++e) { alpar@399: if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) { alpar@399: _flow->set(e, (*_up)[e]); alpar@399: _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]); alpar@399: _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]); alpar@399: } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) { alpar@399: _flow->set(e, (*_lo)[e]); alpar@399: _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]); alpar@399: _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]); alpar@399: } else { alpar@399: Value fc = -(*_excess)[_g.target(e)]; alpar@399: _flow->set(e, fc); alpar@399: _excess->set(_g.target(e), 0); alpar@399: _excess->set(_g.source(e), (*_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: alpar@399: ///Starts the algorithm alpar@399: alpar@399: ///This function starts the algorithm. alpar@399: ///\return This function returns true if it found a feasible circulation. alpar@399: /// alpar@399: ///\sa barrier() 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; alpar@399: Value exc=(*_excess)[act]; alpar@399: alpar@399: for(OutArcIt e(_g,act);e!=INVALID; ++e) { alpar@399: Node v = _g.target(e); alpar@399: Value fc=(*_up)[e]-(*_flow)[e]; alpar@399: if(!_tol.positive(fc)) continue; alpar@399: if((*_level)[v]set(e, (*_flow)[e] + exc); alpar@399: _excess->set(v, (*_excess)[v] + exc); alpar@399: if(!_level->active(v) && _tol.positive((*_excess)[v])) alpar@399: _level->activate(v); alpar@399: _excess->set(act,0); alpar@399: _level->deactivate(act); alpar@399: goto next_l; alpar@399: } alpar@399: else { alpar@399: _flow->set(e, (*_up)[e]); alpar@399: _excess->set(v, (*_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); alpar@399: _excess->set(v, (*_excess)[v] + exc); alpar@399: if(!_level->active(v) && _tol.positive((*_excess)[v])) alpar@399: _level->activate(v); alpar@399: _excess->set(act,0); alpar@399: _level->deactivate(act); alpar@399: goto next_l; alpar@399: } alpar@399: else { alpar@399: _flow->set(e, (*_lo)[e]); alpar@399: _excess->set(v, (*_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(act, exc); alpar@399: if(!_tol.positive(exc)) _level->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: alpar@399: /// Runs the circulation algorithm. alpar@399: alpar@399: /// Runs the circulation algorithm. alpar@399: /// \note fc.run() is just a shortcut of the following code. alpar@399: /// \code alpar@399: /// fc.greedyInit(); alpar@399: /// return fc.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 alpar@399: /// The result of the %Circulation algorithm can be obtained using alpar@399: /// these functions. alpar@399: /// \n alpar@399: /// Before the use of these functions, alpar@399: /// either run() or start() must be called. alpar@399: alpar@399: ///@{ alpar@399: alpar@399: /** alpar@399: \brief Returns a barrier alpar@399: alpar@399: Barrier is a set \e B of nodes for which alpar@399: \f[ \sum_{v\in B}-delta(v)< alpar@399: \sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f] alpar@399: holds. The existence of a set with this property prooves that a feasible alpar@399: flow cannot exists. alpar@399: \sa checkBarrier() alpar@399: \sa run() alpar@399: */ alpar@399: template alpar@399: void barrierMap(GT &bar) 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: ///Returns true if the node is in the barrier alpar@399: alpar@399: ///Returns true if the node is in the barrier alpar@399: ///\sa barrierMap() alpar@399: bool barrier(const Node& node) alpar@399: { alpar@399: return (*_level)[node] >= _el; alpar@399: } alpar@399: alpar@399: /// \brief Returns the flow on the arc. alpar@399: /// alpar@399: /// Sets the \c flowMap to the flow on the arcs. This method can alpar@399: /// be called after the second phase of algorithm. alpar@399: Value flow(const Arc& arc) const { alpar@399: return (*_flow)[arc]; alpar@399: } alpar@399: alpar@399: /// @} alpar@399: alpar@399: /// \name Checker Functions alpar@399: /// The feasibility of the results can be checked using alpar@399: /// these functions. alpar@399: /// \n alpar@399: /// Before the use of these functions, alpar@399: /// either run() or start() must be called. alpar@399: alpar@399: ///@{ alpar@399: alpar@399: ///Check if the \c flow is a feasible circulation alpar@399: bool checkFlow() { 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: { alpar@399: Value dif=-(*_delta)[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: alpar@399: ///Check whether or not the last execution provides a barrier alpar@399: ///\sa barrier() alpar@399: bool checkBarrier() alpar@399: { alpar@399: Value delta=0; alpar@399: for(NodeIt n(_g);n!=INVALID;++n) alpar@399: if(barrier(n)) alpar@399: delta-=(*_delta)[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); alpar@399: if(barrier(s)&&!barrier(t)) delta+=(*_up)[e]; 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