Reimplemented MinMeanCycle to be much more efficient.
authorkpeter
Sun, 13 Jan 2008 10:26:55 +0000
changeset 2555a84e52e99f57
parent 2554 1775aaa02ac4
child 2556 74c2c81055e1
Reimplemented MinMeanCycle to be much more efficient.
The new version implements Howard's algorithm instead of Karp's algorithm and
it is at least 10-20 times faster on all the 40-50 random graphs we have tested.
lemon/min_mean_cycle.h
     1.1 --- a/lemon/min_mean_cycle.h	Mon Jan 07 17:07:40 2008 +0000
     1.2 +++ b/lemon/min_mean_cycle.h	Sun Jan 13 10:26:55 2008 +0000
     1.3 @@ -22,224 +22,242 @@
     1.4  /// \ingroup shortest_path
     1.5  ///
     1.6  /// \file
     1.7 -/// \brief Karp's algorithm for finding a minimum mean (directed) cycle.
     1.8 +/// \brief Howard's algorithm for finding a minimum mean directed cycle.
     1.9  
    1.10  #include <vector>
    1.11  #include <lemon/graph_utils.h>
    1.12  #include <lemon/topology.h>
    1.13  #include <lemon/path.h>
    1.14 +#include <lemon/tolerance.h>
    1.15  
    1.16  namespace lemon {
    1.17  
    1.18    /// \addtogroup shortest_path
    1.19    /// @{
    1.20  
    1.21 -  /// \brief Implementation of Karp's algorithm for finding a
    1.22 -  /// minimum mean (directed) cycle.
    1.23 +  /// \brief Implementation of Howard's algorithm for finding a minimum
    1.24 +  /// mean directed cycle.
    1.25    ///
    1.26 -  /// The \ref lemon::MinMeanCycle "MinMeanCycle" implements Karp's
    1.27 -  /// algorithm for finding a minimum mean (directed) cycle.
    1.28 +  /// \ref MinMeanCycle implements Howard's algorithm for finding a
    1.29 +  /// minimum mean directed cycle.
    1.30    ///
    1.31    /// \param Graph The directed graph type the algorithm runs on.
    1.32    /// \param LengthMap The type of the length (cost) map.
    1.33    ///
    1.34 +  /// \warning \c LengthMap::Value must be convertible to \c double.
    1.35 +  ///
    1.36    /// \author Peter Kovacs
    1.37  
    1.38 -  template <typename Graph,
    1.39 -    typename LengthMap = typename Graph::template EdgeMap<int> >
    1.40 +  template < typename Graph,
    1.41 +             typename LengthMap = typename Graph::template EdgeMap<int> >
    1.42    class MinMeanCycle
    1.43    {
    1.44 -    typedef typename Graph::Node Node;
    1.45 -    typedef typename Graph::NodeIt NodeIt;
    1.46 -    typedef typename Graph::Edge Edge;
    1.47 -    typedef typename Graph::EdgeIt EdgeIt;
    1.48 -    typedef typename Graph::OutEdgeIt OutEdgeIt;
    1.49 +    GRAPH_TYPEDEFS(typename Graph);
    1.50  
    1.51      typedef typename LengthMap::Value Length;
    1.52 -
    1.53 -    typedef typename Graph::template NodeMap<int> IntNodeMap;
    1.54 -    typedef typename Graph::template NodeMap<Edge> PredNodeMap;
    1.55      typedef Path<Graph> Path;
    1.56 -    typedef std::vector<Node> NodeVector;
    1.57  
    1.58    protected:
    1.59  
    1.60 -    /// \brief Data structure for path data.
    1.61 -    struct PathData
    1.62 -    {
    1.63 -      bool found;
    1.64 -      Length dist;
    1.65 -      Edge pred;
    1.66 -      PathData(bool _found = false, Length _dist = 0) :
    1.67 -	found(_found), dist(_dist), pred(INVALID) {}
    1.68 -      PathData(bool _found, Length _dist, Edge _pred) :
    1.69 -	found(_found), dist(_dist), pred(_pred) {}
    1.70 -    };
    1.71 +    typename Graph::template NodeMap<bool> _reached;
    1.72 +    typename Graph::template NodeMap<double> _dist;
    1.73 +    typename Graph::template NodeMap<Edge> _policy;
    1.74  
    1.75 -  private:
    1.76 +    /// The directed graph the algorithm runs on.
    1.77 +    const Graph &_graph;
    1.78 +    /// The length of the edges.
    1.79 +    const LengthMap &_length;
    1.80  
    1.81 -    typedef typename Graph::template NodeMap<std::vector<PathData> >
    1.82 -      PathDataNodeMap;
    1.83 +    /// The total length of the found cycle.
    1.84 +    Length _cycle_length;
    1.85 +    /// The number of edges on the found cycle.
    1.86 +    int _cycle_size;
    1.87 +    /// The found cycle.
    1.88 +    Path *_cycle_path;
    1.89  
    1.90 -  protected:
    1.91 +    bool _local_path;
    1.92 +    bool _cycle_found;
    1.93 +    Node _cycle_node;
    1.94  
    1.95 -    /// \brief Node map for storing path data.
    1.96 -    ///
    1.97 -    /// Node map for storing path data of all nodes in the current
    1.98 -    /// component. dmap[v][k] is the length of a shortest directed walk
    1.99 -    /// to node v from the starting node containing exactly k edges.
   1.100 -    PathDataNodeMap dmap;
   1.101 +    typename Graph::template NodeMap<int> _component;
   1.102 +    int _component_num;
   1.103  
   1.104 -    /// \brief The directed graph the algorithm runs on.
   1.105 -    const Graph &graph;
   1.106 -    /// \brief The length of the edges.
   1.107 -    const LengthMap &length;
   1.108 +    std::vector<Node> _nodes;
   1.109 +    std::vector<Edge> _edges;
   1.110 +    Tolerance<double> _tolerance;
   1.111  
   1.112 -    /// \brief The total length of the found cycle.
   1.113 -    Length cycle_length;
   1.114 -    /// \brief The number of edges in the found cycle.
   1.115 -    int cycle_size;
   1.116 -    /// \brief A node for obtaining a minimum mean cycle.
   1.117 -    Node cycle_node;
   1.118 -
   1.119 -    /// \brief The found cycle.
   1.120 -    Path *cycle_path;
   1.121 -    /// \brief The algorithm uses local \ref lemon::Path "Path"
   1.122 -    /// structure to store the found cycle.
   1.123 -    bool local_path;
   1.124 -
   1.125 -    /// \brief Node map for identifying strongly connected components.
   1.126 -    IntNodeMap comp;
   1.127 -    /// \brief The number of strongly connected components.
   1.128 -    int comp_num;
   1.129 -    /// \brief Counter for identifying the current component.
   1.130 -    int comp_cnt;
   1.131 -    /// \brief Nodes of the current component.
   1.132 -    NodeVector nodes;
   1.133 -    /// \brief The processed nodes in the last round.
   1.134 -    NodeVector process;
   1.135 -
   1.136 -  public :
   1.137 +  public:
   1.138  
   1.139      /// \brief The constructor of the class.
   1.140      ///
   1.141      /// The constructor of the class.
   1.142      ///
   1.143 -    /// \param _graph The directed graph the algorithm runs on.
   1.144 -    /// \param _length The length (cost) of the edges.
   1.145 -    MinMeanCycle( const Graph &_graph,
   1.146 -		  const LengthMap &_length ) :
   1.147 -      graph(_graph), length(_length), dmap(_graph), comp(_graph),
   1.148 -      cycle_length(0), cycle_size(-1), cycle_node(INVALID),
   1.149 -      cycle_path(NULL), local_path(false)
   1.150 +    /// \param graph The directed graph the algorithm runs on.
   1.151 +    /// \param length The length (cost) of the edges.
   1.152 +    MinMeanCycle( const Graph &graph,
   1.153 +                  const LengthMap &length ) :
   1.154 +      _graph(graph), _length(length), _dist(graph), _reached(graph),
   1.155 +      _policy(graph), _component(graph), _cycle_length(0),
   1.156 +      _cycle_size(-1), _cycle_path(NULL), _local_path(false)
   1.157      { }
   1.158  
   1.159 -    /// \brief The destructor of the class.
   1.160 +    /// The destructor of the class.
   1.161      ~MinMeanCycle() {
   1.162 -      if (local_path) delete cycle_path;
   1.163 +      if (_local_path) delete _cycle_path;
   1.164      }
   1.165  
   1.166    protected:
   1.167  
   1.168 -    /// \brief Initializes the internal data structures for the current
   1.169 -    /// component.
   1.170 -    void initCurrent() {
   1.171 -      nodes.clear();
   1.172 +    // Initializes the internal data structures for the current strongly
   1.173 +    // connected component and creating the policy graph.
   1.174 +    // The policy graph can be represented by the _policy map because
   1.175 +    // the out degree of every node is 1.
   1.176 +    bool initCurrentComponent(int comp) {
   1.177        // Finding the nodes of the current component
   1.178 -      for (NodeIt v(graph); v != INVALID; ++v) {
   1.179 -	if (comp[v] == comp_cnt) nodes.push_back(v);
   1.180 +      _nodes.clear();
   1.181 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.182 +        if (_component[n] == comp) _nodes.push_back(n);
   1.183        }
   1.184 -      // Creating vectors for all nodes
   1.185 -      int n = nodes.size();
   1.186 -      for (int i = 0; i < n; ++i) {
   1.187 -	dmap[nodes[i]].resize(n + 1);
   1.188 +      if (_nodes.size() <= 1) return false;
   1.189 +      // Finding the edges of the current component
   1.190 +      _edges.clear();
   1.191 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.192 +        if ( _component[_graph.source(e)] == comp &&
   1.193 +             _component[_graph.target(e)] == comp )
   1.194 +          _edges.push_back(e);
   1.195 +      }
   1.196 +      // Initializing _reached, _dist, _policy maps
   1.197 +      for (int i = 0; i < _nodes.size(); ++i) {
   1.198 +        _reached[_nodes[i]] = false;
   1.199 +        _policy[_nodes[i]] = INVALID;
   1.200 +      }
   1.201 +      Node u; Edge e;
   1.202 +      for (int j = 0; j < _edges.size(); ++j) {
   1.203 +        e = _edges[j];
   1.204 +        u = _graph.source(e);
   1.205 +        if (!_reached[u] || _length[e] < _dist[u]) {
   1.206 +          _dist[u] = _length[e];
   1.207 +          _policy[u] = e;
   1.208 +          _reached[u] = true;
   1.209 +        }
   1.210 +      }
   1.211 +      return true;
   1.212 +    }
   1.213 +
   1.214 +    // Finds all cycles in the policy graph.
   1.215 +    // Sets _cycle_found to true if a cycle is found and sets
   1.216 +    // _cycle_length, _cycle_size, _cycle_node to represent the minimum
   1.217 +    // mean cycle in the policy graph.
   1.218 +    bool findPolicyCycles(int comp) {
   1.219 +      typename Graph::template NodeMap<int> level(_graph, -1);
   1.220 +      _cycle_found = false;
   1.221 +      Length clength;
   1.222 +      int csize;
   1.223 +      int path_cnt = 0;
   1.224 +      Node u, v;
   1.225 +      // Searching for cycles
   1.226 +      for (int i = 0; i < _nodes.size(); ++i) {
   1.227 +        if (level[_nodes[i]] < 0) {
   1.228 +          u = _nodes[i];
   1.229 +          level[u] = path_cnt;
   1.230 +          while (level[u = _graph.target(_policy[u])] < 0)
   1.231 +            level[u] = path_cnt;
   1.232 +          if (level[u] == path_cnt) {
   1.233 +            // A cycle is found
   1.234 +            clength = _length[_policy[u]];
   1.235 +            csize = 1;
   1.236 +            for (v = u; (v = _graph.target(_policy[v])) != u; ) {
   1.237 +              clength += _length[_policy[v]];
   1.238 +              ++csize;
   1.239 +            }
   1.240 +            if ( !_cycle_found ||
   1.241 +                 clength * _cycle_size < _cycle_length * csize ) {
   1.242 +              _cycle_found = true;
   1.243 +              _cycle_length = clength;
   1.244 +              _cycle_size = csize;
   1.245 +              _cycle_node = u;
   1.246 +            }
   1.247 +          }
   1.248 +          ++path_cnt;
   1.249 +        }
   1.250 +      }
   1.251 +      return _cycle_found;
   1.252 +    }
   1.253 +
   1.254 +    // Contracts the policy graph to be connected by cutting all cycles
   1.255 +    // except for the main cycle (i.e. the minimum mean cycle).
   1.256 +    void contractPolicyGraph(int comp) {
   1.257 +      // Finding the component of the main cycle using
   1.258 +      // reverse BFS search
   1.259 +      typename Graph::template NodeMap<int> found(_graph, false);
   1.260 +      std::deque<Node> queue;
   1.261 +      queue.push_back(_cycle_node);
   1.262 +      found[_cycle_node] = true;
   1.263 +      Node u, v;
   1.264 +      while (!queue.empty()) {
   1.265 +        v = queue.front(); queue.pop_front();
   1.266 +        for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
   1.267 +          u = _graph.source(e);
   1.268 +          if (_component[u] == comp && !found[u] && _policy[u] == e) {
   1.269 +            found[u] = true;
   1.270 +            queue.push_back(u);
   1.271 +          }
   1.272 +        }
   1.273 +      }
   1.274 +      // Connecting all other nodes to this component using
   1.275 +      // reverse BFS search
   1.276 +      queue.clear();
   1.277 +      for (int i = 0; i < _nodes.size(); ++i)
   1.278 +        if (found[_nodes[i]]) queue.push_back(_nodes[i]);
   1.279 +      int found_cnt = queue.size();
   1.280 +      while (found_cnt < _nodes.size() && !queue.empty()) {
   1.281 +        v = queue.front(); queue.pop_front();
   1.282 +        for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
   1.283 +          u = _graph.source(e);
   1.284 +          if (_component[u] == comp && !found[u]) {
   1.285 +            found[u] = true;
   1.286 +            ++found_cnt;
   1.287 +            _policy[u] = e;
   1.288 +            queue.push_back(u);
   1.289 +          }
   1.290 +        }
   1.291        }
   1.292      }
   1.293  
   1.294 -    /// \brief Processes all rounds of computing required path data for
   1.295 -    /// the current component.
   1.296 -    void processRounds() {
   1.297 -      dmap[nodes[0]][0] = PathData(true, 0);
   1.298 -      process.clear();
   1.299 -      // Processing the first round
   1.300 -      for (OutEdgeIt e(graph, nodes[0]); e != INVALID; ++e) {
   1.301 -	Node v = graph.target(e);
   1.302 -	if (comp[v] != comp_cnt || v == nodes[0]) continue;
   1.303 -	dmap[v][1] = PathData(true, length[e], e);
   1.304 -	process.push_back(v);
   1.305 +    // Computes node distances in the policy graph and updates the
   1.306 +    // policy graph if the node distances can be improved.
   1.307 +    bool computeNodeDistances(int comp) {
   1.308 +      // Computing node distances using reverse BFS search
   1.309 +      double cycle_mean = (double)_cycle_length / _cycle_size;
   1.310 +      typename Graph::template NodeMap<int> found(_graph, false);
   1.311 +      std::deque<Node> queue;
   1.312 +      queue.push_back(_cycle_node);
   1.313 +      found[_cycle_node] = true;
   1.314 +      Node u, v;
   1.315 +      while (!queue.empty()) {
   1.316 +        v = queue.front(); queue.pop_front();
   1.317 +        for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
   1.318 +          u = _graph.source(e);
   1.319 +          if (_component[u] == comp && !found[u] && _policy[u] == e) {
   1.320 +            found[u] = true;
   1.321 +            _dist[u] = _dist[v] + _length[e] - cycle_mean;
   1.322 +            queue.push_back(u);
   1.323 +          }
   1.324 +        }
   1.325        }
   1.326 -      // Processing other rounds
   1.327 -      int n = nodes.size(), k;
   1.328 -      for (k = 2; k <= n && process.size() < n; ++k)
   1.329 -	processNextBuildRound(k);
   1.330 -      for ( ; k <= n; ++k)
   1.331 -	processNextFullRound(k);
   1.332 -    }
   1.333 -
   1.334 -    /// \brief Processes one round of computing required path data and
   1.335 -    /// rebuilds \ref process vector.
   1.336 -    void processNextBuildRound(int k) {
   1.337 -      NodeVector next;
   1.338 -      for (int i = 0; i < process.size(); ++i) {
   1.339 -	for (OutEdgeIt e(graph, process[i]); e != INVALID; ++e) {
   1.340 -	  Node v = graph.target(e);
   1.341 -	  if (comp[v] != comp_cnt) continue;
   1.342 -	  if (!dmap[v][k].found) {
   1.343 -	    next.push_back(v);
   1.344 -	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);
   1.345 -	  }
   1.346 -	  else if (dmap[process[i]][k-1].dist + length[e] < dmap[v][k].dist) {
   1.347 -	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);
   1.348 -	  }
   1.349 -	}
   1.350 +      // Improving node distances
   1.351 +      bool improved = false;
   1.352 +      for (int j = 0; j < _edges.size(); ++j) {
   1.353 +        Edge e = _edges[j];
   1.354 +        u = _graph.source(e); v = _graph.target(e);
   1.355 +        double delta = _dist[v] + _length[e] - cycle_mean;
   1.356 +        if (_tolerance.less(delta, _dist[u])) {
   1.357 +          improved = true;
   1.358 +          _dist[u] = delta;
   1.359 +          _policy[u] = e;
   1.360 +        }
   1.361        }
   1.362 -      process.swap(next);
   1.363 -    }
   1.364 -
   1.365 -    /// \brief Processes one round of computing required path data
   1.366 -    /// using \ref nodes vector instead of \ref process vector.
   1.367 -    void processNextFullRound(int k) {
   1.368 -      for (int i = 0; i < nodes.size(); ++i) {
   1.369 -	for (OutEdgeIt e(graph, nodes[i]); e != INVALID; ++e) {
   1.370 -	  Node v = graph.target(e);
   1.371 -	  if (comp[v] != comp_cnt) continue;
   1.372 -	  if ( !dmap[v][k].found ||
   1.373 -	       dmap[nodes[i]][k-1].dist + length[e] < dmap[v][k].dist ) {
   1.374 -	    dmap[v][k] = PathData(true, dmap[nodes[i]][k-1].dist + length[e], e);
   1.375 -	  }
   1.376 -	}
   1.377 -      }
   1.378 -    }
   1.379 -
   1.380 -    /// \brief Finds the minimum cycle mean value in the current
   1.381 -    /// component.
   1.382 -    bool findCurrentMin(Length &min_length, int &min_size, Node &min_node) {
   1.383 -      bool found_min = false;
   1.384 -      for (int i = 0; i < nodes.size(); ++i) {
   1.385 -	int n = nodes.size();
   1.386 -	if (!dmap[nodes[i]][n].found) continue;
   1.387 -	Length len;
   1.388 -	int size;
   1.389 -	bool found_one = false;
   1.390 -	for (int k = 0; k < n; ++k) {
   1.391 -	  if (!dmap[nodes[i]][k].found) continue;
   1.392 -	  Length _len = dmap[nodes[i]][n].dist - dmap[nodes[i]][k].dist;
   1.393 -	  int _size = n - k;
   1.394 -	  if (!found_one || len * _size < _len * size) {
   1.395 -	    found_one = true;
   1.396 -	    len = _len;
   1.397 -	    size = _size;
   1.398 -	  }
   1.399 -	}
   1.400 -	if ( found_one &&
   1.401 -	     (!found_min || len * min_size < min_length * size) ) {
   1.402 -	  found_min = true;
   1.403 -	  min_length = len;
   1.404 -	  min_size = size;
   1.405 -	  min_node = nodes[i];
   1.406 -	}
   1.407 -      }
   1.408 -      return found_min;
   1.409 +      return improved;
   1.410      }
   1.411  
   1.412    public:
   1.413 @@ -248,19 +266,18 @@
   1.414      ///
   1.415      /// Runs the algorithm.
   1.416      ///
   1.417 -    /// \return \c true if a cycle exists in the graph.
   1.418 +    /// \return Returns \c true if a directed cycle exists in the graph.
   1.419      ///
   1.420 -    /// \note Apart from the return value, <tt>m.run()</tt> is just a
   1.421 +    /// \note Apart from the return value, <tt>mmc.run()</tt> is just a
   1.422      /// shortcut of the following code.
   1.423      /// \code
   1.424 -    ///   m.init();
   1.425 -    ///   m.findMinMean();
   1.426 -    ///   m.findCycle();
   1.427 +    ///   mmc.init();
   1.428 +    ///   mmc.findMinMean();
   1.429 +    ///   mmc.findCycle();
   1.430      /// \endcode
   1.431      bool run() {
   1.432        init();
   1.433 -      findMinMean();
   1.434 -      return findCycle();
   1.435 +      return findMinMean() && findCycle();
   1.436      }
   1.437  
   1.438      /// \brief Initializes the internal data structures.
   1.439 @@ -269,75 +286,13 @@
   1.440      ///
   1.441      /// \sa reset()
   1.442      void init() {
   1.443 -      comp_num = stronglyConnectedComponents(graph, comp);
   1.444 -      if (!cycle_path) {
   1.445 -	local_path = true;
   1.446 -	cycle_path = new Path;
   1.447 +      _tolerance.epsilon(1e-8);
   1.448 +      if (!_cycle_path) {
   1.449 +        _local_path = true;
   1.450 +        _cycle_path = new Path;
   1.451        }
   1.452 -    }
   1.453 -
   1.454 -    /// \brief Finds the minimum cycle mean value in the graph.
   1.455 -    ///
   1.456 -    /// Computes all the required path data and finds the minimum cycle
   1.457 -    /// mean value in the graph.
   1.458 -    ///
   1.459 -    /// \return \c true if a cycle exists in the graph.
   1.460 -    ///
   1.461 -    /// \pre \ref init() must be called before using this function.
   1.462 -    bool findMinMean() {
   1.463 -      cycle_node = INVALID;
   1.464 -      for (comp_cnt = 0; comp_cnt < comp_num; ++comp_cnt) {
   1.465 -	initCurrent();
   1.466 -	processRounds();
   1.467 -
   1.468 -	Length min_length;
   1.469 -	int min_size;
   1.470 -	Node min_node;
   1.471 -	bool found_min = findCurrentMin(min_length, min_size, min_node);
   1.472 -
   1.473 -	if ( found_min && (cycle_node == INVALID ||
   1.474 -	     min_length * cycle_size < cycle_length * min_size) ) {
   1.475 -	  cycle_length = min_length;
   1.476 -	  cycle_size = min_size;
   1.477 -	  cycle_node = min_node;
   1.478 -	}
   1.479 -      }
   1.480 -      return (cycle_node != INVALID);
   1.481 -    }
   1.482 -
   1.483 -    /// \brief Finds a critical (minimum mean) cycle.
   1.484 -    ///
   1.485 -    /// Finds a critical (minimum mean) cycle using the path data
   1.486 -    /// stored in \ref dmap.
   1.487 -    ///
   1.488 -    /// \return \c true if a cycle exists in the graph.
   1.489 -    ///
   1.490 -    /// \pre \ref init() and \ref findMinMean() must be called before
   1.491 -    /// using this function.
   1.492 -    bool findCycle() {
   1.493 -      if (cycle_node == INVALID) return false;
   1.494 -      cycle_length = 0;
   1.495 -      cycle_size = 0;
   1.496 -      IntNodeMap reached(graph, -1);
   1.497 -      int r = reached[cycle_node] = dmap[cycle_node].size() - 1;
   1.498 -      Node u = graph.source(dmap[cycle_node][r].pred);
   1.499 -      while (reached[u] < 0) {
   1.500 -	reached[u] = --r;
   1.501 -	u = graph.source(dmap[u][r].pred);
   1.502 -      }
   1.503 -      r = reached[u];
   1.504 -      Edge e = dmap[u][r].pred;
   1.505 -      cycle_path->addFront(e);
   1.506 -      cycle_length = cycle_length + length[e];
   1.507 -      ++cycle_size;
   1.508 -      Node v;
   1.509 -      while ((v = graph.source(e)) != u) {
   1.510 -	e = dmap[v][--r].pred;
   1.511 -	cycle_path->addFront(e);
   1.512 -	cycle_length = cycle_length + length[e];
   1.513 -	++cycle_size;
   1.514 -      }
   1.515 -      return true;
   1.516 +      _cycle_found = false;
   1.517 +      _component_num = stronglyConnectedComponents(_graph, _component);
   1.518      }
   1.519  
   1.520      /// \brief Resets the internal data structures.
   1.521 @@ -348,33 +303,68 @@
   1.522      ///
   1.523      /// \sa init()
   1.524      void reset() {
   1.525 -      for (NodeIt u(graph); u != INVALID; ++u)
   1.526 -	dmap[u].clear();
   1.527 -      cycle_node = INVALID;
   1.528 -      if (cycle_path) cycle_path->clear();
   1.529 -      comp_num = stronglyConnectedComponents(graph, comp);
   1.530 +      if (_cycle_path) _cycle_path->clear();
   1.531 +      _cycle_found = false;
   1.532 +      _component_num = stronglyConnectedComponents(_graph, _component);
   1.533 +    }
   1.534 +
   1.535 +    /// \brief Finds the minimum cycle mean length in the graph.
   1.536 +    ///
   1.537 +    /// Computes all the required data and finds the minimum cycle mean
   1.538 +    /// length in the graph.
   1.539 +    ///
   1.540 +    /// \return Returns \c true if a directed cycle exists in the graph.
   1.541 +    ///
   1.542 +    /// \pre \ref init() must be called before using this function.
   1.543 +    bool findMinMean() {
   1.544 +      // Finding the minimum mean cycle in the components
   1.545 +      for (int comp = 0; comp < _component_num; ++comp) {
   1.546 +        if (!initCurrentComponent(comp)) continue;
   1.547 +        while (true) {
   1.548 +          if (!findPolicyCycles(comp)) return false;
   1.549 +          contractPolicyGraph(comp);
   1.550 +          if (!computeNodeDistances(comp)) return true;
   1.551 +        }
   1.552 +      }
   1.553 +    }
   1.554 +
   1.555 +    /// \brief Finds a critical (minimum mean) directed cycle.
   1.556 +    ///
   1.557 +    /// Finds a critical (minimum mean) directed cycle using the data
   1.558 +    /// computed in the \ref findMinMean() function.
   1.559 +    ///
   1.560 +    /// \return Returns \c true if a directed cycle exists in the graph.
   1.561 +    ///
   1.562 +    /// \pre \ref init() and \ref findMinMean() must be called before
   1.563 +    /// using this function.
   1.564 +    bool findCycle() {
   1.565 +      if (!_cycle_found) return false;
   1.566 +      _cycle_path->addBack(_policy[_cycle_node]);
   1.567 +      for ( Node v = _cycle_node;
   1.568 +            (v = _graph.target(_policy[v])) != _cycle_node; ) {
   1.569 +        _cycle_path->addBack(_policy[v]);
   1.570 +      }
   1.571 +      return true;
   1.572      }
   1.573  
   1.574      /// \brief Returns the total length of the found cycle.
   1.575      ///
   1.576      /// Returns the total length of the found cycle.
   1.577      ///
   1.578 -    /// \pre \ref run() or \ref findCycle() must be called before
   1.579 -    /// using this function. If only \ref findMinMean() is called,
   1.580 -    /// the return value is not valid.
   1.581 +    /// \pre \ref run() or \ref findMinMean() must be called before
   1.582 +    /// using this function.
   1.583      Length cycleLength() const {
   1.584 -      return cycle_length;
   1.585 +      return _cycle_length;
   1.586      }
   1.587  
   1.588 -    /// \brief Returns the number of edges in the found cycle.
   1.589 +    /// \brief Returns the number of edges on the found cycle.
   1.590      ///
   1.591 -    /// Returns the number of edges in the found cycle.
   1.592 +    /// Returns the number of edges on the found cycle.
   1.593      ///
   1.594 -    /// \pre \ref run() or \ref findCycle() must be called before
   1.595 -    /// using this function. If only \ref findMinMean() is called,
   1.596 -    /// the return value is not valid.
   1.597 +    /// \pre \ref run() or \ref findMinMean() must be called before
   1.598 +    /// using this function.
   1.599      int cycleEdgeNum() const {
   1.600 -      return cycle_size;
   1.601 +      return _cycle_size;
   1.602      }
   1.603  
   1.604      /// \brief Returns the mean length of the found cycle.
   1.605 @@ -384,52 +374,53 @@
   1.606      /// \pre \ref run() or \ref findMinMean() must be called before
   1.607      /// using this function.
   1.608      ///
   1.609 -    /// \warning LengthMap::Value must be convertible to double.
   1.610 -    ///
   1.611 -    /// \note <tt>m.minMean()</tt> is just a shortcut of the following
   1.612 -    /// code.
   1.613 +    /// \note <tt>mmc.cycleMean()</tt> is just a shortcut of the
   1.614 +    /// following code.
   1.615      /// \code
   1.616 -    ///   return m.cycleLength() / double(m.cycleEdgeNum());
   1.617 +    ///   return double(mmc.cycleLength()) / mmc.cycleEdgeNum();
   1.618      /// \endcode
   1.619 -    /// However if only \ref findMinMean() is called before using this
   1.620 -    /// function, <tt>m.cycleLength()</tt> and <tt>m.cycleEdgeNum()</tt>
   1.621 -    /// are not valid alone, only their ratio is valid.
   1.622 -    double minMean() const {
   1.623 -      return cycle_length / (double)cycle_size;
   1.624 +    double cycleMean() const {
   1.625 +      return (double)_cycle_length / _cycle_size;
   1.626      }
   1.627  
   1.628 -    /// \brief Returns a const reference to the \ref lemon::Path "Path"
   1.629 -    /// structure of the found cycle.
   1.630 +    /// \brief Returns a const reference to the \ref Path "path"
   1.631 +    /// structure storing the found cycle.
   1.632      ///
   1.633 -    /// Returns a const reference to the \ref lemon::Path "Path"
   1.634 -    /// structure of the found cycle.
   1.635 +    /// Returns a const reference to the \ref Path "path"
   1.636 +    /// structure storing the found cycle.
   1.637      ///
   1.638 -    /// \pre \ref run() must be called before using this function.
   1.639 +    /// \pre \ref run() or \ref findCycle() must be called before using
   1.640 +    /// this function.
   1.641      ///
   1.642 -    /// \sa \ref cyclePath()
   1.643 +    /// \sa cyclePath()
   1.644      const Path& cycle() const {
   1.645 -      return *cycle_path;
   1.646 +      return *_cycle_path;
   1.647      }
   1.648  
   1.649 -    /// \brief Sets the \ref lemon::Path "Path" structure storing the
   1.650 +    /// \brief Sets the \ref Path "path" structure for storing the found
   1.651 +    /// cycle.
   1.652 +    ///
   1.653 +    /// Sets an external \ref Path "path" structure for storing the
   1.654      /// found cycle.
   1.655      ///
   1.656 -    /// Sets the \ref lemon::Path "Path" structure storing the found
   1.657 -    /// cycle. If you don't use this function before calling
   1.658 -    /// \ref run(), it will allocate one. The destuctor deallocates
   1.659 -    /// this automatically allocated map, of course.
   1.660 +    /// If you don't call this function before calling \ref run() or
   1.661 +    /// \ref init(), it will allocate a local \ref Path "path"
   1.662 +    /// structure.
   1.663 +    /// The destuctor deallocates this automatically allocated map,
   1.664 +    /// of course.
   1.665      ///
   1.666 -    /// \note The algorithm calls only the \ref lemon::Path::addFront()
   1.667 -    /// "addFront()" method of the given \ref lemon::Path "Path"
   1.668 -    /// structure.
   1.669 +    /// \note The algorithm calls only the \ref lemon::Path::addBack()
   1.670 +    /// "addBack()" function of the given \ref Path "path" structure.
   1.671      ///
   1.672 -    /// \return \c (*this)
   1.673 +    /// \return <tt>(*this)</tt>
   1.674 +    ///
   1.675 +    /// \sa cycle()
   1.676      MinMeanCycle& cyclePath(Path &path) {
   1.677 -      if (local_path) {
   1.678 -	delete cycle_path;
   1.679 -	local_path = false;
   1.680 +      if (_local_path) {
   1.681 +        delete _cycle_path;
   1.682 +        _local_path = false;
   1.683        }
   1.684 -      cycle_path = &path;
   1.685 +      _cycle_path = &path;
   1.686        return *this;
   1.687      }
   1.688