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