diff --git a/.gitignore b/.gitignore index 9253a1c..6befc51 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ /vendor /img /.direnv +mutants.out/ diff --git a/Cargo.lock b/Cargo.lock index 60f6817..ce3cd8f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1063,7 +1063,7 @@ checksum = "e320a6c5ad31d271ad523dcf3ad13e2767ad8b1cb8f047f75a8aeaf8da139da2" [[package]] name = "sliding_features" -version = "7.1.3" +version = "8.0.0" dependencies = [ "criterion", "getset", diff --git a/Cargo.toml b/Cargo.toml index 38c4dce..d251c59 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -2,7 +2,7 @@ cargo-features = ["edition2024"] [package] name = "sliding_features" -version = "7.1.3" +version = "8.0.0" authors = ["MathisWellmann "] edition = "2024" license-file = "LICENSE" diff --git a/benches/vsct.rs b/benches/vsct.rs index 819f5a3..d66b239 100644 --- a/benches/vsct.rs +++ b/benches/vsct.rs @@ -15,7 +15,7 @@ use rand::{ use sliding_features::{ View, pure_functions::Echo, - sliding_windows::Vsct, + sliding_windows::ZScoreStandardization, }; fn criterion_benchmark(c: &mut Criterion) { @@ -26,7 +26,8 @@ fn criterion_benchmark(c: &mut Criterion) { group.bench_function("f64", |b| { let vals = Vec::::from_iter((0..N).map(|_| rng.random())); b.iter(|| { - let mut view = Vsct::::new(Echo::new(), NonZeroUsize::new(1024).unwrap()); + let mut view = + ZScoreStandardization::::new(Echo::new(), NonZeroUsize::new(1024).unwrap()); for v in vals.iter() { view.update(*v); let _ = black_box(view.last()); @@ -36,7 +37,8 @@ fn criterion_benchmark(c: &mut Criterion) { group.bench_function("f32", |b| { let vals = Vec::::from_iter((0..N).map(|_| rng.random())); b.iter(|| { - let mut view = Vsct::::new(Echo::new(), NonZeroUsize::new(1024).unwrap()); + let mut view = + ZScoreStandardization::::new(Echo::new(), NonZeroUsize::new(1024).unwrap()); for v in vals.iter() { view.update(*v); let _ = black_box(view.last()); diff --git a/examples/basic_chainable_view.rs b/examples/basic_chainable_view.rs index dafb341..573f257 100644 --- a/examples/basic_chainable_view.rs +++ b/examples/basic_chainable_view.rs @@ -9,7 +9,7 @@ use sliding_features::{ pure_functions::Echo, sliding_windows::{ Alma, - Vsct, + ZScoreStandardization, }, }; @@ -24,7 +24,7 @@ fn main() { let window_len = NonZeroUsize::new(20).unwrap(); let mut chain = Alma::new( // first, define the last function which gets applied in the chain - Vsct::new(Echo::new(), window_len), // Make the first transformation in the chain a VSCT + ZScoreStandardization::new(Echo::new(), window_len), // Make the first transformation in the chain a VSCT window_len, ); for v in &rands { diff --git a/flake.lock b/flake.lock index 85e7d68..61dbe22 100644 --- a/flake.lock +++ b/flake.lock @@ -20,11 +20,11 @@ }, "nixpkgs": { "locked": { - "lastModified": 1770169770, - "narHash": "sha256-awR8qIwJxJJiOmcEGgP2KUqYmHG4v/z8XpL9z8FnT1A=", + "lastModified": 1780246643, + "narHash": "sha256-4T1KWX7xWGQMs9hNZ24IOY3aYOi8D6+5WtkNRBSttB8=", "owner": "NixOS", "repo": "nixpkgs", - "rev": "aa290c9891fa4ebe88f8889e59633d20cc06a5f2", + "rev": "3109eaae18e09d0b8aef23dc2579e7d94b8d4b4e", "type": "github" }, "original": { @@ -50,11 +50,11 @@ }, "nixpks": { "locked": { - "lastModified": 1770197578, - "narHash": "sha256-AYqlWrX09+HvGs8zM6ebZ1pwUqjkfpnv8mewYwAo+iM=", + "lastModified": 1780243769, + "narHash": "sha256-x5UQuRsH3MqI0U9afaXSNqzTPSeZlRLvFAav2Ux1pNw=", "owner": "NixOS", "repo": "nixpkgs", - "rev": "00c21e4c93d963c50d4c0c89bfa84ed6e0694df2", + "rev": "331800de5053fcebacf6813adb5db9c9dca22a0c", "type": "github" }, "original": { @@ -77,11 +77,11 @@ "nixpkgs": "nixpkgs_2" }, "locked": { - "lastModified": 1770260791, - "narHash": "sha256-ADTBfENFjRVDQMcCycyX/pAy6NFI/Ct6Mrar3gsmXI0=", + "lastModified": 1780284119, + "narHash": "sha256-y2wR4Mk6D/N1ID4FZa2oUMStCUxyIoRzmgOOpLzoWmo=", "owner": "oxalica", "repo": "rust-overlay", - "rev": "42ec85352e419e601775c57256a52f6d48a39906", + "rev": "51390d0bfca0a68a8c337d215a4bbeddc2ca616e", "type": "github" }, "original": { diff --git a/flake.nix b/flake.nix index 4b22467..76664d7 100644 --- a/flake.nix +++ b/flake.nix @@ -57,6 +57,7 @@ taplo cargo-semver-checks cargo_upgrades + cargo-mutants ]; nix_tools = with pkgs; [ alejandra # Nix code formatter diff --git a/img/trend_flex.png b/img/trend_flex.png index 4c6d679..7ead302 100644 Binary files a/img/trend_flex.png and b/img/trend_flex.png differ diff --git a/img/vsct.png b/img/vsct.png index a24ccd9..e71ad49 100644 Binary files a/img/vsct.png and b/img/vsct.png differ diff --git a/img/welford_online_sliding.png b/img/welford_online_sliding.png index 1de445a..0f2a495 100644 Binary files a/img/welford_online_sliding.png and b/img/welford_online_sliding.png differ diff --git a/src/sliding_windows/mod.rs b/src/sliding_windows/mod.rs index ba217e4..e916962 100644 --- a/src/sliding_windows/mod.rs +++ b/src/sliding_windows/mod.rs @@ -26,8 +26,8 @@ mod sma; mod super_smoother; mod trend_flex; mod variance_stabilizing_transformation; -mod vsct; mod welford_online; +mod z_score_standardization; pub use alma::Alma; pub use binary_entropy::BinaryEntropy; @@ -54,5 +54,5 @@ pub use sma::Sma; pub use super_smoother::SuperSmoother; pub use trend_flex::TrendFlex; pub use variance_stabilizing_transformation::Vst; -pub use vsct::Vsct; pub use welford_online::WelfordOnline; +pub use z_score_standardization::ZScoreStandardization; diff --git a/src/sliding_windows/vsct.rs b/src/sliding_windows/vsct.rs deleted file mode 100644 index bc8e7cc..0000000 --- a/src/sliding_windows/vsct.rs +++ /dev/null @@ -1,92 +0,0 @@ -//! Variance Stabilizing Centering Transform Sliding Window - -use std::num::NonZeroUsize; - -use num::Float; - -use super::WelfordOnline; -use crate::{ - View, - pure_functions::Echo, -}; - -/// Variance Stabilizing Centering Transform Sliding Window -#[derive(Debug, Clone)] -pub struct Vsct { - view: V, - welford_online: WelfordOnline>, - last: T, -} - -impl Vsct -where - V: View, - T: Float, -{ - /// Create a new Variance Stabilizing Centering Transform with a chained View - /// and a given sliding window length - #[inline] - pub fn new(view: V, window_len: NonZeroUsize) -> Self { - Vsct { - view, - welford_online: WelfordOnline::new(Echo::new(), window_len), - last: T::zero(), - } - } - - /// The sliding window length. - #[inline(always)] - pub fn window_len(&self) -> NonZeroUsize { - self.welford_online.window_len() - } -} - -impl View for Vsct -where - V: View, - T: Float, -{ - fn update(&mut self, val: T) { - debug_assert!(val.is_finite(), "value must be finite"); - self.view.update(val); - let Some(val) = self.view.last() else { return }; - debug_assert!(val.is_finite(), "value must be finite"); - - self.welford_online.update(val); - self.last = val; - } - - fn last(&self) -> Option { - let std_dev = self.welford_online.last()?; - if std_dev == T::zero() { - return Some(T::zero()); - } - let mean = self.welford_online.mean(); - let out = (self.last - mean) / std_dev; - debug_assert!(out.is_finite(), "value must be finite"); - Some(out) - } -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::{ - plot::plot_values, - test_data::TEST_DATA, - }; - - #[test] - fn vsct_plot() { - let mut vsct = Vsct::new(Echo::new(), NonZeroUsize::new(16).unwrap()); - let mut out: Vec = Vec::with_capacity(TEST_DATA.len()); - for v in &TEST_DATA { - vsct.update(*v); - if let Some(val) = vsct.last() { - out.push(val); - } - } - let filename = "img/vsct.png"; - plot_values(out, filename).unwrap(); - } -} diff --git a/src/sliding_windows/welford_online.rs b/src/sliding_windows/welford_online.rs index 2a8cf3a..7f3e52e 100644 --- a/src/sliding_windows/welford_online.rs +++ b/src/sliding_windows/welford_online.rs @@ -54,10 +54,29 @@ where #[inline] fn update_stats_remove(&mut self, old_value: T) { - let delta = old_value - self.mean; - self.mean = self.mean - (delta / T::from(self.count).unwrap()); - self.m2 = self.m2 - (delta * (old_value - self.mean)); + // Derive sum and sum-of-squares from Welford state: + // mean = sum / n, m2 = sum_sq - n * mean^2 + // => sum = n * mean, sum_sq = m2 + n * mean^2 + let n = T::from(self.count).unwrap(); + let sum = self.mean * n; + let sum_sq = self.m2 + sum * self.mean; + + let sum_new = sum - old_value; + let sum_sq_new = sum_sq - old_value * old_value; self.count -= 1; + + if self.count > 0 { + let n_new = T::from(self.count).unwrap(); + self.mean = sum_new / n_new; + self.m2 = sum_sq_new - n_new * self.mean * self.mean; + // Clamp floating-point noise that could push m2 slightly negative. + if self.m2 < T::zero() { + self.m2 = T::zero(); + } + } else { + self.mean = T::zero(); + self.m2 = T::zero(); + } } /// Return the variance of the sliding window @@ -139,6 +158,47 @@ mod tests { assert_eq!(round(w_std_dev, 4), round(std_dev, 4)); } + #[test] + fn welford_online_sliding_matches_direct_computation() { + // Feed [1, 2, 3, 4, 5, 6] through a window of length 3. + // After each step, verify the std dev matches the direct formula for + // the window contents. + let mut wo = WelfordOnline::new(Echo::new(), NonZeroUsize::new(3).unwrap()); + let all: Vec = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; + + // Direct std-dev helper (sample sd). + let direct_sd = |vs: &[f64]| -> f64 { + let n = vs.len() as f64; + let m = vs.iter().sum::() / n; + let var = vs.iter().map(|v| (v - m).powi(2)).sum::() / (n - 1.0); + var.sqrt() + }; + + for i in 0..all.len() { + wo.update(all[i]); + let start = if i < 2 { 0 } else { i - 2 }; // i+1-window + let window = &all[start..=i]; + let expected = direct_sd(&window); + if let Some(got) = wo.last() { + let diff = (got - expected).abs(); + if diff > 1e-12 { + panic!( + "step {} window {:?}: expected sd={}, got={}, diff={}", + i, window, expected, got, diff + ); + } + } else { + // window not full yet → None is correct for i < 2 + assert!( + i < 2, + "step {}: got None but window is full ({:?})", + i, + window + ); + } + } + } + #[test] fn welford_online_plot() { let mut wo = WelfordOnline::new(Echo::new(), NonZeroUsize::new(16).unwrap()); diff --git a/src/sliding_windows/z_score_standardization.rs b/src/sliding_windows/z_score_standardization.rs new file mode 100644 index 0000000..6e9e310 --- /dev/null +++ b/src/sliding_windows/z_score_standardization.rs @@ -0,0 +1,187 @@ +//! Z-Score Standardization over a sliding window. +//! +//! Computes `(x - μ) / σ` where μ and σ are the sample mean and +//! standard deviation of the previous sliding window. + +use std::num::NonZeroUsize; + +use num::Float; + +use super::WelfordOnline; +use crate::{ + View, + pure_functions::Echo, +}; + +/// Z-Score Standardization over a sliding window. +/// +/// Computes `(x - μ) / σ` where μ and σ are the sample mean and +/// sample standard deviation derived from the previous sliding window. +/// +/// The current value is intentionally excluded from the statistics used to +/// standardize it, so the transformation only depends on past values and does +/// not leak information from the current/future sample into the rolling mean or +/// standard deviation. +#[derive(Debug, Clone)] +pub struct ZScoreStandardization { + view: V, + welford_online: WelfordOnline>, + history_len: usize, + out: Option, +} + +impl ZScoreStandardization +where + V: View, + T: Float, +{ + /// Create a new ZScoreStandardization with a chained View + /// and a given sliding window length + #[inline] + pub fn new(view: V, window_len: NonZeroUsize) -> Self { + ZScoreStandardization { + view, + welford_online: WelfordOnline::new(Echo::new(), window_len), + history_len: 0, + out: None, + } + } + + /// The sliding window length. + #[inline(always)] + pub fn window_len(&self) -> NonZeroUsize { + self.welford_online.window_len() + } +} + +impl View for ZScoreStandardization +where + V: View, + T: Float, +{ + fn update(&mut self, val: T) { + debug_assert!(val.is_finite(), "value must be finite"); + self.view.update(val); + let Some(val) = self.view.last() else { return }; + debug_assert!(val.is_finite(), "value must be finite"); + + self.out = if self.history_len >= self.window_len().get() { + let std_dev = self.welford_online.last().expect("history is warm"); + if std_dev == T::zero() { + Some(T::zero()) + } else { + let mean = self.welford_online.mean(); + let out = (val - mean) / std_dev; + debug_assert!(out.is_finite(), "value must be finite"); + Some(out) + } + } else { + None + }; + + self.welford_online.update(val); + self.history_len = (self.history_len + 1).min(self.window_len().get()); + } + + fn last(&self) -> Option { + self.out + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{ + plot::plot_values, + test_data::TEST_DATA, + }; + + #[test] + fn z_score_plot() { + let mut vsct = ZScoreStandardization::new(Echo::new(), NonZeroUsize::new(16).unwrap()); + let mut out: Vec = Vec::with_capacity(TEST_DATA.len()); + for v in &TEST_DATA { + vsct.update(*v); + if let Some(val) = vsct.last() { + out.push(val); + } + } + let filename = "img/vsct.png"; + plot_values(out, filename).unwrap(); + } + + #[test] + fn z_score_matches_direct_past_window_computation() { + // Feed [1, 2, 3, 4, 5, 6] through a z-score normalizer with window=3. + // After warm-up, verify the z-score of the current value equals the + // direct formula: (x - μ) / σ where μ,σ are sample stats of the three + // previous values. The current value must not be part of the fitted + // normalization window. + let mut zs = ZScoreStandardization::new(Echo::new(), NonZeroUsize::new(3).unwrap()); + let all: Vec = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; + + let direct_z_score = |x: f64, history: &[f64]| -> f64 { + let n = history.len() as f64; + let m = history.iter().sum::() / n; + let var = history.iter().map(|v| (v - m).powi(2)).sum::() / (n - 1.0); + (x - m) / var.sqrt() + }; + + for i in 0..all.len() { + zs.update(all[i]); + if i >= 3 { + let history = &all[i - 3..i]; + let expected = direct_z_score(all[i], history); + let got = zs.last().unwrap(); + let diff = (got - expected).abs(); + assert!( + diff < 1e-12, + "step {} history {:?}: expected z={}, got={}, diff={}", + i, + history, + expected, + got, + diff + ); + } else { + assert!( + zs.last().is_none(), + "step {}: got {:?} before {} past samples were available", + i, + zs.last(), + zs.window_len() + ); + } + } + } + + #[test] + fn z_score_does_not_leak_current_value_into_statistics() { + let mut zs = ZScoreStandardization::new(Echo::new(), NonZeroUsize::new(3).unwrap()); + for val in [1.0, 2.0, 3.0] { + zs.update(val); + } + + zs.update(100.0); + + // If the current value leaked into the statistics, this would be the + // z-score against [2, 3, 100] (roughly 1.1545). The causal value uses + // only the previous window [1, 2, 3]. + let expected = (100.0 - 2.0) / 1.0; + let got = zs.last().unwrap(); + assert!( + (got - expected).abs() < 1e-12, + "expected {expected}, got {got}" + ); + } + + #[test] + fn z_score_of_identical_values_is_zero() { + // When all values in the window are the same, std_dev = 0 → z = 0. + let mut zs = ZScoreStandardization::new(Echo::new(), NonZeroUsize::new(4).unwrap()); + for _ in 0..10 { + zs.update(5.0); + } + assert_eq!(zs.last(), Some(0.0)); + } +}