NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mapping_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 <nvBowtie/bowtie2/cuda/defs.h>
31 #include <nvBowtie/bowtie2/cuda/utils.h>
32 #include <nvBowtie/bowtie2/cuda/params.h>
33 #include <nvBowtie/bowtie2/cuda/seed_hit.h>
34 #include <nvBowtie/bowtie2/cuda/seed_hit_deque_array.h>
35 #include <nvbio/io/utils.h>
36 #include <nvbio/basic/cuda/arch.h>
37 #include <nvbio/basic/cuda/pingpong_queues.h>
38 #include <nvbio/basic/cached_iterator.h>
39 #include <nvbio/basic/packedstream.h>
40 #include <nvbio/basic/priority_deque.h>
41 #include <nvbio/basic/strided_iterator.h>
42 #include <nvbio/basic/transform_iterator.h>
43 #include <nvbio/basic/index_transform_iterator.h>
44 #include <nvbio/basic/algorithms.h>
45 
46 namespace nvbio {
47 namespace bowtie2 {
48 namespace cuda {
49 
50 namespace detail {
51 
54 
57 
60 
61 // Checks if the seed's Ns are compatible with the search, i.e.
62 // if there are 0 Ns in the exactly-matched region, and at most 1 N
63 // in the inexact matching region.
64 template <typename ReadStream>
65 NVBIO_DEVICE NVBIO_FORCEINLINE
66 bool check_N(const ReadStream& seed, const uint32 exact_len, const uint32 seed_len)
67 {
68  for (uint32 i = 0; i < exact_len; ++i)
69  if (seed[i] >= 4)
70  return false;
71 
72  uint32 N_count = 0;
73  for (uint32 i = exact_len; i < seed_len; ++i)
74  {
75  if (seed[i] >= 4 && ++N_count > 1)
76  return false;
77  }
78  return true;
79 }
80 
81 // This function is a lot like match_reverse,
82 // except that it accepts lower and upper bounds on the sequence stream
83 template< typename FMType, typename StreamType > NVBIO_DEVICE NVBIO_FORCEINLINE
84 uint2 match_range(uint2 range, const FMType index, StreamType query,
85  uint32 begin, uint32 end)
86 {
87  for (uint32 i = begin; i < end && range.x <= range.y; ++i)
88  {
89  const uint8 c = query[i];
90  if (c > 3) return make_uint2(1, 0);
91 
92  const uint2 bwt_cnts = rank(index, make_uint2(range.x-1, range.y), c);
93  range.x = index.L2(c) + bwt_cnts.x + 1;
94  range.y = index.L2(c) + bwt_cnts.y;
95  }
96  return range;
97 }
98 
99 // Store an array of hits into a SeedHitDequeArray
100 //
101 NVBIO_DEVICE NVBIO_FORCEINLINE
102 void store_deque(
103  SeedHitDequeArrayDeviceView hit_deques,
104  const uint32 read_id,
105  const uint32 n_hits,
106  const SeedHit* hitstorage)
107 {
108  SeedHit* hit_data = hit_deques.alloc_deque( read_id, n_hits );
109  if (hit_data != NULL)
110  {
111  hit_deques.resize_deque( read_id, n_hits );
112  for (uint32 i = 0; i < n_hits; ++i)
113  hit_data[i] = hitstorage[i];
114  }
115 }
116 
117 enum { CHECK_EXACT=1u, IGNORE_EXACT=0u };
118 enum MappingAlgorithm
119 {
120  EXACT_MAPPING = 0u,
121  APPROX_MAPPING = 1u,
122  CASE_PRUNING_MAPPING = 2u,
123 };
124 
128 template< bool find_exact,
129  typename Stream,
130  typename FMType,
131  typename HitType> NVBIO_DEVICE NVBIO_FORCEINLINE
132 void map(const Stream query,
133  /*const*/ uint32 len1,
134  const uint32 len2,
135  const FMType index,
136  const SeedHit::Flags hit_flags,
137  HitType& hitheap, uint32 max_hits,
138  uint32& range_sum,
139  uint32& range_count)
140 {
141  // TODO:
142  // if we know there is one (and only one) N in the second part of the seed,
143  // we could avoid doing the mismatch search at all other positions.
144  // However, due to the control divergence it would generate, it might be
145  // better to filter these seeds away and leave them for a second kernel
146  // call.
147  //
148 
149  // if there is an N, we can do exact matching till the N
150  uint32 N_pos = 0;
151  uint32 N_cnt = 0;
152  for (uint32 i = 0; i < len2; ++i)
153  {
154  const uint8 c = query[i];
155  if (c > 3)
156  {
157  // if the N is in our exact-matching region, or it is the second one, we cannot align this
158  if (i < len1 || N_cnt) return;
159 
160  N_pos = i;
161  N_cnt++;
162  }
163  }
164  // do exact-matching till the N, if there is one
165  len1 = N_cnt ? N_pos : len1;
166 
167  // Exact Match first half of read
168  uint2 base_range = make_uint2(0, index.length());
169  base_range = match_range(base_range, index, query, 0, len1);
170 
171  // Allow an error in the next batch
172  for(uint32 i = len1; i < len2 && base_range.x <= base_range.y; ++i)
173  {
174  const uint8 c = query[i];
175 
176  // Get ranges for all chars at pos i
177  uint4 occ_lo, occ_hi;
178  rank4(index, make_uint2(base_range.x-1, base_range.y), &occ_lo, &occ_hi);
179 
180  // Try each character
181  for (uint8 sub = 0; sub < 4; ++sub)
182  {
183  if (sub != c && comp(occ_hi, sub) > comp(occ_lo, sub))
184  {
185  // Manually insert sub
186  uint2 range = make_uint2(index.L2(sub) + comp(occ_lo, sub) + 1,
187  index.L2(sub) + comp(occ_hi, sub) );
188 
189  // Extend match to end of read
190  range = match_range(range, index, query, i+1, len2);
191 
192  // Store Result
193  if (range.x <= range.y) {
194  if (hitheap.size() == max_hits)
195  hitheap.pop_bottom();
196  hitheap.push(SeedHit(hit_flags, inclusive_to_exclusive(range)));
197 
198  range_sum += range.y - range.x + 1u;
199  range_count++;
200  }
201  }
202  }
203 
204  // Extend base_range
205  if (c < 4) {
206  base_range.x = index.L2(c) + comp(occ_lo, c) + 1;
207  base_range.y = index.L2(c) + comp(occ_hi, c);
208  } else { base_range = make_uint2(1, 0); break; }
209  }
210 
211  if (find_exact && base_range.x <= base_range.y)
212  {
213  if (hitheap.size() == max_hits)
214  hitheap.pop_bottom();
215  hitheap.push(SeedHit(hit_flags, inclusive_to_exclusive(base_range)));
216 
217  range_sum += base_range.y - base_range.x + 1u;
218  range_count++;
219  }
220 }
221 
223 template <MappingAlgorithm ALG>
224 struct seed_mapper {};
225 
229 template <>
230 struct seed_mapper<EXACT_MAPPING>
231 {
232  template<typename SeedIterator, typename FMType, typename rFMType, typename HitType>
233  NVBIO_FORCEINLINE NVBIO_DEVICE
234  static void enact(
235  const FMType fmi, const rFMType rfmi,
236  const SeedIterator seed,
237  const uint2 read_range,
238  const uint32 pos,
239  const uint32 seed_len,
240  HitType& hitheap,
241  uint32& range_sum,
242  uint32& range_count,
243  const ParamsPOD& params,
244  const bool fw,
245  const bool rc)
246  {
247  const OffsetXform <typename SeedIterator::index_type> forward_offset(0);
248  const ReverseXform<typename SeedIterator::index_type> reverse_offset(seed_len);
249 
250  typedef index_transform_iterator< SeedIterator, OffsetXform <typename SeedIterator::index_type> > fSeedReader;
251  typedef index_transform_iterator< SeedIterator, ReverseXform<typename SeedIterator::index_type> > rSeedReader;
252 
253  SeedHit::Flags flags;
254  uint2 range;
255 
256  // forward scan, forward index=0
257  const fSeedReader f_reader(seed, forward_offset);
258  if (util::count_occurrences( f_reader, seed_len, 4u, 1u ))
259  return;
260 
261  if (fw)
262  {
263  range = make_uint2(0, fmi.length());
264  range = match_range(range, fmi, f_reader, 0, seed_len);
265  if (range.x <= range.y)
266  {
267  flags = SeedHit::build_flags(STANDARD, FORWARD, read_range.y-pos-seed_len);
268  if (hitheap.size() == params.max_hits)
269  hitheap.pop_bottom();
270  hitheap.push(SeedHit(flags, inclusive_to_exclusive(range)));
271 
272  range_sum += range.y - range.x + 1u;
273  range_count++;
274  }
275  }
276 
277  if (rc)
278  {
279  #if USE_REVERSE_INDEX
280  // Complement seed=1, forward scan, reverse index=1
281  const transform_iterator<fSeedReader, complement_functor<4> > cf_reader(f_reader, complement_functor<4>());
282 
283  range = make_uint2(0, rfmi.length());
284  range = match_range(range, rfmi, cf_reader, 0, seed_len);
285  if (range.x <= range.y)
286  {
287  flags = SeedHit::build_flags(COMPLEMENT, REVERSE, pos-read_range.x+seed_len-1);
288  if (hitheap.size() == params.max_hits)
289  hitheap.pop_bottom();
290  hitheap.push(SeedHit(flags, inclusive_to_exclusive(range)));
291  }
292  #else
293  // Complement seed=1, reverse scan, forward index=0
294  const rSeedReader r_reader(seed, reverse_offset);
295  const transform_iterator<rSeedReader, complement_functor<4> > cr_reader(r_reader, complement_functor<4>());
296 
297  range = make_uint2(0, fmi.length());
298  range = match_range(range, fmi, cr_reader, 0, seed_len);
299  if (range.x <= range.y)
300  {
301  flags = SeedHit::build_flags(COMPLEMENT, FORWARD, pos-read_range.x);
302  if (hitheap.size() == params.max_hits)
303  hitheap.pop_bottom();
304  hitheap.push(SeedHit(flags, inclusive_to_exclusive(range)));
305 
306  range_sum += range.y - range.x + 1u;
307  range_count++;
308  }
309  #endif
310  }
311  }
312 };
313 
318 template <>
319 struct seed_mapper<APPROX_MAPPING>
320 {
321  template<typename SeedIterator, typename FMType, typename rFMType, typename HitType>
322  NVBIO_FORCEINLINE NVBIO_DEVICE
323  static void enact(
324  const FMType fmi, const rFMType rfmi,
325  const SeedIterator seed,
326  const uint2 read_range,
327  const uint32 pos,
328  const uint32 seed_len,
329  HitType& hitheap,
330  uint32& range_sum,
331  uint32& range_count,
332  const ParamsPOD& params,
333  const bool fw,
334  const bool rc)
335  {
336  const OffsetXform <typename SeedIterator::index_type> forward_offset(0);
337  const ReverseXform<typename SeedIterator::index_type> reverse_offset(seed_len);
338 
339  typedef index_transform_iterator< SeedIterator, OffsetXform <typename SeedIterator::index_type> > fSeedReader;
340  typedef index_transform_iterator< SeedIterator, ReverseXform<typename SeedIterator::index_type> > rSeedReader;
341 
342  SeedHit::Flags flags;
343 
344  // Standard seed=0, forward scan, forward index=0
345  const fSeedReader f_reader(seed, forward_offset);
346  if (util::count_occurrences( f_reader, seed_len, 4u, 2u ))
347  return;
348 
349  if (fw)
350  {
351  flags = SeedHit::build_flags(STANDARD, FORWARD, read_range.y-pos-seed_len);
352  map<CHECK_EXACT> (f_reader, params.subseed_len, seed_len, fmi, flags, hitheap, params.max_hits, range_sum, range_count);
353  }
354 
355  // Complement seed=1, reverse scan, forward index=0
356  if (rc)
357  {
358  const rSeedReader r_reader(seed, reverse_offset);
359  const transform_iterator<rSeedReader, complement_functor<4> > cr_reader(r_reader, complement_functor<4>());
360  flags = SeedHit::build_flags(COMPLEMENT, FORWARD, pos-read_range.x);
361  map<IGNORE_EXACT>(cr_reader, params.subseed_len, seed_len, fmi, flags, hitheap, params.max_hits, range_sum, range_count);
362  }
363  }
364 };
365 
372 template <>
373 struct seed_mapper<CASE_PRUNING_MAPPING>
374 {
375  template<typename SeedIterator, typename FMType, typename rFMType, typename HitType>
376  NVBIO_FORCEINLINE NVBIO_DEVICE
377  static void enact(
378  const FMType fmi, const rFMType rfmi,
379  const SeedIterator seed,
380  const uint2 read_range,
381  const uint32 pos,
382  const uint32 seed_len,
383  HitType& hitheap,
384  uint32& range_sum,
385  uint32& range_count,
386  const ParamsPOD& params,
387  const bool fw,
388  const bool rc)
389  {
390  const OffsetXform <typename SeedIterator::index_type> forward_offset(0);
391  const ReverseXform<typename SeedIterator::index_type> reverse_offset(seed_len);
392 
393  typedef index_transform_iterator< SeedIterator, OffsetXform <typename SeedIterator::index_type> > fSeedReader;
394  typedef index_transform_iterator< SeedIterator, ReverseXform<typename SeedIterator::index_type> > rSeedReader;
395 
396  SeedHit::Flags flags;
397 
398  // Standard seed=0, forward scan, forward index=0
399  const fSeedReader f_reader(seed, forward_offset);
400  if (fw)
401  {
402  flags = SeedHit::build_flags(STANDARD, FORWARD, read_range.y-pos-seed_len);
403  map<CHECK_EXACT> (f_reader, (seed_len )/2, seed_len, fmi, flags, hitheap, params.max_hits, range_sum, range_count);
404  }
405  // Standard seed=0, reverse scan, reverse index=1
406  const rSeedReader r_reader(seed, reverse_offset);
407  if (fw)
408  {
409  flags = SeedHit::build_flags(STANDARD, REVERSE, read_range.y-pos-1);
410  map<IGNORE_EXACT>(r_reader, (seed_len+1)/2, seed_len, rfmi, flags, hitheap, params.max_hits, range_sum, range_count);
411  }
412 
413  // Complement seed=1, forward scan, reverse index=1
414  if (rc)
415  {
416  const transform_iterator<fSeedReader, complement_functor<4> > cf_reader(f_reader, complement_functor<4>());
417  flags = SeedHit::build_flags(COMPLEMENT, REVERSE, pos-read_range.x+seed_len-1);
418  map<CHECK_EXACT> (cf_reader, (seed_len )/2, seed_len, rfmi, flags, hitheap, params.max_hits, range_sum, range_count);
419  }
420 
421  // Complement seed=1, reverse scan, forward index=0
422  if (rc)
423  {
424  const transform_iterator<rSeedReader, complement_functor<4> > cr_reader(r_reader, complement_functor<4>());
425  flags = SeedHit::build_flags(COMPLEMENT, FORWARD, pos-read_range.x);
426  map<IGNORE_EXACT>(cr_reader, (seed_len+1)/2, seed_len, fmi, flags, hitheap, params.max_hits, range_sum, range_count);
427  }
428  }
429 };
430 
434 template<typename BatchType, typename FMType, typename rFMType> __global__
435 void map_whole_read_kernel(
436  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
437  const nvbio::cuda::PingPongQueuesView<uint32> queues,
438  uint8* reseed,
439  SeedHitDequeArrayDeviceView hits,
440  const ParamsPOD params,
441  const bool fw,
442  const bool rc)
443 {
444  typedef PackedStringLoader<
445  typename BatchType::sequence_storage_iterator,
446  BatchType::SEQUENCE_BITS,
447  BatchType::SEQUENCE_BIG_ENDIAN,
448  uncached_tag > packed_loader_type;
449  typedef typename packed_loader_type::iterator seed_iterator;
450 
451  // instantiate a packed loader to cache the seed in local memory
452  packed_loader_type packer_loader;
453 
454  // instantiate storage for a local memory hits queue
455  SeedHit local_hits[512];
456 
457  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
458  const uint32 grid_size = BLOCKDIM * gridDim.x;
459 
460  for (uint32 id = thread_id; id < queues.in_size; id += grid_size)
461  {
462  const uint32 read_id = queues.in_queue[ id ];
463  NVBIO_CUDA_ASSERT( read_id < read_batch.size() );
464  const uint2 read_range = read_batch.get_range( read_id );
465 
466  // filter away short strings
467  if (read_range.y - read_range.x < params.min_read_len)
468  {
469  hits.resize_deque( read_id, 0 );
470  return;
471  }
472 
473  typedef SeedHitDequeArrayDeviceView::hit_vector_type hit_storage_type;
474  typedef SeedHitDequeArrayDeviceView::hit_deque_type hit_deque_type;
475  hit_deque_type hitheap( hit_storage_type( 0, local_hits ), true );
476 
477  uint32 range_sum = 0;
478  uint32 range_count = 0;
479 
480  // load the whole read
481  const seed_iterator read = packer_loader.load( read_batch.sequence_stream() + read_range.x, read_range.y - read_range.x );
482 
483  // map the read
484  seed_mapper<EXACT_MAPPING>::enact(
485  fmi, rfmi,
486  read,
487  read_range,
488  read_range.x,
489  read_range.y - read_range.x,
490  hitheap,
491  range_sum,
492  range_count,
493  params,
494  fw,
495  rc );
496 
497  // save the hits
498  store_deque( hits, read_id, hitheap.size(), local_hits );
499 
500  // enqueue for re-seeding if there are no hits
501  if (reseed)
502  reseed[ id ] = (range_count == 0);
503  }
504 }
505 
511 template<MappingAlgorithm ALGO, typename BatchType, typename FMType, typename rFMType> __global__
512 void map_queues_kernel(
513  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
514  const uint32 retry,
515  const nvbio::cuda::PingPongQueuesView<uint32> queues,
516  uint8* reseed,
517  SeedHitDequeArrayDeviceView hits,
518  const ParamsPOD params,
519  const bool fw,
520  const bool rc)
521 {
522  typedef PackedStringLoader<
523  typename BatchType::sequence_storage_iterator,
524  BatchType::SEQUENCE_BITS,
525  BatchType::SEQUENCE_BIG_ENDIAN,
526  uncached_tag > packed_loader_type;
527  //lmem_cache_tag<128> > packed_loader_type; // TODO: this produces a crashing kernel, for no good reason
528  typedef typename packed_loader_type::iterator seed_iterator;
529 
530  // instantiate a packed loader to cache the seed in local memory
531  packed_loader_type packer_loader;
532 
533  // instantiate storage for a local memory hits queue
534  SeedHit local_hits[512];
535 
536  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
537  const uint32 grid_size = BLOCKDIM * gridDim.x;
538 
539  for (uint32 id = thread_id; id < queues.in_size; id += grid_size)
540  {
541  const uint32 read_id = queues.in_queue[ id ];
542  NVBIO_CUDA_ASSERT( read_id < read_batch.size() );
543  const uint2 read_range = read_batch.get_range( read_id );
544  const uint32 read_len = read_range.y - read_range.x;
545 
546  // filter away short strings
547  if (read_len < params.min_read_len)
548  {
549  hits.resize_deque( read_id, 0 );
550  continue;
551  }
552 
553  const uint32 seed_len = nvbio::min( params.seed_len, read_len );
554  const uint32 seed_freq = (uint32)params.seed_freq( read_len );
555  const uint32 retry_stride = seed_freq/(params.max_reseed+1);
556 
557  typedef SeedHitDequeArrayDeviceView::hit_vector_type hit_storage_type;
558  typedef SeedHitDequeArrayDeviceView::hit_deque_type hit_deque_type;
559  hit_deque_type hitheap( hit_storage_type( 0, local_hits ), true );
560 
561  uint32 range_sum = 0;
562  uint32 range_count = 0;
563 
564  // loop over seeds
565  for (uint32 pos = read_range.x + retry * retry_stride; pos + seed_len <= read_range.y; pos += seed_freq)
566  {
567  const seed_iterator seed = packer_loader.load( read_batch.sequence_stream() + pos, seed_len );
568 
569  seed_mapper<ALGO>::enact(
570  fmi, rfmi,
571  seed,
572  read_range,
573  pos,
574  seed_len,
575  hitheap,
576  range_sum,
577  range_count,
578  params,
579  fw,
580  rc );
581  }
582 
583  // save the hits
584  store_deque( hits, read_id, hitheap.size(), local_hits );
585 
586  // enqueue for re-seeding if there are too many hits
587  if (reseed)
588  reseed[ id ] = (range_count == 0 || range_sum >= params.rep_seeds * range_count);
589 
590  NVBIO_CUDA_DEBUG_PRINT_IF( params.debug.read_id == read_id, "map: ranges[%u], rows[%u]\n", read_id, range_count, range_sum );
591  }
592 }
593 
599 template<MappingAlgorithm ALGO, typename BatchType, typename FMType, typename rFMType> __global__
600 void map_kernel(
601  const BatchType read_batch, const FMType fmi, const rFMType rfmi,
602  SeedHitDequeArrayDeviceView hits,
603  const uint2 seed_range,
604  const ParamsPOD params,
605  const bool fw,
606  const bool rc)
607 {
608  typedef PackedStringLoader<
609  typename BatchType::sequence_storage_iterator,
610  BatchType::SEQUENCE_BITS,
611  BatchType::SEQUENCE_BIG_ENDIAN,
612  uncached_tag > packed_loader_type;
613  //lmem_cache_tag<128> > packed_loader_type; // TODO: this produces a crashing kernel, for no good reason
614  typedef typename packed_loader_type::iterator seed_iterator;
615 
616  // instantiate a packed loader to cache the seed in local memory
617  packed_loader_type packer_loader;
618 
619  // instantiate storage for a local memory hits queue
620  SeedHit local_hits[512];
621 
622  const uint32 thread_id = threadIdx.x + BLOCKDIM*blockIdx.x;
623  const uint32 grid_size = BLOCKDIM * gridDim.x;
624 
625  for (uint32 read_id = thread_id; read_id < read_batch.size(); read_id += grid_size)
626  {
627  const uint2 read_range = read_batch.get_range( read_id );
628  const uint32 read_len = read_range.y - read_range.x;
629 
630  // filter away short strings
631  if (read_len < params.min_read_len)
632  {
633  hits.resize_deque( read_id, 0 );
634  return;
635  }
636 
637  const uint32 seed_len = nvbio::min( params.seed_len, read_len );
638  const uint32 seed_freq = (uint32)params.seed_freq( read_len );
639  const uint32 retry_stride = seed_freq/(params.max_reseed+1);
640 
641  typedef SeedHitDequeArrayDeviceView::hit_vector_type hit_storage_type;
642  typedef SeedHitDequeArrayDeviceView::hit_deque_type hit_deque_type;
643  hit_deque_type hitheap( hit_storage_type( 0, local_hits ), true );
644 
645  bool active = true;
646 
647  //for (uint32 retry = 0; retry <= params.max_reseed; ++retry) // TODO: this makes the kernel crash even when
648  // max_reseed is 0, for apparently no good reason...
649  for (uint32 retry = 0; retry <= 0; ++retry)
650  {
651  if (active)
652  {
653  // restore the result heap
654  hitheap.clear();
655 
656  uint32 range_sum = 0;
657  uint32 range_count = 0;
658 
659  for (uint32 i = seed_range.x; i < seed_range.y; ++i)
660  {
661  const uint32 pos = read_range.x + retry * retry_stride + i * seed_freq;
662  if (pos + seed_len <= read_range.y)
663  {
664  const seed_iterator seed = packer_loader.load( read_batch.sequence_stream() + pos, seed_len );
665 
666  seed_mapper<ALGO>::enact(
667  fmi, rfmi,
668  seed,
669  read_range,
670  pos,
671  seed_len,
672  hitheap,
673  range_sum,
674  range_count,
675  params,
676  fw,
677  rc );
678  }
679  }
680 
681  // check whether we could stop here or if we need re-seeding
682  if (retry == params.max_reseed || range_count == 0 && range_sum < params.rep_seeds * range_count)
683  active = false;
684  }
685  }
686 
687  // save the hits
688  store_deque( hits, read_id, hitheap.size(), local_hits );
689  }
690 }
691 
695 
696 } // namespace detail
697 
698 //
699 // call the map-sub kernel
700 //
701 template <typename BatchType, typename FMType, typename rFMType>
702 void map_case_pruning_t(
703  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
704  const uint32 retry,
705  const nvbio::cuda::PingPongQueuesView<uint32> queues,
706  uint8* reseed,
707  SeedHitDequeArrayDeviceView hits,
708  const ParamsPOD params,
709  const bool fw,
710  const bool rc)
711 {
712  const int blocks = nvbio::cuda::max_active_blocks( detail::map_queues_kernel<detail::CASE_PRUNING_MAPPING,BatchType,FMType,rFMType>, BLOCKDIM, 0 );
713  detail::map_queues_kernel<detail::CASE_PRUNING_MAPPING> <<<blocks, BLOCKDIM>>>(
714  read_batch, fmi, rfmi, retry, queues, reseed, hits, params, fw, rc );
715 }
716 
717 //
718 // call the map-hybrid kernel
719 //
720 template <typename BatchType, typename FMType, typename rFMType>
721 void map_approx_t(
722  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
723  const uint32 retry,
724  const nvbio::cuda::PingPongQueuesView<uint32> queues,
725  uint8* reseed,
726  SeedHitDequeArrayDeviceView hits,
727  const ParamsPOD params,
728  const bool fw,
729  const bool rc)
730 {
731  const int blocks = nvbio::cuda::max_active_blocks( detail::map_queues_kernel<detail::APPROX_MAPPING,BatchType,FMType,rFMType>, BLOCKDIM, 0 );
732  detail::map_queues_kernel<detail::APPROX_MAPPING> <<<blocks, BLOCKDIM>>>(
733  read_batch, fmi, rfmi, retry, queues, reseed, hits, params, fw, rc );
734 }
735 
736 //
737 // call the map-hybrid kernel
738 //
739 template <typename BatchType, typename FMType, typename rFMType>
740 void map_approx_t(
741  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
742  SeedHitDequeArrayDeviceView hits,
743  const uint2 seed_range,
744  const ParamsPOD params,
745  const bool fw,
746  const bool rc)
747 {
748  const int blocks = nvbio::cuda::max_active_blocks( detail::map_kernel<detail::APPROX_MAPPING,BatchType,FMType,rFMType>, BLOCKDIM, 0 );
749  detail::map_kernel<detail::APPROX_MAPPING> <<<blocks, BLOCKDIM>>>(
750  read_batch, fmi, rfmi, hits, seed_range, params, fw, rc );
751 }
752 
753 //
754 // call the map-exact kernel
755 //
756 template <typename BatchType, typename FMType, typename rFMType>
757 void map_whole_read_t(
758  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
759  const nvbio::cuda::PingPongQueuesView<uint32> queues,
760  uint8* reseed,
761  SeedHitDequeArrayDeviceView hits,
762  const ParamsPOD params,
763  const bool fw,
764  const bool rc)
765 {
766  const uint32 blocks = nvbio::cuda::max_active_blocks( detail::map_whole_read_kernel<BatchType,FMType,rFMType>, BLOCKDIM, 0 );
767  detail::map_whole_read_kernel<<<blocks, BLOCKDIM>>> (
768  read_batch, fmi, rfmi, queues, reseed, hits, params, fw, rc );
769 }
770 
771 //
772 // call the map-exact kernel
773 //
774 template <typename BatchType, typename FMType, typename rFMType>
775 void map_exact_t(
776  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
777  const uint32 retry,
778  const nvbio::cuda::PingPongQueuesView<uint32> queues,
779  uint8* reseed,
780  SeedHitDequeArrayDeviceView hits,
781  const ParamsPOD params,
782  const bool fw,
783  const bool rc)
784 {
785  const uint32 blocks = nvbio::cuda::max_active_blocks( detail::map_queues_kernel<detail::EXACT_MAPPING,BatchType,FMType,rFMType>, BLOCKDIM, 0 );
786  detail::map_queues_kernel<detail::EXACT_MAPPING> <<<blocks, BLOCKDIM>>>(
787  read_batch, fmi, rfmi, retry, queues, reseed, hits, params, fw, rc );
788 }
789 
790 //
791 // call the map-exact kernel
792 //
793 template <typename BatchType, typename FMType, typename rFMType>
794 void map_exact_t(
795  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
796  SeedHitDequeArrayDeviceView hits,
797  const uint2 seed_range,
798  const ParamsPOD params,
799  const bool fw,
800  const bool rc)
801 {
802  const uint32 blocks = nvbio::cuda::max_active_blocks( detail::map_kernel<detail::EXACT_MAPPING,BatchType,FMType,rFMType>, BLOCKDIM, 0 );
803  detail::map_kernel<detail::EXACT_MAPPING> <<<blocks, BLOCKDIM>>>(
804  read_batch, fmi, rfmi, hits, seed_range, params, fw, rc );
805 }
806 
807 //
808 // call the appropriate mapping kernel
809 //
810 template <typename BatchType, typename FMType, typename rFMType>
811 void map_t(
812  const BatchType& read_batch, const FMType fmi, const rFMType rfmi,
813  const uint32 retry,
814  const nvbio::cuda::PingPongQueuesView<uint32> queues,
815  uint8* reseed,
816  SeedHitDequeArrayDeviceView hits,
817  const ParamsPOD params,
818  const bool fw,
819  const bool rc)
820 {
821  // check whether we allow substitutions
822  if (params.allow_sub)
823  {
824  // check whether we have to perform exact matching on a subseed
825  if (params.subseed_len == 0)
826  {
827  // call the hybrid fuzzy mapping kernel
828  map_case_pruning_t(
829  read_batch, fmi, rfmi, retry, queues, reseed, hits, params, fw, rc );
830  }
831  else
832  {
833  // call the hybrid exact-fuzzy mapping kernel
834  map_approx_t(
835  read_batch, fmi, rfmi, retry, queues, reseed, hits, params, fw, rc );
836  }
837  }
838  else
839  {
840  // call the hybrid exact mapping kernel
841  map_exact_t(
842  read_batch, fmi, rfmi, retry, queues, reseed, hits, params, fw, rc );
843  }
844 }
845 
846 } // namespace cuda
847 } // namespace bowtie2
848 } // namespace nvbio
Generated on Wed Feb 25 2015 08:33:05 for NVBIO by   doxygen 1.8.4