Changeset 807:83ce7ce39f21 in lemon
- Timestamp:
- 08/06/09 20:12:43 (15 years ago)
- Branch:
- default
- Phase:
- public
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/min_mean_cycle.h
r806 r807 75 75 const LengthMap &_length; 76 76 77 // The total length of the found cycle 78 Value _cycle_length; 79 // The number of arcs on the found cycle 80 int _cycle_size; 81 // The found cycle 77 // Data for the found cycles 78 bool _curr_found, _best_found; 79 Value _curr_length, _best_length; 80 int _curr_size, _best_size; 81 Node _curr_node, _best_node; 82 82 83 Path *_cycle_path; 83 84 84 bool _local_path; 85 bool _cycle_found; 86 Node _cycle_node;87 85 86 // Internal data used by the algorithm 87 typename Digraph::template NodeMap<Arc> _policy; 88 88 typename Digraph::template NodeMap<bool> _reached; 89 typename Digraph::template NodeMap<int> _level; 89 90 typename Digraph::template NodeMap<double> _dist; 90 typename Digraph::template NodeMap<Arc> _policy; 91 91 92 // Data for storing the strongly connected components 93 int _comp_num; 92 94 typename Digraph::template NodeMap<int> _comp; 93 int _comp_num; 94 95 std::vector<Node> _nodes; 96 std::vector<Arc> _arcs; 95 std::vector<std::vector<Node> > _comp_nodes; 96 std::vector<Node>* _nodes; 97 typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs; 98 99 // Queue used for BFS search 100 std::vector<Node> _queue; 101 int _qfront, _qback; 102 97 103 Tolerance<double> _tol; 98 104 … … 107 113 MinMeanCycle( const Digraph &digraph, 108 114 const LengthMap &length ) : 109 _gr(digraph), _length(length), _cycle_ length(0), _cycle_size(-1),110 _ cycle_path(NULL), _local_path(false), _reached(digraph),111 _ dist(digraph), _policy(digraph), _comp(digraph)115 _gr(digraph), _length(length), _cycle_path(NULL), _local_path(false), 116 _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph), 117 _comp(digraph), _in_arcs(digraph) 112 118 {} 113 119 … … 173 179 /// \return \c true if a directed cycle exists in the digraph. 174 180 bool findMinMean() { 175 // Initialize 181 // Initialize and find strongly connected components 182 init(); 183 findComponents(); 184 185 // Find the minimum cycle mean in the components 186 for (int comp = 0; comp < _comp_num; ++comp) { 187 // Find the minimum mean cycle in the current component 188 if (!buildPolicyGraph(comp)) continue; 189 while (true) { 190 findPolicyCycle(); 191 if (!computeNodeDistances()) break; 192 } 193 // Update the best cycle (global minimum mean cycle) 194 if ( !_best_found || (_curr_found && 195 _curr_length * _best_size < _best_length * _curr_size) ) { 196 _best_found = true; 197 _best_length = _curr_length; 198 _best_size = _curr_size; 199 _best_node = _curr_node; 200 } 201 } 202 return _best_found; 203 } 204 205 /// \brief Find a minimum mean directed cycle. 206 /// 207 /// This function finds a directed cycle of minimum mean length 208 /// in the digraph using the data computed by findMinMean(). 209 /// 210 /// \return \c true if a directed cycle exists in the digraph. 211 /// 212 /// \pre \ref findMinMean() must be called before using this function. 213 bool findCycle() { 214 if (!_best_found) return false; 215 _cycle_path->addBack(_policy[_best_node]); 216 for ( Node v = _best_node; 217 (v = _gr.target(_policy[v])) != _best_node; ) { 218 _cycle_path->addBack(_policy[v]); 219 } 220 return true; 221 } 222 223 /// @} 224 225 /// \name Query Functions 226 /// The results of the algorithm can be obtained using these 227 /// functions.\n 228 /// The algorithm should be executed before using them. 229 230 /// @{ 231 232 /// \brief Return the total length of the found cycle. 233 /// 234 /// This function returns the total length of the found cycle. 235 /// 236 /// \pre \ref run() or \ref findMinMean() must be called before 237 /// using this function. 238 Value cycleLength() const { 239 return _best_length; 240 } 241 242 /// \brief Return the number of arcs on the found cycle. 243 /// 244 /// This function returns the number of arcs on the found cycle. 245 /// 246 /// \pre \ref run() or \ref findMinMean() must be called before 247 /// using this function. 248 int cycleArcNum() const { 249 return _best_size; 250 } 251 252 /// \brief Return the mean length of the found cycle. 253 /// 254 /// This function returns the mean length of the found cycle. 255 /// 256 /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the 257 /// following code. 258 /// \code 259 /// return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum(); 260 /// \endcode 261 /// 262 /// \pre \ref run() or \ref findMinMean() must be called before 263 /// using this function. 264 double cycleMean() const { 265 return static_cast<double>(_best_length) / _best_size; 266 } 267 268 /// \brief Return the found cycle. 269 /// 270 /// This function returns a const reference to the path structure 271 /// storing the found cycle. 272 /// 273 /// \pre \ref run() or \ref findCycle() must be called before using 274 /// this function. 275 /// 276 /// \sa cyclePath() 277 const Path& cycle() const { 278 return *_cycle_path; 279 } 280 281 ///@} 282 283 private: 284 285 // Initialize 286 void init() { 176 287 _tol.epsilon(1e-6); 177 288 if (!_cycle_path) { … … 179 290 _cycle_path = new Path; 180 291 } 292 _queue.resize(countNodes(_gr)); 293 _best_found = false; 294 _best_length = 0; 295 _best_size = 1; 181 296 _cycle_path->clear(); 182 _cycle_found = false; 183 184 // Find the minimum cycle mean in the components 297 } 298 299 // Find strongly connected components and initialize _comp_nodes 300 // and _in_arcs 301 void findComponents() { 185 302 _comp_num = stronglyConnectedComponents(_gr, _comp); 186 for (int comp = 0; comp < _comp_num; ++comp) { 187 if (!initCurrentComponent(comp)) continue; 188 while (true) { 189 if (!findPolicyCycles()) break; 190 contractPolicyGraph(comp); 191 if (!computeNodeDistances()) break; 192 } 193 } 194 return _cycle_found; 195 } 196 197 /// \brief Find a minimum mean directed cycle. 198 /// 199 /// This function finds a directed cycle of minimum mean length 200 /// in the digraph using the data computed by findMinMean(). 201 /// 202 /// \return \c true if a directed cycle exists in the digraph. 203 /// 204 /// \pre \ref findMinMean() must be called before using this function. 205 bool findCycle() { 206 if (!_cycle_found) return false; 207 _cycle_path->addBack(_policy[_cycle_node]); 208 for ( Node v = _cycle_node; 209 (v = _gr.target(_policy[v])) != _cycle_node; ) { 210 _cycle_path->addBack(_policy[v]); 303 _comp_nodes.resize(_comp_num); 304 if (_comp_num == 1) { 305 _comp_nodes[0].clear(); 306 for (NodeIt n(_gr); n != INVALID; ++n) { 307 _comp_nodes[0].push_back(n); 308 _in_arcs[n].clear(); 309 for (InArcIt a(_gr, n); a != INVALID; ++a) { 310 _in_arcs[n].push_back(a); 311 } 312 } 313 } else { 314 for (int i = 0; i < _comp_num; ++i) 315 _comp_nodes[i].clear(); 316 for (NodeIt n(_gr); n != INVALID; ++n) { 317 int k = _comp[n]; 318 _comp_nodes[k].push_back(n); 319 _in_arcs[n].clear(); 320 for (InArcIt a(_gr, n); a != INVALID; ++a) { 321 if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a); 322 } 323 } 324 } 325 } 326 327 // Build the policy graph in the given strongly connected component 328 // (the out-degree of every node is 1) 329 bool buildPolicyGraph(int comp) { 330 _nodes = &(_comp_nodes[comp]); 331 if (_nodes->size() < 1 || 332 (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) { 333 return false; 334 } 335 for (int i = 0; i < int(_nodes->size()); ++i) { 336 _dist[(*_nodes)[i]] = std::numeric_limits<double>::max(); 337 } 338 Node u, v; 339 Arc e; 340 for (int i = 0; i < int(_nodes->size()); ++i) { 341 v = (*_nodes)[i]; 342 for (int j = 0; j < int(_in_arcs[v].size()); ++j) { 343 e = _in_arcs[v][j]; 344 u = _gr.source(e); 345 if (_length[e] < _dist[u]) { 346 _dist[u] = _length[e]; 347 _policy[u] = e; 348 } 349 } 211 350 } 212 351 return true; 213 352 } 214 353 215 /// @} 216 217 /// \name Query Functions 218 /// The results of the algorithm can be obtained using these 219 /// functions.\n 220 /// The algorithm should be executed before using them. 221 222 /// @{ 223 224 /// \brief Return the total length of the found cycle. 225 /// 226 /// This function returns the total length of the found cycle. 227 /// 228 /// \pre \ref run() or \ref findCycle() must be called before 229 /// using this function. 230 Value cycleLength() const { 231 return _cycle_length; 232 } 233 234 /// \brief Return the number of arcs on the found cycle. 235 /// 236 /// This function returns the number of arcs on the found cycle. 237 /// 238 /// \pre \ref run() or \ref findCycle() must be called before 239 /// using this function. 240 int cycleArcNum() const { 241 return _cycle_size; 242 } 243 244 /// \brief Return the mean length of the found cycle. 245 /// 246 /// This function returns the mean length of the found cycle. 247 /// 248 /// \note <tt>mmc.cycleMean()</tt> is just a shortcut of the 249 /// following code. 250 /// \code 251 /// return double(mmc.cycleLength()) / mmc.cycleArcNum(); 252 /// \endcode 253 /// 254 /// \pre \ref run() or \ref findMinMean() must be called before 255 /// using this function. 256 double cycleMean() const { 257 return double(_cycle_length) / _cycle_size; 258 } 259 260 /// \brief Return the found cycle. 261 /// 262 /// This function returns a const reference to the path structure 263 /// storing the found cycle. 264 /// 265 /// \pre \ref run() or \ref findCycle() must be called before using 266 /// this function. 267 /// 268 /// \sa cyclePath() 269 const Path& cycle() const { 270 return *_cycle_path; 271 } 272 273 ///@} 274 275 private: 276 277 // Initialize the internal data structures for the current strongly 278 // connected component and create the policy graph. 279 // The policy graph can be represented by the _policy map because 280 // the out-degree of every node is 1. 281 bool initCurrentComponent(int comp) { 282 // Find the nodes of the current component 283 _nodes.clear(); 284 for (NodeIt n(_gr); n != INVALID; ++n) { 285 if (_comp[n] == comp) _nodes.push_back(n); 286 } 287 if (_nodes.size() <= 1) return false; 288 // Find the arcs of the current component 289 _arcs.clear(); 290 for (ArcIt e(_gr); e != INVALID; ++e) { 291 if ( _comp[_gr.source(e)] == comp && 292 _comp[_gr.target(e)] == comp ) 293 _arcs.push_back(e); 294 } 295 // Initialize _reached, _dist, _policy maps 296 for (int i = 0; i < int(_nodes.size()); ++i) { 297 _reached[_nodes[i]] = false; 298 _policy[_nodes[i]] = INVALID; 299 } 300 Node u; Arc e; 301 for (int j = 0; j < int(_arcs.size()); ++j) { 302 e = _arcs[j]; 303 u = _gr.source(e); 304 if (!_reached[u] || _length[e] < _dist[u]) { 305 _dist[u] = _length[e]; 306 _policy[u] = e; 307 _reached[u] = true; 308 } 309 } 310 return true; 311 } 312 313 // Find all cycles in the policy graph. 314 // Set _cycle_found to true if a cycle is found and set 315 // _cycle_length, _cycle_size, _cycle_node to represent the minimum 316 // mean cycle in the policy graph. 317 bool findPolicyCycles() { 318 typename Digraph::template NodeMap<int> level(_gr, -1); 319 bool curr_cycle_found = false; 354 // Find the minimum mean cycle in the policy graph 355 void findPolicyCycle() { 356 for (int i = 0; i < int(_nodes->size()); ++i) { 357 _level[(*_nodes)[i]] = -1; 358 } 320 359 Value clength; 321 360 int csize; 322 int path_cnt = 0;323 361 Node u, v; 324 // Searching for cycles325 for (int i = 0; i < int(_nodes .size()); ++i) {326 if (level[_nodes[i]] < 0) {327 u = _nodes[i];328 level[u] = path_cnt;329 while (level[u = _gr.target(_policy[u])] < 0)330 level[u] = path_cnt;331 if (level[u] == path_cnt) {332 333 curr_cycle_found = true;334 clength = _length[_policy[u]];335 csize = 1;336 for (v = u; (v = _gr.target(_policy[v])) != u; ) {337 clength += _length[_policy[v]];338 ++csize;339 }340 if ( !_cycle_found ||341 clength * _cycle_size < _cycle_length * csize ) {342 _cycle_found = true;343 _cycle_length = clength;344 _cycle_size = csize;345 _cycle_node = u;346 347 348 ++path_cnt;349 } 350 }351 return curr_cycle_found;352 }353 354 // Contract the policy graph to be connected by cutting all cycles355 // except for the main cycle (i.e. the minimum mean cycle).356 void contractPolicyGraph(int comp) {357 // Find the component of the main cycle using reverse BFS search358 typename Digraph::template NodeMap<int> found(_gr, false);359 std::deque<Node> queue;360 queue.push_back(_cycle_node);361 found[_cycle_node] = true;362 _curr_found = false; 363 for (int i = 0; i < int(_nodes->size()); ++i) { 364 u = (*_nodes)[i]; 365 if (_level[u] >= 0) continue; 366 for (; _level[u] < 0; u = _gr.target(_policy[u])) { 367 _level[u] = i; 368 } 369 if (_level[u] == i) { 370 // A cycle is found 371 clength = _length[_policy[u]]; 372 csize = 1; 373 for (v = u; (v = _gr.target(_policy[v])) != u; ) { 374 clength += _length[_policy[v]]; 375 ++csize; 376 } 377 if ( !_curr_found || 378 (clength * _curr_size < _curr_length * csize) ) { 379 _curr_found = true; 380 _curr_length = clength; 381 _curr_size = csize; 382 _curr_node = u; 383 } 384 } 385 } 386 } 387 388 // Contract the policy graph and compute node distances 389 bool computeNodeDistances() { 390 // Find the component of the main cycle and compute node distances 391 // using reverse BFS 392 for (int i = 0; i < int(_nodes->size()); ++i) { 393 _reached[(*_nodes)[i]] = false; 394 } 395 double curr_mean = double(_curr_length) / _curr_size; 396 _qfront = _qback = 0; 397 _queue[0] = _curr_node; 398 _reached[_curr_node] = true; 399 _dist[_curr_node] = 0; 362 400 Node u, v; 363 while (!queue.empty()) { 364 v = queue.front(); queue.pop_front(); 365 for (InArcIt e(_gr, v); e != INVALID; ++e) { 401 Arc e; 402 while (_qfront <= _qback) { 403 v = _queue[_qfront++]; 404 for (int j = 0; j < int(_in_arcs[v].size()); ++j) { 405 e = _in_arcs[v][j]; 366 406 u = _gr.source(e); 367 if (_policy[u] == e && !found[u]) { 368 found[u] = true; 369 queue.push_back(u); 370 } 371 } 372 } 373 // Connect all other nodes to this component using reverse BFS search 374 queue.clear(); 375 for (int i = 0; i < int(_nodes.size()); ++i) 376 if (found[_nodes[i]]) queue.push_back(_nodes[i]); 377 int found_cnt = queue.size(); 378 while (found_cnt < int(_nodes.size())) { 379 v = queue.front(); queue.pop_front(); 380 for (InArcIt e(_gr, v); e != INVALID; ++e) { 407 if (_policy[u] == e && !_reached[u]) { 408 _reached[u] = true; 409 _dist[u] = _dist[v] + _length[e] - curr_mean; 410 _queue[++_qback] = u; 411 } 412 } 413 } 414 415 // Connect all other nodes to this component and compute node 416 // distances using reverse BFS 417 _qfront = 0; 418 while (_qback < int(_nodes->size())-1) { 419 v = _queue[_qfront++]; 420 for (int j = 0; j < int(_in_arcs[v].size()); ++j) { 421 e = _in_arcs[v][j]; 381 422 u = _gr.source(e); 382 if (_comp[u] == comp && !found[u]) { 383 found[u] = true; 384 ++found_cnt; 423 if (!_reached[u]) { 424 _reached[u] = true; 385 425 _policy[u] = e; 386 queue.push_back(u); 387 } 388 } 389 } 390 } 391 392 // Compute node distances in the policy graph and update the 393 // policy graph if the node distances can be improved. 394 bool computeNodeDistances() { 395 // Compute node distances using reverse BFS search 396 double cycle_mean = double(_cycle_length) / _cycle_size; 397 typename Digraph::template NodeMap<int> found(_gr, false); 398 std::deque<Node> queue; 399 queue.push_back(_cycle_node); 400 found[_cycle_node] = true; 401 _dist[_cycle_node] = 0; 402 Node u, v; 403 while (!queue.empty()) { 404 v = queue.front(); queue.pop_front(); 405 for (InArcIt e(_gr, v); e != INVALID; ++e) { 426 _dist[u] = _dist[v] + _length[e] - curr_mean; 427 _queue[++_qback] = u; 428 } 429 } 430 } 431 432 // Improve node distances 433 bool improved = false; 434 for (int i = 0; i < int(_nodes->size()); ++i) { 435 v = (*_nodes)[i]; 436 for (int j = 0; j < int(_in_arcs[v].size()); ++j) { 437 e = _in_arcs[v][j]; 406 438 u = _gr.source(e); 407 if (_policy[u] == e && !found[u]) { 408 found[u] = true; 409 _dist[u] = _dist[v] + _length[e] - cycle_mean; 410 queue.push_back(u); 411 } 412 } 413 } 414 // Improving node distances 415 bool improved = false; 416 for (int j = 0; j < int(_arcs.size()); ++j) { 417 Arc e = _arcs[j]; 418 u = _gr.source(e); v = _gr.target(e); 419 double delta = _dist[v] + _length[e] - cycle_mean; 420 if (_tol.less(delta, _dist[u])) { 421 improved = true; 422 _dist[u] = delta; 423 _policy[u] = e; 439 double delta = _dist[v] + _length[e] - curr_mean; 440 if (_tol.less(delta, _dist[u])) { 441 _dist[u] = delta; 442 _policy[u] = e; 443 improved = true; 444 } 424 445 } 425 446 }
Note: See TracChangeset
for help on using the changeset viewer.