alpar@2375: /* -*- C++ -*- alpar@2375: * alpar@2391: * This file is a part of LEMON, a generic C++ optimization library alpar@2391: * alpar@2391: * Copyright (C) 2003-2007 alpar@2391: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport alpar@2375: * (Egervary Research Group on Combinatorial Optimization, EGRES). alpar@2375: * alpar@2375: * Permission to use, modify and distribute this software is granted alpar@2375: * provided that this copyright notice appears in all copies. For alpar@2375: * precise terms see the accompanying LICENSE file. alpar@2375: * alpar@2375: * This software is provided "AS IS" with no warranty of any kind, alpar@2375: * express or implied, and with no claim as to its suitability for any alpar@2375: * purpose. alpar@2375: * alpar@2375: */ alpar@2375: alpar@2375: #ifndef LEMON_CIRCULATION_H alpar@2375: #define LEMON_CIRCULATION_H alpar@2375: alpar@2375: #include alpar@2375: #include alpar@2375: #include alpar@2375: #include alpar@2375: #include alpar@2375: deba@2376: ///\ingroup max_flow alpar@2375: ///\file alpar@2375: ///\brief Push-prelabel algorithm for finding a feasible circulation. alpar@2375: /// alpar@2375: namespace lemon { alpar@2375: deba@2526: /// \brief Default traits class of Circulation class. deba@2526: /// deba@2526: /// Default traits class of Circulation class. deba@2526: /// \param _Graph Graph type. deba@2526: /// \param _CapacityMap Type of capacity map. deba@2526: template deba@2526: struct CirculationDefaultTraits { deba@2526: deba@2526: /// \brief The graph type the algorithm runs on. deba@2526: typedef _Graph Graph; deba@2526: deba@2526: /// \brief The type of the map that stores the circulation lower deba@2526: /// bound. deba@2526: /// deba@2526: /// The type of the map that stores the circulation lower bound. deba@2526: /// It must meet the \ref concepts::ReadMap "ReadMap" concept. deba@2526: typedef _LCapMap LCapMap; deba@2526: deba@2526: /// \brief The type of the map that stores the circulation upper deba@2526: /// bound. deba@2526: /// deba@2526: /// The type of the map that stores the circulation upper bound. deba@2526: /// It must meet the \ref concepts::ReadMap "ReadMap" concept. deba@2526: typedef _UCapMap UCapMap; deba@2526: deba@2526: /// \brief The type of the map that stores the upper bound of deba@2526: /// node excess. deba@2526: /// deba@2526: /// The type of the map that stores the lower bound of node deba@2526: /// excess. It must meet the \ref concepts::ReadMap "ReadMap" deba@2526: /// concept. deba@2526: typedef _DeltaMap DeltaMap; deba@2526: deba@2526: /// \brief The type of the length of the edges. deba@2526: typedef typename DeltaMap::Value Value; deba@2526: deba@2526: /// \brief The map type that stores the flow values. deba@2526: /// deba@2526: /// The map type that stores the flow values. deba@2526: /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. deba@2526: typedef typename Graph::template EdgeMap FlowMap; deba@2526: deba@2526: /// \brief Instantiates a FlowMap. deba@2526: /// deba@2526: /// This function instantiates a \ref FlowMap. deba@2526: /// \param graph The graph, to which we would like to define the flow map. deba@2526: static FlowMap* createFlowMap(const Graph& graph) { deba@2526: return new FlowMap(graph); deba@2526: } deba@2526: deba@2526: /// \brief The eleavator type used by Circulation algorithm. deba@2526: /// deba@2526: /// The elevator type used by Circulation algorithm. deba@2526: /// deba@2526: /// \sa Elevator deba@2526: /// \sa LinkedElevator deba@2526: typedef Elevator Elevator; deba@2526: deba@2526: /// \brief Instantiates an Elevator. deba@2526: /// deba@2526: /// This function instantiates a \ref Elevator. deba@2526: /// \param graph The graph, to which we would like to define the elevator. deba@2526: /// \param max_level The maximum level of the elevator. deba@2526: static Elevator* createElevator(const Graph& graph, int max_level) { deba@2526: return new Elevator(graph, max_level); deba@2526: } deba@2526: deba@2526: /// \brief The tolerance used by the algorithm deba@2526: /// deba@2526: /// The tolerance used by the algorithm to handle inexact computation. deba@2526: typedef Tolerance Tolerance; deba@2526: deba@2526: }; deba@2526: deba@2526: ///Push-relabel algorithm for the Network Circulation Problem. alpar@2375: deba@2378: ///\ingroup max_flow deba@2526: ///This class implements a push-relabel algorithm alpar@2375: ///for the Network Circulation Problem. alpar@2375: ///The exact formulation of this problem is the following. alpar@2450: /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq -delta(v)\quad \forall v\in V \f] alpar@2375: /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f] alpar@2375: /// deba@2526: template, deba@2526: class _UCapMap=_LCapMap, deba@2526: class _DeltaMap=typename _Graph::template NodeMap< deba@2526: typename _UCapMap::Value>, deba@2526: class _Traits=CirculationDefaultTraits<_Graph, _LCapMap, deba@2526: _UCapMap, _DeltaMap> > alpar@2375: class Circulation { deba@2526: deba@2526: typedef _Traits Traits; deba@2526: typedef typename Traits::Graph Graph; deba@2526: GRAPH_TYPEDEFS(typename Graph); deba@2526: deba@2526: typedef typename Traits::Value Value; deba@2526: deba@2526: typedef typename Traits::LCapMap LCapMap; deba@2526: typedef typename Traits::UCapMap UCapMap; deba@2526: typedef typename Traits::DeltaMap DeltaMap; deba@2526: typedef typename Traits::FlowMap FlowMap; deba@2526: typedef typename Traits::Elevator Elevator; deba@2526: typedef typename Traits::Tolerance Tolerance; deba@2526: deba@2526: typedef typename Graph::template NodeMap ExcessMap; alpar@2375: alpar@2375: const Graph &_g; alpar@2375: int _node_num; alpar@2375: deba@2526: const LCapMap *_lo; deba@2526: const UCapMap *_up; deba@2526: const DeltaMap *_delta; deba@2526: deba@2526: FlowMap *_flow; deba@2526: bool _local_flow; deba@2526: deba@2526: Elevator* _level; deba@2526: bool _local_level; deba@2526: deba@2526: ExcessMap* _excess; deba@2526: deba@2526: Tolerance _tol; deba@2526: int _el; alpar@2375: alpar@2375: public: deba@2526: deba@2526: typedef Circulation Create; deba@2526: deba@2526: ///\name Named template parameters deba@2526: deba@2526: ///@{ deba@2526: deba@2526: template deba@2526: struct DefFlowMapTraits : public Traits { deba@2526: typedef _FlowMap FlowMap; deba@2526: static FlowMap *createFlowMap(const Graph&) { deba@2526: throw UninitializedParameter(); deba@2526: } deba@2526: }; deba@2526: deba@2526: /// \brief \ref named-templ-param "Named parameter" for setting deba@2526: /// FlowMap type deba@2526: /// deba@2526: /// \ref named-templ-param "Named parameter" for setting FlowMap deba@2526: /// type deba@2526: template deba@2526: struct DefFlowMap deba@2526: : public Circulation > { deba@2526: typedef Circulation > Create; deba@2526: }; deba@2526: deba@2526: template deba@2526: struct DefElevatorTraits : public Traits { deba@2526: typedef _Elevator Elevator; deba@2526: static Elevator *createElevator(const Graph&, int) { deba@2526: throw UninitializedParameter(); deba@2526: } deba@2526: }; deba@2526: deba@2526: /// \brief \ref named-templ-param "Named parameter" for setting deba@2526: /// Elevator type deba@2526: /// deba@2526: /// \ref named-templ-param "Named parameter" for setting Elevator deba@2526: /// type deba@2526: template deba@2526: struct DefElevator deba@2526: : public Circulation > { deba@2526: typedef Circulation > Create; deba@2526: }; deba@2526: deba@2526: template deba@2526: struct DefStandardElevatorTraits : public Traits { deba@2526: typedef _Elevator Elevator; deba@2526: static Elevator *createElevator(const Graph& graph, int max_level) { deba@2526: return new Elevator(graph, max_level); deba@2526: } deba@2526: }; deba@2526: deba@2526: /// \brief \ref named-templ-param "Named parameter" for setting deba@2526: /// Elevator type deba@2526: /// deba@2526: /// \ref named-templ-param "Named parameter" for setting Elevator deba@2526: /// type. The Elevator should be standard constructor interface, ie. deba@2526: /// the graph and the maximum level should be passed to it. deba@2526: template deba@2526: struct DefStandardElevator deba@2526: : public Circulation > { deba@2526: typedef Circulation > Create; deba@2526: }; deba@2526: deba@2526: /// @} deba@2526: deba@2527: protected: deba@2527: deba@2527: Circulation() {} deba@2527: deba@2527: public: deba@2527: deba@2526: /// The constructor of the class. deba@2526: deba@2526: /// The constructor of the class. deba@2526: /// \param g The directed graph the algorithm runs on. deba@2526: /// \param lo The lower bound capacity of the edges. deba@2526: /// \param up The upper bound capacity of the edges. deba@2526: /// \param delta The lower bound on node excess. deba@2526: Circulation(const Graph &g,const LCapMap &lo, deba@2526: const UCapMap &up,const DeltaMap &delta) deba@2526: : _g(g), _node_num(), deba@2526: _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false), deba@2526: _level(0), _local_level(false), _excess(0), _el() {} deba@2526: deba@2526: /// Destrcutor. deba@2526: ~Circulation() { deba@2526: destroyStructures(); alpar@2375: } deba@2526: alpar@2375: private: alpar@2375: deba@2526: void createStructures() { deba@2526: _node_num = _el = countNodes(_g); deba@2526: deba@2526: if (!_flow) { deba@2526: _flow = Traits::createFlowMap(_g); deba@2526: _local_flow = true; deba@2526: } deba@2526: if (!_level) { deba@2526: _level = Traits::createElevator(_g, _node_num); deba@2526: _local_level = true; deba@2526: } deba@2526: if (!_excess) { deba@2526: _excess = new ExcessMap(_g); deba@2526: } alpar@2375: } alpar@2375: deba@2526: void destroyStructures() { deba@2526: if (_local_flow) { deba@2526: delete _flow; deba@2526: } deba@2526: if (_local_level) { deba@2526: delete _level; deba@2526: } deba@2526: if (_excess) { deba@2526: delete _excess; deba@2526: } alpar@2375: } alpar@2375: alpar@2375: public: deba@2526: deba@2526: /// Sets the lower bound capacity map. deba@2526: deba@2526: /// Sets the lower bound capacity map. deba@2526: /// \return \c (*this) deba@2526: Circulation& lowerCapMap(const LCapMap& map) { deba@2526: _lo = ↦ deba@2526: return *this; deba@2526: } deba@2526: deba@2526: /// Sets the upper bound capacity map. deba@2526: deba@2526: /// Sets the upper bound capacity map. deba@2526: /// \return \c (*this) deba@2526: Circulation& upperCapMap(const LCapMap& map) { deba@2526: _up = ↦ deba@2526: return *this; deba@2526: } deba@2526: deba@2526: /// Sets the lower bound map on excess. deba@2526: deba@2526: /// Sets the lower bound map on excess. deba@2526: /// \return \c (*this) deba@2526: Circulation& deltaMap(const DeltaMap& map) { deba@2526: _delta = ↦ deba@2526: return *this; deba@2526: } deba@2526: deba@2526: /// Sets the flow map. deba@2526: deba@2526: /// Sets the flow map. deba@2526: /// \return \c (*this) deba@2526: Circulation& flowMap(FlowMap& map) { deba@2526: if (_local_flow) { deba@2526: delete _flow; deba@2526: _local_flow = false; deba@2526: } deba@2526: _flow = ↦ deba@2526: return *this; deba@2526: } deba@2526: deba@2526: /// Returns the flow map. deba@2526: deba@2526: /// \return The flow map. deba@2526: /// deba@2526: const FlowMap& flowMap() { deba@2526: return *_flow; deba@2526: } deba@2526: deba@2526: /// Sets the elevator. deba@2526: deba@2526: /// Sets the elevator. deba@2526: /// \return \c (*this) deba@2526: Circulation& elevator(Elevator& elevator) { deba@2526: if (_local_level) { deba@2526: delete _level; deba@2526: _local_level = false; deba@2526: } deba@2526: _level = &elevator; deba@2526: return *this; deba@2526: } deba@2526: deba@2526: /// Returns the elevator. deba@2526: deba@2526: /// \return The elevator. deba@2526: /// deba@2526: const Elevator& elevator() { deba@2526: return *_level; deba@2526: } deba@2526: deba@2526: /// Sets the tolerance used by algorithm. deba@2526: deba@2526: /// Sets the tolerance used by algorithm. deba@2526: /// deba@2526: Circulation& tolerance(const Tolerance& tolerance) const { deba@2526: _tol = tolerance; deba@2526: return *this; deba@2526: } deba@2526: deba@2526: /// Returns the tolerance used by algorithm. deba@2526: deba@2526: /// Returns the tolerance used by algorithm. deba@2526: /// deba@2526: const Tolerance& tolerance() const { deba@2526: return tolerance; deba@2526: } deba@2526: kpeter@2533: /// \name Execution control kpeter@2533: /// The simplest way to execute the algorithm is to use one of the kpeter@2533: /// member functions called \c run(). deba@2526: /// \n deba@2526: /// If you need more control on initial solution or execution then deba@2526: /// you have to call one \ref init() function and then the start() deba@2526: /// function. deba@2526: deba@2526: ///@{ deba@2526: deba@2526: /// Initializes the internal data structures. deba@2526: deba@2526: /// Initializes the internal data structures. This function sets deba@2526: /// all flow values to the lower bound. deba@2526: /// \return This function returns false if the initialization deba@2526: /// process found a barrier. deba@2526: void init() deba@2526: { deba@2526: createStructures(); deba@2526: deba@2526: for(NodeIt n(_g);n!=INVALID;++n) { deba@2526: _excess->set(n, (*_delta)[n]); deba@2526: } deba@2526: deba@2526: for (EdgeIt e(_g);e!=INVALID;++e) { deba@2526: _flow->set(e, (*_lo)[e]); deba@2526: _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]); deba@2526: _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]); deba@2526: } deba@2526: deba@2526: typename Graph::template NodeMap reached(_g, false); deba@2526: deba@2526: deba@2526: // global relabeling tested, but in general case it provides deba@2526: // worse performance for random graphs deba@2526: _level->initStart(); deba@2526: for(NodeIt n(_g);n!=INVALID;++n) deba@2526: _level->initAddItem(n); deba@2526: _level->initFinish(); deba@2526: for(NodeIt n(_g);n!=INVALID;++n) deba@2526: if(_tol.positive((*_excess)[n])) deba@2526: _level->activate(n); deba@2526: } deba@2526: deba@2526: /// Initializes the internal data structures. deba@2526: deba@2526: /// Initializes the internal data structures. This functions uses deba@2526: /// greedy approach to construct the initial solution. deba@2526: void greedyInit() deba@2526: { deba@2526: createStructures(); deba@2526: deba@2526: for(NodeIt n(_g);n!=INVALID;++n) { deba@2526: _excess->set(n, (*_delta)[n]); deba@2526: } deba@2526: deba@2526: for (EdgeIt e(_g);e!=INVALID;++e) { deba@2526: if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) { deba@2526: _flow->set(e, (*_up)[e]); deba@2526: _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]); deba@2526: _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]); deba@2526: } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) { deba@2526: _flow->set(e, (*_lo)[e]); deba@2526: _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]); deba@2526: _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]); deba@2526: } else { deba@2526: Value fc = -(*_excess)[_g.target(e)]; deba@2526: _flow->set(e, fc); deba@2526: _excess->set(_g.target(e), 0); deba@2526: _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc); deba@2526: } deba@2526: } deba@2526: deba@2526: _level->initStart(); deba@2526: for(NodeIt n(_g);n!=INVALID;++n) deba@2526: _level->initAddItem(n); deba@2526: _level->initFinish(); deba@2526: for(NodeIt n(_g);n!=INVALID;++n) deba@2526: if(_tol.positive((*_excess)[n])) deba@2526: _level->activate(n); deba@2526: } deba@2526: deba@2526: ///Starts the algorithm deba@2526: deba@2526: ///This function starts the algorithm. deba@2526: ///\return This function returns true if it found a feasible circulation. deba@2526: /// deba@2526: ///\sa barrier() deba@2526: bool start() deba@2526: { deba@2526: deba@2526: Node act; deba@2526: Node bact=INVALID; deba@2526: Node last_activated=INVALID; deba@2526: while((act=_level->highestActive())!=INVALID) { deba@2526: int actlevel=(*_level)[act]; deba@2526: int mlevel=_node_num; deba@2526: Value exc=(*_excess)[act]; deba@2526: deba@2526: for(OutEdgeIt e(_g,act);e!=INVALID; ++e) { deba@2526: Node v = _g.target(e); deba@2526: Value fc=(*_up)[e]-(*_flow)[e]; deba@2526: if(!_tol.positive(fc)) continue; deba@2526: if((*_level)[v]set(e, (*_flow)[e] + exc); deba@2526: _excess->set(v, (*_excess)[v] + exc); deba@2526: if(!_level->active(v) && _tol.positive((*_excess)[v])) deba@2526: _level->activate(v); deba@2526: _excess->set(act,0); deba@2526: _level->deactivate(act); deba@2526: goto next_l; deba@2526: } deba@2526: else { deba@2526: _flow->set(e, (*_up)[e]); deba@2526: _excess->set(v, (*_excess)[v] + exc); deba@2526: if(!_level->active(v) && _tol.positive((*_excess)[v])) deba@2526: _level->activate(v); deba@2526: exc-=fc; deba@2526: } deba@2526: } deba@2526: else if((*_level)[v]set(e, (*_flow)[e] - exc); deba@2526: _excess->set(v, (*_excess)[v] + exc); deba@2526: if(!_level->active(v) && _tol.positive((*_excess)[v])) deba@2526: _level->activate(v); deba@2526: _excess->set(act,0); deba@2526: _level->deactivate(act); deba@2526: goto next_l; deba@2526: } deba@2526: else { deba@2526: _flow->set(e, (*_lo)[e]); deba@2526: _excess->set(v, (*_excess)[v] + fc); deba@2526: if(!_level->active(v) && _tol.positive((*_excess)[v])) deba@2526: _level->activate(v); deba@2526: exc-=fc; deba@2526: } deba@2526: } deba@2526: else if((*_level)[v]set(act, exc); deba@2526: if(!_tol.positive(exc)) _level->deactivate(act); deba@2526: else if(mlevel==_node_num) { deba@2526: _level->liftHighestActiveToTop(); deba@2526: _el = _node_num; deba@2526: return false; deba@2526: } deba@2526: else { deba@2526: _level->liftHighestActive(mlevel+1); deba@2526: if(_level->onLevel(actlevel)==0) { deba@2526: _el = actlevel; deba@2526: return false; deba@2526: } deba@2526: } deba@2526: next_l: deba@2526: ; deba@2526: } deba@2526: return true; deba@2526: } deba@2526: deba@2526: /// Runs the circulation algorithm. deba@2526: deba@2526: /// Runs the circulation algorithm. deba@2526: /// \note fc.run() is just a shortcut of the following code. deba@2526: /// \code deba@2526: /// fc.greedyInit(); deba@2526: /// return fc.start(); deba@2526: /// \endcode deba@2526: bool run() { deba@2526: greedyInit(); deba@2526: return start(); deba@2526: } deba@2526: deba@2526: /// @} deba@2526: deba@2526: /// \name Query Functions deba@2526: /// The result of the %Circulation algorithm can be obtained using deba@2526: /// these functions. deba@2526: /// \n deba@2526: /// Before the use of these functions, deba@2526: /// either run() or start() must be called. deba@2526: deba@2526: ///@{ deba@2526: deba@2526: ///Returns a barrier deba@2526: deba@2526: ///Barrier is a set \e B of nodes for which deba@2526: /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f] deba@2526: ///holds. The existence of a set with this property prooves that a feasible deba@2526: ///flow cannot exists. deba@2526: ///\sa checkBarrier() deba@2526: ///\sa run() deba@2526: template deba@2526: void barrierMap(GT &bar) deba@2526: { deba@2526: for(NodeIt n(_g);n!=INVALID;++n) deba@2526: bar.set(n, (*_level)[n] >= _el); deba@2526: } deba@2526: deba@2526: ///Returns true if the node is in the barrier deba@2526: deba@2526: ///Returns true if the node is in the barrier deba@2526: ///\sa barrierMap() deba@2526: bool barrier(const Node& node) deba@2526: { deba@2526: return (*_level)[node] >= _el; deba@2526: } deba@2526: deba@2526: /// \brief Returns the flow on the edge. deba@2526: /// deba@2526: /// Sets the \c flowMap to the flow on the edges. This method can deba@2526: /// be called after the second phase of algorithm. deba@2526: Value flow(const Edge& edge) const { deba@2526: return (*_flow)[edge]; deba@2526: } deba@2526: deba@2526: /// @} deba@2526: deba@2526: /// \name Checker Functions deba@2526: /// The feasibility of the results can be checked using deba@2526: /// these functions. deba@2526: /// \n deba@2526: /// Before the use of these functions, deba@2526: /// either run() or start() must be called. deba@2526: deba@2526: ///@{ deba@2526: deba@2526: ///Check if the \c flow is a feasible circulation deba@2526: bool checkFlow() { alpar@2375: for(EdgeIt e(_g);e!=INVALID;++e) deba@2526: if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false; alpar@2375: for(NodeIt n(_g);n!=INVALID;++n) alpar@2375: { deba@2526: Value dif=-(*_delta)[n]; deba@2526: for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e]; deba@2526: for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e]; alpar@2375: if(_tol.negative(dif)) return false; alpar@2375: } alpar@2375: return true; deba@2526: } alpar@2375: alpar@2408: ///Check whether or not the last execution provides a barrier alpar@2375: alpar@2408: ///Check whether or not the last execution provides a barrier alpar@2375: ///\sa barrier() alpar@2375: bool checkBarrier() alpar@2375: { deba@2526: Value delta=0; deba@2526: for(NodeIt n(_g);n!=INVALID;++n) deba@2526: if(barrier(n)) deba@2526: delta-=(*_delta)[n]; deba@2526: for(EdgeIt e(_g);e!=INVALID;++e) deba@2526: { deba@2526: Node s=_g.source(e); deba@2526: Node t=_g.target(e); deba@2526: if(barrier(s)&&!barrier(t)) delta+=(*_up)[e]; deba@2526: else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e]; deba@2526: } deba@2526: return _tol.negative(delta); alpar@2375: } alpar@2375: deba@2526: /// @} alpar@2375: alpar@2375: }; alpar@2375: alpar@2375: } alpar@2375: alpar@2375: #endif