NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
utils_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/packedstream.h>
31 #include <nvbio/alignment/sink.h>
32 
33 namespace nvbio {
34 namespace aln {
35 
36 //
37 // Calculate the maximum possible number of pattern gaps that could occur in a
38 // given score boundary
39 //
40 template <AlignmentType TYPE, typename algorithm_tag>
41 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
42 uint32 max_pattern_gaps(
43  const EditDistanceAligner<TYPE,algorithm_tag>& aligner,
44  int32 min_score,
45  int32 pattern_len)
46 {
47  return -min_score;
48 }
49 
50 //
51 // Calculate the maximum possible number of reference gaps that could occur in a
52 // given score boundary
53 //
54 template <AlignmentType TYPE, typename algorithm_tag>
55 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
56 uint32 max_text_gaps(
57  const EditDistanceAligner<TYPE,algorithm_tag>& aligner,
58  int32 min_score,
59  int32 pattern_len)
60 {
61  return -min_score;
62 }
63 
64 //
65 // Calculate the maximum possible number of pattern gaps that could occur in a
66 // given score boundary
67 //
68 template <AlignmentType TYPE, typename scoring_scheme_type, typename algorithm_tag>
69 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
70 uint32 max_pattern_gaps(
71  const HammingDistanceAligner<TYPE,scoring_scheme_type,algorithm_tag>& aligner,
72  int32 min_score,
73  int32 pattern_len)
74 {
75  return 0;
76 }
77 
78 //
79 // Calculate the maximum possible number of reference gaps that could occur in a
80 // given score boundary
81 //
82 template <AlignmentType TYPE, typename scoring_scheme_type, typename algorithm_tag>
83 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
84 uint32 max_text_gaps(
85  const HammingDistanceAligner<TYPE,scoring_scheme_type,algorithm_tag>& aligner,
86  int32 min_score,
87  int32 pattern_len)
88 {
89  return 0;
90 }
91 
92 //
93 // Calculate the maximum possible number of pattern gaps that could occur in a
94 // given score boundary
95 //
96 template <AlignmentType TYPE, typename scoring_scheme_type, typename algorithm_tag>
97 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
98 uint32 max_pattern_gaps(
99  const SmithWatermanAligner<TYPE,scoring_scheme_type,algorithm_tag>& scoring,
100  int32 min_score,
101  int32 pattern_len)
102 {
103  // compute the optimal score
104  int32 score = pattern_len * scoring.scheme.match(30);
105  if (score < min_score)
106  return 0u;
107 
108  uint32 n = 0;
109  while (score >= min_score && n < pattern_len)
110  {
111  // subtract just the extension penalty
112  score += scoring.scheme.deletion();
113 
114  ++n;
115  }
116  return n-1;
117 }
118 
119 //
120 // Calculate the maximum possible number of reference gaps that could occur in a
121 // given score boundary
122 //
123 template <AlignmentType TYPE, typename scoring_scheme_type, typename algorithm_tag>
124 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
125 uint32 max_text_gaps(
126  const SmithWatermanAligner<TYPE,scoring_scheme_type,algorithm_tag>& scoring,
127  int32 min_score,
128  int32 pattern_len)
129 {
130  // compute the optimal score
131  int32 score = pattern_len * scoring.scheme.match(30);
132  if (score < min_score)
133  return 0u;
134 
135  uint32 n = 0;
136  while (score >= min_score && n < pattern_len)
137  {
138  // subtract just the extension penalty
139  score += scoring.scheme.insertion();
140 
141  ++n;
142  }
143  return n-1;
144 }
145 
146 //
147 // Calculate the maximum possible number of pattern gaps that could occur in a
148 // given score boundary
149 //
150 template <AlignmentType TYPE, typename scoring_scheme_type, typename algorithm_tag>
151 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
152 uint32 max_pattern_gaps(
153  const GotohAligner<TYPE,scoring_scheme_type,algorithm_tag>& scoring,
154  int32 min_score,
155  int32 pattern_len)
156 {
157  // compute the optimal score
158  int32 score = pattern_len * scoring.scheme.match(30);
159  if (score < min_score)
160  return 0u;
161 
162  // subtract the gap open penalty
163  score += scoring.scheme.pattern_gap_open();
164 
165  uint32 n = 0;
166  while (score >= min_score && n < pattern_len)
167  {
168  // subtract just the extension penalty
169  score += scoring.scheme.pattern_gap_extension();
170 
171  ++n;
172  }
173  return n-1;
174 }
175 
176 //
177 // Calculate the maximum possible number of reference gaps that could occur in a
178 // given score boundary
179 //
180 template <AlignmentType TYPE, typename scoring_scheme_type, typename algorithm_tag>
181 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
182 uint32 max_text_gaps(
183  const GotohAligner<TYPE,scoring_scheme_type,algorithm_tag>& scoring,
184  int32 min_score,
185  int32 pattern_len)
186 {
187  // compute the optimal score
188  int32 score = pattern_len * scoring.scheme.match(30);
189  if (score < min_score)
190  return 0u;
191 
192  // subtract the gap open penalty
193  score += scoring.scheme.text_gap_open();
194 
195  uint32 n = 0;
196  while (score >= min_score && n < pattern_len)
197  {
198  // subtract just the extension penalty
199  score += scoring.scheme.text_gap_extension();
200 
201  ++n;
202  }
203  return n-1;
204 }
205 
206 template <uint32 BAND_LEN, typename score_type>
207 struct select_dispatch
208 {
209  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
210  static score_type enact(const uint32 M, const score_type* band)
211  {
212  // get the highest score along the long edge of the path graph
213  score_type r = 0;
214 
215  const uint32 M_mod = (M-1) & (BAND_LEN-1);
216 
217  #pragma unroll
218  for (uint32 j = 1; j <= BAND_LEN; ++j)
219  {
220  const score_type value = band[j];
221  if (j == M_mod+1)
222  r = value;
223  }
224  return r;
225  }
226 };
227 template <uint32 BAND_LEN>
228 struct select_dispatch<BAND_LEN, simd4u8>
229 {
230  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
231  static simd4u8 enact(const uint4 M, const simd4u8* band)
232  {
233  // get the highest score along the long edge of the path graph
234  uint4 r;
235 
236  const uint4 M_mod = make_uint4(
237  (M.x-1) & (BAND_LEN-1),
238  (M.y-1) & (BAND_LEN-1),
239  (M.z-1) & (BAND_LEN-1),
240  (M.w-1) & (BAND_LEN-1) );
241 
242  #pragma unroll
243  for (uint32 j = 1; j <= BAND_LEN; ++j)
244  {
245  const simd4u8 value = band[j];
246  if (j == M_mod.x+1) r.x = get<0>( value );
247  if (j == M_mod.y+1) r.y = get<1>( value );
248  if (j == M_mod.z+1) r.z = get<2>( value );
249  if (j == M_mod.w+1) r.w = get<3>( value );
250  }
251  return simd4u8( r );
252  }
253 };
254 
255 template <uint32 BAND_LEN>
256 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
257 void save_Mth(const uint4 M, const simd4u8* band, simd4u8& best_score)
258 {
259  const simd4u8 m = select_dispatch<BAND_LEN,simd4u8>::enact( M, band );
260  best_score = nvbio::max( best_score, m );
261 }
262 
263 template <uint32 BAND_LEN, typename sink_type>
264 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
265 void save_Mth(const uint4 M, const simd4u8* band, const uint32 i, sink_type& sink, const simd4u8 mask)
266 {
267  const simd4u8 m = select_dispatch<BAND_LEN,simd4u8>::enact( M, band );
268  sink.report( and_op( m, mask ) );
269 }
270 
271 template <uint32 BAND_LEN, typename score_type>
272 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
273 void save_Mth(const uint32 M, const score_type* band, score_type& best_score)
274 {
275  const score_type m = select_dispatch<BAND_LEN,score_type>::enact( M, band );
276  best_score = nvbio::max( best_score, m );
277 }
278 
279 template <uint32 BAND_LEN, typename score_type, typename sink_type>
280 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
281 void save_Mth(const uint32 M, const score_type* band, const uint32 i, sink_type& sink)
282 {
283  const score_type m = select_dispatch<BAND_LEN,score_type>::enact( M, band );
284  sink.report( m, make_uint2( i+1, M ) );
285 }
286 template <uint32 BAND_LEN, typename score_type>
287 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
288 void save_boundary(const uint32 block, const uint32 M, const score_type* band, score_type& best_score)
289 {
290  if (block + BAND_LEN >= M)
291  save_Mth<BAND_LEN>( M, band, best_score );
292 }
293 template <uint32 BAND_LEN, typename score_type, typename sink_type>
294 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
295 void save_boundary(const uint32 block, const uint32 M, const score_type* band, const uint32 i, sink_type& sink)
296 {
297  if (block + BAND_LEN >= M)
298  save_Mth<BAND_LEN>( M, band, i, sink );
299 }
300 
301 template <uint32 BAND_LEN, typename sink_type>
302 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
303 void save_boundary(const uint32 block, const uint4 M, const simd4u8* band, const uint32 i, sink_type& sink, const simd4u8 active_mask)
304 {
305  const simd4u8 mask = and_op( simd4u8(
306  block < M.x && block + BAND_LEN >= M.x ? 0xFFu : 0u ,
307  block < M.y && block + BAND_LEN >= M.y ? 0xFFu : 0u ,
308  block < M.z && block + BAND_LEN >= M.z ? 0xFFu : 0u ,
309  block < M.w && block + BAND_LEN >= M.w ? 0xFFu : 0u ),
310  active_mask );
311 
312  if (any( mask ))
313  save_Mth<BAND_LEN>( M, band, i, sink, mask );
314 }
315 
316 } // namespace aln
317 } // namespace nvbio
Generated on Wed Feb 25 2015 08:32:54 for NVBIO by   doxygen 1.8.4