harmos_signal/waveform/

summary.rs

//! What a run of samples came to, in a form that survives being split up.

/// The statistics of a run of samples, as accumulators that merge.
///
/// Every field is additive or extremal, so the summary of two runs is a
/// function of their two summaries and never of their samples: a window folded
/// from ten sub-windows is the same value as the window folded from every
/// sample at once. That is what makes this usable beneath a tiered store, where
/// the coarse level is folded from the fine one and the samples themselves are
/// long gone.
///
/// The accumulators are deliberately the ones harmos's `Bucket` carries —
/// `count`, `sum`, `sum_sq`, and the two extremes — so the two merge by one
/// algebra and a summary's numbers are a bucket's, minus the axis coordinates a
/// kernel has no notion of. Nothing here depends on harmos; the compatibility
/// is in the arithmetic.
///
/// Variance is the exception, and the reason this is a reshape rather than a
/// port: `sum_sq / n − mean²` loses every significant digit when the mean is
/// large compared with the spread, and returns garbage — often negative. A
/// Welford sum of squared deviations is carried alongside instead, merged by
/// Chan's formula, and the variance is read off that.
///
/// ```
/// use harmos_signal::Summary;
///
/// let whole = Summary::of(&[1.0, 2.0, 3.0, 4.0]);
/// let split = Summary::of(&[1.0, 2.0]).merge(&Summary::of(&[3.0, 4.0]));
///
/// assert_eq!(whole, split);
/// assert_eq!(whole.mean(), Some(2.5));
/// ```
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Summary {
    count: u64,
    sum: f64,
    sum_sq: f64,
    lowest: f64,
    highest: f64,
    deviations: f64,
}

impl Default for Summary {
    fn default() -> Self {
        Self {
            count: 0,
            sum: 0.0,
            sum_sq: 0.0,
            lowest: f64::INFINITY,
            highest: f64::NEG_INFINITY,
            deviations: 0.0,
        }
    }
}

impl Summary {
    /// A summary of nothing, which is the identity of [`merge`](Self::merge).
    #[must_use]
    pub fn new() -> Self {
        Self::default()
    }

    /// The summary of one run of samples.
    #[must_use]
    pub fn of(values: &[f64]) -> Self {
        let mut summary = Self::new();
        for value in values {
            summary.push(*value);
        }
        summary
    }

    /// Folds one more sample in.
    ///
    /// The running mean is what the squared deviation is measured against, both
    /// before and after this sample — Welford's update, which is stable where
    /// subtracting two large sums is not.
    pub fn push(&mut self, value: f64) {
        let before = self.mean().unwrap_or(0.0);

        self.count += 1;
        self.sum += value;
        self.sum_sq += value * value;
        self.lowest = self.lowest.min(value);
        self.highest = self.highest.max(value);

        let after = self.sum / self.count as f64;
        self.deviations += (value - before) * (value - after);
    }

    /// The summary of both runs together.
    ///
    /// Associative, and that is the contract: which order a tree of summaries
    /// is folded in cannot change the answer. Chan's parallel formula carries
    /// the squared deviations across, which needs both means and both counts
    /// and nothing else.
    #[must_use]
    pub fn merge(&self, other: &Self) -> Self {
        if self.count == 0 {
            return *other;
        }
        if other.count == 0 {
            return *self;
        }

        let count = self.count + other.count;
        let apart = (other.sum / other.count as f64) - (self.sum / self.count as f64);
        let weight = (self.count as f64 * other.count as f64) / count as f64;

        Self {
            count,
            sum: self.sum + other.sum,
            sum_sq: self.sum_sq + other.sum_sq,
            lowest: self.lowest.min(other.lowest),
            highest: self.highest.max(other.highest),
            deviations: self.deviations + other.deviations + apart * apart * weight,
        }
    }

    /// How many samples this summary folded.
    #[must_use]
    pub fn count(&self) -> u64 {
        self.count
    }

    /// Whether this summary folded no samples at all.
    #[must_use]
    pub fn is_empty(&self) -> bool {
        self.count == 0
    }

    /// Every sample added together.
    #[must_use]
    pub fn sum(&self) -> f64 {
        self.sum
    }

    /// Every sample squared and added together.
    #[must_use]
    pub fn sum_sq(&self) -> f64 {
        self.sum_sq
    }

    /// The lowest sample.
    #[must_use]
    pub fn min(&self) -> Option<f64> {
        (self.count > 0).then_some(self.lowest)
    }

    /// The highest sample.
    #[must_use]
    pub fn max(&self) -> Option<f64> {
        (self.count > 0).then_some(self.highest)
    }

    /// The distance between the highest sample and the lowest.
    #[must_use]
    pub fn span(&self) -> Option<f64> {
        (self.count > 0).then_some(self.highest - self.lowest)
    }

    /// The arithmetic mean.
    #[must_use]
    pub fn mean(&self) -> Option<f64> {
        (self.count > 0).then(|| self.sum / self.count as f64)
    }

    /// The population variance, read off the accumulated squared deviations.
    ///
    /// Population and not sample: a summary describes the samples it folded,
    /// which are the whole of what it knows.
    #[must_use]
    pub fn variance(&self) -> Option<f64> {
        (self.count > 0).then(|| self.deviations / self.count as f64)
    }

    /// The population standard deviation.
    #[must_use]
    pub fn deviation(&self) -> Option<f64> {
        self.variance().map(f64::sqrt)
    }

    /// The root mean square: `sqrt(sum_sq / count)`.
    #[must_use]
    pub fn rms(&self) -> Option<f64> {
        (self.count > 0).then(|| (self.sum_sq / self.count as f64).sqrt())
    }
}

#[cfg(test)]
mod tests {
    use std::f64::consts::TAU;

    use super::*;

    #[test]
    fn a_summary_of_nothing_answers_nothing() {
        let empty = Summary::new();

        assert!(empty.is_empty());
        assert_eq!(empty.count(), 0);
        assert_eq!(empty.sum(), 0.0);
        assert_eq!(empty.mean(), None);
        assert_eq!(empty.min(), None);
        assert_eq!(empty.max(), None);
        assert_eq!(empty.span(), None);
        assert_eq!(empty.variance(), None);
        assert_eq!(empty.deviation(), None);
        assert_eq!(empty.rms(), None);
    }

    #[test]
    fn a_summary_reports_the_run_it_folded() {
        let summary = Summary::of(&[3.0, 1.0, 4.0, 1.0, 5.0, 9.0, 2.0, 6.0]);

        assert_eq!(summary.count(), 8);
        assert_eq!(summary.sum(), 31.0);
        assert_eq!(summary.min(), Some(1.0));
        assert_eq!(summary.max(), Some(9.0));
        assert_eq!(summary.span(), Some(8.0));
        assert_eq!(summary.mean(), Some(31.0 / 8.0));
    }

    #[test]
    fn the_rms_of_a_unit_sine_is_one_over_root_two() {
        let samples: Vec<f64> = (0..10_000)
            .map(|index| (TAU * f64::from(index) / 10_000.0).sin())
            .collect();

        let rms = Summary::of(&samples)
            .rms()
            .expect("a full period has an rms");

        assert!((rms - 1.0 / 2.0_f64.sqrt()).abs() < 1e-4, "rms was {rms}");
    }

    #[test]
    fn a_constant_run_has_no_spread() {
        let summary = Summary::of(&[5.0; 100]);

        assert_eq!(summary.variance(), Some(0.0));
        assert_eq!(summary.deviation(), Some(0.0));
        assert_eq!(summary.rms(), Some(5.0));
    }

    #[test]
    fn variance_survives_a_mean_far_larger_than_the_spread() {
        // The reshape's reason for existing. Naive sum_sq/n − mean² over these
        // three is the difference of two numbers near 1e18: every digit of the
        // answer is lost, and the result is routinely negative.
        let summary = Summary::of(&[1e9, 1e9 + 1.0, 1e9 + 2.0]);

        let variance = summary.variance().expect("three samples have a variance");
        assert!(
            (variance - 2.0 / 3.0).abs() < 1e-9,
            "variance was {variance}"
        );

        let naive = summary.sum_sq() / 3.0 - (summary.sum() / 3.0).powi(2);
        assert!(
            naive < 0.0 || (naive - 2.0 / 3.0).abs() > 1e-3,
            "the naive form was accurate"
        );
    }

    #[test]
    fn a_split_run_summarizes_to_the_whole_run() {
        let samples: Vec<f64> = (0..64).map(f64::from).collect();

        let whole = Summary::of(&samples);
        let split = Summary::of(&samples[..17]).merge(&Summary::of(&samples[17..]));

        assert_eq!(whole.count(), split.count());
        assert_eq!(whole.sum(), split.sum());
        assert_eq!(whole.sum_sq(), split.sum_sq());
        assert_eq!(whole.min(), split.min());
        assert_eq!(whole.max(), split.max());
        assert!(
            (whole.variance().unwrap() - split.variance().unwrap()).abs() < 1e-12,
            "the split variance drifted"
        );
    }

    #[test]
    fn merging_is_associative() {
        // Exactly representable samples, so association is testable as equality
        // rather than as a tolerance.
        let (left, middle, right) = (
            Summary::of(&[1.0, 2.0, 3.0, 4.0]),
            Summary::of(&[8.0, 16.0]),
            Summary::of(&[32.0, 64.0, 128.0]),
        );

        let (one, other) = (
            left.merge(&middle).merge(&right),
            left.merge(&middle.merge(&right)),
        );

        assert_eq!(one.count(), other.count());
        assert_eq!(one.sum(), other.sum());
        assert_eq!(one.sum_sq(), other.sum_sq());
        assert_eq!(one.min(), other.min());
        assert_eq!(one.max(), other.max());
        // The squared deviations cross a merge through a division by the
        // combined count, so association holds to the arithmetic's precision
        // rather than to the bit.
        assert!(
            (one.variance().unwrap() - other.variance().unwrap()).abs() < 1e-9,
            "{one:?} against {other:?}"
        );
    }

    #[test]
    fn an_empty_summary_is_the_identity_of_a_merge() {
        let summary = Summary::of(&[1.0, 2.0, 3.0]);
        let empty = Summary::new();

        assert_eq!(summary.merge(&empty), summary);
        assert_eq!(empty.merge(&summary), summary);
        assert_eq!(empty.merge(&empty), empty);
    }

    #[test]
    fn pushing_one_sample_at_a_time_is_summarizing_the_run() {
        let samples = [2.5, -1.0, 7.25, 0.0];

        let mut pushed = Summary::new();
        for sample in samples {
            pushed.push(sample);
        }

        assert_eq!(pushed, Summary::of(&samples));
    }
}