1//! Descriptive statistics and bootstrap confidence intervals for timing samples.2//! Floats are used here only for statistics (never for field arithmetic).34use rand::{RngCore, SeedableRng};5use rand_chacha::ChaCha20Rng;6use serde::Serialize;78/// Linear-interpolated quantile of a sorted slice (same convention as numpy `linear`).9pub fn quantile_sorted(sorted: &[f64], q: f64) -> f64 {10 if sorted.is_empty() {11 return f64::NAN;12 }13 let pos = q * (sorted.len() as f64 - 1.0);14 let lo = pos.floor() as usize;15 let hi = pos.ceil() as usize;16 if lo == hi {17 sorted[lo]18 } else {19 sorted[lo] + (sorted[hi] - sorted[lo]) * (pos - lo as f64)20 }21}2223pub fn median(xs: &[f64]) -> f64 {24 let mut v = xs.to_vec();25 v.sort_by(|a, b| a.partial_cmp(b).unwrap());26 quantile_sorted(&v, 0.5)27}2829#[derive(Serialize, Clone, Debug)]30pub struct Timing {31 pub reps: usize,32 pub ops_per_rep: usize,33 pub ns_per_op_samples: Vec<f64>,34 pub median_ns: f64,35 pub q1_ns: f64,36 pub q3_ns: f64,37 pub iqr_ns: f64,38 pub min_ns: f64,39 pub mean_ns: f64,40 pub sd_ns: f64,41 pub ops_per_s: f64,42}4344impl Timing {45 pub fn from_samples(ops_per_rep: usize, samples: Vec<f64>) -> Timing {46 let mut v = samples.clone();47 v.sort_by(|a, b| a.partial_cmp(b).unwrap());48 let n = v.len() as f64;49 let mean = v.iter().sum::<f64>() / n;50 let sd = if v.len() > 1 {51 (v.iter().map(|x| (x - mean) * (x - mean)).sum::<f64>() / (n - 1.0)).sqrt()52 } else {53 0.054 };55 let med = quantile_sorted(&v, 0.5);56 let q1 = quantile_sorted(&v, 0.25);57 let q3 = quantile_sorted(&v, 0.75);58 Timing {59 reps: samples.len(),60 ops_per_rep,61 ns_per_op_samples: samples,62 median_ns: med,63 q1_ns: q1,64 q3_ns: q3,65 iqr_ns: q3 - q1,66 min_ns: v[0],67 mean_ns: mean,68 sd_ns: sd,69 ops_per_s: 1e9 / med,70 }71 }72}7374#[derive(Serialize, Clone, Debug)]75pub struct RatioCI {76 pub point: f64,77 pub ci95_low: f64,78 pub ci95_high: f64,79 pub bootstrap_resamples: usize,80}8182/// Bootstrap 95 % CI for the ratio median(num) / median(den), resampling each sample set with83/// replacement (percentile method). Deterministic given `seed`.84pub fn bootstrap_ratio_of_medians(num: &[f64], den: &[f64], resamples: usize, seed: u64) -> RatioCI {85 let mut rng = ChaCha20Rng::seed_from_u64(seed);86 let mut ratios = Vec::with_capacity(resamples);87 let mut bn = vec![0.0; num.len()];88 let mut bd = vec![0.0; den.len()];89 for _ in 0..resamples {90 for x in bn.iter_mut() {91 *x = num[(rng.next_u64() % num.len() as u64) as usize];92 }93 for x in bd.iter_mut() {94 *x = den[(rng.next_u64() % den.len() as u64) as usize];95 }96 ratios.push(median(&bn) / median(&bd));97 }98 ratios.sort_by(|a, b| a.partial_cmp(b).unwrap());99 RatioCI {100 point: median(num) / median(den),101 ci95_low: quantile_sorted(&ratios, 0.025),102 ci95_high: quantile_sorted(&ratios, 0.975),103 bootstrap_resamples: resamples,104 }105}106