... | ... |
@@ -74,26 +74,32 @@ |
74 | 74 |
// The length of the arcs |
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 |
// |
|
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; |
84 |
bool _local_path; |
|
83 | 85 |
|
84 |
bool _local_path; |
|
85 |
bool _cycle_found; |
|
86 |
|
|
86 |
// Internal data used by the algorithm |
|
87 |
typename Digraph::template NodeMap<Arc> _policy; |
|
88 |
typename Digraph::template NodeMap<bool> _reached; |
|
89 |
typename Digraph::template NodeMap<int> _level; |
|
90 |
typename Digraph::template NodeMap<double> _dist; |
|
87 | 91 |
|
88 |
typename Digraph::template NodeMap<bool> _reached; |
|
89 |
typename Digraph::template NodeMap<double> _dist; |
|
90 |
typename Digraph::template NodeMap<Arc> _policy; |
|
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 |
|
99 | 105 |
public: |
... | ... |
@@ -106,9 +112,9 @@ |
106 | 112 |
/// \param length The lengths (costs) of the arcs. |
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 |
|
|
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 |
|
114 | 120 |
/// Destructor. |
... | ... |
@@ -172,26 +178,28 @@ |
172 | 178 |
/// |
173 | 179 |
/// \return \c true if a directed cycle exists in the digraph. |
174 | 180 |
bool findMinMean() { |
175 |
// Initialize |
|
176 |
_tol.epsilon(1e-6); |
|
177 |
if (!_cycle_path) { |
|
178 |
_local_path = true; |
|
179 |
_cycle_path = new Path; |
|
180 |
} |
|
181 |
_cycle_path->clear(); |
|
182 |
_cycle_found = false; |
|
183 |
|
|
181 |
// Initialize and find strongly connected components |
|
182 |
init(); |
|
183 |
findComponents(); |
|
184 |
|
|
184 | 185 |
// Find the minimum cycle mean in the components |
185 |
_comp_num = stronglyConnectedComponents(_gr, _comp); |
|
186 | 186 |
for (int comp = 0; comp < _comp_num; ++comp) { |
187 |
|
|
187 |
// Find the minimum mean cycle in the current component |
|
188 |
if (!buildPolicyGraph(comp)) continue; |
|
188 | 189 |
while (true) { |
189 |
if (!findPolicyCycles()) break; |
|
190 |
contractPolicyGraph(comp); |
|
190 |
findPolicyCycle(); |
|
191 | 191 |
if (!computeNodeDistances()) break; |
192 | 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 |
} |
|
193 | 201 |
} |
194 |
return |
|
202 |
return _best_found; |
|
195 | 203 |
} |
196 | 204 |
|
197 | 205 |
/// \brief Find a minimum mean directed cycle. |
... | ... |
@@ -203,10 +211,10 @@ |
203 | 211 |
/// |
204 | 212 |
/// \pre \ref findMinMean() must be called before using this function. |
205 | 213 |
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; ) { |
|
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; ) { |
|
210 | 218 |
_cycle_path->addBack(_policy[v]); |
211 | 219 |
} |
212 | 220 |
return true; |
... | ... |
@@ -225,36 +233,36 @@ |
225 | 233 |
/// |
226 | 234 |
/// This function returns the total length of the found cycle. |
227 | 235 |
/// |
228 |
/// \pre \ref run() or \ref |
|
236 |
/// \pre \ref run() or \ref findMinMean() must be called before |
|
229 | 237 |
/// using this function. |
230 | 238 |
Value cycleLength() const { |
231 |
return |
|
239 |
return _best_length; |
|
232 | 240 |
} |
233 | 241 |
|
234 | 242 |
/// \brief Return the number of arcs on the found cycle. |
235 | 243 |
/// |
236 | 244 |
/// This function returns the number of arcs on the found cycle. |
237 | 245 |
/// |
238 |
/// \pre \ref run() or \ref |
|
246 |
/// \pre \ref run() or \ref findMinMean() must be called before |
|
239 | 247 |
/// using this function. |
240 | 248 |
int cycleArcNum() const { |
241 |
return |
|
249 |
return _best_size; |
|
242 | 250 |
} |
243 | 251 |
|
244 | 252 |
/// \brief Return the mean length of the found cycle. |
245 | 253 |
/// |
246 | 254 |
/// This function returns the mean length of the found cycle. |
247 | 255 |
/// |
248 |
/// \note <tt> |
|
256 |
/// \note <tt>alg.cycleMean()</tt> is just a shortcut of the |
|
249 | 257 |
/// following code. |
250 | 258 |
/// \code |
251 |
/// return double( |
|
259 |
/// return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum(); |
|
252 | 260 |
/// \endcode |
253 | 261 |
/// |
254 | 262 |
/// \pre \ref run() or \ref findMinMean() must be called before |
255 | 263 |
/// using this function. |
256 | 264 |
double cycleMean() const { |
257 |
return double( |
|
265 |
return static_cast<double>(_best_length) / _best_size; |
|
258 | 266 |
} |
259 | 267 |
|
260 | 268 |
/// \brief Return the found cycle. |
... | ... |
@@ -274,153 +282,166 @@ |
274 | 282 |
|
275 | 283 |
private: |
276 | 284 |
|
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 |
|
|
285 |
// Initialize |
|
286 |
void init() { |
|
287 |
_tol.epsilon(1e-6); |
|
288 |
if (!_cycle_path) { |
|
289 |
_local_path = true; |
|
290 |
_cycle_path = new Path; |
|
286 | 291 |
} |
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 |
|
|
292 |
_queue.resize(countNodes(_gr)); |
|
293 |
_best_found = false; |
|
294 |
_best_length = 0; |
|
295 |
_best_size = 1; |
|
296 |
_cycle_path->clear(); |
|
297 |
} |
|
298 |
|
|
299 |
// Find strongly connected components and initialize _comp_nodes |
|
300 |
// and _in_arcs |
|
301 |
void findComponents() { |
|
302 |
_comp_num = stronglyConnectedComponents(_gr, _comp); |
|
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 |
} |
|
294 | 324 |
} |
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; |
|
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; |
|
299 | 334 |
} |
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; |
|
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 |
} |
|
308 | 349 |
} |
309 | 350 |
} |
310 | 351 |
return true; |
311 | 352 |
} |
312 | 353 |
|
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 |
|
|
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 cycles |
|
325 |
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 |
// A cycle is found |
|
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 |
|
|
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; |
|
347 | 376 |
} |
348 |
++path_cnt; |
|
349 |
} |
|
350 |
} |
|
351 |
return curr_cycle_found; |
|
352 |
} |
|
353 |
|
|
354 |
// Contract the policy graph to be connected by cutting all cycles |
|
355 |
// 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 search |
|
358 |
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 |
Node u, v; |
|
363 |
while (!queue.empty()) { |
|
364 |
v = queue.front(); queue.pop_front(); |
|
365 |
for (InArcIt e(_gr, v); e != INVALID; ++e) { |
|
366 |
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) { |
|
381 |
u = _gr.source(e); |
|
382 |
if (_comp[u] == comp && !found[u]) { |
|
383 |
found[u] = true; |
|
384 |
++found_cnt; |
|
385 |
_policy[u] = e; |
|
386 |
|
|
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; |
|
387 | 383 |
} |
388 | 384 |
} |
389 | 385 |
} |
390 | 386 |
} |
391 | 387 |
|
392 |
// Compute node distances in the policy graph and update the |
|
393 |
// policy graph if the node distances can be improved. |
|
388 |
// Contract the policy graph and compute node distances |
|
394 | 389 |
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 |
|
|
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; |
|
402 | 400 |
Node u, v; |
403 |
while (!queue.empty()) { |
|
404 |
v = queue.front(); queue.pop_front(); |
|
405 |
|
|
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]; |
|
406 | 406 |
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); |
|
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 | 411 |
} |
412 | 412 |
} |
413 | 413 |
} |
414 |
|
|
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]; |
|
422 |
u = _gr.source(e); |
|
423 |
if (!_reached[u]) { |
|
424 |
_reached[u] = true; |
|
425 |
_policy[u] = e; |
|
426 |
_dist[u] = _dist[v] + _length[e] - curr_mean; |
|
427 |
_queue[++_qback] = u; |
|
428 |
} |
|
429 |
} |
|
430 |
} |
|
431 |
|
|
432 |
// Improve node distances |
|
415 | 433 |
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; |
|
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]; |
|
438 |
u = _gr.source(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 |
} |
426 | 447 |
return improved; |
0 comments (0 inline)