Index: doc/groups.dox
===================================================================
--- doc/groups.dox (revision 802)
+++ doc/groups.dox (revision 817)
@@ -454,4 +454,38 @@
/**
+@defgroup min_mean_cycle Minimum Mean Cycle Algorithms
+@ingroup algs
+\brief Algorithms for finding minimum mean cycles.
+
+This group contains the algorithms for finding minimum mean cycles.
+
+The \e minimum \e mean \e cycle \e problem is to find a directed cycle
+of minimum mean length (cost) in a digraph.
+The mean length of a cycle is the average length of its arcs, i.e. the
+ratio between the total length of the cycle and the number of arcs on it.
+
+This problem has an important connection to \e conservative \e length
+\e functions, too. A length function on the arcs of a digraph is called
+conservative if and only if there is no directed cycle of negative total
+length. For an arbitrary length function, the negative of the minimum
+cycle mean is the smallest \f$\epsilon\f$ value so that increasing the
+arc lengths uniformly by \f$\epsilon\f$ results in a conservative length
+function.
+
+LEMON contains three algorithms for solving the minimum mean cycle problem:
+- \ref Karp "Karp"'s original algorithm.
+- \ref HartmannOrlin "Hartmann-Orlin"'s algorithm, which is an improved
+ version of Karp's algorithm.
+- \ref Howard "Howard"'s policy iteration algorithm.
+
+In practice, the Howard algorithm proved to be by far the most efficient
+one, though the best known theoretical bound on its running time is
+exponential.
+Both Karp and HartmannOrlin algorithms run in time O(ne) and use space
+O(n2+e), but the latter one is typically faster due to the
+applied early termination scheme.
+*/
+
+/**
@defgroup matching Matching Algorithms
@ingroup algs
Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am (revision 755)
+++ lemon/Makefile.am (revision 817)
@@ -87,5 +87,8 @@
lemon/graph_to_eps.h \
lemon/grid_graph.h \
+ lemon/hartmann_orlin.h \
+ lemon/howard.h \
lemon/hypercube_graph.h \
+ lemon/karp.h \
lemon/kary_heap.h \
lemon/kruskal.h \
Index: lemon/hartmann_orlin.h
===================================================================
--- lemon/hartmann_orlin.h (revision 816)
+++ lemon/hartmann_orlin.h (revision 816)
@@ -0,0 +1,639 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_HARTMANN_ORLIN_H
+#define LEMON_HARTMANN_ORLIN_H
+
+/// \ingroup min_mean_cycle
+///
+/// \file
+/// \brief Hartmann-Orlin's algorithm for finding a minimum mean cycle.
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace lemon {
+
+ /// \brief Default traits class of HartmannOrlin algorithm.
+ ///
+ /// Default traits class of HartmannOrlin algorithm.
+ /// \tparam GR The type of the digraph.
+ /// \tparam LEN The type of the length map.
+ /// It must conform to the \ref concepts::Rea_data "Rea_data" concept.
+#ifdef DOXYGEN
+ template
+#else
+ template ::is_integer>
+#endif
+ struct HartmannOrlinDefaultTraits
+ {
+ /// The type of the digraph
+ typedef GR Digraph;
+ /// The type of the length map
+ typedef LEN LengthMap;
+ /// The type of the arc lengths
+ typedef typename LengthMap::Value Value;
+
+ /// \brief The large value type used for internal computations
+ ///
+ /// The large value type used for internal computations.
+ /// It is \c long \c long if the \c Value type is integer,
+ /// otherwise it is \c double.
+ /// \c Value must be convertible to \c LargeValue.
+ typedef double LargeValue;
+
+ /// The tolerance type used for internal computations
+ typedef lemon::Tolerance Tolerance;
+
+ /// \brief The path type of the found cycles
+ ///
+ /// The path type of the found cycles.
+ /// It must conform to the \ref lemon::concepts::Path "Path" concept
+ /// and it must have an \c addBack() function.
+ typedef lemon::Path Path;
+ };
+
+ // Default traits class for integer value types
+ template
+ struct HartmannOrlinDefaultTraits
+ {
+ typedef GR Digraph;
+ typedef LEN LengthMap;
+ typedef typename LengthMap::Value Value;
+#ifdef LEMON_HAVE_LONG_LONG
+ typedef long long LargeValue;
+#else
+ typedef long LargeValue;
+#endif
+ typedef lemon::Tolerance Tolerance;
+ typedef lemon::Path Path;
+ };
+
+
+ /// \addtogroup min_mean_cycle
+ /// @{
+
+ /// \brief Implementation of the Hartmann-Orlin algorithm for finding
+ /// a minimum mean cycle.
+ ///
+ /// This class implements the Hartmann-Orlin algorithm for finding
+ /// a directed cycle of minimum mean length (cost) in a digraph.
+ /// It is an improved version of \ref Karp "Karp"'s original algorithm,
+ /// it applies an efficient early termination scheme.
+ /// It runs in time O(ne) and uses space O(n2+e).
+ ///
+ /// \tparam GR The type of the digraph the algorithm runs on.
+ /// \tparam LEN The type of the length map. The default
+ /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap".
+#ifdef DOXYGEN
+ template
+#else
+ template < typename GR,
+ typename LEN = typename GR::template ArcMap,
+ typename TR = HartmannOrlinDefaultTraits >
+#endif
+ class HartmannOrlin
+ {
+ public:
+
+ /// The type of the digraph
+ typedef typename TR::Digraph Digraph;
+ /// The type of the length map
+ typedef typename TR::LengthMap LengthMap;
+ /// The type of the arc lengths
+ typedef typename TR::Value Value;
+
+ /// \brief The large value type
+ ///
+ /// The large value type used for internal computations.
+ /// Using the \ref HartmannOrlinDefaultTraits "default traits class",
+ /// it is \c long \c long if the \c Value type is integer,
+ /// otherwise it is \c double.
+ typedef typename TR::LargeValue LargeValue;
+
+ /// The tolerance type
+ typedef typename TR::Tolerance Tolerance;
+
+ /// \brief The path type of the found cycles
+ ///
+ /// The path type of the found cycles.
+ /// Using the \ref HartmannOrlinDefaultTraits "default traits class",
+ /// it is \ref lemon::Path "Path".
+ typedef typename TR::Path Path;
+
+ /// The \ref HartmannOrlinDefaultTraits "traits class" of the algorithm
+ typedef TR Traits;
+
+ private:
+
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+ // Data sturcture for path data
+ struct PathData
+ {
+ LargeValue dist;
+ Arc pred;
+ PathData(LargeValue d, Arc p = INVALID) :
+ dist(d), pred(p) {}
+ };
+
+ typedef typename Digraph::template NodeMap >
+ PathDataNodeMap;
+
+ private:
+
+ // The digraph the algorithm runs on
+ const Digraph &_gr;
+ // The length of the arcs
+ const LengthMap &_length;
+
+ // Data for storing the strongly connected components
+ int _comp_num;
+ typename Digraph::template NodeMap _comp;
+ std::vector > _comp_nodes;
+ std::vector* _nodes;
+ typename Digraph::template NodeMap > _out_arcs;
+
+ // Data for the found cycles
+ bool _curr_found, _best_found;
+ LargeValue _curr_length, _best_length;
+ int _curr_size, _best_size;
+ Node _curr_node, _best_node;
+ int _curr_level, _best_level;
+
+ Path *_cycle_path;
+ bool _local_path;
+
+ // Node map for storing path data
+ PathDataNodeMap _data;
+ // The processed nodes in the last round
+ std::vector _process;
+
+ Tolerance _tolerance;
+
+ // Infinite constant
+ const LargeValue INF;
+
+ public:
+
+ /// \name Named Template Parameters
+ /// @{
+
+ template
+ struct SetLargeValueTraits : public Traits {
+ typedef T LargeValue;
+ typedef lemon::Tolerance Tolerance;
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// \c LargeValue type.
+ ///
+ /// \ref named-templ-param "Named parameter" for setting \c LargeValue
+ /// type. It is used for internal computations in the algorithm.
+ template
+ struct SetLargeValue
+ : public HartmannOrlin > {
+ typedef HartmannOrlin > Create;
+ };
+
+ template
+ struct SetPathTraits : public Traits {
+ typedef T Path;
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// \c %Path type.
+ ///
+ /// \ref named-templ-param "Named parameter" for setting the \c %Path
+ /// type of the found cycles.
+ /// It must conform to the \ref lemon::concepts::Path "Path" concept
+ /// and it must have an \c addFront() function.
+ template
+ struct SetPath
+ : public HartmannOrlin > {
+ typedef HartmannOrlin > Create;
+ };
+
+ /// @}
+
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// The constructor of the class.
+ ///
+ /// \param digraph The digraph the algorithm runs on.
+ /// \param length The lengths (costs) of the arcs.
+ HartmannOrlin( const Digraph &digraph,
+ const LengthMap &length ) :
+ _gr(digraph), _length(length), _comp(digraph), _out_arcs(digraph),
+ _best_found(false), _best_length(0), _best_size(1),
+ _cycle_path(NULL), _local_path(false), _data(digraph),
+ INF(std::numeric_limits::has_infinity ?
+ std::numeric_limits::infinity() :
+ std::numeric_limits::max())
+ {}
+
+ /// Destructor.
+ ~HartmannOrlin() {
+ if (_local_path) delete _cycle_path;
+ }
+
+ /// \brief Set the path structure for storing the found cycle.
+ ///
+ /// This function sets an external path structure for storing the
+ /// found cycle.
+ ///
+ /// If you don't call this function before calling \ref run() or
+ /// \ref findMinMean(), it will allocate a local \ref Path "path"
+ /// structure. The destuctor deallocates this automatically
+ /// allocated object, of course.
+ ///
+ /// \note The algorithm calls only the \ref lemon::Path::addFront()
+ /// "addFront()" function of the given path structure.
+ ///
+ /// \return (*this)
+ HartmannOrlin& cycle(Path &path) {
+ if (_local_path) {
+ delete _cycle_path;
+ _local_path = false;
+ }
+ _cycle_path = &path;
+ return *this;
+ }
+
+ /// \brief Set the tolerance used by the algorithm.
+ ///
+ /// This function sets the tolerance object used by the algorithm.
+ ///
+ /// \return (*this)
+ HartmannOrlin& tolerance(const Tolerance& tolerance) {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ /// \brief Return a const reference to the tolerance.
+ ///
+ /// This function returns a const reference to the tolerance object
+ /// used by the algorithm.
+ const Tolerance& tolerance() const {
+ return _tolerance;
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to call the \ref run()
+ /// function.\n
+ /// If you only need the minimum mean length, you may call
+ /// \ref findMinMean().
+
+ /// @{
+
+ /// \brief Run the algorithm.
+ ///
+ /// This function runs the algorithm.
+ /// It can be called more than once (e.g. if the underlying digraph
+ /// and/or the arc lengths have been modified).
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ ///
+ /// \note mmc.run() is just a shortcut of the following code.
+ /// \code
+ /// return mmc.findMinMean() && mmc.findCycle();
+ /// \endcode
+ bool run() {
+ return findMinMean() && findCycle();
+ }
+
+ /// \brief Find the minimum cycle mean.
+ ///
+ /// This function finds the minimum mean length of the directed
+ /// cycles in the digraph.
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ bool findMinMean() {
+ // Initialization and find strongly connected components
+ init();
+ findComponents();
+
+ // Find the minimum cycle mean in the components
+ for (int comp = 0; comp < _comp_num; ++comp) {
+ if (!initComponent(comp)) continue;
+ processRounds();
+
+ // Update the best cycle (global minimum mean cycle)
+ if ( _curr_found && (!_best_found ||
+ _curr_length * _best_size < _best_length * _curr_size) ) {
+ _best_found = true;
+ _best_length = _curr_length;
+ _best_size = _curr_size;
+ _best_node = _curr_node;
+ _best_level = _curr_level;
+ }
+ }
+ return _best_found;
+ }
+
+ /// \brief Find a minimum mean directed cycle.
+ ///
+ /// This function finds a directed cycle of minimum mean length
+ /// in the digraph using the data computed by findMinMean().
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ ///
+ /// \pre \ref findMinMean() must be called before using this function.
+ bool findCycle() {
+ if (!_best_found) return false;
+ IntNodeMap reached(_gr, -1);
+ int r = _best_level + 1;
+ Node u = _best_node;
+ while (reached[u] < 0) {
+ reached[u] = --r;
+ u = _gr.source(_data[u][r].pred);
+ }
+ r = reached[u];
+ Arc e = _data[u][r].pred;
+ _cycle_path->addFront(e);
+ _best_length = _length[e];
+ _best_size = 1;
+ Node v;
+ while ((v = _gr.source(e)) != u) {
+ e = _data[v][--r].pred;
+ _cycle_path->addFront(e);
+ _best_length += _length[e];
+ ++_best_size;
+ }
+ return true;
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The results of the algorithm can be obtained using these
+ /// functions.\n
+ /// The algorithm should be executed before using them.
+
+ /// @{
+
+ /// \brief Return the total length of the found cycle.
+ ///
+ /// This function returns the total length of the found cycle.
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ LargeValue cycleLength() const {
+ return _best_length;
+ }
+
+ /// \brief Return the number of arcs on the found cycle.
+ ///
+ /// This function returns the number of arcs on the found cycle.
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ int cycleArcNum() const {
+ return _best_size;
+ }
+
+ /// \brief Return the mean length of the found cycle.
+ ///
+ /// This function returns the mean length of the found cycle.
+ ///
+ /// \note alg.cycleMean() is just a shortcut of the
+ /// following code.
+ /// \code
+ /// return static_cast(alg.cycleLength()) / alg.cycleArcNum();
+ /// \endcode
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ double cycleMean() const {
+ return static_cast(_best_length) / _best_size;
+ }
+
+ /// \brief Return the found cycle.
+ ///
+ /// This function returns a const reference to the path structure
+ /// storing the found cycle.
+ ///
+ /// \pre \ref run() or \ref findCycle() must be called before using
+ /// this function.
+ const Path& cycle() const {
+ return *_cycle_path;
+ }
+
+ ///@}
+
+ private:
+
+ // Initialization
+ void init() {
+ if (!_cycle_path) {
+ _local_path = true;
+ _cycle_path = new Path;
+ }
+ _cycle_path->clear();
+ _best_found = false;
+ _best_length = 0;
+ _best_size = 1;
+ _cycle_path->clear();
+ for (NodeIt u(_gr); u != INVALID; ++u)
+ _data[u].clear();
+ }
+
+ // Find strongly connected components and initialize _comp_nodes
+ // and _out_arcs
+ void findComponents() {
+ _comp_num = stronglyConnectedComponents(_gr, _comp);
+ _comp_nodes.resize(_comp_num);
+ if (_comp_num == 1) {
+ _comp_nodes[0].clear();
+ for (NodeIt n(_gr); n != INVALID; ++n) {
+ _comp_nodes[0].push_back(n);
+ _out_arcs[n].clear();
+ for (OutArcIt a(_gr, n); a != INVALID; ++a) {
+ _out_arcs[n].push_back(a);
+ }
+ }
+ } else {
+ for (int i = 0; i < _comp_num; ++i)
+ _comp_nodes[i].clear();
+ for (NodeIt n(_gr); n != INVALID; ++n) {
+ int k = _comp[n];
+ _comp_nodes[k].push_back(n);
+ _out_arcs[n].clear();
+ for (OutArcIt a(_gr, n); a != INVALID; ++a) {
+ if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
+ }
+ }
+ }
+ }
+
+ // Initialize path data for the current component
+ bool initComponent(int comp) {
+ _nodes = &(_comp_nodes[comp]);
+ int n = _nodes->size();
+ if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
+ return false;
+ }
+ for (int i = 0; i < n; ++i) {
+ _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
+ }
+ return true;
+ }
+
+ // Process all rounds of computing path data for the current component.
+ // _data[v][k] is the length of a shortest directed walk from the root
+ // node to node v containing exactly k arcs.
+ void processRounds() {
+ Node start = (*_nodes)[0];
+ _data[start][0] = PathData(0);
+ _process.clear();
+ _process.push_back(start);
+
+ int k, n = _nodes->size();
+ int next_check = 4;
+ bool terminate = false;
+ for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {
+ processNextBuildRound(k);
+ if (k == next_check || k == n) {
+ terminate = checkTermination(k);
+ next_check = next_check * 3 / 2;
+ }
+ }
+ for ( ; k <= n && !terminate; ++k) {
+ processNextFullRound(k);
+ if (k == next_check || k == n) {
+ terminate = checkTermination(k);
+ next_check = next_check * 3 / 2;
+ }
+ }
+ }
+
+ // Process one round and rebuild _process
+ void processNextBuildRound(int k) {
+ std::vector next;
+ Node u, v;
+ Arc e;
+ LargeValue d;
+ for (int i = 0; i < int(_process.size()); ++i) {
+ u = _process[i];
+ for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
+ e = _out_arcs[u][j];
+ v = _gr.target(e);
+ d = _data[u][k-1].dist + _length[e];
+ if (_tolerance.less(d, _data[v][k].dist)) {
+ if (_data[v][k].dist == INF) next.push_back(v);
+ _data[v][k] = PathData(d, e);
+ }
+ }
+ }
+ _process.swap(next);
+ }
+
+ // Process one round using _nodes instead of _process
+ void processNextFullRound(int k) {
+ Node u, v;
+ Arc e;
+ LargeValue d;
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ u = (*_nodes)[i];
+ for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
+ e = _out_arcs[u][j];
+ v = _gr.target(e);
+ d = _data[u][k-1].dist + _length[e];
+ if (_tolerance.less(d, _data[v][k].dist)) {
+ _data[v][k] = PathData(d, e);
+ }
+ }
+ }
+ }
+
+ // Check early termination
+ bool checkTermination(int k) {
+ typedef std::pair Pair;
+ typename GR::template NodeMap level(_gr, Pair(-1, 0));
+ typename GR::template NodeMap pi(_gr);
+ int n = _nodes->size();
+ LargeValue length;
+ int size;
+ Node u;
+
+ // Search for cycles that are already found
+ _curr_found = false;
+ for (int i = 0; i < n; ++i) {
+ u = (*_nodes)[i];
+ if (_data[u][k].dist == INF) continue;
+ for (int j = k; j >= 0; --j) {
+ if (level[u].first == i && level[u].second > 0) {
+ // A cycle is found
+ length = _data[u][level[u].second].dist - _data[u][j].dist;
+ size = level[u].second - j;
+ if (!_curr_found || length * _curr_size < _curr_length * size) {
+ _curr_length = length;
+ _curr_size = size;
+ _curr_node = u;
+ _curr_level = level[u].second;
+ _curr_found = true;
+ }
+ }
+ level[u] = Pair(i, j);
+ u = _gr.source(_data[u][j].pred);
+ }
+ }
+
+ // If at least one cycle is found, check the optimality condition
+ LargeValue d;
+ if (_curr_found && k < n) {
+ // Find node potentials
+ for (int i = 0; i < n; ++i) {
+ u = (*_nodes)[i];
+ pi[u] = INF;
+ for (int j = 0; j <= k; ++j) {
+ if (_data[u][j].dist < INF) {
+ d = _data[u][j].dist * _curr_size - j * _curr_length;
+ if (_tolerance.less(d, pi[u])) pi[u] = d;
+ }
+ }
+ }
+
+ // Check the optimality condition for all arcs
+ bool done = true;
+ for (ArcIt a(_gr); a != INVALID; ++a) {
+ if (_tolerance.less(_length[a] * _curr_size - _curr_length,
+ pi[_gr.target(a)] - pi[_gr.source(a)]) ) {
+ done = false;
+ break;
+ }
+ }
+ return done;
+ }
+ return (k == n);
+ }
+
+ }; //class HartmannOrlin
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_HARTMANN_ORLIN_H
Index: lemon/howard.h
===================================================================
--- lemon/howard.h (revision 816)
+++ lemon/howard.h (revision 816)
@@ -0,0 +1,596 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_HOWARD_H
+#define LEMON_HOWARD_H
+
+/// \ingroup min_mean_cycle
+///
+/// \file
+/// \brief Howard's algorithm for finding a minimum mean cycle.
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace lemon {
+
+ /// \brief Default traits class of Howard class.
+ ///
+ /// Default traits class of Howard class.
+ /// \tparam GR The type of the digraph.
+ /// \tparam LEN The type of the length map.
+ /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
+#ifdef DOXYGEN
+ template
+#else
+ template ::is_integer>
+#endif
+ struct HowardDefaultTraits
+ {
+ /// The type of the digraph
+ typedef GR Digraph;
+ /// The type of the length map
+ typedef LEN LengthMap;
+ /// The type of the arc lengths
+ typedef typename LengthMap::Value Value;
+
+ /// \brief The large value type used for internal computations
+ ///
+ /// The large value type used for internal computations.
+ /// It is \c long \c long if the \c Value type is integer,
+ /// otherwise it is \c double.
+ /// \c Value must be convertible to \c LargeValue.
+ typedef double LargeValue;
+
+ /// The tolerance type used for internal computations
+ typedef lemon::Tolerance Tolerance;
+
+ /// \brief The path type of the found cycles
+ ///
+ /// The path type of the found cycles.
+ /// It must conform to the \ref lemon::concepts::Path "Path" concept
+ /// and it must have an \c addBack() function.
+ typedef lemon::Path Path;
+ };
+
+ // Default traits class for integer value types
+ template
+ struct HowardDefaultTraits
+ {
+ typedef GR Digraph;
+ typedef LEN LengthMap;
+ typedef typename LengthMap::Value Value;
+#ifdef LEMON_HAVE_LONG_LONG
+ typedef long long LargeValue;
+#else
+ typedef long LargeValue;
+#endif
+ typedef lemon::Tolerance Tolerance;
+ typedef lemon::Path Path;
+ };
+
+
+ /// \addtogroup min_mean_cycle
+ /// @{
+
+ /// \brief Implementation of Howard's algorithm for finding a minimum
+ /// mean cycle.
+ ///
+ /// This class implements Howard's policy iteration algorithm for finding
+ /// a directed cycle of minimum mean length (cost) in a digraph.
+ /// This class provides the most efficient algorithm for the
+ /// minimum mean cycle problem, though the best known theoretical
+ /// bound on its running time is exponential.
+ ///
+ /// \tparam GR The type of the digraph the algorithm runs on.
+ /// \tparam LEN The type of the length map. The default
+ /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap".
+#ifdef DOXYGEN
+ template
+#else
+ template < typename GR,
+ typename LEN = typename GR::template ArcMap,
+ typename TR = HowardDefaultTraits >
+#endif
+ class Howard
+ {
+ public:
+
+ /// The type of the digraph
+ typedef typename TR::Digraph Digraph;
+ /// The type of the length map
+ typedef typename TR::LengthMap LengthMap;
+ /// The type of the arc lengths
+ typedef typename TR::Value Value;
+
+ /// \brief The large value type
+ ///
+ /// The large value type used for internal computations.
+ /// Using the \ref HowardDefaultTraits "default traits class",
+ /// it is \c long \c long if the \c Value type is integer,
+ /// otherwise it is \c double.
+ typedef typename TR::LargeValue LargeValue;
+
+ /// The tolerance type
+ typedef typename TR::Tolerance Tolerance;
+
+ /// \brief The path type of the found cycles
+ ///
+ /// The path type of the found cycles.
+ /// Using the \ref HowardDefaultTraits "default traits class",
+ /// it is \ref lemon::Path "Path".
+ typedef typename TR::Path Path;
+
+ /// The \ref HowardDefaultTraits "traits class" of the algorithm
+ typedef TR Traits;
+
+ private:
+
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+ // The digraph the algorithm runs on
+ const Digraph &_gr;
+ // The length of the arcs
+ const LengthMap &_length;
+
+ // Data for the found cycles
+ bool _curr_found, _best_found;
+ LargeValue _curr_length, _best_length;
+ int _curr_size, _best_size;
+ Node _curr_node, _best_node;
+
+ Path *_cycle_path;
+ bool _local_path;
+
+ // Internal data used by the algorithm
+ typename Digraph::template NodeMap _policy;
+ typename Digraph::template NodeMap _reached;
+ typename Digraph::template NodeMap _level;
+ typename Digraph::template NodeMap _dist;
+
+ // Data for storing the strongly connected components
+ int _comp_num;
+ typename Digraph::template NodeMap _comp;
+ std::vector > _comp_nodes;
+ std::vector* _nodes;
+ typename Digraph::template NodeMap > _in_arcs;
+
+ // Queue used for BFS search
+ std::vector _queue;
+ int _qfront, _qback;
+
+ Tolerance _tolerance;
+
+ // Infinite constant
+ const LargeValue INF;
+
+ public:
+
+ /// \name Named Template Parameters
+ /// @{
+
+ template
+ struct SetLargeValueTraits : public Traits {
+ typedef T LargeValue;
+ typedef lemon::Tolerance Tolerance;
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// \c LargeValue type.
+ ///
+ /// \ref named-templ-param "Named parameter" for setting \c LargeValue
+ /// type. It is used for internal computations in the algorithm.
+ template
+ struct SetLargeValue
+ : public Howard > {
+ typedef Howard > Create;
+ };
+
+ template
+ struct SetPathTraits : public Traits {
+ typedef T Path;
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// \c %Path type.
+ ///
+ /// \ref named-templ-param "Named parameter" for setting the \c %Path
+ /// type of the found cycles.
+ /// It must conform to the \ref lemon::concepts::Path "Path" concept
+ /// and it must have an \c addBack() function.
+ template
+ struct SetPath
+ : public Howard > {
+ typedef Howard > Create;
+ };
+
+ /// @}
+
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// The constructor of the class.
+ ///
+ /// \param digraph The digraph the algorithm runs on.
+ /// \param length The lengths (costs) of the arcs.
+ Howard( const Digraph &digraph,
+ const LengthMap &length ) :
+ _gr(digraph), _length(length), _best_found(false),
+ _best_length(0), _best_size(1), _cycle_path(NULL), _local_path(false),
+ _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),
+ _comp(digraph), _in_arcs(digraph),
+ INF(std::numeric_limits::has_infinity ?
+ std::numeric_limits::infinity() :
+ std::numeric_limits::max())
+ {}
+
+ /// Destructor.
+ ~Howard() {
+ if (_local_path) delete _cycle_path;
+ }
+
+ /// \brief Set the path structure for storing the found cycle.
+ ///
+ /// This function sets an external path structure for storing the
+ /// found cycle.
+ ///
+ /// If you don't call this function before calling \ref run() or
+ /// \ref findMinMean(), it will allocate a local \ref Path "path"
+ /// structure. The destuctor deallocates this automatically
+ /// allocated object, of course.
+ ///
+ /// \note The algorithm calls only the \ref lemon::Path::addBack()
+ /// "addBack()" function of the given path structure.
+ ///
+ /// \return (*this)
+ Howard& cycle(Path &path) {
+ if (_local_path) {
+ delete _cycle_path;
+ _local_path = false;
+ }
+ _cycle_path = &path;
+ return *this;
+ }
+
+ /// \brief Set the tolerance used by the algorithm.
+ ///
+ /// This function sets the tolerance object used by the algorithm.
+ ///
+ /// \return (*this)
+ Howard& tolerance(const Tolerance& tolerance) {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ /// \brief Return a const reference to the tolerance.
+ ///
+ /// This function returns a const reference to the tolerance object
+ /// used by the algorithm.
+ const Tolerance& tolerance() const {
+ return _tolerance;
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to call the \ref run()
+ /// function.\n
+ /// If you only need the minimum mean length, you may call
+ /// \ref findMinMean().
+
+ /// @{
+
+ /// \brief Run the algorithm.
+ ///
+ /// This function runs the algorithm.
+ /// It can be called more than once (e.g. if the underlying digraph
+ /// and/or the arc lengths have been modified).
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ ///
+ /// \note mmc.run() is just a shortcut of the following code.
+ /// \code
+ /// return mmc.findMinMean() && mmc.findCycle();
+ /// \endcode
+ bool run() {
+ return findMinMean() && findCycle();
+ }
+
+ /// \brief Find the minimum cycle mean.
+ ///
+ /// This function finds the minimum mean length of the directed
+ /// cycles in the digraph.
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ bool findMinMean() {
+ // Initialize and find strongly connected components
+ init();
+ findComponents();
+
+ // Find the minimum cycle mean in the components
+ for (int comp = 0; comp < _comp_num; ++comp) {
+ // Find the minimum mean cycle in the current component
+ if (!buildPolicyGraph(comp)) continue;
+ while (true) {
+ findPolicyCycle();
+ if (!computeNodeDistances()) break;
+ }
+ // Update the best cycle (global minimum mean cycle)
+ if ( _curr_found && (!_best_found ||
+ _curr_length * _best_size < _best_length * _curr_size) ) {
+ _best_found = true;
+ _best_length = _curr_length;
+ _best_size = _curr_size;
+ _best_node = _curr_node;
+ }
+ }
+ return _best_found;
+ }
+
+ /// \brief Find a minimum mean directed cycle.
+ ///
+ /// This function finds a directed cycle of minimum mean length
+ /// in the digraph using the data computed by findMinMean().
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ ///
+ /// \pre \ref findMinMean() must be called before using this function.
+ bool findCycle() {
+ if (!_best_found) return false;
+ _cycle_path->addBack(_policy[_best_node]);
+ for ( Node v = _best_node;
+ (v = _gr.target(_policy[v])) != _best_node; ) {
+ _cycle_path->addBack(_policy[v]);
+ }
+ return true;
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The results of the algorithm can be obtained using these
+ /// functions.\n
+ /// The algorithm should be executed before using them.
+
+ /// @{
+
+ /// \brief Return the total length of the found cycle.
+ ///
+ /// This function returns the total length of the found cycle.
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ LargeValue cycleLength() const {
+ return _best_length;
+ }
+
+ /// \brief Return the number of arcs on the found cycle.
+ ///
+ /// This function returns the number of arcs on the found cycle.
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ int cycleArcNum() const {
+ return _best_size;
+ }
+
+ /// \brief Return the mean length of the found cycle.
+ ///
+ /// This function returns the mean length of the found cycle.
+ ///
+ /// \note alg.cycleMean() is just a shortcut of the
+ /// following code.
+ /// \code
+ /// return static_cast(alg.cycleLength()) / alg.cycleArcNum();
+ /// \endcode
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ double cycleMean() const {
+ return static_cast(_best_length) / _best_size;
+ }
+
+ /// \brief Return the found cycle.
+ ///
+ /// This function returns a const reference to the path structure
+ /// storing the found cycle.
+ ///
+ /// \pre \ref run() or \ref findCycle() must be called before using
+ /// this function.
+ const Path& cycle() const {
+ return *_cycle_path;
+ }
+
+ ///@}
+
+ private:
+
+ // Initialize
+ void init() {
+ if (!_cycle_path) {
+ _local_path = true;
+ _cycle_path = new Path;
+ }
+ _queue.resize(countNodes(_gr));
+ _best_found = false;
+ _best_length = 0;
+ _best_size = 1;
+ _cycle_path->clear();
+ }
+
+ // Find strongly connected components and initialize _comp_nodes
+ // and _in_arcs
+ void findComponents() {
+ _comp_num = stronglyConnectedComponents(_gr, _comp);
+ _comp_nodes.resize(_comp_num);
+ if (_comp_num == 1) {
+ _comp_nodes[0].clear();
+ for (NodeIt n(_gr); n != INVALID; ++n) {
+ _comp_nodes[0].push_back(n);
+ _in_arcs[n].clear();
+ for (InArcIt a(_gr, n); a != INVALID; ++a) {
+ _in_arcs[n].push_back(a);
+ }
+ }
+ } else {
+ for (int i = 0; i < _comp_num; ++i)
+ _comp_nodes[i].clear();
+ for (NodeIt n(_gr); n != INVALID; ++n) {
+ int k = _comp[n];
+ _comp_nodes[k].push_back(n);
+ _in_arcs[n].clear();
+ for (InArcIt a(_gr, n); a != INVALID; ++a) {
+ if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);
+ }
+ }
+ }
+ }
+
+ // Build the policy graph in the given strongly connected component
+ // (the out-degree of every node is 1)
+ bool buildPolicyGraph(int comp) {
+ _nodes = &(_comp_nodes[comp]);
+ if (_nodes->size() < 1 ||
+ (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {
+ return false;
+ }
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ _dist[(*_nodes)[i]] = INF;
+ }
+ Node u, v;
+ Arc e;
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ v = (*_nodes)[i];
+ for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+ e = _in_arcs[v][j];
+ u = _gr.source(e);
+ if (_length[e] < _dist[u]) {
+ _dist[u] = _length[e];
+ _policy[u] = e;
+ }
+ }
+ }
+ return true;
+ }
+
+ // Find the minimum mean cycle in the policy graph
+ void findPolicyCycle() {
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ _level[(*_nodes)[i]] = -1;
+ }
+ LargeValue clength;
+ int csize;
+ Node u, v;
+ _curr_found = false;
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ u = (*_nodes)[i];
+ if (_level[u] >= 0) continue;
+ for (; _level[u] < 0; u = _gr.target(_policy[u])) {
+ _level[u] = i;
+ }
+ if (_level[u] == i) {
+ // A cycle is found
+ clength = _length[_policy[u]];
+ csize = 1;
+ for (v = u; (v = _gr.target(_policy[v])) != u; ) {
+ clength += _length[_policy[v]];
+ ++csize;
+ }
+ if ( !_curr_found ||
+ (clength * _curr_size < _curr_length * csize) ) {
+ _curr_found = true;
+ _curr_length = clength;
+ _curr_size = csize;
+ _curr_node = u;
+ }
+ }
+ }
+ }
+
+ // Contract the policy graph and compute node distances
+ bool computeNodeDistances() {
+ // Find the component of the main cycle and compute node distances
+ // using reverse BFS
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ _reached[(*_nodes)[i]] = false;
+ }
+ _qfront = _qback = 0;
+ _queue[0] = _curr_node;
+ _reached[_curr_node] = true;
+ _dist[_curr_node] = 0;
+ Node u, v;
+ Arc e;
+ while (_qfront <= _qback) {
+ v = _queue[_qfront++];
+ for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+ e = _in_arcs[v][j];
+ u = _gr.source(e);
+ if (_policy[u] == e && !_reached[u]) {
+ _reached[u] = true;
+ _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
+ _queue[++_qback] = u;
+ }
+ }
+ }
+
+ // Connect all other nodes to this component and compute node
+ // distances using reverse BFS
+ _qfront = 0;
+ while (_qback < int(_nodes->size())-1) {
+ v = _queue[_qfront++];
+ for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+ e = _in_arcs[v][j];
+ u = _gr.source(e);
+ if (!_reached[u]) {
+ _reached[u] = true;
+ _policy[u] = e;
+ _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
+ _queue[++_qback] = u;
+ }
+ }
+ }
+
+ // Improve node distances
+ bool improved = false;
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ v = (*_nodes)[i];
+ for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+ e = _in_arcs[v][j];
+ u = _gr.source(e);
+ LargeValue delta = _dist[v] + _length[e] * _curr_size - _curr_length;
+ if (_tolerance.less(delta, _dist[u])) {
+ _dist[u] = delta;
+ _policy[u] = e;
+ improved = true;
+ }
+ }
+ }
+ return improved;
+ }
+
+ }; //class Howard
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_HOWARD_H
Index: lemon/karp.h
===================================================================
--- lemon/karp.h (revision 816)
+++ lemon/karp.h (revision 816)
@@ -0,0 +1,581 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_KARP_H
+#define LEMON_KARP_H
+
+/// \ingroup min_mean_cycle
+///
+/// \file
+/// \brief Karp's algorithm for finding a minimum mean cycle.
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace lemon {
+
+ /// \brief Default traits class of Karp algorithm.
+ ///
+ /// Default traits class of Karp algorithm.
+ /// \tparam GR The type of the digraph.
+ /// \tparam LEN The type of the length map.
+ /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
+#ifdef DOXYGEN
+ template
+#else
+ template ::is_integer>
+#endif
+ struct KarpDefaultTraits
+ {
+ /// The type of the digraph
+ typedef GR Digraph;
+ /// The type of the length map
+ typedef LEN LengthMap;
+ /// The type of the arc lengths
+ typedef typename LengthMap::Value Value;
+
+ /// \brief The large value type used for internal computations
+ ///
+ /// The large value type used for internal computations.
+ /// It is \c long \c long if the \c Value type is integer,
+ /// otherwise it is \c double.
+ /// \c Value must be convertible to \c LargeValue.
+ typedef double LargeValue;
+
+ /// The tolerance type used for internal computations
+ typedef lemon::Tolerance Tolerance;
+
+ /// \brief The path type of the found cycles
+ ///
+ /// The path type of the found cycles.
+ /// It must conform to the \ref lemon::concepts::Path "Path" concept
+ /// and it must have an \c addBack() function.
+ typedef lemon::Path Path;
+ };
+
+ // Default traits class for integer value types
+ template
+ struct KarpDefaultTraits
+ {
+ typedef GR Digraph;
+ typedef LEN LengthMap;
+ typedef typename LengthMap::Value Value;
+#ifdef LEMON_HAVE_LONG_LONG
+ typedef long long LargeValue;
+#else
+ typedef long LargeValue;
+#endif
+ typedef lemon::Tolerance Tolerance;
+ typedef lemon::Path Path;
+ };
+
+
+ /// \addtogroup min_mean_cycle
+ /// @{
+
+ /// \brief Implementation of Karp's algorithm for finding a minimum
+ /// mean cycle.
+ ///
+ /// This class implements Karp's algorithm for finding a directed
+ /// cycle of minimum mean length (cost) in a digraph.
+ /// It runs in time O(ne) and uses space O(n2+e).
+ ///
+ /// \tparam GR The type of the digraph the algorithm runs on.
+ /// \tparam LEN The type of the length map. The default
+ /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap".
+#ifdef DOXYGEN
+ template
+#else
+ template < typename GR,
+ typename LEN = typename GR::template ArcMap,
+ typename TR = KarpDefaultTraits >
+#endif
+ class Karp
+ {
+ public:
+
+ /// The type of the digraph
+ typedef typename TR::Digraph Digraph;
+ /// The type of the length map
+ typedef typename TR::LengthMap LengthMap;
+ /// The type of the arc lengths
+ typedef typename TR::Value Value;
+
+ /// \brief The large value type
+ ///
+ /// The large value type used for internal computations.
+ /// Using the \ref KarpDefaultTraits "default traits class",
+ /// it is \c long \c long if the \c Value type is integer,
+ /// otherwise it is \c double.
+ typedef typename TR::LargeValue LargeValue;
+
+ /// The tolerance type
+ typedef typename TR::Tolerance Tolerance;
+
+ /// \brief The path type of the found cycles
+ ///
+ /// The path type of the found cycles.
+ /// Using the \ref KarpDefaultTraits "default traits class",
+ /// it is \ref lemon::Path "Path".
+ typedef typename TR::Path Path;
+
+ /// The \ref KarpDefaultTraits "traits class" of the algorithm
+ typedef TR Traits;
+
+ private:
+
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+ // Data sturcture for path data
+ struct PathData
+ {
+ LargeValue dist;
+ Arc pred;
+ PathData(LargeValue d, Arc p = INVALID) :
+ dist(d), pred(p) {}
+ };
+
+ typedef typename Digraph::template NodeMap >
+ PathDataNodeMap;
+
+ private:
+
+ // The digraph the algorithm runs on
+ const Digraph &_gr;
+ // The length of the arcs
+ const LengthMap &_length;
+
+ // Data for storing the strongly connected components
+ int _comp_num;
+ typename Digraph::template NodeMap _comp;
+ std::vector > _comp_nodes;
+ std::vector* _nodes;
+ typename Digraph::template NodeMap > _out_arcs;
+
+ // Data for the found cycle
+ LargeValue _cycle_length;
+ int _cycle_size;
+ Node _cycle_node;
+
+ Path *_cycle_path;
+ bool _local_path;
+
+ // Node map for storing path data
+ PathDataNodeMap _data;
+ // The processed nodes in the last round
+ std::vector _process;
+
+ Tolerance _tolerance;
+
+ // Infinite constant
+ const LargeValue INF;
+
+ public:
+
+ /// \name Named Template Parameters
+ /// @{
+
+ template
+ struct SetLargeValueTraits : public Traits {
+ typedef T LargeValue;
+ typedef lemon::Tolerance Tolerance;
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// \c LargeValue type.
+ ///
+ /// \ref named-templ-param "Named parameter" for setting \c LargeValue
+ /// type. It is used for internal computations in the algorithm.
+ template
+ struct SetLargeValue
+ : public Karp > {
+ typedef Karp > Create;
+ };
+
+ template
+ struct SetPathTraits : public Traits {
+ typedef T Path;
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// \c %Path type.
+ ///
+ /// \ref named-templ-param "Named parameter" for setting the \c %Path
+ /// type of the found cycles.
+ /// It must conform to the \ref lemon::concepts::Path "Path" concept
+ /// and it must have an \c addFront() function.
+ template
+ struct SetPath
+ : public Karp > {
+ typedef Karp > Create;
+ };
+
+ /// @}
+
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// The constructor of the class.
+ ///
+ /// \param digraph The digraph the algorithm runs on.
+ /// \param length The lengths (costs) of the arcs.
+ Karp( const Digraph &digraph,
+ const LengthMap &length ) :
+ _gr(digraph), _length(length), _comp(digraph), _out_arcs(digraph),
+ _cycle_length(0), _cycle_size(1), _cycle_node(INVALID),
+ _cycle_path(NULL), _local_path(false), _data(digraph),
+ INF(std::numeric_limits::has_infinity ?
+ std::numeric_limits::infinity() :
+ std::numeric_limits::max())
+ {}
+
+ /// Destructor.
+ ~Karp() {
+ if (_local_path) delete _cycle_path;
+ }
+
+ /// \brief Set the path structure for storing the found cycle.
+ ///
+ /// This function sets an external path structure for storing the
+ /// found cycle.
+ ///
+ /// If you don't call this function before calling \ref run() or
+ /// \ref findMinMean(), it will allocate a local \ref Path "path"
+ /// structure. The destuctor deallocates this automatically
+ /// allocated object, of course.
+ ///
+ /// \note The algorithm calls only the \ref lemon::Path::addFront()
+ /// "addFront()" function of the given path structure.
+ ///
+ /// \return (*this)
+ Karp& cycle(Path &path) {
+ if (_local_path) {
+ delete _cycle_path;
+ _local_path = false;
+ }
+ _cycle_path = &path;
+ return *this;
+ }
+
+ /// \brief Set the tolerance used by the algorithm.
+ ///
+ /// This function sets the tolerance object used by the algorithm.
+ ///
+ /// \return (*this)
+ Karp& tolerance(const Tolerance& tolerance) {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ /// \brief Return a const reference to the tolerance.
+ ///
+ /// This function returns a const reference to the tolerance object
+ /// used by the algorithm.
+ const Tolerance& tolerance() const {
+ return _tolerance;
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to call the \ref run()
+ /// function.\n
+ /// If you only need the minimum mean length, you may call
+ /// \ref findMinMean().
+
+ /// @{
+
+ /// \brief Run the algorithm.
+ ///
+ /// This function runs the algorithm.
+ /// It can be called more than once (e.g. if the underlying digraph
+ /// and/or the arc lengths have been modified).
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ ///
+ /// \note mmc.run() is just a shortcut of the following code.
+ /// \code
+ /// return mmc.findMinMean() && mmc.findCycle();
+ /// \endcode
+ bool run() {
+ return findMinMean() && findCycle();
+ }
+
+ /// \brief Find the minimum cycle mean.
+ ///
+ /// This function finds the minimum mean length of the directed
+ /// cycles in the digraph.
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ bool findMinMean() {
+ // Initialization and find strongly connected components
+ init();
+ findComponents();
+
+ // Find the minimum cycle mean in the components
+ for (int comp = 0; comp < _comp_num; ++comp) {
+ if (!initComponent(comp)) continue;
+ processRounds();
+ updateMinMean();
+ }
+ return (_cycle_node != INVALID);
+ }
+
+ /// \brief Find a minimum mean directed cycle.
+ ///
+ /// This function finds a directed cycle of minimum mean length
+ /// in the digraph using the data computed by findMinMean().
+ ///
+ /// \return \c true if a directed cycle exists in the digraph.
+ ///
+ /// \pre \ref findMinMean() must be called before using this function.
+ bool findCycle() {
+ if (_cycle_node == INVALID) return false;
+ IntNodeMap reached(_gr, -1);
+ int r = _data[_cycle_node].size();
+ Node u = _cycle_node;
+ while (reached[u] < 0) {
+ reached[u] = --r;
+ u = _gr.source(_data[u][r].pred);
+ }
+ r = reached[u];
+ Arc e = _data[u][r].pred;
+ _cycle_path->addFront(e);
+ _cycle_length = _length[e];
+ _cycle_size = 1;
+ Node v;
+ while ((v = _gr.source(e)) != u) {
+ e = _data[v][--r].pred;
+ _cycle_path->addFront(e);
+ _cycle_length += _length[e];
+ ++_cycle_size;
+ }
+ return true;
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The results of the algorithm can be obtained using these
+ /// functions.\n
+ /// The algorithm should be executed before using them.
+
+ /// @{
+
+ /// \brief Return the total length of the found cycle.
+ ///
+ /// This function returns the total length of the found cycle.
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ LargeValue cycleLength() const {
+ return _cycle_length;
+ }
+
+ /// \brief Return the number of arcs on the found cycle.
+ ///
+ /// This function returns the number of arcs on the found cycle.
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ int cycleArcNum() const {
+ return _cycle_size;
+ }
+
+ /// \brief Return the mean length of the found cycle.
+ ///
+ /// This function returns the mean length of the found cycle.
+ ///
+ /// \note alg.cycleMean() is just a shortcut of the
+ /// following code.
+ /// \code
+ /// return static_cast(alg.cycleLength()) / alg.cycleArcNum();
+ /// \endcode
+ ///
+ /// \pre \ref run() or \ref findMinMean() must be called before
+ /// using this function.
+ double cycleMean() const {
+ return static_cast(_cycle_length) / _cycle_size;
+ }
+
+ /// \brief Return the found cycle.
+ ///
+ /// This function returns a const reference to the path structure
+ /// storing the found cycle.
+ ///
+ /// \pre \ref run() or \ref findCycle() must be called before using
+ /// this function.
+ const Path& cycle() const {
+ return *_cycle_path;
+ }
+
+ ///@}
+
+ private:
+
+ // Initialization
+ void init() {
+ if (!_cycle_path) {
+ _local_path = true;
+ _cycle_path = new Path;
+ }
+ _cycle_path->clear();
+ _cycle_length = 0;
+ _cycle_size = 1;
+ _cycle_node = INVALID;
+ for (NodeIt u(_gr); u != INVALID; ++u)
+ _data[u].clear();
+ }
+
+ // Find strongly connected components and initialize _comp_nodes
+ // and _out_arcs
+ void findComponents() {
+ _comp_num = stronglyConnectedComponents(_gr, _comp);
+ _comp_nodes.resize(_comp_num);
+ if (_comp_num == 1) {
+ _comp_nodes[0].clear();
+ for (NodeIt n(_gr); n != INVALID; ++n) {
+ _comp_nodes[0].push_back(n);
+ _out_arcs[n].clear();
+ for (OutArcIt a(_gr, n); a != INVALID; ++a) {
+ _out_arcs[n].push_back(a);
+ }
+ }
+ } else {
+ for (int i = 0; i < _comp_num; ++i)
+ _comp_nodes[i].clear();
+ for (NodeIt n(_gr); n != INVALID; ++n) {
+ int k = _comp[n];
+ _comp_nodes[k].push_back(n);
+ _out_arcs[n].clear();
+ for (OutArcIt a(_gr, n); a != INVALID; ++a) {
+ if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
+ }
+ }
+ }
+ }
+
+ // Initialize path data for the current component
+ bool initComponent(int comp) {
+ _nodes = &(_comp_nodes[comp]);
+ int n = _nodes->size();
+ if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
+ return false;
+ }
+ for (int i = 0; i < n; ++i) {
+ _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
+ }
+ return true;
+ }
+
+ // Process all rounds of computing path data for the current component.
+ // _data[v][k] is the length of a shortest directed walk from the root
+ // node to node v containing exactly k arcs.
+ void processRounds() {
+ Node start = (*_nodes)[0];
+ _data[start][0] = PathData(0);
+ _process.clear();
+ _process.push_back(start);
+
+ int k, n = _nodes->size();
+ for (k = 1; k <= n && int(_process.size()) < n; ++k) {
+ processNextBuildRound(k);
+ }
+ for ( ; k <= n; ++k) {
+ processNextFullRound(k);
+ }
+ }
+
+ // Process one round and rebuild _process
+ void processNextBuildRound(int k) {
+ std::vector next;
+ Node u, v;
+ Arc e;
+ LargeValue d;
+ for (int i = 0; i < int(_process.size()); ++i) {
+ u = _process[i];
+ for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
+ e = _out_arcs[u][j];
+ v = _gr.target(e);
+ d = _data[u][k-1].dist + _length[e];
+ if (_tolerance.less(d, _data[v][k].dist)) {
+ if (_data[v][k].dist == INF) next.push_back(v);
+ _data[v][k] = PathData(d, e);
+ }
+ }
+ }
+ _process.swap(next);
+ }
+
+ // Process one round using _nodes instead of _process
+ void processNextFullRound(int k) {
+ Node u, v;
+ Arc e;
+ LargeValue d;
+ for (int i = 0; i < int(_nodes->size()); ++i) {
+ u = (*_nodes)[i];
+ for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
+ e = _out_arcs[u][j];
+ v = _gr.target(e);
+ d = _data[u][k-1].dist + _length[e];
+ if (_tolerance.less(d, _data[v][k].dist)) {
+ _data[v][k] = PathData(d, e);
+ }
+ }
+ }
+ }
+
+ // Update the minimum cycle mean
+ void updateMinMean() {
+ int n = _nodes->size();
+ for (int i = 0; i < n; ++i) {
+ Node u = (*_nodes)[i];
+ if (_data[u][n].dist == INF) continue;
+ LargeValue length, max_length = 0;
+ int size, max_size = 1;
+ bool found_curr = false;
+ for (int k = 0; k < n; ++k) {
+ if (_data[u][k].dist == INF) continue;
+ length = _data[u][n].dist - _data[u][k].dist;
+ size = n - k;
+ if (!found_curr || length * max_size > max_length * size) {
+ found_curr = true;
+ max_length = length;
+ max_size = size;
+ }
+ }
+ if ( found_curr && (_cycle_node == INVALID ||
+ max_length * _cycle_size < _cycle_length * max_size) ) {
+ _cycle_length = max_length;
+ _cycle_size = max_size;
+ _cycle_node = u;
+ }
+ }
+ }
+
+ }; //class Karp
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_KARP_H
Index: test/CMakeLists.txt
===================================================================
--- test/CMakeLists.txt (revision 745)
+++ test/CMakeLists.txt (revision 817)
@@ -33,4 +33,5 @@
min_cost_arborescence_test
min_cost_flow_test
+ min_mean_cycle_test
path_test
preflow_test
Index: test/Makefile.am
===================================================================
--- test/Makefile.am (revision 745)
+++ test/Makefile.am (revision 817)
@@ -31,4 +31,5 @@
test/min_cost_arborescence_test \
test/min_cost_flow_test \
+ test/min_mean_cycle_test \
test/path_test \
test/preflow_test \
@@ -79,4 +80,5 @@
test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
+test_min_mean_cycle_test_SOURCES = test/min_mean_cycle_test.cc
test_path_test_SOURCES = test/path_test.cc
test_preflow_test_SOURCES = test/preflow_test.cc
Index: test/min_mean_cycle_test.cc
===================================================================
--- test/min_mean_cycle_test.cc (revision 816)
+++ test/min_mean_cycle_test.cc (revision 816)
@@ -0,0 +1,216 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2009
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include "test_tools.h"
+
+using namespace lemon;
+
+char test_lgf[] =
+ "@nodes\n"
+ "label\n"
+ "1\n"
+ "2\n"
+ "3\n"
+ "4\n"
+ "5\n"
+ "6\n"
+ "7\n"
+ "@arcs\n"
+ " len1 len2 len3 len4 c1 c2 c3 c4\n"
+ "1 2 1 1 1 1 0 0 0 0\n"
+ "2 4 5 5 5 5 1 0 0 0\n"
+ "2 3 8 8 8 8 0 0 0 0\n"
+ "3 2 -2 0 0 0 1 0 0 0\n"
+ "3 4 4 4 4 4 0 0 0 0\n"
+ "3 7 -4 -4 -4 -4 0 0 0 0\n"
+ "4 1 2 2 2 2 0 0 0 0\n"
+ "4 3 3 3 3 3 1 0 0 0\n"
+ "4 4 3 3 0 0 0 0 1 0\n"
+ "5 2 4 4 4 4 0 0 0 0\n"
+ "5 6 3 3 3 3 0 1 0 0\n"
+ "6 5 2 2 2 2 0 1 0 0\n"
+ "6 4 -1 -1 -1 -1 0 0 0 0\n"
+ "6 7 1 1 1 1 0 0 0 0\n"
+ "7 7 4 4 4 -1 0 0 0 1\n";
+
+
+// Check the interface of an MMC algorithm
+template
+struct MmcClassConcept
+{
+ template
+ struct Constraints {
+ void constraints() {
+ const Constraints& me = *this;
+
+ typedef typename MMC
+ ::template SetPath >
+ ::template SetLargeValue
+ ::Create MmcAlg;
+ MmcAlg mmc(me.g, me.length);
+ const MmcAlg& const_mmc = mmc;
+
+ typename MmcAlg::Tolerance tol = const_mmc.tolerance();
+ mmc.tolerance(tol);
+
+ b = mmc.cycle(p).run();
+ b = mmc.findMinMean();
+ b = mmc.findCycle();
+
+ v = const_mmc.cycleLength();
+ i = const_mmc.cycleArcNum();
+ d = const_mmc.cycleMean();
+ p = const_mmc.cycle();
+ }
+
+ typedef concepts::ReadMap LM;
+
+ GR g;
+ LM length;
+ ListPath p;
+ Value v;
+ int i;
+ double d;
+ bool b;
+ };
+};
+
+// Perform a test with the given parameters
+template
+void checkMmcAlg(const SmartDigraph& gr,
+ const SmartDigraph::ArcMap& lm,
+ const SmartDigraph::ArcMap& cm,
+ int length, int size) {
+ MMC alg(gr, lm);
+ alg.findMinMean();
+ check(alg.cycleMean() == static_cast(length) / size,
+ "Wrong cycle mean");
+ alg.findCycle();
+ check(alg.cycleLength() == length && alg.cycleArcNum() == size,
+ "Wrong path");
+ SmartDigraph::ArcMap cycle(gr, 0);
+ for (typename MMC::Path::ArcIt a(alg.cycle()); a != INVALID; ++a) {
+ ++cycle[a];
+ }
+ for (SmartDigraph::ArcIt a(gr); a != INVALID; ++a) {
+ check(cm[a] == cycle[a], "Wrong path");
+ }
+}
+
+// Class for comparing types
+template
+struct IsSameType {
+ static const int result = 0;
+};
+
+template
+struct IsSameType {
+ static const int result = 1;
+};
+
+
+int main() {
+ #ifdef LEMON_HAVE_LONG_LONG
+ typedef long long long_int;
+ #else
+ typedef long long_int;
+ #endif
+
+ // Check the interface
+ {
+ typedef concepts::Digraph GR;
+
+ // Karp
+ checkConcept< MmcClassConcept,
+ Karp > >();
+ checkConcept< MmcClassConcept,
+ Karp > >();
+
+ // HartmannOrlin
+ checkConcept< MmcClassConcept,
+ HartmannOrlin > >();
+ checkConcept< MmcClassConcept,
+ HartmannOrlin > >();
+
+ // Howard
+ checkConcept< MmcClassConcept,
+ Howard > >();
+ checkConcept< MmcClassConcept,
+ Howard > >();
+
+ if (IsSameType >::LargeValue,
+ long_int>::result == 0) check(false, "Wrong LargeValue type");
+ if (IsSameType >::LargeValue,
+ double>::result == 0) check(false, "Wrong LargeValue type");
+ }
+
+ // Run various tests
+ {
+ typedef SmartDigraph GR;
+ DIGRAPH_TYPEDEFS(GR);
+
+ GR gr;
+ IntArcMap l1(gr), l2(gr), l3(gr), l4(gr);
+ IntArcMap c1(gr), c2(gr), c3(gr), c4(gr);
+
+ std::istringstream input(test_lgf);
+ digraphReader(gr, input).
+ arcMap("len1", l1).
+ arcMap("len2", l2).
+ arcMap("len3", l3).
+ arcMap("len4", l4).
+ arcMap("c1", c1).
+ arcMap("c2", c2).
+ arcMap("c3", c3).
+ arcMap("c4", c4).
+ run();
+
+ // Karp
+ checkMmcAlg >(gr, l1, c1, 6, 3);
+ checkMmcAlg >(gr, l2, c2, 5, 2);
+ checkMmcAlg >(gr, l3, c3, 0, 1);
+ checkMmcAlg >(gr, l4, c4, -1, 1);
+
+ // HartmannOrlin
+ checkMmcAlg >(gr, l1, c1, 6, 3);
+ checkMmcAlg >(gr, l2, c2, 5, 2);
+ checkMmcAlg >(gr, l3, c3, 0, 1);
+ checkMmcAlg >(gr, l4, c4, -1, 1);
+
+ // Howard
+ checkMmcAlg >(gr, l1, c1, 6, 3);
+ checkMmcAlg >(gr, l2, c2, 5, 2);
+ checkMmcAlg >(gr, l3, c3, 0, 1);
+ checkMmcAlg >(gr, l4, c4, -1, 1);
+ }
+
+ return 0;
+}
Index: test/test_tools.h
===================================================================
--- test/test_tools.h (revision 463)
+++ test/test_tools.h (revision 810)
@@ -38,9 +38,13 @@
///print something like this (and then exits).
///\verbatim file_name.cc:123: error: This is obviously false. \endverbatim
-#define check(rc, msg) \
- if(!(rc)) { \
- std::cerr << __FILE__ ":" << __LINE__ << ": error: " << msg << std::endl; \
- abort(); \
- } else { } \
+#define check(rc, msg) \
+ { \
+ if(!(rc)) { \
+ std::cerr << __FILE__ ":" << __LINE__ << ": error: " \
+ << msg << std::endl; \
+ abort(); \
+ } else { } \
+ } \
+
#endif