henad_core/explore/
summary.rs

1//! Statistics over the replicates of each config.
2//!
3//! The 95% confidence interval for a mean is `mean ± t * sd / sqrt(n)`, where `t` is the 97.5% quantile of
4//! Student's t distribution with `n - 1` degrees of freedom.
5
6use std::collections::BTreeMap;
7
8use crate::explore::outcome::RunStatus;
9
10/// Quantiles of Student's t distribution at 97.5% for 1 to 30 degrees of freedom.
11const T_975: [f64; 30] = [
12    12.706_204_736,
13    4.302_652_730,
14    3.182_446_305,
15    2.776_445_105,
16    2.570_581_836,
17    2.446_911_851,
18    2.364_624_252,
19    2.306_004_135,
20    2.262_157_163,
21    2.228_138_852,
22    2.200_985_160,
23    2.178_812_830,
24    2.160_368_656,
25    2.144_786_688,
26    2.131_449_546,
27    2.119_905_299,
28    2.109_815_578,
29    2.100_922_040,
30    2.093_024_054,
31    2.085_963_447,
32    2.079_613_845,
33    2.073_873_068,
34    2.068_657_610,
35    2.063_898_562,
36    2.059_538_553,
37    2.055_529_439,
38    2.051_830_516,
39    2.048_407_142,
40    2.045_229_642,
41    2.042_272_456,
42];
43
44/// Quantile of the standard normal distribution at 97.5%.
45const Z_975: f64 = 1.959_963_984_540_054;
46
47/// Returns the 97.5% quantile of Student's t distribution with `degrees_of_freedom` degrees of freedom.
48///
49/// Up to 30 degrees of freedom the value comes from a table, and past it from the Cornish-Fisher expansion in
50/// `1 / degrees_of_freedom`. Note that the result is infinity for 0 degrees of freedom.
51pub fn student_t_975(degrees_of_freedom: u64) -> f64 {
52    if degrees_of_freedom == 0 {
53        return f64::INFINITY;
54    }
55    if let Some(&quantile) = usize::try_from(degrees_of_freedom - 1)
56        .ok()
57        .and_then(|index| T_975.get(index))
58    {
59        return quantile;
60    }
61    let z = Z_975;
62    let z2 = z * z;
63    let z3 = z2 * z;
64    let z5 = z3 * z2;
65    let z7 = z5 * z2;
66    let z9 = z7 * z2;
67    let g1 = (z3 + z) / 4.0;
68    let g2 = (5.0 * z5 + 16.0 * z3 + 3.0 * z) / 96.0;
69    let g3 = (3.0 * z7 + 19.0 * z5 + 17.0 * z3 - 15.0 * z) / 384.0;
70    let g4 = (79.0 * z9 + 776.0 * z7 + 1482.0 * z5 - 1920.0 * z3 - 945.0 * z) / 92160.0;
71    let inverse = 1.0 / degrees_of_freedom as f64;
72    z + inverse * (g1 + inverse * (g2 + inverse * (g3 + inverse * g4)))
73}
74
75/// Count, mean and spread of a stream of values, updated one value at a time with Welford's method.
76#[derive(Debug, Clone, Copy, Default, PartialEq)]
77pub struct RunningMoments {
78    count: u64,
79    mean: f64,
80    /// Sum of squared distances from the mean.
81    squares: f64,
82}
83
84impl RunningMoments {
85    /// Adds `value`. A value that is not finite is left out.
86    pub fn push(&mut self, value: f64) {
87        if !value.is_finite() {
88            return;
89        }
90        self.count += 1;
91        let delta = value - self.mean;
92        self.mean += delta / self.count as f64;
93        self.squares += delta * (value - self.mean);
94    }
95
96    /// Number of finite values pushed.
97    pub fn count(&self) -> u64 {
98        self.count
99    }
100
101    /// Returns the mean, or `None` before any value.
102    pub fn mean(&self) -> Option<f64> {
103        (self.count > 0).then_some(self.mean)
104    }
105
106    /// Returns the sample standard deviation, with `n - 1` in the denominator, or `None` below two values.
107    pub fn standard_deviation(&self) -> Option<f64> {
108        (self.count > 1).then(|| (self.squares / (self.count - 1) as f64).sqrt())
109    }
110
111    /// Returns the count, mean, standard deviation and 95% confidence interval of the values so far.
112    pub fn summary(&self) -> ReplicateSummary {
113        let standard_deviation = self.standard_deviation();
114        let ci95 = standard_deviation.map(|standard_deviation| {
115            let half_width = student_t_975(self.count - 1) * standard_deviation / (self.count as f64).sqrt();
116            (self.mean - half_width, self.mean + half_width)
117        });
118        ReplicateSummary {
119            n: self.count,
120            mean: self.mean(),
121            standard_deviation,
122            ci95,
123        }
124    }
125}
126
127/// Statistics of one output over the replicates of a config.
128#[derive(Debug, Clone, Copy, PartialEq)]
129pub struct ReplicateSummary {
130    /// Number of finite values.
131    pub n: u64,
132    /// Mean of the finite values, `None` when there are no finite values.
133    pub mean: Option<f64>,
134    /// Sample standard deviation, `None` below two values.
135    pub standard_deviation: Option<f64>,
136    /// Low and high ends of the 95% confidence interval for the mean, `None` below two values.
137    pub ci95: Option<(f64, f64)>,
138}
139
140/// Returns the summary of the finite values in `values`.
141pub fn summarize(values: impl IntoIterator<Item = f64>) -> ReplicateSummary {
142    let mut moments = RunningMoments::default();
143    for value in values {
144        moments.push(value);
145    }
146    moments.summary()
147}
148
149/// Run counts and statistics of one config.
150#[derive(Debug, Clone, PartialEq)]
151pub struct SummaryRow {
152    /// Id of the config.
153    pub config_id: u64,
154    /// Number of runs of the config, whatever their status.
155    pub runs: u64,
156    /// Number of runs whose status is [`RunStatus::Ok`].
157    pub ok: u64,
158    /// Number of runs that ended on a fault or a timeout, as [`RunStatus::is_failure`] decides.
159    pub failed: u64,
160    /// Tick each run ended on.
161    pub ticks: ReplicateSummary,
162    /// One summary per reducer.
163    pub reducers: Vec<ReplicateSummary>,
164}
165
166/// Counts and running statistics of every config, fed one run at a time.
167///
168/// Tick and reducer statistics come from the runs that did not fail.
169#[derive(Debug, Clone, Default, PartialEq)]
170pub struct SummaryAccumulator {
171    reducer_count: usize,
172    configs: BTreeMap<u64, ConfigTally>,
173}
174
175#[derive(Debug, Clone, PartialEq)]
176struct ConfigTally {
177    runs: u64,
178    ok: u64,
179    failed: u64,
180    ticks: RunningMoments,
181    reducers: Vec<RunningMoments>,
182}
183
184impl SummaryAccumulator {
185    /// Returns an accumulator for runs with `reducer_count` reducer values each.
186    pub fn new(reducer_count: usize) -> Self {
187        Self {
188            reducer_count,
189            configs: BTreeMap::new(),
190        }
191    }
192
193    /// Adds a run of config `config_id` that ended with `status` at tick `ticks`, with the values of its reducers.
194    ///
195    /// `reducers` holds one value per reducer. Note that the order in which runs are added affects the statistics down
196    /// to the last bit.
197    pub fn push(&mut self, config_id: u64, status: RunStatus, ticks: u64, reducers: &[Option<f64>]) {
198        let reducer_count = self.reducer_count;
199        let tally = self.configs.entry(config_id).or_insert_with(|| ConfigTally {
200            runs: 0,
201            ok: 0,
202            failed: 0,
203            ticks: RunningMoments::default(),
204            reducers: vec![RunningMoments::default(); reducer_count],
205        });
206        tally.runs += 1;
207        if status == RunStatus::Ok {
208            tally.ok += 1;
209        }
210        if status.is_failure() {
211            tally.failed += 1;
212            return;
213        }
214        tally.ticks.push(ticks as f64);
215        for (moments, value) in tally.reducers.iter_mut().zip(reducers) {
216            if let Some(value) = value {
217                moments.push(*value);
218            }
219        }
220    }
221
222    /// Returns a row for each config with at least one run, in config order.
223    pub fn rows(&self) -> impl Iterator<Item = SummaryRow> + '_ {
224        self.configs.iter().map(|(&config_id, tally)| SummaryRow {
225            config_id,
226            runs: tally.runs,
227            ok: tally.ok,
228            failed: tally.failed,
229            ticks: tally.ticks.summary(),
230            reducers: tally.reducers.iter().map(RunningMoments::summary).collect(),
231        })
232    }
233}
234
235#[cfg(test)]
236mod tests {
237    use super::{RunningMoments, SummaryAccumulator, student_t_975, summarize};
238    use crate::explore::outcome::RunStatus;
239
240    fn close(actual: f64, expected: f64, tolerance: f64) -> bool {
241        (actual - expected).abs() <= tolerance
242    }
243
244    #[test]
245    fn t_quantiles_match_the_table() {
246        for (df, expected) in [(1, 12.706), (2, 4.303), (10, 2.228), (30, 2.042), (120, 1.980)] {
247            let quantile = student_t_975(df);
248            assert!(close(quantile, expected, 1e-3), "df {df}: {quantile}");
249        }
250        assert!(
251            close(student_t_975(31), 2.039_513, 1e-6),
252            "the expansion joins the table"
253        );
254        assert!(
255            close(student_t_975(1_000_000), 1.959_966, 1e-6),
256            "and tends to the normal quantile"
257        );
258        assert!(student_t_975(0).is_infinite());
259    }
260
261    #[test]
262    fn a_replicate_summary_matches_hand_computed_statistics() {
263        let summary = summarize([1.0, 2.0, 3.0, 4.0]);
264        assert_eq!(summary.n, 4);
265        assert_eq!(summary.mean, Some(2.5));
266        let sd = summary.standard_deviation.expect("four values have a spread");
267        assert!(close(sd, 1.2910, 1e-4), "sd {sd}");
268        let (low, high) = summary.ci95.expect("four values have an interval");
269        assert!(close(high - 2.5, 2.054, 1e-3), "half-width {}", high - 2.5);
270        assert!(close(2.5 - low, high - 2.5, 1e-12), "the interval is symmetric");
271    }
272
273    #[test]
274    fn a_single_replicate_has_no_spread() {
275        let summary = summarize([7.0]);
276        assert_eq!(summary.n, 1);
277        assert_eq!(summary.mean, Some(7.0));
278        assert_eq!(summary.standard_deviation, None);
279        assert_eq!(summary.ci95, None);
280        let empty = summarize([]);
281        assert_eq!((empty.n, empty.mean), (0, None));
282    }
283
284    #[test]
285    fn values_that_are_not_finite_are_left_out() {
286        let mut moments = RunningMoments::default();
287        for value in [1.0, f64::NAN, 3.0, f64::INFINITY] {
288            moments.push(value);
289        }
290        assert_eq!(moments.count(), 2);
291        assert_eq!(moments.mean(), Some(2.0));
292    }
293
294    #[test]
295    fn configs_are_summarized_apart_and_failed_runs_are_counted_only() {
296        let mut accumulator = SummaryAccumulator::new(1);
297        accumulator.push(1, RunStatus::Ok, 100, &[Some(4.0)]);
298        accumulator.push(0, RunStatus::Ok, 100, &[Some(1.0)]);
299        accumulator.push(0, RunStatus::NonFinite, 100, &[Some(3.0)]);
300        accumulator.push(0, RunStatus::Panicked, 40, &[Some(99.0)]);
301        accumulator.push(0, RunStatus::Ok, 100, &[None]);
302        let rows: Vec<_> = accumulator.rows().collect();
303        assert_eq!(rows.len(), 2);
304        let first = &rows[0];
305        assert_eq!(first.config_id, 0, "rows come in config order");
306        assert_eq!((first.runs, first.ok, first.failed), (4, 2, 1));
307        assert_eq!(first.ticks.mean, Some(100.0), "the failed run's tick is left out");
308        assert_eq!(first.reducers[0].n, 2, "an empty value is left out");
309        assert_eq!(first.reducers[0].mean, Some(2.0));
310        assert_eq!(rows[1].reducers[0].mean, Some(4.0));
311    }
312}
henad_core/explore/summary.rs

henad_core/explore/
summary.rs
