NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
myers_banded_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 #include <nvbio/basic/types.h>
31 #include <nvbio/basic/numbers.h>
32 
33 namespace nvbio {
34 namespace aln {
35 
36 namespace priv {
37 
40 
41 
44 template <uint32 ALPHABET_SIZE>
45 struct MyersBitVectors
46 {
47  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
48  MyersBitVectors()
49  {
50  #pragma unroll ALPHABET_SIZE
51  for (uint32 c = 0; c < ALPHABET_SIZE; c++)
52  B[c] = 0u;
53  }
54 
55  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
56  void load(const uint8 c, const uint32 v) { B[c] |= v; }
57 
58  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
59  uint32 get(const uint8 c) const { return B[c]; }
60 
61  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
62  void shift()
63  {
64  #pragma unroll ALPHABET_SIZE
65  for(uint32 c = 0; c < ALPHABET_SIZE; c++)
66  B[c] >>= 1;
67  }
68 
69  uint32 B[ALPHABET_SIZE];
70 };
71 
74 template <>
75 struct MyersBitVectors<2>
76 {
77  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
78  MyersBitVectors()
79  {
80  #pragma unroll 2
81  for (uint32 c = 0; c < 2; c++)
82  B[c] = 0u;
83  }
84 
85  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
86  void load(const uint8 c, const uint32 v)
87  {
88  const uint8 cc = c & 1u;
89  if (cc == 0) B[0] |= v;
90  if (cc == 1) B[1] |= v;
91  }
92 
93  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
94  uint32 get(const uint8 c) const
95  {
96  const uint8 cc = c & 1u;
97  return cc == 0 ? B[0] : B[1];
98  }
99 
100  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
101  void shift()
102  {
103  #pragma unroll 2
104  for(uint32 c = 0; c < 2; c++)
105  B[c] >>= 1;
106  }
107 
108  uint32 B[4];
109 };
110 
113 template <>
114 struct MyersBitVectors<4>
115 {
116  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
117  MyersBitVectors()
118  {
119  #pragma unroll 4
120  for (uint32 c = 0; c < 4; c++)
121  B[c] = 0u;
122  }
123 
124  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
125  void load(const uint8 c, const uint32 v)
126  {
127  const uint8 cc = c & 3u;
128  if (cc == 0) B[0] |= v;
129  if (cc == 1) B[1] |= v;
130  if (cc == 2) B[2] |= v;
131  if (cc == 3) B[3] |= v;
132  }
133 
134  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
135  uint32 get(const uint8 c) const
136  {
137  const uint8 cc = c & 3u;
138  return cc <= 1 ?
139  cc == 0 ? B[0] : B[1] :
140  cc == 2 ? B[2] : B[3];
141  }
142 
143  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
144  void shift()
145  {
146  #pragma unroll 4
147  for(uint32 c = 0; c < 4; c++)
148  B[c] >>= 1;
149  }
150 
151  uint32 B[4];
152 };
153 
156 template <>
157 struct MyersBitVectors<5>
158 {
159  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
160  MyersBitVectors()
161  {
162  #pragma unroll 5
163  for (uint32 c = 0; c < 4; c++)
164  B[c] = 0u;
165  }
166 
167  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
168  void load(const uint8 c, const uint32 v)
169  {
170  if (c == 0) B[0] |= v;
171  if (c == 1) B[1] |= v;
172  if (c == 2) B[2] |= v;
173  if (c == 3) B[3] |= v;
174  if (c == 4) B[4] |= v;
175  }
176 
177  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
178  uint32 get(const uint8 c) const
179  {
180  return c <= 1 ? (c == 0 ? B[0] : B[1]) :
181  c <= 3 ? (c == 2 ? B[2] : B[3]) :
182  B[4];
183  }
184 
185  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
186  void shift()
187  {
188  #pragma unroll 5
189  for(uint32 c = 0; c < 5; c++)
190  B[c] >>= 1;
191  }
192 
193  uint32 B[5];
194 };
195 
196 template <uint32 BAND_WIDTH>
197 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
198 int diagonal_column(const uint32 B_c_ref, uint32& VP, uint32& VN)
199 {
200  uint32 X = B_c_ref | VN;
201  const uint32 D0 = ((VP + (X & VP)) ^ VP) | X;
202  const uint32 HN = VP & D0;
203  const uint32 HP = VN | ~ (VP | D0);
204  X = D0 >> 1;
205  VN = X & HP;
206  VP = HN | ~ (X | HP);
207 
208  return 1 - ((D0 >> (BAND_WIDTH - 1)) & 1u);
209 }
210 
211 template <uint32 BAND_WIDTH>
212 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
213 int horizontal_column(const uint32 B_c_ref, uint32& VP, uint32& VN, const int s)
214 {
215  uint32 X = B_c_ref | VN;
216  const uint32 D0 = ((VP + (X & VP)) ^ VP) | X;
217  const uint32 HN = VP & D0;
218  const uint32 HP = VN | ~ (VP | D0);
219  X = D0 >> 1;
220  VN = X & HP;
221  VP = HN | ~ (X | HP);
222 
223  return ((HP >> s) & 1) - ((HN >> s) & 1);
224 }
225 
226 template <
227  uint32 BAND_WIDTH,
228  uint32 C,
229  AlignmentType TYPE,
230  uint32 ALPHABET_SIZE,
231  typename pattern_string,
232  typename text_string,
233  typename sink_type>
234 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
235 bool banded_myers(
236  const pattern_string pattern,
237  const text_string text,
238  const int16 min_score,
239  sink_type& sink)
240 {
241  const uint32 pattern_len = pattern.length();
242  const uint32 text_len = text.length();
243 
244  // check whether the text is long enough
245  if (text_len < pattern_len)
246  return false;
247 
248  MyersBitVectors<ALPHABET_SIZE> B;
249 
250  uint32 VP = 0xFFFFFFFFu;
251  uint32 VN = 0x0u;
252 
253  int dist = -int( C );
254 
255  // initialize the column with C bases from the pattern
256  #pragma unroll 16
257  for (int j = 0; j < C; j++)
258  B.load( pattern[j], 1u << (BAND_WIDTH - C + j) );
259 
260  // phase 1 (diagonal)
261  const uint32 last = nvbio::min( text_len - 1u, pattern_len - C );
262  for(uint32 i = 0; i < last; ++i)
263  {
264  B.shift();
265 
266  B.load( pattern[ C + i ], 1u << (BAND_WIDTH - 1) );
267 
268  dist -= diagonal_column<BAND_WIDTH>( B.get( text[i] ), VP, VN );
269  }
270 
271  // phase 2 (horizontal)
272  int s = BAND_WIDTH - 1u + pattern_len - C - last;
273  for (uint32 i = last; i < text_len && s >= 0; ++i)
274  {
275  B.shift();
276 
277  dist -= horizontal_column<BAND_WIDTH>( B.get( text[i] ), VP, VN, s );
278 
279  // report a potential hit
280  if (TYPE == SEMI_GLOBAL && dist >= min_score)
281  sink.report( dist, make_uint2( i+1, pattern_len ) );
282 
283  --s;
284  }
285  if (TYPE == GLOBAL && dist >= min_score)
286  sink.report( dist, make_uint2( text_len, pattern_len ) );
287  return true;
288 }
289 
302 template <
303  uint32 BAND_LEN,
304  AlignmentType TYPE,
305  uint32 ALPHABET_SIZE,
306  typename pattern_type,
307  typename qual_type,
308  typename text_type,
309  typename sink_type>
310 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
311 bool banded_alignment_score(
312  const EditDistanceAligner<TYPE, MyersTag<ALPHABET_SIZE> >& aligner,
313  pattern_type pattern,
314  qual_type quals,
315  text_type text,
316  const int32 min_score,
317  sink_type& sink)
318 {
319  return banded_myers<BAND_LEN,0u,TYPE,ALPHABET_SIZE>(
320  pattern,
321  text,
322  min_score,
323  sink );
324 }
325 
339 template <
340  uint32 BAND_LEN,
341  AlignmentType TYPE,
342  uint32 ALPHABET_SIZE,
343  typename pattern_type,
344  typename qual_type,
345  typename text_type,
346  typename sink_type,
347  typename checkpoint_type>
348 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
349 bool banded_alignment_score(
350  const EditDistanceAligner<TYPE, MyersTag<ALPHABET_SIZE> >& aligner,
351  pattern_type pattern,
352  qual_type quals,
353  text_type text,
354  const int32 min_score,
355  const uint32 window_begin,
356  const uint32 window_end,
357  sink_type& sink,
358  checkpoint_type checkpoint)
359 {
360  return banded_myers<BAND_LEN,0u,TYPE,ALPHABET_SIZE>(
361  pattern,
362  text,
363  min_score,
364  sink );
365 }
366 
368 
369 } // namespace priv
370 
371 } // namespace aln
372 } // namespace nvbio
373 
Generated on Wed Feb 25 2015 08:32:53 for NVBIO by   doxygen 1.8.4