conv_bench/
conv_bench.rs

1//! Isolated microbenchmark for the GC bias-corrected effective-length
2//! convolution (`corrected_effective_length_full`), so it can be profiled and
3//! A/B'd without the surrounding 130s quant run (where it's only ~5% of wall).
4//!
5//! Generates a realistic set of transcripts (human-cDNA-like length spread,
6//! random ACGT), a 3×25 GC ratio model, and a fragment-length CDF, then times
7//! the convolution over every transcript. Run:
8//!   cargo run --release --example conv_bench -- [num_transcripts] [passes]
9
10use std::hint::black_box;
11use std::time::Instant;
12
13use salmon_model::gcbias::{DEFAULT_COND_BINS, DEFAULT_GC_BINS};
14use salmon_model::{
15    corrected_effective_length_full, gc_prefix, gc_ratio, BiasInputs, GcFragModel, GC_SAMP_STRIDE,
16};
17
18/// Tiny deterministic LCG so the bench is reproducible without a rand dep.
19struct Lcg(u64);
20impl Lcg {
21    fn next_u32(&mut self) -> u32 {
22        self.0 = self
23            .0
24            .wrapping_mul(6364136223846793005)
25            .wrapping_add(1442695040888963407);
26        (self.0 >> 33) as u32
27    }
28    fn next_f64(&mut self) -> f64 {
29        self.next_u32() as f64 / u32::MAX as f64
30    }
31}
32
33fn main() {
34    let args: Vec<String> = std::env::args().collect();
35    let n: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(100_000);
36    let passes: usize = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(3);
37
38    let mut rng = Lcg(0x9E3779B97F4A7C15);
39
40    // Build a realistic transcript-length distribution: log-normal-ish around
41    // ~1.5 kb with a long tail, clamped to [80, 15000] — close to GRCh38 cDNA.
42    let bases = [b'A', b'C', b'G', b'T'];
43    let mut seqs: Vec<Vec<u8>> = Vec::with_capacity(n);
44    let mut prefixes: Vec<Vec<u32>> = Vec::with_capacity(n);
45    for _ in 0..n {
46        let u = rng.next_f64();
47        // exp of a normal-ish variate -> heavy right tail
48        let len = (7.3 + 0.9 * (u - 0.5) * 4.0).exp() as usize;
49        let len = len.clamp(80, 15_000);
50        let seq: Vec<u8> = (0..len)
51            .map(|_| bases[(rng.next_u32() & 3) as usize])
52            .collect();
53        prefixes.push(gc_prefix(&seq));
54        seqs.push(seq);
55    }
56
57    // A normalized 3×25 GC ratio model (the shape doesn't matter for timing; we
58    // exercise the same binning + lookup the real run does).
59    let mut obs = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
60    let mut exp = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
61    for ctx in 0..=100 {
62        for gc in 0..=100 {
63            obs.inc(gc, ctx, 1.0 + 0.3 * ((gc + ctx) as f64).sin());
64            exp.inc(gc, ctx, 1.0);
65        }
66    }
67    let gc_model = gc_ratio(&mut obs, &mut exp, 1000.0);
68
69    // Fragment-length CDF: gaussian-ish pmf around mean 250, sd 40, over 0..1000.
70    let fld_max = 1000usize;
71    let mean = 250.0f64;
72    let sd = 40.0f64;
73    let pmf: Vec<f64> = (0..=fld_max)
74        .map(|l| {
75            let z = (l as f64 - mean) / sd;
76            (-0.5 * z * z).exp()
77        })
78        .collect();
79    let (cdf, fld_low, fld_high) = salmon_model::seqbias::fld_cdf_and_bounds(&pmf);
80
81    let mut acc = 0.0f64;
82    let mut best = f64::INFINITY;
83    for p in 0..passes {
84        let t = Instant::now();
85        for (seq, prefix) in seqs.iter().zip(&prefixes) {
86            let ref_len = seq.len() as f64;
87            let elen = (ref_len - 200.0).max(1.0); // ensure unprocessed > 0
88            let bias = BiasInputs {
89                seq: None,
90                gc: Some((&gc_model, salmon_model::GcView::Dense(prefix.as_slice()))),
91                pos: None,
92            };
93            acc += corrected_effective_length_full(
94                seq,
95                &cdf,
96                fld_low,
97                fld_high,
98                &bias,
99                elen,
100                GC_SAMP_STRIDE,
101                false,
102            );
103        }
104        let dt = t.elapsed().as_secs_f64();
105        best = best.min(dt);
106        eprintln!(
107            "pass {p}: {:.3}s  ({:.2} µs/transcript)  acc={}",
108            dt,
109            dt * 1e6 / n as f64,
110            black_box(acc)
111        );
112    }
113    eprintln!(
114        "BEST: {:.3}s over {} transcripts ({:.2} µs/transcript)",
115        best,
116        n,
117        best * 1e6 / n as f64
118    );
119}
conv_bench/conv_bench.rs

conv_bench/
conv_bench.rs
