Add HartmannOrlin algorithm class (#179)
authorPeter Kovacs <kpeter@inf.elte.hu>
Tue, 11 Aug 2009 21:53:39 +0200 (2009-08-11)
changeset 76697744b6dabf8
parent 765 3b544a9c92db
child 767 11c946fa8d13
Add HartmannOrlin algorithm class (#179)
This algorithm is an improved version of Karp's original method,
it applies an efficient early termination scheme.
The interface is the same as Karp's and Howard's interface.
lemon/Makefile.am
lemon/hartmann_orlin.h
test/min_mean_cycle_test.cc
     1.1 --- a/lemon/Makefile.am	Tue Aug 11 20:55:40 2009 +0200
     1.2 +++ b/lemon/Makefile.am	Tue Aug 11 21:53:39 2009 +0200
     1.3 @@ -83,6 +83,7 @@
     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/hypercube_graph.h \
    1.10  	lemon/karp.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/hartmann_orlin.h	Tue Aug 11 21:53:39 2009 +0200
     2.3 @@ -0,0 +1,618 @@
     2.4 +/* -*- C++ -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library
     2.7 + *
     2.8 + * Copyright (C) 2003-2008
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_HARTMANN_ORLIN_H
    2.23 +#define LEMON_HARTMANN_ORLIN_H
    2.24 +
    2.25 +/// \ingroup shortest_path
    2.26 +///
    2.27 +/// \file
    2.28 +/// \brief Hartmann-Orlin's algorithm for finding a minimum mean cycle.
    2.29 +
    2.30 +#include <vector>
    2.31 +#include <limits>
    2.32 +#include <lemon/core.h>
    2.33 +#include <lemon/path.h>
    2.34 +#include <lemon/tolerance.h>
    2.35 +#include <lemon/connectivity.h>
    2.36 +
    2.37 +namespace lemon {
    2.38 +
    2.39 +  /// \brief Default traits class of HartmannOrlin algorithm.
    2.40 +  ///
    2.41 +  /// Default traits class of HartmannOrlin algorithm.
    2.42 +  /// \tparam GR The type of the digraph.
    2.43 +  /// \tparam LEN The type of the length map.
    2.44 +  /// It must conform to the \ref concepts::Rea_data "Rea_data" concept.
    2.45 +#ifdef DOXYGEN
    2.46 +  template <typename GR, typename LEN>
    2.47 +#else
    2.48 +  template <typename GR, typename LEN,
    2.49 +    bool integer = std::numeric_limits<typename LEN::Value>::is_integer>
    2.50 +#endif
    2.51 +  struct HartmannOrlinDefaultTraits
    2.52 +  {
    2.53 +    /// The type of the digraph
    2.54 +    typedef GR Digraph;
    2.55 +    /// The type of the length map
    2.56 +    typedef LEN LengthMap;
    2.57 +    /// The type of the arc lengths
    2.58 +    typedef typename LengthMap::Value Value;
    2.59 +
    2.60 +    /// \brief The large value type used for internal computations
    2.61 +    ///
    2.62 +    /// The large value type used for internal computations.
    2.63 +    /// It is \c long \c long if the \c Value type is integer,
    2.64 +    /// otherwise it is \c double.
    2.65 +    /// \c Value must be convertible to \c LargeValue.
    2.66 +    typedef double LargeValue;
    2.67 +
    2.68 +    /// The tolerance type used for internal computations
    2.69 +    typedef lemon::Tolerance<LargeValue> Tolerance;
    2.70 +
    2.71 +    /// \brief The path type of the found cycles
    2.72 +    ///
    2.73 +    /// The path type of the found cycles.
    2.74 +    /// It must conform to the \ref lemon::concepts::Path "Path" concept
    2.75 +    /// and it must have an \c addBack() function.
    2.76 +    typedef lemon::Path<Digraph> Path;
    2.77 +  };
    2.78 +
    2.79 +  // Default traits class for integer value types
    2.80 +  template <typename GR, typename LEN>
    2.81 +  struct HartmannOrlinDefaultTraits<GR, LEN, true>
    2.82 +  {
    2.83 +    typedef GR Digraph;
    2.84 +    typedef LEN LengthMap;
    2.85 +    typedef typename LengthMap::Value Value;
    2.86 +#ifdef LEMON_HAVE_LONG_LONG
    2.87 +    typedef long long LargeValue;
    2.88 +#else
    2.89 +    typedef long LargeValue;
    2.90 +#endif
    2.91 +    typedef lemon::Tolerance<LargeValue> Tolerance;
    2.92 +    typedef lemon::Path<Digraph> Path;
    2.93 +  };
    2.94 +
    2.95 +
    2.96 +  /// \addtogroup shortest_path
    2.97 +  /// @{
    2.98 +
    2.99 +  /// \brief Implementation of the Hartmann-Orlin algorithm for finding
   2.100 +  /// a minimum mean cycle.
   2.101 +  ///
   2.102 +  /// This class implements the Hartmann-Orlin algorithm for finding
   2.103 +  /// a directed cycle of minimum mean length (cost) in a digraph.
   2.104 +  /// It is an improved version of \ref Karp "Karp's original algorithm",
   2.105 +  /// it applies an efficient early termination scheme.
   2.106 +  ///
   2.107 +  /// \tparam GR The type of the digraph the algorithm runs on.
   2.108 +  /// \tparam LEN The type of the length map. The default
   2.109 +  /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
   2.110 +#ifdef DOXYGEN
   2.111 +  template <typename GR, typename LEN, typename TR>
   2.112 +#else
   2.113 +  template < typename GR,
   2.114 +             typename LEN = typename GR::template ArcMap<int>,
   2.115 +             typename TR = HartmannOrlinDefaultTraits<GR, LEN> >
   2.116 +#endif
   2.117 +  class HartmannOrlin
   2.118 +  {
   2.119 +  public:
   2.120 +
   2.121 +    /// The type of the digraph
   2.122 +    typedef typename TR::Digraph Digraph;
   2.123 +    /// The type of the length map
   2.124 +    typedef typename TR::LengthMap LengthMap;
   2.125 +    /// The type of the arc lengths
   2.126 +    typedef typename TR::Value Value;
   2.127 +
   2.128 +    /// \brief The large value type
   2.129 +    ///
   2.130 +    /// The large value type used for internal computations.
   2.131 +    /// Using the \ref HartmannOrlinDefaultTraits "default traits class",
   2.132 +    /// it is \c long \c long if the \c Value type is integer,
   2.133 +    /// otherwise it is \c double.
   2.134 +    typedef typename TR::LargeValue LargeValue;
   2.135 +
   2.136 +    /// The tolerance type
   2.137 +    typedef typename TR::Tolerance Tolerance;
   2.138 +
   2.139 +    /// \brief The path type of the found cycles
   2.140 +    ///
   2.141 +    /// The path type of the found cycles.
   2.142 +    /// Using the \ref HartmannOrlinDefaultTraits "default traits class",
   2.143 +    /// it is \ref lemon::Path "Path<Digraph>".
   2.144 +    typedef typename TR::Path Path;
   2.145 +
   2.146 +    /// The \ref HartmannOrlinDefaultTraits "traits class" of the algorithm
   2.147 +    typedef TR Traits;
   2.148 +
   2.149 +  private:
   2.150 +
   2.151 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   2.152 +
   2.153 +    // Data sturcture for path data
   2.154 +    struct PathData
   2.155 +    {
   2.156 +      bool found;
   2.157 +      LargeValue dist;
   2.158 +      Arc pred;
   2.159 +      PathData(bool f = false, LargeValue d = 0, Arc p = INVALID) :
   2.160 +        found(f), dist(d), pred(p) {}
   2.161 +    };
   2.162 +
   2.163 +    typedef typename Digraph::template NodeMap<std::vector<PathData> >
   2.164 +      PathDataNodeMap;
   2.165 +
   2.166 +  private:
   2.167 +
   2.168 +    // The digraph the algorithm runs on
   2.169 +    const Digraph &_gr;
   2.170 +    // The length of the arcs
   2.171 +    const LengthMap &_length;
   2.172 +
   2.173 +    // Data for storing the strongly connected components
   2.174 +    int _comp_num;
   2.175 +    typename Digraph::template NodeMap<int> _comp;
   2.176 +    std::vector<std::vector<Node> > _comp_nodes;
   2.177 +    std::vector<Node>* _nodes;
   2.178 +    typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
   2.179 +
   2.180 +    // Data for the found cycles
   2.181 +    bool _curr_found, _best_found;
   2.182 +    LargeValue _curr_length, _best_length;
   2.183 +    int _curr_size, _best_size;
   2.184 +    Node _curr_node, _best_node;
   2.185 +    int _curr_level, _best_level;
   2.186 +
   2.187 +    Path *_cycle_path;
   2.188 +    bool _local_path;
   2.189 +
   2.190 +    // Node map for storing path data
   2.191 +    PathDataNodeMap _data;
   2.192 +    // The processed nodes in the last round
   2.193 +    std::vector<Node> _process;
   2.194 +
   2.195 +    Tolerance _tolerance;
   2.196 +
   2.197 +  public:
   2.198 +
   2.199 +    /// \name Named Template Parameters
   2.200 +    /// @{
   2.201 +
   2.202 +    template <typename T>
   2.203 +    struct SetLargeValueTraits : public Traits {
   2.204 +      typedef T LargeValue;
   2.205 +      typedef lemon::Tolerance<T> Tolerance;
   2.206 +    };
   2.207 +
   2.208 +    /// \brief \ref named-templ-param "Named parameter" for setting
   2.209 +    /// \c LargeValue type.
   2.210 +    ///
   2.211 +    /// \ref named-templ-param "Named parameter" for setting \c LargeValue
   2.212 +    /// type. It is used for internal computations in the algorithm.
   2.213 +    template <typename T>
   2.214 +    struct SetLargeValue
   2.215 +      : public HartmannOrlin<GR, LEN, SetLargeValueTraits<T> > {
   2.216 +      typedef HartmannOrlin<GR, LEN, SetLargeValueTraits<T> > Create;
   2.217 +    };
   2.218 +
   2.219 +    template <typename T>
   2.220 +    struct SetPathTraits : public Traits {
   2.221 +      typedef T Path;
   2.222 +    };
   2.223 +
   2.224 +    /// \brief \ref named-templ-param "Named parameter" for setting
   2.225 +    /// \c %Path type.
   2.226 +    ///
   2.227 +    /// \ref named-templ-param "Named parameter" for setting the \c %Path
   2.228 +    /// type of the found cycles.
   2.229 +    /// It must conform to the \ref lemon::concepts::Path "Path" concept
   2.230 +    /// and it must have an \c addFront() function.
   2.231 +    template <typename T>
   2.232 +    struct SetPath
   2.233 +      : public HartmannOrlin<GR, LEN, SetPathTraits<T> > {
   2.234 +      typedef HartmannOrlin<GR, LEN, SetPathTraits<T> > Create;
   2.235 +    };
   2.236 +
   2.237 +    /// @}
   2.238 +
   2.239 +  public:
   2.240 +
   2.241 +    /// \brief Constructor.
   2.242 +    ///
   2.243 +    /// The constructor of the class.
   2.244 +    ///
   2.245 +    /// \param digraph The digraph the algorithm runs on.
   2.246 +    /// \param length The lengths (costs) of the arcs.
   2.247 +    HartmannOrlin( const Digraph &digraph,
   2.248 +                   const LengthMap &length ) :
   2.249 +      _gr(digraph), _length(length), _comp(digraph), _out_arcs(digraph),
   2.250 +      _best_found(false), _best_length(0), _best_size(1),
   2.251 +      _cycle_path(NULL), _local_path(false), _data(digraph)
   2.252 +    {}
   2.253 +
   2.254 +    /// Destructor.
   2.255 +    ~HartmannOrlin() {
   2.256 +      if (_local_path) delete _cycle_path;
   2.257 +    }
   2.258 +
   2.259 +    /// \brief Set the path structure for storing the found cycle.
   2.260 +    ///
   2.261 +    /// This function sets an external path structure for storing the
   2.262 +    /// found cycle.
   2.263 +    ///
   2.264 +    /// If you don't call this function before calling \ref run() or
   2.265 +    /// \ref findMinMean(), it will allocate a local \ref Path "path"
   2.266 +    /// structure. The destuctor deallocates this automatically
   2.267 +    /// allocated object, of course.
   2.268 +    ///
   2.269 +    /// \note The algorithm calls only the \ref lemon::Path::addFront()
   2.270 +    /// "addFront()" function of the given path structure.
   2.271 +    ///
   2.272 +    /// \return <tt>(*this)</tt>
   2.273 +    HartmannOrlin& cycle(Path &path) {
   2.274 +      if (_local_path) {
   2.275 +        delete _cycle_path;
   2.276 +        _local_path = false;
   2.277 +      }
   2.278 +      _cycle_path = &path;
   2.279 +      return *this;
   2.280 +    }
   2.281 +
   2.282 +    /// \name Execution control
   2.283 +    /// The simplest way to execute the algorithm is to call the \ref run()
   2.284 +    /// function.\n
   2.285 +    /// If you only need the minimum mean length, you may call
   2.286 +    /// \ref findMinMean().
   2.287 +
   2.288 +    /// @{
   2.289 +
   2.290 +    /// \brief Run the algorithm.
   2.291 +    ///
   2.292 +    /// This function runs the algorithm.
   2.293 +    /// It can be called more than once (e.g. if the underlying digraph
   2.294 +    /// and/or the arc lengths have been modified).
   2.295 +    ///
   2.296 +    /// \return \c true if a directed cycle exists in the digraph.
   2.297 +    ///
   2.298 +    /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
   2.299 +    /// \code
   2.300 +    ///   return mmc.findMinMean() && mmc.findCycle();
   2.301 +    /// \endcode
   2.302 +    bool run() {
   2.303 +      return findMinMean() && findCycle();
   2.304 +    }
   2.305 +
   2.306 +    /// \brief Find the minimum cycle mean.
   2.307 +    ///
   2.308 +    /// This function finds the minimum mean length of the directed
   2.309 +    /// cycles in the digraph.
   2.310 +    ///
   2.311 +    /// \return \c true if a directed cycle exists in the digraph.
   2.312 +    bool findMinMean() {
   2.313 +      // Initialization and find strongly connected components
   2.314 +      init();
   2.315 +      findComponents();
   2.316 +      
   2.317 +      // Find the minimum cycle mean in the components
   2.318 +      for (int comp = 0; comp < _comp_num; ++comp) {
   2.319 +        if (!initComponent(comp)) continue;
   2.320 +        processRounds();
   2.321 +        
   2.322 +        // Update the best cycle (global minimum mean cycle)
   2.323 +        if ( _curr_found && (!_best_found || 
   2.324 +             _curr_length * _best_size < _best_length * _curr_size) ) {
   2.325 +          _best_found = true;
   2.326 +          _best_length = _curr_length;
   2.327 +          _best_size = _curr_size;
   2.328 +          _best_node = _curr_node;
   2.329 +          _best_level = _curr_level;
   2.330 +        }
   2.331 +      }
   2.332 +      return _best_found;
   2.333 +    }
   2.334 +
   2.335 +    /// \brief Find a minimum mean directed cycle.
   2.336 +    ///
   2.337 +    /// This function finds a directed cycle of minimum mean length
   2.338 +    /// in the digraph using the data computed by findMinMean().
   2.339 +    ///
   2.340 +    /// \return \c true if a directed cycle exists in the digraph.
   2.341 +    ///
   2.342 +    /// \pre \ref findMinMean() must be called before using this function.
   2.343 +    bool findCycle() {
   2.344 +      if (!_best_found) return false;
   2.345 +      IntNodeMap reached(_gr, -1);
   2.346 +      int r = _best_level + 1;
   2.347 +      Node u = _best_node;
   2.348 +      while (reached[u] < 0) {
   2.349 +        reached[u] = --r;
   2.350 +        u = _gr.source(_data[u][r].pred);
   2.351 +      }
   2.352 +      r = reached[u];
   2.353 +      Arc e = _data[u][r].pred;
   2.354 +      _cycle_path->addFront(e);
   2.355 +      _best_length = _length[e];
   2.356 +      _best_size = 1;
   2.357 +      Node v;
   2.358 +      while ((v = _gr.source(e)) != u) {
   2.359 +        e = _data[v][--r].pred;
   2.360 +        _cycle_path->addFront(e);
   2.361 +        _best_length += _length[e];
   2.362 +        ++_best_size;
   2.363 +      }
   2.364 +      return true;
   2.365 +    }
   2.366 +
   2.367 +    /// @}
   2.368 +
   2.369 +    /// \name Query Functions
   2.370 +    /// The results of the algorithm can be obtained using these
   2.371 +    /// functions.\n
   2.372 +    /// The algorithm should be executed before using them.
   2.373 +
   2.374 +    /// @{
   2.375 +
   2.376 +    /// \brief Return the total length of the found cycle.
   2.377 +    ///
   2.378 +    /// This function returns the total length of the found cycle.
   2.379 +    ///
   2.380 +    /// \pre \ref run() or \ref findMinMean() must be called before
   2.381 +    /// using this function.
   2.382 +    LargeValue cycleLength() const {
   2.383 +      return _best_length;
   2.384 +    }
   2.385 +
   2.386 +    /// \brief Return the number of arcs on the found cycle.
   2.387 +    ///
   2.388 +    /// This function returns the number of arcs on the found cycle.
   2.389 +    ///
   2.390 +    /// \pre \ref run() or \ref findMinMean() must be called before
   2.391 +    /// using this function.
   2.392 +    int cycleArcNum() const {
   2.393 +      return _best_size;
   2.394 +    }
   2.395 +
   2.396 +    /// \brief Return the mean length of the found cycle.
   2.397 +    ///
   2.398 +    /// This function returns the mean length of the found cycle.
   2.399 +    ///
   2.400 +    /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
   2.401 +    /// following code.
   2.402 +    /// \code
   2.403 +    ///   return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
   2.404 +    /// \endcode
   2.405 +    ///
   2.406 +    /// \pre \ref run() or \ref findMinMean() must be called before
   2.407 +    /// using this function.
   2.408 +    double cycleMean() const {
   2.409 +      return static_cast<double>(_best_length) / _best_size;
   2.410 +    }
   2.411 +
   2.412 +    /// \brief Return the found cycle.
   2.413 +    ///
   2.414 +    /// This function returns a const reference to the path structure
   2.415 +    /// storing the found cycle.
   2.416 +    ///
   2.417 +    /// \pre \ref run() or \ref findCycle() must be called before using
   2.418 +    /// this function.
   2.419 +    const Path& cycle() const {
   2.420 +      return *_cycle_path;
   2.421 +    }
   2.422 +
   2.423 +    ///@}
   2.424 +
   2.425 +  private:
   2.426 +
   2.427 +    // Initialization
   2.428 +    void init() {
   2.429 +      if (!_cycle_path) {
   2.430 +        _local_path = true;
   2.431 +        _cycle_path = new Path;
   2.432 +      }
   2.433 +      _cycle_path->clear();
   2.434 +      _best_found = false;
   2.435 +      _best_length = 0;
   2.436 +      _best_size = 1;
   2.437 +      _cycle_path->clear();
   2.438 +      for (NodeIt u(_gr); u != INVALID; ++u)
   2.439 +        _data[u].clear();
   2.440 +    }
   2.441 +
   2.442 +    // Find strongly connected components and initialize _comp_nodes
   2.443 +    // and _out_arcs
   2.444 +    void findComponents() {
   2.445 +      _comp_num = stronglyConnectedComponents(_gr, _comp);
   2.446 +      _comp_nodes.resize(_comp_num);
   2.447 +      if (_comp_num == 1) {
   2.448 +        _comp_nodes[0].clear();
   2.449 +        for (NodeIt n(_gr); n != INVALID; ++n) {
   2.450 +          _comp_nodes[0].push_back(n);
   2.451 +          _out_arcs[n].clear();
   2.452 +          for (OutArcIt a(_gr, n); a != INVALID; ++a) {
   2.453 +            _out_arcs[n].push_back(a);
   2.454 +          }
   2.455 +        }
   2.456 +      } else {
   2.457 +        for (int i = 0; i < _comp_num; ++i)
   2.458 +          _comp_nodes[i].clear();
   2.459 +        for (NodeIt n(_gr); n != INVALID; ++n) {
   2.460 +          int k = _comp[n];
   2.461 +          _comp_nodes[k].push_back(n);
   2.462 +          _out_arcs[n].clear();
   2.463 +          for (OutArcIt a(_gr, n); a != INVALID; ++a) {
   2.464 +            if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
   2.465 +          }
   2.466 +        }
   2.467 +      }
   2.468 +    }
   2.469 +
   2.470 +    // Initialize path data for the current component
   2.471 +    bool initComponent(int comp) {
   2.472 +      _nodes = &(_comp_nodes[comp]);
   2.473 +      int n = _nodes->size();
   2.474 +      if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
   2.475 +        return false;
   2.476 +      }      
   2.477 +      for (int i = 0; i < n; ++i) {
   2.478 +        _data[(*_nodes)[i]].resize(n + 1);
   2.479 +      }
   2.480 +      return true;
   2.481 +    }
   2.482 +
   2.483 +    // Process all rounds of computing path data for the current component.
   2.484 +    // _data[v][k] is the length of a shortest directed walk from the root
   2.485 +    // node to node v containing exactly k arcs.
   2.486 +    void processRounds() {
   2.487 +      Node start = (*_nodes)[0];
   2.488 +      _data[start][0] = PathData(true, 0);
   2.489 +      _process.clear();
   2.490 +      _process.push_back(start);
   2.491 +
   2.492 +      int k, n = _nodes->size();
   2.493 +      int next_check = 4;
   2.494 +      bool terminate = false;
   2.495 +      for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {
   2.496 +        processNextBuildRound(k);
   2.497 +        if (k == next_check || k == n) {
   2.498 +          terminate = checkTermination(k);
   2.499 +          next_check = next_check * 3 / 2;
   2.500 +        }
   2.501 +      }
   2.502 +      for ( ; k <= n && !terminate; ++k) {
   2.503 +        processNextFullRound(k);
   2.504 +        if (k == next_check || k == n) {
   2.505 +          terminate = checkTermination(k);
   2.506 +          next_check = next_check * 3 / 2;
   2.507 +        }
   2.508 +      }
   2.509 +    }
   2.510 +
   2.511 +    // Process one round and rebuild _process
   2.512 +    void processNextBuildRound(int k) {
   2.513 +      std::vector<Node> next;
   2.514 +      Node u, v;
   2.515 +      Arc e;
   2.516 +      LargeValue d;
   2.517 +      for (int i = 0; i < int(_process.size()); ++i) {
   2.518 +        u = _process[i];
   2.519 +        for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
   2.520 +          e = _out_arcs[u][j];
   2.521 +          v = _gr.target(e);
   2.522 +          d = _data[u][k-1].dist + _length[e];
   2.523 +          if (!_data[v][k].found) {
   2.524 +            next.push_back(v);
   2.525 +            _data[v][k] = PathData(true, _data[u][k-1].dist + _length[e], e);
   2.526 +          }
   2.527 +          else if (_tolerance.less(d, _data[v][k].dist)) {
   2.528 +            _data[v][k] = PathData(true, d, e);
   2.529 +          }
   2.530 +        }
   2.531 +      }
   2.532 +      _process.swap(next);
   2.533 +    }
   2.534 +
   2.535 +    // Process one round using _nodes instead of _process
   2.536 +    void processNextFullRound(int k) {
   2.537 +      Node u, v;
   2.538 +      Arc e;
   2.539 +      LargeValue d;
   2.540 +      for (int i = 0; i < int(_nodes->size()); ++i) {
   2.541 +        u = (*_nodes)[i];
   2.542 +        for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
   2.543 +          e = _out_arcs[u][j];
   2.544 +          v = _gr.target(e);
   2.545 +          d = _data[u][k-1].dist + _length[e];
   2.546 +          if (!_data[v][k].found || _tolerance.less(d, _data[v][k].dist)) {
   2.547 +            _data[v][k] = PathData(true, d, e);
   2.548 +          }
   2.549 +        }
   2.550 +      }
   2.551 +    }
   2.552 +    
   2.553 +    // Check early termination
   2.554 +    bool checkTermination(int k) {
   2.555 +      typedef std::pair<int, int> Pair;
   2.556 +      typename GR::template NodeMap<Pair> level(_gr, Pair(-1, 0));
   2.557 +      typename GR::template NodeMap<LargeValue> pi(_gr);
   2.558 +      int n = _nodes->size();
   2.559 +      LargeValue length;
   2.560 +      int size;
   2.561 +      Node u;
   2.562 +      
   2.563 +      // Search for cycles that are already found
   2.564 +      _curr_found = false;
   2.565 +      for (int i = 0; i < n; ++i) {
   2.566 +        u = (*_nodes)[i];
   2.567 +        if (!_data[u][k].found) continue;
   2.568 +        for (int j = k; j >= 0; --j) {
   2.569 +          if (level[u].first == i && level[u].second > 0) {
   2.570 +            // A cycle is found
   2.571 +            length = _data[u][level[u].second].dist - _data[u][j].dist;
   2.572 +            size = level[u].second - j;
   2.573 +            if (!_curr_found || length * _curr_size < _curr_length * size) {
   2.574 +              _curr_length = length;
   2.575 +              _curr_size = size;
   2.576 +              _curr_node = u;
   2.577 +              _curr_level = level[u].second;
   2.578 +              _curr_found = true;
   2.579 +            }
   2.580 +          }
   2.581 +          level[u] = Pair(i, j);
   2.582 +          u = _gr.source(_data[u][j].pred);
   2.583 +        }
   2.584 +      }
   2.585 +
   2.586 +      // If at least one cycle is found, check the optimality condition
   2.587 +      LargeValue d;
   2.588 +      if (_curr_found && k < n) {
   2.589 +        // Find node potentials
   2.590 +        for (int i = 0; i < n; ++i) {
   2.591 +          u = (*_nodes)[i];
   2.592 +          pi[u] = std::numeric_limits<LargeValue>::max();
   2.593 +          for (int j = 0; j <= k; ++j) {
   2.594 +            d = _data[u][j].dist * _curr_size - j * _curr_length;
   2.595 +            if (_data[u][j].found && _tolerance.less(d, pi[u])) {
   2.596 +              pi[u] = d;
   2.597 +            }
   2.598 +          }
   2.599 +        }
   2.600 +
   2.601 +        // Check the optimality condition for all arcs
   2.602 +        bool done = true;
   2.603 +        for (ArcIt a(_gr); a != INVALID; ++a) {
   2.604 +          if (_tolerance.less(_length[a] * _curr_size - _curr_length,
   2.605 +                              pi[_gr.target(a)] - pi[_gr.source(a)]) ) {
   2.606 +            done = false;
   2.607 +            break;
   2.608 +          }
   2.609 +        }
   2.610 +        return done;
   2.611 +      }
   2.612 +      return (k == n);
   2.613 +    }
   2.614 +
   2.615 +  }; //class HartmannOrlin
   2.616 +
   2.617 +  ///@}
   2.618 +
   2.619 +} //namespace lemon
   2.620 +
   2.621 +#endif //LEMON_HARTMANN_ORLIN_H
     3.1 --- a/test/min_mean_cycle_test.cc	Tue Aug 11 20:55:40 2009 +0200
     3.2 +++ b/test/min_mean_cycle_test.cc	Tue Aug 11 21:53:39 2009 +0200
     3.3 @@ -26,6 +26,7 @@
     3.4  #include <lemon/concept_check.h>
     3.5  
     3.6  #include <lemon/karp.h>
     3.7 +#include <lemon/hartmann_orlin.h>
     3.8  #include <lemon/howard.h>
     3.9  
    3.10  #include "test_tools.h"
    3.11 @@ -150,6 +151,12 @@
    3.12      checkConcept< MmcClassConcept<GR, float>,
    3.13                    Karp<GR, concepts::ReadMap<GR::Arc, float> > >();
    3.14      
    3.15 +    // HartmannOrlin
    3.16 +    checkConcept< MmcClassConcept<GR, int>,
    3.17 +                  HartmannOrlin<GR, concepts::ReadMap<GR::Arc, int> > >();
    3.18 +    checkConcept< MmcClassConcept<GR, float>,
    3.19 +                  HartmannOrlin<GR, concepts::ReadMap<GR::Arc, float> > >();
    3.20 +    
    3.21      // Howard
    3.22      checkConcept< MmcClassConcept<GR, int>,
    3.23                    Howard<GR, concepts::ReadMap<GR::Arc, int> > >();
    3.24 @@ -189,6 +196,12 @@
    3.25      checkMmcAlg<Karp<GR, IntArcMap> >(gr, l3, c3,  0, 1);
    3.26      checkMmcAlg<Karp<GR, IntArcMap> >(gr, l4, c4, -1, 1);
    3.27  
    3.28 +    // HartmannOrlin
    3.29 +    checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l1, c1,  6, 3);
    3.30 +    checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l2, c2,  5, 2);
    3.31 +    checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l3, c3,  0, 1);
    3.32 +    checkMmcAlg<HartmannOrlin<GR, IntArcMap> >(gr, l4, c4, -1, 1);
    3.33 +
    3.34      // Howard
    3.35      checkMmcAlg<Howard<GR, IntArcMap> >(gr, l1, c1,  6, 3);
    3.36      checkMmcAlg<Howard<GR, IntArcMap> >(gr, l2, c2,  5, 2);