gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Rework and fix the implementation of MinMeanCycle (#179) - Fix the handling of the cycle means. - Many implementation improvements: - More efficient data storage for the strongly connected components. - Better handling of BFS queues. - Merge consecutive BFS searches (perform two BFS searches instead of three). This version is about two times faster on average and an order of magnitude faster if there are a lot of strongly connected components.
0 1 0
default
1 file changed with 188 insertions and 167 deletions:
↑ Collapse diff ↑
Ignore white space 24 line context
... ...
@@ -65,59 +65,65 @@
65 65
    /// The type of the paths
66 66
    typedef lemon::Path<Digraph> Path;
67 67

	
68 68
  private:
69 69

	
70 70
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
71 71
  
72 72
    // The digraph the algorithm runs on
73 73
    const Digraph &_gr;
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
    // 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;
84
    bool _local_path;
83 85

	
84
    bool _local_path;
85
    bool _cycle_found;
86
    Node _cycle_node;
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:
100 106

	
101 107
    /// \brief Constructor.
102 108
    ///
103 109
    /// The constructor of the class.
104 110
    ///
105 111
    /// \param digraph The digraph the algorithm runs on.
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
      _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

	
114 120
    /// Destructor.
115 121
    ~MinMeanCycle() {
116 122
      if (_local_path) delete _cycle_path;
117 123
    }
118 124

	
119 125
    /// \brief Set the path structure for storing the found cycle.
120 126
    ///
121 127
    /// This function sets an external path structure for storing the
122 128
    /// found cycle.
123 129
    ///
... ...
@@ -163,273 +169,288 @@
163 169
    /// \endcode
164 170
    bool run() {
165 171
      return findMinMean() && findCycle();
166 172
    }
167 173

	
168 174
    /// \brief Find the minimum cycle mean.
169 175
    ///
170 176
    /// This function finds the minimum mean length of the directed
171 177
    /// cycles in the digraph.
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
        if (!initCurrentComponent(comp)) continue;
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 _cycle_found;
202
      return _best_found;
195 203
    }
196 204

	
197 205
    /// \brief Find a minimum mean directed cycle.
198 206
    ///
199 207
    /// This function finds a directed cycle of minimum mean length
200 208
    /// in the digraph using the data computed by findMinMean().
201 209
    ///
202 210
    /// \return \c true if a directed cycle exists in the digraph.
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;
213 221
    }
214 222

	
215 223
    /// @}
216 224

	
217 225
    /// \name Query Functions
218 226
    /// The results of the algorithm can be obtained using these
219 227
    /// functions.\n
220 228
    /// The algorithm should be executed before using them.
221 229

	
222 230
    /// @{
223 231

	
224 232
    /// \brief Return the total length of the found cycle.
225 233
    ///
226 234
    /// This function returns the total length of the found cycle.
227 235
    ///
228
    /// \pre \ref run() or \ref findCycle() must be called before
236
    /// \pre \ref run() or \ref findMinMean() must be called before
229 237
    /// using this function.
230 238
    Value cycleLength() const {
231
      return _cycle_length;
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 findCycle() must be called before
246
    /// \pre \ref run() or \ref findMinMean() must be called before
239 247
    /// using this function.
240 248
    int cycleArcNum() const {
241
      return _cycle_size;
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>mmc.cycleMean()</tt> is just a shortcut of the
256
    /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
249 257
    /// following code.
250 258
    /// \code
251
    ///   return double(mmc.cycleLength()) / mmc.cycleArcNum();
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(_cycle_length) / _cycle_size;
265
      return static_cast<double>(_best_length) / _best_size;
258 266
    }
259 267

	
260 268
    /// \brief Return the found cycle.
261 269
    ///
262 270
    /// This function returns a const reference to the path structure
263 271
    /// storing the found cycle.
264 272
    ///
265 273
    /// \pre \ref run() or \ref findCycle() must be called before using
266 274
    /// this function.
267 275
    ///
268 276
    /// \sa cyclePath()
269 277
    const Path& cycle() const {
270 278
      return *_cycle_path;
271 279
    }
272 280

	
273 281
    ///@}
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
        if (_comp[n] == comp) _nodes.push_back(n);
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
          _arcs.push_back(e);
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
      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 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
            queue.push_back(u);
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
      _dist[_cycle_node] = 0;
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
        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];
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
      // Improving node distances
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;
427 448
    }
428 449

	
429 450
  }; //class MinMeanCycle
430 451

	
431 452
  ///@}
432 453

	
433 454
} //namespace lemon
434 455

	
435 456
#endif //LEMON_MIN_MEAN_CYCLE_H
0 comments (0 inline)