alpar@2409: /* -*- C++ -*- alpar@2409: * alpar@2409: * This file is a part of LEMON, a generic C++ optimization library alpar@2409: * alpar@2553: * Copyright (C) 2003-2008 alpar@2409: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport alpar@2409: * (Egervary Research Group on Combinatorial Optimization, EGRES). alpar@2409: * alpar@2409: * Permission to use, modify and distribute this software is granted alpar@2409: * provided that this copyright notice appears in all copies. For alpar@2409: * precise terms see the accompanying LICENSE file. alpar@2409: * alpar@2409: * This software is provided "AS IS" with no warranty of any kind, alpar@2409: * express or implied, and with no claim as to its suitability for any alpar@2409: * purpose. alpar@2409: * alpar@2409: */ alpar@2409: alpar@2409: #ifndef LEMON_MIN_MEAN_CYCLE_H alpar@2409: #define LEMON_MIN_MEAN_CYCLE_H alpar@2409: deba@2529: /// \ingroup shortest_path alpar@2409: /// deba@2437: /// \file kpeter@2555: /// \brief Howard's algorithm for finding a minimum mean directed cycle. alpar@2409: kpeter@2509: #include alpar@2409: #include alpar@2409: #include alpar@2409: #include kpeter@2555: #include alpar@2409: alpar@2409: namespace lemon { alpar@2409: deba@2529: /// \addtogroup shortest_path alpar@2409: /// @{ alpar@2409: kpeter@2555: /// \brief Implementation of Howard's algorithm for finding a minimum kpeter@2555: /// mean directed cycle. alpar@2409: /// kpeter@2555: /// \ref MinMeanCycle implements Howard's algorithm for finding a kpeter@2555: /// minimum mean directed cycle. alpar@2409: /// alpar@2409: /// \param Graph The directed graph type the algorithm runs on. alpar@2409: /// \param LengthMap The type of the length (cost) map. alpar@2409: /// kpeter@2555: /// \warning \c LengthMap::Value must be convertible to \c double. kpeter@2555: /// alpar@2409: /// \author Peter Kovacs alpar@2409: kpeter@2555: template < typename Graph, kpeter@2555: typename LengthMap = typename Graph::template EdgeMap > alpar@2409: class MinMeanCycle alpar@2409: { kpeter@2555: GRAPH_TYPEDEFS(typename Graph); alpar@2409: alpar@2409: typedef typename LengthMap::Value Length; alpar@2409: typedef Path Path; alpar@2409: alpar@2409: protected: deba@2413: kpeter@2555: typename Graph::template NodeMap _reached; kpeter@2555: typename Graph::template NodeMap _dist; kpeter@2555: typename Graph::template NodeMap _policy; deba@2413: kpeter@2555: /// The directed graph the algorithm runs on. kpeter@2555: const Graph &_graph; kpeter@2555: /// The length of the edges. kpeter@2555: const LengthMap &_length; deba@2413: kpeter@2555: /// The total length of the found cycle. kpeter@2555: Length _cycle_length; kpeter@2555: /// The number of edges on the found cycle. kpeter@2555: int _cycle_size; kpeter@2555: /// The found cycle. kpeter@2555: Path *_cycle_path; deba@2413: kpeter@2555: bool _local_path; kpeter@2555: bool _cycle_found; kpeter@2555: Node _cycle_node; alpar@2409: kpeter@2555: typename Graph::template NodeMap _component; kpeter@2555: int _component_num; deba@2437: kpeter@2555: std::vector _nodes; kpeter@2555: std::vector _edges; kpeter@2555: Tolerance _tolerance; deba@2437: kpeter@2555: public: alpar@2409: alpar@2409: /// \brief The constructor of the class. alpar@2409: /// alpar@2409: /// The constructor of the class. alpar@2409: /// kpeter@2555: /// \param graph The directed graph the algorithm runs on. kpeter@2555: /// \param length The length (cost) of the edges. kpeter@2555: MinMeanCycle( const Graph &graph, kpeter@2555: const LengthMap &length ) : kpeter@2555: _graph(graph), _length(length), _dist(graph), _reached(graph), kpeter@2555: _policy(graph), _component(graph), _cycle_length(0), kpeter@2555: _cycle_size(-1), _cycle_path(NULL), _local_path(false) alpar@2409: { } alpar@2409: kpeter@2555: /// The destructor of the class. deba@2437: ~MinMeanCycle() { kpeter@2555: if (_local_path) delete _cycle_path; alpar@2409: } alpar@2409: alpar@2409: protected: alpar@2409: kpeter@2555: // Initializes the internal data structures for the current strongly kpeter@2555: // connected component and creating the policy graph. kpeter@2555: // The policy graph can be represented by the _policy map because kpeter@2555: // the out degree of every node is 1. kpeter@2555: bool initCurrentComponent(int comp) { alpar@2409: // Finding the nodes of the current component kpeter@2555: _nodes.clear(); kpeter@2555: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@2555: if (_component[n] == comp) _nodes.push_back(n); alpar@2409: } kpeter@2555: if (_nodes.size() <= 1) return false; kpeter@2555: // Finding the edges of the current component kpeter@2555: _edges.clear(); kpeter@2555: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2555: if ( _component[_graph.source(e)] == comp && kpeter@2555: _component[_graph.target(e)] == comp ) kpeter@2555: _edges.push_back(e); kpeter@2555: } kpeter@2555: // Initializing _reached, _dist, _policy maps kpeter@2555: for (int i = 0; i < _nodes.size(); ++i) { kpeter@2555: _reached[_nodes[i]] = false; kpeter@2555: _policy[_nodes[i]] = INVALID; kpeter@2555: } kpeter@2555: Node u; Edge e; kpeter@2555: for (int j = 0; j < _edges.size(); ++j) { kpeter@2555: e = _edges[j]; kpeter@2555: u = _graph.source(e); kpeter@2555: if (!_reached[u] || _length[e] < _dist[u]) { kpeter@2555: _dist[u] = _length[e]; kpeter@2555: _policy[u] = e; kpeter@2555: _reached[u] = true; kpeter@2555: } kpeter@2555: } kpeter@2555: return true; kpeter@2555: } kpeter@2555: kpeter@2555: // Finds all cycles in the policy graph. kpeter@2555: // Sets _cycle_found to true if a cycle is found and sets kpeter@2555: // _cycle_length, _cycle_size, _cycle_node to represent the minimum kpeter@2555: // mean cycle in the policy graph. kpeter@2555: bool findPolicyCycles(int comp) { kpeter@2555: typename Graph::template NodeMap level(_graph, -1); kpeter@2555: _cycle_found = false; kpeter@2555: Length clength; kpeter@2555: int csize; kpeter@2555: int path_cnt = 0; kpeter@2555: Node u, v; kpeter@2555: // Searching for cycles kpeter@2555: for (int i = 0; i < _nodes.size(); ++i) { kpeter@2555: if (level[_nodes[i]] < 0) { kpeter@2555: u = _nodes[i]; kpeter@2555: level[u] = path_cnt; kpeter@2555: while (level[u = _graph.target(_policy[u])] < 0) kpeter@2555: level[u] = path_cnt; kpeter@2555: if (level[u] == path_cnt) { kpeter@2555: // A cycle is found kpeter@2555: clength = _length[_policy[u]]; kpeter@2555: csize = 1; kpeter@2555: for (v = u; (v = _graph.target(_policy[v])) != u; ) { kpeter@2555: clength += _length[_policy[v]]; kpeter@2555: ++csize; kpeter@2555: } kpeter@2555: if ( !_cycle_found || kpeter@2555: clength * _cycle_size < _cycle_length * csize ) { kpeter@2555: _cycle_found = true; kpeter@2555: _cycle_length = clength; kpeter@2555: _cycle_size = csize; kpeter@2555: _cycle_node = u; kpeter@2555: } kpeter@2555: } kpeter@2555: ++path_cnt; kpeter@2555: } kpeter@2555: } kpeter@2555: return _cycle_found; kpeter@2555: } kpeter@2555: kpeter@2555: // Contracts the policy graph to be connected by cutting all cycles kpeter@2555: // except for the main cycle (i.e. the minimum mean cycle). kpeter@2555: void contractPolicyGraph(int comp) { kpeter@2555: // Finding the component of the main cycle using kpeter@2555: // reverse BFS search kpeter@2555: typename Graph::template NodeMap found(_graph, false); kpeter@2555: std::deque queue; kpeter@2555: queue.push_back(_cycle_node); kpeter@2555: found[_cycle_node] = true; kpeter@2555: Node u, v; kpeter@2555: while (!queue.empty()) { kpeter@2555: v = queue.front(); queue.pop_front(); kpeter@2555: for (InEdgeIt e(_graph, v); e != INVALID; ++e) { kpeter@2555: u = _graph.source(e); kpeter@2555: if (_component[u] == comp && !found[u] && _policy[u] == e) { kpeter@2555: found[u] = true; kpeter@2555: queue.push_back(u); kpeter@2555: } kpeter@2555: } kpeter@2555: } kpeter@2555: // Connecting all other nodes to this component using kpeter@2555: // reverse BFS search kpeter@2555: queue.clear(); kpeter@2555: for (int i = 0; i < _nodes.size(); ++i) kpeter@2555: if (found[_nodes[i]]) queue.push_back(_nodes[i]); kpeter@2555: int found_cnt = queue.size(); kpeter@2555: while (found_cnt < _nodes.size() && !queue.empty()) { kpeter@2555: v = queue.front(); queue.pop_front(); kpeter@2555: for (InEdgeIt e(_graph, v); e != INVALID; ++e) { kpeter@2555: u = _graph.source(e); kpeter@2555: if (_component[u] == comp && !found[u]) { kpeter@2555: found[u] = true; kpeter@2555: ++found_cnt; kpeter@2555: _policy[u] = e; kpeter@2555: queue.push_back(u); kpeter@2555: } kpeter@2555: } alpar@2409: } alpar@2409: } deba@2437: kpeter@2555: // Computes node distances in the policy graph and updates the kpeter@2555: // policy graph if the node distances can be improved. kpeter@2555: bool computeNodeDistances(int comp) { kpeter@2555: // Computing node distances using reverse BFS search kpeter@2555: double cycle_mean = (double)_cycle_length / _cycle_size; kpeter@2555: typename Graph::template NodeMap found(_graph, false); kpeter@2555: std::deque queue; kpeter@2555: queue.push_back(_cycle_node); kpeter@2555: found[_cycle_node] = true; kpeter@2555: Node u, v; kpeter@2555: while (!queue.empty()) { kpeter@2555: v = queue.front(); queue.pop_front(); kpeter@2555: for (InEdgeIt e(_graph, v); e != INVALID; ++e) { kpeter@2555: u = _graph.source(e); kpeter@2555: if (_component[u] == comp && !found[u] && _policy[u] == e) { kpeter@2555: found[u] = true; kpeter@2555: _dist[u] = _dist[v] + _length[e] - cycle_mean; kpeter@2555: queue.push_back(u); kpeter@2555: } kpeter@2555: } alpar@2409: } kpeter@2555: // Improving node distances kpeter@2555: bool improved = false; kpeter@2555: for (int j = 0; j < _edges.size(); ++j) { kpeter@2555: Edge e = _edges[j]; kpeter@2555: u = _graph.source(e); v = _graph.target(e); kpeter@2555: double delta = _dist[v] + _length[e] - cycle_mean; kpeter@2555: if (_tolerance.less(delta, _dist[u])) { kpeter@2555: improved = true; kpeter@2555: _dist[u] = delta; kpeter@2555: _policy[u] = e; kpeter@2555: } alpar@2409: } kpeter@2555: return improved; alpar@2409: } deba@2437: deba@2437: public: deba@2437: alpar@2409: /// \brief Runs the algorithm. alpar@2409: /// alpar@2409: /// Runs the algorithm. alpar@2409: /// kpeter@2555: /// \return Returns \c true if a directed cycle exists in the graph. deba@2413: /// kpeter@2555: /// \note Apart from the return value, mmc.run() is just a kpeter@2517: /// shortcut of the following code. deba@2413: /// \code kpeter@2555: /// mmc.init(); kpeter@2555: /// mmc.findMinMean(); kpeter@2555: /// mmc.findCycle(); deba@2413: /// \endcode alpar@2409: bool run() { alpar@2409: init(); kpeter@2555: return findMinMean() && findCycle(); deba@2413: } deba@2437: deba@2413: /// \brief Initializes the internal data structures. kpeter@2517: /// kpeter@2517: /// Initializes the internal data structures. kpeter@2517: /// kpeter@2517: /// \sa reset() deba@2413: void init() { kpeter@2555: _tolerance.epsilon(1e-8); kpeter@2555: if (!_cycle_path) { kpeter@2555: _local_path = true; kpeter@2555: _cycle_path = new Path; deba@2413: } kpeter@2555: _cycle_found = false; kpeter@2555: _component_num = stronglyConnectedComponents(_graph, _component); deba@2413: } deba@2413: deba@2413: /// \brief Resets the internal data structures. deba@2413: /// deba@2437: /// Resets the internal data structures so that \ref findMinMean() deba@2437: /// and \ref findCycle() can be called again (e.g. when the deba@2413: /// underlaying graph has been modified). kpeter@2517: /// kpeter@2517: /// \sa init() deba@2413: void reset() { kpeter@2555: if (_cycle_path) _cycle_path->clear(); kpeter@2555: _cycle_found = false; kpeter@2555: _component_num = stronglyConnectedComponents(_graph, _component); kpeter@2555: } kpeter@2555: kpeter@2555: /// \brief Finds the minimum cycle mean length in the graph. kpeter@2555: /// kpeter@2555: /// Computes all the required data and finds the minimum cycle mean kpeter@2555: /// length in the graph. kpeter@2555: /// kpeter@2555: /// \return Returns \c true if a directed cycle exists in the graph. kpeter@2555: /// kpeter@2555: /// \pre \ref init() must be called before using this function. kpeter@2555: bool findMinMean() { kpeter@2555: // Finding the minimum mean cycle in the components kpeter@2555: for (int comp = 0; comp < _component_num; ++comp) { kpeter@2555: if (!initCurrentComponent(comp)) continue; kpeter@2555: while (true) { kpeter@2555: if (!findPolicyCycles(comp)) return false; kpeter@2555: contractPolicyGraph(comp); kpeter@2555: if (!computeNodeDistances(comp)) return true; kpeter@2555: } kpeter@2555: } kpeter@2555: } kpeter@2555: kpeter@2555: /// \brief Finds a critical (minimum mean) directed cycle. kpeter@2555: /// kpeter@2555: /// Finds a critical (minimum mean) directed cycle using the data kpeter@2555: /// computed in the \ref findMinMean() function. kpeter@2555: /// kpeter@2555: /// \return Returns \c true if a directed cycle exists in the graph. kpeter@2555: /// kpeter@2555: /// \pre \ref init() and \ref findMinMean() must be called before kpeter@2555: /// using this function. kpeter@2555: bool findCycle() { kpeter@2555: if (!_cycle_found) return false; kpeter@2555: _cycle_path->addBack(_policy[_cycle_node]); kpeter@2555: for ( Node v = _cycle_node; kpeter@2555: (v = _graph.target(_policy[v])) != _cycle_node; ) { kpeter@2555: _cycle_path->addBack(_policy[v]); kpeter@2555: } kpeter@2555: return true; deba@2413: } deba@2437: deba@2413: /// \brief Returns the total length of the found cycle. deba@2413: /// deba@2413: /// Returns the total length of the found cycle. alpar@2409: /// kpeter@2555: /// \pre \ref run() or \ref findMinMean() must be called before kpeter@2555: /// using this function. deba@2437: Length cycleLength() const { kpeter@2555: return _cycle_length; alpar@2409: } deba@2437: kpeter@2555: /// \brief Returns the number of edges on the found cycle. alpar@2409: /// kpeter@2555: /// Returns the number of edges on the found cycle. alpar@2409: /// kpeter@2555: /// \pre \ref run() or \ref findMinMean() must be called before kpeter@2555: /// using this function. deba@2437: int cycleEdgeNum() const { kpeter@2555: return _cycle_size; alpar@2409: } deba@2437: deba@2413: /// \brief Returns the mean length of the found cycle. alpar@2409: /// deba@2413: /// Returns the mean length of the found cycle. alpar@2409: /// kpeter@2517: /// \pre \ref run() or \ref findMinMean() must be called before kpeter@2517: /// using this function. alpar@2409: /// kpeter@2555: /// \note mmc.cycleMean() is just a shortcut of the kpeter@2555: /// following code. deba@2413: /// \code kpeter@2555: /// return double(mmc.cycleLength()) / mmc.cycleEdgeNum(); deba@2413: /// \endcode kpeter@2555: double cycleMean() const { kpeter@2555: return (double)_cycle_length / _cycle_size; alpar@2409: } alpar@2409: kpeter@2555: /// \brief Returns a const reference to the \ref Path "path" kpeter@2555: /// structure storing the found cycle. alpar@2409: /// kpeter@2555: /// Returns a const reference to the \ref Path "path" kpeter@2555: /// structure storing the found cycle. alpar@2409: /// kpeter@2555: /// \pre \ref run() or \ref findCycle() must be called before using kpeter@2555: /// this function. alpar@2409: /// kpeter@2555: /// \sa cyclePath() deba@2437: const Path& cycle() const { kpeter@2555: return *_cycle_path; alpar@2409: } deba@2437: kpeter@2555: /// \brief Sets the \ref Path "path" structure for storing the found kpeter@2555: /// cycle. kpeter@2555: /// kpeter@2555: /// Sets an external \ref Path "path" structure for storing the alpar@2409: /// found cycle. deba@2437: /// kpeter@2555: /// If you don't call this function before calling \ref run() or kpeter@2555: /// \ref init(), it will allocate a local \ref Path "path" kpeter@2555: /// structure. kpeter@2555: /// The destuctor deallocates this automatically allocated map, kpeter@2555: /// of course. alpar@2409: /// kpeter@2555: /// \note The algorithm calls only the \ref lemon::Path::addBack() kpeter@2555: /// "addBack()" function of the given \ref Path "path" structure. deba@2437: /// kpeter@2555: /// \return (*this) kpeter@2555: /// kpeter@2555: /// \sa cycle() deba@2437: MinMeanCycle& cyclePath(Path &path) { kpeter@2555: if (_local_path) { kpeter@2555: delete _cycle_path; kpeter@2555: _local_path = false; alpar@2409: } kpeter@2555: _cycle_path = &path; alpar@2409: return *this; alpar@2409: } alpar@2409: alpar@2409: }; //class MinMeanCycle alpar@2409: alpar@2409: ///@} alpar@2409: alpar@2409: } //namespace lemon alpar@2409: alpar@2409: #endif //LEMON_MIN_MEAN_CYCLE_H