A toolkit for working with phylogenetic data.
v0.18.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
range_minimum_query.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2017 Lucas Czech
4 
5  This program is free software: you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation, either version 3 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 
18  Contact:
19  Lucas Czech <lucas.czech@h-its.org>
20  Exelixis Lab, Heidelberg Institute for Theoretical Studies
21  Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
22 */
23 
32 
33 #include <assert.h>
34 #include <limits>
35 #include <stdexcept>
36 #include <string>
37 
38 namespace genesis {
39 namespace utils {
40 
41 // ================================================================================================
42 // Lookup Tables
43 // ================================================================================================
44 
48 const size_t RangeMinimumQuery::catalan_numbers_[17][17] = {
49  { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
50  { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 },
51  { 0, 0, 2, 5, 9, 14, 20, 27, 35, 44, 54, 65, 77, 90, 104, 119, 135 },
52  { 0, 0, 0, 5, 14, 28, 48, 75, 110, 154, 208, 273, 350, 440, 544, 663, 798 },
53  { 0, 0, 0, 0, 14, 42, 90, 165, 275, 429, 637, 910, 1260, 1700, 2244, 2907, 3705 },
54  { 0, 0, 0, 0, 0, 42, 132, 297, 572, 1001, 1638, 2548, 3808, 5508, 7752, 10659, 14364 },
55  { 0, 0, 0, 0, 0, 0, 132, 429, 1001, 2002, 3640, 6188, 9996, 15504, 23256, 33915, 48279 },
56  { 0, 0, 0, 0, 0, 0, 0, 429, 1430, 3432, 7072, 13260, 23256, 38760, 62016, 95931, 144210 },
57  { 0, 0, 0, 0, 0, 0, 0, 0, 1430, 4862, 11934, 25194, 48450, 87210, 149226, 245157, 389367 },
58  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 4862, 16796, 41990, 90440, 177650, 326876, 572033, 961400 },
59  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16796, 58786, 149226, 326876, 653752, 1225785, 2187185 },
60  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 58786, 208012, 534888, 1188640, 2414425, 4601610 },
61  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 208012, 742900, 1931540, 4345965, 8947575 },
62  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 742900, 2674440, 7020405, 15967980 },
63  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2674440, 9694845, 25662825 },
64  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9694845, 35357670 },
65  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 35357670 }
66 };
67 
72 const unsigned char RangeMinimumQuery::log_table_256_[256] =
73 {
74  0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
75  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
76  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
77  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
78  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
79  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
80  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
81  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
82  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
83  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
84  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
85  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
86  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
87  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
88  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
89  7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7
90 };
91 
95 const unsigned char RangeMinimumQuery::lsb_table_256_[256] =
96 {
97  0, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
98  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
99  5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
100  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
101  6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
102  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
103  5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
104  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
105  7, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
106  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
107  5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
108  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
109  6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
110  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
111  5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
112  4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
113 };
114 
115 const RangeMinimumQuery::SuccinctType RangeMinimumQuery::highest_bits_set_[8] =
116 {
117  // ~1, ~3, ~7, ~15, ~31, ~63, ~127
118  static_cast< RangeMinimumQuery::SuccinctType >( ~0 ),
119  static_cast< RangeMinimumQuery::SuccinctType >( ~1 ),
120  static_cast< RangeMinimumQuery::SuccinctType >( ~3 ),
121  static_cast< RangeMinimumQuery::SuccinctType >( ~7 ),
122  static_cast< RangeMinimumQuery::SuccinctType >( ~15 ),
123  static_cast< RangeMinimumQuery::SuccinctType >( ~31 ),
124  static_cast< RangeMinimumQuery::SuccinctType >( ~63 ),
125  static_cast< RangeMinimumQuery::SuccinctType >( ~127 )
126 };
127 
128 // ================================================================================================
129 // Constructor and Rule of Five
130 // ================================================================================================
131 
135 RangeMinimumQuery::RangeMinimumQuery( std::vector<RangeMinimumQuery::IntType> const& array )
136  : array_( array )
137 {
138  init_();
139 }
140 
144 RangeMinimumQuery::RangeMinimumQuery( std::vector<RangeMinimumQuery::IntType>&& array )
145  : array_( std::move( array ))
146 {
147  init_();
148 }
149 
150 // ================================================================================================
151 // Member Functions
152 // ================================================================================================
153 
163 size_t RangeMinimumQuery::query( size_t i, size_t j ) const
164 {
165  // Safety checks.
166  if( i >= array_.size() || j >= array_.size() || j < i ) {
167  throw std::invalid_argument(
168  "Invalid range minimum query with indices i==" +
169  std::to_string(i) + " and j==" + std::to_string(j) + "."
170  );
171  }
172 
173  // If the array is too small, we do not precompute the internal data structures,
174  // but want to use the naive approach.
175  if( naive_ ) {
176  size_t naive_min = j;
177  for( size_t x = i; x < j; x++ ) {
178  if( array_[x] < array_[naive_min] ) {
179  naive_min = x;
180  }
181  }
182  return naive_min;
183  }
184 
185  size_t mb_i = microblock_(i); // i's microblock
186  size_t mb_j = microblock_(j); // j's microblock
187  size_t min, min_i, min_j; // min: to be returned
188  size_t s_mi = mb_i * micro_size_; // start of i's microblock
189  size_t i_pos = i - s_mi; // pos. of i in its microblock
190 
191  if( mb_i == mb_j ) { // only one microblock-query
192  min_i = clearbits_(precomputed_queries_(block_types_[mb_i], j-s_mi), i_pos);
193  min = min_i == 0 ? j : s_mi + lsb_(min_i);
194  } else {
195  size_t b_i = block_(i); // i's block
196  size_t b_j = block_(j); // j's block
197  size_t s_mj = mb_j * micro_size_; // start of j's microblock
198  size_t j_pos = j - s_mj; // position of j in its microblock
199  min_i = clearbits_(precomputed_queries_(block_types_[mb_i], micro_size_-1), i_pos);
200  min = min_i == 0 ? s_mi + micro_size_ - 1 : s_mi + lsb_(min_i); // left in-microblock-query
201 
202  // right in-microblock-query
203  if( precomputed_queries_(block_types_[mb_j], j_pos) == 0 ) {
204  min_j = j;
205  } else {
206  min_j = s_mj + lsb_(precomputed_queries_(block_types_[mb_j], j_pos));
207  }
208  if( array_[min_j] < array_[min] ) {
209  min = min_j;
210  }
211 
212  if( mb_j > mb_i + 1 ) { // otherwise we're done!
213  size_t s_bi = b_i * block_size_; // start of block i
214  size_t s_bj = b_j * block_size_; // start of block j
215  if( s_bi+micro_size_ > i ) { // do another microblock-query!
216  mb_i++; // go one microblock to the right
217 
218  // right in-block-query
219  if( precomputed_queries_(block_types_[mb_i], micro_size_-1) == 0 ) {
220  min_i = s_bi + block_size_ - 1;
221  } else {
222  min_i = s_mi + micro_size_ + lsb_(
223  precomputed_queries_(block_types_[mb_i], micro_size_-1)
224  );
225  }
226  if( array_[min_i] < array_[min] ) {
227  min = min_i;
228  }
229  }
230  if( j >= s_bj+micro_size_ ) { // and yet another microblock-query!
231  mb_j--; // go one microblock to the left
232 
233  // right in-block-query
234  if( precomputed_queries_(block_types_[mb_j], micro_size_-1) == 0 ) {
235  min_j = s_mj - 1;
236  } else {
237  min_j = s_bj + lsb_(precomputed_queries_(block_types_[mb_j], micro_size_-1));
238  }
239  if( array_[min_j] < array_[min] ) {
240  min = min_j;
241  }
242  }
243 
244  size_t block_difference = b_j - b_i;
245  if( block_difference > 1 ) { // otherwise we're done!
246  size_t k, twotothek, block_tmp; // for index calculations in M and M'
247  b_i++; // block where out-of-block-query starts
248  if( s_bj - s_bi - block_size_ <= super_size_ ) { // just one out-of-block-query
249  k = log2fast_(block_difference - 2);
250  twotothek = 1 << k; // 2^k
251  i = m_(k, b_i);
252  j = m_(k, b_j-twotothek);
253  min_i = array_[i] <= array_[j] ? i : j;
254  }
255  else { // here we have to answer a superblock-query:
256  size_t sb_i = superblock_(i); // i's superblock
257  size_t sb_j = superblock_(j); // j's superblock
258 
259  block_tmp = block_((sb_i+1)*super_size_); // end of left out-of-block-query
260  k = log2fast_(block_tmp - b_i);
261  twotothek = 1 << k; // 2^k
262  i = m_(k, b_i);
263  j = m_(k, block_tmp+1-twotothek);
264  min_i = array_[i] <= array_[j] ? i : j;
265 
266  block_tmp = block_(sb_j*super_size_); // start of right out-of-block-query
267  k = log2fast_(b_j - block_tmp);
268  twotothek = 1 << k; // 2^k
269  block_tmp--; // going one block to the left doesn't harm and saves some tests
270  i = m_(k, block_tmp);
271  j = m_(k, b_j-twotothek);
272  min_j = array_[i] <= array_[j] ? i : j;
273 
274  if( array_[min_j] < array_[min_i] ) {
275  min_i = min_j;
276  }
277 
278  if( sb_j > sb_i + 1 ) { // finally, the superblock-query:
279  k = log2fast_(sb_j - sb_i - 2);
280  twotothek = 1 << k;
281  i = m_prime_(k, sb_i+1);
282  j = m_prime_(k, sb_j-twotothek);
283  min_j = array_[i] <= array_[j] ? i : j;
284  if (array_[min_j] < array_[min_i]) {
285  min_i = min_j; // does NOT always return leftmost min!!!
286  }
287  }
288 
289  }
290  if( array_[min_i] < array_[min] ) {
291  min = min_i; // does NOT always return leftmost min!!!
292  }
293  }
294  }
295  }
296 
297  return min;
298 }
299 
300 // ================================================================================================
301 // Internal Functions
302 // ================================================================================================
303 
307 void RangeMinimumQuery::init_()
308 {
309  micro_size_ = 1 << 3; // microblock-size
310  block_size_ = 1 << 4; // block-size
311  super_size_ = 1 << 8; // superblock-size
312 
313  auto const n = array_.size();
314 
315  // number of microblocks = array_.size() / micro_size_
316  // number of blocks = array_.size() / block_size_
317  // number of superblocks = array_.size() / super_size_
318  auto const micro_count = microblock_(n-1)+1;
319  auto const block_count = block_(n-1)+1;
320  auto const super_count = superblock_(n-1)+1;
321 
322  // The following is necessary because we've fixed s, s' and s'' according to the computer's
323  // word size and NOT according to the input size. This may cause the (super-)block-size
324  // to be too big, or, in other words, the array too small. If this code is compiled on
325  // a 32-bit computer, this happens iff n < 113. For such small instances it isn't
326  // advisable anyway to use the whole data structure, because simpler methods are faster and
327  // less space consuming.
328  // Thus, in those cases, we default to a naive implementation.
329  // We do not need the precomputation. Instead, query() will simply use the array directly.
330  if( block_count < super_size_ / (2*block_size_) ) {
331  naive_ = true;
332  return;
333  }
334 
335  // Type-calculation for the microblocks and pre-computation of in-microblock-queries:
336  block_types_ = std::vector<BlockTypeType>(micro_count);
337 
338  precomputed_queries_ = Matrix<SuccinctType>(
339  catalan_numbers_[micro_size_][micro_size_], micro_size_
340  );
341  for( size_t i = 0; i < catalan_numbers_[micro_size_][micro_size_]; i++ ) {
342  precomputed_queries_(i, 0) = 1; // init with impossible value
343  }
344 
345  std::vector<IntType> rp(micro_size_+1); // rp: rightmost path in Cart. tree
346  size_t z = 0; // index in array a
347  size_t start; // start of current block
348  size_t end; // end of current block
349  size_t q; // position in catalan_numbers_ triangle
350  size_t p; // --------- " ----------------
351  rp[0] = std::numeric_limits<IntType>::min(); // stopper (minus infinity)
352 
353  // prec[i]: the jth bit is 1 iff j is 1. pos. to the left of i where array_[j] < array_[i]
354  std::vector<size_t> gstack( micro_size_ );
355  size_t gstacksize;
356  size_t g; // first position to the left of i where array_[g[i]] < array_[i]
357 
358  for( size_t i = 0; i < micro_count; i++ ) { // step through microblocks
359  start = z; // init start
360  end = start + micro_size_; // end of block (not inclusive!)
361  if( end > n ) {
362  end = n; // last block could be smaller than s!
363  }
364 
365  // compute block type as in Fischer/Heun CPM'06:
366  q = micro_size_; // init q
367  p = micro_size_-1; // init p
368  block_types_[i] = 0; // init type (will be increased!)
369  rp[1] = array_[z]; // init rightmost path
370 
371  while( ++z < end ) { // step through current block:
372  p--;
373  while( rp[q-p-1] > array_[z] ) {
374  block_types_[i] += catalan_numbers_[p][q]; // update type
375  q--;
376  }
377  rp[q-p] = array_[z]; // add last element to rightmost path
378  }
379 
380  // precompute in-block-queries for this microblock (if necessary)
381  // as in Alstrup et al. SPAA'02:
382  if( precomputed_queries_(block_types_[i], 0) == 1 ) {
383  precomputed_queries_(block_types_[i], 0) = 0;
384  gstacksize = 0;
385  for( size_t j = start; j < end; j++ ) {
386  while( gstacksize > 0 && (array_[j] < array_[gstack[gstacksize-1]]) ) {
387  gstacksize--;
388  }
389  if( gstacksize > 0 ) {
390  g = gstack[gstacksize-1];
391  precomputed_queries_(block_types_[i], j-start)
392  = precomputed_queries_(block_types_[i], g-start) | (1 << (g % micro_size_));
393  } else {
394  precomputed_queries_(block_types_[i], j-start) = 0;
395  }
396  gstack[gstacksize++] = j;
397  }
398  }
399  }
400 
401  // space for out-of-block- and out-of-superblock-queries:
402  auto M_depth = static_cast<size_t>(
403  floor( log2( static_cast<double>( super_size_ ) / static_cast<double>( block_size_ ) ))
404  );
405  m_matrix_ = Matrix<SuccinctType>(M_depth, block_count);
406 
407  auto Mprime_depth = static_cast<size_t>( floor(log2(super_count)) ) + 1;
408  m_prime_ = Matrix<size_t>(Mprime_depth, super_count);
409 
410  // fill 0'th rows of M and Mprime:
411  z = 0; // minimum in current block
412  q = 0; // pos. of min in current superblock
413  g = 0; // number of current superblock
414  for( size_t i = 0; i < block_count; i++ ) { // step through blocks
415  start = z; // init start
416  p = start; // init minimum
417  end = start + block_size_; // end of block (not inclusive!)
418  if( end > n ) {
419  end = n; // last block could be smaller than block_size_!
420  }
421  if( array_[z] < array_[q] ) {
422  q = z; // update minimum in superblock
423  }
424 
425  while( ++z < end ) { // step through current block:
426  if( array_[z] < array_[p] ) {
427  p = z; // update minimum in block
428  }
429  if( array_[z] < array_[q] ) {
430  q = z; // update minimum in superblock
431  }
432  }
433  m_matrix_(0, i) = p-start; // store index of block-minimum (offset!)
434  if( z % super_size_ == 0 || z == n ) { // reached end of superblock?
435  m_prime_(0, g++) = q; // store index of superblock-minimum
436  q = z;
437  }
438  }
439 
440  // fill M:
441  size_t dist = 1; // always 2^(j-1)
442  for( size_t j = 1; j < M_depth; j++ ) {
443  for( size_t i = 0; i < block_count - dist; i++ ) { // be careful: loop may go too far
444  // add 'skipped' elements in a
445  if( array_[m_(j-1, i)] <= array_[m_(j-1,i+dist)] ) {
446  m_matrix_(j, i) = m_matrix_(j-1, i);
447  } else {
448  m_matrix_(j, i) = m_matrix_(j-1, i+dist) + (dist*block_size_);
449  }
450  }
451  for (size_t i = block_count - dist; i < block_count; i++) {
452  m_matrix_(j, i) = m_matrix_(j-1, i); // fill overhang
453  }
454  dist *= 2;
455  }
456 
457  // fill M':
458  dist = 1; // always 2^(j-1)
459  for( size_t j = 1; j < Mprime_depth; j++ ) {
460  for( size_t i = 0; i < super_count - dist; i++ ) {
461  // overhang
462  if( array_[m_prime_(j-1, i)] <= array_[m_prime_(j-1, i+dist)] ) {
463  m_prime_(j, i) = m_prime_(j-1, i);
464  } else {
465  m_prime_(j, i) = m_prime_(j-1, i+dist);
466  }
467  }
468  for (size_t i = super_count - dist; i < super_count; i++) {
469  m_prime_(j, i) = m_prime_(j-1, i);
470  }
471  dist *= 2;
472  }
473 }
474 
478 size_t RangeMinimumQuery::log2fast_( size_t v ) const
479 {
480  size_t c = 0; // c will be lg(v)
481  size_t t, tt; // temporaries
482 
483  if(( tt = v >> 16 )) {
484  c = (t = v >> 24) ? 24 + log_table_256_[t] : 16 + log_table_256_[tt & 0xFF];
485  } else {
486  c = (t = v >> 8) ? 8 + log_table_256_[t] : log_table_256_[v];
487  }
488  return c;
489 }
490 
494 RangeMinimumQuery::SuccinctType RangeMinimumQuery::clearbits_( SuccinctType n, size_t x ) const
495 {
496  return n & highest_bits_set_[x];
497 }
498 
499 } // namespace utils
500 } // namespace genesis
std::string to_string(T const &v)
Return a string representation of a given value.
Definition: string.hpp:300
size_t query(size_t i, size_t j) const
Main function for the Range Minimum Query.