NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dcs.h
Go to the documentation of this file.
1 /*
2  * nvbio
3  * Copyright (c) 2011-2014, NVIDIA CORPORATION. All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  * * Redistributions of source code must retain the above copyright
8  * notice, this list of conditions and the following disclaimer.
9  * * Redistributions in binary form must reproduce the above copyright
10  * notice, this list of conditions and the following disclaimer in the
11  * documentation and/or other materials provided with the distribution.
12  * * Neither the name of the NVIDIA CORPORATION nor the
13  * names of its contributors may be used to endorse or promote products
14  * derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL NVIDIA CORPORATION BE LIABLE FOR ANY
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  */
27 
28 #pragma once
29 
30 #include <nvbio/sufsort/sufsort_priv.h>
31 #include <nvbio/strings/string_set.h>
32 #include <nvbio/basic/thrust_view.h>
33 #include <thrust/host_vector.h>
34 #include <thrust/device_vector.h>
35 #include <thrust/iterator/constant_iterator.h>
36 #include <thrust/iterator/counting_iterator.h>
37 
38 namespace nvbio {
39 
40 
41 // Precomputed Difference Covers
42 template <uint32 Q> struct DCTable {};
43 
44 // Precomputed DC-64
45 template <> struct DCTable<64>
46 {
47  static const uint32 N = 9; // DC quorum
48 
49  static const uint32* S()
50  {
51  static const uint32 dc[9] = { 1, 2, 3, 6, 15, 17, 35, 43, 60 };
52  return dc;
53  }
54 };
55 // Precomputed DC-128
56 template <> struct DCTable<128>
57 {
58  static const uint32 N = 16; // DC quorum
59 
60  static const uint32* S()
61  {
62  static const uint32 dc[16] = { 0, 1, 2, 5, 10, 15, 26, 37, 48, 59, 70, 76, 82, 88, 89, 90 };
63  return dc;
64  }
65 };
66 // Precomputed DC-256
67 template <> struct DCTable<256>
68 {
69  static const uint32 N = 22; // DC quorum
70 
71  static const uint32* S()
72  {
73  static const uint32 dc[22] = { 0, 1, 2, 3, 7, 14, 21, 28, 43, 58, 73, 88, 103, 118, 133, 141, 149, 157, 165, 166, 167, 168 };
74  return dc;
75  }
76 };
77 // Precomputed DC-512
78 template <> struct DCTable<512>
79 {
80  static const uint32 N = 28; // DC quorum
81 
82  static const uint32* S()
83  {
84  static const uint32 dc[28] = { 0, 1, 2, 3, 4, 9, 18, 27, 36, 45, 64, 83, 102, 121, 140, 159, 178, 197, 216, 226, 236, 246, 256, 266, 267, 268, 269, 270 };
85  return dc;
86  }
87 };
88 // Precomputed DC-1024
89 template <> struct DCTable<1024>
90 {
91  static const uint32 N = 40; // DC quorum
92 
93  static const uint32* S()
94  {
95  static const uint32 dc[40] = { 0, 1, 2, 3, 4, 5, 6, 13, 26, 39, 52, 65, 78, 91, 118, 145, 172, 199, 226, 253, 280, 307, 334, 361, 388, 415, 442, 456, 470, 484, 498, 512, 526, 540, 541, 542, 543, 544, 545, 546 };
96  return dc;
97  }
98 };
99 // Precomputed DC-2048
100 template <> struct DCTable<2048>
101 {
102  static const uint32 N = 58; // DC quorum
103 
104  static const uint32* S()
105  {
106  static const uint32 dc[58] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 19, 38, 57, 76, 95, 114, 133, 152, 171, 190, 229, 268, 307, 346, 385, 424, 463, 502, 541, 580, 619, 658, 697, 736, 775, 814, 853, 892, 931, 951, 971, 991, 1011, 1031, 1051, 1071, 1091, 1111, 1131, 1132, 1133, 1134, 1135, 1136, 1137, 1138, 1139, 1140 };
107  return dc;
108  }
109 };
110 
113 struct DCSView
114 {
115  DCSView(
116  const uint32 _Q = 0,
117  const uint32 _N = 0,
118  const uint32 _size = 0,
119  uint32* _dc = NULL,
120  uint32* _lut = NULL,
121  uint32* _pos = NULL,
122  uint8* _bitmask = NULL,
123  uint32* _ranks = NULL) :
124  Q ( _Q ),
125  N ( _N ),
126  dc ( _dc ),
127  lut ( _lut ),
128  pos ( _pos ),
129  bitmask ( _bitmask ),
130  ranks ( _ranks ),
131  size ( _size ) {}
132 
133 
136  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
137  uint32 index(const uint32 i) const;
138 
139  const uint32 Q;
140  const uint32 N;
141  uint32* dc;
142  uint32* lut;
143  uint32* pos;
144  uint8* bitmask;
145  uint32* ranks;
146  const uint32 size;
147 };
148 
151 struct DCS
152 {
153  typedef DCSView plain_view_type;
154 
157  template <uint32 QT>
158  void init();
159 
162  template <uint32 QT>
163  static uint32 estimated_sample_size(const uint64 string_len) { return uint32( util::divide_ri( string_len * DCTable<QT>::N, QT ) + 1u ); }
164 
167  uint32 estimate_sample_size(const uint64 string_len) const { return uint32( util::divide_ri( string_len * N, Q ) + 1u ); }
168 
169  uint32 Q;
170  uint32 N;
171 
172  thrust::device_vector<uint32> d_dc;
173  thrust::device_vector<uint32> d_lut;
174  thrust::device_vector<uint32> d_pos;
175  thrust::device_vector<uint8> d_bitmask;
176  thrust::device_vector<uint32> d_ranks;
177 };
178 
181 inline DCSView plain_view(DCS& dcs)
182 {
183  return DCSView(
184  dcs.Q,
185  dcs.N,
186  uint32( dcs.d_ranks.size() ),
187  nvbio::plain_view( dcs.d_dc ),
188  nvbio::plain_view( dcs.d_lut ),
189  nvbio::plain_view( dcs.d_pos ),
190  nvbio::plain_view( dcs.d_bitmask ),
191  nvbio::plain_view( dcs.d_ranks ) );
192 }
193 
196 inline DCSView plain_view(const DCS& dcs)
197 {
198  return plain_view( const_cast<DCS&>( dcs ) );
199 }
200 
201 namespace priv {
202 
205 struct DCS_predicate
206 {
207  typedef uint32 argument_type;
208  typedef uint32 result_type;
209 
212  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
213  DCS_predicate(const uint32 _Q, const uint8* _dc_bitmask) : Q(_Q), dc_bitmask(_dc_bitmask) {}
214 
217  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
218  result_type operator() (const uint32 suffix) const { return (dc_bitmask[ suffix & (Q-1) ] != 0); }
219 
220  const uint32 Q;
221  const uint8* dc_bitmask;
222 };
223 
226 struct DCS_string_suffix_index
227 {
228  typedef uint32 argument_type;
229  typedef uint32 result_type;
230 
231  DCS_string_suffix_index(const DCSView _dcs) : dcs( _dcs ) {}
232 
235  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
236  result_type operator() (const uint32 suffix_idx) const
237  {
238  return nvbio::min( dcs.index( suffix_idx ), dcs.size );
239  }
240 
241  const DCSView dcs;
242 };
243 
246 template <uint32 SYMBOL_SIZE, typename string_type>
247 struct DCS_string_suffix_less
248 {
249  typedef uint32 first_argument_type;
250  typedef uint32 second_argument_type;
251  typedef uint32 result_type;
252 
255  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
256  DCS_string_suffix_less(
257  const uint64 _string_len,
258  const string_type _string,
259  const DCSView _dcs) :
260  string_len(_string_len),
261  string(_string),
262  dcs( _dcs ) {}
263 
266  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
267  result_type operator() (const uint64 suffix_idx1, const uint64 suffix_idx2) const
268  {
269  const uint32 Q = dcs.Q;
270 
271  //#define DCS_CHECKS
272  #if defined(DCS_CHECKS)
273  const string_suffix_less<SYMBOL_SIZE,string_type> less( string_len, string );
274  const bool r = less( suffix_idx1, suffix_idx2 );
275  #endif
276 
277  const uint32 WORD_BITS = 32u; // use 32-bit words
278  const uint32 DOLLAR_BITS = 4u; // 4 is the minimum number needed to encode up to 16 symbols per word
279  const uint32 SYMBOLS_PER_WORD = symbols_per_word<SYMBOL_SIZE,WORD_BITS,DOLLAR_BITS>();
280 
281  const uint32 suffix_len1 = string_len - suffix_idx1;
282  const uint32 suffix_len2 = string_len - suffix_idx2;
283  const uint32 q_words = (Q + SYMBOLS_PER_WORD-1) / SYMBOLS_PER_WORD;
284  const uint32 n_words = nvbio::min(
285  uint32( nvbio::min(
286  suffix_len1,
287  suffix_len2 ) + SYMBOLS_PER_WORD-1 ) / SYMBOLS_PER_WORD,
288  q_words );
289 
290  // loop through all string-words
291  for (uint32 w = 0; w < n_words; ++w)
292  {
293  string_suffix_word_functor<SYMBOL_SIZE,WORD_BITS,DOLLAR_BITS,string_type,uint32> word_functor( string_len, string, w );
294 
295  const uint32 w1 = word_functor( suffix_idx1 );
296  const uint32 w2 = word_functor( suffix_idx2 );
297  if (w1 < w2) return true;
298  if (w1 > w2) return false;
299  }
300 
301  // check whether the suffixes are shorter than Q - this should never happen...
302  if (suffix_len1 < Q ||
303  suffix_len2 < Q)
304  {
305  #if defined(DCS_CHECKS)
306  const bool r2 = suffix_len1 < suffix_len2;
307  if (r != r2)
308  printf("short suffixes %u, %u\n", suffix_len1, suffix_len2 );
309  #endif
310 
311  return suffix_len1 < suffix_len2;
312  }
313 
314  // compare the DCS ranks
315  {
316  const uint32 i_mod_Q = suffix_idx1 & (Q-1);
317  const uint32 j_mod_Q = suffix_idx2 & (Q-1);
318 
319  // lookup the smallest number l such that (i + l) and (j + l) are in the DCS
320  const uint32 l = dcs.lut[ i_mod_Q * Q + j_mod_Q ];
321 
322  // by construction (suffix_idx1 + l) and (suffix_idx2 + l) are both in the DCS,
323  // we just need to find exactly where...
324  const uint32 pos_i = dcs.index( suffix_idx1 + l );
325  const uint32 pos_j = dcs.index( suffix_idx2 + l );
326 
327  // now we can lookup the ranks of the suffixes in the DCS
328  const uint32 rank_i = dcs.ranks[ pos_i ];
329  const uint32 rank_j = dcs.ranks[ pos_j ];
330 
331  #if defined(DCS_CHECKS)
332  const bool r2 = rank_i < rank_j;
333  if (r != r2)
334  {
335  printf("(%u,%u) : %u != %u, l[%u], pos[%u,%u], rank[%u,%u]\n",
336  (uint32)suffix_idx1, (uint32)suffix_idx2,
337  (uint32)r,
338  (uint32)r2,
339  (uint32)l,
340  (uint32)pos_i,
341  (uint32)pos_j,
342  (uint32)rank_i,
343  (uint32)rank_j);
344  }
345  #endif
346 
347  return rank_i < rank_j;
348  }
349  }
350 
351  const uint64 string_len;
352  const string_type string;
353  const DCSView dcs;
354 };
355 
356 } // namespace priv
357 } // namespace nvbio
358 
359 #include <nvbio/sufsort/dcs_inl.h>
Generated on Wed Feb 25 2015 08:33:03 for NVBIO by   doxygen 1.8.4