1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/circulation.h Mon Dec 01 14:18:40 2008 +0000
1.3 @@ -0,0 +1,755 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2008
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_CIRCULATION_H
1.23 +#define LEMON_CIRCULATION_H
1.24 +
1.25 +#include <lemon/tolerance.h>
1.26 +#include <lemon/elevator.h>
1.27 +
1.28 +///\ingroup max_flow
1.29 +///\file
1.30 +///\brief Push-relabel algorithm for finding a feasible circulation.
1.31 +///
1.32 +namespace lemon {
1.33 +
1.34 + /// \brief Default traits class of Circulation class.
1.35 + ///
1.36 + /// Default traits class of Circulation class.
1.37 + /// \tparam _Diraph Digraph type.
1.38 + /// \tparam _LCapMap Lower bound capacity map type.
1.39 + /// \tparam _UCapMap Upper bound capacity map type.
1.40 + /// \tparam _DeltaMap Delta map type.
1.41 + template <typename _Diraph, typename _LCapMap,
1.42 + typename _UCapMap, typename _DeltaMap>
1.43 + struct CirculationDefaultTraits {
1.44 +
1.45 + /// \brief The type of the digraph the algorithm runs on.
1.46 + typedef _Diraph Digraph;
1.47 +
1.48 + /// \brief The type of the map that stores the circulation lower
1.49 + /// bound.
1.50 + ///
1.51 + /// The type of the map that stores the circulation lower bound.
1.52 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.53 + typedef _LCapMap LCapMap;
1.54 +
1.55 + /// \brief The type of the map that stores the circulation upper
1.56 + /// bound.
1.57 + ///
1.58 + /// The type of the map that stores the circulation upper bound.
1.59 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.60 + typedef _UCapMap UCapMap;
1.61 +
1.62 + /// \brief The type of the map that stores the lower bound for
1.63 + /// the supply of the nodes.
1.64 + ///
1.65 + /// The type of the map that stores the lower bound for the supply
1.66 + /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
1.67 + /// concept.
1.68 + typedef _DeltaMap DeltaMap;
1.69 +
1.70 + /// \brief The type of the flow values.
1.71 + typedef typename DeltaMap::Value Value;
1.72 +
1.73 + /// \brief The type of the map that stores the flow values.
1.74 + ///
1.75 + /// The type of the map that stores the flow values.
1.76 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.77 + typedef typename Digraph::template ArcMap<Value> FlowMap;
1.78 +
1.79 + /// \brief Instantiates a FlowMap.
1.80 + ///
1.81 + /// This function instantiates a \ref FlowMap.
1.82 + /// \param digraph The digraph, to which we would like to define
1.83 + /// the flow map.
1.84 + static FlowMap* createFlowMap(const Digraph& digraph) {
1.85 + return new FlowMap(digraph);
1.86 + }
1.87 +
1.88 + /// \brief The elevator type used by the algorithm.
1.89 + ///
1.90 + /// The elevator type used by the algorithm.
1.91 + ///
1.92 + /// \sa Elevator
1.93 + /// \sa LinkedElevator
1.94 + typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
1.95 +
1.96 + /// \brief Instantiates an Elevator.
1.97 + ///
1.98 + /// This function instantiates an \ref Elevator.
1.99 + /// \param digraph The digraph, to which we would like to define
1.100 + /// the elevator.
1.101 + /// \param max_level The maximum level of the elevator.
1.102 + static Elevator* createElevator(const Digraph& digraph, int max_level) {
1.103 + return new Elevator(digraph, max_level);
1.104 + }
1.105 +
1.106 + /// \brief The tolerance used by the algorithm
1.107 + ///
1.108 + /// The tolerance used by the algorithm to handle inexact computation.
1.109 + typedef lemon::Tolerance<Value> Tolerance;
1.110 +
1.111 + };
1.112 +
1.113 + /**
1.114 + \brief Push-relabel algorithm for the network circulation problem.
1.115 +
1.116 + \ingroup max_flow
1.117 + This class implements a push-relabel algorithm for the network
1.118 + circulation problem.
1.119 + It is to find a feasible circulation when lower and upper bounds
1.120 + are given for the flow values on the arcs and lower bounds
1.121 + are given for the supply values of the nodes.
1.122 +
1.123 + The exact formulation of this problem is the following.
1.124 + Let \f$G=(V,A)\f$ be a digraph,
1.125 + \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
1.126 + \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
1.127 + \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
1.128 + \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
1.129 + \geq delta(v) \quad \forall v\in V, \f]
1.130 + \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
1.131 + \note \f$delta(v)\f$ specifies a lower bound for the supply of node
1.132 + \f$v\f$. It can be either positive or negative, however note that
1.133 + \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
1.134 + have a feasible solution.
1.135 +
1.136 + \note A special case of this problem is when
1.137 + \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
1.138 + will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
1.139 + Thus a feasible solution for the
1.140 + \ref min_cost_flow "minimum cost flow" problem can be calculated
1.141 + in this way.
1.142 +
1.143 + \tparam _Digraph The type of the digraph the algorithm runs on.
1.144 + \tparam _LCapMap The type of the lower bound capacity map. The default
1.145 + map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
1.146 + \tparam _UCapMap The type of the upper bound capacity map. The default
1.147 + map type is \c _LCapMap.
1.148 + \tparam _DeltaMap The type of the map that stores the lower bound
1.149 + for the supply of the nodes. The default map type is
1.150 + \c _Digraph::ArcMap<_UCapMap::Value>.
1.151 + */
1.152 +#ifdef DOXYGEN
1.153 +template< typename _Digraph,
1.154 + typename _LCapMap,
1.155 + typename _UCapMap,
1.156 + typename _DeltaMap,
1.157 + typename _Traits >
1.158 +#else
1.159 +template< typename _Digraph,
1.160 + typename _LCapMap = typename _Digraph::template ArcMap<int>,
1.161 + typename _UCapMap = _LCapMap,
1.162 + typename _DeltaMap = typename _Digraph::
1.163 + template NodeMap<typename _UCapMap::Value>,
1.164 + typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
1.165 + _UCapMap, _DeltaMap> >
1.166 +#endif
1.167 + class Circulation {
1.168 + public:
1.169 +
1.170 + ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
1.171 + typedef _Traits Traits;
1.172 + ///The type of the digraph the algorithm runs on.
1.173 + typedef typename Traits::Digraph Digraph;
1.174 + ///The type of the flow values.
1.175 + typedef typename Traits::Value Value;
1.176 +
1.177 + /// The type of the lower bound capacity map.
1.178 + typedef typename Traits::LCapMap LCapMap;
1.179 + /// The type of the upper bound capacity map.
1.180 + typedef typename Traits::UCapMap UCapMap;
1.181 + /// \brief The type of the map that stores the lower bound for
1.182 + /// the supply of the nodes.
1.183 + typedef typename Traits::DeltaMap DeltaMap;
1.184 + ///The type of the flow map.
1.185 + typedef typename Traits::FlowMap FlowMap;
1.186 +
1.187 + ///The type of the elevator.
1.188 + typedef typename Traits::Elevator Elevator;
1.189 + ///The type of the tolerance.
1.190 + typedef typename Traits::Tolerance Tolerance;
1.191 +
1.192 + private:
1.193 +
1.194 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.195 +
1.196 + const Digraph &_g;
1.197 + int _node_num;
1.198 +
1.199 + const LCapMap *_lo;
1.200 + const UCapMap *_up;
1.201 + const DeltaMap *_delta;
1.202 +
1.203 + FlowMap *_flow;
1.204 + bool _local_flow;
1.205 +
1.206 + Elevator* _level;
1.207 + bool _local_level;
1.208 +
1.209 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
1.210 + ExcessMap* _excess;
1.211 +
1.212 + Tolerance _tol;
1.213 + int _el;
1.214 +
1.215 + public:
1.216 +
1.217 + typedef Circulation Create;
1.218 +
1.219 + ///\name Named Template Parameters
1.220 +
1.221 + ///@{
1.222 +
1.223 + template <typename _FlowMap>
1.224 + struct SetFlowMapTraits : public Traits {
1.225 + typedef _FlowMap FlowMap;
1.226 + static FlowMap *createFlowMap(const Digraph&) {
1.227 + LEMON_ASSERT(false, "FlowMap is not initialized");
1.228 + return 0; // ignore warnings
1.229 + }
1.230 + };
1.231 +
1.232 + /// \brief \ref named-templ-param "Named parameter" for setting
1.233 + /// FlowMap type
1.234 + ///
1.235 + /// \ref named-templ-param "Named parameter" for setting FlowMap
1.236 + /// type.
1.237 + template <typename _FlowMap>
1.238 + struct SetFlowMap
1.239 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.240 + SetFlowMapTraits<_FlowMap> > {
1.241 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.242 + SetFlowMapTraits<_FlowMap> > Create;
1.243 + };
1.244 +
1.245 + template <typename _Elevator>
1.246 + struct SetElevatorTraits : public Traits {
1.247 + typedef _Elevator Elevator;
1.248 + static Elevator *createElevator(const Digraph&, int) {
1.249 + LEMON_ASSERT(false, "Elevator is not initialized");
1.250 + return 0; // ignore warnings
1.251 + }
1.252 + };
1.253 +
1.254 + /// \brief \ref named-templ-param "Named parameter" for setting
1.255 + /// Elevator type
1.256 + ///
1.257 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.258 + /// type. If this named parameter is used, then an external
1.259 + /// elevator object must be passed to the algorithm using the
1.260 + /// \ref elevator(Elevator&) "elevator()" function before calling
1.261 + /// \ref run() or \ref init().
1.262 + /// \sa SetStandardElevator
1.263 + template <typename _Elevator>
1.264 + struct SetElevator
1.265 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.266 + SetElevatorTraits<_Elevator> > {
1.267 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.268 + SetElevatorTraits<_Elevator> > Create;
1.269 + };
1.270 +
1.271 + template <typename _Elevator>
1.272 + struct SetStandardElevatorTraits : public Traits {
1.273 + typedef _Elevator Elevator;
1.274 + static Elevator *createElevator(const Digraph& digraph, int max_level) {
1.275 + return new Elevator(digraph, max_level);
1.276 + }
1.277 + };
1.278 +
1.279 + /// \brief \ref named-templ-param "Named parameter" for setting
1.280 + /// Elevator type with automatic allocation
1.281 + ///
1.282 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.283 + /// type with automatic allocation.
1.284 + /// The Elevator should have standard constructor interface to be
1.285 + /// able to automatically created by the algorithm (i.e. the
1.286 + /// digraph and the maximum level should be passed to it).
1.287 + /// However an external elevator object could also be passed to the
1.288 + /// algorithm with the \ref elevator(Elevator&) "elevator()" function
1.289 + /// before calling \ref run() or \ref init().
1.290 + /// \sa SetElevator
1.291 + template <typename _Elevator>
1.292 + struct SetStandardElevator
1.293 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.294 + SetStandardElevatorTraits<_Elevator> > {
1.295 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.296 + SetStandardElevatorTraits<_Elevator> > Create;
1.297 + };
1.298 +
1.299 + /// @}
1.300 +
1.301 + protected:
1.302 +
1.303 + Circulation() {}
1.304 +
1.305 + public:
1.306 +
1.307 + /// The constructor of the class.
1.308 +
1.309 + /// The constructor of the class.
1.310 + /// \param g The digraph the algorithm runs on.
1.311 + /// \param lo The lower bound capacity of the arcs.
1.312 + /// \param up The upper bound capacity of the arcs.
1.313 + /// \param delta The lower bound for the supply of the nodes.
1.314 + Circulation(const Digraph &g,const LCapMap &lo,
1.315 + const UCapMap &up,const DeltaMap &delta)
1.316 + : _g(g), _node_num(),
1.317 + _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
1.318 + _level(0), _local_level(false), _excess(0), _el() {}
1.319 +
1.320 + /// Destructor.
1.321 + ~Circulation() {
1.322 + destroyStructures();
1.323 + }
1.324 +
1.325 +
1.326 + private:
1.327 +
1.328 + void createStructures() {
1.329 + _node_num = _el = countNodes(_g);
1.330 +
1.331 + if (!_flow) {
1.332 + _flow = Traits::createFlowMap(_g);
1.333 + _local_flow = true;
1.334 + }
1.335 + if (!_level) {
1.336 + _level = Traits::createElevator(_g, _node_num);
1.337 + _local_level = true;
1.338 + }
1.339 + if (!_excess) {
1.340 + _excess = new ExcessMap(_g);
1.341 + }
1.342 + }
1.343 +
1.344 + void destroyStructures() {
1.345 + if (_local_flow) {
1.346 + delete _flow;
1.347 + }
1.348 + if (_local_level) {
1.349 + delete _level;
1.350 + }
1.351 + if (_excess) {
1.352 + delete _excess;
1.353 + }
1.354 + }
1.355 +
1.356 + public:
1.357 +
1.358 + /// Sets the lower bound capacity map.
1.359 +
1.360 + /// Sets the lower bound capacity map.
1.361 + /// \return <tt>(*this)</tt>
1.362 + Circulation& lowerCapMap(const LCapMap& map) {
1.363 + _lo = ↦
1.364 + return *this;
1.365 + }
1.366 +
1.367 + /// Sets the upper bound capacity map.
1.368 +
1.369 + /// Sets the upper bound capacity map.
1.370 + /// \return <tt>(*this)</tt>
1.371 + Circulation& upperCapMap(const LCapMap& map) {
1.372 + _up = ↦
1.373 + return *this;
1.374 + }
1.375 +
1.376 + /// Sets the lower bound map for the supply of the nodes.
1.377 +
1.378 + /// Sets the lower bound map for the supply of the nodes.
1.379 + /// \return <tt>(*this)</tt>
1.380 + Circulation& deltaMap(const DeltaMap& map) {
1.381 + _delta = ↦
1.382 + return *this;
1.383 + }
1.384 +
1.385 + /// \brief Sets the flow map.
1.386 + ///
1.387 + /// Sets the flow map.
1.388 + /// If you don't use this function before calling \ref run() or
1.389 + /// \ref init(), an instance will be allocated automatically.
1.390 + /// The destructor deallocates this automatically allocated map,
1.391 + /// of course.
1.392 + /// \return <tt>(*this)</tt>
1.393 + Circulation& flowMap(FlowMap& map) {
1.394 + if (_local_flow) {
1.395 + delete _flow;
1.396 + _local_flow = false;
1.397 + }
1.398 + _flow = ↦
1.399 + return *this;
1.400 + }
1.401 +
1.402 + /// \brief Sets the elevator used by algorithm.
1.403 + ///
1.404 + /// Sets the elevator used by algorithm.
1.405 + /// If you don't use this function before calling \ref run() or
1.406 + /// \ref init(), an instance will be allocated automatically.
1.407 + /// The destructor deallocates this automatically allocated elevator,
1.408 + /// of course.
1.409 + /// \return <tt>(*this)</tt>
1.410 + Circulation& elevator(Elevator& elevator) {
1.411 + if (_local_level) {
1.412 + delete _level;
1.413 + _local_level = false;
1.414 + }
1.415 + _level = &elevator;
1.416 + return *this;
1.417 + }
1.418 +
1.419 + /// \brief Returns a const reference to the elevator.
1.420 + ///
1.421 + /// Returns a const reference to the elevator.
1.422 + ///
1.423 + /// \pre Either \ref run() or \ref init() must be called before
1.424 + /// using this function.
1.425 + const Elevator& elevator() {
1.426 + return *_level;
1.427 + }
1.428 +
1.429 + /// \brief Sets the tolerance used by algorithm.
1.430 + ///
1.431 + /// Sets the tolerance used by algorithm.
1.432 + Circulation& tolerance(const Tolerance& tolerance) const {
1.433 + _tol = tolerance;
1.434 + return *this;
1.435 + }
1.436 +
1.437 + /// \brief Returns a const reference to the tolerance.
1.438 + ///
1.439 + /// Returns a const reference to the tolerance.
1.440 + const Tolerance& tolerance() const {
1.441 + return tolerance;
1.442 + }
1.443 +
1.444 + /// \name Execution Control
1.445 + /// The simplest way to execute the algorithm is to call \ref run().\n
1.446 + /// If you need more control on the initial solution or the execution,
1.447 + /// first you have to call one of the \ref init() functions, then
1.448 + /// the \ref start() function.
1.449 +
1.450 + ///@{
1.451 +
1.452 + /// Initializes the internal data structures.
1.453 +
1.454 + /// Initializes the internal data structures and sets all flow values
1.455 + /// to the lower bound.
1.456 + void init()
1.457 + {
1.458 + createStructures();
1.459 +
1.460 + for(NodeIt n(_g);n!=INVALID;++n) {
1.461 + _excess->set(n, (*_delta)[n]);
1.462 + }
1.463 +
1.464 + for (ArcIt e(_g);e!=INVALID;++e) {
1.465 + _flow->set(e, (*_lo)[e]);
1.466 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
1.467 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
1.468 + }
1.469 +
1.470 + // global relabeling tested, but in general case it provides
1.471 + // worse performance for random digraphs
1.472 + _level->initStart();
1.473 + for(NodeIt n(_g);n!=INVALID;++n)
1.474 + _level->initAddItem(n);
1.475 + _level->initFinish();
1.476 + for(NodeIt n(_g);n!=INVALID;++n)
1.477 + if(_tol.positive((*_excess)[n]))
1.478 + _level->activate(n);
1.479 + }
1.480 +
1.481 + /// Initializes the internal data structures using a greedy approach.
1.482 +
1.483 + /// Initializes the internal data structures using a greedy approach
1.484 + /// to construct the initial solution.
1.485 + void greedyInit()
1.486 + {
1.487 + createStructures();
1.488 +
1.489 + for(NodeIt n(_g);n!=INVALID;++n) {
1.490 + _excess->set(n, (*_delta)[n]);
1.491 + }
1.492 +
1.493 + for (ArcIt e(_g);e!=INVALID;++e) {
1.494 + if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
1.495 + _flow->set(e, (*_up)[e]);
1.496 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
1.497 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
1.498 + } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
1.499 + _flow->set(e, (*_lo)[e]);
1.500 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
1.501 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
1.502 + } else {
1.503 + Value fc = -(*_excess)[_g.target(e)];
1.504 + _flow->set(e, fc);
1.505 + _excess->set(_g.target(e), 0);
1.506 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
1.507 + }
1.508 + }
1.509 +
1.510 + _level->initStart();
1.511 + for(NodeIt n(_g);n!=INVALID;++n)
1.512 + _level->initAddItem(n);
1.513 + _level->initFinish();
1.514 + for(NodeIt n(_g);n!=INVALID;++n)
1.515 + if(_tol.positive((*_excess)[n]))
1.516 + _level->activate(n);
1.517 + }
1.518 +
1.519 + ///Executes the algorithm
1.520 +
1.521 + ///This function executes the algorithm.
1.522 + ///
1.523 + ///\return \c true if a feasible circulation is found.
1.524 + ///
1.525 + ///\sa barrier()
1.526 + ///\sa barrierMap()
1.527 + bool start()
1.528 + {
1.529 +
1.530 + Node act;
1.531 + Node bact=INVALID;
1.532 + Node last_activated=INVALID;
1.533 + while((act=_level->highestActive())!=INVALID) {
1.534 + int actlevel=(*_level)[act];
1.535 + int mlevel=_node_num;
1.536 + Value exc=(*_excess)[act];
1.537 +
1.538 + for(OutArcIt e(_g,act);e!=INVALID; ++e) {
1.539 + Node v = _g.target(e);
1.540 + Value fc=(*_up)[e]-(*_flow)[e];
1.541 + if(!_tol.positive(fc)) continue;
1.542 + if((*_level)[v]<actlevel) {
1.543 + if(!_tol.less(fc, exc)) {
1.544 + _flow->set(e, (*_flow)[e] + exc);
1.545 + _excess->set(v, (*_excess)[v] + exc);
1.546 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.547 + _level->activate(v);
1.548 + _excess->set(act,0);
1.549 + _level->deactivate(act);
1.550 + goto next_l;
1.551 + }
1.552 + else {
1.553 + _flow->set(e, (*_up)[e]);
1.554 + _excess->set(v, (*_excess)[v] + fc);
1.555 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.556 + _level->activate(v);
1.557 + exc-=fc;
1.558 + }
1.559 + }
1.560 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
1.561 + }
1.562 + for(InArcIt e(_g,act);e!=INVALID; ++e) {
1.563 + Node v = _g.source(e);
1.564 + Value fc=(*_flow)[e]-(*_lo)[e];
1.565 + if(!_tol.positive(fc)) continue;
1.566 + if((*_level)[v]<actlevel) {
1.567 + if(!_tol.less(fc, exc)) {
1.568 + _flow->set(e, (*_flow)[e] - exc);
1.569 + _excess->set(v, (*_excess)[v] + exc);
1.570 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.571 + _level->activate(v);
1.572 + _excess->set(act,0);
1.573 + _level->deactivate(act);
1.574 + goto next_l;
1.575 + }
1.576 + else {
1.577 + _flow->set(e, (*_lo)[e]);
1.578 + _excess->set(v, (*_excess)[v] + fc);
1.579 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.580 + _level->activate(v);
1.581 + exc-=fc;
1.582 + }
1.583 + }
1.584 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
1.585 + }
1.586 +
1.587 + _excess->set(act, exc);
1.588 + if(!_tol.positive(exc)) _level->deactivate(act);
1.589 + else if(mlevel==_node_num) {
1.590 + _level->liftHighestActiveToTop();
1.591 + _el = _node_num;
1.592 + return false;
1.593 + }
1.594 + else {
1.595 + _level->liftHighestActive(mlevel+1);
1.596 + if(_level->onLevel(actlevel)==0) {
1.597 + _el = actlevel;
1.598 + return false;
1.599 + }
1.600 + }
1.601 + next_l:
1.602 + ;
1.603 + }
1.604 + return true;
1.605 + }
1.606 +
1.607 + /// Runs the algorithm.
1.608 +
1.609 + /// This function runs the algorithm.
1.610 + ///
1.611 + /// \return \c true if a feasible circulation is found.
1.612 + ///
1.613 + /// \note Apart from the return value, c.run() is just a shortcut of
1.614 + /// the following code.
1.615 + /// \code
1.616 + /// c.greedyInit();
1.617 + /// c.start();
1.618 + /// \endcode
1.619 + bool run() {
1.620 + greedyInit();
1.621 + return start();
1.622 + }
1.623 +
1.624 + /// @}
1.625 +
1.626 + /// \name Query Functions
1.627 + /// The results of the circulation algorithm can be obtained using
1.628 + /// these functions.\n
1.629 + /// Either \ref run() or \ref start() should be called before
1.630 + /// using them.
1.631 +
1.632 + ///@{
1.633 +
1.634 + /// \brief Returns the flow on the given arc.
1.635 + ///
1.636 + /// Returns the flow on the given arc.
1.637 + ///
1.638 + /// \pre Either \ref run() or \ref init() must be called before
1.639 + /// using this function.
1.640 + Value flow(const Arc& arc) const {
1.641 + return (*_flow)[arc];
1.642 + }
1.643 +
1.644 + /// \brief Returns a const reference to the flow map.
1.645 + ///
1.646 + /// Returns a const reference to the arc map storing the found flow.
1.647 + ///
1.648 + /// \pre Either \ref run() or \ref init() must be called before
1.649 + /// using this function.
1.650 + const FlowMap& flowMap() {
1.651 + return *_flow;
1.652 + }
1.653 +
1.654 + /**
1.655 + \brief Returns \c true if the given node is in a barrier.
1.656 +
1.657 + Barrier is a set \e B of nodes for which
1.658 +
1.659 + \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
1.660 + \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
1.661 +
1.662 + holds. The existence of a set with this property prooves that a
1.663 + feasible circualtion cannot exist.
1.664 +
1.665 + This function returns \c true if the given node is in the found
1.666 + barrier. If a feasible circulation is found, the function
1.667 + gives back \c false for every node.
1.668 +
1.669 + \pre Either \ref run() or \ref init() must be called before
1.670 + using this function.
1.671 +
1.672 + \sa barrierMap()
1.673 + \sa checkBarrier()
1.674 + */
1.675 + bool barrier(const Node& node)
1.676 + {
1.677 + return (*_level)[node] >= _el;
1.678 + }
1.679 +
1.680 + /// \brief Gives back a barrier.
1.681 + ///
1.682 + /// This function sets \c bar to the characteristic vector of the
1.683 + /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
1.684 + /// node map with \c bool (or convertible) value type.
1.685 + ///
1.686 + /// If a feasible circulation is found, the function gives back an
1.687 + /// empty set, so \c bar[v] will be \c false for all nodes \c v.
1.688 + ///
1.689 + /// \note This function calls \ref barrier() for each node,
1.690 + /// so it runs in \f$O(n)\f$ time.
1.691 + ///
1.692 + /// \pre Either \ref run() or \ref init() must be called before
1.693 + /// using this function.
1.694 + ///
1.695 + /// \sa barrier()
1.696 + /// \sa checkBarrier()
1.697 + template<class BarrierMap>
1.698 + void barrierMap(BarrierMap &bar)
1.699 + {
1.700 + for(NodeIt n(_g);n!=INVALID;++n)
1.701 + bar.set(n, (*_level)[n] >= _el);
1.702 + }
1.703 +
1.704 + /// @}
1.705 +
1.706 + /// \name Checker Functions
1.707 + /// The feasibility of the results can be checked using
1.708 + /// these functions.\n
1.709 + /// Either \ref run() or \ref start() should be called before
1.710 + /// using them.
1.711 +
1.712 + ///@{
1.713 +
1.714 + ///Check if the found flow is a feasible circulation
1.715 +
1.716 + ///Check if the found flow is a feasible circulation,
1.717 + ///
1.718 + bool checkFlow() {
1.719 + for(ArcIt e(_g);e!=INVALID;++e)
1.720 + if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
1.721 + for(NodeIt n(_g);n!=INVALID;++n)
1.722 + {
1.723 + Value dif=-(*_delta)[n];
1.724 + for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
1.725 + for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
1.726 + if(_tol.negative(dif)) return false;
1.727 + }
1.728 + return true;
1.729 + }
1.730 +
1.731 + ///Check whether or not the last execution provides a barrier
1.732 +
1.733 + ///Check whether or not the last execution provides a barrier.
1.734 + ///\sa barrier()
1.735 + ///\sa barrierMap()
1.736 + bool checkBarrier()
1.737 + {
1.738 + Value delta=0;
1.739 + for(NodeIt n(_g);n!=INVALID;++n)
1.740 + if(barrier(n))
1.741 + delta-=(*_delta)[n];
1.742 + for(ArcIt e(_g);e!=INVALID;++e)
1.743 + {
1.744 + Node s=_g.source(e);
1.745 + Node t=_g.target(e);
1.746 + if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
1.747 + else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
1.748 + }
1.749 + return _tol.negative(delta);
1.750 + }
1.751 +
1.752 + /// @}
1.753 +
1.754 + };
1.755 +
1.756 +}
1.757 +
1.758 +#endif