# HG changeset patch # User kpeter # Date 1200220015 0 # Node ID a84e52e99f5726c81b9e86fc8a95375c13ed502b # Parent 1775aaa02ac4dd3b102006564e7a85bf3c88aa40 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. diff -r 1775aaa02ac4 -r a84e52e99f57 lemon/min_mean_cycle.h --- a/lemon/min_mean_cycle.h Mon Jan 07 17:07:40 2008 +0000 +++ b/lemon/min_mean_cycle.h Sun Jan 13 10:26:55 2008 +0000 @@ -22,224 +22,242 @@ /// \ingroup shortest_path /// /// \file -/// \brief Karp's algorithm for finding a minimum mean (directed) cycle. +/// \brief Howard's algorithm for finding a minimum mean directed cycle. #include #include #include #include +#include namespace lemon { /// \addtogroup shortest_path /// @{ - /// \brief Implementation of Karp's algorithm for finding a - /// minimum mean (directed) cycle. + /// \brief Implementation of Howard's algorithm for finding a minimum + /// mean directed cycle. /// - /// The \ref lemon::MinMeanCycle "MinMeanCycle" implements Karp's - /// algorithm for finding a minimum mean (directed) cycle. + /// \ref MinMeanCycle implements Howard's algorithm for finding a + /// minimum mean directed cycle. /// /// \param Graph The directed graph type the algorithm runs on. /// \param LengthMap The type of the length (cost) map. /// + /// \warning \c LengthMap::Value must be convertible to \c double. + /// /// \author Peter Kovacs - template > + template < typename Graph, + typename LengthMap = typename Graph::template EdgeMap > class MinMeanCycle { - typedef typename Graph::Node Node; - typedef typename Graph::NodeIt NodeIt; - typedef typename Graph::Edge Edge; - typedef typename Graph::EdgeIt EdgeIt; - typedef typename Graph::OutEdgeIt OutEdgeIt; + GRAPH_TYPEDEFS(typename Graph); typedef typename LengthMap::Value Length; - - typedef typename Graph::template NodeMap IntNodeMap; - typedef typename Graph::template NodeMap PredNodeMap; typedef Path Path; - typedef std::vector NodeVector; protected: - /// \brief Data structure for path data. - struct PathData - { - bool found; - Length dist; - Edge pred; - PathData(bool _found = false, Length _dist = 0) : - found(_found), dist(_dist), pred(INVALID) {} - PathData(bool _found, Length _dist, Edge _pred) : - found(_found), dist(_dist), pred(_pred) {} - }; + typename Graph::template NodeMap _reached; + typename Graph::template NodeMap _dist; + typename Graph::template NodeMap _policy; - private: + /// The directed graph the algorithm runs on. + const Graph &_graph; + /// The length of the edges. + const LengthMap &_length; - typedef typename Graph::template NodeMap > - PathDataNodeMap; + /// The total length of the found cycle. + Length _cycle_length; + /// The number of edges on the found cycle. + int _cycle_size; + /// The found cycle. + Path *_cycle_path; - protected: + bool _local_path; + bool _cycle_found; + Node _cycle_node; - /// \brief Node map for storing path data. - /// - /// Node map for storing path data of all nodes in the current - /// component. dmap[v][k] is the length of a shortest directed walk - /// to node v from the starting node containing exactly k edges. - PathDataNodeMap dmap; + typename Graph::template NodeMap _component; + int _component_num; - /// \brief The directed graph the algorithm runs on. - const Graph &graph; - /// \brief The length of the edges. - const LengthMap &length; + std::vector _nodes; + std::vector _edges; + Tolerance _tolerance; - /// \brief The total length of the found cycle. - Length cycle_length; - /// \brief The number of edges in the found cycle. - int cycle_size; - /// \brief A node for obtaining a minimum mean cycle. - Node cycle_node; - - /// \brief The found cycle. - Path *cycle_path; - /// \brief The algorithm uses local \ref lemon::Path "Path" - /// structure to store the found cycle. - bool local_path; - - /// \brief Node map for identifying strongly connected components. - IntNodeMap comp; - /// \brief The number of strongly connected components. - int comp_num; - /// \brief Counter for identifying the current component. - int comp_cnt; - /// \brief Nodes of the current component. - NodeVector nodes; - /// \brief The processed nodes in the last round. - NodeVector process; - - public : + public: /// \brief The constructor of the class. /// /// The constructor of the class. /// - /// \param _graph The directed graph the algorithm runs on. - /// \param _length The length (cost) of the edges. - MinMeanCycle( const Graph &_graph, - const LengthMap &_length ) : - graph(_graph), length(_length), dmap(_graph), comp(_graph), - cycle_length(0), cycle_size(-1), cycle_node(INVALID), - cycle_path(NULL), local_path(false) + /// \param graph The directed graph the algorithm runs on. + /// \param length The length (cost) of the edges. + MinMeanCycle( const Graph &graph, + const LengthMap &length ) : + _graph(graph), _length(length), _dist(graph), _reached(graph), + _policy(graph), _component(graph), _cycle_length(0), + _cycle_size(-1), _cycle_path(NULL), _local_path(false) { } - /// \brief The destructor of the class. + /// The destructor of the class. ~MinMeanCycle() { - if (local_path) delete cycle_path; + if (_local_path) delete _cycle_path; } protected: - /// \brief Initializes the internal data structures for the current - /// component. - void initCurrent() { - nodes.clear(); + // Initializes the internal data structures for the current strongly + // connected component and creating the policy graph. + // The policy graph can be represented by the _policy map because + // the out degree of every node is 1. + bool initCurrentComponent(int comp) { // Finding the nodes of the current component - for (NodeIt v(graph); v != INVALID; ++v) { - if (comp[v] == comp_cnt) nodes.push_back(v); + _nodes.clear(); + for (NodeIt n(_graph); n != INVALID; ++n) { + if (_component[n] == comp) _nodes.push_back(n); } - // Creating vectors for all nodes - int n = nodes.size(); - for (int i = 0; i < n; ++i) { - dmap[nodes[i]].resize(n + 1); + if (_nodes.size() <= 1) return false; + // Finding the edges of the current component + _edges.clear(); + for (EdgeIt e(_graph); e != INVALID; ++e) { + if ( _component[_graph.source(e)] == comp && + _component[_graph.target(e)] == comp ) + _edges.push_back(e); + } + // Initializing _reached, _dist, _policy maps + for (int i = 0; i < _nodes.size(); ++i) { + _reached[_nodes[i]] = false; + _policy[_nodes[i]] = INVALID; + } + Node u; Edge e; + for (int j = 0; j < _edges.size(); ++j) { + e = _edges[j]; + u = _graph.source(e); + if (!_reached[u] || _length[e] < _dist[u]) { + _dist[u] = _length[e]; + _policy[u] = e; + _reached[u] = true; + } + } + return true; + } + + // Finds all cycles in the policy graph. + // Sets _cycle_found to true if a cycle is found and sets + // _cycle_length, _cycle_size, _cycle_node to represent the minimum + // mean cycle in the policy graph. + bool findPolicyCycles(int comp) { + typename Graph::template NodeMap level(_graph, -1); + _cycle_found = false; + Length clength; + int csize; + int path_cnt = 0; + Node u, v; + // Searching for cycles + for (int i = 0; i < _nodes.size(); ++i) { + if (level[_nodes[i]] < 0) { + u = _nodes[i]; + level[u] = path_cnt; + while (level[u = _graph.target(_policy[u])] < 0) + level[u] = path_cnt; + if (level[u] == path_cnt) { + // A cycle is found + clength = _length[_policy[u]]; + csize = 1; + for (v = u; (v = _graph.target(_policy[v])) != u; ) { + clength += _length[_policy[v]]; + ++csize; + } + if ( !_cycle_found || + clength * _cycle_size < _cycle_length * csize ) { + _cycle_found = true; + _cycle_length = clength; + _cycle_size = csize; + _cycle_node = u; + } + } + ++path_cnt; + } + } + return _cycle_found; + } + + // Contracts the policy graph to be connected by cutting all cycles + // except for the main cycle (i.e. the minimum mean cycle). + void contractPolicyGraph(int comp) { + // Finding the component of the main cycle using + // reverse BFS search + typename Graph::template NodeMap found(_graph, false); + std::deque queue; + queue.push_back(_cycle_node); + found[_cycle_node] = true; + Node u, v; + while (!queue.empty()) { + v = queue.front(); queue.pop_front(); + for (InEdgeIt e(_graph, v); e != INVALID; ++e) { + u = _graph.source(e); + if (_component[u] == comp && !found[u] && _policy[u] == e) { + found[u] = true; + queue.push_back(u); + } + } + } + // Connecting all other nodes to this component using + // reverse BFS search + queue.clear(); + for (int i = 0; i < _nodes.size(); ++i) + if (found[_nodes[i]]) queue.push_back(_nodes[i]); + int found_cnt = queue.size(); + while (found_cnt < _nodes.size() && !queue.empty()) { + v = queue.front(); queue.pop_front(); + for (InEdgeIt e(_graph, v); e != INVALID; ++e) { + u = _graph.source(e); + if (_component[u] == comp && !found[u]) { + found[u] = true; + ++found_cnt; + _policy[u] = e; + queue.push_back(u); + } + } } } - /// \brief Processes all rounds of computing required path data for - /// the current component. - void processRounds() { - dmap[nodes[0]][0] = PathData(true, 0); - process.clear(); - // Processing the first round - for (OutEdgeIt e(graph, nodes[0]); e != INVALID; ++e) { - Node v = graph.target(e); - if (comp[v] != comp_cnt || v == nodes[0]) continue; - dmap[v][1] = PathData(true, length[e], e); - process.push_back(v); + // Computes node distances in the policy graph and updates the + // policy graph if the node distances can be improved. + bool computeNodeDistances(int comp) { + // Computing node distances using reverse BFS search + double cycle_mean = (double)_cycle_length / _cycle_size; + typename Graph::template NodeMap found(_graph, false); + std::deque queue; + queue.push_back(_cycle_node); + found[_cycle_node] = true; + Node u, v; + while (!queue.empty()) { + v = queue.front(); queue.pop_front(); + for (InEdgeIt e(_graph, v); e != INVALID; ++e) { + u = _graph.source(e); + if (_component[u] == comp && !found[u] && _policy[u] == e) { + found[u] = true; + _dist[u] = _dist[v] + _length[e] - cycle_mean; + queue.push_back(u); + } + } } - // Processing other rounds - int n = nodes.size(), k; - for (k = 2; k <= n && process.size() < n; ++k) - processNextBuildRound(k); - for ( ; k <= n; ++k) - processNextFullRound(k); - } - - /// \brief Processes one round of computing required path data and - /// rebuilds \ref process vector. - void processNextBuildRound(int k) { - NodeVector next; - for (int i = 0; i < process.size(); ++i) { - for (OutEdgeIt e(graph, process[i]); e != INVALID; ++e) { - Node v = graph.target(e); - if (comp[v] != comp_cnt) continue; - if (!dmap[v][k].found) { - next.push_back(v); - dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e); - } - else if (dmap[process[i]][k-1].dist + length[e] < dmap[v][k].dist) { - dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e); - } - } + // Improving node distances + bool improved = false; + for (int j = 0; j < _edges.size(); ++j) { + Edge e = _edges[j]; + u = _graph.source(e); v = _graph.target(e); + double delta = _dist[v] + _length[e] - cycle_mean; + if (_tolerance.less(delta, _dist[u])) { + improved = true; + _dist[u] = delta; + _policy[u] = e; + } } - process.swap(next); - } - - /// \brief Processes one round of computing required path data - /// using \ref nodes vector instead of \ref process vector. - void processNextFullRound(int k) { - for (int i = 0; i < nodes.size(); ++i) { - for (OutEdgeIt e(graph, nodes[i]); e != INVALID; ++e) { - Node v = graph.target(e); - if (comp[v] != comp_cnt) continue; - if ( !dmap[v][k].found || - dmap[nodes[i]][k-1].dist + length[e] < dmap[v][k].dist ) { - dmap[v][k] = PathData(true, dmap[nodes[i]][k-1].dist + length[e], e); - } - } - } - } - - /// \brief Finds the minimum cycle mean value in the current - /// component. - bool findCurrentMin(Length &min_length, int &min_size, Node &min_node) { - bool found_min = false; - for (int i = 0; i < nodes.size(); ++i) { - int n = nodes.size(); - if (!dmap[nodes[i]][n].found) continue; - Length len; - int size; - bool found_one = false; - for (int k = 0; k < n; ++k) { - if (!dmap[nodes[i]][k].found) continue; - Length _len = dmap[nodes[i]][n].dist - dmap[nodes[i]][k].dist; - int _size = n - k; - if (!found_one || len * _size < _len * size) { - found_one = true; - len = _len; - size = _size; - } - } - if ( found_one && - (!found_min || len * min_size < min_length * size) ) { - found_min = true; - min_length = len; - min_size = size; - min_node = nodes[i]; - } - } - return found_min; + return improved; } public: @@ -248,19 +266,18 @@ /// /// Runs the algorithm. /// - /// \return \c true if a cycle exists in the graph. + /// \return Returns \c true if a directed cycle exists in the graph. /// - /// \note Apart from the return value, m.run() is just a + /// \note Apart from the return value, mmc.run() is just a /// shortcut of the following code. /// \code - /// m.init(); - /// m.findMinMean(); - /// m.findCycle(); + /// mmc.init(); + /// mmc.findMinMean(); + /// mmc.findCycle(); /// \endcode bool run() { init(); - findMinMean(); - return findCycle(); + return findMinMean() && findCycle(); } /// \brief Initializes the internal data structures. @@ -269,75 +286,13 @@ /// /// \sa reset() void init() { - comp_num = stronglyConnectedComponents(graph, comp); - if (!cycle_path) { - local_path = true; - cycle_path = new Path; + _tolerance.epsilon(1e-8); + if (!_cycle_path) { + _local_path = true; + _cycle_path = new Path; } - } - - /// \brief Finds the minimum cycle mean value in the graph. - /// - /// Computes all the required path data and finds the minimum cycle - /// mean value in the graph. - /// - /// \return \c true if a cycle exists in the graph. - /// - /// \pre \ref init() must be called before using this function. - bool findMinMean() { - cycle_node = INVALID; - for (comp_cnt = 0; comp_cnt < comp_num; ++comp_cnt) { - initCurrent(); - processRounds(); - - Length min_length; - int min_size; - Node min_node; - bool found_min = findCurrentMin(min_length, min_size, min_node); - - if ( found_min && (cycle_node == INVALID || - min_length * cycle_size < cycle_length * min_size) ) { - cycle_length = min_length; - cycle_size = min_size; - cycle_node = min_node; - } - } - return (cycle_node != INVALID); - } - - /// \brief Finds a critical (minimum mean) cycle. - /// - /// Finds a critical (minimum mean) cycle using the path data - /// stored in \ref dmap. - /// - /// \return \c true if a cycle exists in the graph. - /// - /// \pre \ref init() and \ref findMinMean() must be called before - /// using this function. - bool findCycle() { - if (cycle_node == INVALID) return false; - cycle_length = 0; - cycle_size = 0; - IntNodeMap reached(graph, -1); - int r = reached[cycle_node] = dmap[cycle_node].size() - 1; - Node u = graph.source(dmap[cycle_node][r].pred); - while (reached[u] < 0) { - reached[u] = --r; - u = graph.source(dmap[u][r].pred); - } - r = reached[u]; - Edge e = dmap[u][r].pred; - cycle_path->addFront(e); - cycle_length = cycle_length + length[e]; - ++cycle_size; - Node v; - while ((v = graph.source(e)) != u) { - e = dmap[v][--r].pred; - cycle_path->addFront(e); - cycle_length = cycle_length + length[e]; - ++cycle_size; - } - return true; + _cycle_found = false; + _component_num = stronglyConnectedComponents(_graph, _component); } /// \brief Resets the internal data structures. @@ -348,33 +303,68 @@ /// /// \sa init() void reset() { - for (NodeIt u(graph); u != INVALID; ++u) - dmap[u].clear(); - cycle_node = INVALID; - if (cycle_path) cycle_path->clear(); - comp_num = stronglyConnectedComponents(graph, comp); + if (_cycle_path) _cycle_path->clear(); + _cycle_found = false; + _component_num = stronglyConnectedComponents(_graph, _component); + } + + /// \brief Finds the minimum cycle mean length in the graph. + /// + /// Computes all the required data and finds the minimum cycle mean + /// length in the graph. + /// + /// \return Returns \c true if a directed cycle exists in the graph. + /// + /// \pre \ref init() must be called before using this function. + bool findMinMean() { + // Finding the minimum mean cycle in the components + for (int comp = 0; comp < _component_num; ++comp) { + if (!initCurrentComponent(comp)) continue; + while (true) { + if (!findPolicyCycles(comp)) return false; + contractPolicyGraph(comp); + if (!computeNodeDistances(comp)) return true; + } + } + } + + /// \brief Finds a critical (minimum mean) directed cycle. + /// + /// Finds a critical (minimum mean) directed cycle using the data + /// computed in the \ref findMinMean() function. + /// + /// \return Returns \c true if a directed cycle exists in the graph. + /// + /// \pre \ref init() and \ref findMinMean() must be called before + /// using this function. + bool findCycle() { + if (!_cycle_found) return false; + _cycle_path->addBack(_policy[_cycle_node]); + for ( Node v = _cycle_node; + (v = _graph.target(_policy[v])) != _cycle_node; ) { + _cycle_path->addBack(_policy[v]); + } + return true; } /// \brief Returns the total length of the found cycle. /// /// Returns the total length of the found cycle. /// - /// \pre \ref run() or \ref findCycle() must be called before - /// using this function. If only \ref findMinMean() is called, - /// the return value is not valid. + /// \pre \ref run() or \ref findMinMean() must be called before + /// using this function. Length cycleLength() const { - return cycle_length; + return _cycle_length; } - /// \brief Returns the number of edges in the found cycle. + /// \brief Returns the number of edges on the found cycle. /// - /// Returns the number of edges in the found cycle. + /// Returns the number of edges on the found cycle. /// - /// \pre \ref run() or \ref findCycle() must be called before - /// using this function. If only \ref findMinMean() is called, - /// the return value is not valid. + /// \pre \ref run() or \ref findMinMean() must be called before + /// using this function. int cycleEdgeNum() const { - return cycle_size; + return _cycle_size; } /// \brief Returns the mean length of the found cycle. @@ -384,52 +374,53 @@ /// \pre \ref run() or \ref findMinMean() must be called before /// using this function. /// - /// \warning LengthMap::Value must be convertible to double. - /// - /// \note m.minMean() is just a shortcut of the following - /// code. + /// \note mmc.cycleMean() is just a shortcut of the + /// following code. /// \code - /// return m.cycleLength() / double(m.cycleEdgeNum()); + /// return double(mmc.cycleLength()) / mmc.cycleEdgeNum(); /// \endcode - /// However if only \ref findMinMean() is called before using this - /// function, m.cycleLength() and m.cycleEdgeNum() - /// are not valid alone, only their ratio is valid. - double minMean() const { - return cycle_length / (double)cycle_size; + double cycleMean() const { + return (double)_cycle_length / _cycle_size; } - /// \brief Returns a const reference to the \ref lemon::Path "Path" - /// structure of the found cycle. + /// \brief Returns a const reference to the \ref Path "path" + /// structure storing the found cycle. /// - /// Returns a const reference to the \ref lemon::Path "Path" - /// structure of the found cycle. + /// Returns a const reference to the \ref Path "path" + /// structure storing the found cycle. /// - /// \pre \ref run() must be called before using this function. + /// \pre \ref run() or \ref findCycle() must be called before using + /// this function. /// - /// \sa \ref cyclePath() + /// \sa cyclePath() const Path& cycle() const { - return *cycle_path; + return *_cycle_path; } - /// \brief Sets the \ref lemon::Path "Path" structure storing the + /// \brief Sets the \ref Path "path" structure for storing the found + /// cycle. + /// + /// Sets an external \ref Path "path" structure for storing the /// found cycle. /// - /// Sets the \ref lemon::Path "Path" structure storing the found - /// cycle. If you don't use this function before calling - /// \ref run(), it will allocate one. The destuctor deallocates - /// this automatically allocated map, of course. + /// If you don't call this function before calling \ref run() or + /// \ref init(), it will allocate a local \ref Path "path" + /// structure. + /// The destuctor deallocates this automatically allocated map, + /// of course. /// - /// \note The algorithm calls only the \ref lemon::Path::addFront() - /// "addFront()" method of the given \ref lemon::Path "Path" - /// structure. + /// \note The algorithm calls only the \ref lemon::Path::addBack() + /// "addBack()" function of the given \ref Path "path" structure. /// - /// \return \c (*this) + /// \return (*this) + /// + /// \sa cycle() MinMeanCycle& cyclePath(Path &path) { - if (local_path) { - delete cycle_path; - local_path = false; + if (_local_path) { + delete _cycle_path; + _local_path = false; } - cycle_path = &path; + _cycle_path = &path; return *this; }