1.1 --- a/lemon/circulation.h Wed Nov 28 17:40:41 2007 +0000
1.2 +++ b/lemon/circulation.h Wed Nov 28 17:51:02 2007 +0000
1.3 @@ -31,268 +31,612 @@
1.4 ///
1.5 namespace lemon {
1.6
1.7 - ///Preflow algorithm for the Network Circulation Problem.
1.8 + /// \brief Default traits class of Circulation class.
1.9 + ///
1.10 + /// Default traits class of Circulation class.
1.11 + /// \param _Graph Graph type.
1.12 + /// \param _CapacityMap Type of capacity map.
1.13 + template <typename _Graph, typename _LCapMap,
1.14 + typename _UCapMap, typename _DeltaMap>
1.15 + struct CirculationDefaultTraits {
1.16 +
1.17 + /// \brief The graph type the algorithm runs on.
1.18 + typedef _Graph Graph;
1.19 +
1.20 + /// \brief The type of the map that stores the circulation lower
1.21 + /// bound.
1.22 + ///
1.23 + /// The type of the map that stores the circulation lower bound.
1.24 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.25 + typedef _LCapMap LCapMap;
1.26 +
1.27 + /// \brief The type of the map that stores the circulation upper
1.28 + /// bound.
1.29 + ///
1.30 + /// The type of the map that stores the circulation upper bound.
1.31 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.32 + typedef _UCapMap UCapMap;
1.33 +
1.34 + /// \brief The type of the map that stores the upper bound of
1.35 + /// node excess.
1.36 + ///
1.37 + /// The type of the map that stores the lower bound of node
1.38 + /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
1.39 + /// concept.
1.40 + typedef _DeltaMap DeltaMap;
1.41 +
1.42 + /// \brief The type of the length of the edges.
1.43 + typedef typename DeltaMap::Value Value;
1.44 +
1.45 + /// \brief The map type that stores the flow values.
1.46 + ///
1.47 + /// The map type that stores the flow values.
1.48 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.49 + typedef typename Graph::template EdgeMap<Value> FlowMap;
1.50 +
1.51 + /// \brief Instantiates a FlowMap.
1.52 + ///
1.53 + /// This function instantiates a \ref FlowMap.
1.54 + /// \param graph The graph, to which we would like to define the flow map.
1.55 + static FlowMap* createFlowMap(const Graph& graph) {
1.56 + return new FlowMap(graph);
1.57 + }
1.58 +
1.59 + /// \brief The eleavator type used by Circulation algorithm.
1.60 + ///
1.61 + /// The elevator type used by Circulation algorithm.
1.62 + ///
1.63 + /// \sa Elevator
1.64 + /// \sa LinkedElevator
1.65 + typedef Elevator<Graph, typename Graph::Node> Elevator;
1.66 +
1.67 + /// \brief Instantiates an Elevator.
1.68 + ///
1.69 + /// This function instantiates a \ref Elevator.
1.70 + /// \param graph The graph, to which we would like to define the elevator.
1.71 + /// \param max_level The maximum level of the elevator.
1.72 + static Elevator* createElevator(const Graph& graph, int max_level) {
1.73 + return new Elevator(graph, max_level);
1.74 + }
1.75 +
1.76 + /// \brief The tolerance used by the algorithm
1.77 + ///
1.78 + /// The tolerance used by the algorithm to handle inexact computation.
1.79 + typedef Tolerance<Value> Tolerance;
1.80 +
1.81 + };
1.82 +
1.83 + ///Push-relabel algorithm for the Network Circulation Problem.
1.84
1.85 ///\ingroup max_flow
1.86 - ///This class implements a preflow algorithm
1.87 + ///This class implements a push-relabel algorithm
1.88 ///for the Network Circulation Problem.
1.89 ///The exact formulation of this problem is the following.
1.90 /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq -delta(v)\quad \forall v\in V \f]
1.91 /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
1.92 ///
1.93 - template<class Graph,
1.94 - class Value,
1.95 - class FlowMap=typename Graph::template EdgeMap<Value>,
1.96 - class LCapMap=typename Graph::template EdgeMap<Value>,
1.97 - class UCapMap=LCapMap,
1.98 - class DeltaMap=typename Graph::template NodeMap<Value>
1.99 - >
1.100 + template<class _Graph,
1.101 + class _LCapMap=typename _Graph::template EdgeMap<int>,
1.102 + class _UCapMap=_LCapMap,
1.103 + class _DeltaMap=typename _Graph::template NodeMap<
1.104 + typename _UCapMap::Value>,
1.105 + class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
1.106 + _UCapMap, _DeltaMap> >
1.107 class Circulation {
1.108 - typedef typename Graph::Node Node;
1.109 - typedef typename Graph::NodeIt NodeIt;
1.110 - typedef typename Graph::Edge Edge;
1.111 - typedef typename Graph::EdgeIt EdgeIt;
1.112 - typedef typename Graph::InEdgeIt InEdgeIt;
1.113 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.114 -
1.115 +
1.116 + typedef _Traits Traits;
1.117 + typedef typename Traits::Graph Graph;
1.118 + GRAPH_TYPEDEFS(typename Graph);
1.119 +
1.120 + typedef typename Traits::Value Value;
1.121 +
1.122 + typedef typename Traits::LCapMap LCapMap;
1.123 + typedef typename Traits::UCapMap UCapMap;
1.124 + typedef typename Traits::DeltaMap DeltaMap;
1.125 + typedef typename Traits::FlowMap FlowMap;
1.126 + typedef typename Traits::Elevator Elevator;
1.127 + typedef typename Traits::Tolerance Tolerance;
1.128 +
1.129 + typedef typename Graph::template NodeMap<Value> ExcessMap;
1.130
1.131 const Graph &_g;
1.132 int _node_num;
1.133
1.134 - const LCapMap &_lo;
1.135 - const UCapMap &_up;
1.136 - const DeltaMap &_delta;
1.137 - FlowMap &_x;
1.138 - Tolerance<Value> _tol;
1.139 - Elevator<Graph,typename Graph::Node> _levels;
1.140 - typename Graph::template NodeMap<Value> _excess;
1.141 + const LCapMap *_lo;
1.142 + const UCapMap *_up;
1.143 + const DeltaMap *_delta;
1.144 +
1.145 + FlowMap *_flow;
1.146 + bool _local_flow;
1.147 +
1.148 + Elevator* _level;
1.149 + bool _local_level;
1.150 +
1.151 + ExcessMap* _excess;
1.152 +
1.153 + Tolerance _tol;
1.154 + int _el;
1.155
1.156 public:
1.157 - ///\e
1.158 - Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
1.159 - const DeltaMap &delta,FlowMap &x) :
1.160 - _g(g),
1.161 - _node_num(countNodes(g)),
1.162 - _lo(lo),_up(up),_delta(delta),_x(x),
1.163 - _levels(g,_node_num),
1.164 - _excess(g)
1.165 - {
1.166 +
1.167 + typedef Circulation Create;
1.168 +
1.169 + ///\name Named template parameters
1.170 +
1.171 + ///@{
1.172 +
1.173 + template <typename _FlowMap>
1.174 + struct DefFlowMapTraits : public Traits {
1.175 + typedef _FlowMap FlowMap;
1.176 + static FlowMap *createFlowMap(const Graph&) {
1.177 + throw UninitializedParameter();
1.178 + }
1.179 + };
1.180 +
1.181 + /// \brief \ref named-templ-param "Named parameter" for setting
1.182 + /// FlowMap type
1.183 + ///
1.184 + /// \ref named-templ-param "Named parameter" for setting FlowMap
1.185 + /// type
1.186 + template <typename _FlowMap>
1.187 + struct DefFlowMap
1.188 + : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
1.189 + DefFlowMapTraits<_FlowMap> > {
1.190 + typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
1.191 + DefFlowMapTraits<_FlowMap> > Create;
1.192 + };
1.193 +
1.194 + template <typename _Elevator>
1.195 + struct DefElevatorTraits : public Traits {
1.196 + typedef _Elevator Elevator;
1.197 + static Elevator *createElevator(const Graph&, int) {
1.198 + throw UninitializedParameter();
1.199 + }
1.200 + };
1.201 +
1.202 + /// \brief \ref named-templ-param "Named parameter" for setting
1.203 + /// Elevator type
1.204 + ///
1.205 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.206 + /// type
1.207 + template <typename _Elevator>
1.208 + struct DefElevator
1.209 + : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
1.210 + DefElevatorTraits<_Elevator> > {
1.211 + typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
1.212 + DefElevatorTraits<_Elevator> > Create;
1.213 + };
1.214 +
1.215 + template <typename _Elevator>
1.216 + struct DefStandardElevatorTraits : public Traits {
1.217 + typedef _Elevator Elevator;
1.218 + static Elevator *createElevator(const Graph& graph, int max_level) {
1.219 + return new Elevator(graph, max_level);
1.220 + }
1.221 + };
1.222 +
1.223 + /// \brief \ref named-templ-param "Named parameter" for setting
1.224 + /// Elevator type
1.225 + ///
1.226 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.227 + /// type. The Elevator should be standard constructor interface, ie.
1.228 + /// the graph and the maximum level should be passed to it.
1.229 + template <typename _Elevator>
1.230 + struct DefStandardElevator
1.231 + : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
1.232 + DefStandardElevatorTraits<_Elevator> > {
1.233 + typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
1.234 + DefStandardElevatorTraits<_Elevator> > Create;
1.235 + };
1.236 +
1.237 + /// @}
1.238 +
1.239 + /// The constructor of the class.
1.240 +
1.241 + /// The constructor of the class.
1.242 + /// \param g The directed graph the algorithm runs on.
1.243 + /// \param lo The lower bound capacity of the edges.
1.244 + /// \param up The upper bound capacity of the edges.
1.245 + /// \param delta The lower bound on node excess.
1.246 + Circulation(const Graph &g,const LCapMap &lo,
1.247 + const UCapMap &up,const DeltaMap &delta)
1.248 + : _g(g), _node_num(),
1.249 + _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
1.250 + _level(0), _local_level(false), _excess(0), _el() {}
1.251 +
1.252 + /// Destrcutor.
1.253 + ~Circulation() {
1.254 + destroyStructures();
1.255 }
1.256 -
1.257 +
1.258 private:
1.259
1.260 - void addExcess(Node n,Value v)
1.261 - {
1.262 - if(_tol.positive(_excess[n]+=v))
1.263 - {
1.264 - if(!_levels.active(n)) _levels.activate(n);
1.265 - }
1.266 - else if(_levels.active(n)) _levels.deactivate(n);
1.267 + void createStructures() {
1.268 + _node_num = _el = countNodes(_g);
1.269 +
1.270 + if (!_flow) {
1.271 + _flow = Traits::createFlowMap(_g);
1.272 + _local_flow = true;
1.273 + }
1.274 + if (!_level) {
1.275 + _level = Traits::createElevator(_g, _node_num);
1.276 + _local_level = true;
1.277 + }
1.278 + if (!_excess) {
1.279 + _excess = new ExcessMap(_g);
1.280 + }
1.281 }
1.282 -
1.283 - void init()
1.284 - {
1.285 -
1.286 - _x=_lo;
1.287
1.288 - for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=_delta[n];
1.289 -
1.290 - for(EdgeIt e(_g);e!=INVALID;++e)
1.291 - {
1.292 - _excess[_g.target(e)]+=_x[e];
1.293 - _excess[_g.source(e)]-=_x[e];
1.294 - }
1.295 -
1.296 - _levels.initStart();
1.297 - for(NodeIt n(_g);n!=INVALID;++n)
1.298 - _levels.initAddItem(n);
1.299 - _levels.initFinish();
1.300 - for(NodeIt n(_g);n!=INVALID;++n)
1.301 - if(_tol.positive(_excess[n]))
1.302 - _levels.activate(n);
1.303 + void destroyStructures() {
1.304 + if (_local_flow) {
1.305 + delete _flow;
1.306 + }
1.307 + if (_local_level) {
1.308 + delete _level;
1.309 + }
1.310 + if (_excess) {
1.311 + delete _excess;
1.312 + }
1.313 }
1.314
1.315 public:
1.316 - ///Check if \c x is a feasible circulation
1.317 - template<class FT>
1.318 - bool checkX(FT &x) {
1.319 +
1.320 + /// Sets the lower bound capacity map.
1.321 +
1.322 + /// Sets the lower bound capacity map.
1.323 + /// \return \c (*this)
1.324 + Circulation& lowerCapMap(const LCapMap& map) {
1.325 + _lo = ↦
1.326 + return *this;
1.327 + }
1.328 +
1.329 + /// Sets the upper bound capacity map.
1.330 +
1.331 + /// Sets the upper bound capacity map.
1.332 + /// \return \c (*this)
1.333 + Circulation& upperCapMap(const LCapMap& map) {
1.334 + _up = ↦
1.335 + return *this;
1.336 + }
1.337 +
1.338 + /// Sets the lower bound map on excess.
1.339 +
1.340 + /// Sets the lower bound map on excess.
1.341 + /// \return \c (*this)
1.342 + Circulation& deltaMap(const DeltaMap& map) {
1.343 + _delta = ↦
1.344 + return *this;
1.345 + }
1.346 +
1.347 + /// Sets the flow map.
1.348 +
1.349 + /// Sets the flow map.
1.350 + /// \return \c (*this)
1.351 + Circulation& flowMap(FlowMap& map) {
1.352 + if (_local_flow) {
1.353 + delete _flow;
1.354 + _local_flow = false;
1.355 + }
1.356 + _flow = ↦
1.357 + return *this;
1.358 + }
1.359 +
1.360 + /// Returns the flow map.
1.361 +
1.362 + /// \return The flow map.
1.363 + ///
1.364 + const FlowMap& flowMap() {
1.365 + return *_flow;
1.366 + }
1.367 +
1.368 + /// Sets the elevator.
1.369 +
1.370 + /// Sets the elevator.
1.371 + /// \return \c (*this)
1.372 + Circulation& elevator(Elevator& elevator) {
1.373 + if (_local_level) {
1.374 + delete _level;
1.375 + _local_level = false;
1.376 + }
1.377 + _level = &elevator;
1.378 + return *this;
1.379 + }
1.380 +
1.381 + /// Returns the elevator.
1.382 +
1.383 + /// \return The elevator.
1.384 + ///
1.385 + const Elevator& elevator() {
1.386 + return *_level;
1.387 + }
1.388 +
1.389 + /// Sets the tolerance used by algorithm.
1.390 +
1.391 + /// Sets the tolerance used by algorithm.
1.392 + ///
1.393 + Circulation& tolerance(const Tolerance& tolerance) const {
1.394 + _tol = tolerance;
1.395 + return *this;
1.396 + }
1.397 +
1.398 + /// Returns the tolerance used by algorithm.
1.399 +
1.400 + /// Returns the tolerance used by algorithm.
1.401 + ///
1.402 + const Tolerance& tolerance() const {
1.403 + return tolerance;
1.404 + }
1.405 +
1.406 + /// \name Execution control The simplest way to execute the
1.407 + /// algorithm is to use one of the member functions called \c
1.408 + /// run().
1.409 + /// \n
1.410 + /// If you need more control on initial solution or execution then
1.411 + /// you have to call one \ref init() function and then the start()
1.412 + /// function.
1.413 +
1.414 + ///@{
1.415 +
1.416 + /// Initializes the internal data structures.
1.417 +
1.418 + /// Initializes the internal data structures. This function sets
1.419 + /// all flow values to the lower bound.
1.420 + /// \return This function returns false if the initialization
1.421 + /// process found a barrier.
1.422 + void init()
1.423 + {
1.424 + createStructures();
1.425 +
1.426 + for(NodeIt n(_g);n!=INVALID;++n) {
1.427 + _excess->set(n, (*_delta)[n]);
1.428 + }
1.429 +
1.430 + for (EdgeIt e(_g);e!=INVALID;++e) {
1.431 + _flow->set(e, (*_lo)[e]);
1.432 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
1.433 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
1.434 + }
1.435 +
1.436 + typename Graph::template NodeMap<bool> reached(_g, false);
1.437 +
1.438 +
1.439 + // global relabeling tested, but in general case it provides
1.440 + // worse performance for random graphs
1.441 + _level->initStart();
1.442 + for(NodeIt n(_g);n!=INVALID;++n)
1.443 + _level->initAddItem(n);
1.444 + _level->initFinish();
1.445 + for(NodeIt n(_g);n!=INVALID;++n)
1.446 + if(_tol.positive((*_excess)[n]))
1.447 + _level->activate(n);
1.448 + }
1.449 +
1.450 + /// Initializes the internal data structures.
1.451 +
1.452 + /// Initializes the internal data structures. This functions uses
1.453 + /// greedy approach to construct the initial solution.
1.454 + void greedyInit()
1.455 + {
1.456 + createStructures();
1.457 +
1.458 + for(NodeIt n(_g);n!=INVALID;++n) {
1.459 + _excess->set(n, (*_delta)[n]);
1.460 + }
1.461 +
1.462 + for (EdgeIt e(_g);e!=INVALID;++e) {
1.463 + if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
1.464 + _flow->set(e, (*_up)[e]);
1.465 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
1.466 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
1.467 + } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
1.468 + _flow->set(e, (*_lo)[e]);
1.469 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
1.470 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
1.471 + } else {
1.472 + Value fc = -(*_excess)[_g.target(e)];
1.473 + _flow->set(e, fc);
1.474 + _excess->set(_g.target(e), 0);
1.475 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
1.476 + }
1.477 + }
1.478 +
1.479 + _level->initStart();
1.480 + for(NodeIt n(_g);n!=INVALID;++n)
1.481 + _level->initAddItem(n);
1.482 + _level->initFinish();
1.483 + for(NodeIt n(_g);n!=INVALID;++n)
1.484 + if(_tol.positive((*_excess)[n]))
1.485 + _level->activate(n);
1.486 + }
1.487 +
1.488 + ///Starts the algorithm
1.489 +
1.490 + ///This function starts the algorithm.
1.491 + ///\return This function returns true if it found a feasible circulation.
1.492 + ///
1.493 + ///\sa barrier()
1.494 + bool start()
1.495 + {
1.496 +
1.497 + Node act;
1.498 + Node bact=INVALID;
1.499 + Node last_activated=INVALID;
1.500 + while((act=_level->highestActive())!=INVALID) {
1.501 + int actlevel=(*_level)[act];
1.502 + int mlevel=_node_num;
1.503 + Value exc=(*_excess)[act];
1.504 +
1.505 + for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
1.506 + Node v = _g.target(e);
1.507 + Value fc=(*_up)[e]-(*_flow)[e];
1.508 + if(!_tol.positive(fc)) continue;
1.509 + if((*_level)[v]<actlevel) {
1.510 + if(!_tol.less(fc, exc)) {
1.511 + _flow->set(e, (*_flow)[e] + exc);
1.512 + _excess->set(v, (*_excess)[v] + exc);
1.513 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.514 + _level->activate(v);
1.515 + _excess->set(act,0);
1.516 + _level->deactivate(act);
1.517 + goto next_l;
1.518 + }
1.519 + else {
1.520 + _flow->set(e, (*_up)[e]);
1.521 + _excess->set(v, (*_excess)[v] + exc);
1.522 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.523 + _level->activate(v);
1.524 + exc-=fc;
1.525 + }
1.526 + }
1.527 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
1.528 + }
1.529 + for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
1.530 + Node v = _g.source(e);
1.531 + Value fc=(*_flow)[e]-(*_lo)[e];
1.532 + if(!_tol.positive(fc)) continue;
1.533 + if((*_level)[v]<actlevel) {
1.534 + if(!_tol.less(fc, exc)) {
1.535 + _flow->set(e, (*_flow)[e] - exc);
1.536 + _excess->set(v, (*_excess)[v] + exc);
1.537 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.538 + _level->activate(v);
1.539 + _excess->set(act,0);
1.540 + _level->deactivate(act);
1.541 + goto next_l;
1.542 + }
1.543 + else {
1.544 + _flow->set(e, (*_lo)[e]);
1.545 + _excess->set(v, (*_excess)[v] + fc);
1.546 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.547 + _level->activate(v);
1.548 + exc-=fc;
1.549 + }
1.550 + }
1.551 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
1.552 + }
1.553 +
1.554 + _excess->set(act, exc);
1.555 + if(!_tol.positive(exc)) _level->deactivate(act);
1.556 + else if(mlevel==_node_num) {
1.557 + _level->liftHighestActiveToTop();
1.558 + _el = _node_num;
1.559 + return false;
1.560 + }
1.561 + else {
1.562 + _level->liftHighestActive(mlevel+1);
1.563 + if(_level->onLevel(actlevel)==0) {
1.564 + _el = actlevel;
1.565 + return false;
1.566 + }
1.567 + }
1.568 + next_l:
1.569 + ;
1.570 + }
1.571 + return true;
1.572 + }
1.573 +
1.574 + /// Runs the circulation algorithm.
1.575 +
1.576 + /// Runs the circulation algorithm.
1.577 + /// \note fc.run() is just a shortcut of the following code.
1.578 + /// \code
1.579 + /// fc.greedyInit();
1.580 + /// return fc.start();
1.581 + /// \endcode
1.582 + bool run() {
1.583 + greedyInit();
1.584 + return start();
1.585 + }
1.586 +
1.587 + /// @}
1.588 +
1.589 + /// \name Query Functions
1.590 + /// The result of the %Circulation algorithm can be obtained using
1.591 + /// these functions.
1.592 + /// \n
1.593 + /// Before the use of these functions,
1.594 + /// either run() or start() must be called.
1.595 +
1.596 + ///@{
1.597 +
1.598 + ///Returns a barrier
1.599 +
1.600 + ///Barrier is a set \e B of nodes for which
1.601 + /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
1.602 + ///holds. The existence of a set with this property prooves that a feasible
1.603 + ///flow cannot exists.
1.604 + ///\sa checkBarrier()
1.605 + ///\sa run()
1.606 + template<class GT>
1.607 + void barrierMap(GT &bar)
1.608 + {
1.609 + for(NodeIt n(_g);n!=INVALID;++n)
1.610 + bar.set(n, (*_level)[n] >= _el);
1.611 + }
1.612 +
1.613 + ///Returns true if the node is in the barrier
1.614 +
1.615 + ///Returns true if the node is in the barrier
1.616 + ///\sa barrierMap()
1.617 + bool barrier(const Node& node)
1.618 + {
1.619 + return (*_level)[node] >= _el;
1.620 + }
1.621 +
1.622 + /// \brief Returns the flow on the edge.
1.623 + ///
1.624 + /// Sets the \c flowMap to the flow on the edges. This method can
1.625 + /// be called after the second phase of algorithm.
1.626 + Value flow(const Edge& edge) const {
1.627 + return (*_flow)[edge];
1.628 + }
1.629 +
1.630 + /// @}
1.631 +
1.632 + /// \name Checker Functions
1.633 + /// The feasibility of the results can be checked using
1.634 + /// these functions.
1.635 + /// \n
1.636 + /// Before the use of these functions,
1.637 + /// either run() or start() must be called.
1.638 +
1.639 + ///@{
1.640 +
1.641 + ///Check if the \c flow is a feasible circulation
1.642 + bool checkFlow() {
1.643 for(EdgeIt e(_g);e!=INVALID;++e)
1.644 - if(x[e]<_lo[e]||x[e]>_up[e]) return false;
1.645 + if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
1.646 for(NodeIt n(_g);n!=INVALID;++n)
1.647 {
1.648 - Value dif=-_delta[n];
1.649 - for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];
1.650 - for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];
1.651 + Value dif=-(*_delta)[n];
1.652 + for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
1.653 + for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
1.654 if(_tol.negative(dif)) return false;
1.655 }
1.656 return true;
1.657 - };
1.658 - ///Check if the default \c x is a feasible circulation
1.659 - bool checkX() { return checkX(_x); }
1.660 + }
1.661
1.662 - ///Check if \c bar is a real barrier
1.663 -
1.664 - ///Check if \c bar is a real barrier
1.665 - ///\sa barrier()
1.666 - template<class GT>
1.667 - bool checkBarrier(GT &bar)
1.668 - {
1.669 - Value delta=0;
1.670 - for(NodeIt n(_g);n!=INVALID;++n)
1.671 - if(bar[n])
1.672 - delta-=_delta[n];
1.673 - for(EdgeIt e(_g);e!=INVALID;++e)
1.674 - {
1.675 - Node s=_g.source(e);
1.676 - Node t=_g.target(e);
1.677 - if(bar[s]&&!bar[t]) delta+=_up[e];
1.678 - else if(bar[t]&&!bar[s]) delta-=_lo[e];
1.679 - }
1.680 - return _tol.negative(delta);
1.681 - }
1.682 ///Check whether or not the last execution provides a barrier
1.683
1.684 ///Check whether or not the last execution provides a barrier
1.685 ///\sa barrier()
1.686 bool checkBarrier()
1.687 {
1.688 - typename Graph:: template NodeMap<bool> bar(_g);
1.689 - barrier(bar);
1.690 - return checkBarrier(bar);
1.691 + Value delta=0;
1.692 + for(NodeIt n(_g);n!=INVALID;++n)
1.693 + if(barrier(n))
1.694 + delta-=(*_delta)[n];
1.695 + for(EdgeIt e(_g);e!=INVALID;++e)
1.696 + {
1.697 + Node s=_g.source(e);
1.698 + Node t=_g.target(e);
1.699 + if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
1.700 + else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
1.701 + }
1.702 + return _tol.negative(delta);
1.703 }
1.704 - ///Run the algorithm
1.705
1.706 - ///This function runs the algorithm.
1.707 - ///\return This function returns -1 if it found a feasible circulation.
1.708 - /// nonnegative values (including 0) mean that no feasible solution is
1.709 - /// found. In this case the return value means an "empty level".
1.710 - ///
1.711 - ///\sa barrier()
1.712 - int run()
1.713 - {
1.714 - init();
1.715 -
1.716 -#ifdef LEMON_CIRCULATION_DEBUG
1.717 - for(NodeIt n(_g);n!=INVALID;++n)
1.718 - std::cerr<< _levels[n] << ' ';
1.719 - std::cerr << std::endl;
1.720 -#endif
1.721 - Node act;
1.722 - Node bact=INVALID;
1.723 - Node last_activated=INVALID;
1.724 - while((act=_levels.highestActive())!=INVALID) {
1.725 - int actlevel=_levels[act];
1.726 - int tlevel;
1.727 - int mlevel=_node_num;
1.728 - Value exc=_excess[act];
1.729 - Value fc;
1.730 -
1.731 -#ifdef LEMON_CIRCULATION_DEBUG
1.732 - for(NodeIt n(_g);n!=INVALID;++n)
1.733 - std::cerr<< _levels[n] << ' ';
1.734 - std::cerr << std::endl;
1.735 - std::cerr << "Process node " << _g.id(act)
1.736 - << " on level " << actlevel
1.737 - << " with excess " << exc
1.738 - << std::endl;
1.739 -#endif
1.740 - for(OutEdgeIt e(_g,act);e!=INVALID; ++e)
1.741 - if((fc=_up[e]-_x[e])>0)
1.742 - if((tlevel=_levels[_g.target(e)])<actlevel)
1.743 - if(fc<=exc) {
1.744 - _x[e]=_up[e];
1.745 - addExcess(_g.target(e),fc);
1.746 - exc-=fc;
1.747 -#ifdef LEMON_CIRCULATION_DEBUG
1.748 - std::cerr << " Push " << fc
1.749 - << " toward " << _g.id(_g.target(e)) << std::endl;
1.750 -#endif
1.751 - }
1.752 - else {
1.753 - _x[e]+=exc;
1.754 - addExcess(_g.target(e),exc);
1.755 - //exc=0;
1.756 - _excess[act]=0;
1.757 - _levels.deactivate(act);
1.758 -#ifdef LEMON_CIRCULATION_DEBUG
1.759 - std::cerr << " Push " << exc
1.760 - << " toward " << _g.id(_g.target(e)) << std::endl;
1.761 - std::cerr << " Deactivate\n";
1.762 -#endif
1.763 - goto next_l;
1.764 - }
1.765 - else if(tlevel<mlevel) mlevel=tlevel;
1.766 -
1.767 - for(InEdgeIt e(_g,act);e!=INVALID; ++e)
1.768 - if((fc=_x[e]-_lo[e])>0)
1.769 - if((tlevel=_levels[_g.source(e)])<actlevel)
1.770 - if(fc<=exc) {
1.771 - _x[e]=_lo[e];
1.772 - addExcess(_g.source(e),fc);
1.773 - exc-=fc;
1.774 -#ifdef LEMON_CIRCULATION_DEBUG
1.775 - std::cerr << " Push " << fc
1.776 - << " toward " << _g.id(_g.source(e)) << std::endl;
1.777 -#endif
1.778 - }
1.779 - else {
1.780 - _x[e]-=exc;
1.781 - addExcess(_g.source(e),exc);
1.782 - //exc=0;
1.783 - _excess[act]=0;
1.784 - _levels.deactivate(act);
1.785 -#ifdef LEMON_CIRCULATION_DEBUG
1.786 - std::cerr << " Push " << exc
1.787 - << " toward " << _g.id(_g.source(e)) << std::endl;
1.788 - std::cerr << " Deactivate\n";
1.789 -#endif
1.790 - goto next_l;
1.791 - }
1.792 - else if(tlevel<mlevel) mlevel=tlevel;
1.793 -
1.794 - _excess[act]=exc;
1.795 - if(!_tol.positive(exc)) _levels.deactivate(act);
1.796 - else if(mlevel==_node_num) {
1.797 - _levels.liftHighestActiveToTop();
1.798 -#ifdef LEMON_CIRCULATION_DEBUG
1.799 - std::cerr << " Lift to level " << _node_num << std::endl;
1.800 -#endif
1.801 - return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
1.802 - }
1.803 - else {
1.804 - _levels.liftHighestActive(mlevel+1);
1.805 -#ifdef LEMON_CIRCULATION_DEBUG
1.806 - std::cerr << " Lift to level " << mlevel+1 << std::endl;
1.807 -#endif
1.808 - if(_levels.onLevel(actlevel)==0)
1.809 - return actlevel;
1.810 - }
1.811 - next_l:
1.812 - ;
1.813 - }
1.814 -#ifdef LEMON_CIRCULATION_DEBUG
1.815 - std::cerr << "Feasible flow found.\n";
1.816 -#endif
1.817 - return -1;
1.818 - }
1.819 -
1.820 - ///Return a barrier
1.821 -
1.822 - ///Barrier is a set \e B of nodes for which
1.823 - /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
1.824 - ///holds. The existence of a set with this property prooves that a feasible
1.825 - ///flow cannot exists.
1.826 - ///\pre The run() must have been executed, and its return value was -1.
1.827 - ///\sa checkBarrier()
1.828 - ///\sa run()
1.829 - template<class GT>
1.830 - void barrier(GT &bar,int empty_level=-1)
1.831 - {
1.832 - if(empty_level==-1)
1.833 - for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
1.834 - for(NodeIt n(_g);n!=INVALID;++n)
1.835 - bar[n] = _levels[n]>empty_level;
1.836 - }
1.837 + /// @}
1.838
1.839 };
1.840