|
| 1 | +// This Source Code Form is subject to the terms of the Mozilla Public |
| 2 | +// License, v. 2.0. If a copy of the MPL was not distributed with this |
| 3 | +// file, You can obtain one at https://mozilla.org/MPL/2.0/. |
| 4 | + |
| 5 | +// This code implements Welford's method for calculating mean and variance |
| 6 | +// from streaming data as described here: |
| 7 | +// https://jonisalonen.com/2013/deriving-welfords-method-for-computing-variance/ |
| 8 | +// The TLDR is: |
| 9 | +// variance(samples): |
| 10 | +// M := 0 |
| 11 | +// S := 0 |
| 12 | +// for k from 1 to N: |
| 13 | +// x := samples[k] |
| 14 | +// oldM := M |
| 15 | +// M := M + (x-M)/k |
| 16 | +// S := S + (x-M)*(x-oldM) |
| 17 | +// return S/(N-1) |
| 18 | +// |
| 19 | +// where: |
| 20 | +// - M is the mean |
| 21 | +// - S is the squared distance from the mean |
| 22 | + |
| 23 | +use anyhow::{Context, Result, anyhow}; |
| 24 | +use clap::Parser; |
| 25 | +use std::{ |
| 26 | + fs::File, |
| 27 | + io::{self, BufRead, BufReader, Read}, |
| 28 | + path::PathBuf, |
| 29 | +}; |
| 30 | + |
| 31 | +/// Read input samples formatted as a single sample per line. Each sample is |
| 32 | +/// an positive integer that can fit in a u32. This is the same format emitted |
| 33 | +/// by `attest-time`. |
| 34 | +#[derive(Debug, Parser)] |
| 35 | +#[clap(author, version, about, long_about = None)] |
| 36 | +struct Args { |
| 37 | + /// An optional input file. If omitted stdin will be read. |
| 38 | + input: Option<PathBuf>, |
| 39 | +} |
| 40 | + |
| 41 | +/// This strcture holds data used to generate some measures of central |
| 42 | +/// tendency and dispersion |
| 43 | +struct Data { |
| 44 | + /// `count` is a float because it's mostly used as the denominator in the |
| 45 | + /// mean and variance calculation |
| 46 | + count: u32, |
| 47 | + /// collection used to collect input data |
| 48 | + /// TODO: we can get rid of this once we're confident in our impl of the |
| 49 | + /// Welford's algorithms |
| 50 | + durations: Vec<u32>, |
| 51 | + /// `max` is the running max for the dataset |
| 52 | + max: u32, |
| 53 | + /// `mean` holds the running mean calculated w/ Welford's method |
| 54 | + mean: f64, |
| 55 | + /// the running min for the dataset |
| 56 | + min: u32, |
| 57 | + /// `distance_2` is the running squared distance from the mean calculated |
| 58 | + /// w/ Welford's method |
| 59 | + distance_2: f64, |
| 60 | +} |
| 61 | + |
| 62 | +/// impl `Default` manually to set initial value for `min` |
| 63 | +impl Default for Data { |
| 64 | + fn default() -> Self { |
| 65 | + Self { |
| 66 | + count: 0, |
| 67 | + distance_2: 0.0, |
| 68 | + max: 0, |
| 69 | + mean: 0.0, |
| 70 | + min: u32::MAX, |
| 71 | + durations: Vec::new(), |
| 72 | + } |
| 73 | + } |
| 74 | +} |
| 75 | + |
| 76 | +/// my naive and expensive mean calculation |
| 77 | +fn mean(durations: &Vec<u32>, count: u32) -> Result<u32> { |
| 78 | + let mut total: u128 = 0; |
| 79 | + for v in durations { |
| 80 | + // data in the durations collection is u32, this conversion is safe |
| 81 | + total += *v as u128; |
| 82 | + } |
| 83 | + |
| 84 | + let mean = total / count as u128; |
| 85 | + u32::try_from(mean).context("mean u128 to u32") |
| 86 | +} |
| 87 | + |
| 88 | +/// my naive and expensive variance calculation |
| 89 | +fn variance(durations: &Vec<u32>, count: u32, mean: u32) -> Result<u32> { |
| 90 | + // accumulate sum of the squared difference between each sample and the mean |
| 91 | + // this is the numerator in the classic variance equation |
| 92 | + let mut variance: u128 = 0; |
| 93 | + for v in durations { |
| 94 | + let diff = *v as i128 - mean as i128; |
| 95 | + let square = i128::pow(diff, 2); |
| 96 | + variance += square as u128; |
| 97 | + } |
| 98 | + |
| 99 | + let variance = variance / count as u128 - 1; |
| 100 | + Ok(variance as u32) |
| 101 | +} |
| 102 | + |
| 103 | +fn main() -> Result<()> { |
| 104 | + let args = Args::parse(); |
| 105 | + |
| 106 | + // if args.input file not provided use stdin |
| 107 | + let reader: Box<dyn Read> = match args.input { |
| 108 | + Some(i) => Box::new( |
| 109 | + File::open(&i) |
| 110 | + .with_context(|| format!("open file: {}", i.display()))?, |
| 111 | + ), |
| 112 | + None => Box::new(io::stdin()), |
| 113 | + }; |
| 114 | + let reader = BufReader::new(reader); |
| 115 | + |
| 116 | + let mut data = Data::default(); |
| 117 | + |
| 118 | + for line in reader.lines() { |
| 119 | + let line = line.context("read line")?; |
| 120 | + let words: Vec<&str> = line.split_whitespace().collect(); |
| 121 | + if words.len() != 1 { |
| 122 | + return Err(anyhow!("malformed line")); |
| 123 | + } |
| 124 | + |
| 125 | + let micros: u32 = words[0].parse().context("parse u32 from str")?; |
| 126 | + |
| 127 | + if micros < data.min { |
| 128 | + data.min = micros; |
| 129 | + } |
| 130 | + |
| 131 | + if micros > data.max { |
| 132 | + data.max = micros; |
| 133 | + } |
| 134 | + |
| 135 | + data.count = data |
| 136 | + .count |
| 137 | + .checked_add(1) |
| 138 | + .ok_or(anyhow!("too many samples: count overflow"))?; |
| 139 | + |
| 140 | + data.durations.push(micros); |
| 141 | + |
| 142 | + let micros = f64::from(micros); |
| 143 | + let old_mean = data.mean; |
| 144 | + data.mean = data.mean + (micros - data.mean) / (data.count as f64); |
| 145 | + data.distance_2 += (micros - data.mean) * (micros - old_mean); |
| 146 | + } |
| 147 | + |
| 148 | + println!("sample count: {}", data.count); |
| 149 | + println!("min: {}", data.min); |
| 150 | + println!("max: {}", data.max); |
| 151 | + |
| 152 | + // streaming mean, variance, & standard deviation |
| 153 | + { |
| 154 | + println!("welford's mean: {}", data.mean); |
| 155 | + |
| 156 | + // final calculation from welford's method for variance |
| 157 | + // return S/(n-1) |
| 158 | + let variance = data.distance_2 / f64::from(data.count - 1); |
| 159 | + println!("welford's variance: {}", variance); |
| 160 | + |
| 161 | + println!("welford's standard deviation: {}", variance.sqrt()); |
| 162 | + } |
| 163 | + |
| 164 | + // classic mean, variance, & standard deviation |
| 165 | + { |
| 166 | + let mean = mean(&data.durations, data.count) |
| 167 | + .context("calculate mean from dataset")?; |
| 168 | + println!("mean: {mean}"); |
| 169 | + |
| 170 | + let variance = variance(&data.durations, data.count, mean) |
| 171 | + .context("calculate variance from dataset")?; |
| 172 | + println!("variance: {variance}"); |
| 173 | + |
| 174 | + let sqrt = f64::from(variance).sqrt(); |
| 175 | + println!("standard deviation: {}", sqrt); |
| 176 | + } |
| 177 | + |
| 178 | + Ok(()) |
| 179 | +} |
0 commit comments