NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hamming_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 
31 #include <nvbio/basic/packedstream.h>
32 #include <nvbio/alignment/sink.h>
33 #include <nvbio/alignment/utils.h>
34 #include <nvbio/alignment/alignment_base_inl.h>
35 #include <nvbio/basic/iterator.h>
36 
37 
38 namespace nvbio {
39 namespace aln {
40 
41 namespace priv
42 {
43 
46 
55 template <uint32 BAND_LEN, AlignmentType TYPE, typename algorithm_tag = PatternBlockingTag>
56 struct HammingScoringContext
57 {
58 
66  template <typename column_type, typename scoring_type, typename score_type>
67  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
68  void init(
69  const uint32 j,
70  const uint32 N,
71  column_type& column,
72  const scoring_type& scoring,
73  const score_type zero)
74  {
75  if (j == 0) // ensure this context can be used for multi-pass scoring
76  {
77  for (uint32 i = 0; i < N; ++i)
78  column[i] = zero;
79  }
80  }
81 
87  template <typename column_type>
88  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
89  void previous_column(
90  const uint32 j,
91  const uint32 N,
92  const column_type column) {}
93 
99  template <typename column_type>
100  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
101  void last_column(
102  const uint32 j,
103  const uint32 M,
104  const uint32 N,
105  const column_type column) {}
106 
115  template <typename score_type>
116  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
117  void new_cell(
118  const uint32 i,
119  const uint32 N,
120  const uint32 j,
121  const uint32 M,
122  const score_type score,
123  const DirectionVector dir) {}
124 };
125 
130 template <uint32 BAND_LEN, AlignmentType TYPE, typename algorithm_tag, typename checkpoint_type = PatternBlockingTag>
131 struct HammingCheckpointedScoringContext
132 {
133  // constructor
134  //
135  // \param checkpoints input checkpoints array
136  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
137  HammingCheckpointedScoringContext(checkpoint_type checkpoint) :
138  m_checkpoint( checkpoint ) {}
139 
147  template <typename column_type, typename scoring_type, typename score_type>
148  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
149  void init(
150  const uint32 j,
151  const uint32 N,
152  column_type& column,
153  const scoring_type& scoring,
154  const score_type zero)
155  {
156  if (j == 0)
157  {
158  for (uint32 i = 0; i < N; ++i)
159  column[i] = zero;
160  }
161  else
162  {
163  for (uint32 i = 0; i < N; ++i)
164  column[i] = m_checkpoint[i];
165  }
166  }
167 
173  template <typename column_type>
174  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
175  void previous_column(
176  const uint32 j,
177  const uint32 N,
178  const column_type column) {}
179 
185  template <typename column_type>
186  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
187  void last_column(
188  const uint32 j,
189  const uint32 M,
190  const uint32 N,
191  const column_type column)
192  {
193  for (uint32 i = 0; i < N; ++i)
194  {
195  m_checkpoint[i].x = column[i].x;
196  m_checkpoint[i].y = column[i].y;
197  }
198  }
199 
208  template <typename score_type>
209  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
210  void new_cell(
211  const uint32 i,
212  const uint32 N,
213  const uint32 j,
214  const uint32 M,
215  const score_type score,
216  const DirectionVector dir) {}
217 
218  checkpoint_type m_checkpoint;
219 };
220 
224 template <uint32 BAND_LEN, AlignmentType TYPE, uint32 CHECKPOINTS, typename checkpoint_type>
225 struct HammingCheckpointContext
226 {
227  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
228  HammingCheckpointContext(checkpoint_type checkpoints) :
229  m_checkpoints( checkpoints ) {}
230 
238  template <typename column_type, typename scoring_type, typename score_type>
239  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
240  void init(
241  const uint32 j,
242  const uint32 N,
243  column_type& column,
244  const scoring_type& scoring,
245  const score_type zero)
246  {
247  for (uint32 i = 0; i < N; ++i)
248  column[i] = zero;
249  }
250 
256  template <typename column_type>
257  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
258  void previous_column(
259  const uint32 j,
260  const uint32 N,
261  const column_type column)
262  {
263  // save checkpoint
264  if ((j & (CHECKPOINTS-1)) == 0u)
265  {
266  const uint32 checkpoint_id = j / CHECKPOINTS;
267 
268  for (uint32 i = 0; i < N; ++i)
269  m_checkpoints[ checkpoint_id*N + i ] = column[i];
270  }
271  }
272 
278  template <typename column_type>
279  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
280  void last_column(
281  const uint32 j,
282  const uint32 M,
283  const uint32 N,
284  const column_type column) {}
285 
294  template <typename score_type>
295  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
296  void new_cell(
297  const uint32 i,
298  const uint32 N,
299  const uint32 j,
300  const uint32 M,
301  const score_type score,
302  const DirectionVector dir) {}
303 
304  checkpoint_type m_checkpoints;
305 };
306 
311 template <uint32 BAND_LEN, AlignmentType TYPE, uint32 CHECKPOINTS, typename checkpoint_type, typename submatrix_type>
312 struct HammingSubmatrixContext
313 {
314  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
315  HammingSubmatrixContext(
316  const checkpoint_type checkpoints,
317  const uint32 checkpoint_id,
318  const submatrix_type submatrix) :
319  m_checkpoints( checkpoints ),
320  m_checkpoint_id( checkpoint_id ),
321  m_submatrix( submatrix ) {}
322 
330  template <typename column_type, typename scoring_type, typename score_type>
331  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
332  void init(
333  const uint32 j,
334  const uint32 N,
335  column_type& column,
336  const scoring_type& scoring,
337  const score_type zero)
338  {
339  // restore the checkpoint
340  for (uint32 i = 0; i < N; ++i)
341  column[i] = m_checkpoints[ m_checkpoint_id*N + i ];
342  }
343 
349  template <typename column_type>
350  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
351  void previous_column(
352  const uint32 j,
353  const uint32 N,
354  const column_type column) {}
355 
361  template <typename column_type>
362  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
363  void last_column(
364  const uint32 j,
365  const uint32 M,
366  const uint32 N,
367  const column_type column) {}
368 
369  // do something with the newly computed cell
370  //
371  // \param i row index
372  // \param N number of rows (column size)
373  // \param j column index
374  // \param M number of columns (row size)
375  // \param score computed score
376  // \param dir direction flow
377  template <typename score_type>
378  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
379  void new_cell(
380  const uint32 i,
381  const uint32 N,
382  const uint32 j,
383  const uint32 M,
384  const score_type score,
385  const DirectionVector dir)
386  {
387  // save the direction vector
388  const uint32 offset = m_checkpoint_id * CHECKPOINTS;
389 
390  if (TYPE == LOCAL)
391  m_submatrix[ i * CHECKPOINTS + (j - offset) ] = score ? dir : SINK;
392  else
393  m_submatrix[ i * CHECKPOINTS + (j - offset) ] = dir;
394  }
395 
396  checkpoint_type m_checkpoints;
397  uint32 m_checkpoint_id;
398  submatrix_type m_submatrix;
399 };
400 
401 // -------------------------- Basic Smith-Waterman functions ---------------------------- //
402 
403 template <uint32 BAND_LEN, AlignmentType TYPE, typename algorithm_tag, typename symbol_type>
404 struct hamming_alignment_score_dispatch {};
405 
413 template <uint32 BAND_LEN, AlignmentType TYPE, typename symbol_type>
414 struct hamming_alignment_score_dispatch<BAND_LEN,TYPE,PatternBlockingTag,symbol_type>
415 {
416  template <
417  bool CHECK_M,
418  typename context_type,
419  typename query_cache,
420  typename score_type,
421  typename temp_iterator,
422  typename sink_type,
423  typename scoring_type>
424  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
425  static void update_row(
426  context_type& context,
427  const uint32 block,
428  const uint32 M,
429  const uint32 i,
430  const uint32 N,
431  const symbol_type r_i,
432  const query_cache q_cache,
433  temp_iterator temp,
434  score_type& temp_i,
435  score_type* band,
436  sink_type& sink,
437  score_type& max_score,
438  const score_type zero,
439  const scoring_type scoring)
440  {
441  // set the 0-th coefficient in the band to be equal to the (i-1)-th row of the left column (diagonal term)
442  score_type prev = temp_i;
443 
444  // set the 0-th coefficient in the band to be equal to the i-th row of the left column (left term)
445  band[0] = temp_i = temp[i];
446 
447  // update all cells
448  #pragma unroll
449  for (uint32 j = 1; j <= BAND_LEN; ++j)
450  {
451  const symbol_type q_j = q_cache[ j-1 ].x;
452  const uint8 qq_j = q_cache[ j-1 ].y;
453  const score_type S_ij = (r_i == q_j) ? scoring.match(qq_j) : scoring.mismatch( r_i, q_j, qq_j );
454  score_type hi = prev + S_ij;
455  hi = TYPE == LOCAL ? nvbio::max( hi, zero ) : hi; // clamp to zero
456  prev = band[j];
457  band[j] = hi;
458 
459  if ((CHECK_M == false) || (block + j <= M))
460  {
461  context.new_cell(
462  i, N,
463  block + j - 1, M,
464  hi,
465  SUBSTITUTION );
466  }
467  }
468 
469  // save the last entry of the band
470  temp[i] = band[ BAND_LEN ];
471 
472  max_score = nvbio::max( max_score, band[ BAND_LEN ] );
473 
474  if (TYPE == LOCAL)
475  {
476  if (CHECK_M)
477  {
478  // during local alignment we save the best score across all bands
479  for (uint32 j = 1; j <= BAND_LEN; ++j)
480  {
481  if (block + j <= M)
482  sink.report( band[j], make_uint2( i+1, block+j ) );
483  }
484  }
485  else
486  {
487  // during local alignment we save the best score across all bands
488  for (uint32 j = 1; j <= BAND_LEN; ++j)
489  sink.report( band[j], make_uint2( i+1, block+j ) );
490  }
491  }
492  else if (CHECK_M)
493  {
494  if (TYPE == SEMI_GLOBAL)
495  {
496  // during semi-global alignment we save the best score across the last column H[*][M], at each row
497  save_boundary<BAND_LEN>( block, M, band, i, sink );
498  }
499  }
500  }
501 
574  template <
575  typename context_type,
576  typename string_type,
577  typename qual_type,
578  typename ref_type,
579  typename scoring_type,
580  typename sink_type,
581  typename column_type>
582  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
583  static
584  bool run(
585  const scoring_type& scoring,
586  context_type& context,
587  string_type query,
588  qual_type quals,
589  ref_type ref,
590  const int32 min_score,
591  sink_type& sink,
592  const uint32 window_begin,
593  const uint32 window_end,
594  column_type temp)
595  {
596  const uint32 M = query.length();
597  const uint32 N = ref.length();
598 
599  typedef int32 score_type;
600  score_type band[BAND_LEN+1];
601 
602  const score_type zero = score_type(0);
603 
604  // we are going to process the matrix in stripes, where each partial row within each stripe is processed
605  // in registers: hence, we need to save the right-boundary of each stripe (a full matrix column) in some
606  // temporary storage.
607 
608  // fill row 0 with zeros, this will never be updated
609  uint2 q_cache[ BAND_LEN ];
610 
611  // initialize the first column
612  context.init( window_begin, N, temp, scoring, zero );
613 
614  const uint32 end_block = (window_end == M) ?
615  nvbio::max( BAND_LEN, BAND_LEN * ((M + BAND_LEN-1) / BAND_LEN) ) :
616  window_end + BAND_LEN;
617 
618  // loop across the short edge of the DP matrix (i.e. the columns)
619  for (uint32 block = window_begin; block + BAND_LEN < end_block; block += BAND_LEN)
620  {
621  // save the previous column
622  context.previous_column( block, N, temp );
623 
624  // initialize the first band (corresponding to the 0-th row of the DP matrix)
625  #pragma unroll
626  for (uint32 j = 0; j <= BAND_LEN; ++j)
627  band[j] = zero;
628 
629  // load a block of entries from each query
630  #pragma unroll
631  for (uint32 t = 0; t < BAND_LEN; ++t)
632  q_cache[ t ] = make_uint2( query[ block + t ], quals[ block + t ] );
633 
634  score_type max_score = Field_traits<score_type>::min();
635 
636  score_type temp_i = band[0];
637 
638  // loop across the long edge of the DP matrix (i.e. the rows)
639  for (uint32 i = 0; i < N; ++i)
640  {
641  // load the new character from the reference
642  const uint8 r_i = ref[i];
643 
644  update_row<false>(
645  context,
646  block, M,
647  i, N,
648  r_i,
649  q_cache,
650  temp,
651  temp_i,
652  band,
653  sink,
654  max_score,
655  zero,
656  scoring );
657  }
658 
659  // we are now (M - block - BAND_LEN) columns from the last one: check whether
660  // we could theoretically reach the minimum score
661  const score_type missing_cols = score_type(M - block - BAND_LEN);
662  if (max_score + missing_cols * scoring.match(255) < score_type( min_score ))
663  return false;
664  }
665 
666  if (window_end == M)
667  {
668  const uint32 block = end_block - BAND_LEN;
669 
670  // save the previous column
671  context.previous_column( block, N, temp );
672 
673  // initialize the first band (corresponding to the 0-th row of the DP matrix)
674  for (uint32 j = 0; j <= BAND_LEN; ++j)
675  band[j] = zero;
676 
677  // load a block of entries from each query
678  const uint32 block_end = nvbio::min( block + BAND_LEN, M );
679  #pragma unroll
680  for (uint32 t = 0; t < BAND_LEN; ++t)
681  {
682  if (block + t < block_end)
683  q_cache[ t ] = make_uint2( query[ block + t ], quals[ block + t ] );
684  }
685 
686  score_type max_score = Field_traits<score_type>::min();
687 
688  score_type temp_i = band[0];
689 
690  // loop across the long edge of the DP matrix (i.e. the rows)
691  for (uint32 i = 0; i < N; ++i)
692  {
693  // load the new character from the reference
694  const uint8 r_i = ref[i];
695 
696  update_row<true>(
697  context,
698  block, M,
699  i, N,
700  r_i,
701  q_cache,
702  temp,
703  temp_i,
704  band,
705  sink,
706  max_score,
707  zero,
708  scoring );
709  }
710  }
711 
712  // save the last column
713  context.last_column( window_end, M, N, temp );
714 
715  if (TYPE == GLOBAL)
716  save_Mth<BAND_LEN>( M, band, N-1, sink );
717 
718  return true;
719  }
720 
795  template <
796  uint32 MAX_REF_LEN,
797  typename context_type,
798  typename string_type,
799  typename qual_type,
800  typename ref_type,
801  typename scoring_type,
802  typename sink_type>
803  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
804  static
805  bool run(
806  const scoring_type& scoring,
807  context_type& context,
808  string_type query,
809  qual_type quals,
810  ref_type ref,
811  const int32 min_score,
812  sink_type& sink,
813  const uint32 window_begin,
814  const uint32 window_end)
815  {
816  // instantiated a local memory array
817  int16 temp[MAX_REF_LEN];
818  int16* temp_ptr = temp;
819 
820  return run(
821  scoring,
822  context,
823  query,
824  quals,
825  ref,
826  min_score,
827  sink,
828  window_begin,
829  window_end,
830  temp_ptr );
831  }
832 };
833 
841 template <uint32 BAND_LEN, AlignmentType TYPE, typename symbol_type>
842 struct hamming_alignment_score_dispatch<BAND_LEN,TYPE,TextBlockingTag,symbol_type>
843 {
844  template <
845  bool CHECK_N,
846  typename context_type,
847  typename ref_cache,
848  typename score_type,
849  typename temp_iterator,
850  typename sink_type,
851  typename scoring_type>
852  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
853  static void update_row(
854  context_type& context,
855  const uint32 block,
856  const uint32 N,
857  const uint32 i,
858  const uint32 M,
859  const symbol_type q_i,
860  const uint8 qq_i,
861  const score_type m_i,
862  const ref_cache r_cache,
863  temp_iterator temp,
864  score_type& temp_i,
865  score_type* band,
866  sink_type& sink,
867  score_type& max_score,
868  const score_type zero,
869  const scoring_type scoring)
870  {
871  // set the 0-th coefficient in the band to be equal to the (i-1)-th row of the left column (diagonal term)
872  score_type prev = temp_i;
873 
874  // set the 0-th coefficient in the band to be equal to the i-th row of the left column (left term)
875  band[0] = temp_i = temp[i];
876 
877  // update all cells
878  #pragma unroll
879  for (uint32 j = 1; j <= BAND_LEN; ++j)
880  {
881  const symbol_type r_j = r_cache[ j-1 ];
882  const score_type S_ij = (r_j == q_i) ? m_i : scoring.mismatch( r_j, q_i, qq_i );
883  score_type hi = prev + S_ij;
884  hi = TYPE == LOCAL ? nvbio::max( hi, zero ) : hi; // clamp to zero
885  prev = band[j];
886  band[j] = hi;
887 
888  if ((CHECK_N == false) || (block + j <= N))
889  {
890  context.new_cell(
891  i, M,
892  block + j - 1, N,
893  hi,
894  SUBSTITUTION );
895  }
896  }
897 
898  // save the last entry of the band
899  temp[i] = band[ BAND_LEN ];
900 
901  max_score = nvbio::max( max_score, band[ BAND_LEN ] );
902 
903  if (TYPE == LOCAL)
904  {
905  if (CHECK_N)
906  {
907  // during local alignment we save the best score across all bands
908  for (uint32 j = 1; j <= BAND_LEN; ++j)
909  {
910  if (block + j <= N)
911  sink.report( band[j], make_uint2( block+j, i+1 ) );
912  }
913  }
914  else
915  {
916  // during local alignment we save the best score across all bands
917  for (uint32 j = 1; j <= BAND_LEN; ++j)
918  sink.report( band[j], make_uint2( block+j, i+1 ) );
919  }
920  }
921  }
922 
995  template <
996  typename context_type,
997  typename string_type,
998  typename qual_type,
999  typename ref_type,
1000  typename scoring_type,
1001  typename sink_type,
1002  typename column_type>
1003  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1004  static
1005  bool run(
1006  const scoring_type& scoring,
1007  context_type& context,
1008  string_type query,
1009  qual_type quals,
1010  ref_type ref,
1011  const int32 min_score,
1012  sink_type& sink,
1013  const uint32 window_begin,
1014  const uint32 window_end,
1015  column_type temp)
1016  {
1017  const uint32 M = query.length();
1018  const uint32 N = ref.length();
1019 
1020  typedef int32 score_type;
1021  score_type band[BAND_LEN+1];
1022 
1023  const score_type zero = score_type(0);
1024 
1025  // we are going to process the matrix in stripes, where each partial row within each stripe is processed
1026  // in registers: hence, we need to save the right-boundary of each stripe (a full matrix column) in some
1027  // temporary storage.
1028 
1029  uint8 r_cache[ BAND_LEN ];
1030 
1031  // initialize the first column
1032  context.init( window_begin, M, temp, scoring, zero );
1033 
1034  const uint32 end_block = (window_end == N) ?
1035  nvbio::max( BAND_LEN, BAND_LEN * ((N + BAND_LEN-1) / BAND_LEN) ) :
1036  window_end + BAND_LEN;
1037 
1038  // loop across the short edge of the DP matrix (i.e. the columns)
1039  for (uint32 block = window_begin; block + BAND_LEN < end_block; block += BAND_LEN)
1040  {
1041  // save the previous column
1042  context.previous_column( block, M, temp );
1043 
1044  // initialize the first band (corresponding to the 0-th row of the DP matrix)
1045  #pragma unroll
1046  for (uint32 j = 0; j <= BAND_LEN; ++j)
1047  band[j] = zero;
1048 
1049  // load a block of entries from the reference
1050  #pragma unroll
1051  for (uint32 t = 0; t < BAND_LEN; ++t)
1052  r_cache[ t ] = ref[ block + t ];
1053 
1054  score_type max_score = Field_traits<score_type>::min();
1055 
1056  score_type temp_i = band[0];
1057 
1058  // loop across the long edge of the DP matrix (i.e. the rows)
1059  for (uint32 i = 0; i < M; ++i)
1060  {
1061  // load the new character from the query
1062  const uint8 q_i = query[i];
1063  const uint8 qq_i = quals[i];
1064 
1065  const int32 m_i = scoring.match(qq_i);
1066 
1067  update_row<false>(
1068  context,
1069  block, N,
1070  i, M,
1071  q_i,
1072  qq_i,
1073  m_i,
1074  r_cache,
1075  temp,
1076  temp_i,
1077  band,
1078  sink,
1079  max_score,
1080  zero,
1081  scoring );
1082  }
1083  if (TYPE == SEMI_GLOBAL)
1084  {
1085  // during semi-global alignment we save the best score across the last row
1086  for (uint32 j = 1; j <= BAND_LEN; ++j)
1087  sink.report( band[j], make_uint2( block+j, M ) );
1088  }
1089  }
1090 
1091  if (window_end == N)
1092  {
1093  const uint32 block = end_block - BAND_LEN;
1094 
1095  // save the previous column
1096  context.previous_column( block, M, temp );
1097 
1098  // initialize the first band (corresponding to the 0-th row of the DP matrix)
1099  #pragma unroll
1100  for (uint32 j = 0; j <= BAND_LEN; ++j)
1101  band[j] = zero;
1102 
1103  // load a block of entries from each query
1104  const uint32 block_end = nvbio::min( block + BAND_LEN, N );
1105  #pragma unroll
1106  for (uint32 t = 0; t < BAND_LEN; ++t)
1107  {
1108  if (block + t < block_end)
1109  r_cache[ t ] = ref[ block + t ];
1110  }
1111 
1112  score_type max_score = Field_traits<score_type>::min();
1113 
1114  score_type temp_i = band[0];
1115 
1116  // loop across the long edge of the DP matrix (i.e. the rows)
1117  for (uint32 i = 0; i < M; ++i)
1118  {
1119  // load the new character from the reference
1120  const uint8 q_i = query[i];
1121  const uint8 qq_i = quals[i];
1122 
1123  const int32 m_i = scoring.match(qq_i);
1124 
1125  update_row<true>(
1126  context,
1127  block, N,
1128  i, M,
1129  q_i,
1130  qq_i,
1131  m_i,
1132  r_cache,
1133  temp,
1134  temp_i,
1135  band,
1136  sink,
1137  max_score,
1138  zero,
1139  scoring );
1140  }
1141 
1142  if (TYPE == SEMI_GLOBAL)
1143  {
1144  // during semi-global alignment we save the best score across the last row
1145  for (uint32 j = 1; j <= BAND_LEN; ++j)
1146  {
1147  if (block + j <= N)
1148  sink.report( band[j], make_uint2( block+j, M ) );
1149  }
1150  }
1151  else if (TYPE == GLOBAL)
1152  {
1153  // during global alignment we save the best score at cell [N][M]
1154  for (uint32 j = 1; j <= BAND_LEN; ++j)
1155  {
1156  if (block + j == N)
1157  sink.report( band[j], make_uint2( block+j, M ) );
1158  }
1159  }
1160  }
1161 
1162  // save the last column
1163  context.last_column( window_end, N, M, temp );
1164  return true;
1165  }
1166 
1241  template <
1242  uint32 MAX_PATTERN_LEN,
1243  typename context_type,
1244  typename string_type,
1245  typename qual_type,
1246  typename ref_type,
1247  typename scoring_type,
1248  typename sink_type>
1249  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1250  static
1251  bool run(
1252  const scoring_type& scoring,
1253  context_type& context,
1254  string_type query,
1255  qual_type quals,
1256  ref_type ref,
1257  const int32 min_score,
1258  sink_type& sink,
1259  const uint32 window_begin,
1260  const uint32 window_end)
1261  {
1262  // instantiated a local memory array
1263  int16 temp[MAX_PATTERN_LEN];
1264  int16* temp_ptr = temp;
1265 
1266  return run(
1267  scoring,
1268  context,
1269  query,
1270  quals,
1271  ref,
1272  min_score,
1273  sink,
1274  window_begin,
1275  window_end,
1276  temp_ptr );
1277  }
1278 };
1279 
1280 template <AlignmentType TYPE, uint32 DIM, typename symbol_type>
1281 struct hamming_bandlen_selector
1282 {
1283  static const uint32 BAND_LEN = 16u / DIM;
1284 };
1285 
1286 template <AlignmentType TYPE, uint32 DIM>
1287 struct hamming_bandlen_selector<TYPE,DIM,simd4u8>
1288 {
1289 #if __CUDA_ARCH__ >= 300
1290  static const uint32 BAND_LEN = 8u;
1291 #else
1292  static const uint32 BAND_LEN = 1u;
1293 #endif
1294 };
1295 
1305 template <
1306  AlignmentType TYPE,
1307  typename scoring_type,
1308  typename algorithm_tag,
1309  typename pattern_string,
1310  typename qual_string,
1311  typename text_string,
1312  typename column_type>
1313 struct alignment_score_dispatch<
1314  HammingDistanceAligner<TYPE,scoring_type,algorithm_tag>,
1315  pattern_string,
1316  qual_string,
1317  text_string,
1318  column_type>
1319 {
1320  typedef HammingDistanceAligner<TYPE,scoring_type,algorithm_tag> aligner_type;
1321 
1334  template <typename sink_type>
1335  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1336  static bool dispatch(
1337  const aligner_type aligner,
1338  const pattern_string pattern,
1339  const qual_string quals,
1340  const text_string text,
1341  const int32 min_score,
1342  sink_type& sink,
1343  column_type column)
1344  {
1345  typedef typename pattern_string::value_type symbol_type;
1346 
1347  NVBIO_VAR_UNUSED const uint32 BAND_LEN = hamming_bandlen_selector<TYPE,1u,symbol_type>::BAND_LEN;
1348 
1349  priv::HammingScoringContext<BAND_LEN,TYPE,algorithm_tag> context;
1350 
1351  const uint32 length = equal<algorithm_tag,PatternBlockingTag>() ? pattern.length() : text.length();
1352 
1353  return hamming_alignment_score_dispatch<BAND_LEN,TYPE,algorithm_tag,symbol_type>::run( aligner.scheme, context, pattern, quals, text, min_score, sink, 0, length, column );
1354  }
1355 
1371  template <
1372  typename sink_type,
1373  typename checkpoint_type>
1374  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1375  static bool dispatch(
1376  const aligner_type aligner,
1377  const pattern_string pattern,
1378  const qual_string quals,
1379  const text_string text,
1380  const int32 min_score,
1381  const uint32 window_begin,
1382  const uint32 window_end,
1383  sink_type& sink,
1384  checkpoint_type checkpoint,
1385  column_type column)
1386  {
1387  typedef typename pattern_string::value_type symbol_type;
1388 
1389  NVBIO_VAR_UNUSED const uint32 BAND_LEN = hamming_bandlen_selector<TYPE,1u,symbol_type>::BAND_LEN;
1390 
1391  priv::HammingCheckpointedScoringContext<BAND_LEN,TYPE,algorithm_tag,checkpoint_type> context( checkpoint );
1392 
1393  return hamming_alignment_score_dispatch<BAND_LEN,TYPE,algorithm_tag,symbol_type>::run( aligner.scheme, context, pattern, quals, text, min_score, sink, window_begin, window_end, column );
1394  }
1395 
1409  template <typename sink_type>
1410  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1411  static bool dispatch(
1412  const aligner_type aligner,
1413  const pattern_string pattern,
1414  const qual_string quals,
1415  const text_string text,
1416  const int32 min_score,
1417  const uint32 window_begin,
1418  const uint32 window_end,
1419  sink_type& sink,
1420  column_type column)
1421  {
1422  typedef typename pattern_string::value_type symbol_type;
1423 
1424  NVBIO_VAR_UNUSED const uint32 BAND_LEN = hamming_bandlen_selector<TYPE,1u,symbol_type>::BAND_LEN;
1425 
1426  priv::HammingScoringContext<BAND_LEN,TYPE,algorithm_tag> context;
1427 
1428  return hamming_alignment_score_dispatch<BAND_LEN,TYPE,algorithm_tag,symbol_type>::run( aligner.scheme, context, pattern, quals, text, min_score, sink, window_begin, window_end, column );
1429  }
1430 };
1431 
1442 template <
1443  uint32 CHECKPOINTS,
1444  AlignmentType TYPE,
1445  typename scoring_type,
1446  typename pattern_string,
1447  typename qual_string,
1448  typename text_string,
1449  typename column_type>
1450 struct alignment_checkpointed_dispatch<
1451  CHECKPOINTS,
1452  HammingDistanceAligner<TYPE,scoring_type>,
1453  pattern_string,
1454  qual_string,
1455  text_string,
1456  column_type>
1457 {
1458  typedef HammingDistanceAligner<TYPE,scoring_type> aligner_type;
1459 
1479  template <
1480  typename sink_type,
1481  typename checkpoint_type>
1482  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1483  static
1484  void dispatch_checkpoints(
1485  const aligner_type aligner,
1486  const pattern_string pattern,
1487  const qual_string quals,
1488  const text_string text,
1489  const int32 min_score,
1490  sink_type& sink,
1491  checkpoint_type checkpoints,
1492  column_type column)
1493  {
1494  typedef typename pattern_string::value_type symbol_type;
1495 
1496  NVBIO_VAR_UNUSED const uint32 BAND_LEN = hamming_bandlen_selector<TYPE,1u,symbol_type>::BAND_LEN;
1497 
1498  priv::HammingCheckpointContext<BAND_LEN,TYPE,CHECKPOINTS,checkpoint_type> context( checkpoints );
1499 
1500  hamming_alignment_score_dispatch<BAND_LEN,TYPE,PatternBlockingTag,symbol_type>::run( aligner.scheme, context, pattern, quals, text, min_score, sink, 0, pattern.length(), column );
1501  }
1502 
1527  template <
1528  typename checkpoint_type,
1529  typename submatrix_type>
1530  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1531  static
1532  uint32 dispatch_submatrix(
1533  const aligner_type aligner,
1534  const pattern_string pattern,
1535  const qual_string quals,
1536  const text_string text,
1537  const int32 min_score,
1538  checkpoint_type checkpoints,
1539  const uint32 checkpoint_id,
1540  submatrix_type submatrix,
1541  column_type column)
1542  {
1543  typedef typename pattern_string::value_type symbol_type;
1544 
1545  NVBIO_VAR_UNUSED const uint32 BAND_LEN = hamming_bandlen_selector<TYPE,1u,symbol_type>::BAND_LEN;
1546 
1547  priv::HammingSubmatrixContext<BAND_LEN,TYPE,CHECKPOINTS,checkpoint_type,submatrix_type>
1548  context( checkpoints, checkpoint_id, submatrix );
1549 
1550  const uint32 window_begin = checkpoint_id * CHECKPOINTS;
1551  const uint32 window_end = nvbio::min( window_begin + CHECKPOINTS, uint32(pattern.length()) );
1552 
1553  NullSink null_sink;
1554  hamming_alignment_score_dispatch<BAND_LEN,TYPE,PatternBlockingTag,symbol_type>::run( aligner.scheme, context, pattern, quals, text, min_score, null_sink, window_begin, window_end, column );
1555 
1556  return window_end - window_begin;
1557  }
1558 };
1559 
1595 template <
1596  uint32 CHECKPOINTS,
1597  AlignmentType TYPE,
1598  typename scoring_type,
1599  typename checkpoint_type,
1600  typename submatrix_type,
1601  typename output_type>
1602 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
1603 bool alignment_traceback(
1604  const HammingDistanceAligner<TYPE,scoring_type> aligner,
1605  checkpoint_type checkpoints,
1606  const uint32 checkpoint_id,
1607  submatrix_type submatrix,
1608  const uint32 submatrix_width,
1609  const uint32 submatrix_height,
1610  uint8& state,
1611  uint2& sink,
1612  output_type& output)
1613 {
1614  //
1615  // Backtrack from the second checkpoint to the first looking up the flow submatrix.
1616  //
1617  int32 current_row = sink.x;
1618  int32 current_col = sink.y - checkpoint_id*CHECKPOINTS - 1u;
1619 
1620  NVBIO_CUDA_DEBUG_ASSERT( current_row > 0 &&
1621  current_row <= submatrix_height, "hamming::alignment_backtrack(): sink (%u,%u) -> local x coordinate %d not in [0,%d[\n", sink.x, sink.y, current_row, submatrix_height );
1622  NVBIO_CUDA_DEBUG_ASSERT( current_col >= 0, "hamming::alignment_backtrack(): sink (%u,%u) -> local y coordinate %d not in [0,%u[ (checkpt %u)\n", sink.x, sink.y, current_col, submatrix_width, checkpoint_id );
1623  NVBIO_CUDA_DEBUG_ASSERT( current_col < submatrix_width, "hamming::alignment_backtrack(): sink (%u,%u) -> local y coordinate %d not in [0,%u[ (checkpt %u)\n", sink.x, sink.y, current_col, submatrix_width, checkpoint_id );
1624 
1625  while (current_row > 0 &&
1626  current_col >= 0)
1627  {
1628  const uint32 submatrix_cell = (current_row-1u) * CHECKPOINTS + current_col;
1629  const uint8 op = submatrix[ submatrix_cell ];
1630 
1631  if (TYPE == LOCAL)
1632  {
1633  if (op == SINK)
1634  {
1635  sink.x = current_row;
1636  sink.y = current_col + checkpoint_id*CHECKPOINTS + 1u;
1637  return true;
1638  }
1639  }
1640 
1641  output.push( op );
1642  }
1643  sink.x = current_row;
1644  sink.y = current_col + checkpoint_id*CHECKPOINTS + 1u;
1645  return current_row ? false : true; // if current_row == 0 we reached the end of the alignment along the text
1646 }
1647 
1649 
1650 } // namespace priv
1651 
1652 } // namespace aln
1653 } // namespace nvbio
1654 
Generated on Wed Feb 25 2015 08:32:53 for NVBIO by   doxygen 1.8.4