block_aligner/
scores.rs

1//! Structs for representing match/mismatch scoring matrices.
2
3#[cfg(feature = "simd_sse2")]
4use crate::sse2::*;
5
6#[cfg(feature = "simd_avx2")]
7use crate::avx2::*;
8
9#[cfg(feature = "simd_wasm")]
10use crate::simd128::*;
11
12#[cfg(feature = "simd_neon")]
13use crate::neon::*;
14
15use std::i8;
16
17pub trait Matrix {
18    /// Byte to use as padding.
19    const NULL: u8;
20    /// Create a new matrix with default (usually nonsense) values.
21    ///
22    /// Use `new_simple` to create a sensible scoring matrix.
23    fn new() -> Self;
24    /// Set the score for a pair of bytes.
25    fn set(&mut self, a: u8, b: u8, score: i8);
26    /// Get the score for a pair of bytes.
27    fn get(&self, a: u8, b: u8) -> i8;
28    /// Get the pointer for a specific index.
29    fn as_ptr(&self, i: usize) -> *const i8;
30    /// Get the scores for a certain byte and a certain SIMD vector of bytes.
31    unsafe fn get_scores(&self, c: u8, v: HalfSimd, right: bool) -> Simd;
32    /// Convert a byte to a better storage format that makes retrieving scores
33    /// easier.
34    fn convert_char(c: u8) -> u8;
35}
36
37/// Amino acid scoring matrix.
38///
39/// Supports characters `A` to `Z`. Lowercase characters are uppercased.
40#[repr(C, align(32))]
41#[derive(Clone, PartialEq, Debug)]
42pub struct AAMatrix {
43    scores: [i8; 27 * 32]
44}
45
46impl AAMatrix {
47    /// Create a simple matrix with a certain match and mismatch score.
48    pub const fn new_simple(match_score: i8, mismatch_score: i8) -> Self {
49        let mut scores = [i8::MIN; 27 * 32];
50        let mut i = b'A';
51        while i <= b'Z' {
52            let mut j = b'A';
53            while j <= b'Z' {
54                let idx = ((i - b'A') as usize) * 32 + ((j - b'A') as usize);
55                scores[idx] = if i == j { match_score } else { mismatch_score };
56                j += 1;
57            }
58            i += 1;
59        }
60        Self { scores }
61    }
62
63    /// Create an AAMatrix from a tab-separated table with no headers.
64    ///
65    /// Use `aa_order` to pass in the amino acids in order.
66    pub fn from_tsv(tsv: &str, aa_order: &str) -> Self {
67        let tsv = tsv.trim();
68        let aa_order = aa_order.split_ascii_whitespace().map(|s| s.as_bytes()[0]).collect::<Vec<_>>();
69        let mut res = Self::new();
70
71        for (line, &a) in tsv.split("\n").zip(&aa_order) {
72            for (score, &b) in line.split_ascii_whitespace().zip(&aa_order) {
73                let score = score.parse::<i8>().unwrap();
74                res.set(a, b, score);
75            }
76        }
77
78        res
79    }
80}
81
82impl Matrix for AAMatrix {
83    const NULL: u8 = b'A' + 26u8;
84
85    fn new() -> Self {
86        Self { scores: [i8::MIN; 27 * 32] }
87    }
88
89    fn set(&mut self, a: u8, b: u8, score: i8) {
90        let a = a.to_ascii_uppercase();
91        let b = b.to_ascii_uppercase();
92        assert!(b'A' <= a && a <= b'Z' + 1);
93        assert!(b'A' <= b && b <= b'Z' + 1);
94        let idx = ((a - b'A') as usize) * 32 + ((b - b'A') as usize);
95        self.scores[idx] = score;
96        let idx = ((b - b'A') as usize) * 32 + ((a - b'A') as usize);
97        self.scores[idx] = score;
98    }
99
100    fn get(&self, a: u8, b: u8) -> i8 {
101        let a = a.to_ascii_uppercase();
102        let b = b.to_ascii_uppercase();
103        assert!(b'A' <= a && a <= b'Z' + 1);
104        assert!(b'A' <= b && b <= b'Z' + 1);
105        let idx = ((a - b'A') as usize) * 32 + ((b - b'A') as usize);
106        self.scores[idx]
107    }
108
109    #[inline]
110    fn as_ptr(&self, i: usize) -> *const i8 {
111        debug_assert!(i < 27);
112        unsafe { self.scores.as_ptr().add(i * 32) }
113    }
114
115    // TODO: get rid of lookup for around half of the shifts by constructing position specific scoring matrix?
116    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
117    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
118    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
119    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
120    #[inline]
121    unsafe fn get_scores(&self, c: u8, v: HalfSimd, _right: bool) -> Simd {
122        // efficiently lookup scores for each character in v
123        let matrix_ptr = self.as_ptr(c as usize);
124        let scores1 = lutsimd_load(matrix_ptr as *const LutSimd);
125        let scores2 = lutsimd_load((matrix_ptr as *const LutSimd).add(1));
126        halfsimd_lookup2_i16(scores1, scores2, v)
127    }
128
129    #[inline]
130    fn convert_char(c: u8) -> u8 {
131        let c = c.to_ascii_uppercase();
132        assert!(c >= b'A' && c <= Self::NULL);
133        c - b'A'
134    }
135}
136
137/// Nucleotide scoring matrix.
138///
139/// Supports characters `A`, `C`, `G`, `N`, and `T`. Lowercase characters are uppercased.
140///
141/// If a larger alphabet is needed (for example, with IUPAC characters), use `AAMatrix` instead.
142#[repr(C, align(32))]
143#[derive(Clone, PartialEq, Debug)]
144pub struct NucMatrix {
145    scores: [i8; 8 * 16]
146}
147
148impl NucMatrix {
149    /// Create a simple matrix with a certain match and mismatch score.
150    pub const fn new_simple(match_score: i8, mismatch_score: i8) -> Self {
151        let mut scores = [i8::MIN; 8 * 16];
152        let alpha = [b'A', b'T', b'C', b'G', b'N'];
153        let mut i = 0;
154        while i < alpha.len() {
155            let mut j = 0;
156            while j < alpha.len() {
157                let idx = ((alpha[i] & 0b111) as usize) * 16 + ((alpha[j] & 0b1111) as usize);
158                scores[idx] = if i == j { match_score } else { mismatch_score };
159                j += 1;
160            }
161            i += 1;
162        }
163        Self { scores }
164    }
165}
166
167impl Matrix for NucMatrix {
168    const NULL: u8 = b'Z';
169
170    fn new() -> Self {
171        Self { scores: [i8::MIN; 8 * 16] }
172    }
173
174    fn set(&mut self, a: u8, b: u8, score: i8) {
175        let a = a.to_ascii_uppercase();
176        let b = b.to_ascii_uppercase();
177        assert!(b'A' <= a && a <= b'Z');
178        assert!(b'A' <= b && b <= b'Z');
179        let idx = ((a & 0b111) as usize) * 16 + ((b & 0b1111) as usize);
180        self.scores[idx] = score;
181        let idx = ((b & 0b111) as usize) * 16 + ((a & 0b1111) as usize);
182        self.scores[idx] = score;
183    }
184
185    fn get(&self, a: u8, b: u8) -> i8 {
186        let a = a.to_ascii_uppercase();
187        let b = b.to_ascii_uppercase();
188        assert!(b'A' <= a && a <= b'Z');
189        assert!(b'A' <= b && b <= b'Z');
190        let idx = ((a & 0b111) as usize) * 16 + ((b & 0b1111) as usize);
191        self.scores[idx]
192    }
193
194    #[inline]
195    fn as_ptr(&self, i: usize) -> *const i8 {
196        unsafe { self.scores.as_ptr().add((i & 0b111) * 16) }
197    }
198
199    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
200    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
201    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
202    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
203    #[inline]
204    unsafe fn get_scores(&self, c: u8, v: HalfSimd, _right: bool) -> Simd {
205        // efficiently lookup scores for each character in v
206        let matrix_ptr = self.as_ptr(c as usize);
207        let scores = lutsimd_load(matrix_ptr as *const LutSimd);
208        halfsimd_lookup1_i16(scores, v)
209    }
210
211    #[inline]
212    fn convert_char(c: u8) -> u8 {
213        let c = c.to_ascii_uppercase();
214        assert!(c >= b'A' && c <= Self::NULL);
215        c
216    }
217}
218
219/// Arbitrary bytes scoring matrix.
220#[repr(C)]
221#[derive(Clone, PartialEq, Debug)]
222pub struct ByteMatrix {
223    match_score: i8,
224    mismatch_score: i8
225}
226
227impl ByteMatrix {
228    /// Create a simple matrix with a certain match and mismatch score.
229    pub const fn new_simple(match_score: i8, mismatch_score: i8) -> Self {
230        Self { match_score, mismatch_score }
231    }
232}
233
234impl Matrix for ByteMatrix {
235    /// May lead to inaccurate results with x drop alignment,
236    /// if the block reaches the ends of the strings.
237    ///
238    /// Avoid using `ByteMatrix` with x drop alignment.
239    const NULL: u8 = b'\0';
240
241    fn new() -> Self {
242        Self { match_score: i8::MIN, mismatch_score: i8::MIN }
243    }
244
245    fn set(&mut self, _a: u8, _b: u8, _score: i8) {
246        unimplemented!();
247    }
248
249    fn get(&self, a: u8, b: u8) -> i8 {
250        if a == b { self.match_score } else { self.mismatch_score }
251    }
252
253    #[inline]
254    fn as_ptr(&self, _i: usize) -> *const i8 {
255        unimplemented!()
256    }
257
258    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
259    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
260    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
261    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
262    #[inline]
263    unsafe fn get_scores(&self, c: u8, v: HalfSimd, _right: bool) -> Simd {
264        let match_scores = halfsimd_set1_i8(self.match_score);
265        let mismatch_scores = halfsimd_set1_i8(self.mismatch_score);
266        halfsimd_lookup_bytes_i16(match_scores, mismatch_scores, halfsimd_set1_i8(c as i8), v)
267    }
268
269    #[inline]
270    fn convert_char(c: u8) -> u8 {
271        c
272    }
273}
274
275/// Match = 1, mismatch = -1.
276#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
277pub static NW1: NucMatrix = NucMatrix::new_simple(1, -1);
278
279#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
280pub static BLOSUM45: AAMatrix = AAMatrix { scores: include!("../matrices/BLOSUM45") };
281
282#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
283pub static BLOSUM50: AAMatrix = AAMatrix { scores: include!("../matrices/BLOSUM50") };
284
285#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
286pub static BLOSUM62: AAMatrix = AAMatrix { scores: include!("../matrices/BLOSUM62") };
287
288#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
289pub static BLOSUM80: AAMatrix = AAMatrix { scores: include!("../matrices/BLOSUM80") };
290
291#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
292pub static BLOSUM90: AAMatrix = AAMatrix { scores: include!("../matrices/BLOSUM90") };
293
294#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
295pub static PAM100: AAMatrix = AAMatrix { scores: include!("../matrices/PAM100") };
296
297#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
298pub static PAM120: AAMatrix = AAMatrix { scores: include!("../matrices/PAM120") };
299
300#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
301pub static PAM160: AAMatrix = AAMatrix { scores: include!("../matrices/PAM160") };
302
303#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
304pub static PAM200: AAMatrix = AAMatrix { scores: include!("../matrices/PAM200") };
305
306#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
307pub static PAM250: AAMatrix = AAMatrix { scores: include!("../matrices/PAM250") };
308
309/// Match = 1, mismatch = -1.
310#[cfg_attr(not(target_arch = "wasm32"), no_mangle)]
311pub static BYTES1: ByteMatrix = ByteMatrix::new_simple(1, -1);
312
313/*pub trait ScoreParams {
314    const GAP_OPEN: i8;
315    const GAP_EXTEND: i8;
316    const I: usize;
317}
318
319pub struct Params<const GAP_OPEN: i8, const GAP_EXTEND: i8, const I: usize>;
320
321impl<const GAP_OPEN: i8, const GAP_EXTEND: i8, const I: usize> ScoreParams for Params<{ GAP_OPEN }, { GAP_EXTEND }, { I }> {
322    const GAP_OPEN: i8 = GAP_OPEN;
323    const GAP_EXTEND: i8 = GAP_EXTEND;
324    const I: usize = I;
325}
326
327pub type GapParams<const GAP_OPEN: i8, const GAP_EXTEND: i8> = Params<{ GAP_OPEN }, { GAP_EXTEND }, 0>;*/
328
329/// Open and extend gap costs.
330///
331/// Open cost must include the extend cost. For example, with `Gaps { open: -11, extend: -1 }`,
332/// a gap of length 1 costs -11, and a gap of length 2 costs -12.
333#[derive(Copy, Clone, PartialEq, Debug)]
334#[repr(C)]
335pub struct Gaps {
336    pub open: i8,
337    pub extend: i8
338}
339
340#[allow(non_snake_case)]
341pub trait Profile {
342    /// Byte to use as padding.
343    const NULL: u8;
344
345    /// Create a new profile of a specific length, with default (large negative) values.
346    ///
347    /// Note that internally, the created profile is longer than a conventional position-specific scoring
348    /// matrix (and `str_len`) by 1, so the profile will have the same length as the number of
349    /// columns in the DP matrix.
350    /// The first column of scores in the profile should be large negative values (padding).
351    /// This allows gap open costs to be specified for the first column of the DP matrix.
352    fn new(str_len: usize, block_size: usize, gap_extend: i8) -> Self;
353    /// Create a new profile from a byte string.
354    fn from_bytes(b: &[u8], block_size: usize, match_score: i8, mismatch_score: i8, gap_open_C: i8, gap_close_C: i8, gap_open_R: i8, gap_extend: i8) -> Self;
355
356    /// Get the length of the profile.
357    fn len(&self) -> usize;
358    /// Clear the profile so it can be reused for profile lengths less than or equal
359    /// to the length this struct was created with.
360    fn clear(&mut self, str_len: usize, block_size: usize);
361
362    /// Set the score for a position and byte.
363    ///
364    /// The profile should be first `clear`ed before it is reused with different lengths.
365    ///
366    /// The first column (`i = 0`) should be padded with large negative values.
367    /// Therefore, set values starting from `i = 1`.
368    fn set(&mut self, i: usize, b: u8, score: i8);
369    /// Set the scores for all positions in the position specific scoring matrix.
370    ///
371    /// The profile should be first `clear`ed before it is reused with different lengths.
372    ///
373    /// Use `order` to specify the order of bytes that is used in the `scores` matrix.
374    /// Scores (in `scores`) should be stored in row-major order, where each row is a different position
375    /// and each column is a different byte.
376    ///
377    /// Use `left_shift` and `right_shift` to scale all the scores.
378    fn set_all(&mut self, order: &[u8], scores: &[i8], left_shift: usize, right_shift: usize);
379    /// Set the scores for all positions in reverse in the position specific scoring matrix.
380    ///
381    /// The profile should be first `clear`ed before it is reused with different lengths.
382    ///
383    /// Use `order` to specify the order of bytes that is used in the `scores` matrix.
384    /// Scores (in `scores`) should be stored in row-major order, where each row is a different position
385    /// and each column is a different byte.
386    ///
387    /// Use `left_shift` and `right_shift` to scale all the scores.
388    fn set_all_rev(&mut self, order: &[u8], scores: &[i8], left_shift: usize, right_shift: usize);
389
390    /// Set the gap open cost for a column.
391    ///
392    /// When aligning a sequence `q` to a profile `r`, this is the gap open cost at column `i` for a
393    /// column transition in the DP matrix with `|q| + 1` rows and `|r| + 1` columns.
394    /// This represents starting a gap in `q`.
395    fn set_gap_open_C(&mut self, i: usize, gap: i8);
396    /// Set the gap close cost for a column.
397    ///
398    /// When aligning a sequence `q` to a profile `r`, this is the gap close cost at column `i` for
399    /// ending column transitions in the DP matrix with `|q| + 1` rows and `|r| + 1` columns.
400    /// This represents ending a gap in `q`.
401    fn set_gap_close_C(&mut self, i: usize, gap: i8);
402    /// Set the gap open cost for a row.
403    ///
404    /// When aligning a sequence `q` to a profile `r`, this is the gap open cost at column `i` for
405    /// a row transition in the DP matrix with `|q| + 1` rows and `|r| + 1` columns.
406    /// This represents starting a gap in `r`.
407    fn set_gap_open_R(&mut self, i: usize, gap: i8);
408
409    /// Set the gap open cost for all column transitions.
410    fn set_all_gap_open_C(&mut self, gap: i8);
411    /// Set the gap close cost for all column transitions.
412    fn set_all_gap_close_C(&mut self, gap: i8);
413    /// Set the gap open cost for all row transitions.
414    fn set_all_gap_open_R(&mut self, gap: i8);
415
416    /// Get the score for a position and byte.
417    fn get(&self, i: usize, b: u8) -> i8;
418    /// Get the gap extend cost.
419    fn get_gap_extend(&self) -> i8;
420    /// Get the pointer for a specific index.
421    fn as_ptr_pos(&self, i: usize) -> *const i8;
422    /// Get the pointer for a specific amino acid.
423    fn as_ptr_aa(&self, a: usize) -> *const i16;
424
425    /// Get the scores for a certain SIMD vector of bytes at a specific position in the profile.
426    unsafe fn get_scores_pos(&self, i: usize, v: HalfSimd, right: bool) -> Simd;
427    /// Get the scores for a certain byte starting at a specific position in the profile.
428    unsafe fn get_scores_aa(&self, i: usize, c: u8, right: bool) -> Simd;
429
430    /// Get the gap open cost for a column.
431    unsafe fn get_gap_open_right_C(&self, i: usize) -> Simd;
432    /// Get the gap close cost for a column.
433    unsafe fn get_gap_close_right_C(&self, i: usize) -> Simd;
434    /// Get the gap open cost for a row.
435    unsafe fn get_gap_open_right_R(&self, i: usize) -> Simd;
436
437    /// Get the gap open cost for a column.
438    unsafe fn get_gap_open_down_C(&self, i: usize) -> Simd;
439    /// Get the gap close cost for a column.
440    unsafe fn get_gap_close_down_C(&self, i: usize) -> Simd;
441    /// Get the gap open cost for a row.
442    unsafe fn get_gap_open_down_R(&self, i: usize) -> Simd;
443
444    /// Convert a byte to a better storage format that makes retrieving scores
445    /// easier.
446    fn convert_char(c: u8) -> u8;
447}
448
449/// Amino acid position specific scoring matrix.
450///
451/// Supports characters `A` to `Z`. Lowercase characters are uppercased.
452#[allow(non_snake_case)]
453#[derive(Clone, PartialEq, Debug)]
454pub struct AAProfile {
455    aa_pos: Vec<i16>,
456    pos_aa: Vec<i8>,
457    gap_extend: i8,
458    pos_gap_open_C: Vec<i16>,
459    pos_gap_close_C: Vec<i16>,
460    pos_gap_open_R: Vec<i16>,
461    // length used for underlying allocated vectors
462    max_len: usize,
463    // length used for the current padded profile
464    curr_len: usize,
465    // length of the profile without padding (same length as the consensus sequence of the position
466    // specific scoring matrix)
467    str_len: usize
468}
469
470impl Profile for AAProfile {
471    const NULL: u8 = b'A' + 26u8;
472
473    fn new(str_len: usize, block_size: usize, gap_extend: i8) -> Self {
474        let max_len = str_len + block_size + 1;
475        Self {
476            aa_pos: vec![i8::MIN as i16; 32 * max_len],
477            pos_aa: vec![i8::MIN; max_len * 32],
478            gap_extend,
479            pos_gap_open_C: vec![i8::MIN as i16; max_len],
480            pos_gap_close_C: vec![i8::MIN as i16; max_len],
481            pos_gap_open_R: vec![i8::MIN as i16; max_len],
482            max_len,
483            curr_len: max_len,
484            str_len
485        }
486    }
487
488    #[allow(non_snake_case)]
489    fn from_bytes(b: &[u8], block_size: usize, match_score: i8, mismatch_score: i8, gap_open_C: i8, gap_close_C: i8, gap_open_R: i8, gap_extend: i8) -> Self {
490        let mut res = Self::new(b.len(), block_size, gap_extend);
491
492        for i in 0..b.len() {
493            for c in b'A'..=b'Z' {
494                res.set(i + 1, c, if c == b[i] { match_score } else { mismatch_score });
495            }
496        }
497
498        for i in 0..b.len() + 1 {
499            res.set_gap_open_C(i, gap_open_C);
500            res.set_gap_close_C(i, gap_close_C);
501            res.set_gap_open_R(i, gap_open_R);
502        }
503
504        res
505    }
506
507    fn len(&self) -> usize {
508        self.str_len
509    }
510
511    fn clear(&mut self, str_len: usize, block_size: usize) {
512        let curr_len = str_len + block_size + 1;
513        assert!(curr_len <= self.max_len);
514        self.aa_pos[..32 * curr_len].fill(i8::MIN as i16);
515        self.pos_aa[..curr_len * 32].fill(i8::MIN);
516        self.pos_gap_open_C[..curr_len].fill(i8::MIN as i16);
517        self.pos_gap_close_C[..curr_len].fill(i8::MIN as i16);
518        self.pos_gap_open_R[..curr_len].fill(i8::MIN as i16);
519        self.str_len = str_len;
520        self.curr_len = curr_len;
521    }
522
523    fn set(&mut self, i: usize, b: u8, score: i8) {
524        let b = b.to_ascii_uppercase();
525        assert!(b'A' <= b && b <= b'Z' + 1);
526        let idx = i * 32 + ((b - b'A') as usize);
527        self.pos_aa[idx] = score;
528        let idx = ((b - b'A') as usize) * self.curr_len + i;
529        self.aa_pos[idx] = score as i16;
530    }
531
532    fn set_all(&mut self, order: &[u8], scores: &[i8], left_shift: usize, right_shift: usize) {
533        self.set_all_core::<false>(order, scores, left_shift, right_shift);
534    }
535
536    fn set_all_rev(&mut self, order: &[u8], scores: &[i8], left_shift: usize, right_shift: usize) {
537        self.set_all_core::<true>(order, scores, left_shift, right_shift);
538    }
539
540    fn set_gap_open_C(&mut self, i: usize, gap: i8) {
541        assert!(gap < 0, "Gap open cost must be negative!");
542        self.pos_gap_open_C[i] = gap as i16;
543    }
544
545    fn set_gap_close_C(&mut self, i: usize, gap: i8) {
546        self.pos_gap_close_C[i] = gap as i16;
547    }
548
549    fn set_gap_open_R(&mut self, i: usize, gap: i8) {
550        assert!(gap < 0, "Gap open cost must be negative!");
551        self.pos_gap_open_R[i] = gap as i16;
552    }
553
554    fn set_all_gap_open_C(&mut self, gap: i8) {
555        assert!(gap < 0, "Gap open cost must be negative!");
556        self.pos_gap_open_C[..self.str_len + 1].fill(gap as i16);
557    }
558
559    fn set_all_gap_close_C(&mut self, gap: i8) {
560        self.pos_gap_close_C[..self.str_len + 1].fill(gap as i16);
561    }
562
563    fn set_all_gap_open_R(&mut self, gap: i8) {
564        assert!(gap < 0, "Gap open cost must be negative!");
565        self.pos_gap_open_R[..self.str_len + 1].fill(gap as i16);
566    }
567
568    fn get(&self, i: usize, b: u8) -> i8 {
569        let b = b.to_ascii_uppercase();
570        assert!(b'A' <= b && b <= b'Z' + 1);
571        let idx = i * 32 + ((b - b'A') as usize);
572        self.pos_aa[idx]
573    }
574
575    fn get_gap_extend(&self) -> i8 {
576        self.gap_extend
577    }
578
579    #[inline]
580    fn as_ptr_pos(&self, i: usize) -> *const i8 {
581        debug_assert!(i < self.curr_len);
582        unsafe { self.pos_aa.as_ptr().add(i * 32) }
583    }
584
585    #[inline]
586    fn as_ptr_aa(&self, a: usize) -> *const i16 {
587        debug_assert!(a < 27);
588        unsafe { self.aa_pos.as_ptr().add(a * self.curr_len) }
589    }
590
591    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
592    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
593    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
594    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
595    #[inline]
596    unsafe fn get_scores_pos(&self, i: usize, v: HalfSimd, _right: bool) -> Simd {
597        // efficiently lookup scores for each character in v
598        let matrix_ptr = self.as_ptr_pos(i);
599        let scores1 = lutsimd_loadu(matrix_ptr as *const LutSimd);
600        let scores2 = lutsimd_loadu((matrix_ptr as *const LutSimd).add(1));
601        halfsimd_lookup2_i16(scores1, scores2, v)
602    }
603
604    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
605    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
606    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
607    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
608    #[inline]
609    unsafe fn get_scores_aa(&self, i: usize, c: u8, _right: bool) -> Simd {
610        let matrix_ptr = self.as_ptr_aa(c as usize);
611        simd_loadu(matrix_ptr.add(i) as *const Simd)
612    }
613
614    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
615    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
616    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
617    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
618    #[inline]
619    unsafe fn get_gap_open_right_C(&self, i: usize) -> Simd {
620        simd_set1_i16(*self.pos_gap_open_C.as_ptr().add(i))
621    }
622
623    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
624    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
625    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
626    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
627    #[inline]
628    unsafe fn get_gap_close_right_C(&self, i: usize) -> Simd {
629        simd_set1_i16(*self.pos_gap_close_C.as_ptr().add(i))
630    }
631
632    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
633    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
634    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
635    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
636    #[inline]
637    unsafe fn get_gap_open_right_R(&self, i: usize) -> Simd {
638        simd_set1_i16(*self.pos_gap_open_R.as_ptr().add(i))
639    }
640
641    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
642    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
643    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
644    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
645    #[inline]
646    unsafe fn get_gap_open_down_C(&self, i: usize) -> Simd {
647        simd_loadu(self.pos_gap_open_C.as_ptr().add(i) as *const Simd)
648    }
649
650    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
651    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
652    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
653    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
654    #[inline]
655    unsafe fn get_gap_close_down_C(&self, i: usize) -> Simd {
656        simd_loadu(self.pos_gap_close_C.as_ptr().add(i) as *const Simd)
657    }
658
659    #[cfg_attr(feature = "simd_sse2", target_feature(enable = "sse2"))]
660    #[cfg_attr(feature = "simd_avx2", target_feature(enable = "avx2"))]
661    #[cfg_attr(feature = "simd_wasm", target_feature(enable = "simd128"))]
662    #[cfg_attr(feature = "simd_neon", target_feature(enable = "neon"))]
663    #[inline]
664    unsafe fn get_gap_open_down_R(&self, i: usize) -> Simd {
665        simd_loadu(self.pos_gap_open_R.as_ptr().add(i) as *const Simd)
666    }
667
668    #[inline]
669    fn convert_char(c: u8) -> u8 {
670        let c = c.to_ascii_uppercase();
671        assert!(c >= b'A' && c <= Self::NULL);
672        c - b'A'
673    }
674}
675
676impl AAProfile {
677    fn set_all_core<const REV: bool>(&mut self, order: &[u8], scores: &[i8], left_shift: usize, right_shift: usize) {
678        #[repr(align(32))]
679        struct A([u8; 32]);
680        let mut o = A([Self::NULL - b'A'; 32]);
681        assert!(order.len() <= 32);
682
683        for (i, &b) in order.iter().enumerate() {
684            let b = b.to_ascii_uppercase();
685            assert!(b'A' <= b && b <= b'Z' + 1);
686            o.0[i] = b - b'A';
687        }
688        assert_eq!(scores.len() / order.len(), self.str_len);
689
690        let mut i = if REV { self.str_len } else { 1 };
691        let mut score_idx = 0;
692
693        while if REV { i >= 1 } else { i <= self.str_len } {
694            let mut j = 0;
695
696            while j < order.len() {
697                unsafe {
698                    let score = ((*scores.as_ptr().add(score_idx)) << left_shift) >> right_shift;
699                    let b = *o.0.as_ptr().add(j) as usize;
700                    *self.pos_aa.as_mut_ptr().add(i * 32 + b) = score;
701                    *self.aa_pos.as_mut_ptr().add(b * self.curr_len + i) = score as i16;
702                }
703
704                score_idx += 1;
705                j += 1;
706            }
707
708            if REV {
709                i -= 1;
710            } else {
711                i += 1;
712            }
713        }
714    }
715}
block_aligner/scores.rs

block_aligner/
scores.rs
