1.1 --- a/lemon/min_mean_cycle.h Thu Feb 28 02:55:23 2008 +0000
1.2 +++ b/lemon/min_mean_cycle.h Thu Feb 28 02:57:36 2008 +0000
1.3 @@ -26,9 +26,9 @@
1.4
1.5 #include <vector>
1.6 #include <lemon/graph_utils.h>
1.7 -#include <lemon/topology.h>
1.8 #include <lemon/path.h>
1.9 #include <lemon/tolerance.h>
1.10 +#include <lemon/topology.h>
1.11
1.12 namespace lemon {
1.13
1.14 @@ -41,8 +41,8 @@
1.15 /// \ref MinMeanCycle implements Howard's algorithm for finding a
1.16 /// minimum mean directed cycle.
1.17 ///
1.18 - /// \param Graph The directed graph type the algorithm runs on.
1.19 - /// \param LengthMap The type of the length (cost) map.
1.20 + /// \tparam Graph The directed graph type the algorithm runs on.
1.21 + /// \tparam LengthMap The type of the length (cost) map.
1.22 ///
1.23 /// \warning \c LengthMap::Value must be convertible to \c double.
1.24 ///
1.25 @@ -57,28 +57,28 @@
1.26 typedef typename LengthMap::Value Length;
1.27 typedef Path<Graph> Path;
1.28
1.29 - protected:
1.30 + private:
1.31
1.32 - typename Graph::template NodeMap<bool> _reached;
1.33 - typename Graph::template NodeMap<double> _dist;
1.34 - typename Graph::template NodeMap<Edge> _policy;
1.35 -
1.36 - /// The directed graph the algorithm runs on.
1.37 + // The directed graph the algorithm runs on
1.38 const Graph &_graph;
1.39 - /// The length of the edges.
1.40 + // The length of the edges
1.41 const LengthMap &_length;
1.42
1.43 - /// The total length of the found cycle.
1.44 + // The total length of the found cycle
1.45 Length _cycle_length;
1.46 - /// The number of edges on the found cycle.
1.47 + // The number of edges on the found cycle
1.48 int _cycle_size;
1.49 - /// The found cycle.
1.50 + // The found cycle
1.51 Path *_cycle_path;
1.52
1.53 bool _local_path;
1.54 bool _cycle_found;
1.55 Node _cycle_node;
1.56
1.57 + typename Graph::template NodeMap<bool> _reached;
1.58 + typename Graph::template NodeMap<double> _dist;
1.59 + typename Graph::template NodeMap<Edge> _policy;
1.60 +
1.61 typename Graph::template NodeMap<int> _component;
1.62 int _component_num;
1.63
1.64 @@ -96,171 +96,55 @@
1.65 /// \param length The length (cost) of the edges.
1.66 MinMeanCycle( const Graph &graph,
1.67 const LengthMap &length ) :
1.68 - _graph(graph), _length(length), _dist(graph), _reached(graph),
1.69 - _policy(graph), _component(graph), _cycle_length(0),
1.70 - _cycle_size(-1), _cycle_path(NULL), _local_path(false)
1.71 - { }
1.72 + _graph(graph), _length(length), _cycle_length(0), _cycle_size(-1),
1.73 + _cycle_path(NULL), _local_path(false), _reached(graph),
1.74 + _dist(graph), _policy(graph), _component(graph)
1.75 + {}
1.76
1.77 /// The destructor of the class.
1.78 ~MinMeanCycle() {
1.79 if (_local_path) delete _cycle_path;
1.80 }
1.81
1.82 - protected:
1.83 -
1.84 - // Initializes the internal data structures for the current strongly
1.85 - // connected component and creating the policy graph.
1.86 - // The policy graph can be represented by the _policy map because
1.87 - // the out degree of every node is 1.
1.88 - bool initCurrentComponent(int comp) {
1.89 - // Finding the nodes of the current component
1.90 - _nodes.clear();
1.91 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.92 - if (_component[n] == comp) _nodes.push_back(n);
1.93 + /// \brief Sets the \ref Path "path" structure for storing the found
1.94 + /// cycle.
1.95 + ///
1.96 + /// Sets an external \ref Path "path" structure for storing the
1.97 + /// found cycle.
1.98 + ///
1.99 + /// If you don't call this function before calling \ref run() or
1.100 + /// \ref init(), it will allocate a local \ref Path "path"
1.101 + /// structure.
1.102 + /// The destuctor deallocates this automatically allocated map,
1.103 + /// of course.
1.104 + ///
1.105 + /// \note The algorithm calls only the \ref lemon::Path::addBack()
1.106 + /// "addBack()" function of the given \ref Path "path" structure.
1.107 + ///
1.108 + /// \return <tt>(*this)</tt>
1.109 + ///
1.110 + /// \sa cycle()
1.111 + MinMeanCycle& cyclePath(Path &path) {
1.112 + if (_local_path) {
1.113 + delete _cycle_path;
1.114 + _local_path = false;
1.115 }
1.116 - if (_nodes.size() <= 1) return false;
1.117 - // Finding the edges of the current component
1.118 - _edges.clear();
1.119 - for (EdgeIt e(_graph); e != INVALID; ++e) {
1.120 - if ( _component[_graph.source(e)] == comp &&
1.121 - _component[_graph.target(e)] == comp )
1.122 - _edges.push_back(e);
1.123 - }
1.124 - // Initializing _reached, _dist, _policy maps
1.125 - for (int i = 0; i < _nodes.size(); ++i) {
1.126 - _reached[_nodes[i]] = false;
1.127 - _policy[_nodes[i]] = INVALID;
1.128 - }
1.129 - Node u; Edge e;
1.130 - for (int j = 0; j < _edges.size(); ++j) {
1.131 - e = _edges[j];
1.132 - u = _graph.source(e);
1.133 - if (!_reached[u] || _length[e] < _dist[u]) {
1.134 - _dist[u] = _length[e];
1.135 - _policy[u] = e;
1.136 - _reached[u] = true;
1.137 - }
1.138 - }
1.139 - return true;
1.140 + _cycle_path = &path;
1.141 + return *this;
1.142 }
1.143
1.144 - // Finds all cycles in the policy graph.
1.145 - // Sets _cycle_found to true if a cycle is found and sets
1.146 - // _cycle_length, _cycle_size, _cycle_node to represent the minimum
1.147 - // mean cycle in the policy graph.
1.148 - bool findPolicyCycles(int comp) {
1.149 - typename Graph::template NodeMap<int> level(_graph, -1);
1.150 - _cycle_found = false;
1.151 - Length clength;
1.152 - int csize;
1.153 - int path_cnt = 0;
1.154 - Node u, v;
1.155 - // Searching for cycles
1.156 - for (int i = 0; i < _nodes.size(); ++i) {
1.157 - if (level[_nodes[i]] < 0) {
1.158 - u = _nodes[i];
1.159 - level[u] = path_cnt;
1.160 - while (level[u = _graph.target(_policy[u])] < 0)
1.161 - level[u] = path_cnt;
1.162 - if (level[u] == path_cnt) {
1.163 - // A cycle is found
1.164 - clength = _length[_policy[u]];
1.165 - csize = 1;
1.166 - for (v = u; (v = _graph.target(_policy[v])) != u; ) {
1.167 - clength += _length[_policy[v]];
1.168 - ++csize;
1.169 - }
1.170 - if ( !_cycle_found ||
1.171 - clength * _cycle_size < _cycle_length * csize ) {
1.172 - _cycle_found = true;
1.173 - _cycle_length = clength;
1.174 - _cycle_size = csize;
1.175 - _cycle_node = u;
1.176 - }
1.177 - }
1.178 - ++path_cnt;
1.179 - }
1.180 - }
1.181 - return _cycle_found;
1.182 - }
1.183 + /// \name Execution control
1.184 + /// The simplest way to execute the algorithm is to call run().
1.185 + /// \n
1.186 + /// If you only need the minimum mean value, you may call init()
1.187 + /// and findMinMean().
1.188 + /// \n
1.189 + /// If you would like to run the algorithm again (e.g. the
1.190 + /// underlaying graph and/or the edge costs were modified), you may
1.191 + /// not create a new instance of the class, rather call reset(),
1.192 + /// findMinMean(), and findCycle() instead.
1.193
1.194 - // Contracts the policy graph to be connected by cutting all cycles
1.195 - // except for the main cycle (i.e. the minimum mean cycle).
1.196 - void contractPolicyGraph(int comp) {
1.197 - // Finding the component of the main cycle using
1.198 - // reverse BFS search
1.199 - typename Graph::template NodeMap<int> found(_graph, false);
1.200 - std::deque<Node> queue;
1.201 - queue.push_back(_cycle_node);
1.202 - found[_cycle_node] = true;
1.203 - Node u, v;
1.204 - while (!queue.empty()) {
1.205 - v = queue.front(); queue.pop_front();
1.206 - for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.207 - u = _graph.source(e);
1.208 - if (_component[u] == comp && !found[u] && _policy[u] == e) {
1.209 - found[u] = true;
1.210 - queue.push_back(u);
1.211 - }
1.212 - }
1.213 - }
1.214 - // Connecting all other nodes to this component using
1.215 - // reverse BFS search
1.216 - queue.clear();
1.217 - for (int i = 0; i < _nodes.size(); ++i)
1.218 - if (found[_nodes[i]]) queue.push_back(_nodes[i]);
1.219 - int found_cnt = queue.size();
1.220 - while (found_cnt < _nodes.size() && !queue.empty()) {
1.221 - v = queue.front(); queue.pop_front();
1.222 - for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.223 - u = _graph.source(e);
1.224 - if (_component[u] == comp && !found[u]) {
1.225 - found[u] = true;
1.226 - ++found_cnt;
1.227 - _policy[u] = e;
1.228 - queue.push_back(u);
1.229 - }
1.230 - }
1.231 - }
1.232 - }
1.233 -
1.234 - // Computes node distances in the policy graph and updates the
1.235 - // policy graph if the node distances can be improved.
1.236 - bool computeNodeDistances(int comp) {
1.237 - // Computing node distances using reverse BFS search
1.238 - double cycle_mean = (double)_cycle_length / _cycle_size;
1.239 - typename Graph::template NodeMap<int> found(_graph, false);
1.240 - std::deque<Node> queue;
1.241 - queue.push_back(_cycle_node);
1.242 - found[_cycle_node] = true;
1.243 - Node u, v;
1.244 - while (!queue.empty()) {
1.245 - v = queue.front(); queue.pop_front();
1.246 - for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.247 - u = _graph.source(e);
1.248 - if (_component[u] == comp && !found[u] && _policy[u] == e) {
1.249 - found[u] = true;
1.250 - _dist[u] = _dist[v] + _length[e] - cycle_mean;
1.251 - queue.push_back(u);
1.252 - }
1.253 - }
1.254 - }
1.255 - // Improving node distances
1.256 - bool improved = false;
1.257 - for (int j = 0; j < _edges.size(); ++j) {
1.258 - Edge e = _edges[j];
1.259 - u = _graph.source(e); v = _graph.target(e);
1.260 - double delta = _dist[v] + _length[e] - cycle_mean;
1.261 - if (_tolerance.less(delta, _dist[u])) {
1.262 - improved = true;
1.263 - _dist[u] = delta;
1.264 - _policy[u] = e;
1.265 - }
1.266 - }
1.267 - return improved;
1.268 - }
1.269 -
1.270 - public:
1.271 + /// @{
1.272
1.273 /// \brief Runs the algorithm.
1.274 ///
1.275 @@ -286,7 +170,7 @@
1.276 ///
1.277 /// \sa reset()
1.278 void init() {
1.279 - _tolerance.epsilon(1e-8);
1.280 + _tolerance.epsilon(1e-6);
1.281 if (!_cycle_path) {
1.282 _local_path = true;
1.283 _cycle_path = new Path;
1.284 @@ -321,11 +205,12 @@
1.285 for (int comp = 0; comp < _component_num; ++comp) {
1.286 if (!initCurrentComponent(comp)) continue;
1.287 while (true) {
1.288 - if (!findPolicyCycles(comp)) return false;
1.289 + if (!findPolicyCycles()) break;
1.290 contractPolicyGraph(comp);
1.291 - if (!computeNodeDistances(comp)) return true;
1.292 + if (!computeNodeDistances(comp)) break;
1.293 }
1.294 }
1.295 + return _cycle_found;
1.296 }
1.297
1.298 /// \brief Finds a critical (minimum mean) directed cycle.
1.299 @@ -346,7 +231,16 @@
1.300 }
1.301 return true;
1.302 }
1.303 +
1.304 + /// @}
1.305
1.306 + /// \name Query Functions
1.307 + /// The result of the algorithm can be obtained using these
1.308 + /// functions.
1.309 + /// \n run() must be called before using them.
1.310 +
1.311 + /// @{
1.312 +
1.313 /// \brief Returns the total length of the found cycle.
1.314 ///
1.315 /// Returns the total length of the found cycle.
1.316 @@ -380,7 +274,7 @@
1.317 /// return double(mmc.cycleLength()) / mmc.cycleEdgeNum();
1.318 /// \endcode
1.319 double cycleMean() const {
1.320 - return (double)_cycle_length / _cycle_size;
1.321 + return double(_cycle_length) / _cycle_size;
1.322 }
1.323
1.324 /// \brief Returns a const reference to the \ref Path "path"
1.325 @@ -396,32 +290,163 @@
1.326 const Path& cycle() const {
1.327 return *_cycle_path;
1.328 }
1.329 +
1.330 + ///@}
1.331 +
1.332 + private:
1.333
1.334 - /// \brief Sets the \ref Path "path" structure for storing the found
1.335 - /// cycle.
1.336 - ///
1.337 - /// Sets an external \ref Path "path" structure for storing the
1.338 - /// found cycle.
1.339 - ///
1.340 - /// If you don't call this function before calling \ref run() or
1.341 - /// \ref init(), it will allocate a local \ref Path "path"
1.342 - /// structure.
1.343 - /// The destuctor deallocates this automatically allocated map,
1.344 - /// of course.
1.345 - ///
1.346 - /// \note The algorithm calls only the \ref lemon::Path::addBack()
1.347 - /// "addBack()" function of the given \ref Path "path" structure.
1.348 - ///
1.349 - /// \return <tt>(*this)</tt>
1.350 - ///
1.351 - /// \sa cycle()
1.352 - MinMeanCycle& cyclePath(Path &path) {
1.353 - if (_local_path) {
1.354 - delete _cycle_path;
1.355 - _local_path = false;
1.356 + // Initializes the internal data structures for the current strongly
1.357 + // connected component and creating the policy graph.
1.358 + // The policy graph can be represented by the _policy map because
1.359 + // the out degree of every node is 1.
1.360 + bool initCurrentComponent(int comp) {
1.361 + // Finding the nodes of the current component
1.362 + _nodes.clear();
1.363 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.364 + if (_component[n] == comp) _nodes.push_back(n);
1.365 }
1.366 - _cycle_path = &path;
1.367 - return *this;
1.368 + if (_nodes.size() <= 1) return false;
1.369 + // Finding the edges of the current component
1.370 + _edges.clear();
1.371 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.372 + if ( _component[_graph.source(e)] == comp &&
1.373 + _component[_graph.target(e)] == comp )
1.374 + _edges.push_back(e);
1.375 + }
1.376 + // Initializing _reached, _dist, _policy maps
1.377 + for (int i = 0; i < int(_nodes.size()); ++i) {
1.378 + _reached[_nodes[i]] = false;
1.379 + _policy[_nodes[i]] = INVALID;
1.380 + }
1.381 + Node u; Edge e;
1.382 + for (int j = 0; j < int(_edges.size()); ++j) {
1.383 + e = _edges[j];
1.384 + u = _graph.source(e);
1.385 + if (!_reached[u] || _length[e] < _dist[u]) {
1.386 + _dist[u] = _length[e];
1.387 + _policy[u] = e;
1.388 + _reached[u] = true;
1.389 + }
1.390 + }
1.391 + return true;
1.392 + }
1.393 +
1.394 + // Finds all cycles in the policy graph.
1.395 + // Sets _cycle_found to true if a cycle is found and sets
1.396 + // _cycle_length, _cycle_size, _cycle_node to represent the minimum
1.397 + // mean cycle in the policy graph.
1.398 + bool findPolicyCycles() {
1.399 + typename Graph::template NodeMap<int> level(_graph, -1);
1.400 + bool curr_cycle_found = false;
1.401 + Length clength;
1.402 + int csize;
1.403 + int path_cnt = 0;
1.404 + Node u, v;
1.405 + // Searching for cycles
1.406 + for (int i = 0; i < int(_nodes.size()); ++i) {
1.407 + if (level[_nodes[i]] < 0) {
1.408 + u = _nodes[i];
1.409 + level[u] = path_cnt;
1.410 + while (level[u = _graph.target(_policy[u])] < 0)
1.411 + level[u] = path_cnt;
1.412 + if (level[u] == path_cnt) {
1.413 + // A cycle is found
1.414 + curr_cycle_found = true;
1.415 + clength = _length[_policy[u]];
1.416 + csize = 1;
1.417 + for (v = u; (v = _graph.target(_policy[v])) != u; ) {
1.418 + clength += _length[_policy[v]];
1.419 + ++csize;
1.420 + }
1.421 + if ( !_cycle_found ||
1.422 + clength * _cycle_size < _cycle_length * csize ) {
1.423 + _cycle_found = true;
1.424 + _cycle_length = clength;
1.425 + _cycle_size = csize;
1.426 + _cycle_node = u;
1.427 + }
1.428 + }
1.429 + ++path_cnt;
1.430 + }
1.431 + }
1.432 + return curr_cycle_found;
1.433 + }
1.434 +
1.435 + // Contracts the policy graph to be connected by cutting all cycles
1.436 + // except for the main cycle (i.e. the minimum mean cycle).
1.437 + void contractPolicyGraph(int comp) {
1.438 + // Finding the component of the main cycle using
1.439 + // reverse BFS search
1.440 + typename Graph::template NodeMap<int> found(_graph, false);
1.441 + std::deque<Node> queue;
1.442 + queue.push_back(_cycle_node);
1.443 + found[_cycle_node] = true;
1.444 + Node u, v;
1.445 + while (!queue.empty()) {
1.446 + v = queue.front(); queue.pop_front();
1.447 + for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.448 + u = _graph.source(e);
1.449 + if (_component[u] == comp && !found[u] && _policy[u] == e) {
1.450 + found[u] = true;
1.451 + queue.push_back(u);
1.452 + }
1.453 + }
1.454 + }
1.455 + // Connecting all other nodes to this component using
1.456 + // reverse BFS search
1.457 + queue.clear();
1.458 + for (int i = 0; i < int(_nodes.size()); ++i)
1.459 + if (found[_nodes[i]]) queue.push_back(_nodes[i]);
1.460 + int found_cnt = queue.size();
1.461 + while (found_cnt < int(_nodes.size()) && !queue.empty()) {
1.462 + v = queue.front(); queue.pop_front();
1.463 + for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.464 + u = _graph.source(e);
1.465 + if (_component[u] == comp && !found[u]) {
1.466 + found[u] = true;
1.467 + ++found_cnt;
1.468 + _policy[u] = e;
1.469 + queue.push_back(u);
1.470 + }
1.471 + }
1.472 + }
1.473 + }
1.474 +
1.475 + // Computes node distances in the policy graph and updates the
1.476 + // policy graph if the node distances can be improved.
1.477 + bool computeNodeDistances(int comp) {
1.478 + // Computing node distances using reverse BFS search
1.479 + double cycle_mean = double(_cycle_length) / _cycle_size;
1.480 + typename Graph::template NodeMap<int> found(_graph, false);
1.481 + std::deque<Node> queue;
1.482 + queue.push_back(_cycle_node);
1.483 + found[_cycle_node] = true;
1.484 + _dist[_cycle_node] = 0;
1.485 + Node u, v;
1.486 + while (!queue.empty()) {
1.487 + v = queue.front(); queue.pop_front();
1.488 + for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.489 + u = _graph.source(e);
1.490 + if (_component[u] == comp && !found[u] && _policy[u] == e) {
1.491 + found[u] = true;
1.492 + _dist[u] = _dist[v] + _length[e] - cycle_mean;
1.493 + queue.push_back(u);
1.494 + }
1.495 + }
1.496 + }
1.497 + // Improving node distances
1.498 + bool improved = false;
1.499 + for (int j = 0; j < int(_edges.size()); ++j) {
1.500 + Edge e = _edges[j];
1.501 + u = _graph.source(e); v = _graph.target(e);
1.502 + double delta = _dist[v] + _length[e] - cycle_mean;
1.503 + if (_tolerance.less(delta, _dist[u])) {
1.504 + improved = true;
1.505 + _dist[u] = delta;
1.506 + _policy[u] = e;
1.507 + }
1.508 + }
1.509 + return improved;
1.510 }
1.511
1.512 }; //class MinMeanCycle