1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/howard.h Wed Dec 09 11:14:06 2009 +0100
1.3 @@ -0,0 +1,597 @@
1.4 +/* -*- C++ -*-
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_HOWARD_H
1.23 +#define LEMON_HOWARD_H
1.24 +
1.25 +/// \ingroup min_mean_cycle
1.26 +///
1.27 +/// \file
1.28 +/// \brief Howard's algorithm for finding a minimum mean cycle.
1.29 +
1.30 +#include <vector>
1.31 +#include <limits>
1.32 +#include <lemon/core.h>
1.33 +#include <lemon/path.h>
1.34 +#include <lemon/tolerance.h>
1.35 +#include <lemon/connectivity.h>
1.36 +
1.37 +namespace lemon {
1.38 +
1.39 + /// \brief Default traits class of Howard class.
1.40 + ///
1.41 + /// Default traits class of Howard class.
1.42 + /// \tparam GR The type of the digraph.
1.43 + /// \tparam LEN The type of the length map.
1.44 + /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
1.45 +#ifdef DOXYGEN
1.46 + template <typename GR, typename LEN>
1.47 +#else
1.48 + template <typename GR, typename LEN,
1.49 + bool integer = std::numeric_limits<typename LEN::Value>::is_integer>
1.50 +#endif
1.51 + struct HowardDefaultTraits
1.52 + {
1.53 + /// The type of the digraph
1.54 + typedef GR Digraph;
1.55 + /// The type of the length map
1.56 + typedef LEN LengthMap;
1.57 + /// The type of the arc lengths
1.58 + typedef typename LengthMap::Value Value;
1.59 +
1.60 + /// \brief The large value type used for internal computations
1.61 + ///
1.62 + /// The large value type used for internal computations.
1.63 + /// It is \c long \c long if the \c Value type is integer,
1.64 + /// otherwise it is \c double.
1.65 + /// \c Value must be convertible to \c LargeValue.
1.66 + typedef double LargeValue;
1.67 +
1.68 + /// The tolerance type used for internal computations
1.69 + typedef lemon::Tolerance<LargeValue> Tolerance;
1.70 +
1.71 + /// \brief The path type of the found cycles
1.72 + ///
1.73 + /// The path type of the found cycles.
1.74 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
1.75 + /// and it must have an \c addBack() function.
1.76 + typedef lemon::Path<Digraph> Path;
1.77 + };
1.78 +
1.79 + // Default traits class for integer value types
1.80 + template <typename GR, typename LEN>
1.81 + struct HowardDefaultTraits<GR, LEN, true>
1.82 + {
1.83 + typedef GR Digraph;
1.84 + typedef LEN LengthMap;
1.85 + typedef typename LengthMap::Value Value;
1.86 +#ifdef LEMON_HAVE_LONG_LONG
1.87 + typedef long long LargeValue;
1.88 +#else
1.89 + typedef long LargeValue;
1.90 +#endif
1.91 + typedef lemon::Tolerance<LargeValue> Tolerance;
1.92 + typedef lemon::Path<Digraph> Path;
1.93 + };
1.94 +
1.95 +
1.96 + /// \addtogroup min_mean_cycle
1.97 + /// @{
1.98 +
1.99 + /// \brief Implementation of Howard's algorithm for finding a minimum
1.100 + /// mean cycle.
1.101 + ///
1.102 + /// This class implements Howard's policy iteration algorithm for finding
1.103 + /// a directed cycle of minimum mean length (cost) in a digraph
1.104 + /// \ref amo93networkflows, \ref dasdan98minmeancycle.
1.105 + /// This class provides the most efficient algorithm for the
1.106 + /// minimum mean cycle problem, though the best known theoretical
1.107 + /// bound on its running time is exponential.
1.108 + ///
1.109 + /// \tparam GR The type of the digraph the algorithm runs on.
1.110 + /// \tparam LEN The type of the length map. The default
1.111 + /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
1.112 +#ifdef DOXYGEN
1.113 + template <typename GR, typename LEN, typename TR>
1.114 +#else
1.115 + template < typename GR,
1.116 + typename LEN = typename GR::template ArcMap<int>,
1.117 + typename TR = HowardDefaultTraits<GR, LEN> >
1.118 +#endif
1.119 + class Howard
1.120 + {
1.121 + public:
1.122 +
1.123 + /// The type of the digraph
1.124 + typedef typename TR::Digraph Digraph;
1.125 + /// The type of the length map
1.126 + typedef typename TR::LengthMap LengthMap;
1.127 + /// The type of the arc lengths
1.128 + typedef typename TR::Value Value;
1.129 +
1.130 + /// \brief The large value type
1.131 + ///
1.132 + /// The large value type used for internal computations.
1.133 + /// Using the \ref HowardDefaultTraits "default traits class",
1.134 + /// it is \c long \c long if the \c Value type is integer,
1.135 + /// otherwise it is \c double.
1.136 + typedef typename TR::LargeValue LargeValue;
1.137 +
1.138 + /// The tolerance type
1.139 + typedef typename TR::Tolerance Tolerance;
1.140 +
1.141 + /// \brief The path type of the found cycles
1.142 + ///
1.143 + /// The path type of the found cycles.
1.144 + /// Using the \ref HowardDefaultTraits "default traits class",
1.145 + /// it is \ref lemon::Path "Path<Digraph>".
1.146 + typedef typename TR::Path Path;
1.147 +
1.148 + /// The \ref HowardDefaultTraits "traits class" of the algorithm
1.149 + typedef TR Traits;
1.150 +
1.151 + private:
1.152 +
1.153 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.154 +
1.155 + // The digraph the algorithm runs on
1.156 + const Digraph &_gr;
1.157 + // The length of the arcs
1.158 + const LengthMap &_length;
1.159 +
1.160 + // Data for the found cycles
1.161 + bool _curr_found, _best_found;
1.162 + LargeValue _curr_length, _best_length;
1.163 + int _curr_size, _best_size;
1.164 + Node _curr_node, _best_node;
1.165 +
1.166 + Path *_cycle_path;
1.167 + bool _local_path;
1.168 +
1.169 + // Internal data used by the algorithm
1.170 + typename Digraph::template NodeMap<Arc> _policy;
1.171 + typename Digraph::template NodeMap<bool> _reached;
1.172 + typename Digraph::template NodeMap<int> _level;
1.173 + typename Digraph::template NodeMap<LargeValue> _dist;
1.174 +
1.175 + // Data for storing the strongly connected components
1.176 + int _comp_num;
1.177 + typename Digraph::template NodeMap<int> _comp;
1.178 + std::vector<std::vector<Node> > _comp_nodes;
1.179 + std::vector<Node>* _nodes;
1.180 + typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;
1.181 +
1.182 + // Queue used for BFS search
1.183 + std::vector<Node> _queue;
1.184 + int _qfront, _qback;
1.185 +
1.186 + Tolerance _tolerance;
1.187 +
1.188 + // Infinite constant
1.189 + const LargeValue INF;
1.190 +
1.191 + public:
1.192 +
1.193 + /// \name Named Template Parameters
1.194 + /// @{
1.195 +
1.196 + template <typename T>
1.197 + struct SetLargeValueTraits : public Traits {
1.198 + typedef T LargeValue;
1.199 + typedef lemon::Tolerance<T> Tolerance;
1.200 + };
1.201 +
1.202 + /// \brief \ref named-templ-param "Named parameter" for setting
1.203 + /// \c LargeValue type.
1.204 + ///
1.205 + /// \ref named-templ-param "Named parameter" for setting \c LargeValue
1.206 + /// type. It is used for internal computations in the algorithm.
1.207 + template <typename T>
1.208 + struct SetLargeValue
1.209 + : public Howard<GR, LEN, SetLargeValueTraits<T> > {
1.210 + typedef Howard<GR, LEN, SetLargeValueTraits<T> > Create;
1.211 + };
1.212 +
1.213 + template <typename T>
1.214 + struct SetPathTraits : public Traits {
1.215 + typedef T Path;
1.216 + };
1.217 +
1.218 + /// \brief \ref named-templ-param "Named parameter" for setting
1.219 + /// \c %Path type.
1.220 + ///
1.221 + /// \ref named-templ-param "Named parameter" for setting the \c %Path
1.222 + /// type of the found cycles.
1.223 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
1.224 + /// and it must have an \c addBack() function.
1.225 + template <typename T>
1.226 + struct SetPath
1.227 + : public Howard<GR, LEN, SetPathTraits<T> > {
1.228 + typedef Howard<GR, LEN, SetPathTraits<T> > Create;
1.229 + };
1.230 +
1.231 + /// @}
1.232 +
1.233 + public:
1.234 +
1.235 + /// \brief Constructor.
1.236 + ///
1.237 + /// The constructor of the class.
1.238 + ///
1.239 + /// \param digraph The digraph the algorithm runs on.
1.240 + /// \param length The lengths (costs) of the arcs.
1.241 + Howard( const Digraph &digraph,
1.242 + const LengthMap &length ) :
1.243 + _gr(digraph), _length(length), _best_found(false),
1.244 + _best_length(0), _best_size(1), _cycle_path(NULL), _local_path(false),
1.245 + _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),
1.246 + _comp(digraph), _in_arcs(digraph),
1.247 + INF(std::numeric_limits<LargeValue>::has_infinity ?
1.248 + std::numeric_limits<LargeValue>::infinity() :
1.249 + std::numeric_limits<LargeValue>::max())
1.250 + {}
1.251 +
1.252 + /// Destructor.
1.253 + ~Howard() {
1.254 + if (_local_path) delete _cycle_path;
1.255 + }
1.256 +
1.257 + /// \brief Set the path structure for storing the found cycle.
1.258 + ///
1.259 + /// This function sets an external path structure for storing the
1.260 + /// found cycle.
1.261 + ///
1.262 + /// If you don't call this function before calling \ref run() or
1.263 + /// \ref findMinMean(), it will allocate a local \ref Path "path"
1.264 + /// structure. The destuctor deallocates this automatically
1.265 + /// allocated object, of course.
1.266 + ///
1.267 + /// \note The algorithm calls only the \ref lemon::Path::addBack()
1.268 + /// "addBack()" function of the given path structure.
1.269 + ///
1.270 + /// \return <tt>(*this)</tt>
1.271 + Howard& cycle(Path &path) {
1.272 + if (_local_path) {
1.273 + delete _cycle_path;
1.274 + _local_path = false;
1.275 + }
1.276 + _cycle_path = &path;
1.277 + return *this;
1.278 + }
1.279 +
1.280 + /// \brief Set the tolerance used by the algorithm.
1.281 + ///
1.282 + /// This function sets the tolerance object used by the algorithm.
1.283 + ///
1.284 + /// \return <tt>(*this)</tt>
1.285 + Howard& tolerance(const Tolerance& tolerance) {
1.286 + _tolerance = tolerance;
1.287 + return *this;
1.288 + }
1.289 +
1.290 + /// \brief Return a const reference to the tolerance.
1.291 + ///
1.292 + /// This function returns a const reference to the tolerance object
1.293 + /// used by the algorithm.
1.294 + const Tolerance& tolerance() const {
1.295 + return _tolerance;
1.296 + }
1.297 +
1.298 + /// \name Execution control
1.299 + /// The simplest way to execute the algorithm is to call the \ref run()
1.300 + /// function.\n
1.301 + /// If you only need the minimum mean length, you may call
1.302 + /// \ref findMinMean().
1.303 +
1.304 + /// @{
1.305 +
1.306 + /// \brief Run the algorithm.
1.307 + ///
1.308 + /// This function runs the algorithm.
1.309 + /// It can be called more than once (e.g. if the underlying digraph
1.310 + /// and/or the arc lengths have been modified).
1.311 + ///
1.312 + /// \return \c true if a directed cycle exists in the digraph.
1.313 + ///
1.314 + /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
1.315 + /// \code
1.316 + /// return mmc.findMinMean() && mmc.findCycle();
1.317 + /// \endcode
1.318 + bool run() {
1.319 + return findMinMean() && findCycle();
1.320 + }
1.321 +
1.322 + /// \brief Find the minimum cycle mean.
1.323 + ///
1.324 + /// This function finds the minimum mean length of the directed
1.325 + /// cycles in the digraph.
1.326 + ///
1.327 + /// \return \c true if a directed cycle exists in the digraph.
1.328 + bool findMinMean() {
1.329 + // Initialize and find strongly connected components
1.330 + init();
1.331 + findComponents();
1.332 +
1.333 + // Find the minimum cycle mean in the components
1.334 + for (int comp = 0; comp < _comp_num; ++comp) {
1.335 + // Find the minimum mean cycle in the current component
1.336 + if (!buildPolicyGraph(comp)) continue;
1.337 + while (true) {
1.338 + findPolicyCycle();
1.339 + if (!computeNodeDistances()) break;
1.340 + }
1.341 + // Update the best cycle (global minimum mean cycle)
1.342 + if ( _curr_found && (!_best_found ||
1.343 + _curr_length * _best_size < _best_length * _curr_size) ) {
1.344 + _best_found = true;
1.345 + _best_length = _curr_length;
1.346 + _best_size = _curr_size;
1.347 + _best_node = _curr_node;
1.348 + }
1.349 + }
1.350 + return _best_found;
1.351 + }
1.352 +
1.353 + /// \brief Find a minimum mean directed cycle.
1.354 + ///
1.355 + /// This function finds a directed cycle of minimum mean length
1.356 + /// in the digraph using the data computed by findMinMean().
1.357 + ///
1.358 + /// \return \c true if a directed cycle exists in the digraph.
1.359 + ///
1.360 + /// \pre \ref findMinMean() must be called before using this function.
1.361 + bool findCycle() {
1.362 + if (!_best_found) return false;
1.363 + _cycle_path->addBack(_policy[_best_node]);
1.364 + for ( Node v = _best_node;
1.365 + (v = _gr.target(_policy[v])) != _best_node; ) {
1.366 + _cycle_path->addBack(_policy[v]);
1.367 + }
1.368 + return true;
1.369 + }
1.370 +
1.371 + /// @}
1.372 +
1.373 + /// \name Query Functions
1.374 + /// The results of the algorithm can be obtained using these
1.375 + /// functions.\n
1.376 + /// The algorithm should be executed before using them.
1.377 +
1.378 + /// @{
1.379 +
1.380 + /// \brief Return the total length of the found cycle.
1.381 + ///
1.382 + /// This function returns the total length of the found cycle.
1.383 + ///
1.384 + /// \pre \ref run() or \ref findMinMean() must be called before
1.385 + /// using this function.
1.386 + LargeValue cycleLength() const {
1.387 + return _best_length;
1.388 + }
1.389 +
1.390 + /// \brief Return the number of arcs on the found cycle.
1.391 + ///
1.392 + /// This function returns the number of arcs on the found cycle.
1.393 + ///
1.394 + /// \pre \ref run() or \ref findMinMean() must be called before
1.395 + /// using this function.
1.396 + int cycleArcNum() const {
1.397 + return _best_size;
1.398 + }
1.399 +
1.400 + /// \brief Return the mean length of the found cycle.
1.401 + ///
1.402 + /// This function returns the mean length of the found cycle.
1.403 + ///
1.404 + /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
1.405 + /// following code.
1.406 + /// \code
1.407 + /// return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
1.408 + /// \endcode
1.409 + ///
1.410 + /// \pre \ref run() or \ref findMinMean() must be called before
1.411 + /// using this function.
1.412 + double cycleMean() const {
1.413 + return static_cast<double>(_best_length) / _best_size;
1.414 + }
1.415 +
1.416 + /// \brief Return the found cycle.
1.417 + ///
1.418 + /// This function returns a const reference to the path structure
1.419 + /// storing the found cycle.
1.420 + ///
1.421 + /// \pre \ref run() or \ref findCycle() must be called before using
1.422 + /// this function.
1.423 + const Path& cycle() const {
1.424 + return *_cycle_path;
1.425 + }
1.426 +
1.427 + ///@}
1.428 +
1.429 + private:
1.430 +
1.431 + // Initialize
1.432 + void init() {
1.433 + if (!_cycle_path) {
1.434 + _local_path = true;
1.435 + _cycle_path = new Path;
1.436 + }
1.437 + _queue.resize(countNodes(_gr));
1.438 + _best_found = false;
1.439 + _best_length = 0;
1.440 + _best_size = 1;
1.441 + _cycle_path->clear();
1.442 + }
1.443 +
1.444 + // Find strongly connected components and initialize _comp_nodes
1.445 + // and _in_arcs
1.446 + void findComponents() {
1.447 + _comp_num = stronglyConnectedComponents(_gr, _comp);
1.448 + _comp_nodes.resize(_comp_num);
1.449 + if (_comp_num == 1) {
1.450 + _comp_nodes[0].clear();
1.451 + for (NodeIt n(_gr); n != INVALID; ++n) {
1.452 + _comp_nodes[0].push_back(n);
1.453 + _in_arcs[n].clear();
1.454 + for (InArcIt a(_gr, n); a != INVALID; ++a) {
1.455 + _in_arcs[n].push_back(a);
1.456 + }
1.457 + }
1.458 + } else {
1.459 + for (int i = 0; i < _comp_num; ++i)
1.460 + _comp_nodes[i].clear();
1.461 + for (NodeIt n(_gr); n != INVALID; ++n) {
1.462 + int k = _comp[n];
1.463 + _comp_nodes[k].push_back(n);
1.464 + _in_arcs[n].clear();
1.465 + for (InArcIt a(_gr, n); a != INVALID; ++a) {
1.466 + if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);
1.467 + }
1.468 + }
1.469 + }
1.470 + }
1.471 +
1.472 + // Build the policy graph in the given strongly connected component
1.473 + // (the out-degree of every node is 1)
1.474 + bool buildPolicyGraph(int comp) {
1.475 + _nodes = &(_comp_nodes[comp]);
1.476 + if (_nodes->size() < 1 ||
1.477 + (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {
1.478 + return false;
1.479 + }
1.480 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.481 + _dist[(*_nodes)[i]] = INF;
1.482 + }
1.483 + Node u, v;
1.484 + Arc e;
1.485 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.486 + v = (*_nodes)[i];
1.487 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
1.488 + e = _in_arcs[v][j];
1.489 + u = _gr.source(e);
1.490 + if (_length[e] < _dist[u]) {
1.491 + _dist[u] = _length[e];
1.492 + _policy[u] = e;
1.493 + }
1.494 + }
1.495 + }
1.496 + return true;
1.497 + }
1.498 +
1.499 + // Find the minimum mean cycle in the policy graph
1.500 + void findPolicyCycle() {
1.501 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.502 + _level[(*_nodes)[i]] = -1;
1.503 + }
1.504 + LargeValue clength;
1.505 + int csize;
1.506 + Node u, v;
1.507 + _curr_found = false;
1.508 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.509 + u = (*_nodes)[i];
1.510 + if (_level[u] >= 0) continue;
1.511 + for (; _level[u] < 0; u = _gr.target(_policy[u])) {
1.512 + _level[u] = i;
1.513 + }
1.514 + if (_level[u] == i) {
1.515 + // A cycle is found
1.516 + clength = _length[_policy[u]];
1.517 + csize = 1;
1.518 + for (v = u; (v = _gr.target(_policy[v])) != u; ) {
1.519 + clength += _length[_policy[v]];
1.520 + ++csize;
1.521 + }
1.522 + if ( !_curr_found ||
1.523 + (clength * _curr_size < _curr_length * csize) ) {
1.524 + _curr_found = true;
1.525 + _curr_length = clength;
1.526 + _curr_size = csize;
1.527 + _curr_node = u;
1.528 + }
1.529 + }
1.530 + }
1.531 + }
1.532 +
1.533 + // Contract the policy graph and compute node distances
1.534 + bool computeNodeDistances() {
1.535 + // Find the component of the main cycle and compute node distances
1.536 + // using reverse BFS
1.537 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.538 + _reached[(*_nodes)[i]] = false;
1.539 + }
1.540 + _qfront = _qback = 0;
1.541 + _queue[0] = _curr_node;
1.542 + _reached[_curr_node] = true;
1.543 + _dist[_curr_node] = 0;
1.544 + Node u, v;
1.545 + Arc e;
1.546 + while (_qfront <= _qback) {
1.547 + v = _queue[_qfront++];
1.548 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
1.549 + e = _in_arcs[v][j];
1.550 + u = _gr.source(e);
1.551 + if (_policy[u] == e && !_reached[u]) {
1.552 + _reached[u] = true;
1.553 + _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
1.554 + _queue[++_qback] = u;
1.555 + }
1.556 + }
1.557 + }
1.558 +
1.559 + // Connect all other nodes to this component and compute node
1.560 + // distances using reverse BFS
1.561 + _qfront = 0;
1.562 + while (_qback < int(_nodes->size())-1) {
1.563 + v = _queue[_qfront++];
1.564 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
1.565 + e = _in_arcs[v][j];
1.566 + u = _gr.source(e);
1.567 + if (!_reached[u]) {
1.568 + _reached[u] = true;
1.569 + _policy[u] = e;
1.570 + _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
1.571 + _queue[++_qback] = u;
1.572 + }
1.573 + }
1.574 + }
1.575 +
1.576 + // Improve node distances
1.577 + bool improved = false;
1.578 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.579 + v = (*_nodes)[i];
1.580 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
1.581 + e = _in_arcs[v][j];
1.582 + u = _gr.source(e);
1.583 + LargeValue delta = _dist[v] + _length[e] * _curr_size - _curr_length;
1.584 + if (_tolerance.less(delta, _dist[u])) {
1.585 + _dist[u] = delta;
1.586 + _policy[u] = e;
1.587 + improved = true;
1.588 + }
1.589 + }
1.590 + }
1.591 + return improved;
1.592 + }
1.593 +
1.594 + }; //class Howard
1.595 +
1.596 + ///@}
1.597 +
1.598 +} //namespace lemon
1.599 +
1.600 +#endif //LEMON_HOWARD_H