Counts : speed improvment on materliaze Bin

src/scan/scan.rs

@@ -697,48 +697,72 @@ pub fn scan_contig(
 
     // Materialize BinCount rows in the same order/ranges as defs
     let mut out = Vec::with_capacity(defs.len());
-
+    let contig = contig.to_string();
     for (def_idx, def) in defs.into_iter().enumerate() {
         let b = &mut bins[def_idx];
+        let len = b.depths_fwd.len();
+
+        assert_eq!(len, b.depths_rev.len());
+        assert_eq!(len, b.lowq.len());
+        assert_eq!(len, b.fb_invs.len());
 
         let n_reads = b.reps.len() as u32;
-        let n_sa = b.reps.values().filter(|r| r.has_sa).count() as u32;
+        let mut n_sa = 0_u32;
+        let region_len_u32 = def
+            .end
+            .checked_sub(def.start)
+            .and_then(|v| v.checked_add(1))
+            .expect("invalid region");
+        let region_len = region_len_u32 as usize;
+        let mut start_counts = vec![0_u32; region_len];
+        let mut end_counts = vec![0_u32; region_len];
 
-        let mut start_counts: HashMap<u32, u32> = HashMap::new();
-        let mut end_counts: HashMap<u32, u32> = HashMap::new();
         for r in b.reps.values() {
-            if (def.start..=def.end).contains(&r.aln_start) {
-                *start_counts.entry(r.aln_start).or_insert(0) += 1;
+            n_sa += u32::from(r.has_sa);
+
+            if let Some(offset) = r.aln_start.checked_sub(def.start) {
+                if offset < region_len_u32 {
+                    start_counts[offset as usize] += 1;
+                }
             }
-            if (def.start..=def.end).contains(&r.aln_end) {
-                *end_counts.entry(r.aln_end).or_insert(0) += 1;
+
+            if let Some(offset) = r.aln_end.checked_sub(def.start) {
+                if offset < region_len_u32 {
+                    end_counts[offset as usize] += 1;
+                }
             }
         }
-        let max_start = start_counts.values().copied().max().unwrap_or(0);
-        let max_end = end_counts.values().copied().max().unwrap_or(0);
-
-        // total depths + per-position sens ratio, reconstructed from fwd/rev
-        let mut depths: Vec<u32> = Vec::with_capacity(b.depths_fwd.len());
-        let mut depths_sens_ratio: Vec<f64> = Vec::with_capacity(b.depths_fwd.len());
-        let mut sum_fwd: u64 = 0;
-        let mut sum_rev: u64 = 0;
-        for (&fwd, &rev) in b.depths_fwd.iter().zip(b.depths_rev.iter()) {
+
+        let max_start = start_counts.into_iter().max().unwrap_or(0);
+        let max_end = end_counts.into_iter().max().unwrap_or(0);
+
+        let mut depths = vec![0_u32; len];
+        let mut depths_sens_ratio = vec![0.0_f64; len];
+
+        let mut sum_fwd = 0_u64;
+        let mut sum_rev = 0_u64;
+
+        for i in 0..len {
+            let fwd = b.depths_fwd[i];
+            let rev = b.depths_rev[i];
             let total = fwd + rev;
-            depths.push(total);
-            depths_sens_ratio.push(if total == 0 {
+
+            depths[i] = total;
+            depths_sens_ratio[i] = if total == 0 {
                 f64::NAN
             } else {
-                fwd as f64 / total as f64
-            });
-            sum_fwd += fwd as u64;
-            sum_rev += rev as u64;
+                f64::from(fwd) / f64::from(total)
+            };
+
+            sum_fwd += u64::from(fwd);
+            sum_rev += u64::from(rev);
         }
-        let sum_depths = sum_fwd + sum_rev;
 
-        let coverage = if depths.is_empty() {
+        let sum_depths = sum_fwd + sum_rev;
+        let coverage = if len == 0 {
             0.0
         } else {
-            sum_depths as f64 / (depths.len() as f64)
+            sum_depths as f64 / len as f64
         };
 
         let n_reads_f = n_reads as f64;
@@ -752,28 +776,21 @@ pub fn scan_contig(
             )
         };
 
-        // NEW 1A: ratio_lowq = sum(low_qualities) / sum(depths)
-        let sum_lowq: u64 = b.lowq.iter().map(|&v| v as u64).sum();
+        // ratio_lowq = sum(low_qualities) / sum(depths)
+        let sum_lowq: u64 = b.lowq.iter().copied().map(u64::from).sum();
         let ratio_lowq = if sum_depths == 0 {
             f64::NAN
         } else {
             sum_lowq as f64 / sum_depths as f64
         };
 
-        // NEW 2C (depth-normalized): ratio_fb = max(fb_invs[i]) / depths[i]
+        // ratio_fb = max_i(fb_invs[i] / depths[i])
         let ratio_fb = b
             .fb_invs
             .iter()
-            .enumerate()
-            .max_by_key(|&(_, v)| *v)
-            .map(|(idx, &max_val)| {
-                let d = depths[idx];
-                if d == 0 {
-                    f64::NAN
-                } else {
-                    max_val as f64 / d as f64
-                }
-            })
+            .zip(&depths)
+            .filter_map(|(&fb, &depth)| (depth != 0).then_some(f64::from(fb) / f64::from(depth)))
+            .reduce(f64::max)
             .unwrap_or(f64::NAN);
 
         // NEW 3A: ratio_sens = pooled sum(fwd) / sum(fwd+rev), bin-level
@@ -784,7 +801,7 @@ pub fn scan_contig(
         };
 
         out.push(BinCount {
-            contig: contig.to_string(),
+            contig: contig.clone(),
             start: def.start,
             end: def.end,
             n_reads,
@@ -1328,7 +1345,7 @@ mod tests {
         let end0 = 32_323_039;
         let start0 = 32_177_773 - 10;
         let end0 = start0 + 1000;
-        let case = "CHALO";
+        let case = "CALO";
         let time = "diag";
 
         let config = Config::default();