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 169 insertions and 148 deletions:
↑ Collapse diff ↑
Show white space 48 line context
... ...
@@ -53,83 +53,89 @@
53 53
             typename LEN = typename GR::template ArcMap<int> >
54 54
#endif
55 55
  class MinMeanCycle
56 56
  {
57 57
  public:
58 58
  
59 59
    /// The type of the digraph the algorithm runs on
60 60
    typedef GR Digraph;
61 61
    /// The type of the length map
62 62
    typedef LEN LengthMap;
63 63
    /// The type of the arc lengths
64 64
    typedef typename LengthMap::Value Value;
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;
92
    // Data for storing the strongly connected components
93
    int _comp_num;
94
    typename Digraph::template NodeMap<int> _comp;
95
    std::vector<std::vector<Node> > _comp_nodes;
96
    std::vector<Node>* _nodes;
97
    typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;
91 98

	
92
    typename Digraph::template NodeMap<int> _comp;
93
    int _comp_num;
99
    // Queue used for BFS search
100
    std::vector<Node> _queue;
101
    int _qfront, _qback;
94 102

	
95
    std::vector<Node> _nodes;
96
    std::vector<Arc> _arcs;
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
    ///
124 130
    /// If you don't call this function before calling \ref run() or
125 131
    /// \ref findMinMean(), it will allocate a local \ref Path "path"
126 132
    /// structure. The destuctor deallocates this automatically
127 133
    /// allocated object, of course.
128 134
    ///
129 135
    /// \note The algorithm calls only the \ref lemon::Path::addBack()
130 136
    /// "addBack()" function of the given path structure.
131 137
    ///
132 138
    /// \return <tt>(*this)</tt>
133 139
    ///
134 140
    /// \sa cycle()
135 141
    MinMeanCycle& cyclePath(Path &path) {
... ...
@@ -151,285 +157,300 @@
151 157

	
152 158
    /// \brief Run the algorithm.
153 159
    ///
154 160
    /// This function runs the algorithm.
155 161
    /// It can be called more than once (e.g. if the underlying digraph
156 162
    /// and/or the arc lengths have been modified).
157 163
    ///
158 164
    /// \return \c true if a directed cycle exists in the digraph.
159 165
    ///
160 166
    /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
161 167
    /// \code
162 168
    ///   return mmc.findMinMean() && mmc.findCycle();
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;
181
      // Initialize and find strongly connected components
182
      init();
183
      findComponents();
183 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;
193 200
      }
194
      return _cycle_found;
201
      }
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();
285
    // Initialize
286
    void init() {
287
      _tol.epsilon(1e-6);
288
      if (!_cycle_path) {
289
        _local_path = true;
290
        _cycle_path = new Path;
291
      }
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();
284 306
      for (NodeIt n(_gr); n != INVALID; ++n) {
285
        if (_comp[n] == comp) _nodes.push_back(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);
286 311
      }
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 312
      }
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;
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);
299 322
      }
300
      Node u; Arc e;
301
      for (int j = 0; j < int(_arcs.size()); ++j) {
302
        e = _arcs[j];
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];
303 344
        u = _gr.source(e);
304
        if (!_reached[u] || _length[e] < _dist[u]) {
345
          if (_length[e] < _dist[u]) {
305 346
          _dist[u] = _length[e];
306 347
          _policy[u] = e;
307
          _reached[u] = true;
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) {
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) {
332 370
            // A cycle is found
333
            curr_cycle_found = true;
334 371
            clength = _length[_policy[u]];
335 372
            csize = 1;
336 373
            for (v = u; (v = _gr.target(_policy[v])) != u; ) {
337 374
              clength += _length[_policy[v]];
338 375
              ++csize;
339 376
            }
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;
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;
346 383
            }
347 384
          }
348
          ++path_cnt;
349 385
        }
350 386
      }
351
      return curr_cycle_found;
352
    }
353 387

	
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;
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) {
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);
387
          }
407
          if (_policy[u] == e && !_reached[u]) {
408
            _reached[u] = true;
409
            _dist[u] = _dist[v] + _length[e] - curr_mean;
410
            _queue[++_qback] = u;
388 411
        }
389 412
      }
390 413
    }
391 414

	
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) {
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];
406 422
          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);
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;
411 428
          }
412 429
        }
413 430
      }
414
      // Improving node distances
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;
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;
420 440
        if (_tol.less(delta, _dist[u])) {
421
          improved = true;
422 441
          _dist[u] = delta;
423 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)