fix(test): the unmapped-reason test re-implemented the logic it was testing

CHANGELOG.md

@@ -131,3 +131,15 @@ Sections commonly used: Features, Bug fixes, Other changes.
   needs random access to the SA in RAM.
 
 Initial release of Rust rewrite of STAR.
+### Other changes
+
+- The per-read transcript-filter bookkeeping is a fixed set of counters
+  rather than a `HashMap<&str, i32>`, removing a heap allocation per read
+  and a string hash per filtered transcript. Output is unchanged.
+
+- The splice-motif check in the junction scan carries a sliding window
+  and looks the motif up in a table instead of re-reading four genome
+  bases and matching on them at every position. Output is unchanged
+  (SAM byte-identical on 200k reads); the scan itself goes from 8.5 ns
+  to 5.2 ns per iteration, measured by `examples/jrbench.rs`.
+

examples/jrbench.rs

@@ -0,0 +1,97 @@
+//! Is the junction scan stalled on memory, and is it the acceptor stream?
+//!
+//! The stack sampler gives self time, not cache misses, so this asks the
+//! question by construction instead: run the same scan over the same genome,
+//! changing only how far the acceptor sits from the donor. Everything else
+//! (iteration count, branch pattern, instruction mix) is held fixed, because
+//! `del` enters the loop only as an address offset once it is inside the
+//! intron-length range.
+//!
+//! If time per iteration is flat across `del`, the loop is not memory-bound on
+//! the acceptor and prefetching it would be wasted work. If it climbs with
+//! `del`, that stream is the cost.
+//!
+//! Run: cargo run --release --example jrbench
+
+use rustar_aligner::align::score::AlignmentScorer;
+use rustar_aligner::genome::{Genome, GenomeSeq};
+use std::time::Instant;
+
+/// Deterministic pseudo-random bases, so the run is reproducible and the motif
+/// hit rate is the same at every `del`.
+fn synthetic_genome(n: usize) -> Genome {
+    let mut state = 0x2545_F491_4F6C_DD1Du64;
+    let mut seq = Vec::with_capacity(n);
+    for _ in 0..n {
+        state ^= state << 13;
+        state ^= state >> 7;
+        state ^= state << 17;
+        seq.push((state >> 33) as u8 & 3);
+    }
+    Genome {
+        transform_blocks: None,
+        sequence: GenomeSeq::Owned(seq),
+        n_genome: n as u64,
+        n_genome_real: n as u64,
+        n_chr_real: 1,
+        chr_name: vec!["chr1".to_string()],
+        chr_length: vec![n as u64],
+        chr_start: vec![0, n as u64],
+    }
+}
+
+fn main() {
+    // Large enough that the donor and acceptor streams cannot both stay in
+    // cache when they are far apart.
+    const N: usize = 256 << 20; // 256 MB of bases, one byte each
+    const READ_LEN: usize = 150;
+    const SCANS: usize = 20_000;
+
+    let genome = synthetic_genome(N);
+    let mut scorer = AlignmentScorer::from_params_minimal();
+    scorer.align_intron_min = 20;
+    scorer.align_intron_max = u32::MAX;
+
+    let mut read = vec![0u8; READ_LEN];
+    let mut state = 0x9E37_79B9_7F4A_7C15u64;
+    for b in read.iter_mut() {
+        state ^= state << 13;
+        state ^= state >> 7;
+        state ^= state << 17;
+        *b = (state >> 33) as u8 & 3;
+    }
+
+    println!("{:>12}  {:>10}  {:>12}", "del", "wall", "ns/iter");
+    for del in [64u64, 1_024, 16_384, 262_144, 4_194_304, 67_108_864] {
+        // Spread the scans over the genome so no single window stays resident.
+        let stride = (N as u64 - del - 4096) / SCANS as u64;
+        let iters_per_scan = 100usize; // r_gap 0 + next_seed_len 100
+        let t0 = Instant::now();
+        let mut sink = 0i64;
+        for k in 0..SCANS {
+            let g_a_end = 2048 + k as u64 * stride;
+            let (jr, _motif, score, _l, _r) = scorer.find_best_junction_position(
+                &read,
+                75,
+                g_a_end,
+                0,
+                del as i64,
+                &genome,
+                false,
+                N as u64,
+                75,
+                iters_per_scan,
+            );
+            sink += jr as i64 + score as i64;
+        }
+        let dt = t0.elapsed();
+        let iters = (SCANS * iters_per_scan) as f64;
+        println!(
+            "{:>12}  {:>9.3}s  {:>11.2}   (sink {})",
+            del,
+            dt.as_secs_f64(),
+            dt.as_secs_f64() * 1e9 / iters,
+            sink
+        );
+    }
+}

src/align/read_align.rs

@@ -12,6 +12,81 @@ use crate::params::{IntronMotifFilter, IntronStrandFilter, MultimapperOrder, Par
 use crate::stats::UnmappedReason;
 use std::hash::{DefaultHasher, Hash, Hasher};
 
+/// Why a transcript was dropped by the quality filters.
+///
+/// A fixed set of counters rather than a `HashMap<&str, i32>`: this is built
+/// once per read and incremented per filtered transcript, so the map cost a
+/// heap allocation per read and a string hash per increment. The counts feed
+/// two debug logs and the `unmapped_reason` decision below, all of which ask
+/// only whether a count is non-zero, so nothing here can reach the output
+/// except through that decision.
+#[derive(Default, Clone, Copy)]
+struct FilterReasons {
+    counts: [u32; FilterReasons::N],
+}
+
+impl FilterReasons {
+    const N: usize = 9;
+    const NAMES: [&'static str; Self::N] = [
+        "score_min",
+        "score_min_relative",
+        "mismatch_max",
+        "mismatch_rate",
+        "match_min",
+        "match_min_relative",
+        "noncanonical_junction",
+        "noncanonical_unannotated_junction",
+        "inconsistent_strand",
+    ];
+
+    const SCORE_MIN: usize = 0;
+    const SCORE_MIN_RELATIVE: usize = 1;
+    const MISMATCH_MAX: usize = 2;
+    const MISMATCH_RATE: usize = 3;
+    const MATCH_MIN: usize = 4;
+    const MATCH_MIN_RELATIVE: usize = 5;
+    const NONCANONICAL_JUNCTION: usize = 6;
+    const NONCANONICAL_UNANNOTATED_JUNCTION: usize = 7;
+    const INCONSISTENT_STRAND: usize = 8;
+
+    #[inline]
+    fn add(&mut self, reason: usize) {
+        self.counts[reason] += 1;
+    }
+
+    #[inline]
+    fn has(&self, reason: usize) -> bool {
+        self.counts[reason] > 0
+    }
+
+    /// Any reason other than the two mismatch ones.
+    #[inline]
+    fn has_non_mismatch(&self) -> bool {
+        self.counts
+            .iter()
+            .enumerate()
+            .any(|(i, &c)| c > 0 && i != Self::MISMATCH_MAX && i != Self::MISMATCH_RATE)
+    }
+}
+
+impl std::fmt::Debug for FilterReasons {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        let mut first = true;
+        write!(f, "{{")?;
+        for (i, &c) in self.counts.iter().enumerate() {
+            if c == 0 {
+                continue;
+            }
+            if !first {
+                write!(f, ", ")?;
+            }
+            first = false;
+            write!(f, "{}: {}", Self::NAMES[i], c)?;
+        }
+        write!(f, "}}")
+    }
+}
+
 /// Derive a deterministic per-read RNG seed from `run_rng_seed` + the read name.
 ///
 /// STAR seeds `std::mt19937` once per chunk/thread (`runRNGseed*(iChunk+1)`),
@@ -409,24 +484,24 @@ pub fn align_read(
 
     // Log filtering statistics
     let pre_filter_count = transcripts.len();
-    let mut filter_reasons = std::collections::HashMap::new();
+    let mut filter_reasons = FilterReasons::default();
 
     transcripts.retain(|t| {
         // Absolute score threshold
         if t.score < params.out_filter_score_min {
-            *filter_reasons.entry("score_min").or_insert(0) += 1;
+            filter_reasons.add(FilterReasons::SCORE_MIN);
             return false;
         }
 
         // Relative score threshold: STAR casts to intScore (i32)
         if t.score < (params.out_filter_score_min_over_lread * lread_m1) as i32 {
-            *filter_reasons.entry("score_min_relative").or_insert(0) += 1;
+            filter_reasons.add(FilterReasons::SCORE_MIN_RELATIVE);
             return false;
         }
 
         // Absolute mismatch count
         if t.n_mismatch > params.out_filter_mismatch_nmax {
-            *filter_reasons.entry("mismatch_max").or_insert(0) += 1;
+            filter_reasons.add(FilterReasons::MISMATCH_MAX);
             log::debug!(
                 "Filtered {}: {} mismatches > {} max (read_len={}, score={})",
                 read_name,
@@ -441,7 +516,7 @@ pub fn align_read(
         // Relative mismatch count (mismatches / read_length)
         let mismatch_rate = t.n_mismatch as f64 / read_length;
         if mismatch_rate > params.out_filter_mismatch_nover_lmax {
-            *filter_reasons.entry("mismatch_rate").or_insert(0) += 1;
+            filter_reasons.add(FilterReasons::MISMATCH_RATE);
             log::debug!(
                 "Filtered {}: {:.1}% mismatch rate > {:.1}% max ({}/{} bases, score={})",
                 read_name,
@@ -457,13 +532,13 @@ pub fn align_read(
         // Absolute matched bases
         let n_matched = t.n_matched();
         if n_matched < params.out_filter_match_nmin as usize {
-            *filter_reasons.entry("match_min").or_insert(0) += 1;
+            filter_reasons.add(FilterReasons::MATCH_MIN);
             return false;
         }
 
         // Relative matched bases: STAR casts to uint (u32)
         if (n_matched as f64) < params.out_filter_match_nmin_over_lread * lread_m1 {
-            *filter_reasons.entry("match_min_relative").or_insert(0) += 1;
+            filter_reasons.add(FilterReasons::MATCH_MIN_RELATIVE);
             return false;
         }
 
@@ -475,7 +550,7 @@ pub fn align_read(
             IntronMotifFilter::RemoveNoncanonical => {
                 // Reject if any junction is non-canonical
                 if t.junction_motifs.contains(&SpliceMotif::NonCanonical) {
-                    *filter_reasons.entry("noncanonical_junction").or_insert(0) += 1;
+                    filter_reasons.add(FilterReasons::NONCANONICAL_JUNCTION);
                     return false;
                 }
             }
@@ -486,9 +561,7 @@ pub fn align_read(
                     .zip(t.junction_annotated.iter())
                     .any(|(m, annotated)| *m == SpliceMotif::NonCanonical && !annotated)
                 {
-                    *filter_reasons
-                        .entry("noncanonical_unannotated_junction")
-                        .or_insert(0) += 1;
+                    filter_reasons.add(FilterReasons::NONCANONICAL_UNANNOTATED_JUNCTION);
                     return false;
                 }
             }
@@ -512,7 +585,7 @@ pub fn align_read(
                 }
             }
             if has_plus && has_minus {
-                *filter_reasons.entry("inconsistent_strand").or_insert(0) += 1;
+                filter_reasons.add(FilterReasons::INCONSISTENT_STRAND);
                 return false;
             }
         }
@@ -639,11 +712,9 @@ pub fn align_read(
     // if the quality filter emptied the set, the best transcript failed → classify
     // short/mismatch; otherwise, if too many loci survive, it's multi; else mapped.
     let unmapped_reason = if transcripts.is_empty() {
-        let has_mismatch = filter_reasons.contains_key("mismatch_max")
-            || filter_reasons.contains_key("mismatch_rate");
-        let has_other = filter_reasons
-            .keys()
-            .any(|k| *k != "mismatch_max" && *k != "mismatch_rate");
+        let has_mismatch = filter_reasons.has(FilterReasons::MISMATCH_MAX)
+            || filter_reasons.has(FilterReasons::MISMATCH_RATE);
+        let has_other = filter_reasons.has_non_mismatch();
         if has_mismatch && !has_other {
             Some(UnmappedReason::TooManyMismatches)
         } else {
@@ -1708,14 +1779,17 @@ mod tests {
     fn test_unmapped_reason_mismatch_classification() {
         // Verify the filter_reasons logic used to derive unmapped_reason.
         // mismatch-only → TooManyMismatches; anything else (or mixed) → TooShort.
-        let classify = |reasons: &[&str]| -> UnmappedReason {
-            let filter_reasons: std::collections::HashMap<&str, i32> =
-                reasons.iter().map(|k| (*k, 1)).collect();
-            let has_mismatch = filter_reasons.contains_key("mismatch_max")
-                || filter_reasons.contains_key("mismatch_rate");
-            let has_other = filter_reasons
-                .keys()
-                .any(|k| *k != "mismatch_max" && *k != "mismatch_rate");
+        // Drives `FilterReasons` itself rather than re-implementing the
+        // decision: the previous version of this test kept its own copy of the
+        // logic, so it would have passed even if the shipped code had changed.
+        let classify = |reasons: &[usize]| -> UnmappedReason {
+            let mut filter_reasons = FilterReasons::default();
+            for &r in reasons {
+                filter_reasons.add(r);
+            }
+            let has_mismatch = filter_reasons.has(FilterReasons::MISMATCH_MAX)
+                || filter_reasons.has(FilterReasons::MISMATCH_RATE);
+            let has_other = filter_reasons.has_non_mismatch();
             if has_mismatch && !has_other {
                 UnmappedReason::TooManyMismatches
             } else {
@@ -1724,25 +1798,31 @@ mod tests {
         };
 
         assert_eq!(
-            classify(&["mismatch_max"]),
+            classify(&[FilterReasons::MISMATCH_MAX]),
             UnmappedReason::TooManyMismatches
         );
         assert_eq!(
-            classify(&["mismatch_rate"]),
+            classify(&[FilterReasons::MISMATCH_RATE]),
             UnmappedReason::TooManyMismatches
         );
         assert_eq!(
-            classify(&["mismatch_max", "mismatch_rate"]),
+            classify(&[FilterReasons::MISMATCH_MAX, FilterReasons::MISMATCH_RATE]),
             UnmappedReason::TooManyMismatches
         );
         // Mixed with score filter → TooShort
         assert_eq!(
-            classify(&["mismatch_max", "score_min"]),
+            classify(&[FilterReasons::MISMATCH_MAX, FilterReasons::SCORE_MIN]),
             UnmappedReason::TooShort
         );
         // Score-only → TooShort
-        assert_eq!(classify(&["score_min"]), UnmappedReason::TooShort);
-        assert_eq!(classify(&["match_min"]), UnmappedReason::TooShort);
+        assert_eq!(
+            classify(&[FilterReasons::SCORE_MIN]),
+            UnmappedReason::TooShort
+        );
+        assert_eq!(
+            classify(&[FilterReasons::MATCH_MIN]),
+            UnmappedReason::TooShort
+        );
     }
 
     #[test]