NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dcs_inl.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 namespace nvbio {
31 
41 inline
42 void calculateColbournAndLingDC(const uint32 r, uint32& maxv, std::vector<uint32>& samples)
43 {
44  maxv = 24*r*r + 36*r + 13; // Corollary 2.3
45 
46  uint32 numsamp = 6*r + 4;
47 
48  // allocate enough storage
49  samples.resize( numsamp );
50 
51  // initialize samples to 0
52  for (uint32 i = 0; i < numsamp; ++i)
53  samples[i] = 0;
54 
55  // fill in the 1^r part of the B series
56  for (uint32 i = 1; i < r+1; i++)
57  samples[i] = samples[i-1] + 1;
58 
59  // fill in the (r + 1)^1 part
60  samples[r+1] = samples[r] + r + 1;
61 
62  // fill in the (2r + 1)^r part
63  for (uint32 i = r+2; i < r+2+r; i++)
64  samples[i] = samples[i-1] + 2*r + 1;
65 
66  // fill in the (4r + 3)^(2r + 1) part
67  for (uint32 i = r+2+r; i < r+2+r+2*r+1; i++)
68  samples[i] = samples[i-1] + 4*r + 3;
69 
70  // fill in the (2r + 2)^(r + 1) part
71  for (uint32 i = r+2+r+2*r+1; i < r+2+r+2*r+1+r+1; i++)
72  samples[i] = samples[i-1] + 2*r + 2;
73 
74  // fill in the last 1^r part
75  for (uint32 i = r+2+r+2*r+1+r+1; i < r+2+r+2*r+1+r+1+r; i++)
76  samples[i] = samples[i-1] + 1;
77 }
78 
79 // a utility SuffixHandler to rank the sorted suffixes
80 //
82 {
83  // constructor
84  //
85  DCSSuffixRanker(DCSView _dcs) : dcs( _dcs ), n_output(0) {}
86 
87  // process the next batch of suffixes
88  //
90  const uint32 n_suffixes,
91  const uint32* d_suffixes)
92  {
93  // essentially, invert the suffix array
94  thrust::scatter(
95  thrust::make_counting_iterator<uint32>( n_output ),
96  thrust::make_counting_iterator<uint32>( n_output ) + n_suffixes,
98  thrust::device_ptr<const uint32>( d_suffixes ),
99  priv::DCS_string_suffix_index( dcs ) ), // localize the index to the DCS
100  thrust::device_ptr<uint32>( dcs.ranks ) );
101 
102  n_output += n_suffixes;
103  }
104 
105  // process a sparse set of suffixes
106  //
108  const uint32 n_suffixes,
109  const uint32* d_suffixes,
110  const uint32* d_slots)
111  {
112  // essentially, invert the suffix array
113  thrust::scatter(
114  thrust::device_ptr<const uint32>( d_slots ),
115  thrust::device_ptr<const uint32>( d_slots ) + n_suffixes,
117  thrust::device_ptr<const uint32>( d_suffixes ),
118  priv::DCS_string_suffix_index( dcs ) ), // localize the index to the DCS
119  thrust::device_ptr<uint32>( dcs.ranks ) );
120  }
121 
122  const DCSView dcs;
124 };
125 
126 // return the sampled position of a given suffix index
127 //
130 {
131  const uint32 block_i = i / Q;
132  const uint32 mod_i = i & (Q-1);
133 
134  return block_i * N + pos[ mod_i ];
135 }
136 
137 // constructor
138 //
139 template <uint32 QT>
140 void DCS::init()
141 {
142  // build a table for our Difference Cover
143  const uint32* h_dc = DCTable<QT>::S();
144 
145  Q = QT;
146  N = DCTable<QT>::N;
147 
148  thrust::host_vector<uint8> h_bitmask( Q, 0u );
149  thrust::host_vector<uint32> h_lut( Q*Q, 0u );
150  thrust::host_vector<uint32> h_pos( Q, 0u );
151 
152  // build the DC bitmask
153  thrust::scatter(
154  thrust::make_constant_iterator<uint32>(1u),
155  thrust::make_constant_iterator<uint32>(1u) + N,
156  h_dc,
157  h_bitmask.begin() );
158 
159  // build the DC position table, mapping each entry in DC to its position (q -> i | DC[i] = q)
160  thrust::scatter(
161  thrust::make_counting_iterator<uint32>(0u),
162  thrust::make_counting_iterator<uint32>(0u) + N,
163  h_dc,
164  h_pos.begin() );
165 
166  // build the LUT (i,j) -> l | [(i + l) in DC && (j + l) in DC]
167  for (uint32 i = 0; i < Q; ++i)
168  {
169  for (uint32 j = 0; j < Q; ++j)
170  {
171  for (uint32 l = 0; l < Q; ++l)
172  {
173  if (h_bitmask[ (i + l) & (Q-1) ] &&
174  h_bitmask[ (j + l) & (Q-1) ])
175  {
176  h_lut[ i * Q + j ] = l;
177  break;
178  }
179  if (l == Q-1)
180  throw nvbio::logic_error("DCS: could not find a period for (%u,%u)!\n", i, j);
181  }
182  }
183  }
184 
185  // copy all tables to the device
186  d_dc.resize( N );
187  thrust::copy(
188  h_dc,
189  h_dc + N,
190  d_dc.begin() );
191 
192  d_lut = h_lut;
193  d_pos = h_pos;
194  d_bitmask = h_bitmask;
195 }
196 
197 } // namespace nvbio