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
33// Compares the scalar effective-length convolution against the FFT
34// cross-correlation. The two must agree numerically (there are tests for that);
35// this measures where the crossover in cost lies, which is what justifies keeping
36// both paths.
37fn main() {
38    let args: Vec<String> = std::env::args().collect();
39    let n: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(100_000);
40    let passes: usize = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(3);
41
42    let mut rng = Lcg(0x9E3779B97F4A7C15);
43
44    // Build a realistic transcript-length distribution: log-normal-ish around
45    // ~1.5 kb with a long tail, clamped to [80, 15000] — close to GRCh38 cDNA.
46    let bases = *b"ACGT";
47    let mut seqs: Vec<Vec<u8>> = Vec::with_capacity(n);
48    let mut prefixes: Vec<Vec<u32>> = Vec::with_capacity(n);
49    for _ in 0..n {
50        let u = rng.next_f64();
51        // exp of a normal-ish variate -> heavy right tail
52        let len = (7.3 + 0.9 * (u - 0.5) * 4.0).exp() as usize;
53        let len = len.clamp(80, 15_000);
54        let seq: Vec<u8> = (0..len)
55            .map(|_| bases[(rng.next_u32() & 3) as usize])
56            .collect();
57        prefixes.push(gc_prefix(&seq));
58        seqs.push(seq);
59    }
60
61    // A normalized 3×25 GC ratio model (the shape doesn't matter for timing; we
62    // exercise the same binning + lookup the real run does).
63    let mut obs = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
64    let mut exp = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
65    for ctx in 0..=100 {
66        for gc in 0..=100 {
67            obs.inc(gc, ctx, 1.0 + 0.3 * ((gc + ctx) as f64).sin());
68            exp.inc(gc, ctx, 1.0);
69        }
70    }
71    let gc_model = gc_ratio(&mut obs, &mut exp, 1000.0);
72
73    // Fragment-length CDF: gaussian-ish pmf around mean 250, sd 40, over 0..1000.
74    let fld_max = 1000usize;
75    let mean = 250.0f64;
76    let sd = 40.0f64;
77    let pmf: Vec<f64> = (0..=fld_max)
78        .map(|l| {
79            let z = (l as f64 - mean) / sd;
80            (-0.5 * z * z).exp()
81        })
82        .collect();
83    let (cdf, fld_low, fld_high) = salmon_model::seqbias::fld_cdf_and_bounds(&pmf);
84
85    let mut acc = 0.0f64;
86    let mut best = f64::INFINITY;
87    for p in 0..passes {
88        let t = Instant::now();
89        for (seq, prefix) in seqs.iter().zip(&prefixes) {
90            let ref_len = seq.len() as f64;
91            let elen = (ref_len - 200.0).max(1.0); // ensure unprocessed > 0
92            let bias = BiasInputs {
93                seq: None,
94                gc: Some((&gc_model, salmon_model::GcView::Dense(prefix.as_slice()))),
95                pos: None,
96            };
97            acc += corrected_effective_length_full(
98                seq,
99                &cdf,
100                fld_low,
101                fld_high,
102                &bias,
103                elen,
104                GC_SAMP_STRIDE,
105                false,
106            );
107        }
108        let dt = t.elapsed().as_secs_f64();
109        best = best.min(dt);
110        eprintln!(
111            "pass {p}: {:.3}s  ({:.2} µs/transcript)  acc={}",
112            dt,
113            dt * 1e6 / n as f64,
114            black_box(acc)
115        );
116    }
117    eprintln!(
118        "BEST: {:.3}s over {} transcripts ({:.2} µs/transcript)",
119        best,
120        n,
121        best * 1e6 / n as f64
122    );
123}
conv_bench/conv_bench.rs

conv_bench/
conv_bench.rs
