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