1.1 --- a/lemon/Makefile.am Mon Mar 08 08:33:41 2010 +0100
1.2 +++ b/lemon/Makefile.am Sat Mar 13 22:01:38 2010 +0100
1.3 @@ -89,10 +89,10 @@
1.4 lemon/gomory_hu.h \
1.5 lemon/graph_to_eps.h \
1.6 lemon/grid_graph.h \
1.7 - lemon/hartmann_orlin.h \
1.8 - lemon/howard.h \
1.9 + lemon/hartmann_orlin_mmc.h \
1.10 + lemon/howard_mmc.h \
1.11 lemon/hypercube_graph.h \
1.12 - lemon/karp.h \
1.13 + lemon/karp_mmc.h \
1.14 lemon/kruskal.h \
1.15 lemon/hao_orlin.h \
1.16 lemon/lgf_reader.h \
2.1 --- a/lemon/cycle_canceling.h Mon Mar 08 08:33:41 2010 +0100
2.2 +++ b/lemon/cycle_canceling.h Sat Mar 13 22:01:38 2010 +0100
2.3 @@ -34,7 +34,7 @@
2.4 #include <lemon/adaptors.h>
2.5 #include <lemon/circulation.h>
2.6 #include <lemon/bellman_ford.h>
2.7 -#include <lemon/howard.h>
2.8 +#include <lemon/howard_mmc.h>
2.9
2.10 namespace lemon {
2.11
2.12 @@ -924,14 +924,14 @@
2.13 void startMinMeanCycleCanceling() {
2.14 typedef SimplePath<StaticDigraph> SPath;
2.15 typedef typename SPath::ArcIt SPathArcIt;
2.16 - typedef typename Howard<StaticDigraph, CostArcMap>
2.17 + typedef typename HowardMmc<StaticDigraph, CostArcMap>
2.18 ::template SetPath<SPath>::Create MMC;
2.19
2.20 SPath cycle;
2.21 MMC mmc(_sgr, _cost_map);
2.22 mmc.cycle(cycle);
2.23 buildResidualNetwork();
2.24 - while (mmc.findMinMean() && mmc.cycleLength() < 0) {
2.25 + while (mmc.findCycleMean() && mmc.cycleCost() < 0) {
2.26 // Find the cycle
2.27 mmc.findCycle();
2.28
2.29 @@ -1132,17 +1132,17 @@
2.30 }
2.31 }
2.32 } else {
2.33 - typedef Howard<StaticDigraph, CostArcMap> MMC;
2.34 + typedef HowardMmc<StaticDigraph, CostArcMap> MMC;
2.35 typedef typename BellmanFord<StaticDigraph, CostArcMap>
2.36 ::template SetDistMap<CostNodeMap>::Create BF;
2.37
2.38 // Set epsilon to the minimum cycle mean
2.39 buildResidualNetwork();
2.40 MMC mmc(_sgr, _cost_map);
2.41 - mmc.findMinMean();
2.42 + mmc.findCycleMean();
2.43 epsilon = -mmc.cycleMean();
2.44 - Cost cycle_cost = mmc.cycleLength();
2.45 - int cycle_size = mmc.cycleArcNum();
2.46 + Cost cycle_cost = mmc.cycleCost();
2.47 + int cycle_size = mmc.cycleSize();
2.48
2.49 // Compute feasible potentials for the current epsilon
2.50 for (int i = 0; i != int(_cost_vec.size()); ++i) {
3.1 --- a/lemon/hartmann_orlin.h Mon Mar 08 08:33:41 2010 +0100
3.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
3.3 @@ -1,650 +0,0 @@
3.4 -/* -*- C++ -*-
3.5 - *
3.6 - * This file is a part of LEMON, a generic C++ optimization library
3.7 - *
3.8 - * Copyright (C) 2003-2008
3.9 - * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
3.10 - * (Egervary Research Group on Combinatorial Optimization, EGRES).
3.11 - *
3.12 - * Permission to use, modify and distribute this software is granted
3.13 - * provided that this copyright notice appears in all copies. For
3.14 - * precise terms see the accompanying LICENSE file.
3.15 - *
3.16 - * This software is provided "AS IS" with no warranty of any kind,
3.17 - * express or implied, and with no claim as to its suitability for any
3.18 - * purpose.
3.19 - *
3.20 - */
3.21 -
3.22 -#ifndef LEMON_HARTMANN_ORLIN_H
3.23 -#define LEMON_HARTMANN_ORLIN_H
3.24 -
3.25 -/// \ingroup min_mean_cycle
3.26 -///
3.27 -/// \file
3.28 -/// \brief Hartmann-Orlin's algorithm for finding a minimum mean cycle.
3.29 -
3.30 -#include <vector>
3.31 -#include <limits>
3.32 -#include <lemon/core.h>
3.33 -#include <lemon/path.h>
3.34 -#include <lemon/tolerance.h>
3.35 -#include <lemon/connectivity.h>
3.36 -
3.37 -namespace lemon {
3.38 -
3.39 - /// \brief Default traits class of HartmannOrlin algorithm.
3.40 - ///
3.41 - /// Default traits class of HartmannOrlin algorithm.
3.42 - /// \tparam GR The type of the digraph.
3.43 - /// \tparam LEN The type of the length map.
3.44 - /// It must conform to the \ref concepts::Rea_data "Rea_data" concept.
3.45 -#ifdef DOXYGEN
3.46 - template <typename GR, typename LEN>
3.47 -#else
3.48 - template <typename GR, typename LEN,
3.49 - bool integer = std::numeric_limits<typename LEN::Value>::is_integer>
3.50 -#endif
3.51 - struct HartmannOrlinDefaultTraits
3.52 - {
3.53 - /// The type of the digraph
3.54 - typedef GR Digraph;
3.55 - /// The type of the length map
3.56 - typedef LEN LengthMap;
3.57 - /// The type of the arc lengths
3.58 - typedef typename LengthMap::Value Value;
3.59 -
3.60 - /// \brief The large value type used for internal computations
3.61 - ///
3.62 - /// The large value type used for internal computations.
3.63 - /// It is \c long \c long if the \c Value type is integer,
3.64 - /// otherwise it is \c double.
3.65 - /// \c Value must be convertible to \c LargeValue.
3.66 - typedef double LargeValue;
3.67 -
3.68 - /// The tolerance type used for internal computations
3.69 - typedef lemon::Tolerance<LargeValue> Tolerance;
3.70 -
3.71 - /// \brief The path type of the found cycles
3.72 - ///
3.73 - /// The path type of the found cycles.
3.74 - /// It must conform to the \ref lemon::concepts::Path "Path" concept
3.75 - /// and it must have an \c addFront() function.
3.76 - typedef lemon::Path<Digraph> Path;
3.77 - };
3.78 -
3.79 - // Default traits class for integer value types
3.80 - template <typename GR, typename LEN>
3.81 - struct HartmannOrlinDefaultTraits<GR, LEN, true>
3.82 - {
3.83 - typedef GR Digraph;
3.84 - typedef LEN LengthMap;
3.85 - typedef typename LengthMap::Value Value;
3.86 -#ifdef LEMON_HAVE_LONG_LONG
3.87 - typedef long long LargeValue;
3.88 -#else
3.89 - typedef long LargeValue;
3.90 -#endif
3.91 - typedef lemon::Tolerance<LargeValue> Tolerance;
3.92 - typedef lemon::Path<Digraph> Path;
3.93 - };
3.94 -
3.95 -
3.96 - /// \addtogroup min_mean_cycle
3.97 - /// @{
3.98 -
3.99 - /// \brief Implementation of the Hartmann-Orlin algorithm for finding
3.100 - /// a minimum mean cycle.
3.101 - ///
3.102 - /// This class implements the Hartmann-Orlin algorithm for finding
3.103 - /// a directed cycle of minimum mean length (cost) in a digraph
3.104 - /// \ref amo93networkflows, \ref dasdan98minmeancycle.
3.105 - /// It is an improved version of \ref Karp "Karp"'s original algorithm,
3.106 - /// it applies an efficient early termination scheme.
3.107 - /// It runs in time O(ne) and uses space O(n<sup>2</sup>+e).
3.108 - ///
3.109 - /// \tparam GR The type of the digraph the algorithm runs on.
3.110 - /// \tparam LEN The type of the length map. The default
3.111 - /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
3.112 - /// \tparam TR The traits class that defines various types used by the
3.113 - /// algorithm. By default, it is \ref HartmannOrlinDefaultTraits
3.114 - /// "HartmannOrlinDefaultTraits<GR, LEN>".
3.115 - /// In most cases, this parameter should not be set directly,
3.116 - /// consider to use the named template parameters instead.
3.117 -#ifdef DOXYGEN
3.118 - template <typename GR, typename LEN, typename TR>
3.119 -#else
3.120 - template < typename GR,
3.121 - typename LEN = typename GR::template ArcMap<int>,
3.122 - typename TR = HartmannOrlinDefaultTraits<GR, LEN> >
3.123 -#endif
3.124 - class HartmannOrlin
3.125 - {
3.126 - public:
3.127 -
3.128 - /// The type of the digraph
3.129 - typedef typename TR::Digraph Digraph;
3.130 - /// The type of the length map
3.131 - typedef typename TR::LengthMap LengthMap;
3.132 - /// The type of the arc lengths
3.133 - typedef typename TR::Value Value;
3.134 -
3.135 - /// \brief The large value type
3.136 - ///
3.137 - /// The large value type used for internal computations.
3.138 - /// By default, it is \c long \c long if the \c Value type is integer,
3.139 - /// otherwise it is \c double.
3.140 - typedef typename TR::LargeValue LargeValue;
3.141 -
3.142 - /// The tolerance type
3.143 - typedef typename TR::Tolerance Tolerance;
3.144 -
3.145 - /// \brief The path type of the found cycles
3.146 - ///
3.147 - /// The path type of the found cycles.
3.148 - /// Using the \ref HartmannOrlinDefaultTraits "default traits class",
3.149 - /// it is \ref lemon::Path "Path<Digraph>".
3.150 - typedef typename TR::Path Path;
3.151 -
3.152 - /// The \ref HartmannOrlinDefaultTraits "traits class" of the algorithm
3.153 - typedef TR Traits;
3.154 -
3.155 - private:
3.156 -
3.157 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
3.158 -
3.159 - // Data sturcture for path data
3.160 - struct PathData
3.161 - {
3.162 - LargeValue dist;
3.163 - Arc pred;
3.164 - PathData(LargeValue d, Arc p = INVALID) :
3.165 - dist(d), pred(p) {}
3.166 - };
3.167 -
3.168 - typedef typename Digraph::template NodeMap<std::vector<PathData> >
3.169 - PathDataNodeMap;
3.170 -
3.171 - private:
3.172 -
3.173 - // The digraph the algorithm runs on
3.174 - const Digraph &_gr;
3.175 - // The length of the arcs
3.176 - const LengthMap &_length;
3.177 -
3.178 - // Data for storing the strongly connected components
3.179 - int _comp_num;
3.180 - typename Digraph::template NodeMap<int> _comp;
3.181 - std::vector<std::vector<Node> > _comp_nodes;
3.182 - std::vector<Node>* _nodes;
3.183 - typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
3.184 -
3.185 - // Data for the found cycles
3.186 - bool _curr_found, _best_found;
3.187 - LargeValue _curr_length, _best_length;
3.188 - int _curr_size, _best_size;
3.189 - Node _curr_node, _best_node;
3.190 - int _curr_level, _best_level;
3.191 -
3.192 - Path *_cycle_path;
3.193 - bool _local_path;
3.194 -
3.195 - // Node map for storing path data
3.196 - PathDataNodeMap _data;
3.197 - // The processed nodes in the last round
3.198 - std::vector<Node> _process;
3.199 -
3.200 - Tolerance _tolerance;
3.201 -
3.202 - // Infinite constant
3.203 - const LargeValue INF;
3.204 -
3.205 - public:
3.206 -
3.207 - /// \name Named Template Parameters
3.208 - /// @{
3.209 -
3.210 - template <typename T>
3.211 - struct SetLargeValueTraits : public Traits {
3.212 - typedef T LargeValue;
3.213 - typedef lemon::Tolerance<T> Tolerance;
3.214 - };
3.215 -
3.216 - /// \brief \ref named-templ-param "Named parameter" for setting
3.217 - /// \c LargeValue type.
3.218 - ///
3.219 - /// \ref named-templ-param "Named parameter" for setting \c LargeValue
3.220 - /// type. It is used for internal computations in the algorithm.
3.221 - template <typename T>
3.222 - struct SetLargeValue
3.223 - : public HartmannOrlin<GR, LEN, SetLargeValueTraits<T> > {
3.224 - typedef HartmannOrlin<GR, LEN, SetLargeValueTraits<T> > Create;
3.225 - };
3.226 -
3.227 - template <typename T>
3.228 - struct SetPathTraits : public Traits {
3.229 - typedef T Path;
3.230 - };
3.231 -
3.232 - /// \brief \ref named-templ-param "Named parameter" for setting
3.233 - /// \c %Path type.
3.234 - ///
3.235 - /// \ref named-templ-param "Named parameter" for setting the \c %Path
3.236 - /// type of the found cycles.
3.237 - /// It must conform to the \ref lemon::concepts::Path "Path" concept
3.238 - /// and it must have an \c addFront() function.
3.239 - template <typename T>
3.240 - struct SetPath
3.241 - : public HartmannOrlin<GR, LEN, SetPathTraits<T> > {
3.242 - typedef HartmannOrlin<GR, LEN, SetPathTraits<T> > Create;
3.243 - };
3.244 -
3.245 - /// @}
3.246 -
3.247 - protected:
3.248 -
3.249 - HartmannOrlin() {}
3.250 -
3.251 - public:
3.252 -
3.253 - /// \brief Constructor.
3.254 - ///
3.255 - /// The constructor of the class.
3.256 - ///
3.257 - /// \param digraph The digraph the algorithm runs on.
3.258 - /// \param length The lengths (costs) of the arcs.
3.259 - HartmannOrlin( const Digraph &digraph,
3.260 - const LengthMap &length ) :
3.261 - _gr(digraph), _length(length), _comp(digraph), _out_arcs(digraph),
3.262 - _best_found(false), _best_length(0), _best_size(1),
3.263 - _cycle_path(NULL), _local_path(false), _data(digraph),
3.264 - INF(std::numeric_limits<LargeValue>::has_infinity ?
3.265 - std::numeric_limits<LargeValue>::infinity() :
3.266 - std::numeric_limits<LargeValue>::max())
3.267 - {}
3.268 -
3.269 - /// Destructor.
3.270 - ~HartmannOrlin() {
3.271 - if (_local_path) delete _cycle_path;
3.272 - }
3.273 -
3.274 - /// \brief Set the path structure for storing the found cycle.
3.275 - ///
3.276 - /// This function sets an external path structure for storing the
3.277 - /// found cycle.
3.278 - ///
3.279 - /// If you don't call this function before calling \ref run() or
3.280 - /// \ref findMinMean(), it will allocate a local \ref Path "path"
3.281 - /// structure. The destuctor deallocates this automatically
3.282 - /// allocated object, of course.
3.283 - ///
3.284 - /// \note The algorithm calls only the \ref lemon::Path::addFront()
3.285 - /// "addFront()" function of the given path structure.
3.286 - ///
3.287 - /// \return <tt>(*this)</tt>
3.288 - HartmannOrlin& cycle(Path &path) {
3.289 - if (_local_path) {
3.290 - delete _cycle_path;
3.291 - _local_path = false;
3.292 - }
3.293 - _cycle_path = &path;
3.294 - return *this;
3.295 - }
3.296 -
3.297 - /// \brief Set the tolerance used by the algorithm.
3.298 - ///
3.299 - /// This function sets the tolerance object used by the algorithm.
3.300 - ///
3.301 - /// \return <tt>(*this)</tt>
3.302 - HartmannOrlin& tolerance(const Tolerance& tolerance) {
3.303 - _tolerance = tolerance;
3.304 - return *this;
3.305 - }
3.306 -
3.307 - /// \brief Return a const reference to the tolerance.
3.308 - ///
3.309 - /// This function returns a const reference to the tolerance object
3.310 - /// used by the algorithm.
3.311 - const Tolerance& tolerance() const {
3.312 - return _tolerance;
3.313 - }
3.314 -
3.315 - /// \name Execution control
3.316 - /// The simplest way to execute the algorithm is to call the \ref run()
3.317 - /// function.\n
3.318 - /// If you only need the minimum mean length, you may call
3.319 - /// \ref findMinMean().
3.320 -
3.321 - /// @{
3.322 -
3.323 - /// \brief Run the algorithm.
3.324 - ///
3.325 - /// This function runs the algorithm.
3.326 - /// It can be called more than once (e.g. if the underlying digraph
3.327 - /// and/or the arc lengths have been modified).
3.328 - ///
3.329 - /// \return \c true if a directed cycle exists in the digraph.
3.330 - ///
3.331 - /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
3.332 - /// \code
3.333 - /// return mmc.findMinMean() && mmc.findCycle();
3.334 - /// \endcode
3.335 - bool run() {
3.336 - return findMinMean() && findCycle();
3.337 - }
3.338 -
3.339 - /// \brief Find the minimum cycle mean.
3.340 - ///
3.341 - /// This function finds the minimum mean length of the directed
3.342 - /// cycles in the digraph.
3.343 - ///
3.344 - /// \return \c true if a directed cycle exists in the digraph.
3.345 - bool findMinMean() {
3.346 - // Initialization and find strongly connected components
3.347 - init();
3.348 - findComponents();
3.349 -
3.350 - // Find the minimum cycle mean in the components
3.351 - for (int comp = 0; comp < _comp_num; ++comp) {
3.352 - if (!initComponent(comp)) continue;
3.353 - processRounds();
3.354 -
3.355 - // Update the best cycle (global minimum mean cycle)
3.356 - if ( _curr_found && (!_best_found ||
3.357 - _curr_length * _best_size < _best_length * _curr_size) ) {
3.358 - _best_found = true;
3.359 - _best_length = _curr_length;
3.360 - _best_size = _curr_size;
3.361 - _best_node = _curr_node;
3.362 - _best_level = _curr_level;
3.363 - }
3.364 - }
3.365 - return _best_found;
3.366 - }
3.367 -
3.368 - /// \brief Find a minimum mean directed cycle.
3.369 - ///
3.370 - /// This function finds a directed cycle of minimum mean length
3.371 - /// in the digraph using the data computed by findMinMean().
3.372 - ///
3.373 - /// \return \c true if a directed cycle exists in the digraph.
3.374 - ///
3.375 - /// \pre \ref findMinMean() must be called before using this function.
3.376 - bool findCycle() {
3.377 - if (!_best_found) return false;
3.378 - IntNodeMap reached(_gr, -1);
3.379 - int r = _best_level + 1;
3.380 - Node u = _best_node;
3.381 - while (reached[u] < 0) {
3.382 - reached[u] = --r;
3.383 - u = _gr.source(_data[u][r].pred);
3.384 - }
3.385 - r = reached[u];
3.386 - Arc e = _data[u][r].pred;
3.387 - _cycle_path->addFront(e);
3.388 - _best_length = _length[e];
3.389 - _best_size = 1;
3.390 - Node v;
3.391 - while ((v = _gr.source(e)) != u) {
3.392 - e = _data[v][--r].pred;
3.393 - _cycle_path->addFront(e);
3.394 - _best_length += _length[e];
3.395 - ++_best_size;
3.396 - }
3.397 - return true;
3.398 - }
3.399 -
3.400 - /// @}
3.401 -
3.402 - /// \name Query Functions
3.403 - /// The results of the algorithm can be obtained using these
3.404 - /// functions.\n
3.405 - /// The algorithm should be executed before using them.
3.406 -
3.407 - /// @{
3.408 -
3.409 - /// \brief Return the total length of the found cycle.
3.410 - ///
3.411 - /// This function returns the total length of the found cycle.
3.412 - ///
3.413 - /// \pre \ref run() or \ref findMinMean() must be called before
3.414 - /// using this function.
3.415 - Value cycleLength() const {
3.416 - return static_cast<Value>(_best_length);
3.417 - }
3.418 -
3.419 - /// \brief Return the number of arcs on the found cycle.
3.420 - ///
3.421 - /// This function returns the number of arcs on the found cycle.
3.422 - ///
3.423 - /// \pre \ref run() or \ref findMinMean() must be called before
3.424 - /// using this function.
3.425 - int cycleArcNum() const {
3.426 - return _best_size;
3.427 - }
3.428 -
3.429 - /// \brief Return the mean length of the found cycle.
3.430 - ///
3.431 - /// This function returns the mean length of the found cycle.
3.432 - ///
3.433 - /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
3.434 - /// following code.
3.435 - /// \code
3.436 - /// return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
3.437 - /// \endcode
3.438 - ///
3.439 - /// \pre \ref run() or \ref findMinMean() must be called before
3.440 - /// using this function.
3.441 - double cycleMean() const {
3.442 - return static_cast<double>(_best_length) / _best_size;
3.443 - }
3.444 -
3.445 - /// \brief Return the found cycle.
3.446 - ///
3.447 - /// This function returns a const reference to the path structure
3.448 - /// storing the found cycle.
3.449 - ///
3.450 - /// \pre \ref run() or \ref findCycle() must be called before using
3.451 - /// this function.
3.452 - const Path& cycle() const {
3.453 - return *_cycle_path;
3.454 - }
3.455 -
3.456 - ///@}
3.457 -
3.458 - private:
3.459 -
3.460 - // Initialization
3.461 - void init() {
3.462 - if (!_cycle_path) {
3.463 - _local_path = true;
3.464 - _cycle_path = new Path;
3.465 - }
3.466 - _cycle_path->clear();
3.467 - _best_found = false;
3.468 - _best_length = 0;
3.469 - _best_size = 1;
3.470 - _cycle_path->clear();
3.471 - for (NodeIt u(_gr); u != INVALID; ++u)
3.472 - _data[u].clear();
3.473 - }
3.474 -
3.475 - // Find strongly connected components and initialize _comp_nodes
3.476 - // and _out_arcs
3.477 - void findComponents() {
3.478 - _comp_num = stronglyConnectedComponents(_gr, _comp);
3.479 - _comp_nodes.resize(_comp_num);
3.480 - if (_comp_num == 1) {
3.481 - _comp_nodes[0].clear();
3.482 - for (NodeIt n(_gr); n != INVALID; ++n) {
3.483 - _comp_nodes[0].push_back(n);
3.484 - _out_arcs[n].clear();
3.485 - for (OutArcIt a(_gr, n); a != INVALID; ++a) {
3.486 - _out_arcs[n].push_back(a);
3.487 - }
3.488 - }
3.489 - } else {
3.490 - for (int i = 0; i < _comp_num; ++i)
3.491 - _comp_nodes[i].clear();
3.492 - for (NodeIt n(_gr); n != INVALID; ++n) {
3.493 - int k = _comp[n];
3.494 - _comp_nodes[k].push_back(n);
3.495 - _out_arcs[n].clear();
3.496 - for (OutArcIt a(_gr, n); a != INVALID; ++a) {
3.497 - if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
3.498 - }
3.499 - }
3.500 - }
3.501 - }
3.502 -
3.503 - // Initialize path data for the current component
3.504 - bool initComponent(int comp) {
3.505 - _nodes = &(_comp_nodes[comp]);
3.506 - int n = _nodes->size();
3.507 - if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
3.508 - return false;
3.509 - }
3.510 - for (int i = 0; i < n; ++i) {
3.511 - _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
3.512 - }
3.513 - return true;
3.514 - }
3.515 -
3.516 - // Process all rounds of computing path data for the current component.
3.517 - // _data[v][k] is the length of a shortest directed walk from the root
3.518 - // node to node v containing exactly k arcs.
3.519 - void processRounds() {
3.520 - Node start = (*_nodes)[0];
3.521 - _data[start][0] = PathData(0);
3.522 - _process.clear();
3.523 - _process.push_back(start);
3.524 -
3.525 - int k, n = _nodes->size();
3.526 - int next_check = 4;
3.527 - bool terminate = false;
3.528 - for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {
3.529 - processNextBuildRound(k);
3.530 - if (k == next_check || k == n) {
3.531 - terminate = checkTermination(k);
3.532 - next_check = next_check * 3 / 2;
3.533 - }
3.534 - }
3.535 - for ( ; k <= n && !terminate; ++k) {
3.536 - processNextFullRound(k);
3.537 - if (k == next_check || k == n) {
3.538 - terminate = checkTermination(k);
3.539 - next_check = next_check * 3 / 2;
3.540 - }
3.541 - }
3.542 - }
3.543 -
3.544 - // Process one round and rebuild _process
3.545 - void processNextBuildRound(int k) {
3.546 - std::vector<Node> next;
3.547 - Node u, v;
3.548 - Arc e;
3.549 - LargeValue d;
3.550 - for (int i = 0; i < int(_process.size()); ++i) {
3.551 - u = _process[i];
3.552 - for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
3.553 - e = _out_arcs[u][j];
3.554 - v = _gr.target(e);
3.555 - d = _data[u][k-1].dist + _length[e];
3.556 - if (_tolerance.less(d, _data[v][k].dist)) {
3.557 - if (_data[v][k].dist == INF) next.push_back(v);
3.558 - _data[v][k] = PathData(d, e);
3.559 - }
3.560 - }
3.561 - }
3.562 - _process.swap(next);
3.563 - }
3.564 -
3.565 - // Process one round using _nodes instead of _process
3.566 - void processNextFullRound(int k) {
3.567 - Node u, v;
3.568 - Arc e;
3.569 - LargeValue d;
3.570 - for (int i = 0; i < int(_nodes->size()); ++i) {
3.571 - u = (*_nodes)[i];
3.572 - for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
3.573 - e = _out_arcs[u][j];
3.574 - v = _gr.target(e);
3.575 - d = _data[u][k-1].dist + _length[e];
3.576 - if (_tolerance.less(d, _data[v][k].dist)) {
3.577 - _data[v][k] = PathData(d, e);
3.578 - }
3.579 - }
3.580 - }
3.581 - }
3.582 -
3.583 - // Check early termination
3.584 - bool checkTermination(int k) {
3.585 - typedef std::pair<int, int> Pair;
3.586 - typename GR::template NodeMap<Pair> level(_gr, Pair(-1, 0));
3.587 - typename GR::template NodeMap<LargeValue> pi(_gr);
3.588 - int n = _nodes->size();
3.589 - LargeValue length;
3.590 - int size;
3.591 - Node u;
3.592 -
3.593 - // Search for cycles that are already found
3.594 - _curr_found = false;
3.595 - for (int i = 0; i < n; ++i) {
3.596 - u = (*_nodes)[i];
3.597 - if (_data[u][k].dist == INF) continue;
3.598 - for (int j = k; j >= 0; --j) {
3.599 - if (level[u].first == i && level[u].second > 0) {
3.600 - // A cycle is found
3.601 - length = _data[u][level[u].second].dist - _data[u][j].dist;
3.602 - size = level[u].second - j;
3.603 - if (!_curr_found || length * _curr_size < _curr_length * size) {
3.604 - _curr_length = length;
3.605 - _curr_size = size;
3.606 - _curr_node = u;
3.607 - _curr_level = level[u].second;
3.608 - _curr_found = true;
3.609 - }
3.610 - }
3.611 - level[u] = Pair(i, j);
3.612 - if (j != 0) {
3.613 - u = _gr.source(_data[u][j].pred);
3.614 - }
3.615 - }
3.616 - }
3.617 -
3.618 - // If at least one cycle is found, check the optimality condition
3.619 - LargeValue d;
3.620 - if (_curr_found && k < n) {
3.621 - // Find node potentials
3.622 - for (int i = 0; i < n; ++i) {
3.623 - u = (*_nodes)[i];
3.624 - pi[u] = INF;
3.625 - for (int j = 0; j <= k; ++j) {
3.626 - if (_data[u][j].dist < INF) {
3.627 - d = _data[u][j].dist * _curr_size - j * _curr_length;
3.628 - if (_tolerance.less(d, pi[u])) pi[u] = d;
3.629 - }
3.630 - }
3.631 - }
3.632 -
3.633 - // Check the optimality condition for all arcs
3.634 - bool done = true;
3.635 - for (ArcIt a(_gr); a != INVALID; ++a) {
3.636 - if (_tolerance.less(_length[a] * _curr_size - _curr_length,
3.637 - pi[_gr.target(a)] - pi[_gr.source(a)]) ) {
3.638 - done = false;
3.639 - break;
3.640 - }
3.641 - }
3.642 - return done;
3.643 - }
3.644 - return (k == n);
3.645 - }
3.646 -
3.647 - }; //class HartmannOrlin
3.648 -
3.649 - ///@}
3.650 -
3.651 -} //namespace lemon
3.652 -
3.653 -#endif //LEMON_HARTMANN_ORLIN_H
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/lemon/hartmann_orlin_mmc.h Sat Mar 13 22:01:38 2010 +0100
4.3 @@ -0,0 +1,650 @@
4.4 +/* -*- C++ -*-
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 +#ifndef LEMON_HARTMANN_ORLIN_MMC_H
4.23 +#define LEMON_HARTMANN_ORLIN_MMC_H
4.24 +
4.25 +/// \ingroup min_mean_cycle
4.26 +///
4.27 +/// \file
4.28 +/// \brief Hartmann-Orlin's algorithm for finding a minimum mean cycle.
4.29 +
4.30 +#include <vector>
4.31 +#include <limits>
4.32 +#include <lemon/core.h>
4.33 +#include <lemon/path.h>
4.34 +#include <lemon/tolerance.h>
4.35 +#include <lemon/connectivity.h>
4.36 +
4.37 +namespace lemon {
4.38 +
4.39 + /// \brief Default traits class of HartmannOrlinMmc class.
4.40 + ///
4.41 + /// Default traits class of HartmannOrlinMmc class.
4.42 + /// \tparam GR The type of the digraph.
4.43 + /// \tparam CM The type of the cost map.
4.44 + /// It must conform to the \ref concepts::Rea_data "Rea_data" concept.
4.45 +#ifdef DOXYGEN
4.46 + template <typename GR, typename CM>
4.47 +#else
4.48 + template <typename GR, typename CM,
4.49 + bool integer = std::numeric_limits<typename CM::Value>::is_integer>
4.50 +#endif
4.51 + struct HartmannOrlinMmcDefaultTraits
4.52 + {
4.53 + /// The type of the digraph
4.54 + typedef GR Digraph;
4.55 + /// The type of the cost map
4.56 + typedef CM CostMap;
4.57 + /// The type of the arc costs
4.58 + typedef typename CostMap::Value Cost;
4.59 +
4.60 + /// \brief The large cost type used for internal computations
4.61 + ///
4.62 + /// The large cost type used for internal computations.
4.63 + /// It is \c long \c long if the \c Cost type is integer,
4.64 + /// otherwise it is \c double.
4.65 + /// \c Cost must be convertible to \c LargeCost.
4.66 + typedef double LargeCost;
4.67 +
4.68 + /// The tolerance type used for internal computations
4.69 + typedef lemon::Tolerance<LargeCost> Tolerance;
4.70 +
4.71 + /// \brief The path type of the found cycles
4.72 + ///
4.73 + /// The path type of the found cycles.
4.74 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
4.75 + /// and it must have an \c addFront() function.
4.76 + typedef lemon::Path<Digraph> Path;
4.77 + };
4.78 +
4.79 + // Default traits class for integer cost types
4.80 + template <typename GR, typename CM>
4.81 + struct HartmannOrlinMmcDefaultTraits<GR, CM, true>
4.82 + {
4.83 + typedef GR Digraph;
4.84 + typedef CM CostMap;
4.85 + typedef typename CostMap::Value Cost;
4.86 +#ifdef LEMON_HAVE_LONG_LONG
4.87 + typedef long long LargeCost;
4.88 +#else
4.89 + typedef long LargeCost;
4.90 +#endif
4.91 + typedef lemon::Tolerance<LargeCost> Tolerance;
4.92 + typedef lemon::Path<Digraph> Path;
4.93 + };
4.94 +
4.95 +
4.96 + /// \addtogroup min_mean_cycle
4.97 + /// @{
4.98 +
4.99 + /// \brief Implementation of the Hartmann-Orlin algorithm for finding
4.100 + /// a minimum mean cycle.
4.101 + ///
4.102 + /// This class implements the Hartmann-Orlin algorithm for finding
4.103 + /// a directed cycle of minimum mean cost in a digraph
4.104 + /// \ref amo93networkflows, \ref dasdan98minmeancycle.
4.105 + /// It is an improved version of \ref Karp "Karp"'s original algorithm,
4.106 + /// it applies an efficient early termination scheme.
4.107 + /// It runs in time O(ne) and uses space O(n<sup>2</sup>+e).
4.108 + ///
4.109 + /// \tparam GR The type of the digraph the algorithm runs on.
4.110 + /// \tparam CM The type of the cost map. The default
4.111 + /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
4.112 + /// \tparam TR The traits class that defines various types used by the
4.113 + /// algorithm. By default, it is \ref HartmannOrlinMmcDefaultTraits
4.114 + /// "HartmannOrlinMmcDefaultTraits<GR, CM>".
4.115 + /// In most cases, this parameter should not be set directly,
4.116 + /// consider to use the named template parameters instead.
4.117 +#ifdef DOXYGEN
4.118 + template <typename GR, typename CM, typename TR>
4.119 +#else
4.120 + template < typename GR,
4.121 + typename CM = typename GR::template ArcMap<int>,
4.122 + typename TR = HartmannOrlinMmcDefaultTraits<GR, CM> >
4.123 +#endif
4.124 + class HartmannOrlinMmc
4.125 + {
4.126 + public:
4.127 +
4.128 + /// The type of the digraph
4.129 + typedef typename TR::Digraph Digraph;
4.130 + /// The type of the cost map
4.131 + typedef typename TR::CostMap CostMap;
4.132 + /// The type of the arc costs
4.133 + typedef typename TR::Cost Cost;
4.134 +
4.135 + /// \brief The large cost type
4.136 + ///
4.137 + /// The large cost type used for internal computations.
4.138 + /// By default, it is \c long \c long if the \c Cost type is integer,
4.139 + /// otherwise it is \c double.
4.140 + typedef typename TR::LargeCost LargeCost;
4.141 +
4.142 + /// The tolerance type
4.143 + typedef typename TR::Tolerance Tolerance;
4.144 +
4.145 + /// \brief The path type of the found cycles
4.146 + ///
4.147 + /// The path type of the found cycles.
4.148 + /// Using the \ref HartmannOrlinMmcDefaultTraits "default traits class",
4.149 + /// it is \ref lemon::Path "Path<Digraph>".
4.150 + typedef typename TR::Path Path;
4.151 +
4.152 + /// The \ref HartmannOrlinMmcDefaultTraits "traits class" of the algorithm
4.153 + typedef TR Traits;
4.154 +
4.155 + private:
4.156 +
4.157 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
4.158 +
4.159 + // Data sturcture for path data
4.160 + struct PathData
4.161 + {
4.162 + LargeCost dist;
4.163 + Arc pred;
4.164 + PathData(LargeCost d, Arc p = INVALID) :
4.165 + dist(d), pred(p) {}
4.166 + };
4.167 +
4.168 + typedef typename Digraph::template NodeMap<std::vector<PathData> >
4.169 + PathDataNodeMap;
4.170 +
4.171 + private:
4.172 +
4.173 + // The digraph the algorithm runs on
4.174 + const Digraph &_gr;
4.175 + // The cost of the arcs
4.176 + const CostMap &_cost;
4.177 +
4.178 + // Data for storing the strongly connected components
4.179 + int _comp_num;
4.180 + typename Digraph::template NodeMap<int> _comp;
4.181 + std::vector<std::vector<Node> > _comp_nodes;
4.182 + std::vector<Node>* _nodes;
4.183 + typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
4.184 +
4.185 + // Data for the found cycles
4.186 + bool _curr_found, _best_found;
4.187 + LargeCost _curr_cost, _best_cost;
4.188 + int _curr_size, _best_size;
4.189 + Node _curr_node, _best_node;
4.190 + int _curr_level, _best_level;
4.191 +
4.192 + Path *_cycle_path;
4.193 + bool _local_path;
4.194 +
4.195 + // Node map for storing path data
4.196 + PathDataNodeMap _data;
4.197 + // The processed nodes in the last round
4.198 + std::vector<Node> _process;
4.199 +
4.200 + Tolerance _tolerance;
4.201 +
4.202 + // Infinite constant
4.203 + const LargeCost INF;
4.204 +
4.205 + public:
4.206 +
4.207 + /// \name Named Template Parameters
4.208 + /// @{
4.209 +
4.210 + template <typename T>
4.211 + struct SetLargeCostTraits : public Traits {
4.212 + typedef T LargeCost;
4.213 + typedef lemon::Tolerance<T> Tolerance;
4.214 + };
4.215 +
4.216 + /// \brief \ref named-templ-param "Named parameter" for setting
4.217 + /// \c LargeCost type.
4.218 + ///
4.219 + /// \ref named-templ-param "Named parameter" for setting \c LargeCost
4.220 + /// type. It is used for internal computations in the algorithm.
4.221 + template <typename T>
4.222 + struct SetLargeCost
4.223 + : public HartmannOrlinMmc<GR, CM, SetLargeCostTraits<T> > {
4.224 + typedef HartmannOrlinMmc<GR, CM, SetLargeCostTraits<T> > Create;
4.225 + };
4.226 +
4.227 + template <typename T>
4.228 + struct SetPathTraits : public Traits {
4.229 + typedef T Path;
4.230 + };
4.231 +
4.232 + /// \brief \ref named-templ-param "Named parameter" for setting
4.233 + /// \c %Path type.
4.234 + ///
4.235 + /// \ref named-templ-param "Named parameter" for setting the \c %Path
4.236 + /// type of the found cycles.
4.237 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
4.238 + /// and it must have an \c addFront() function.
4.239 + template <typename T>
4.240 + struct SetPath
4.241 + : public HartmannOrlinMmc<GR, CM, SetPathTraits<T> > {
4.242 + typedef HartmannOrlinMmc<GR, CM, SetPathTraits<T> > Create;
4.243 + };
4.244 +
4.245 + /// @}
4.246 +
4.247 + protected:
4.248 +
4.249 + HartmannOrlinMmc() {}
4.250 +
4.251 + public:
4.252 +
4.253 + /// \brief Constructor.
4.254 + ///
4.255 + /// The constructor of the class.
4.256 + ///
4.257 + /// \param digraph The digraph the algorithm runs on.
4.258 + /// \param cost The costs of the arcs.
4.259 + HartmannOrlinMmc( const Digraph &digraph,
4.260 + const CostMap &cost ) :
4.261 + _gr(digraph), _cost(cost), _comp(digraph), _out_arcs(digraph),
4.262 + _best_found(false), _best_cost(0), _best_size(1),
4.263 + _cycle_path(NULL), _local_path(false), _data(digraph),
4.264 + INF(std::numeric_limits<LargeCost>::has_infinity ?
4.265 + std::numeric_limits<LargeCost>::infinity() :
4.266 + std::numeric_limits<LargeCost>::max())
4.267 + {}
4.268 +
4.269 + /// Destructor.
4.270 + ~HartmannOrlinMmc() {
4.271 + if (_local_path) delete _cycle_path;
4.272 + }
4.273 +
4.274 + /// \brief Set the path structure for storing the found cycle.
4.275 + ///
4.276 + /// This function sets an external path structure for storing the
4.277 + /// found cycle.
4.278 + ///
4.279 + /// If you don't call this function before calling \ref run() or
4.280 + /// \ref findCycleMean(), it will allocate a local \ref Path "path"
4.281 + /// structure. The destuctor deallocates this automatically
4.282 + /// allocated object, of course.
4.283 + ///
4.284 + /// \note The algorithm calls only the \ref lemon::Path::addFront()
4.285 + /// "addFront()" function of the given path structure.
4.286 + ///
4.287 + /// \return <tt>(*this)</tt>
4.288 + HartmannOrlinMmc& cycle(Path &path) {
4.289 + if (_local_path) {
4.290 + delete _cycle_path;
4.291 + _local_path = false;
4.292 + }
4.293 + _cycle_path = &path;
4.294 + return *this;
4.295 + }
4.296 +
4.297 + /// \brief Set the tolerance used by the algorithm.
4.298 + ///
4.299 + /// This function sets the tolerance object used by the algorithm.
4.300 + ///
4.301 + /// \return <tt>(*this)</tt>
4.302 + HartmannOrlinMmc& tolerance(const Tolerance& tolerance) {
4.303 + _tolerance = tolerance;
4.304 + return *this;
4.305 + }
4.306 +
4.307 + /// \brief Return a const reference to the tolerance.
4.308 + ///
4.309 + /// This function returns a const reference to the tolerance object
4.310 + /// used by the algorithm.
4.311 + const Tolerance& tolerance() const {
4.312 + return _tolerance;
4.313 + }
4.314 +
4.315 + /// \name Execution control
4.316 + /// The simplest way to execute the algorithm is to call the \ref run()
4.317 + /// function.\n
4.318 + /// If you only need the minimum mean cost, you may call
4.319 + /// \ref findCycleMean().
4.320 +
4.321 + /// @{
4.322 +
4.323 + /// \brief Run the algorithm.
4.324 + ///
4.325 + /// This function runs the algorithm.
4.326 + /// It can be called more than once (e.g. if the underlying digraph
4.327 + /// and/or the arc costs have been modified).
4.328 + ///
4.329 + /// \return \c true if a directed cycle exists in the digraph.
4.330 + ///
4.331 + /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
4.332 + /// \code
4.333 + /// return mmc.findCycleMean() && mmc.findCycle();
4.334 + /// \endcode
4.335 + bool run() {
4.336 + return findCycleMean() && findCycle();
4.337 + }
4.338 +
4.339 + /// \brief Find the minimum cycle mean.
4.340 + ///
4.341 + /// This function finds the minimum mean cost of the directed
4.342 + /// cycles in the digraph.
4.343 + ///
4.344 + /// \return \c true if a directed cycle exists in the digraph.
4.345 + bool findCycleMean() {
4.346 + // Initialization and find strongly connected components
4.347 + init();
4.348 + findComponents();
4.349 +
4.350 + // Find the minimum cycle mean in the components
4.351 + for (int comp = 0; comp < _comp_num; ++comp) {
4.352 + if (!initComponent(comp)) continue;
4.353 + processRounds();
4.354 +
4.355 + // Update the best cycle (global minimum mean cycle)
4.356 + if ( _curr_found && (!_best_found ||
4.357 + _curr_cost * _best_size < _best_cost * _curr_size) ) {
4.358 + _best_found = true;
4.359 + _best_cost = _curr_cost;
4.360 + _best_size = _curr_size;
4.361 + _best_node = _curr_node;
4.362 + _best_level = _curr_level;
4.363 + }
4.364 + }
4.365 + return _best_found;
4.366 + }
4.367 +
4.368 + /// \brief Find a minimum mean directed cycle.
4.369 + ///
4.370 + /// This function finds a directed cycle of minimum mean cost
4.371 + /// in the digraph using the data computed by findCycleMean().
4.372 + ///
4.373 + /// \return \c true if a directed cycle exists in the digraph.
4.374 + ///
4.375 + /// \pre \ref findCycleMean() must be called before using this function.
4.376 + bool findCycle() {
4.377 + if (!_best_found) return false;
4.378 + IntNodeMap reached(_gr, -1);
4.379 + int r = _best_level + 1;
4.380 + Node u = _best_node;
4.381 + while (reached[u] < 0) {
4.382 + reached[u] = --r;
4.383 + u = _gr.source(_data[u][r].pred);
4.384 + }
4.385 + r = reached[u];
4.386 + Arc e = _data[u][r].pred;
4.387 + _cycle_path->addFront(e);
4.388 + _best_cost = _cost[e];
4.389 + _best_size = 1;
4.390 + Node v;
4.391 + while ((v = _gr.source(e)) != u) {
4.392 + e = _data[v][--r].pred;
4.393 + _cycle_path->addFront(e);
4.394 + _best_cost += _cost[e];
4.395 + ++_best_size;
4.396 + }
4.397 + return true;
4.398 + }
4.399 +
4.400 + /// @}
4.401 +
4.402 + /// \name Query Functions
4.403 + /// The results of the algorithm can be obtained using these
4.404 + /// functions.\n
4.405 + /// The algorithm should be executed before using them.
4.406 +
4.407 + /// @{
4.408 +
4.409 + /// \brief Return the total cost of the found cycle.
4.410 + ///
4.411 + /// This function returns the total cost of the found cycle.
4.412 + ///
4.413 + /// \pre \ref run() or \ref findCycleMean() must be called before
4.414 + /// using this function.
4.415 + Cost cycleCost() const {
4.416 + return static_cast<Cost>(_best_cost);
4.417 + }
4.418 +
4.419 + /// \brief Return the number of arcs on the found cycle.
4.420 + ///
4.421 + /// This function returns the number of arcs on the found cycle.
4.422 + ///
4.423 + /// \pre \ref run() or \ref findCycleMean() must be called before
4.424 + /// using this function.
4.425 + int cycleSize() const {
4.426 + return _best_size;
4.427 + }
4.428 +
4.429 + /// \brief Return the mean cost of the found cycle.
4.430 + ///
4.431 + /// This function returns the mean cost of the found cycle.
4.432 + ///
4.433 + /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
4.434 + /// following code.
4.435 + /// \code
4.436 + /// return static_cast<double>(alg.cycleCost()) / alg.cycleSize();
4.437 + /// \endcode
4.438 + ///
4.439 + /// \pre \ref run() or \ref findCycleMean() must be called before
4.440 + /// using this function.
4.441 + double cycleMean() const {
4.442 + return static_cast<double>(_best_cost) / _best_size;
4.443 + }
4.444 +
4.445 + /// \brief Return the found cycle.
4.446 + ///
4.447 + /// This function returns a const reference to the path structure
4.448 + /// storing the found cycle.
4.449 + ///
4.450 + /// \pre \ref run() or \ref findCycle() must be called before using
4.451 + /// this function.
4.452 + const Path& cycle() const {
4.453 + return *_cycle_path;
4.454 + }
4.455 +
4.456 + ///@}
4.457 +
4.458 + private:
4.459 +
4.460 + // Initialization
4.461 + void init() {
4.462 + if (!_cycle_path) {
4.463 + _local_path = true;
4.464 + _cycle_path = new Path;
4.465 + }
4.466 + _cycle_path->clear();
4.467 + _best_found = false;
4.468 + _best_cost = 0;
4.469 + _best_size = 1;
4.470 + _cycle_path->clear();
4.471 + for (NodeIt u(_gr); u != INVALID; ++u)
4.472 + _data[u].clear();
4.473 + }
4.474 +
4.475 + // Find strongly connected components and initialize _comp_nodes
4.476 + // and _out_arcs
4.477 + void findComponents() {
4.478 + _comp_num = stronglyConnectedComponents(_gr, _comp);
4.479 + _comp_nodes.resize(_comp_num);
4.480 + if (_comp_num == 1) {
4.481 + _comp_nodes[0].clear();
4.482 + for (NodeIt n(_gr); n != INVALID; ++n) {
4.483 + _comp_nodes[0].push_back(n);
4.484 + _out_arcs[n].clear();
4.485 + for (OutArcIt a(_gr, n); a != INVALID; ++a) {
4.486 + _out_arcs[n].push_back(a);
4.487 + }
4.488 + }
4.489 + } else {
4.490 + for (int i = 0; i < _comp_num; ++i)
4.491 + _comp_nodes[i].clear();
4.492 + for (NodeIt n(_gr); n != INVALID; ++n) {
4.493 + int k = _comp[n];
4.494 + _comp_nodes[k].push_back(n);
4.495 + _out_arcs[n].clear();
4.496 + for (OutArcIt a(_gr, n); a != INVALID; ++a) {
4.497 + if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
4.498 + }
4.499 + }
4.500 + }
4.501 + }
4.502 +
4.503 + // Initialize path data for the current component
4.504 + bool initComponent(int comp) {
4.505 + _nodes = &(_comp_nodes[comp]);
4.506 + int n = _nodes->size();
4.507 + if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
4.508 + return false;
4.509 + }
4.510 + for (int i = 0; i < n; ++i) {
4.511 + _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
4.512 + }
4.513 + return true;
4.514 + }
4.515 +
4.516 + // Process all rounds of computing path data for the current component.
4.517 + // _data[v][k] is the cost of a shortest directed walk from the root
4.518 + // node to node v containing exactly k arcs.
4.519 + void processRounds() {
4.520 + Node start = (*_nodes)[0];
4.521 + _data[start][0] = PathData(0);
4.522 + _process.clear();
4.523 + _process.push_back(start);
4.524 +
4.525 + int k, n = _nodes->size();
4.526 + int next_check = 4;
4.527 + bool terminate = false;
4.528 + for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {
4.529 + processNextBuildRound(k);
4.530 + if (k == next_check || k == n) {
4.531 + terminate = checkTermination(k);
4.532 + next_check = next_check * 3 / 2;
4.533 + }
4.534 + }
4.535 + for ( ; k <= n && !terminate; ++k) {
4.536 + processNextFullRound(k);
4.537 + if (k == next_check || k == n) {
4.538 + terminate = checkTermination(k);
4.539 + next_check = next_check * 3 / 2;
4.540 + }
4.541 + }
4.542 + }
4.543 +
4.544 + // Process one round and rebuild _process
4.545 + void processNextBuildRound(int k) {
4.546 + std::vector<Node> next;
4.547 + Node u, v;
4.548 + Arc e;
4.549 + LargeCost d;
4.550 + for (int i = 0; i < int(_process.size()); ++i) {
4.551 + u = _process[i];
4.552 + for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
4.553 + e = _out_arcs[u][j];
4.554 + v = _gr.target(e);
4.555 + d = _data[u][k-1].dist + _cost[e];
4.556 + if (_tolerance.less(d, _data[v][k].dist)) {
4.557 + if (_data[v][k].dist == INF) next.push_back(v);
4.558 + _data[v][k] = PathData(d, e);
4.559 + }
4.560 + }
4.561 + }
4.562 + _process.swap(next);
4.563 + }
4.564 +
4.565 + // Process one round using _nodes instead of _process
4.566 + void processNextFullRound(int k) {
4.567 + Node u, v;
4.568 + Arc e;
4.569 + LargeCost d;
4.570 + for (int i = 0; i < int(_nodes->size()); ++i) {
4.571 + u = (*_nodes)[i];
4.572 + for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
4.573 + e = _out_arcs[u][j];
4.574 + v = _gr.target(e);
4.575 + d = _data[u][k-1].dist + _cost[e];
4.576 + if (_tolerance.less(d, _data[v][k].dist)) {
4.577 + _data[v][k] = PathData(d, e);
4.578 + }
4.579 + }
4.580 + }
4.581 + }
4.582 +
4.583 + // Check early termination
4.584 + bool checkTermination(int k) {
4.585 + typedef std::pair<int, int> Pair;
4.586 + typename GR::template NodeMap<Pair> level(_gr, Pair(-1, 0));
4.587 + typename GR::template NodeMap<LargeCost> pi(_gr);
4.588 + int n = _nodes->size();
4.589 + LargeCost cost;
4.590 + int size;
4.591 + Node u;
4.592 +
4.593 + // Search for cycles that are already found
4.594 + _curr_found = false;
4.595 + for (int i = 0; i < n; ++i) {
4.596 + u = (*_nodes)[i];
4.597 + if (_data[u][k].dist == INF) continue;
4.598 + for (int j = k; j >= 0; --j) {
4.599 + if (level[u].first == i && level[u].second > 0) {
4.600 + // A cycle is found
4.601 + cost = _data[u][level[u].second].dist - _data[u][j].dist;
4.602 + size = level[u].second - j;
4.603 + if (!_curr_found || cost * _curr_size < _curr_cost * size) {
4.604 + _curr_cost = cost;
4.605 + _curr_size = size;
4.606 + _curr_node = u;
4.607 + _curr_level = level[u].second;
4.608 + _curr_found = true;
4.609 + }
4.610 + }
4.611 + level[u] = Pair(i, j);
4.612 + if (j != 0) {
4.613 + u = _gr.source(_data[u][j].pred);
4.614 + }
4.615 + }
4.616 + }
4.617 +
4.618 + // If at least one cycle is found, check the optimality condition
4.619 + LargeCost d;
4.620 + if (_curr_found && k < n) {
4.621 + // Find node potentials
4.622 + for (int i = 0; i < n; ++i) {
4.623 + u = (*_nodes)[i];
4.624 + pi[u] = INF;
4.625 + for (int j = 0; j <= k; ++j) {
4.626 + if (_data[u][j].dist < INF) {
4.627 + d = _data[u][j].dist * _curr_size - j * _curr_cost;
4.628 + if (_tolerance.less(d, pi[u])) pi[u] = d;
4.629 + }
4.630 + }
4.631 + }
4.632 +
4.633 + // Check the optimality condition for all arcs
4.634 + bool done = true;
4.635 + for (ArcIt a(_gr); a != INVALID; ++a) {
4.636 + if (_tolerance.less(_cost[a] * _curr_size - _curr_cost,
4.637 + pi[_gr.target(a)] - pi[_gr.source(a)]) ) {
4.638 + done = false;
4.639 + break;
4.640 + }
4.641 + }
4.642 + return done;
4.643 + }
4.644 + return (k == n);
4.645 + }
4.646 +
4.647 + }; //class HartmannOrlinMmc
4.648 +
4.649 + ///@}
4.650 +
4.651 +} //namespace lemon
4.652 +
4.653 +#endif //LEMON_HARTMANN_ORLIN_MMC_H
5.1 --- a/lemon/howard.h Mon Mar 08 08:33:41 2010 +0100
5.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
5.3 @@ -1,605 +0,0 @@
5.4 -/* -*- C++ -*-
5.5 - *
5.6 - * This file is a part of LEMON, a generic C++ optimization library
5.7 - *
5.8 - * Copyright (C) 2003-2008
5.9 - * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5.10 - * (Egervary Research Group on Combinatorial Optimization, EGRES).
5.11 - *
5.12 - * Permission to use, modify and distribute this software is granted
5.13 - * provided that this copyright notice appears in all copies. For
5.14 - * precise terms see the accompanying LICENSE file.
5.15 - *
5.16 - * This software is provided "AS IS" with no warranty of any kind,
5.17 - * express or implied, and with no claim as to its suitability for any
5.18 - * purpose.
5.19 - *
5.20 - */
5.21 -
5.22 -#ifndef LEMON_HOWARD_H
5.23 -#define LEMON_HOWARD_H
5.24 -
5.25 -/// \ingroup min_mean_cycle
5.26 -///
5.27 -/// \file
5.28 -/// \brief Howard's algorithm for finding a minimum mean cycle.
5.29 -
5.30 -#include <vector>
5.31 -#include <limits>
5.32 -#include <lemon/core.h>
5.33 -#include <lemon/path.h>
5.34 -#include <lemon/tolerance.h>
5.35 -#include <lemon/connectivity.h>
5.36 -
5.37 -namespace lemon {
5.38 -
5.39 - /// \brief Default traits class of Howard class.
5.40 - ///
5.41 - /// Default traits class of Howard class.
5.42 - /// \tparam GR The type of the digraph.
5.43 - /// \tparam LEN The type of the length map.
5.44 - /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
5.45 -#ifdef DOXYGEN
5.46 - template <typename GR, typename LEN>
5.47 -#else
5.48 - template <typename GR, typename LEN,
5.49 - bool integer = std::numeric_limits<typename LEN::Value>::is_integer>
5.50 -#endif
5.51 - struct HowardDefaultTraits
5.52 - {
5.53 - /// The type of the digraph
5.54 - typedef GR Digraph;
5.55 - /// The type of the length map
5.56 - typedef LEN LengthMap;
5.57 - /// The type of the arc lengths
5.58 - typedef typename LengthMap::Value Value;
5.59 -
5.60 - /// \brief The large value type used for internal computations
5.61 - ///
5.62 - /// The large value type used for internal computations.
5.63 - /// It is \c long \c long if the \c Value type is integer,
5.64 - /// otherwise it is \c double.
5.65 - /// \c Value must be convertible to \c LargeValue.
5.66 - typedef double LargeValue;
5.67 -
5.68 - /// The tolerance type used for internal computations
5.69 - typedef lemon::Tolerance<LargeValue> Tolerance;
5.70 -
5.71 - /// \brief The path type of the found cycles
5.72 - ///
5.73 - /// The path type of the found cycles.
5.74 - /// It must conform to the \ref lemon::concepts::Path "Path" concept
5.75 - /// and it must have an \c addBack() function.
5.76 - typedef lemon::Path<Digraph> Path;
5.77 - };
5.78 -
5.79 - // Default traits class for integer value types
5.80 - template <typename GR, typename LEN>
5.81 - struct HowardDefaultTraits<GR, LEN, true>
5.82 - {
5.83 - typedef GR Digraph;
5.84 - typedef LEN LengthMap;
5.85 - typedef typename LengthMap::Value Value;
5.86 -#ifdef LEMON_HAVE_LONG_LONG
5.87 - typedef long long LargeValue;
5.88 -#else
5.89 - typedef long LargeValue;
5.90 -#endif
5.91 - typedef lemon::Tolerance<LargeValue> Tolerance;
5.92 - typedef lemon::Path<Digraph> Path;
5.93 - };
5.94 -
5.95 -
5.96 - /// \addtogroup min_mean_cycle
5.97 - /// @{
5.98 -
5.99 - /// \brief Implementation of Howard's algorithm for finding a minimum
5.100 - /// mean cycle.
5.101 - ///
5.102 - /// This class implements Howard's policy iteration algorithm for finding
5.103 - /// a directed cycle of minimum mean length (cost) in a digraph
5.104 - /// \ref amo93networkflows, \ref dasdan98minmeancycle.
5.105 - /// This class provides the most efficient algorithm for the
5.106 - /// minimum mean cycle problem, though the best known theoretical
5.107 - /// bound on its running time is exponential.
5.108 - ///
5.109 - /// \tparam GR The type of the digraph the algorithm runs on.
5.110 - /// \tparam LEN The type of the length map. The default
5.111 - /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
5.112 - /// \tparam TR The traits class that defines various types used by the
5.113 - /// algorithm. By default, it is \ref HowardDefaultTraits
5.114 - /// "HowardDefaultTraits<GR, LEN>".
5.115 - /// In most cases, this parameter should not be set directly,
5.116 - /// consider to use the named template parameters instead.
5.117 -#ifdef DOXYGEN
5.118 - template <typename GR, typename LEN, typename TR>
5.119 -#else
5.120 - template < typename GR,
5.121 - typename LEN = typename GR::template ArcMap<int>,
5.122 - typename TR = HowardDefaultTraits<GR, LEN> >
5.123 -#endif
5.124 - class Howard
5.125 - {
5.126 - public:
5.127 -
5.128 - /// The type of the digraph
5.129 - typedef typename TR::Digraph Digraph;
5.130 - /// The type of the length map
5.131 - typedef typename TR::LengthMap LengthMap;
5.132 - /// The type of the arc lengths
5.133 - typedef typename TR::Value Value;
5.134 -
5.135 - /// \brief The large value type
5.136 - ///
5.137 - /// The large value type used for internal computations.
5.138 - /// By default, it is \c long \c long if the \c Value type is integer,
5.139 - /// otherwise it is \c double.
5.140 - typedef typename TR::LargeValue LargeValue;
5.141 -
5.142 - /// The tolerance type
5.143 - typedef typename TR::Tolerance Tolerance;
5.144 -
5.145 - /// \brief The path type of the found cycles
5.146 - ///
5.147 - /// The path type of the found cycles.
5.148 - /// Using the \ref HowardDefaultTraits "default traits class",
5.149 - /// it is \ref lemon::Path "Path<Digraph>".
5.150 - typedef typename TR::Path Path;
5.151 -
5.152 - /// The \ref HowardDefaultTraits "traits class" of the algorithm
5.153 - typedef TR Traits;
5.154 -
5.155 - private:
5.156 -
5.157 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
5.158 -
5.159 - // The digraph the algorithm runs on
5.160 - const Digraph &_gr;
5.161 - // The length of the arcs
5.162 - const LengthMap &_length;
5.163 -
5.164 - // Data for the found cycles
5.165 - bool _curr_found, _best_found;
5.166 - LargeValue _curr_length, _best_length;
5.167 - int _curr_size, _best_size;
5.168 - Node _curr_node, _best_node;
5.169 -
5.170 - Path *_cycle_path;
5.171 - bool _local_path;
5.172 -
5.173 - // Internal data used by the algorithm
5.174 - typename Digraph::template NodeMap<Arc> _policy;
5.175 - typename Digraph::template NodeMap<bool> _reached;
5.176 - typename Digraph::template NodeMap<int> _level;
5.177 - typename Digraph::template NodeMap<LargeValue> _dist;
5.178 -
5.179 - // Data for storing the strongly connected components
5.180 - int _comp_num;
5.181 - typename Digraph::template NodeMap<int> _comp;
5.182 - std::vector<std::vector<Node> > _comp_nodes;
5.183 - std::vector<Node>* _nodes;
5.184 - typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;
5.185 -
5.186 - // Queue used for BFS search
5.187 - std::vector<Node> _queue;
5.188 - int _qfront, _qback;
5.189 -
5.190 - Tolerance _tolerance;
5.191 -
5.192 - // Infinite constant
5.193 - const LargeValue INF;
5.194 -
5.195 - public:
5.196 -
5.197 - /// \name Named Template Parameters
5.198 - /// @{
5.199 -
5.200 - template <typename T>
5.201 - struct SetLargeValueTraits : public Traits {
5.202 - typedef T LargeValue;
5.203 - typedef lemon::Tolerance<T> Tolerance;
5.204 - };
5.205 -
5.206 - /// \brief \ref named-templ-param "Named parameter" for setting
5.207 - /// \c LargeValue type.
5.208 - ///
5.209 - /// \ref named-templ-param "Named parameter" for setting \c LargeValue
5.210 - /// type. It is used for internal computations in the algorithm.
5.211 - template <typename T>
5.212 - struct SetLargeValue
5.213 - : public Howard<GR, LEN, SetLargeValueTraits<T> > {
5.214 - typedef Howard<GR, LEN, SetLargeValueTraits<T> > Create;
5.215 - };
5.216 -
5.217 - template <typename T>
5.218 - struct SetPathTraits : public Traits {
5.219 - typedef T Path;
5.220 - };
5.221 -
5.222 - /// \brief \ref named-templ-param "Named parameter" for setting
5.223 - /// \c %Path type.
5.224 - ///
5.225 - /// \ref named-templ-param "Named parameter" for setting the \c %Path
5.226 - /// type of the found cycles.
5.227 - /// It must conform to the \ref lemon::concepts::Path "Path" concept
5.228 - /// and it must have an \c addBack() function.
5.229 - template <typename T>
5.230 - struct SetPath
5.231 - : public Howard<GR, LEN, SetPathTraits<T> > {
5.232 - typedef Howard<GR, LEN, SetPathTraits<T> > Create;
5.233 - };
5.234 -
5.235 - /// @}
5.236 -
5.237 - protected:
5.238 -
5.239 - Howard() {}
5.240 -
5.241 - public:
5.242 -
5.243 - /// \brief Constructor.
5.244 - ///
5.245 - /// The constructor of the class.
5.246 - ///
5.247 - /// \param digraph The digraph the algorithm runs on.
5.248 - /// \param length The lengths (costs) of the arcs.
5.249 - Howard( const Digraph &digraph,
5.250 - const LengthMap &length ) :
5.251 - _gr(digraph), _length(length), _best_found(false),
5.252 - _best_length(0), _best_size(1), _cycle_path(NULL), _local_path(false),
5.253 - _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),
5.254 - _comp(digraph), _in_arcs(digraph),
5.255 - INF(std::numeric_limits<LargeValue>::has_infinity ?
5.256 - std::numeric_limits<LargeValue>::infinity() :
5.257 - std::numeric_limits<LargeValue>::max())
5.258 - {}
5.259 -
5.260 - /// Destructor.
5.261 - ~Howard() {
5.262 - if (_local_path) delete _cycle_path;
5.263 - }
5.264 -
5.265 - /// \brief Set the path structure for storing the found cycle.
5.266 - ///
5.267 - /// This function sets an external path structure for storing the
5.268 - /// found cycle.
5.269 - ///
5.270 - /// If you don't call this function before calling \ref run() or
5.271 - /// \ref findMinMean(), it will allocate a local \ref Path "path"
5.272 - /// structure. The destuctor deallocates this automatically
5.273 - /// allocated object, of course.
5.274 - ///
5.275 - /// \note The algorithm calls only the \ref lemon::Path::addBack()
5.276 - /// "addBack()" function of the given path structure.
5.277 - ///
5.278 - /// \return <tt>(*this)</tt>
5.279 - Howard& cycle(Path &path) {
5.280 - if (_local_path) {
5.281 - delete _cycle_path;
5.282 - _local_path = false;
5.283 - }
5.284 - _cycle_path = &path;
5.285 - return *this;
5.286 - }
5.287 -
5.288 - /// \brief Set the tolerance used by the algorithm.
5.289 - ///
5.290 - /// This function sets the tolerance object used by the algorithm.
5.291 - ///
5.292 - /// \return <tt>(*this)</tt>
5.293 - Howard& tolerance(const Tolerance& tolerance) {
5.294 - _tolerance = tolerance;
5.295 - return *this;
5.296 - }
5.297 -
5.298 - /// \brief Return a const reference to the tolerance.
5.299 - ///
5.300 - /// This function returns a const reference to the tolerance object
5.301 - /// used by the algorithm.
5.302 - const Tolerance& tolerance() const {
5.303 - return _tolerance;
5.304 - }
5.305 -
5.306 - /// \name Execution control
5.307 - /// The simplest way to execute the algorithm is to call the \ref run()
5.308 - /// function.\n
5.309 - /// If you only need the minimum mean length, you may call
5.310 - /// \ref findMinMean().
5.311 -
5.312 - /// @{
5.313 -
5.314 - /// \brief Run the algorithm.
5.315 - ///
5.316 - /// This function runs the algorithm.
5.317 - /// It can be called more than once (e.g. if the underlying digraph
5.318 - /// and/or the arc lengths have been modified).
5.319 - ///
5.320 - /// \return \c true if a directed cycle exists in the digraph.
5.321 - ///
5.322 - /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
5.323 - /// \code
5.324 - /// return mmc.findMinMean() && mmc.findCycle();
5.325 - /// \endcode
5.326 - bool run() {
5.327 - return findMinMean() && findCycle();
5.328 - }
5.329 -
5.330 - /// \brief Find the minimum cycle mean.
5.331 - ///
5.332 - /// This function finds the minimum mean length of the directed
5.333 - /// cycles in the digraph.
5.334 - ///
5.335 - /// \return \c true if a directed cycle exists in the digraph.
5.336 - bool findMinMean() {
5.337 - // Initialize and find strongly connected components
5.338 - init();
5.339 - findComponents();
5.340 -
5.341 - // Find the minimum cycle mean in the components
5.342 - for (int comp = 0; comp < _comp_num; ++comp) {
5.343 - // Find the minimum mean cycle in the current component
5.344 - if (!buildPolicyGraph(comp)) continue;
5.345 - while (true) {
5.346 - findPolicyCycle();
5.347 - if (!computeNodeDistances()) break;
5.348 - }
5.349 - // Update the best cycle (global minimum mean cycle)
5.350 - if ( _curr_found && (!_best_found ||
5.351 - _curr_length * _best_size < _best_length * _curr_size) ) {
5.352 - _best_found = true;
5.353 - _best_length = _curr_length;
5.354 - _best_size = _curr_size;
5.355 - _best_node = _curr_node;
5.356 - }
5.357 - }
5.358 - return _best_found;
5.359 - }
5.360 -
5.361 - /// \brief Find a minimum mean directed cycle.
5.362 - ///
5.363 - /// This function finds a directed cycle of minimum mean length
5.364 - /// in the digraph using the data computed by findMinMean().
5.365 - ///
5.366 - /// \return \c true if a directed cycle exists in the digraph.
5.367 - ///
5.368 - /// \pre \ref findMinMean() must be called before using this function.
5.369 - bool findCycle() {
5.370 - if (!_best_found) return false;
5.371 - _cycle_path->addBack(_policy[_best_node]);
5.372 - for ( Node v = _best_node;
5.373 - (v = _gr.target(_policy[v])) != _best_node; ) {
5.374 - _cycle_path->addBack(_policy[v]);
5.375 - }
5.376 - return true;
5.377 - }
5.378 -
5.379 - /// @}
5.380 -
5.381 - /// \name Query Functions
5.382 - /// The results of the algorithm can be obtained using these
5.383 - /// functions.\n
5.384 - /// The algorithm should be executed before using them.
5.385 -
5.386 - /// @{
5.387 -
5.388 - /// \brief Return the total length of the found cycle.
5.389 - ///
5.390 - /// This function returns the total length of the found cycle.
5.391 - ///
5.392 - /// \pre \ref run() or \ref findMinMean() must be called before
5.393 - /// using this function.
5.394 - Value cycleLength() const {
5.395 - return static_cast<Value>(_best_length);
5.396 - }
5.397 -
5.398 - /// \brief Return the number of arcs on the found cycle.
5.399 - ///
5.400 - /// This function returns the number of arcs on the found cycle.
5.401 - ///
5.402 - /// \pre \ref run() or \ref findMinMean() must be called before
5.403 - /// using this function.
5.404 - int cycleArcNum() const {
5.405 - return _best_size;
5.406 - }
5.407 -
5.408 - /// \brief Return the mean length of the found cycle.
5.409 - ///
5.410 - /// This function returns the mean length of the found cycle.
5.411 - ///
5.412 - /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
5.413 - /// following code.
5.414 - /// \code
5.415 - /// return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
5.416 - /// \endcode
5.417 - ///
5.418 - /// \pre \ref run() or \ref findMinMean() must be called before
5.419 - /// using this function.
5.420 - double cycleMean() const {
5.421 - return static_cast<double>(_best_length) / _best_size;
5.422 - }
5.423 -
5.424 - /// \brief Return the found cycle.
5.425 - ///
5.426 - /// This function returns a const reference to the path structure
5.427 - /// storing the found cycle.
5.428 - ///
5.429 - /// \pre \ref run() or \ref findCycle() must be called before using
5.430 - /// this function.
5.431 - const Path& cycle() const {
5.432 - return *_cycle_path;
5.433 - }
5.434 -
5.435 - ///@}
5.436 -
5.437 - private:
5.438 -
5.439 - // Initialize
5.440 - void init() {
5.441 - if (!_cycle_path) {
5.442 - _local_path = true;
5.443 - _cycle_path = new Path;
5.444 - }
5.445 - _queue.resize(countNodes(_gr));
5.446 - _best_found = false;
5.447 - _best_length = 0;
5.448 - _best_size = 1;
5.449 - _cycle_path->clear();
5.450 - }
5.451 -
5.452 - // Find strongly connected components and initialize _comp_nodes
5.453 - // and _in_arcs
5.454 - void findComponents() {
5.455 - _comp_num = stronglyConnectedComponents(_gr, _comp);
5.456 - _comp_nodes.resize(_comp_num);
5.457 - if (_comp_num == 1) {
5.458 - _comp_nodes[0].clear();
5.459 - for (NodeIt n(_gr); n != INVALID; ++n) {
5.460 - _comp_nodes[0].push_back(n);
5.461 - _in_arcs[n].clear();
5.462 - for (InArcIt a(_gr, n); a != INVALID; ++a) {
5.463 - _in_arcs[n].push_back(a);
5.464 - }
5.465 - }
5.466 - } else {
5.467 - for (int i = 0; i < _comp_num; ++i)
5.468 - _comp_nodes[i].clear();
5.469 - for (NodeIt n(_gr); n != INVALID; ++n) {
5.470 - int k = _comp[n];
5.471 - _comp_nodes[k].push_back(n);
5.472 - _in_arcs[n].clear();
5.473 - for (InArcIt a(_gr, n); a != INVALID; ++a) {
5.474 - if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);
5.475 - }
5.476 - }
5.477 - }
5.478 - }
5.479 -
5.480 - // Build the policy graph in the given strongly connected component
5.481 - // (the out-degree of every node is 1)
5.482 - bool buildPolicyGraph(int comp) {
5.483 - _nodes = &(_comp_nodes[comp]);
5.484 - if (_nodes->size() < 1 ||
5.485 - (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {
5.486 - return false;
5.487 - }
5.488 - for (int i = 0; i < int(_nodes->size()); ++i) {
5.489 - _dist[(*_nodes)[i]] = INF;
5.490 - }
5.491 - Node u, v;
5.492 - Arc e;
5.493 - for (int i = 0; i < int(_nodes->size()); ++i) {
5.494 - v = (*_nodes)[i];
5.495 - for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
5.496 - e = _in_arcs[v][j];
5.497 - u = _gr.source(e);
5.498 - if (_length[e] < _dist[u]) {
5.499 - _dist[u] = _length[e];
5.500 - _policy[u] = e;
5.501 - }
5.502 - }
5.503 - }
5.504 - return true;
5.505 - }
5.506 -
5.507 - // Find the minimum mean cycle in the policy graph
5.508 - void findPolicyCycle() {
5.509 - for (int i = 0; i < int(_nodes->size()); ++i) {
5.510 - _level[(*_nodes)[i]] = -1;
5.511 - }
5.512 - LargeValue clength;
5.513 - int csize;
5.514 - Node u, v;
5.515 - _curr_found = false;
5.516 - for (int i = 0; i < int(_nodes->size()); ++i) {
5.517 - u = (*_nodes)[i];
5.518 - if (_level[u] >= 0) continue;
5.519 - for (; _level[u] < 0; u = _gr.target(_policy[u])) {
5.520 - _level[u] = i;
5.521 - }
5.522 - if (_level[u] == i) {
5.523 - // A cycle is found
5.524 - clength = _length[_policy[u]];
5.525 - csize = 1;
5.526 - for (v = u; (v = _gr.target(_policy[v])) != u; ) {
5.527 - clength += _length[_policy[v]];
5.528 - ++csize;
5.529 - }
5.530 - if ( !_curr_found ||
5.531 - (clength * _curr_size < _curr_length * csize) ) {
5.532 - _curr_found = true;
5.533 - _curr_length = clength;
5.534 - _curr_size = csize;
5.535 - _curr_node = u;
5.536 - }
5.537 - }
5.538 - }
5.539 - }
5.540 -
5.541 - // Contract the policy graph and compute node distances
5.542 - bool computeNodeDistances() {
5.543 - // Find the component of the main cycle and compute node distances
5.544 - // using reverse BFS
5.545 - for (int i = 0; i < int(_nodes->size()); ++i) {
5.546 - _reached[(*_nodes)[i]] = false;
5.547 - }
5.548 - _qfront = _qback = 0;
5.549 - _queue[0] = _curr_node;
5.550 - _reached[_curr_node] = true;
5.551 - _dist[_curr_node] = 0;
5.552 - Node u, v;
5.553 - Arc e;
5.554 - while (_qfront <= _qback) {
5.555 - v = _queue[_qfront++];
5.556 - for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
5.557 - e = _in_arcs[v][j];
5.558 - u = _gr.source(e);
5.559 - if (_policy[u] == e && !_reached[u]) {
5.560 - _reached[u] = true;
5.561 - _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
5.562 - _queue[++_qback] = u;
5.563 - }
5.564 - }
5.565 - }
5.566 -
5.567 - // Connect all other nodes to this component and compute node
5.568 - // distances using reverse BFS
5.569 - _qfront = 0;
5.570 - while (_qback < int(_nodes->size())-1) {
5.571 - v = _queue[_qfront++];
5.572 - for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
5.573 - e = _in_arcs[v][j];
5.574 - u = _gr.source(e);
5.575 - if (!_reached[u]) {
5.576 - _reached[u] = true;
5.577 - _policy[u] = e;
5.578 - _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
5.579 - _queue[++_qback] = u;
5.580 - }
5.581 - }
5.582 - }
5.583 -
5.584 - // Improve node distances
5.585 - bool improved = false;
5.586 - for (int i = 0; i < int(_nodes->size()); ++i) {
5.587 - v = (*_nodes)[i];
5.588 - for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
5.589 - e = _in_arcs[v][j];
5.590 - u = _gr.source(e);
5.591 - LargeValue delta = _dist[v] + _length[e] * _curr_size - _curr_length;
5.592 - if (_tolerance.less(delta, _dist[u])) {
5.593 - _dist[u] = delta;
5.594 - _policy[u] = e;
5.595 - improved = true;
5.596 - }
5.597 - }
5.598 - }
5.599 - return improved;
5.600 - }
5.601 -
5.602 - }; //class Howard
5.603 -
5.604 - ///@}
5.605 -
5.606 -} //namespace lemon
5.607 -
5.608 -#endif //LEMON_HOWARD_H
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/lemon/howard_mmc.h Sat Mar 13 22:01:38 2010 +0100
6.3 @@ -0,0 +1,605 @@
6.4 +/* -*- C++ -*-
6.5 + *
6.6 + * This file is a part of LEMON, a generic C++ optimization library
6.7 + *
6.8 + * Copyright (C) 2003-2008
6.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
6.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
6.11 + *
6.12 + * Permission to use, modify and distribute this software is granted
6.13 + * provided that this copyright notice appears in all copies. For
6.14 + * precise terms see the accompanying LICENSE file.
6.15 + *
6.16 + * This software is provided "AS IS" with no warranty of any kind,
6.17 + * express or implied, and with no claim as to its suitability for any
6.18 + * purpose.
6.19 + *
6.20 + */
6.21 +
6.22 +#ifndef LEMON_HOWARD_MMC_H
6.23 +#define LEMON_HOWARD_MMC_H
6.24 +
6.25 +/// \ingroup min_mean_cycle
6.26 +///
6.27 +/// \file
6.28 +/// \brief Howard's algorithm for finding a minimum mean cycle.
6.29 +
6.30 +#include <vector>
6.31 +#include <limits>
6.32 +#include <lemon/core.h>
6.33 +#include <lemon/path.h>
6.34 +#include <lemon/tolerance.h>
6.35 +#include <lemon/connectivity.h>
6.36 +
6.37 +namespace lemon {
6.38 +
6.39 + /// \brief Default traits class of HowardMmc class.
6.40 + ///
6.41 + /// Default traits class of HowardMmc class.
6.42 + /// \tparam GR The type of the digraph.
6.43 + /// \tparam CM The type of the cost map.
6.44 + /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
6.45 +#ifdef DOXYGEN
6.46 + template <typename GR, typename CM>
6.47 +#else
6.48 + template <typename GR, typename CM,
6.49 + bool integer = std::numeric_limits<typename CM::Value>::is_integer>
6.50 +#endif
6.51 + struct HowardMmcDefaultTraits
6.52 + {
6.53 + /// The type of the digraph
6.54 + typedef GR Digraph;
6.55 + /// The type of the cost map
6.56 + typedef CM CostMap;
6.57 + /// The type of the arc costs
6.58 + typedef typename CostMap::Value Cost;
6.59 +
6.60 + /// \brief The large cost type used for internal computations
6.61 + ///
6.62 + /// The large cost type used for internal computations.
6.63 + /// It is \c long \c long if the \c Cost type is integer,
6.64 + /// otherwise it is \c double.
6.65 + /// \c Cost must be convertible to \c LargeCost.
6.66 + typedef double LargeCost;
6.67 +
6.68 + /// The tolerance type used for internal computations
6.69 + typedef lemon::Tolerance<LargeCost> Tolerance;
6.70 +
6.71 + /// \brief The path type of the found cycles
6.72 + ///
6.73 + /// The path type of the found cycles.
6.74 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
6.75 + /// and it must have an \c addBack() function.
6.76 + typedef lemon::Path<Digraph> Path;
6.77 + };
6.78 +
6.79 + // Default traits class for integer cost types
6.80 + template <typename GR, typename CM>
6.81 + struct HowardMmcDefaultTraits<GR, CM, true>
6.82 + {
6.83 + typedef GR Digraph;
6.84 + typedef CM CostMap;
6.85 + typedef typename CostMap::Value Cost;
6.86 +#ifdef LEMON_HAVE_LONG_LONG
6.87 + typedef long long LargeCost;
6.88 +#else
6.89 + typedef long LargeCost;
6.90 +#endif
6.91 + typedef lemon::Tolerance<LargeCost> Tolerance;
6.92 + typedef lemon::Path<Digraph> Path;
6.93 + };
6.94 +
6.95 +
6.96 + /// \addtogroup min_mean_cycle
6.97 + /// @{
6.98 +
6.99 + /// \brief Implementation of Howard's algorithm for finding a minimum
6.100 + /// mean cycle.
6.101 + ///
6.102 + /// This class implements Howard's policy iteration algorithm for finding
6.103 + /// a directed cycle of minimum mean cost in a digraph
6.104 + /// \ref amo93networkflows, \ref dasdan98minmeancycle.
6.105 + /// This class provides the most efficient algorithm for the
6.106 + /// minimum mean cycle problem, though the best known theoretical
6.107 + /// bound on its running time is exponential.
6.108 + ///
6.109 + /// \tparam GR The type of the digraph the algorithm runs on.
6.110 + /// \tparam CM The type of the cost map. The default
6.111 + /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
6.112 + /// \tparam TR The traits class that defines various types used by the
6.113 + /// algorithm. By default, it is \ref HowardMmcDefaultTraits
6.114 + /// "HowardMmcDefaultTraits<GR, CM>".
6.115 + /// In most cases, this parameter should not be set directly,
6.116 + /// consider to use the named template parameters instead.
6.117 +#ifdef DOXYGEN
6.118 + template <typename GR, typename CM, typename TR>
6.119 +#else
6.120 + template < typename GR,
6.121 + typename CM = typename GR::template ArcMap<int>,
6.122 + typename TR = HowardMmcDefaultTraits<GR, CM> >
6.123 +#endif
6.124 + class HowardMmc
6.125 + {
6.126 + public:
6.127 +
6.128 + /// The type of the digraph
6.129 + typedef typename TR::Digraph Digraph;
6.130 + /// The type of the cost map
6.131 + typedef typename TR::CostMap CostMap;
6.132 + /// The type of the arc costs
6.133 + typedef typename TR::Cost Cost;
6.134 +
6.135 + /// \brief The large cost type
6.136 + ///
6.137 + /// The large cost type used for internal computations.
6.138 + /// By default, it is \c long \c long if the \c Cost type is integer,
6.139 + /// otherwise it is \c double.
6.140 + typedef typename TR::LargeCost LargeCost;
6.141 +
6.142 + /// The tolerance type
6.143 + typedef typename TR::Tolerance Tolerance;
6.144 +
6.145 + /// \brief The path type of the found cycles
6.146 + ///
6.147 + /// The path type of the found cycles.
6.148 + /// Using the \ref HowardMmcDefaultTraits "default traits class",
6.149 + /// it is \ref lemon::Path "Path<Digraph>".
6.150 + typedef typename TR::Path Path;
6.151 +
6.152 + /// The \ref HowardMmcDefaultTraits "traits class" of the algorithm
6.153 + typedef TR Traits;
6.154 +
6.155 + private:
6.156 +
6.157 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
6.158 +
6.159 + // The digraph the algorithm runs on
6.160 + const Digraph &_gr;
6.161 + // The cost of the arcs
6.162 + const CostMap &_cost;
6.163 +
6.164 + // Data for the found cycles
6.165 + bool _curr_found, _best_found;
6.166 + LargeCost _curr_cost, _best_cost;
6.167 + int _curr_size, _best_size;
6.168 + Node _curr_node, _best_node;
6.169 +
6.170 + Path *_cycle_path;
6.171 + bool _local_path;
6.172 +
6.173 + // Internal data used by the algorithm
6.174 + typename Digraph::template NodeMap<Arc> _policy;
6.175 + typename Digraph::template NodeMap<bool> _reached;
6.176 + typename Digraph::template NodeMap<int> _level;
6.177 + typename Digraph::template NodeMap<LargeCost> _dist;
6.178 +
6.179 + // Data for storing the strongly connected components
6.180 + int _comp_num;
6.181 + typename Digraph::template NodeMap<int> _comp;
6.182 + std::vector<std::vector<Node> > _comp_nodes;
6.183 + std::vector<Node>* _nodes;
6.184 + typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;
6.185 +
6.186 + // Queue used for BFS search
6.187 + std::vector<Node> _queue;
6.188 + int _qfront, _qback;
6.189 +
6.190 + Tolerance _tolerance;
6.191 +
6.192 + // Infinite constant
6.193 + const LargeCost INF;
6.194 +
6.195 + public:
6.196 +
6.197 + /// \name Named Template Parameters
6.198 + /// @{
6.199 +
6.200 + template <typename T>
6.201 + struct SetLargeCostTraits : public Traits {
6.202 + typedef T LargeCost;
6.203 + typedef lemon::Tolerance<T> Tolerance;
6.204 + };
6.205 +
6.206 + /// \brief \ref named-templ-param "Named parameter" for setting
6.207 + /// \c LargeCost type.
6.208 + ///
6.209 + /// \ref named-templ-param "Named parameter" for setting \c LargeCost
6.210 + /// type. It is used for internal computations in the algorithm.
6.211 + template <typename T>
6.212 + struct SetLargeCost
6.213 + : public HowardMmc<GR, CM, SetLargeCostTraits<T> > {
6.214 + typedef HowardMmc<GR, CM, SetLargeCostTraits<T> > Create;
6.215 + };
6.216 +
6.217 + template <typename T>
6.218 + struct SetPathTraits : public Traits {
6.219 + typedef T Path;
6.220 + };
6.221 +
6.222 + /// \brief \ref named-templ-param "Named parameter" for setting
6.223 + /// \c %Path type.
6.224 + ///
6.225 + /// \ref named-templ-param "Named parameter" for setting the \c %Path
6.226 + /// type of the found cycles.
6.227 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
6.228 + /// and it must have an \c addBack() function.
6.229 + template <typename T>
6.230 + struct SetPath
6.231 + : public HowardMmc<GR, CM, SetPathTraits<T> > {
6.232 + typedef HowardMmc<GR, CM, SetPathTraits<T> > Create;
6.233 + };
6.234 +
6.235 + /// @}
6.236 +
6.237 + protected:
6.238 +
6.239 + HowardMmc() {}
6.240 +
6.241 + public:
6.242 +
6.243 + /// \brief Constructor.
6.244 + ///
6.245 + /// The constructor of the class.
6.246 + ///
6.247 + /// \param digraph The digraph the algorithm runs on.
6.248 + /// \param cost The costs of the arcs.
6.249 + HowardMmc( const Digraph &digraph,
6.250 + const CostMap &cost ) :
6.251 + _gr(digraph), _cost(cost), _best_found(false),
6.252 + _best_cost(0), _best_size(1), _cycle_path(NULL), _local_path(false),
6.253 + _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),
6.254 + _comp(digraph), _in_arcs(digraph),
6.255 + INF(std::numeric_limits<LargeCost>::has_infinity ?
6.256 + std::numeric_limits<LargeCost>::infinity() :
6.257 + std::numeric_limits<LargeCost>::max())
6.258 + {}
6.259 +
6.260 + /// Destructor.
6.261 + ~HowardMmc() {
6.262 + if (_local_path) delete _cycle_path;
6.263 + }
6.264 +
6.265 + /// \brief Set the path structure for storing the found cycle.
6.266 + ///
6.267 + /// This function sets an external path structure for storing the
6.268 + /// found cycle.
6.269 + ///
6.270 + /// If you don't call this function before calling \ref run() or
6.271 + /// \ref findCycleMean(), it will allocate a local \ref Path "path"
6.272 + /// structure. The destuctor deallocates this automatically
6.273 + /// allocated object, of course.
6.274 + ///
6.275 + /// \note The algorithm calls only the \ref lemon::Path::addBack()
6.276 + /// "addBack()" function of the given path structure.
6.277 + ///
6.278 + /// \return <tt>(*this)</tt>
6.279 + HowardMmc& cycle(Path &path) {
6.280 + if (_local_path) {
6.281 + delete _cycle_path;
6.282 + _local_path = false;
6.283 + }
6.284 + _cycle_path = &path;
6.285 + return *this;
6.286 + }
6.287 +
6.288 + /// \brief Set the tolerance used by the algorithm.
6.289 + ///
6.290 + /// This function sets the tolerance object used by the algorithm.
6.291 + ///
6.292 + /// \return <tt>(*this)</tt>
6.293 + HowardMmc& tolerance(const Tolerance& tolerance) {
6.294 + _tolerance = tolerance;
6.295 + return *this;
6.296 + }
6.297 +
6.298 + /// \brief Return a const reference to the tolerance.
6.299 + ///
6.300 + /// This function returns a const reference to the tolerance object
6.301 + /// used by the algorithm.
6.302 + const Tolerance& tolerance() const {
6.303 + return _tolerance;
6.304 + }
6.305 +
6.306 + /// \name Execution control
6.307 + /// The simplest way to execute the algorithm is to call the \ref run()
6.308 + /// function.\n
6.309 + /// If you only need the minimum mean cost, you may call
6.310 + /// \ref findCycleMean().
6.311 +
6.312 + /// @{
6.313 +
6.314 + /// \brief Run the algorithm.
6.315 + ///
6.316 + /// This function runs the algorithm.
6.317 + /// It can be called more than once (e.g. if the underlying digraph
6.318 + /// and/or the arc costs have been modified).
6.319 + ///
6.320 + /// \return \c true if a directed cycle exists in the digraph.
6.321 + ///
6.322 + /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
6.323 + /// \code
6.324 + /// return mmc.findCycleMean() && mmc.findCycle();
6.325 + /// \endcode
6.326 + bool run() {
6.327 + return findCycleMean() && findCycle();
6.328 + }
6.329 +
6.330 + /// \brief Find the minimum cycle mean.
6.331 + ///
6.332 + /// This function finds the minimum mean cost of the directed
6.333 + /// cycles in the digraph.
6.334 + ///
6.335 + /// \return \c true if a directed cycle exists in the digraph.
6.336 + bool findCycleMean() {
6.337 + // Initialize and find strongly connected components
6.338 + init();
6.339 + findComponents();
6.340 +
6.341 + // Find the minimum cycle mean in the components
6.342 + for (int comp = 0; comp < _comp_num; ++comp) {
6.343 + // Find the minimum mean cycle in the current component
6.344 + if (!buildPolicyGraph(comp)) continue;
6.345 + while (true) {
6.346 + findPolicyCycle();
6.347 + if (!computeNodeDistances()) break;
6.348 + }
6.349 + // Update the best cycle (global minimum mean cycle)
6.350 + if ( _curr_found && (!_best_found ||
6.351 + _curr_cost * _best_size < _best_cost * _curr_size) ) {
6.352 + _best_found = true;
6.353 + _best_cost = _curr_cost;
6.354 + _best_size = _curr_size;
6.355 + _best_node = _curr_node;
6.356 + }
6.357 + }
6.358 + return _best_found;
6.359 + }
6.360 +
6.361 + /// \brief Find a minimum mean directed cycle.
6.362 + ///
6.363 + /// This function finds a directed cycle of minimum mean cost
6.364 + /// in the digraph using the data computed by findCycleMean().
6.365 + ///
6.366 + /// \return \c true if a directed cycle exists in the digraph.
6.367 + ///
6.368 + /// \pre \ref findCycleMean() must be called before using this function.
6.369 + bool findCycle() {
6.370 + if (!_best_found) return false;
6.371 + _cycle_path->addBack(_policy[_best_node]);
6.372 + for ( Node v = _best_node;
6.373 + (v = _gr.target(_policy[v])) != _best_node; ) {
6.374 + _cycle_path->addBack(_policy[v]);
6.375 + }
6.376 + return true;
6.377 + }
6.378 +
6.379 + /// @}
6.380 +
6.381 + /// \name Query Functions
6.382 + /// The results of the algorithm can be obtained using these
6.383 + /// functions.\n
6.384 + /// The algorithm should be executed before using them.
6.385 +
6.386 + /// @{
6.387 +
6.388 + /// \brief Return the total cost of the found cycle.
6.389 + ///
6.390 + /// This function returns the total cost of the found cycle.
6.391 + ///
6.392 + /// \pre \ref run() or \ref findCycleMean() must be called before
6.393 + /// using this function.
6.394 + Cost cycleCost() const {
6.395 + return static_cast<Cost>(_best_cost);
6.396 + }
6.397 +
6.398 + /// \brief Return the number of arcs on the found cycle.
6.399 + ///
6.400 + /// This function returns the number of arcs on the found cycle.
6.401 + ///
6.402 + /// \pre \ref run() or \ref findCycleMean() must be called before
6.403 + /// using this function.
6.404 + int cycleSize() const {
6.405 + return _best_size;
6.406 + }
6.407 +
6.408 + /// \brief Return the mean cost of the found cycle.
6.409 + ///
6.410 + /// This function returns the mean cost of the found cycle.
6.411 + ///
6.412 + /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
6.413 + /// following code.
6.414 + /// \code
6.415 + /// return static_cast<double>(alg.cycleCost()) / alg.cycleSize();
6.416 + /// \endcode
6.417 + ///
6.418 + /// \pre \ref run() or \ref findCycleMean() must be called before
6.419 + /// using this function.
6.420 + double cycleMean() const {
6.421 + return static_cast<double>(_best_cost) / _best_size;
6.422 + }
6.423 +
6.424 + /// \brief Return the found cycle.
6.425 + ///
6.426 + /// This function returns a const reference to the path structure
6.427 + /// storing the found cycle.
6.428 + ///
6.429 + /// \pre \ref run() or \ref findCycle() must be called before using
6.430 + /// this function.
6.431 + const Path& cycle() const {
6.432 + return *_cycle_path;
6.433 + }
6.434 +
6.435 + ///@}
6.436 +
6.437 + private:
6.438 +
6.439 + // Initialize
6.440 + void init() {
6.441 + if (!_cycle_path) {
6.442 + _local_path = true;
6.443 + _cycle_path = new Path;
6.444 + }
6.445 + _queue.resize(countNodes(_gr));
6.446 + _best_found = false;
6.447 + _best_cost = 0;
6.448 + _best_size = 1;
6.449 + _cycle_path->clear();
6.450 + }
6.451 +
6.452 + // Find strongly connected components and initialize _comp_nodes
6.453 + // and _in_arcs
6.454 + void findComponents() {
6.455 + _comp_num = stronglyConnectedComponents(_gr, _comp);
6.456 + _comp_nodes.resize(_comp_num);
6.457 + if (_comp_num == 1) {
6.458 + _comp_nodes[0].clear();
6.459 + for (NodeIt n(_gr); n != INVALID; ++n) {
6.460 + _comp_nodes[0].push_back(n);
6.461 + _in_arcs[n].clear();
6.462 + for (InArcIt a(_gr, n); a != INVALID; ++a) {
6.463 + _in_arcs[n].push_back(a);
6.464 + }
6.465 + }
6.466 + } else {
6.467 + for (int i = 0; i < _comp_num; ++i)
6.468 + _comp_nodes[i].clear();
6.469 + for (NodeIt n(_gr); n != INVALID; ++n) {
6.470 + int k = _comp[n];
6.471 + _comp_nodes[k].push_back(n);
6.472 + _in_arcs[n].clear();
6.473 + for (InArcIt a(_gr, n); a != INVALID; ++a) {
6.474 + if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);
6.475 + }
6.476 + }
6.477 + }
6.478 + }
6.479 +
6.480 + // Build the policy graph in the given strongly connected component
6.481 + // (the out-degree of every node is 1)
6.482 + bool buildPolicyGraph(int comp) {
6.483 + _nodes = &(_comp_nodes[comp]);
6.484 + if (_nodes->size() < 1 ||
6.485 + (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {
6.486 + return false;
6.487 + }
6.488 + for (int i = 0; i < int(_nodes->size()); ++i) {
6.489 + _dist[(*_nodes)[i]] = INF;
6.490 + }
6.491 + Node u, v;
6.492 + Arc e;
6.493 + for (int i = 0; i < int(_nodes->size()); ++i) {
6.494 + v = (*_nodes)[i];
6.495 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
6.496 + e = _in_arcs[v][j];
6.497 + u = _gr.source(e);
6.498 + if (_cost[e] < _dist[u]) {
6.499 + _dist[u] = _cost[e];
6.500 + _policy[u] = e;
6.501 + }
6.502 + }
6.503 + }
6.504 + return true;
6.505 + }
6.506 +
6.507 + // Find the minimum mean cycle in the policy graph
6.508 + void findPolicyCycle() {
6.509 + for (int i = 0; i < int(_nodes->size()); ++i) {
6.510 + _level[(*_nodes)[i]] = -1;
6.511 + }
6.512 + LargeCost ccost;
6.513 + int csize;
6.514 + Node u, v;
6.515 + _curr_found = false;
6.516 + for (int i = 0; i < int(_nodes->size()); ++i) {
6.517 + u = (*_nodes)[i];
6.518 + if (_level[u] >= 0) continue;
6.519 + for (; _level[u] < 0; u = _gr.target(_policy[u])) {
6.520 + _level[u] = i;
6.521 + }
6.522 + if (_level[u] == i) {
6.523 + // A cycle is found
6.524 + ccost = _cost[_policy[u]];
6.525 + csize = 1;
6.526 + for (v = u; (v = _gr.target(_policy[v])) != u; ) {
6.527 + ccost += _cost[_policy[v]];
6.528 + ++csize;
6.529 + }
6.530 + if ( !_curr_found ||
6.531 + (ccost * _curr_size < _curr_cost * csize) ) {
6.532 + _curr_found = true;
6.533 + _curr_cost = ccost;
6.534 + _curr_size = csize;
6.535 + _curr_node = u;
6.536 + }
6.537 + }
6.538 + }
6.539 + }
6.540 +
6.541 + // Contract the policy graph and compute node distances
6.542 + bool computeNodeDistances() {
6.543 + // Find the component of the main cycle and compute node distances
6.544 + // using reverse BFS
6.545 + for (int i = 0; i < int(_nodes->size()); ++i) {
6.546 + _reached[(*_nodes)[i]] = false;
6.547 + }
6.548 + _qfront = _qback = 0;
6.549 + _queue[0] = _curr_node;
6.550 + _reached[_curr_node] = true;
6.551 + _dist[_curr_node] = 0;
6.552 + Node u, v;
6.553 + Arc e;
6.554 + while (_qfront <= _qback) {
6.555 + v = _queue[_qfront++];
6.556 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
6.557 + e = _in_arcs[v][j];
6.558 + u = _gr.source(e);
6.559 + if (_policy[u] == e && !_reached[u]) {
6.560 + _reached[u] = true;
6.561 + _dist[u] = _dist[v] + _cost[e] * _curr_size - _curr_cost;
6.562 + _queue[++_qback] = u;
6.563 + }
6.564 + }
6.565 + }
6.566 +
6.567 + // Connect all other nodes to this component and compute node
6.568 + // distances using reverse BFS
6.569 + _qfront = 0;
6.570 + while (_qback < int(_nodes->size())-1) {
6.571 + v = _queue[_qfront++];
6.572 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
6.573 + e = _in_arcs[v][j];
6.574 + u = _gr.source(e);
6.575 + if (!_reached[u]) {
6.576 + _reached[u] = true;
6.577 + _policy[u] = e;
6.578 + _dist[u] = _dist[v] + _cost[e] * _curr_size - _curr_cost;
6.579 + _queue[++_qback] = u;
6.580 + }
6.581 + }
6.582 + }
6.583 +
6.584 + // Improve node distances
6.585 + bool improved = false;
6.586 + for (int i = 0; i < int(_nodes->size()); ++i) {
6.587 + v = (*_nodes)[i];
6.588 + for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
6.589 + e = _in_arcs[v][j];
6.590 + u = _gr.source(e);
6.591 + LargeCost delta = _dist[v] + _cost[e] * _curr_size - _curr_cost;
6.592 + if (_tolerance.less(delta, _dist[u])) {
6.593 + _dist[u] = delta;
6.594 + _policy[u] = e;
6.595 + improved = true;
6.596 + }
6.597 + }
6.598 + }
6.599 + return improved;
6.600 + }
6.601 +
6.602 + }; //class HowardMmc
6.603 +
6.604 + ///@}
6.605 +
6.606 +} //namespace lemon
6.607 +
6.608 +#endif //LEMON_HOWARD_MMC_H
7.1 --- a/lemon/karp.h Mon Mar 08 08:33:41 2010 +0100
7.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
7.3 @@ -1,590 +0,0 @@
7.4 -/* -*- C++ -*-
7.5 - *
7.6 - * This file is a part of LEMON, a generic C++ optimization library
7.7 - *
7.8 - * Copyright (C) 2003-2008
7.9 - * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7.10 - * (Egervary Research Group on Combinatorial Optimization, EGRES).
7.11 - *
7.12 - * Permission to use, modify and distribute this software is granted
7.13 - * provided that this copyright notice appears in all copies. For
7.14 - * precise terms see the accompanying LICENSE file.
7.15 - *
7.16 - * This software is provided "AS IS" with no warranty of any kind,
7.17 - * express or implied, and with no claim as to its suitability for any
7.18 - * purpose.
7.19 - *
7.20 - */
7.21 -
7.22 -#ifndef LEMON_KARP_H
7.23 -#define LEMON_KARP_H
7.24 -
7.25 -/// \ingroup min_mean_cycle
7.26 -///
7.27 -/// \file
7.28 -/// \brief Karp's algorithm for finding a minimum mean cycle.
7.29 -
7.30 -#include <vector>
7.31 -#include <limits>
7.32 -#include <lemon/core.h>
7.33 -#include <lemon/path.h>
7.34 -#include <lemon/tolerance.h>
7.35 -#include <lemon/connectivity.h>
7.36 -
7.37 -namespace lemon {
7.38 -
7.39 - /// \brief Default traits class of Karp algorithm.
7.40 - ///
7.41 - /// Default traits class of Karp algorithm.
7.42 - /// \tparam GR The type of the digraph.
7.43 - /// \tparam LEN The type of the length map.
7.44 - /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
7.45 -#ifdef DOXYGEN
7.46 - template <typename GR, typename LEN>
7.47 -#else
7.48 - template <typename GR, typename LEN,
7.49 - bool integer = std::numeric_limits<typename LEN::Value>::is_integer>
7.50 -#endif
7.51 - struct KarpDefaultTraits
7.52 - {
7.53 - /// The type of the digraph
7.54 - typedef GR Digraph;
7.55 - /// The type of the length map
7.56 - typedef LEN LengthMap;
7.57 - /// The type of the arc lengths
7.58 - typedef typename LengthMap::Value Value;
7.59 -
7.60 - /// \brief The large value type used for internal computations
7.61 - ///
7.62 - /// The large value type used for internal computations.
7.63 - /// It is \c long \c long if the \c Value type is integer,
7.64 - /// otherwise it is \c double.
7.65 - /// \c Value must be convertible to \c LargeValue.
7.66 - typedef double LargeValue;
7.67 -
7.68 - /// The tolerance type used for internal computations
7.69 - typedef lemon::Tolerance<LargeValue> Tolerance;
7.70 -
7.71 - /// \brief The path type of the found cycles
7.72 - ///
7.73 - /// The path type of the found cycles.
7.74 - /// It must conform to the \ref lemon::concepts::Path "Path" concept
7.75 - /// and it must have an \c addFront() function.
7.76 - typedef lemon::Path<Digraph> Path;
7.77 - };
7.78 -
7.79 - // Default traits class for integer value types
7.80 - template <typename GR, typename LEN>
7.81 - struct KarpDefaultTraits<GR, LEN, true>
7.82 - {
7.83 - typedef GR Digraph;
7.84 - typedef LEN LengthMap;
7.85 - typedef typename LengthMap::Value Value;
7.86 -#ifdef LEMON_HAVE_LONG_LONG
7.87 - typedef long long LargeValue;
7.88 -#else
7.89 - typedef long LargeValue;
7.90 -#endif
7.91 - typedef lemon::Tolerance<LargeValue> Tolerance;
7.92 - typedef lemon::Path<Digraph> Path;
7.93 - };
7.94 -
7.95 -
7.96 - /// \addtogroup min_mean_cycle
7.97 - /// @{
7.98 -
7.99 - /// \brief Implementation of Karp's algorithm for finding a minimum
7.100 - /// mean cycle.
7.101 - ///
7.102 - /// This class implements Karp's algorithm for finding a directed
7.103 - /// cycle of minimum mean length (cost) in a digraph
7.104 - /// \ref amo93networkflows, \ref dasdan98minmeancycle.
7.105 - /// It runs in time O(ne) and uses space O(n<sup>2</sup>+e).
7.106 - ///
7.107 - /// \tparam GR The type of the digraph the algorithm runs on.
7.108 - /// \tparam LEN The type of the length map. The default
7.109 - /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
7.110 - /// \tparam TR The traits class that defines various types used by the
7.111 - /// algorithm. By default, it is \ref KarpDefaultTraits
7.112 - /// "KarpDefaultTraits<GR, LEN>".
7.113 - /// In most cases, this parameter should not be set directly,
7.114 - /// consider to use the named template parameters instead.
7.115 -#ifdef DOXYGEN
7.116 - template <typename GR, typename LEN, typename TR>
7.117 -#else
7.118 - template < typename GR,
7.119 - typename LEN = typename GR::template ArcMap<int>,
7.120 - typename TR = KarpDefaultTraits<GR, LEN> >
7.121 -#endif
7.122 - class Karp
7.123 - {
7.124 - public:
7.125 -
7.126 - /// The type of the digraph
7.127 - typedef typename TR::Digraph Digraph;
7.128 - /// The type of the length map
7.129 - typedef typename TR::LengthMap LengthMap;
7.130 - /// The type of the arc lengths
7.131 - typedef typename TR::Value Value;
7.132 -
7.133 - /// \brief The large value type
7.134 - ///
7.135 - /// The large value type used for internal computations.
7.136 - /// By default, it is \c long \c long if the \c Value type is integer,
7.137 - /// otherwise it is \c double.
7.138 - typedef typename TR::LargeValue LargeValue;
7.139 -
7.140 - /// The tolerance type
7.141 - typedef typename TR::Tolerance Tolerance;
7.142 -
7.143 - /// \brief The path type of the found cycles
7.144 - ///
7.145 - /// The path type of the found cycles.
7.146 - /// Using the \ref KarpDefaultTraits "default traits class",
7.147 - /// it is \ref lemon::Path "Path<Digraph>".
7.148 - typedef typename TR::Path Path;
7.149 -
7.150 - /// The \ref KarpDefaultTraits "traits class" of the algorithm
7.151 - typedef TR Traits;
7.152 -
7.153 - private:
7.154 -
7.155 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
7.156 -
7.157 - // Data sturcture for path data
7.158 - struct PathData
7.159 - {
7.160 - LargeValue dist;
7.161 - Arc pred;
7.162 - PathData(LargeValue d, Arc p = INVALID) :
7.163 - dist(d), pred(p) {}
7.164 - };
7.165 -
7.166 - typedef typename Digraph::template NodeMap<std::vector<PathData> >
7.167 - PathDataNodeMap;
7.168 -
7.169 - private:
7.170 -
7.171 - // The digraph the algorithm runs on
7.172 - const Digraph &_gr;
7.173 - // The length of the arcs
7.174 - const LengthMap &_length;
7.175 -
7.176 - // Data for storing the strongly connected components
7.177 - int _comp_num;
7.178 - typename Digraph::template NodeMap<int> _comp;
7.179 - std::vector<std::vector<Node> > _comp_nodes;
7.180 - std::vector<Node>* _nodes;
7.181 - typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
7.182 -
7.183 - // Data for the found cycle
7.184 - LargeValue _cycle_length;
7.185 - int _cycle_size;
7.186 - Node _cycle_node;
7.187 -
7.188 - Path *_cycle_path;
7.189 - bool _local_path;
7.190 -
7.191 - // Node map for storing path data
7.192 - PathDataNodeMap _data;
7.193 - // The processed nodes in the last round
7.194 - std::vector<Node> _process;
7.195 -
7.196 - Tolerance _tolerance;
7.197 -
7.198 - // Infinite constant
7.199 - const LargeValue INF;
7.200 -
7.201 - public:
7.202 -
7.203 - /// \name Named Template Parameters
7.204 - /// @{
7.205 -
7.206 - template <typename T>
7.207 - struct SetLargeValueTraits : public Traits {
7.208 - typedef T LargeValue;
7.209 - typedef lemon::Tolerance<T> Tolerance;
7.210 - };
7.211 -
7.212 - /// \brief \ref named-templ-param "Named parameter" for setting
7.213 - /// \c LargeValue type.
7.214 - ///
7.215 - /// \ref named-templ-param "Named parameter" for setting \c LargeValue
7.216 - /// type. It is used for internal computations in the algorithm.
7.217 - template <typename T>
7.218 - struct SetLargeValue
7.219 - : public Karp<GR, LEN, SetLargeValueTraits<T> > {
7.220 - typedef Karp<GR, LEN, SetLargeValueTraits<T> > Create;
7.221 - };
7.222 -
7.223 - template <typename T>
7.224 - struct SetPathTraits : public Traits {
7.225 - typedef T Path;
7.226 - };
7.227 -
7.228 - /// \brief \ref named-templ-param "Named parameter" for setting
7.229 - /// \c %Path type.
7.230 - ///
7.231 - /// \ref named-templ-param "Named parameter" for setting the \c %Path
7.232 - /// type of the found cycles.
7.233 - /// It must conform to the \ref lemon::concepts::Path "Path" concept
7.234 - /// and it must have an \c addFront() function.
7.235 - template <typename T>
7.236 - struct SetPath
7.237 - : public Karp<GR, LEN, SetPathTraits<T> > {
7.238 - typedef Karp<GR, LEN, SetPathTraits<T> > Create;
7.239 - };
7.240 -
7.241 - /// @}
7.242 -
7.243 - protected:
7.244 -
7.245 - Karp() {}
7.246 -
7.247 - public:
7.248 -
7.249 - /// \brief Constructor.
7.250 - ///
7.251 - /// The constructor of the class.
7.252 - ///
7.253 - /// \param digraph The digraph the algorithm runs on.
7.254 - /// \param length The lengths (costs) of the arcs.
7.255 - Karp( const Digraph &digraph,
7.256 - const LengthMap &length ) :
7.257 - _gr(digraph), _length(length), _comp(digraph), _out_arcs(digraph),
7.258 - _cycle_length(0), _cycle_size(1), _cycle_node(INVALID),
7.259 - _cycle_path(NULL), _local_path(false), _data(digraph),
7.260 - INF(std::numeric_limits<LargeValue>::has_infinity ?
7.261 - std::numeric_limits<LargeValue>::infinity() :
7.262 - std::numeric_limits<LargeValue>::max())
7.263 - {}
7.264 -
7.265 - /// Destructor.
7.266 - ~Karp() {
7.267 - if (_local_path) delete _cycle_path;
7.268 - }
7.269 -
7.270 - /// \brief Set the path structure for storing the found cycle.
7.271 - ///
7.272 - /// This function sets an external path structure for storing the
7.273 - /// found cycle.
7.274 - ///
7.275 - /// If you don't call this function before calling \ref run() or
7.276 - /// \ref findMinMean(), it will allocate a local \ref Path "path"
7.277 - /// structure. The destuctor deallocates this automatically
7.278 - /// allocated object, of course.
7.279 - ///
7.280 - /// \note The algorithm calls only the \ref lemon::Path::addFront()
7.281 - /// "addFront()" function of the given path structure.
7.282 - ///
7.283 - /// \return <tt>(*this)</tt>
7.284 - Karp& cycle(Path &path) {
7.285 - if (_local_path) {
7.286 - delete _cycle_path;
7.287 - _local_path = false;
7.288 - }
7.289 - _cycle_path = &path;
7.290 - return *this;
7.291 - }
7.292 -
7.293 - /// \brief Set the tolerance used by the algorithm.
7.294 - ///
7.295 - /// This function sets the tolerance object used by the algorithm.
7.296 - ///
7.297 - /// \return <tt>(*this)</tt>
7.298 - Karp& tolerance(const Tolerance& tolerance) {
7.299 - _tolerance = tolerance;
7.300 - return *this;
7.301 - }
7.302 -
7.303 - /// \brief Return a const reference to the tolerance.
7.304 - ///
7.305 - /// This function returns a const reference to the tolerance object
7.306 - /// used by the algorithm.
7.307 - const Tolerance& tolerance() const {
7.308 - return _tolerance;
7.309 - }
7.310 -
7.311 - /// \name Execution control
7.312 - /// The simplest way to execute the algorithm is to call the \ref run()
7.313 - /// function.\n
7.314 - /// If you only need the minimum mean length, you may call
7.315 - /// \ref findMinMean().
7.316 -
7.317 - /// @{
7.318 -
7.319 - /// \brief Run the algorithm.
7.320 - ///
7.321 - /// This function runs the algorithm.
7.322 - /// It can be called more than once (e.g. if the underlying digraph
7.323 - /// and/or the arc lengths have been modified).
7.324 - ///
7.325 - /// \return \c true if a directed cycle exists in the digraph.
7.326 - ///
7.327 - /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
7.328 - /// \code
7.329 - /// return mmc.findMinMean() && mmc.findCycle();
7.330 - /// \endcode
7.331 - bool run() {
7.332 - return findMinMean() && findCycle();
7.333 - }
7.334 -
7.335 - /// \brief Find the minimum cycle mean.
7.336 - ///
7.337 - /// This function finds the minimum mean length of the directed
7.338 - /// cycles in the digraph.
7.339 - ///
7.340 - /// \return \c true if a directed cycle exists in the digraph.
7.341 - bool findMinMean() {
7.342 - // Initialization and find strongly connected components
7.343 - init();
7.344 - findComponents();
7.345 -
7.346 - // Find the minimum cycle mean in the components
7.347 - for (int comp = 0; comp < _comp_num; ++comp) {
7.348 - if (!initComponent(comp)) continue;
7.349 - processRounds();
7.350 - updateMinMean();
7.351 - }
7.352 - return (_cycle_node != INVALID);
7.353 - }
7.354 -
7.355 - /// \brief Find a minimum mean directed cycle.
7.356 - ///
7.357 - /// This function finds a directed cycle of minimum mean length
7.358 - /// in the digraph using the data computed by findMinMean().
7.359 - ///
7.360 - /// \return \c true if a directed cycle exists in the digraph.
7.361 - ///
7.362 - /// \pre \ref findMinMean() must be called before using this function.
7.363 - bool findCycle() {
7.364 - if (_cycle_node == INVALID) return false;
7.365 - IntNodeMap reached(_gr, -1);
7.366 - int r = _data[_cycle_node].size();
7.367 - Node u = _cycle_node;
7.368 - while (reached[u] < 0) {
7.369 - reached[u] = --r;
7.370 - u = _gr.source(_data[u][r].pred);
7.371 - }
7.372 - r = reached[u];
7.373 - Arc e = _data[u][r].pred;
7.374 - _cycle_path->addFront(e);
7.375 - _cycle_length = _length[e];
7.376 - _cycle_size = 1;
7.377 - Node v;
7.378 - while ((v = _gr.source(e)) != u) {
7.379 - e = _data[v][--r].pred;
7.380 - _cycle_path->addFront(e);
7.381 - _cycle_length += _length[e];
7.382 - ++_cycle_size;
7.383 - }
7.384 - return true;
7.385 - }
7.386 -
7.387 - /// @}
7.388 -
7.389 - /// \name Query Functions
7.390 - /// The results of the algorithm can be obtained using these
7.391 - /// functions.\n
7.392 - /// The algorithm should be executed before using them.
7.393 -
7.394 - /// @{
7.395 -
7.396 - /// \brief Return the total length of the found cycle.
7.397 - ///
7.398 - /// This function returns the total length of the found cycle.
7.399 - ///
7.400 - /// \pre \ref run() or \ref findMinMean() must be called before
7.401 - /// using this function.
7.402 - Value cycleLength() const {
7.403 - return static_cast<Value>(_cycle_length);
7.404 - }
7.405 -
7.406 - /// \brief Return the number of arcs on the found cycle.
7.407 - ///
7.408 - /// This function returns the number of arcs on the found cycle.
7.409 - ///
7.410 - /// \pre \ref run() or \ref findMinMean() must be called before
7.411 - /// using this function.
7.412 - int cycleArcNum() const {
7.413 - return _cycle_size;
7.414 - }
7.415 -
7.416 - /// \brief Return the mean length of the found cycle.
7.417 - ///
7.418 - /// This function returns the mean length of the found cycle.
7.419 - ///
7.420 - /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
7.421 - /// following code.
7.422 - /// \code
7.423 - /// return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
7.424 - /// \endcode
7.425 - ///
7.426 - /// \pre \ref run() or \ref findMinMean() must be called before
7.427 - /// using this function.
7.428 - double cycleMean() const {
7.429 - return static_cast<double>(_cycle_length) / _cycle_size;
7.430 - }
7.431 -
7.432 - /// \brief Return the found cycle.
7.433 - ///
7.434 - /// This function returns a const reference to the path structure
7.435 - /// storing the found cycle.
7.436 - ///
7.437 - /// \pre \ref run() or \ref findCycle() must be called before using
7.438 - /// this function.
7.439 - const Path& cycle() const {
7.440 - return *_cycle_path;
7.441 - }
7.442 -
7.443 - ///@}
7.444 -
7.445 - private:
7.446 -
7.447 - // Initialization
7.448 - void init() {
7.449 - if (!_cycle_path) {
7.450 - _local_path = true;
7.451 - _cycle_path = new Path;
7.452 - }
7.453 - _cycle_path->clear();
7.454 - _cycle_length = 0;
7.455 - _cycle_size = 1;
7.456 - _cycle_node = INVALID;
7.457 - for (NodeIt u(_gr); u != INVALID; ++u)
7.458 - _data[u].clear();
7.459 - }
7.460 -
7.461 - // Find strongly connected components and initialize _comp_nodes
7.462 - // and _out_arcs
7.463 - void findComponents() {
7.464 - _comp_num = stronglyConnectedComponents(_gr, _comp);
7.465 - _comp_nodes.resize(_comp_num);
7.466 - if (_comp_num == 1) {
7.467 - _comp_nodes[0].clear();
7.468 - for (NodeIt n(_gr); n != INVALID; ++n) {
7.469 - _comp_nodes[0].push_back(n);
7.470 - _out_arcs[n].clear();
7.471 - for (OutArcIt a(_gr, n); a != INVALID; ++a) {
7.472 - _out_arcs[n].push_back(a);
7.473 - }
7.474 - }
7.475 - } else {
7.476 - for (int i = 0; i < _comp_num; ++i)
7.477 - _comp_nodes[i].clear();
7.478 - for (NodeIt n(_gr); n != INVALID; ++n) {
7.479 - int k = _comp[n];
7.480 - _comp_nodes[k].push_back(n);
7.481 - _out_arcs[n].clear();
7.482 - for (OutArcIt a(_gr, n); a != INVALID; ++a) {
7.483 - if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
7.484 - }
7.485 - }
7.486 - }
7.487 - }
7.488 -
7.489 - // Initialize path data for the current component
7.490 - bool initComponent(int comp) {
7.491 - _nodes = &(_comp_nodes[comp]);
7.492 - int n = _nodes->size();
7.493 - if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
7.494 - return false;
7.495 - }
7.496 - for (int i = 0; i < n; ++i) {
7.497 - _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
7.498 - }
7.499 - return true;
7.500 - }
7.501 -
7.502 - // Process all rounds of computing path data for the current component.
7.503 - // _data[v][k] is the length of a shortest directed walk from the root
7.504 - // node to node v containing exactly k arcs.
7.505 - void processRounds() {
7.506 - Node start = (*_nodes)[0];
7.507 - _data[start][0] = PathData(0);
7.508 - _process.clear();
7.509 - _process.push_back(start);
7.510 -
7.511 - int k, n = _nodes->size();
7.512 - for (k = 1; k <= n && int(_process.size()) < n; ++k) {
7.513 - processNextBuildRound(k);
7.514 - }
7.515 - for ( ; k <= n; ++k) {
7.516 - processNextFullRound(k);
7.517 - }
7.518 - }
7.519 -
7.520 - // Process one round and rebuild _process
7.521 - void processNextBuildRound(int k) {
7.522 - std::vector<Node> next;
7.523 - Node u, v;
7.524 - Arc e;
7.525 - LargeValue d;
7.526 - for (int i = 0; i < int(_process.size()); ++i) {
7.527 - u = _process[i];
7.528 - for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
7.529 - e = _out_arcs[u][j];
7.530 - v = _gr.target(e);
7.531 - d = _data[u][k-1].dist + _length[e];
7.532 - if (_tolerance.less(d, _data[v][k].dist)) {
7.533 - if (_data[v][k].dist == INF) next.push_back(v);
7.534 - _data[v][k] = PathData(d, e);
7.535 - }
7.536 - }
7.537 - }
7.538 - _process.swap(next);
7.539 - }
7.540 -
7.541 - // Process one round using _nodes instead of _process
7.542 - void processNextFullRound(int k) {
7.543 - Node u, v;
7.544 - Arc e;
7.545 - LargeValue d;
7.546 - for (int i = 0; i < int(_nodes->size()); ++i) {
7.547 - u = (*_nodes)[i];
7.548 - for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
7.549 - e = _out_arcs[u][j];
7.550 - v = _gr.target(e);
7.551 - d = _data[u][k-1].dist + _length[e];
7.552 - if (_tolerance.less(d, _data[v][k].dist)) {
7.553 - _data[v][k] = PathData(d, e);
7.554 - }
7.555 - }
7.556 - }
7.557 - }
7.558 -
7.559 - // Update the minimum cycle mean
7.560 - void updateMinMean() {
7.561 - int n = _nodes->size();
7.562 - for (int i = 0; i < n; ++i) {
7.563 - Node u = (*_nodes)[i];
7.564 - if (_data[u][n].dist == INF) continue;
7.565 - LargeValue length, max_length = 0;
7.566 - int size, max_size = 1;
7.567 - bool found_curr = false;
7.568 - for (int k = 0; k < n; ++k) {
7.569 - if (_data[u][k].dist == INF) continue;
7.570 - length = _data[u][n].dist - _data[u][k].dist;
7.571 - size = n - k;
7.572 - if (!found_curr || length * max_size > max_length * size) {
7.573 - found_curr = true;
7.574 - max_length = length;
7.575 - max_size = size;
7.576 - }
7.577 - }
7.578 - if ( found_curr && (_cycle_node == INVALID ||
7.579 - max_length * _cycle_size < _cycle_length * max_size) ) {
7.580 - _cycle_length = max_length;
7.581 - _cycle_size = max_size;
7.582 - _cycle_node = u;
7.583 - }
7.584 - }
7.585 - }
7.586 -
7.587 - }; //class Karp
7.588 -
7.589 - ///@}
7.590 -
7.591 -} //namespace lemon
7.592 -
7.593 -#endif //LEMON_KARP_H
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
8.2 +++ b/lemon/karp_mmc.h Sat Mar 13 22:01:38 2010 +0100
8.3 @@ -0,0 +1,590 @@
8.4 +/* -*- C++ -*-
8.5 + *
8.6 + * This file is a part of LEMON, a generic C++ optimization library
8.7 + *
8.8 + * Copyright (C) 2003-2008
8.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
8.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
8.11 + *
8.12 + * Permission to use, modify and distribute this software is granted
8.13 + * provided that this copyright notice appears in all copies. For
8.14 + * precise terms see the accompanying LICENSE file.
8.15 + *
8.16 + * This software is provided "AS IS" with no warranty of any kind,
8.17 + * express or implied, and with no claim as to its suitability for any
8.18 + * purpose.
8.19 + *
8.20 + */
8.21 +
8.22 +#ifndef LEMON_KARP_MMC_H
8.23 +#define LEMON_KARP_MMC_H
8.24 +
8.25 +/// \ingroup min_mean_cycle
8.26 +///
8.27 +/// \file
8.28 +/// \brief Karp's algorithm for finding a minimum mean cycle.
8.29 +
8.30 +#include <vector>
8.31 +#include <limits>
8.32 +#include <lemon/core.h>
8.33 +#include <lemon/path.h>
8.34 +#include <lemon/tolerance.h>
8.35 +#include <lemon/connectivity.h>
8.36 +
8.37 +namespace lemon {
8.38 +
8.39 + /// \brief Default traits class of KarpMmc class.
8.40 + ///
8.41 + /// Default traits class of KarpMmc class.
8.42 + /// \tparam GR The type of the digraph.
8.43 + /// \tparam CM The type of the cost map.
8.44 + /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
8.45 +#ifdef DOXYGEN
8.46 + template <typename GR, typename CM>
8.47 +#else
8.48 + template <typename GR, typename CM,
8.49 + bool integer = std::numeric_limits<typename CM::Value>::is_integer>
8.50 +#endif
8.51 + struct KarpMmcDefaultTraits
8.52 + {
8.53 + /// The type of the digraph
8.54 + typedef GR Digraph;
8.55 + /// The type of the cost map
8.56 + typedef CM CostMap;
8.57 + /// The type of the arc costs
8.58 + typedef typename CostMap::Value Cost;
8.59 +
8.60 + /// \brief The large cost type used for internal computations
8.61 + ///
8.62 + /// The large cost type used for internal computations.
8.63 + /// It is \c long \c long if the \c Cost type is integer,
8.64 + /// otherwise it is \c double.
8.65 + /// \c Cost must be convertible to \c LargeCost.
8.66 + typedef double LargeCost;
8.67 +
8.68 + /// The tolerance type used for internal computations
8.69 + typedef lemon::Tolerance<LargeCost> Tolerance;
8.70 +
8.71 + /// \brief The path type of the found cycles
8.72 + ///
8.73 + /// The path type of the found cycles.
8.74 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
8.75 + /// and it must have an \c addFront() function.
8.76 + typedef lemon::Path<Digraph> Path;
8.77 + };
8.78 +
8.79 + // Default traits class for integer cost types
8.80 + template <typename GR, typename CM>
8.81 + struct KarpMmcDefaultTraits<GR, CM, true>
8.82 + {
8.83 + typedef GR Digraph;
8.84 + typedef CM CostMap;
8.85 + typedef typename CostMap::Value Cost;
8.86 +#ifdef LEMON_HAVE_LONG_LONG
8.87 + typedef long long LargeCost;
8.88 +#else
8.89 + typedef long LargeCost;
8.90 +#endif
8.91 + typedef lemon::Tolerance<LargeCost> Tolerance;
8.92 + typedef lemon::Path<Digraph> Path;
8.93 + };
8.94 +
8.95 +
8.96 + /// \addtogroup min_mean_cycle
8.97 + /// @{
8.98 +
8.99 + /// \brief Implementation of Karp's algorithm for finding a minimum
8.100 + /// mean cycle.
8.101 + ///
8.102 + /// This class implements Karp's algorithm for finding a directed
8.103 + /// cycle of minimum mean cost in a digraph
8.104 + /// \ref amo93networkflows, \ref dasdan98minmeancycle.
8.105 + /// It runs in time O(ne) and uses space O(n<sup>2</sup>+e).
8.106 + ///
8.107 + /// \tparam GR The type of the digraph the algorithm runs on.
8.108 + /// \tparam CM The type of the cost map. The default
8.109 + /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
8.110 + /// \tparam TR The traits class that defines various types used by the
8.111 + /// algorithm. By default, it is \ref KarpMmcDefaultTraits
8.112 + /// "KarpMmcDefaultTraits<GR, CM>".
8.113 + /// In most cases, this parameter should not be set directly,
8.114 + /// consider to use the named template parameters instead.
8.115 +#ifdef DOXYGEN
8.116 + template <typename GR, typename CM, typename TR>
8.117 +#else
8.118 + template < typename GR,
8.119 + typename CM = typename GR::template ArcMap<int>,
8.120 + typename TR = KarpMmcDefaultTraits<GR, CM> >
8.121 +#endif
8.122 + class KarpMmc
8.123 + {
8.124 + public:
8.125 +
8.126 + /// The type of the digraph
8.127 + typedef typename TR::Digraph Digraph;
8.128 + /// The type of the cost map
8.129 + typedef typename TR::CostMap CostMap;
8.130 + /// The type of the arc costs
8.131 + typedef typename TR::Cost Cost;
8.132 +
8.133 + /// \brief The large cost type
8.134 + ///
8.135 + /// The large cost type used for internal computations.
8.136 + /// By default, it is \c long \c long if the \c Cost type is integer,
8.137 + /// otherwise it is \c double.
8.138 + typedef typename TR::LargeCost LargeCost;
8.139 +
8.140 + /// The tolerance type
8.141 + typedef typename TR::Tolerance Tolerance;
8.142 +
8.143 + /// \brief The path type of the found cycles
8.144 + ///
8.145 + /// The path type of the found cycles.
8.146 + /// Using the \ref KarpMmcDefaultTraits "default traits class",
8.147 + /// it is \ref lemon::Path "Path<Digraph>".
8.148 + typedef typename TR::Path Path;
8.149 +
8.150 + /// The \ref KarpMmcDefaultTraits "traits class" of the algorithm
8.151 + typedef TR Traits;
8.152 +
8.153 + private:
8.154 +
8.155 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
8.156 +
8.157 + // Data sturcture for path data
8.158 + struct PathData
8.159 + {
8.160 + LargeCost dist;
8.161 + Arc pred;
8.162 + PathData(LargeCost d, Arc p = INVALID) :
8.163 + dist(d), pred(p) {}
8.164 + };
8.165 +
8.166 + typedef typename Digraph::template NodeMap<std::vector<PathData> >
8.167 + PathDataNodeMap;
8.168 +
8.169 + private:
8.170 +
8.171 + // The digraph the algorithm runs on
8.172 + const Digraph &_gr;
8.173 + // The cost of the arcs
8.174 + const CostMap &_cost;
8.175 +
8.176 + // Data for storing the strongly connected components
8.177 + int _comp_num;
8.178 + typename Digraph::template NodeMap<int> _comp;
8.179 + std::vector<std::vector<Node> > _comp_nodes;
8.180 + std::vector<Node>* _nodes;
8.181 + typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
8.182 +
8.183 + // Data for the found cycle
8.184 + LargeCost _cycle_cost;
8.185 + int _cycle_size;
8.186 + Node _cycle_node;
8.187 +
8.188 + Path *_cycle_path;
8.189 + bool _local_path;
8.190 +
8.191 + // Node map for storing path data
8.192 + PathDataNodeMap _data;
8.193 + // The processed nodes in the last round
8.194 + std::vector<Node> _process;
8.195 +
8.196 + Tolerance _tolerance;
8.197 +
8.198 + // Infinite constant
8.199 + const LargeCost INF;
8.200 +
8.201 + public:
8.202 +
8.203 + /// \name Named Template Parameters
8.204 + /// @{
8.205 +
8.206 + template <typename T>
8.207 + struct SetLargeCostTraits : public Traits {
8.208 + typedef T LargeCost;
8.209 + typedef lemon::Tolerance<T> Tolerance;
8.210 + };
8.211 +
8.212 + /// \brief \ref named-templ-param "Named parameter" for setting
8.213 + /// \c LargeCost type.
8.214 + ///
8.215 + /// \ref named-templ-param "Named parameter" for setting \c LargeCost
8.216 + /// type. It is used for internal computations in the algorithm.
8.217 + template <typename T>
8.218 + struct SetLargeCost
8.219 + : public KarpMmc<GR, CM, SetLargeCostTraits<T> > {
8.220 + typedef KarpMmc<GR, CM, SetLargeCostTraits<T> > Create;
8.221 + };
8.222 +
8.223 + template <typename T>
8.224 + struct SetPathTraits : public Traits {
8.225 + typedef T Path;
8.226 + };
8.227 +
8.228 + /// \brief \ref named-templ-param "Named parameter" for setting
8.229 + /// \c %Path type.
8.230 + ///
8.231 + /// \ref named-templ-param "Named parameter" for setting the \c %Path
8.232 + /// type of the found cycles.
8.233 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
8.234 + /// and it must have an \c addFront() function.
8.235 + template <typename T>
8.236 + struct SetPath
8.237 + : public KarpMmc<GR, CM, SetPathTraits<T> > {
8.238 + typedef KarpMmc<GR, CM, SetPathTraits<T> > Create;
8.239 + };
8.240 +
8.241 + /// @}
8.242 +
8.243 + protected:
8.244 +
8.245 + KarpMmc() {}
8.246 +
8.247 + public:
8.248 +
8.249 + /// \brief Constructor.
8.250 + ///
8.251 + /// The constructor of the class.
8.252 + ///
8.253 + /// \param digraph The digraph the algorithm runs on.
8.254 + /// \param cost The costs of the arcs.
8.255 + KarpMmc( const Digraph &digraph,
8.256 + const CostMap &cost ) :
8.257 + _gr(digraph), _cost(cost), _comp(digraph), _out_arcs(digraph),
8.258 + _cycle_cost(0), _cycle_size(1), _cycle_node(INVALID),
8.259 + _cycle_path(NULL), _local_path(false), _data(digraph),
8.260 + INF(std::numeric_limits<LargeCost>::has_infinity ?
8.261 + std::numeric_limits<LargeCost>::infinity() :
8.262 + std::numeric_limits<LargeCost>::max())
8.263 + {}
8.264 +
8.265 + /// Destructor.
8.266 + ~KarpMmc() {
8.267 + if (_local_path) delete _cycle_path;
8.268 + }
8.269 +
8.270 + /// \brief Set the path structure for storing the found cycle.
8.271 + ///
8.272 + /// This function sets an external path structure for storing the
8.273 + /// found cycle.
8.274 + ///
8.275 + /// If you don't call this function before calling \ref run() or
8.276 + /// \ref findCycleMean(), it will allocate a local \ref Path "path"
8.277 + /// structure. The destuctor deallocates this automatically
8.278 + /// allocated object, of course.
8.279 + ///
8.280 + /// \note The algorithm calls only the \ref lemon::Path::addFront()
8.281 + /// "addFront()" function of the given path structure.
8.282 + ///
8.283 + /// \return <tt>(*this)</tt>
8.284 + KarpMmc& cycle(Path &path) {
8.285 + if (_local_path) {
8.286 + delete _cycle_path;
8.287 + _local_path = false;
8.288 + }
8.289 + _cycle_path = &path;
8.290 + return *this;
8.291 + }
8.292 +
8.293 + /// \brief Set the tolerance used by the algorithm.
8.294 + ///
8.295 + /// This function sets the tolerance object used by the algorithm.
8.296 + ///
8.297 + /// \return <tt>(*this)</tt>
8.298 + KarpMmc& tolerance(const Tolerance& tolerance) {
8.299 + _tolerance = tolerance;
8.300 + return *this;
8.301 + }
8.302 +
8.303 + /// \brief Return a const reference to the tolerance.
8.304 + ///
8.305 + /// This function returns a const reference to the tolerance object
8.306 + /// used by the algorithm.
8.307 + const Tolerance& tolerance() const {
8.308 + return _tolerance;
8.309 + }
8.310 +
8.311 + /// \name Execution control
8.312 + /// The simplest way to execute the algorithm is to call the \ref run()
8.313 + /// function.\n
8.314 + /// If you only need the minimum mean cost, you may call
8.315 + /// \ref findCycleMean().
8.316 +
8.317 + /// @{
8.318 +
8.319 + /// \brief Run the algorithm.
8.320 + ///
8.321 + /// This function runs the algorithm.
8.322 + /// It can be called more than once (e.g. if the underlying digraph
8.323 + /// and/or the arc costs have been modified).
8.324 + ///
8.325 + /// \return \c true if a directed cycle exists in the digraph.
8.326 + ///
8.327 + /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
8.328 + /// \code
8.329 + /// return mmc.findCycleMean() && mmc.findCycle();
8.330 + /// \endcode
8.331 + bool run() {
8.332 + return findCycleMean() && findCycle();
8.333 + }
8.334 +
8.335 + /// \brief Find the minimum cycle mean.
8.336 + ///
8.337 + /// This function finds the minimum mean cost of the directed
8.338 + /// cycles in the digraph.
8.339 + ///
8.340 + /// \return \c true if a directed cycle exists in the digraph.
8.341 + bool findCycleMean() {
8.342 + // Initialization and find strongly connected components
8.343 + init();
8.344 + findComponents();
8.345 +
8.346 + // Find the minimum cycle mean in the components
8.347 + for (int comp = 0; comp < _comp_num; ++comp) {
8.348 + if (!initComponent(comp)) continue;
8.349 + processRounds();
8.350 + updateMinMean();
8.351 + }
8.352 + return (_cycle_node != INVALID);
8.353 + }
8.354 +
8.355 + /// \brief Find a minimum mean directed cycle.
8.356 + ///
8.357 + /// This function finds a directed cycle of minimum mean cost
8.358 + /// in the digraph using the data computed by findCycleMean().
8.359 + ///
8.360 + /// \return \c true if a directed cycle exists in the digraph.
8.361 + ///
8.362 + /// \pre \ref findCycleMean() must be called before using this function.
8.363 + bool findCycle() {
8.364 + if (_cycle_node == INVALID) return false;
8.365 + IntNodeMap reached(_gr, -1);
8.366 + int r = _data[_cycle_node].size();
8.367 + Node u = _cycle_node;
8.368 + while (reached[u] < 0) {
8.369 + reached[u] = --r;
8.370 + u = _gr.source(_data[u][r].pred);
8.371 + }
8.372 + r = reached[u];
8.373 + Arc e = _data[u][r].pred;
8.374 + _cycle_path->addFront(e);
8.375 + _cycle_cost = _cost[e];
8.376 + _cycle_size = 1;
8.377 + Node v;
8.378 + while ((v = _gr.source(e)) != u) {
8.379 + e = _data[v][--r].pred;
8.380 + _cycle_path->addFront(e);
8.381 + _cycle_cost += _cost[e];
8.382 + ++_cycle_size;
8.383 + }
8.384 + return true;
8.385 + }
8.386 +
8.387 + /// @}
8.388 +
8.389 + /// \name Query Functions
8.390 + /// The results of the algorithm can be obtained using these
8.391 + /// functions.\n
8.392 + /// The algorithm should be executed before using them.
8.393 +
8.394 + /// @{
8.395 +
8.396 + /// \brief Return the total cost of the found cycle.
8.397 + ///
8.398 + /// This function returns the total cost of the found cycle.
8.399 + ///
8.400 + /// \pre \ref run() or \ref findCycleMean() must be called before
8.401 + /// using this function.
8.402 + Cost cycleCost() const {
8.403 + return static_cast<Cost>(_cycle_cost);
8.404 + }
8.405 +
8.406 + /// \brief Return the number of arcs on the found cycle.
8.407 + ///
8.408 + /// This function returns the number of arcs on the found cycle.
8.409 + ///
8.410 + /// \pre \ref run() or \ref findCycleMean() must be called before
8.411 + /// using this function.
8.412 + int cycleSize() const {
8.413 + return _cycle_size;
8.414 + }
8.415 +
8.416 + /// \brief Return the mean cost of the found cycle.
8.417 + ///
8.418 + /// This function returns the mean cost of the found cycle.
8.419 + ///
8.420 + /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
8.421 + /// following code.
8.422 + /// \code
8.423 + /// return static_cast<double>(alg.cycleCost()) / alg.cycleSize();
8.424 + /// \endcode
8.425 + ///
8.426 + /// \pre \ref run() or \ref findCycleMean() must be called before
8.427 + /// using this function.
8.428 + double cycleMean() const {
8.429 + return static_cast<double>(_cycle_cost) / _cycle_size;
8.430 + }
8.431 +
8.432 + /// \brief Return the found cycle.
8.433 + ///
8.434 + /// This function returns a const reference to the path structure
8.435 + /// storing the found cycle.
8.436 + ///
8.437 + /// \pre \ref run() or \ref findCycle() must be called before using
8.438 + /// this function.
8.439 + const Path& cycle() const {
8.440 + return *_cycle_path;
8.441 + }
8.442 +
8.443 + ///@}
8.444 +
8.445 + private:
8.446 +
8.447 + // Initialization
8.448 + void init() {
8.449 + if (!_cycle_path) {
8.450 + _local_path = true;
8.451 + _cycle_path = new Path;
8.452 + }
8.453 + _cycle_path->clear();
8.454 + _cycle_cost = 0;
8.455 + _cycle_size = 1;
8.456 + _cycle_node = INVALID;
8.457 + for (NodeIt u(_gr); u != INVALID; ++u)
8.458 + _data[u].clear();
8.459 + }
8.460 +
8.461 + // Find strongly connected components and initialize _comp_nodes
8.462 + // and _out_arcs
8.463 + void findComponents() {
8.464 + _comp_num = stronglyConnectedComponents(_gr, _comp);
8.465 + _comp_nodes.resize(_comp_num);
8.466 + if (_comp_num == 1) {
8.467 + _comp_nodes[0].clear();
8.468 + for (NodeIt n(_gr); n != INVALID; ++n) {
8.469 + _comp_nodes[0].push_back(n);
8.470 + _out_arcs[n].clear();
8.471 + for (OutArcIt a(_gr, n); a != INVALID; ++a) {
8.472 + _out_arcs[n].push_back(a);
8.473 + }
8.474 + }
8.475 + } else {
8.476 + for (int i = 0; i < _comp_num; ++i)
8.477 + _comp_nodes[i].clear();
8.478 + for (NodeIt n(_gr); n != INVALID; ++n) {
8.479 + int k = _comp[n];
8.480 + _comp_nodes[k].push_back(n);
8.481 + _out_arcs[n].clear();
8.482 + for (OutArcIt a(_gr, n); a != INVALID; ++a) {
8.483 + if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
8.484 + }
8.485 + }
8.486 + }
8.487 + }
8.488 +
8.489 + // Initialize path data for the current component
8.490 + bool initComponent(int comp) {
8.491 + _nodes = &(_comp_nodes[comp]);
8.492 + int n = _nodes->size();
8.493 + if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
8.494 + return false;
8.495 + }
8.496 + for (int i = 0; i < n; ++i) {
8.497 + _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
8.498 + }
8.499 + return true;
8.500 + }
8.501 +
8.502 + // Process all rounds of computing path data for the current component.
8.503 + // _data[v][k] is the cost of a shortest directed walk from the root
8.504 + // node to node v containing exactly k arcs.
8.505 + void processRounds() {
8.506 + Node start = (*_nodes)[0];
8.507 + _data[start][0] = PathData(0);
8.508 + _process.clear();
8.509 + _process.push_back(start);
8.510 +
8.511 + int k, n = _nodes->size();
8.512 + for (k = 1; k <= n && int(_process.size()) < n; ++k) {
8.513 + processNextBuildRound(k);
8.514 + }
8.515 + for ( ; k <= n; ++k) {
8.516 + processNextFullRound(k);
8.517 + }
8.518 + }
8.519 +
8.520 + // Process one round and rebuild _process
8.521 + void processNextBuildRound(int k) {
8.522 + std::vector<Node> next;
8.523 + Node u, v;
8.524 + Arc e;
8.525 + LargeCost d;
8.526 + for (int i = 0; i < int(_process.size()); ++i) {
8.527 + u = _process[i];
8.528 + for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
8.529 + e = _out_arcs[u][j];
8.530 + v = _gr.target(e);
8.531 + d = _data[u][k-1].dist + _cost[e];
8.532 + if (_tolerance.less(d, _data[v][k].dist)) {
8.533 + if (_data[v][k].dist == INF) next.push_back(v);
8.534 + _data[v][k] = PathData(d, e);
8.535 + }
8.536 + }
8.537 + }
8.538 + _process.swap(next);
8.539 + }
8.540 +
8.541 + // Process one round using _nodes instead of _process
8.542 + void processNextFullRound(int k) {
8.543 + Node u, v;
8.544 + Arc e;
8.545 + LargeCost d;
8.546 + for (int i = 0; i < int(_nodes->size()); ++i) {
8.547 + u = (*_nodes)[i];
8.548 + for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
8.549 + e = _out_arcs[u][j];
8.550 + v = _gr.target(e);
8.551 + d = _data[u][k-1].dist + _cost[e];
8.552 + if (_tolerance.less(d, _data[v][k].dist)) {
8.553 + _data[v][k] = PathData(d, e);
8.554 + }
8.555 + }
8.556 + }
8.557 + }
8.558 +
8.559 + // Update the minimum cycle mean
8.560 + void updateMinMean() {
8.561 + int n = _nodes->size();
8.562 + for (int i = 0; i < n; ++i) {
8.563 + Node u = (*_nodes)[i];
8.564 + if (_data[u][n].dist == INF) continue;
8.565 + LargeCost cost, max_cost = 0;
8.566 + int size, max_size = 1;
8.567 + bool found_curr = false;
8.568 + for (int k = 0; k < n; ++k) {
8.569 + if (_data[u][k].dist == INF) continue;
8.570 + cost = _data[u][n].dist - _data[u][k].dist;
8.571 + size = n - k;
8.572 + if (!found_curr || cost * max_size > max_cost * size) {
8.573 + found_curr = true;
8.574 + max_cost = cost;
8.575 + max_size = size;
8.576 + }
8.577 + }
8.578 + if ( found_curr && (_cycle_node == INVALID ||
8.579 + max_cost * _cycle_size < _cycle_cost * max_size) ) {
8.580 + _cycle_cost = max_cost;
8.581 + _cycle_size = max_size;
8.582 + _cycle_node = u;
8.583 + }
8.584 + }
8.585 + }
8.586 +
8.587 + }; //class KarpMmc
8.588 +
8.589 + ///@}
8.590 +
8.591 +} //namespace lemon
8.592 +
8.593 +#endif //LEMON_KARP_MMC_H
9.1 --- a/test/min_mean_cycle_test.cc Mon Mar 08 08:33:41 2010 +0100
9.2 +++ b/test/min_mean_cycle_test.cc Sat Mar 13 22:01:38 2010 +0100
9.3 @@ -25,9 +25,9 @@
9.4 #include <lemon/concepts/digraph.h>
9.5 #include <lemon/concept_check.h>
9.6
9.7 -#include <lemon/karp.h>
9.8 -#include <lemon/hartmann_orlin.h>
9.9 -#include <lemon/howard.h>
9.10 +#include <lemon/karp_mmc.h>
9.11 +#include <lemon/hartmann_orlin_mmc.h>
9.12 +#include <lemon/howard_mmc.h>
9.13
9.14 #include "test_tools.h"
9.15
9.16 @@ -63,7 +63,7 @@
9.17
9.18
9.19 // Check the interface of an MMC algorithm
9.20 -template <typename GR, typename Value>
9.21 +template <typename GR, typename Cost>
9.22 struct MmcClassConcept
9.23 {
9.24 template <typename MMC>
9.25 @@ -73,30 +73,30 @@
9.26
9.27 typedef typename MMC
9.28 ::template SetPath<ListPath<GR> >
9.29 - ::template SetLargeValue<Value>
9.30 + ::template SetLargeCost<Cost>
9.31 ::Create MmcAlg;
9.32 - MmcAlg mmc(me.g, me.length);
9.33 + MmcAlg mmc(me.g, me.cost);
9.34 const MmcAlg& const_mmc = mmc;
9.35
9.36 typename MmcAlg::Tolerance tol = const_mmc.tolerance();
9.37 mmc.tolerance(tol);
9.38
9.39 b = mmc.cycle(p).run();
9.40 - b = mmc.findMinMean();
9.41 + b = mmc.findCycleMean();
9.42 b = mmc.findCycle();
9.43
9.44 - v = const_mmc.cycleLength();
9.45 - i = const_mmc.cycleArcNum();
9.46 + v = const_mmc.cycleCost();
9.47 + i = const_mmc.cycleSize();
9.48 d = const_mmc.cycleMean();
9.49 p = const_mmc.cycle();
9.50 }
9.51
9.52 - typedef concepts::ReadMap<typename GR::Arc, Value> LM;
9.53 + typedef concepts::ReadMap<typename GR::Arc, Cost> CM;
9.54
9.55 GR g;
9.56 - LM length;
9.57 + CM cost;
9.58 ListPath<GR> p;
9.59 - Value v;
9.60 + Cost v;
9.61 int i;
9.62 double d;
9.63 bool b;
9.64 @@ -108,13 +108,13 @@
9.65 void checkMmcAlg(const SmartDigraph& gr,
9.66 const SmartDigraph::ArcMap<int>& lm,
9.67 const SmartDigraph::ArcMap<int>& cm,
9.68 - int length, int size) {
9.69 + int cost, int size) {
9.70 MMC alg(gr, lm);
9.71 - alg.findMinMean();
9.72 - check(alg.cycleMean() == static_cast<double>(length) / size,
9.73 + alg.findCycleMean();
9.74 + check(alg.cycleMean() == static_cast<double>(cost) / size,
9.75 "Wrong cycle mean");
9.76 alg.findCycle();
9.77 - check(alg.cycleLength() == length && alg.cycleArcNum() == size,
9.78 + check(alg.cycleCost() == cost && alg.cycleSize() == size,
9.79 "Wrong path");
9.80 SmartDigraph::ArcMap<int> cycle(gr, 0);
9.81 for (typename MMC::Path::ArcIt a(alg.cycle()); a != INVALID; ++a) {
9.82 @@ -148,28 +148,28 @@
9.83 {
9.84 typedef concepts::Digraph GR;
9.85
9.86 - // Karp
9.87 + // KarpMmc
9.88 checkConcept< MmcClassConcept<GR, int>,
9.89 - Karp<GR, concepts::ReadMap<GR::Arc, int> > >();
9.90 + KarpMmc<GR, concepts::ReadMap<GR::Arc, int> > >();
9.91 checkConcept< MmcClassConcept<GR, float>,
9.92 - Karp<GR, concepts::ReadMap<GR::Arc, float> > >();
9.93 + KarpMmc<GR, concepts::ReadMap<GR::Arc, float> > >();
9.94
9.95 - // HartmannOrlin
9.96 + // HartmannOrlinMmc
9.97 checkConcept< MmcClassConcept<GR, int>,
9.98 - HartmannOrlin<GR, concepts::ReadMap<GR::Arc, int> > >();
9.99 + HartmannOrlinMmc<GR, concepts::ReadMap<GR::Arc, int> > >();
9.100 checkConcept< MmcClassConcept<GR, float>,
9.101 - HartmannOrlin<GR, concepts::ReadMap<GR::Arc, float> > >();
9.102 + HartmannOrlinMmc<GR, concepts::ReadMap<GR::Arc, float> > >();
9.103
9.104 - // Howard
9.105 + // HowardMmc
9.106 checkConcept< MmcClassConcept<GR, int>,
9.107 - Howard<GR, concepts::ReadMap<GR::Arc, int> > >();
9.108 + HowardMmc<GR, concepts::ReadMap<GR::Arc, int> > >();
9.109 checkConcept< MmcClassConcept<GR, float>,
9.110 - Howard<GR, concepts::ReadMap<GR::Arc, float> > >();
9.111 + HowardMmc<GR, concepts::ReadMap<GR::Arc, float> > >();
9.112
9.113 - if (IsSameType<Howard<GR, concepts::ReadMap<GR::Arc, int> >::LargeValue,
9.114 - long_int>::result == 0) check(false, "Wrong LargeValue type");
9.115 - if (IsSameType<Howard<GR, concepts::ReadMap<GR::Arc, float> >::LargeValue,
9.116 - double>::result == 0) check(false, "Wrong LargeValue type");
9.117 + check((IsSameType<HowardMmc<GR, concepts::ReadMap<GR::Arc, int> >
9.118 + ::LargeCost, long_int>::result == 1), "Wrong LargeCost type");
9.119 + check((IsSameType<HowardMmc<GR, concepts::ReadMap<GR::Arc, float> >
9.120 + ::LargeCost, double>::result == 1), "Wrong LargeCost type");
9.121 }
9.122
9.123 // Run various tests
9.124 @@ -194,22 +194,22 @@
9.125 run();
9.126
9.127 // Karp
9.128 - checkMmcAlg<Karp<GR, IntArcMap> >(gr, l1, c1, 6, 3);
9.129 - checkMmcAlg<Karp<GR, IntArcMap> >(gr, l2, c2, 5, 2);
9.130 - checkMmcAlg<Karp<GR, IntArcMap> >(gr, l3, c3, 0, 1);
9.131 - checkMmcAlg<Karp<GR, IntArcMap> >(gr, l4, c4, -1, 1);
9.132 + checkMmcAlg<KarpMmc<GR, IntArcMap> >(gr, l1, c1, 6, 3);
9.133 + checkMmcAlg<KarpMmc<GR, IntArcMap> >(gr, l2, c2, 5, 2);
9.134 + checkMmcAlg<KarpMmc<GR, IntArcMap> >(gr, l3, c3, 0, 1);
9.135 + checkMmcAlg<KarpMmc<GR, IntArcMap> >(gr, l4, c4, -1, 1);
9.136
9.137 // HartmannOrlin
9.138 - checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l1, c1, 6, 3);
9.139 - checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l2, c2, 5, 2);
9.140 - checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l3, c3, 0, 1);
9.141 - checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l4, c4, -1, 1);
9.142 + checkMmcAlg<HartmannOrlinMmc<GR, IntArcMap> >(gr, l1, c1, 6, 3);
9.143 + checkMmcAlg<HartmannOrlinMmc<GR, IntArcMap> >(gr, l2, c2, 5, 2);
9.144 + checkMmcAlg<HartmannOrlinMmc<GR, IntArcMap> >(gr, l3, c3, 0, 1);
9.145 + checkMmcAlg<HartmannOrlinMmc<GR, IntArcMap> >(gr, l4, c4, -1, 1);
9.146
9.147 // Howard
9.148 - checkMmcAlg<Howard<GR, IntArcMap> >(gr, l1, c1, 6, 3);
9.149 - checkMmcAlg<Howard<GR, IntArcMap> >(gr, l2, c2, 5, 2);
9.150 - checkMmcAlg<Howard<GR, IntArcMap> >(gr, l3, c3, 0, 1);
9.151 - checkMmcAlg<Howard<GR, IntArcMap> >(gr, l4, c4, -1, 1);
9.152 + checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l1, c1, 6, 3);
9.153 + checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l2, c2, 5, 2);
9.154 + checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l3, c3, 0, 1);
9.155 + checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l4, c4, -1, 1);
9.156 }
9.157
9.158 return 0;