Skip to main content

ta_benchmarks/
catalogue_statistics.rs

1//! Deterministic timing statistics for Catalogue measurements.
2
3use std::cmp::Ordering;
4
5#[derive(Clone, Debug, PartialEq)]
6pub struct TimingStats {
7    pub median_ns: f64,
8    pub ci95_lower_ns: f64,
9    pub ci95_upper_ns: f64,
10    pub throughput_observations_per_second: f64,
11    pub sample_count: usize,
12    pub outlier_count: usize,
13    pub outlier_low_count: usize,
14    pub outlier_high_count: usize,
15}
16
17pub fn timing_stats(samples_ns: &[f64], input_length: usize) -> Result<TimingStats, String> {
18    if samples_ns.len() < 2 {
19        return Err("at least two timing samples are required".to_owned());
20    }
21    if samples_ns
22        .iter()
23        .any(|sample| !sample.is_finite() || *sample <= 0.0)
24    {
25        return Err("timing samples must be positive and finite".to_owned());
26    }
27    let mut sorted = samples_ns.to_vec();
28    sorted.sort_by(|left, right| left.partial_cmp(right).unwrap_or(Ordering::Equal));
29    let median_ns = quantile(&sorted, 0.5);
30    let q1 = quantile(&sorted, 0.25);
31    let q3 = quantile(&sorted, 0.75);
32    let iqr = q3 - q1;
33    let low_fence = q1 - 1.5 * iqr;
34    let high_fence = q3 + 1.5 * iqr;
35    let outlier_low_count = sorted
36        .iter()
37        .take_while(|sample| **sample < low_fence)
38        .count();
39    let outlier_high_count = sorted
40        .iter()
41        .rev()
42        .take_while(|sample| **sample > high_fence)
43        .count();
44    let (ci95_lower_ns, ci95_upper_ns) = bootstrap_median_ci(samples_ns);
45    Ok(TimingStats {
46        median_ns,
47        ci95_lower_ns,
48        ci95_upper_ns,
49        throughput_observations_per_second: input_length as f64 * 1.0e9 / median_ns,
50        sample_count: samples_ns.len(),
51        outlier_count: outlier_low_count + outlier_high_count,
52        outlier_low_count,
53        outlier_high_count,
54    })
55}
56
57fn bootstrap_median_ci(samples: &[f64]) -> (f64, f64) {
58    const REPLICATES: usize = 10_000;
59    let mut state = 0x636174616c6f6755_u64 ^ samples.len() as u64;
60    let mut resample = vec![0.0; samples.len()];
61    let mut medians = Vec::with_capacity(REPLICATES);
62    for _ in 0..REPLICATES {
63        for value in &mut resample {
64            state ^= state << 13;
65            state ^= state >> 7;
66            state ^= state << 17;
67            *value = samples[(state as usize) % samples.len()];
68        }
69        resample.sort_by(|left, right| left.partial_cmp(right).unwrap_or(Ordering::Equal));
70        medians.push(quantile(&resample, 0.5));
71    }
72    medians.sort_by(|left, right| left.partial_cmp(right).unwrap_or(Ordering::Equal));
73    (quantile(&medians, 0.025), quantile(&medians, 0.975))
74}
75
76fn quantile(sorted: &[f64], probability: f64) -> f64 {
77    let position = probability * (sorted.len() - 1) as f64;
78    let lower = position.floor() as usize;
79    let upper = position.ceil() as usize;
80    let fraction = position - lower as f64;
81    sorted[lower] + (sorted[upper] - sorted[lower]) * fraction
82}
83pub(crate) fn validate_positive_timing_evidence(
84    median_ns: f64,
85    ci95_lower_ns: f64,
86    ci95_upper_ns: f64,
87    throughput_observations_per_second: f64,
88    sample_count: usize,
89    input_length: usize,
90) -> Result<(), String> {
91    if [
92        median_ns,
93        ci95_lower_ns,
94        ci95_upper_ns,
95        throughput_observations_per_second,
96    ]
97    .into_iter()
98    .any(|value| !value.is_finite() || value <= 0.0)
99    {
100        return Err("timing evidence must be positive and finite".to_owned());
101    }
102    if ci95_lower_ns > median_ns || median_ns > ci95_upper_ns {
103        return Err("95% confidence interval must contain the median".to_owned());
104    }
105    if sample_count == 0 {
106        return Err("sample_count must be positive".to_owned());
107    }
108    if input_length == 0 {
109        return Err("input_length must be positive".to_owned());
110    }
111    let expected_throughput = input_length as f64 * 1.0e9 / median_ns;
112    let relative_error =
113        (throughput_observations_per_second - expected_throughput).abs() / expected_throughput;
114    if relative_error > 1.0e-4 {
115        return Err(format!(
116            "throughput is incoherent with input_length and median_ns (relative error {relative_error:.6})"
117        ));
118    }
119    Ok(())
120}