NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
locate_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 
31 #include <nvbio/io/alignments.h>
34 
35 namespace nvbio {
36 namespace bowtie2 {
37 namespace cuda {
38 
39 namespace detail {
40 
43 
46 
49 
53 template <typename FMType, typename rFMType>
55 uint32 locate(const FMType& fmi, const rFMType& rfmi, const uint32 dir, const uint32 r)
56 {
58  {
59  const FMType fmi_ref = dir == FORWARD ? fmi : rfmi;
60  const uint32 g_pos = locate( fmi_ref, r );
61  NVBIO_CUDA_ASSERT( g_pos < fmi_ref.length() );
62  return (dir == REVERSE) ?
63  rfmi.length()-1 - g_pos : g_pos;
64  }
65  else
66  {
67  if (dir == FORWARD)
68  return locate(fmi, r);
69  else
70  return rfmi.length()-1 - locate(rfmi, r);
71  }
72 }
76 template <typename FMType, typename rFMType>
78 uint2 locate_init(const FMType& fmi, const rFMType& rfmi, const uint32 dir, const uint32 r)
79 {
81  {
82  const FMType fmi_ref = dir == FORWARD ? fmi : rfmi;
83  return locate_ssa_iterator( fmi_ref, r );
84  }
85  else
86  {
87  if (dir == FORWARD)
88  return locate_ssa_iterator( fmi, r );
89  else
90  return locate_ssa_iterator( rfmi, r );
91  }
92 }
93 
97 template <typename FMType, typename rFMType>
99 uint32 locate_lookup(const FMType& fmi, const rFMType& rfmi, const uint32 dir, const uint2 it)
100 {
102  {
103  const FMType fmi_ref = dir == FORWARD ? fmi : rfmi;
104  const uint32 g_pos = lookup_ssa_iterator( fmi_ref, it );
105  return (dir == REVERSE) ?
106  rfmi.length()-1 - g_pos : g_pos;
107  }
108  else
109  {
110  if (dir == FORWARD)
111  return lookup_ssa_iterator( fmi, it );
112  else
113  return rfmi.length()-1 - lookup_ssa_iterator( rfmi, it );
114  }
115 }
116 
122 template <typename BatchType, typename FMType, typename rFMType> __global__
124  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
125  const uint32 in_count,
126  const uint32* idx_queue,
127  HitQueuesDeviceView hits,
128  const ParamsPOD params)
129 {
130  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
131  if (thread_id >= in_count) return;
132 
133  const uint32 sa_idx = idx_queue[ thread_id ]; // fetch the sorting queue index
134 
135  HitReference<HitQueuesDeviceView> hit( hits, sa_idx );
136  const uint32 sa_pos = hit.loc; // fetch the SA coordinate
137  const packed_seed info = hit.seed; // fetch the attached info
138  const uint32 index_dir = info.index_dir; // decode the index direction
139  const uint32 pos_in_read = info.pos_in_read; // decode the seed's position in the read
140 
141  // locate the SA row and calculate the global position
142  const uint32 g_pos = locate( fmi, rfmi, index_dir, sa_pos ) - pos_in_read;
143 
144  // overwrite the locate queue with the final position
145  hit.loc = g_pos;
146 }
147 
153 template <typename BatchType, typename FMType, typename rFMType> __global__
155  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
156  const uint32 in_count,
157  const uint32* idx_queue,
158  HitQueuesDeviceView hits,
159  const ParamsPOD params)
160 {
161  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
162  if (thread_id >= in_count) return;
163 
164  const uint32 sa_idx = idx_queue[ thread_id ]; // fetch the sorting queue index
165 
166  HitReference<HitQueuesDeviceView> hit( hits, sa_idx );
167  const uint32 sa_pos = hit.loc; // fetch the SA coordinate
168  const packed_seed info = hit.seed; // fetch the attached info
169  const uint32 index_dir = info.index_dir; // decode the index direction
170 
171  // locate the SA row and calculate the global position
172  const uint2 ssa = locate_init( fmi, rfmi, index_dir, sa_pos );
173 
174  // overwrite the locate queue with the final position
175  hit.loc = ssa.x;
176  hit.ssa = ssa.y;
177 }
183 template <typename BatchType, typename FMType, typename rFMType> __global__
185  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
186  const uint32 in_count,
187  const uint32* idx_queue,
188  HitQueuesDeviceView hits,
189  const ParamsPOD params)
190 {
191  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
192  if (thread_id >= in_count) return;
193 
194  const uint32 sa_idx = idx_queue[ thread_id ]; // fetch the sorting queue index
195 
196  HitReference<HitQueuesDeviceView> hit( hits, sa_idx );
197  const uint32 sa_pos = hit.loc; // fetch the SA coordinate
198  const uint32 sa_off = hit.ssa; // fetch the SSA offset
199  const packed_seed info = hit.seed; // fetch the attached info
200  const uint32 index_dir = info.index_dir; // decode the index direction
201  const uint32 pos_in_read = info.pos_in_read; // decode the seed's position in the read
202 
203  // locate the SA row and calculate the global position
204  const uint32 g_pos = locate_lookup( fmi, rfmi, index_dir, make_uint2( sa_pos, sa_off ) ) - pos_in_read;
205 
206  // overwrite the locate queue with the final position
207  hit.loc = g_pos;
208 }
209 
213 template <typename index_iterator, typename flags_iterator> __global__
215  const uint32 in_count,
216  const uint32* idx_queue,
217  const uint32 reference_count,
218  const index_iterator reference_index,
219  HitQueuesDeviceView hits,
220  flags_iterator flags,
221  const ParamsPOD params)
222 {
223  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
224  if (thread_id >= in_count) return;
225 
226  const uint32 sa_idx = idx_queue[ thread_id ]; // fetch the sorting queue index
227 
228  HitReference<HitQueuesDeviceView> hit( hits, sa_idx );
229  const uint32 g_pos = hit.loc; // fetch the global reference coordinate
230 
231  // find the sequence
232  const uint32 seq_begin = upper_bound_index(
233  g_pos,
234  reference_index,
235  reference_count+1u ) - 1u;
236 
237  // find the sequence
238  const uint32 seq_end = upper_bound_index(
239  g_pos + params.seed_len,
240  reference_index,
241  reference_count+1u ) - 1u;
242 
243  if (seq_begin != seq_end)
244  flags[ thread_id ] = 0;
245 }
246 
250 
251 } // namespace detail
252 
253 //
254 // Locate the next SA row in the queue.
255 // Since the input loc_queue might have been sorted to gather locality, the
256 // corresponding entry in seed_queue is now specified by an index (idx_queue).
257 //
258 template <typename BatchType, typename FMType, typename rFMType>
259 void locate(
260  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
261  const uint32 in_count,
262  const uint32* idx_queue,
263  HitQueuesDeviceView hits,
264  const ParamsPOD params)
265 {
266  const int blocks = (in_count + BLOCKDIM-1) / BLOCKDIM;
267 
268  detail::locate_kernel<<<blocks, BLOCKDIM>>>(
269  read_batch, fmi, rfmi,
270  in_count,
271  idx_queue,
272  hits,
273  params );
274 }
275 
276 //
277 // Locate the next SA row in the queue.
278 // Since the input loc_queue might have been sorted to gather locality, the
279 // corresponding entry in seed_queue is now specified by an index (idx_queue).
280 //
281 template <typename BatchType, typename FMType, typename rFMType>
283  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
284  const uint32 in_count,
285  const uint32* idx_queue,
286  HitQueuesDeviceView hits,
287  const ParamsPOD params)
288 {
289  const int blocks = (in_count + BLOCKDIM-1) / BLOCKDIM;
290 
291  detail::locate_init_kernel<<<blocks, BLOCKDIM>>>(
292  read_batch, fmi, rfmi,
293  in_count,
294  idx_queue,
295  hits,
296  params );
297 }
298 
299 //
300 // Locate the next SA row in the queue.
301 // Since the input loc_queue might have been sorted to gather locality, the
302 // corresponding entry in seed_queue is now specified by an index (idx_queue).
303 //
304 template <typename ScoringScheme>
307  const ParamsPOD params)
308 {
309  locate_init(
310  pipeline.reads,
311  pipeline.fmi,
312  pipeline.rfmi,
313  pipeline.hits_queue_size,
314  pipeline.idx_queue,
315  pipeline.scoring_queues.hits,
316  params );
317 }
318 
319 //
320 // Locate the next SA row in the queue.
321 // Since the input loc_queue might have been sorted to gather locality, the
322 // corresponding entry in seed_queue is now specified by an index (idx_queue).
323 //
324 template <typename BatchType, typename FMType, typename rFMType>
326  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
327  const uint32 in_count,
328  const uint32* idx_queue,
329  HitQueuesDeviceView hits,
330  const ParamsPOD params)
331 {
332  const int blocks = (in_count + BLOCKDIM-1) / BLOCKDIM;
333 
334  detail::locate_lookup_kernel<<<blocks, BLOCKDIM>>>(
335  read_batch, fmi, rfmi,
336  in_count,
337  idx_queue,
338  hits,
339  params );
340 }
341 
342 //
343 // Locate the next SA row in the queue.
344 // Since the input loc_queue might have been sorted to gather locality, the
345 // corresponding entry in seed_queue is now specified by an index (idx_queue).
346 //
347 template <typename ScoringScheme>
350  const ParamsPOD params)
351 {
352  const int blocks = (pipeline.hits_queue_size + BLOCKDIM-1) / BLOCKDIM;
353 
354  detail::locate_lookup_kernel<<<blocks, BLOCKDIM>>>(
355  pipeline.reads,
356  pipeline.fmi,
357  pipeline.rfmi,
358  pipeline.hits_queue_size,
359  pipeline.idx_queue,
360  pipeline.scoring_queues.hits,
361  params );
362 }
363 
364 //
365 // mark seeds straddling the reference boundaries.
366 //
367 template <typename index_iterator, typename flags_iterator>
369  const uint32 in_count,
370  const uint32* idx_queue,
371  const uint32 reference_count,
372  const index_iterator reference_index,
373  HitQueuesDeviceView hits,
374  flags_iterator flags,
375  const ParamsPOD params)
376 {
377  const int blocks = (in_count + BLOCKDIM-1) / BLOCKDIM;
378 
379  detail::mark_straddling_kernel<<<blocks, BLOCKDIM>>>(
380  in_count,
381  idx_queue,
382  reference_count,
383  reference_index,
384  hits,
385  flags,
386  params );
387 }
388 
389 } // namespace cuda
390 } // namespace bowtie2
391 } // namespace nvbio