1.1 --- a/lemon/Makefile.am Mon Dec 01 13:49:55 2008 +0000
1.2 +++ b/lemon/Makefile.am Mon Dec 01 14:18:40 2008 +0000
1.3 @@ -20,6 +20,7 @@
1.4 lemon/assert.h \
1.5 lemon/bfs.h \
1.6 lemon/bin_heap.h \
1.7 + lemon/circulation.h \
1.8 lemon/color.h \
1.9 lemon/concept_check.h \
1.10 lemon/counter.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/circulation.h Mon Dec 01 14:18:40 2008 +0000
2.3 @@ -0,0 +1,755 @@
2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library.
2.7 + *
2.8 + * Copyright (C) 2003-2008
2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.11 + *
2.12 + * Permission to use, modify and distribute this software is granted
2.13 + * provided that this copyright notice appears in all copies. For
2.14 + * precise terms see the accompanying LICENSE file.
2.15 + *
2.16 + * This software is provided "AS IS" with no warranty of any kind,
2.17 + * express or implied, and with no claim as to its suitability for any
2.18 + * purpose.
2.19 + *
2.20 + */
2.21 +
2.22 +#ifndef LEMON_CIRCULATION_H
2.23 +#define LEMON_CIRCULATION_H
2.24 +
2.25 +#include <lemon/tolerance.h>
2.26 +#include <lemon/elevator.h>
2.27 +
2.28 +///\ingroup max_flow
2.29 +///\file
2.30 +///\brief Push-relabel algorithm for finding a feasible circulation.
2.31 +///
2.32 +namespace lemon {
2.33 +
2.34 + /// \brief Default traits class of Circulation class.
2.35 + ///
2.36 + /// Default traits class of Circulation class.
2.37 + /// \tparam _Diraph Digraph type.
2.38 + /// \tparam _LCapMap Lower bound capacity map type.
2.39 + /// \tparam _UCapMap Upper bound capacity map type.
2.40 + /// \tparam _DeltaMap Delta map type.
2.41 + template <typename _Diraph, typename _LCapMap,
2.42 + typename _UCapMap, typename _DeltaMap>
2.43 + struct CirculationDefaultTraits {
2.44 +
2.45 + /// \brief The type of the digraph the algorithm runs on.
2.46 + typedef _Diraph Digraph;
2.47 +
2.48 + /// \brief The type of the map that stores the circulation lower
2.49 + /// bound.
2.50 + ///
2.51 + /// The type of the map that stores the circulation lower bound.
2.52 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
2.53 + typedef _LCapMap LCapMap;
2.54 +
2.55 + /// \brief The type of the map that stores the circulation upper
2.56 + /// bound.
2.57 + ///
2.58 + /// The type of the map that stores the circulation upper bound.
2.59 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
2.60 + typedef _UCapMap UCapMap;
2.61 +
2.62 + /// \brief The type of the map that stores the lower bound for
2.63 + /// the supply of the nodes.
2.64 + ///
2.65 + /// The type of the map that stores the lower bound for the supply
2.66 + /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
2.67 + /// concept.
2.68 + typedef _DeltaMap DeltaMap;
2.69 +
2.70 + /// \brief The type of the flow values.
2.71 + typedef typename DeltaMap::Value Value;
2.72 +
2.73 + /// \brief The type of the map that stores the flow values.
2.74 + ///
2.75 + /// The type of the map that stores the flow values.
2.76 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
2.77 + typedef typename Digraph::template ArcMap<Value> FlowMap;
2.78 +
2.79 + /// \brief Instantiates a FlowMap.
2.80 + ///
2.81 + /// This function instantiates a \ref FlowMap.
2.82 + /// \param digraph The digraph, to which we would like to define
2.83 + /// the flow map.
2.84 + static FlowMap* createFlowMap(const Digraph& digraph) {
2.85 + return new FlowMap(digraph);
2.86 + }
2.87 +
2.88 + /// \brief The elevator type used by the algorithm.
2.89 + ///
2.90 + /// The elevator type used by the algorithm.
2.91 + ///
2.92 + /// \sa Elevator
2.93 + /// \sa LinkedElevator
2.94 + typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
2.95 +
2.96 + /// \brief Instantiates an Elevator.
2.97 + ///
2.98 + /// This function instantiates an \ref Elevator.
2.99 + /// \param digraph The digraph, to which we would like to define
2.100 + /// the elevator.
2.101 + /// \param max_level The maximum level of the elevator.
2.102 + static Elevator* createElevator(const Digraph& digraph, int max_level) {
2.103 + return new Elevator(digraph, max_level);
2.104 + }
2.105 +
2.106 + /// \brief The tolerance used by the algorithm
2.107 + ///
2.108 + /// The tolerance used by the algorithm to handle inexact computation.
2.109 + typedef lemon::Tolerance<Value> Tolerance;
2.110 +
2.111 + };
2.112 +
2.113 + /**
2.114 + \brief Push-relabel algorithm for the network circulation problem.
2.115 +
2.116 + \ingroup max_flow
2.117 + This class implements a push-relabel algorithm for the network
2.118 + circulation problem.
2.119 + It is to find a feasible circulation when lower and upper bounds
2.120 + are given for the flow values on the arcs and lower bounds
2.121 + are given for the supply values of the nodes.
2.122 +
2.123 + The exact formulation of this problem is the following.
2.124 + Let \f$G=(V,A)\f$ be a digraph,
2.125 + \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
2.126 + \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
2.127 + \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
2.128 + \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
2.129 + \geq delta(v) \quad \forall v\in V, \f]
2.130 + \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
2.131 + \note \f$delta(v)\f$ specifies a lower bound for the supply of node
2.132 + \f$v\f$. It can be either positive or negative, however note that
2.133 + \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
2.134 + have a feasible solution.
2.135 +
2.136 + \note A special case of this problem is when
2.137 + \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
2.138 + will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
2.139 + Thus a feasible solution for the
2.140 + \ref min_cost_flow "minimum cost flow" problem can be calculated
2.141 + in this way.
2.142 +
2.143 + \tparam _Digraph The type of the digraph the algorithm runs on.
2.144 + \tparam _LCapMap The type of the lower bound capacity map. The default
2.145 + map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
2.146 + \tparam _UCapMap The type of the upper bound capacity map. The default
2.147 + map type is \c _LCapMap.
2.148 + \tparam _DeltaMap The type of the map that stores the lower bound
2.149 + for the supply of the nodes. The default map type is
2.150 + \c _Digraph::ArcMap<_UCapMap::Value>.
2.151 + */
2.152 +#ifdef DOXYGEN
2.153 +template< typename _Digraph,
2.154 + typename _LCapMap,
2.155 + typename _UCapMap,
2.156 + typename _DeltaMap,
2.157 + typename _Traits >
2.158 +#else
2.159 +template< typename _Digraph,
2.160 + typename _LCapMap = typename _Digraph::template ArcMap<int>,
2.161 + typename _UCapMap = _LCapMap,
2.162 + typename _DeltaMap = typename _Digraph::
2.163 + template NodeMap<typename _UCapMap::Value>,
2.164 + typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
2.165 + _UCapMap, _DeltaMap> >
2.166 +#endif
2.167 + class Circulation {
2.168 + public:
2.169 +
2.170 + ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
2.171 + typedef _Traits Traits;
2.172 + ///The type of the digraph the algorithm runs on.
2.173 + typedef typename Traits::Digraph Digraph;
2.174 + ///The type of the flow values.
2.175 + typedef typename Traits::Value Value;
2.176 +
2.177 + /// The type of the lower bound capacity map.
2.178 + typedef typename Traits::LCapMap LCapMap;
2.179 + /// The type of the upper bound capacity map.
2.180 + typedef typename Traits::UCapMap UCapMap;
2.181 + /// \brief The type of the map that stores the lower bound for
2.182 + /// the supply of the nodes.
2.183 + typedef typename Traits::DeltaMap DeltaMap;
2.184 + ///The type of the flow map.
2.185 + typedef typename Traits::FlowMap FlowMap;
2.186 +
2.187 + ///The type of the elevator.
2.188 + typedef typename Traits::Elevator Elevator;
2.189 + ///The type of the tolerance.
2.190 + typedef typename Traits::Tolerance Tolerance;
2.191 +
2.192 + private:
2.193 +
2.194 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
2.195 +
2.196 + const Digraph &_g;
2.197 + int _node_num;
2.198 +
2.199 + const LCapMap *_lo;
2.200 + const UCapMap *_up;
2.201 + const DeltaMap *_delta;
2.202 +
2.203 + FlowMap *_flow;
2.204 + bool _local_flow;
2.205 +
2.206 + Elevator* _level;
2.207 + bool _local_level;
2.208 +
2.209 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
2.210 + ExcessMap* _excess;
2.211 +
2.212 + Tolerance _tol;
2.213 + int _el;
2.214 +
2.215 + public:
2.216 +
2.217 + typedef Circulation Create;
2.218 +
2.219 + ///\name Named Template Parameters
2.220 +
2.221 + ///@{
2.222 +
2.223 + template <typename _FlowMap>
2.224 + struct SetFlowMapTraits : public Traits {
2.225 + typedef _FlowMap FlowMap;
2.226 + static FlowMap *createFlowMap(const Digraph&) {
2.227 + LEMON_ASSERT(false, "FlowMap is not initialized");
2.228 + return 0; // ignore warnings
2.229 + }
2.230 + };
2.231 +
2.232 + /// \brief \ref named-templ-param "Named parameter" for setting
2.233 + /// FlowMap type
2.234 + ///
2.235 + /// \ref named-templ-param "Named parameter" for setting FlowMap
2.236 + /// type.
2.237 + template <typename _FlowMap>
2.238 + struct SetFlowMap
2.239 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
2.240 + SetFlowMapTraits<_FlowMap> > {
2.241 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
2.242 + SetFlowMapTraits<_FlowMap> > Create;
2.243 + };
2.244 +
2.245 + template <typename _Elevator>
2.246 + struct SetElevatorTraits : public Traits {
2.247 + typedef _Elevator Elevator;
2.248 + static Elevator *createElevator(const Digraph&, int) {
2.249 + LEMON_ASSERT(false, "Elevator is not initialized");
2.250 + return 0; // ignore warnings
2.251 + }
2.252 + };
2.253 +
2.254 + /// \brief \ref named-templ-param "Named parameter" for setting
2.255 + /// Elevator type
2.256 + ///
2.257 + /// \ref named-templ-param "Named parameter" for setting Elevator
2.258 + /// type. If this named parameter is used, then an external
2.259 + /// elevator object must be passed to the algorithm using the
2.260 + /// \ref elevator(Elevator&) "elevator()" function before calling
2.261 + /// \ref run() or \ref init().
2.262 + /// \sa SetStandardElevator
2.263 + template <typename _Elevator>
2.264 + struct SetElevator
2.265 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
2.266 + SetElevatorTraits<_Elevator> > {
2.267 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
2.268 + SetElevatorTraits<_Elevator> > Create;
2.269 + };
2.270 +
2.271 + template <typename _Elevator>
2.272 + struct SetStandardElevatorTraits : public Traits {
2.273 + typedef _Elevator Elevator;
2.274 + static Elevator *createElevator(const Digraph& digraph, int max_level) {
2.275 + return new Elevator(digraph, max_level);
2.276 + }
2.277 + };
2.278 +
2.279 + /// \brief \ref named-templ-param "Named parameter" for setting
2.280 + /// Elevator type with automatic allocation
2.281 + ///
2.282 + /// \ref named-templ-param "Named parameter" for setting Elevator
2.283 + /// type with automatic allocation.
2.284 + /// The Elevator should have standard constructor interface to be
2.285 + /// able to automatically created by the algorithm (i.e. the
2.286 + /// digraph and the maximum level should be passed to it).
2.287 + /// However an external elevator object could also be passed to the
2.288 + /// algorithm with the \ref elevator(Elevator&) "elevator()" function
2.289 + /// before calling \ref run() or \ref init().
2.290 + /// \sa SetElevator
2.291 + template <typename _Elevator>
2.292 + struct SetStandardElevator
2.293 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
2.294 + SetStandardElevatorTraits<_Elevator> > {
2.295 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
2.296 + SetStandardElevatorTraits<_Elevator> > Create;
2.297 + };
2.298 +
2.299 + /// @}
2.300 +
2.301 + protected:
2.302 +
2.303 + Circulation() {}
2.304 +
2.305 + public:
2.306 +
2.307 + /// The constructor of the class.
2.308 +
2.309 + /// The constructor of the class.
2.310 + /// \param g The digraph the algorithm runs on.
2.311 + /// \param lo The lower bound capacity of the arcs.
2.312 + /// \param up The upper bound capacity of the arcs.
2.313 + /// \param delta The lower bound for the supply of the nodes.
2.314 + Circulation(const Digraph &g,const LCapMap &lo,
2.315 + const UCapMap &up,const DeltaMap &delta)
2.316 + : _g(g), _node_num(),
2.317 + _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
2.318 + _level(0), _local_level(false), _excess(0), _el() {}
2.319 +
2.320 + /// Destructor.
2.321 + ~Circulation() {
2.322 + destroyStructures();
2.323 + }
2.324 +
2.325 +
2.326 + private:
2.327 +
2.328 + void createStructures() {
2.329 + _node_num = _el = countNodes(_g);
2.330 +
2.331 + if (!_flow) {
2.332 + _flow = Traits::createFlowMap(_g);
2.333 + _local_flow = true;
2.334 + }
2.335 + if (!_level) {
2.336 + _level = Traits::createElevator(_g, _node_num);
2.337 + _local_level = true;
2.338 + }
2.339 + if (!_excess) {
2.340 + _excess = new ExcessMap(_g);
2.341 + }
2.342 + }
2.343 +
2.344 + void destroyStructures() {
2.345 + if (_local_flow) {
2.346 + delete _flow;
2.347 + }
2.348 + if (_local_level) {
2.349 + delete _level;
2.350 + }
2.351 + if (_excess) {
2.352 + delete _excess;
2.353 + }
2.354 + }
2.355 +
2.356 + public:
2.357 +
2.358 + /// Sets the lower bound capacity map.
2.359 +
2.360 + /// Sets the lower bound capacity map.
2.361 + /// \return <tt>(*this)</tt>
2.362 + Circulation& lowerCapMap(const LCapMap& map) {
2.363 + _lo = ↦
2.364 + return *this;
2.365 + }
2.366 +
2.367 + /// Sets the upper bound capacity map.
2.368 +
2.369 + /// Sets the upper bound capacity map.
2.370 + /// \return <tt>(*this)</tt>
2.371 + Circulation& upperCapMap(const LCapMap& map) {
2.372 + _up = ↦
2.373 + return *this;
2.374 + }
2.375 +
2.376 + /// Sets the lower bound map for the supply of the nodes.
2.377 +
2.378 + /// Sets the lower bound map for the supply of the nodes.
2.379 + /// \return <tt>(*this)</tt>
2.380 + Circulation& deltaMap(const DeltaMap& map) {
2.381 + _delta = ↦
2.382 + return *this;
2.383 + }
2.384 +
2.385 + /// \brief Sets the flow map.
2.386 + ///
2.387 + /// Sets the flow map.
2.388 + /// If you don't use this function before calling \ref run() or
2.389 + /// \ref init(), an instance will be allocated automatically.
2.390 + /// The destructor deallocates this automatically allocated map,
2.391 + /// of course.
2.392 + /// \return <tt>(*this)</tt>
2.393 + Circulation& flowMap(FlowMap& map) {
2.394 + if (_local_flow) {
2.395 + delete _flow;
2.396 + _local_flow = false;
2.397 + }
2.398 + _flow = ↦
2.399 + return *this;
2.400 + }
2.401 +
2.402 + /// \brief Sets the elevator used by algorithm.
2.403 + ///
2.404 + /// Sets the elevator used by algorithm.
2.405 + /// If you don't use this function before calling \ref run() or
2.406 + /// \ref init(), an instance will be allocated automatically.
2.407 + /// The destructor deallocates this automatically allocated elevator,
2.408 + /// of course.
2.409 + /// \return <tt>(*this)</tt>
2.410 + Circulation& elevator(Elevator& elevator) {
2.411 + if (_local_level) {
2.412 + delete _level;
2.413 + _local_level = false;
2.414 + }
2.415 + _level = &elevator;
2.416 + return *this;
2.417 + }
2.418 +
2.419 + /// \brief Returns a const reference to the elevator.
2.420 + ///
2.421 + /// Returns a const reference to the elevator.
2.422 + ///
2.423 + /// \pre Either \ref run() or \ref init() must be called before
2.424 + /// using this function.
2.425 + const Elevator& elevator() {
2.426 + return *_level;
2.427 + }
2.428 +
2.429 + /// \brief Sets the tolerance used by algorithm.
2.430 + ///
2.431 + /// Sets the tolerance used by algorithm.
2.432 + Circulation& tolerance(const Tolerance& tolerance) const {
2.433 + _tol = tolerance;
2.434 + return *this;
2.435 + }
2.436 +
2.437 + /// \brief Returns a const reference to the tolerance.
2.438 + ///
2.439 + /// Returns a const reference to the tolerance.
2.440 + const Tolerance& tolerance() const {
2.441 + return tolerance;
2.442 + }
2.443 +
2.444 + /// \name Execution Control
2.445 + /// The simplest way to execute the algorithm is to call \ref run().\n
2.446 + /// If you need more control on the initial solution or the execution,
2.447 + /// first you have to call one of the \ref init() functions, then
2.448 + /// the \ref start() function.
2.449 +
2.450 + ///@{
2.451 +
2.452 + /// Initializes the internal data structures.
2.453 +
2.454 + /// Initializes the internal data structures and sets all flow values
2.455 + /// to the lower bound.
2.456 + void init()
2.457 + {
2.458 + createStructures();
2.459 +
2.460 + for(NodeIt n(_g);n!=INVALID;++n) {
2.461 + _excess->set(n, (*_delta)[n]);
2.462 + }
2.463 +
2.464 + for (ArcIt e(_g);e!=INVALID;++e) {
2.465 + _flow->set(e, (*_lo)[e]);
2.466 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
2.467 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
2.468 + }
2.469 +
2.470 + // global relabeling tested, but in general case it provides
2.471 + // worse performance for random digraphs
2.472 + _level->initStart();
2.473 + for(NodeIt n(_g);n!=INVALID;++n)
2.474 + _level->initAddItem(n);
2.475 + _level->initFinish();
2.476 + for(NodeIt n(_g);n!=INVALID;++n)
2.477 + if(_tol.positive((*_excess)[n]))
2.478 + _level->activate(n);
2.479 + }
2.480 +
2.481 + /// Initializes the internal data structures using a greedy approach.
2.482 +
2.483 + /// Initializes the internal data structures using a greedy approach
2.484 + /// to construct the initial solution.
2.485 + void greedyInit()
2.486 + {
2.487 + createStructures();
2.488 +
2.489 + for(NodeIt n(_g);n!=INVALID;++n) {
2.490 + _excess->set(n, (*_delta)[n]);
2.491 + }
2.492 +
2.493 + for (ArcIt e(_g);e!=INVALID;++e) {
2.494 + if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
2.495 + _flow->set(e, (*_up)[e]);
2.496 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
2.497 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
2.498 + } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
2.499 + _flow->set(e, (*_lo)[e]);
2.500 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
2.501 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
2.502 + } else {
2.503 + Value fc = -(*_excess)[_g.target(e)];
2.504 + _flow->set(e, fc);
2.505 + _excess->set(_g.target(e), 0);
2.506 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
2.507 + }
2.508 + }
2.509 +
2.510 + _level->initStart();
2.511 + for(NodeIt n(_g);n!=INVALID;++n)
2.512 + _level->initAddItem(n);
2.513 + _level->initFinish();
2.514 + for(NodeIt n(_g);n!=INVALID;++n)
2.515 + if(_tol.positive((*_excess)[n]))
2.516 + _level->activate(n);
2.517 + }
2.518 +
2.519 + ///Executes the algorithm
2.520 +
2.521 + ///This function executes the algorithm.
2.522 + ///
2.523 + ///\return \c true if a feasible circulation is found.
2.524 + ///
2.525 + ///\sa barrier()
2.526 + ///\sa barrierMap()
2.527 + bool start()
2.528 + {
2.529 +
2.530 + Node act;
2.531 + Node bact=INVALID;
2.532 + Node last_activated=INVALID;
2.533 + while((act=_level->highestActive())!=INVALID) {
2.534 + int actlevel=(*_level)[act];
2.535 + int mlevel=_node_num;
2.536 + Value exc=(*_excess)[act];
2.537 +
2.538 + for(OutArcIt e(_g,act);e!=INVALID; ++e) {
2.539 + Node v = _g.target(e);
2.540 + Value fc=(*_up)[e]-(*_flow)[e];
2.541 + if(!_tol.positive(fc)) continue;
2.542 + if((*_level)[v]<actlevel) {
2.543 + if(!_tol.less(fc, exc)) {
2.544 + _flow->set(e, (*_flow)[e] + exc);
2.545 + _excess->set(v, (*_excess)[v] + exc);
2.546 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.547 + _level->activate(v);
2.548 + _excess->set(act,0);
2.549 + _level->deactivate(act);
2.550 + goto next_l;
2.551 + }
2.552 + else {
2.553 + _flow->set(e, (*_up)[e]);
2.554 + _excess->set(v, (*_excess)[v] + fc);
2.555 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.556 + _level->activate(v);
2.557 + exc-=fc;
2.558 + }
2.559 + }
2.560 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
2.561 + }
2.562 + for(InArcIt e(_g,act);e!=INVALID; ++e) {
2.563 + Node v = _g.source(e);
2.564 + Value fc=(*_flow)[e]-(*_lo)[e];
2.565 + if(!_tol.positive(fc)) continue;
2.566 + if((*_level)[v]<actlevel) {
2.567 + if(!_tol.less(fc, exc)) {
2.568 + _flow->set(e, (*_flow)[e] - exc);
2.569 + _excess->set(v, (*_excess)[v] + exc);
2.570 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.571 + _level->activate(v);
2.572 + _excess->set(act,0);
2.573 + _level->deactivate(act);
2.574 + goto next_l;
2.575 + }
2.576 + else {
2.577 + _flow->set(e, (*_lo)[e]);
2.578 + _excess->set(v, (*_excess)[v] + fc);
2.579 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.580 + _level->activate(v);
2.581 + exc-=fc;
2.582 + }
2.583 + }
2.584 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
2.585 + }
2.586 +
2.587 + _excess->set(act, exc);
2.588 + if(!_tol.positive(exc)) _level->deactivate(act);
2.589 + else if(mlevel==_node_num) {
2.590 + _level->liftHighestActiveToTop();
2.591 + _el = _node_num;
2.592 + return false;
2.593 + }
2.594 + else {
2.595 + _level->liftHighestActive(mlevel+1);
2.596 + if(_level->onLevel(actlevel)==0) {
2.597 + _el = actlevel;
2.598 + return false;
2.599 + }
2.600 + }
2.601 + next_l:
2.602 + ;
2.603 + }
2.604 + return true;
2.605 + }
2.606 +
2.607 + /// Runs the algorithm.
2.608 +
2.609 + /// This function runs the algorithm.
2.610 + ///
2.611 + /// \return \c true if a feasible circulation is found.
2.612 + ///
2.613 + /// \note Apart from the return value, c.run() is just a shortcut of
2.614 + /// the following code.
2.615 + /// \code
2.616 + /// c.greedyInit();
2.617 + /// c.start();
2.618 + /// \endcode
2.619 + bool run() {
2.620 + greedyInit();
2.621 + return start();
2.622 + }
2.623 +
2.624 + /// @}
2.625 +
2.626 + /// \name Query Functions
2.627 + /// The results of the circulation algorithm can be obtained using
2.628 + /// these functions.\n
2.629 + /// Either \ref run() or \ref start() should be called before
2.630 + /// using them.
2.631 +
2.632 + ///@{
2.633 +
2.634 + /// \brief Returns the flow on the given arc.
2.635 + ///
2.636 + /// Returns the flow on the given arc.
2.637 + ///
2.638 + /// \pre Either \ref run() or \ref init() must be called before
2.639 + /// using this function.
2.640 + Value flow(const Arc& arc) const {
2.641 + return (*_flow)[arc];
2.642 + }
2.643 +
2.644 + /// \brief Returns a const reference to the flow map.
2.645 + ///
2.646 + /// Returns a const reference to the arc map storing the found flow.
2.647 + ///
2.648 + /// \pre Either \ref run() or \ref init() must be called before
2.649 + /// using this function.
2.650 + const FlowMap& flowMap() {
2.651 + return *_flow;
2.652 + }
2.653 +
2.654 + /**
2.655 + \brief Returns \c true if the given node is in a barrier.
2.656 +
2.657 + Barrier is a set \e B of nodes for which
2.658 +
2.659 + \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
2.660 + \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
2.661 +
2.662 + holds. The existence of a set with this property prooves that a
2.663 + feasible circualtion cannot exist.
2.664 +
2.665 + This function returns \c true if the given node is in the found
2.666 + barrier. If a feasible circulation is found, the function
2.667 + gives back \c false for every node.
2.668 +
2.669 + \pre Either \ref run() or \ref init() must be called before
2.670 + using this function.
2.671 +
2.672 + \sa barrierMap()
2.673 + \sa checkBarrier()
2.674 + */
2.675 + bool barrier(const Node& node)
2.676 + {
2.677 + return (*_level)[node] >= _el;
2.678 + }
2.679 +
2.680 + /// \brief Gives back a barrier.
2.681 + ///
2.682 + /// This function sets \c bar to the characteristic vector of the
2.683 + /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
2.684 + /// node map with \c bool (or convertible) value type.
2.685 + ///
2.686 + /// If a feasible circulation is found, the function gives back an
2.687 + /// empty set, so \c bar[v] will be \c false for all nodes \c v.
2.688 + ///
2.689 + /// \note This function calls \ref barrier() for each node,
2.690 + /// so it runs in \f$O(n)\f$ time.
2.691 + ///
2.692 + /// \pre Either \ref run() or \ref init() must be called before
2.693 + /// using this function.
2.694 + ///
2.695 + /// \sa barrier()
2.696 + /// \sa checkBarrier()
2.697 + template<class BarrierMap>
2.698 + void barrierMap(BarrierMap &bar)
2.699 + {
2.700 + for(NodeIt n(_g);n!=INVALID;++n)
2.701 + bar.set(n, (*_level)[n] >= _el);
2.702 + }
2.703 +
2.704 + /// @}
2.705 +
2.706 + /// \name Checker Functions
2.707 + /// The feasibility of the results can be checked using
2.708 + /// these functions.\n
2.709 + /// Either \ref run() or \ref start() should be called before
2.710 + /// using them.
2.711 +
2.712 + ///@{
2.713 +
2.714 + ///Check if the found flow is a feasible circulation
2.715 +
2.716 + ///Check if the found flow is a feasible circulation,
2.717 + ///
2.718 + bool checkFlow() {
2.719 + for(ArcIt e(_g);e!=INVALID;++e)
2.720 + if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
2.721 + for(NodeIt n(_g);n!=INVALID;++n)
2.722 + {
2.723 + Value dif=-(*_delta)[n];
2.724 + for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
2.725 + for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
2.726 + if(_tol.negative(dif)) return false;
2.727 + }
2.728 + return true;
2.729 + }
2.730 +
2.731 + ///Check whether or not the last execution provides a barrier
2.732 +
2.733 + ///Check whether or not the last execution provides a barrier.
2.734 + ///\sa barrier()
2.735 + ///\sa barrierMap()
2.736 + bool checkBarrier()
2.737 + {
2.738 + Value delta=0;
2.739 + for(NodeIt n(_g);n!=INVALID;++n)
2.740 + if(barrier(n))
2.741 + delta-=(*_delta)[n];
2.742 + for(ArcIt e(_g);e!=INVALID;++e)
2.743 + {
2.744 + Node s=_g.source(e);
2.745 + Node t=_g.target(e);
2.746 + if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
2.747 + else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
2.748 + }
2.749 + return _tol.negative(delta);
2.750 + }
2.751 +
2.752 + /// @}
2.753 +
2.754 + };
2.755 +
2.756 +}
2.757 +
2.758 +#endif
3.1 --- a/test/Makefile.am Mon Dec 01 13:49:55 2008 +0000
3.2 +++ b/test/Makefile.am Mon Dec 01 14:18:40 2008 +0000
3.3 @@ -9,6 +9,7 @@
3.4
3.5 check_PROGRAMS += \
3.6 test/bfs_test \
3.7 + test/circulation_test \
3.8 test/counter_test \
3.9 test/dfs_test \
3.10 test/digraph_test \
3.11 @@ -35,6 +36,7 @@
3.12 XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
3.13
3.14 test_bfs_test_SOURCES = test/bfs_test.cc
3.15 +test_circulation_test_SOURCES = test/circulation_test.cc
3.16 test_counter_test_SOURCES = test/counter_test.cc
3.17 test_dfs_test_SOURCES = test/dfs_test.cc
3.18 test_digraph_test_SOURCES = test/digraph_test.cc
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/test/circulation_test.cc Mon Dec 01 14:18:40 2008 +0000
4.3 @@ -0,0 +1,156 @@
4.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
4.5 + *
4.6 + * This file is a part of LEMON, a generic C++ optimization library.
4.7 + *
4.8 + * Copyright (C) 2003-2008
4.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
4.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
4.11 + *
4.12 + * Permission to use, modify and distribute this software is granted
4.13 + * provided that this copyright notice appears in all copies. For
4.14 + * precise terms see the accompanying LICENSE file.
4.15 + *
4.16 + * This software is provided "AS IS" with no warranty of any kind,
4.17 + * express or implied, and with no claim as to its suitability for any
4.18 + * purpose.
4.19 + *
4.20 + */
4.21 +
4.22 +#include <iostream>
4.23 +
4.24 +#include "test_tools.h"
4.25 +#include <lemon/list_graph.h>
4.26 +#include <lemon/circulation.h>
4.27 +#include <lemon/lgf_reader.h>
4.28 +#include <lemon/concepts/digraph.h>
4.29 +#include <lemon/concepts/maps.h>
4.30 +
4.31 +using namespace lemon;
4.32 +
4.33 +char test_lgf[] =
4.34 + "@nodes\n"
4.35 + "label\n"
4.36 + "0\n"
4.37 + "1\n"
4.38 + "2\n"
4.39 + "3\n"
4.40 + "4\n"
4.41 + "5\n"
4.42 + "@arcs\n"
4.43 + " lcap ucap\n"
4.44 + "0 1 2 10\n"
4.45 + "0 2 2 6\n"
4.46 + "1 3 4 7\n"
4.47 + "1 4 0 5\n"
4.48 + "2 4 1 3\n"
4.49 + "3 5 3 8\n"
4.50 + "4 5 3 7\n"
4.51 + "@attributes\n"
4.52 + "source 0\n"
4.53 + "sink 5\n";
4.54 +
4.55 +void checkCirculationCompile()
4.56 +{
4.57 + typedef int VType;
4.58 + typedef concepts::Digraph Digraph;
4.59 +
4.60 + typedef Digraph::Node Node;
4.61 + typedef Digraph::Arc Arc;
4.62 + typedef concepts::ReadMap<Arc,VType> CapMap;
4.63 + typedef concepts::ReadMap<Node,VType> DeltaMap;
4.64 + typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
4.65 + typedef concepts::WriteMap<Node,bool> BarrierMap;
4.66 +
4.67 + typedef Elevator<Digraph, Digraph::Node> Elev;
4.68 + typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;
4.69 +
4.70 + Digraph g;
4.71 + Node n;
4.72 + Arc a;
4.73 + CapMap lcap, ucap;
4.74 + DeltaMap delta;
4.75 + FlowMap flow;
4.76 + BarrierMap bar;
4.77 +
4.78 + Circulation<Digraph, CapMap, CapMap, DeltaMap>
4.79 + ::SetFlowMap<FlowMap>
4.80 + ::SetElevator<Elev>
4.81 + ::SetStandardElevator<LinkedElev>
4.82 + ::Create circ_test(g,lcap,ucap,delta);
4.83 +
4.84 + circ_test.lowerCapMap(lcap);
4.85 + circ_test.upperCapMap(ucap);
4.86 + circ_test.deltaMap(delta);
4.87 + flow = circ_test.flowMap();
4.88 + circ_test.flowMap(flow);
4.89 +
4.90 + circ_test.init();
4.91 + circ_test.greedyInit();
4.92 + circ_test.start();
4.93 + circ_test.run();
4.94 +
4.95 + circ_test.barrier(n);
4.96 + circ_test.barrierMap(bar);
4.97 + circ_test.flow(a);
4.98 +}
4.99 +
4.100 +template <class G, class LM, class UM, class DM>
4.101 +void checkCirculation(const G& g, const LM& lm, const UM& um,
4.102 + const DM& dm, bool find)
4.103 +{
4.104 + Circulation<G, LM, UM, DM> circ(g, lm, um, dm);
4.105 + bool ret = circ.run();
4.106 + if (find) {
4.107 + check(ret, "A feasible solution should have been found.");
4.108 + check(circ.checkFlow(), "The found flow is corrupt.");
4.109 + check(!circ.checkBarrier(), "A barrier should not have been found.");
4.110 + } else {
4.111 + check(!ret, "A feasible solution should not have been found.");
4.112 + check(circ.checkBarrier(), "The found barrier is corrupt.");
4.113 + }
4.114 +}
4.115 +
4.116 +int main (int, char*[])
4.117 +{
4.118 + typedef ListDigraph Digraph;
4.119 + DIGRAPH_TYPEDEFS(Digraph);
4.120 +
4.121 + Digraph g;
4.122 + IntArcMap lo(g), up(g);
4.123 + IntNodeMap delta(g, 0);
4.124 + Node s, t;
4.125 +
4.126 + std::istringstream input(test_lgf);
4.127 + DigraphReader<Digraph>(g,input).
4.128 + arcMap("lcap", lo).
4.129 + arcMap("ucap", up).
4.130 + node("source",s).
4.131 + node("sink",t).
4.132 + run();
4.133 +
4.134 + delta[s] = 7; delta[t] = -7;
4.135 + checkCirculation(g, lo, up, delta, true);
4.136 +
4.137 + delta[s] = 13; delta[t] = -13;
4.138 + checkCirculation(g, lo, up, delta, true);
4.139 +
4.140 + delta[s] = 6; delta[t] = -6;
4.141 + checkCirculation(g, lo, up, delta, false);
4.142 +
4.143 + delta[s] = 14; delta[t] = -14;
4.144 + checkCirculation(g, lo, up, delta, false);
4.145 +
4.146 + delta[s] = 7; delta[t] = -13;
4.147 + checkCirculation(g, lo, up, delta, true);
4.148 +
4.149 + delta[s] = 5; delta[t] = -15;
4.150 + checkCirculation(g, lo, up, delta, true);
4.151 +
4.152 + delta[s] = 10; delta[t] = -11;
4.153 + checkCirculation(g, lo, up, delta, true);
4.154 +
4.155 + delta[s] = 11; delta[t] = -10;
4.156 + checkCirculation(g, lo, up, delta, false);
4.157 +
4.158 + return 0;
4.159 +}