# Changeset 2555:a84e52e99f57 in lemon-0.x

Ignore:
Timestamp:
01/13/08 11:26:55 (11 years ago)
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3436
Message:

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.

File:
1 edited

Unmodified
Removed
• ## lemon/min_mean_cycle.h

 r2553 /// /// \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 namespace lemon { /// @{ /// \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) {} }; private: typedef typename Graph::template NodeMap > PathDataNodeMap; typename Graph::template NodeMap _reached; typename Graph::template NodeMap _dist; typename Graph::template NodeMap _policy; /// The directed graph the algorithm runs on. const Graph &_graph; /// The length of the edges. const LengthMap &_length; /// 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; bool _local_path; bool _cycle_found; Node _cycle_node; typename Graph::template NodeMap _component; int _component_num; std::vector _nodes; std::vector _edges; Tolerance _tolerance; 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), _dist(graph), _reached(graph), _policy(graph), _component(graph), _cycle_length(0), _cycle_size(-1), _cycle_path(NULL), _local_path(false) { } /// The destructor of the class. ~MinMeanCycle() { if (_local_path) delete _cycle_path; } protected: /// \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; /// \brief The directed graph the algorithm runs on. const Graph &graph; /// \brief The length of the edges. const LengthMap &length; /// \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 : /// \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) { } /// \brief The destructor of the class. ~MinMeanCycle() { 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); } // Creating vectors for all nodes int n = nodes.size(); for (int i = 0; i < n; ++i) { dmap[nodes[i]].resize(n + 1); } } /// \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); } // 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); } } } 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; _nodes.clear(); for (NodeIt n(_graph); n != INVALID; ++n) { if (_component[n] == comp) _nodes.push_back(n); } 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); } } } } // 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); } } } // 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; } } return improved; } /// Runs the algorithm. /// /// \return \c true if a cycle exists in the graph. /// /// \note Apart from the return value, m.run() is just a /// \return Returns \c true if a directed cycle exists in the graph. /// /// \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(); } /// \sa reset() void init() { comp_num = stronglyConnectedComponents(graph, comp); 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. _tolerance.epsilon(1e-8); if (!_cycle_path) { _local_path = true; _cycle_path = new Path; } _cycle_found = false; _component_num = stronglyConnectedComponents(_graph, _component); } /// \brief Resets the internal data structures. /// /// Resets the internal data structures so that \ref findMinMean() /// and \ref findCycle() can be called again (e.g. when the /// underlaying graph has been modified). /// /// \sa init() void reset() { 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() { 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. // 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_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; 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 Resets the internal data structures. /// /// Resets the internal data structures so that \ref findMinMean() /// and \ref findCycle() can be called again (e.g. when the /// underlaying graph has been modified). /// /// \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); } /// \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. Length cycleLength() const { return cycle_length; } /// \brief Returns the number of edges in the found cycle. /// /// Returns the number of edges in 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. int cycleEdgeNum() const { return cycle_size; } /// \brief Returns the mean length of the found cycle. /// /// Returns the mean length of the found cycle. /// /// \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. Length cycleLength() const { return _cycle_length; } /// \brief Returns the number of edges on the found cycle. /// /// Returns the number of edges on the found cycle. /// /// \pre \ref run() or \ref findMinMean() must be called before /// using this function. int cycleEdgeNum() const { return _cycle_size; } /// \brief Returns the mean length of the found cycle. /// /// Returns the mean length of the found cycle. /// /// \pre \ref run() or \ref findMinMean() must be called before /// using this function. /// /// \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; } /// \brief Returns a const reference to the \ref lemon::Path "Path" /// structure of the found cycle. /// /// Returns a const reference to the \ref lemon::Path "Path" /// structure of the found cycle. /// /// \pre \ref run() must be called before using this function. /// /// \sa \ref cyclePath() double cycleMean() const { return (double)_cycle_length / _cycle_size; } /// \brief Returns a const reference to the \ref Path "path" /// structure storing the found cycle. /// /// Returns a const reference to the \ref Path "path" /// structure storing the found cycle. /// /// \pre \ref run() or \ref findCycle() must be called before using /// this function. /// /// \sa cyclePath() const Path& cycle() const { return *cycle_path; } /// \brief Sets the \ref lemon::Path "Path" structure storing the return *_cycle_path; } /// \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. /// /// \note The algorithm calls only the \ref lemon::Path::addFront() /// "addFront()" method of the given \ref lemon::Path "Path" /// If you don't call this function before calling \ref run() or /// \ref init(), it will allocate a local \ref Path "path" /// structure. /// /// \return \c (*this) /// The destuctor deallocates this automatically allocated map, /// of course. /// /// \note The algorithm calls only the \ref lemon::Path::addBack() /// "addBack()" function of the given \ref Path "path" structure. /// /// \return (*this) /// /// \sa cycle() MinMeanCycle& cyclePath(Path &path) { if (local_path) { delete cycle_path; local_path = false; } cycle_path = &path; if (_local_path) { delete _cycle_path; _local_path = false; } _cycle_path = &path; return *this; }
Note: See TracChangeset for help on using the changeset viewer.