From 03f6f27c4bbf62885f5e8b5107150bb557f6ec85 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 18:16:12 +0000 Subject: [PATCH 1/9] hdr: harvest the rolling floor from lance-graph onto exact moments Adds hpc::rolling_floor, the adaptive half of the HDR exposure meter, harvested from lance-graph graph/blasgraph/hdr.rs (the behavioural reference). Constants, cadence, drift rule and reset semantics carry over unchanged; the arithmetic underneath is #327's exact MomentsU32 instead of the reference's approximate integer Welford. - isqrt_u32: the reference integer Newton square root. - ReservoirU32: deterministic Algorithm-R reservoir (splitmix replacement hash keyed on the observation count), u32 empirical quantile, Pearson second skewness, kurtosis x100. Order-defined; no merge law. - RollingFloor: calibrated (mu, sigma), sigma floors mu - k sigma, empirical floors at the reference percentiles, shape evaluation every 1000 observations after the first 1000 (reservoir >= 100), normal iff |skew| < 2 and 200 < kurt < 500, drift at |dmu| > sigma/2 or |dsigma| > sigma/4, recalibration resets moments, reservoir and shape. - observe_batch folds moments_u32 between checkpoints and stops after the first shift, so any batching reproduces the scalar observe/recalibrate loop exactly. - MomentsU32::observe: scalar fold equal to merging a singleton batch. A test keeps the legacy integer Welford as an oracle: across 95 checkpoints on five streams the (mu, sigma) the drift rule sees are identical. examples/hdr_rolling_floor_bench.rs separates the hot path from the periodic shape path. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- examples/hdr_rolling_floor_bench.rs | 124 ++++ src/hpc/mod.rs | 1 + src/hpc/rolling_floor.rs | 911 ++++++++++++++++++++++++++++ src/hpc/statistics.rs | 10 + 4 files changed, 1046 insertions(+) create mode 100644 examples/hdr_rolling_floor_bench.rs create mode 100644 src/hpc/rolling_floor.rs diff --git a/examples/hdr_rolling_floor_bench.rs b/examples/hdr_rolling_floor_bench.rs new file mode 100644 index 00000000..c44da012 --- /dev/null +++ b/examples/hdr_rolling_floor_bench.rs @@ -0,0 +1,124 @@ +//! Cost split of the HDR rolling floor: the per-observation hot path versus +//! the periodic shape path. +//! +//! ```sh +//! cargo run --release --example hdr_rolling_floor_bench +//! ``` +//! +//! Hot path, per observation: +//! 1. popcount / Hamming distance of two 2048-byte vectors +//! 2. exact moments update (`MomentsU32::observe`, and `moments_u32` batch) +//! 3. reservoir update +//! 4. full rolling-floor update (moments + reservoir + checkpoint test) +//! +//! Periodic path, once per 1000 observations: +//! 5. shape evaluation (sort 1000 samples, median, kurtosis) +//! 6. empirical quantiles (12 lookups into the sorted reservoir) + +use ndarray::hpc::bitwise::hamming_distance_raw; +use ndarray::hpc::rolling_floor::{quantile_of_sorted, ReservoirU32, RollingFloor}; +use ndarray::hpc::statistics::{moments_u32, MomentsU32}; +use std::hint::black_box; +use std::time::Instant; + +fn xorshift(n: usize, mut s: u64) -> Vec { + (0..n) + .map(|_| { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + s + }) + .collect() +} + +fn ns_per(label: &str, n: usize, f: impl FnOnce()) { + let t = Instant::now(); + f(); + let dt = t.elapsed(); + println!("{label:<44} {:>9.2} ns/op ({n} ops)", dt.as_nanos() as f64 / n as f64); +} + +fn main() { + println!("avx512f={} avx2={}", cfg!(target_feature = "avx512f"), cfg!(target_feature = "avx2")); + + const VBYTES: usize = 2048; // 16384-bit vectors + const N: usize = 2_000_000; + let words = xorshift(VBYTES / 8 * 65, 1); + let bytes: Vec = words.iter().flat_map(|w| w.to_le_bytes()).collect(); + let query = &bytes[..VBYTES]; + let db = &bytes[VBYTES..]; + let dists: Vec = (0..N) + .map(|i| { + let j = i % 64; + hamming_distance_raw(query, &db[j * VBYTES..(j + 1) * VBYTES]) as u32 + }) + .collect(); + + println!("-- hot path, per observation --"); + ns_per("1. popcount/Hamming (2048 B)", N, || { + let mut acc = 0u64; + for i in 0..N { + let j = i % 64; + acc += hamming_distance_raw(black_box(query), &db[j * VBYTES..(j + 1) * VBYTES]); + } + black_box(acc); + }); + ns_per("2a. MomentsU32::observe (scalar)", N, || { + let mut m = MomentsU32::default(); + for &d in &dists { + m.observe(black_box(d)); + } + black_box(m); + }); + ns_per("2b. moments_u32 (batch)", N, || { + black_box(moments_u32(black_box(&dists))); + }); + ns_per("3. ReservoirU32::observe (cap 1000)", N, || { + let mut r = ReservoirU32::new(1000); + for &d in &dists { + r.observe(black_box(d)); + } + black_box(r.len()); + }); + ns_per("4a. RollingFloor::observe (incl. checkpoints)", N, || { + let mut f = RollingFloor::for_width(16384); + for &d in &dists { + if let Some(s) = f.observe(black_box(d)) { + f.recalibrate(&s); + } + } + black_box(f.mu()); + }); + ns_per("4b. RollingFloor::observe_batch (incl. checkpoints)", N, || { + let mut f = RollingFloor::for_width(16384); + let mut rest: &[u32] = &dists; + while !rest.is_empty() { + let (used, s) = f.observe_batch(rest); + rest = &rest[used..]; + if let Some(s) = s { + f.recalibrate(&s); + } + } + black_box(f.mu()); + }); + + println!("-- periodic path, per checkpoint (every 1000 observations) --"); + let mut r = ReservoirU32::new(1000); + dists[..5000].iter().for_each(|&d| r.observe(d)); + const K: usize = 20_000; + ns_per("5. shape: sort + median + kurtosis", K, || { + for _ in 0..K { + let sorted = black_box(&r).sorted(); + black_box(quantile_of_sorted(&sorted, 0.5)); + black_box(r.kurtosis(8192, 64)); + } + }); + let sorted = r.sorted(); + ns_per("6. empirical floors: 12 quantiles (sorted)", K, || { + for _ in 0..K { + black_box(RollingFloor::FLOOR_PERCENTILES.map(|p| quantile_of_sorted(black_box(&sorted), p))); + black_box(RollingFloor::CASCADE_PERCENTILES.map(|p| quantile_of_sorted(black_box(&sorted), p))); + } + }); +} diff --git a/src/hpc/mod.rs b/src/hpc/mod.rs index 074921b6..92cd345f 100644 --- a/src/hpc/mod.rs +++ b/src/hpc/mod.rs @@ -25,6 +25,7 @@ pub mod blas_level2; pub mod blas_level3; pub mod reductions; pub mod statistics; +pub mod rolling_floor; /// Reliability & validity statistics: Pearson r, Spearman ρ, Cronbach α, ICC. pub mod reliability; /// Entropy ladder: Staunen↔Wisdom coordinate over NARS truth + Pearl-2³ SPO. diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs new file mode 100644 index 00000000..ca94599f --- /dev/null +++ b/src/hpc/rolling_floor.rs @@ -0,0 +1,911 @@ +//! HDR rolling floor: an online distribution floor for popcount / Hamming +//! observations. +//! +//! This is the adaptive half of the HDR exposure meter, harvested from +//! lance-graph's `graph/blasgraph/hdr.rs`. That file stays the behavioural +//! reference: constants, cadence, drift rule and reset semantics are carried +//! over unchanged. What changed is the arithmetic underneath. The old +//! approximate integer Welford is replaced by the exact, mergeable +//! [`MomentsU32`]. +//! +//! # Two layers, two cadences +//! +//! * **Parameters (continuous, cheap).** Every observation folds into exact +//! `(n, Σx, Σx²)`. Location and spread are available at any time, so the +//! floor can be used from the first observation and only gets better as the +//! population grows. There is no training barrier. +//! * **Shape (periodic, amortised).** A deterministic Algorithm-R reservoir +//! feeds median, empirical quantiles, skewness and kurtosis. They are +//! evaluated once every [`RollingFloor::EVAL_CADENCE`] observations, never on +//! the per-observation path. +//! +//! # Floors +//! +//! While the shape reads as normal, the floors are the analytical quantiles of +//! the calibrated `(μ, σ)`: `μ − kσ`, i.e. `μ + σ·Φ⁻¹(p)` at the percentiles +//! the empirical tables name. When the shape does not read as normal, the +//! floors are the reservoir's empirical quantiles at those same percentiles. +//! +//! # What is and is not mergeable +//! +//! The moments are exact and partition-independent: shards can be merged in +//! any order. The reservoir is not. It is deterministic for a given +//! observation *order* and has no merge law. [`RollingFloor::observe_batch`] +//! therefore reproduces the scalar stream exactly rather than merging shards. +//! +//! # Reference behaviour this preserves, including its limits +//! +//! * Floors move only when a drift alert is acted on +//! ([`RollingFloor::recalibrate`]), not continuously with the running +//! parameters. +//! * A parameter drift resets the shape layer too (reservoir, empirical mode, +//! skewness, kurtosis). The reference does not separate parameter drift from +//! shape drift; neither does this port. + +use super::statistics::{moments_u32, MomentsU32}; + +/// Floor of `√n`, integer Newton iteration. Exact for every `u32`. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::isqrt_u32; +/// assert_eq!(isqrt_u32(4095), 63); +/// assert_eq!(isqrt_u32(4096), 64); +/// assert_eq!(isqrt_u32(u32::MAX), 65535); +/// ``` +pub fn isqrt_u32(n: u32) -> u32 { + if n == 0 { + return 0; + } + // The start value must be >= floor(√n) so Newton descends monotonically. + let mut x = 1u32 << ((33 - n.leading_zeros()) / 2); + loop { + let x1 = (x + n / x) / 2; + if x1 >= x { + return x; + } + x = x1; + } +} + +/// Deterministic reservoir sample of a `u32` stream (Vitter's Algorithm R). +/// +/// Every element seen so far has the same chance of being held. The +/// replacement decision for the `k`-th element is a fixed hash of `k`, so an +/// identical stream always produces an identical reservoir. The result depends +/// on the order of the stream; there is no merge of two reservoirs. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::ReservoirU32; +/// let mut r = ReservoirU32::new(4); +/// for d in [10, 20, 30] { +/// r.observe(d); +/// } +/// assert_eq!(r.len(), 3); +/// assert_eq!(r.quantile(0.5), 20); +/// ``` +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct ReservoirU32 { + samples: Vec, + capacity: usize, + seen: u64, +} + +impl ReservoirU32 { + /// An empty reservoir holding at most `capacity` samples. + pub fn new(capacity: usize) -> Self { + Self { + samples: Vec::with_capacity(capacity), + capacity, + seen: 0, + } + } + + /// Offer one value to the reservoir. + #[inline] + pub fn observe(&mut self, value: u32) { + self.seen += 1; + if self.samples.len() < self.capacity { + self.samples.push(value); + } else { + // Replace slot j with probability capacity / seen. + let j = Self::replacement_hash(self.seen) % self.seen; + if (j as usize) < self.capacity { + self.samples[j as usize] = value; + } + } + } + + /// Maximum number of samples held. + pub fn capacity(&self) -> usize { + self.capacity + } + + /// Number of values offered so far. + pub fn seen(&self) -> u64 { + self.seen + } + + /// Number of samples currently held. + pub fn len(&self) -> usize { + self.samples.len() + } + + /// `true` when no value has been offered yet. + pub fn is_empty(&self) -> bool { + self.samples.is_empty() + } + + /// The held samples, in slot order. + pub fn samples(&self) -> &[u32] { + &self.samples + } + + /// The held samples, sorted ascending. Sort once and use + /// [`quantile_of_sorted`] when several quantiles are needed. + pub fn sorted(&self) -> Vec { + let mut s = self.samples.clone(); + s.sort_unstable(); + s + } + + /// Empirical quantile: the sorted sample at index `⌊q·len⌋`, clamped to + /// the last element. `0` for an empty reservoir. + /// + /// O(len · log len); meant for the periodic shape path, not per + /// observation. + pub fn quantile(&self, q: f32) -> u32 { + quantile_of_sorted(&self.sorted(), q) + } + + /// Pearson's second skewness `3(μ − median) / σ`, in integer division. + /// Positive means right-skewed, `0` symmetric. `0` when `σ = 0` or the + /// reservoir is empty. + pub fn skewness(&self, mu: u32, sigma: u32) -> i32 { + if sigma == 0 || self.samples.is_empty() { + return 0; + } + skewness_from_median(mu, sigma, self.quantile(0.5)) + } + + /// Kurtosis ×100: `100 · E[(X − μ)⁴] / σ⁴` over the reservoir, with the + /// normal distribution at 300. Returns 300 when `σ = 0` or fewer than 4 + /// samples are held. + pub fn kurtosis(&self, mu: u32, sigma: u32) -> u32 { + if sigma == 0 || self.samples.len() < 4 { + return 300; + } + let n = self.samples.len() as u128; + // u128: a fourth power of a u32 difference is below 2^128. + let m4: u128 = self + .samples + .iter() + .map(|&d| { + let diff = u128::from(d.abs_diff(mu)); + diff * diff * diff * diff + }) + .sum::() + / n; + let s4 = u128::from(sigma).pow(4); + u32::try_from(m4 * 100 / s4).unwrap_or(u32::MAX) + } + + /// Deterministic splitmix64-style hash that drives replacement. + fn replacement_hash(seed: u64) -> u64 { + let mut z = seed.wrapping_add(0x9e3779b97f4a7c15); + z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9); + z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb); + z ^ (z >> 31) + } +} + +/// Empirical quantile of an ascending slice: the element at `⌊q·len⌋`, +/// clamped to the last element. `0` for an empty slice. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::quantile_of_sorted; +/// let s = [1, 2, 3, 4]; +/// assert_eq!(quantile_of_sorted(&s, 0.0), 1); +/// assert_eq!(quantile_of_sorted(&s, 0.5), 3); +/// assert_eq!(quantile_of_sorted(&s, 1.0), 4); +/// ``` +pub fn quantile_of_sorted(sorted: &[u32], q: f32) -> u32 { + if sorted.is_empty() { + return 0; + } + let idx = ((q * sorted.len() as f32) as usize).min(sorted.len() - 1); + sorted[idx] +} + +fn skewness_from_median(mu: u32, sigma: u32, median: u32) -> i32 { + // i64 so no difference of two u32 values can overflow. + let s = 3 * (i64::from(mu) - i64::from(median)) / i64::from(sigma); + s.clamp(i64::from(i32::MIN), i64::from(i32::MAX)) as i32 +} + +/// A detected drift of the running parameters away from the calibrated ones. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct FloorShift { + /// Calibrated mean before the shift. + pub old_mu: u32, + /// Running mean at the checkpoint that raised the shift. + pub new_mu: u32, + /// Calibrated standard deviation before the shift. + pub old_sigma: u32, + /// Running standard deviation at that checkpoint. + pub new_sigma: u32, + /// Observation count at that checkpoint. + pub observations: u64, +} + +/// The HDR rolling floor for a stream of `u32` distances. +/// +/// It holds the calibrated `(μ, σ)`, the sigma floors derived from them, and +/// empirical floors from a reservoir. Running parameters come from exact +/// [`MomentsU32`]. Every [`EVAL_CADENCE`](Self::EVAL_CADENCE) observations +/// (after the first cadence) it re-reads the distribution shape, selects sigma +/// or empirical floors, and checks for drift. +/// +/// The floors are lower-is-better thresholds, from strictest to loosest: +/// four band floors at `[μ−3σ, μ−2σ, μ−σ, μ]` and eight cascade floors at +/// quarter-σ steps from `μ−σ` to `μ−3σ`. How a caller turns them into bands +/// is its own policy. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::RollingFloor; +/// let mut floor = RollingFloor::for_width(16384); +/// assert_eq!(floor.mu(), 8192); +/// assert_eq!(floor.sigma(), 64); +/// assert_eq!(floor.active_floors(), [8000, 8064, 8128, 8192]); +/// // Usable immediately; observations refine it as they arrive. +/// for d in 0..3000u32 { +/// if let Some(shift) = floor.observe(8192 + (d % 7)) { +/// floor.recalibrate(&shift); +/// } +/// } +/// ``` +#[derive(Debug, Clone)] +pub struct RollingFloor { + mu: u32, + sigma: u32, + sigma_floors: [u32; 4], + sigma_cascade: [u32; 8], + reservoir: ReservoirU32, + empirical_floors: [u32; 4], + empirical_cascade: [u32; 8], + use_empirical: bool, + skewness: i32, + kurtosis: u32, + moments: MomentsU32, +} + +impl RollingFloor { + /// Reservoir capacity. + pub const RESERVOIR_CAP: usize = 1000; + /// Shape and drift are evaluated when the observation count is a multiple + /// of this, and larger than it. + pub const EVAL_CADENCE: u64 = 1000; + /// Minimum reservoir population before shape is evaluated. + pub const MIN_SHAPE_SAMPLES: usize = 100; + /// Kurtosis ×100 of the normal distribution. + pub const NORMAL_KURTOSIS: u32 = 300; + /// Percentiles of the empirical band floors: ≈3σ, 2σ, 1σ, median. + pub const FLOOR_PERCENTILES: [f32; 4] = [0.001, 0.023, 0.159, 0.500]; + /// Percentiles of the empirical cascade floors: 1σ, 1.5σ, 1.75σ, 2σ, + /// 2.25σ, 2.5σ, 2.75σ, 3σ below the mean. + pub const CASCADE_PERCENTILES: [f32; 8] = [0.1587, 0.0668, 0.0401, 0.0228, 0.0122, 0.0062, 0.0030, 0.0013]; + + /// Band floors `[μ−3σ, μ−2σ, μ−σ, μ]`, saturating at zero. + pub fn sigma_floors_of(mu: u32, sigma: u32) -> [u32; 4] { + [mu.saturating_sub(3 * sigma), mu.saturating_sub(2 * sigma), mu.saturating_sub(sigma), mu] + } + + /// Cascade floors at `μ − kσ/4` for `k = 4, 6, 7, 8, 9, 10, 11, 12`, + /// saturating at zero. + pub fn sigma_cascade_of(mu: u32, sigma: u32) -> [u32; 8] { + [4, 6, 7, 8, 9, 10, 11, 12].map(|k: u32| mu.saturating_sub(k * sigma / 4)) + } + + /// Floors that assume only a prior `(μ, σ)`, with no observations yet. + pub fn from_params(mu: u32, sigma: u32) -> Self { + let sigma_floors = Self::sigma_floors_of(mu, sigma); + let sigma_cascade = Self::sigma_cascade_of(mu, sigma); + Self { + mu, + sigma, + sigma_floors, + sigma_cascade, + reservoir: ReservoirU32::new(Self::RESERVOIR_CAP), + empirical_floors: sigma_floors, + empirical_cascade: sigma_cascade, + use_empirical: false, + skewness: 0, + kurtosis: Self::NORMAL_KURTOSIS, + moments: MomentsU32::default(), + } + } + + /// Resume from calibrated `(μ, σ)` and already-accumulated running + /// moments, with an empty reservoir. The next checkpoint follows from + /// `moments.n`. + pub fn from_params_and_moments(mu: u32, sigma: u32, moments: MomentsU32) -> Self { + let mut floor = Self::from_params(mu, sigma); + floor.moments = moments; + floor + } + + /// Binomial prior for the Hamming distance of two random + /// `total_bits`-bit vectors: `μ = bits/2`, `σ = max(1, ⌊√(bits/4)⌋)`. + pub fn for_width(total_bits: u32) -> Self { + Self::from_params(total_bits / 2, isqrt_u32(total_bits / 4).max(1)) + } + + /// Calibrate from a warm-up sample of at least two distances. + /// + /// `μ` is the floor of the sample mean and `σ` is + /// `max(1, ⌊√⌊Σ(x − μ)² / n⌋⌋)`, spread measured around that integer + /// mean, exactly as the reference does. The sample seeds the reservoir and + /// the running moments; the floors start in sigma mode. + /// + /// # Panics + /// + /// If `sample` has fewer than two values. + pub fn calibrate(sample: &[u32]) -> Self { + assert!(sample.len() > 1, "need at least 2 samples to calibrate"); + let moments = moments_u32(sample); + let mu = saturate_u32(moments.sum / u128::from(moments.n)); + let sigma = isqrt_u32(saturate_u32(centred_on_floor_mean(&moments) / u128::from(moments.n))).max(1); + let mut floor = Self::from_params(mu, sigma); + for &d in sample { + floor.reservoir.observe(d); + } + let sorted = floor.reservoir.sorted(); + floor.empirical_floors = Self::FLOOR_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); + floor.empirical_cascade = Self::CASCADE_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); + floor.moments = moments; + floor + } + + /// Fold one observation in. Returns a shift when this observation lands + /// on a checkpoint and the running parameters have drifted: + /// `|μ_run − μ| > σ/2` or `|σ_run − σ| > σ/4`, against the calibrated + /// `(μ, σ)`. + #[inline] + pub fn observe(&mut self, distance: u32) -> Option { + self.moments.observe(distance); + self.reservoir.observe(distance); + if self.at_checkpoint() { + self.evaluate() + } else { + None + } + } + + /// Fold a batch in, stopping right after the first checkpoint that + /// raises a shift. + /// + /// Returns how many values were consumed and the shift, if any. Feeding + /// the unconsumed rest in after acting on the shift reproduces the scalar + /// loop `if let Some(s) = observe(d) { recalibrate(&s) }` exactly, for any + /// batching: the moments between checkpoints go through + /// [`moments_u32`], the reservoir sees every value in stream order, and + /// every checkpoint is evaluated on the same state the scalar loop sees. + pub fn observe_batch(&mut self, distances: &[u32]) -> (usize, Option) { + let mut consumed = 0; + while consumed < distances.len() { + let to_checkpoint = (Self::EVAL_CADENCE - self.moments.n % Self::EVAL_CADENCE) as usize; + let take = to_checkpoint.min(distances.len() - consumed); + let chunk = &distances[consumed..consumed + take]; + self.moments = self.moments.merge(moments_u32(chunk)); + for &d in chunk { + self.reservoir.observe(d); + } + consumed += take; + if take == to_checkpoint && self.at_checkpoint() { + if let Some(shift) = self.evaluate() { + return (consumed, Some(shift)); + } + } + } + (consumed, None) + } + + /// Adopt the shifted parameters as the new calibration and restart + /// observation from scratch: running moments, reservoir, empirical mode + /// and shape diagnostics are all reset. `σ` is floored at 1. + pub fn recalibrate(&mut self, shift: &FloorShift) { + let capacity = self.reservoir.capacity(); + *self = Self::from_params(shift.new_mu, shift.new_sigma.max(1)); + self.reservoir = ReservoirU32::new(capacity); + } + + fn at_checkpoint(&self) -> bool { + let n = self.moments.n; + n.is_multiple_of(Self::EVAL_CADENCE) && n > Self::EVAL_CADENCE + } + + /// The periodic path: shape evaluation, floor selection, drift check. + fn evaluate(&mut self) -> Option { + let run_mu = saturate_u32(self.moments.sum / u128::from(self.moments.n)); + let run_sigma = isqrt_u32(saturate_u32(variance_floor(&self.moments))).max(1); + + if self.reservoir.len() >= Self::MIN_SHAPE_SAMPLES { + let sorted = self.reservoir.sorted(); + self.skewness = skewness_from_median(run_mu, run_sigma, quantile_of_sorted(&sorted, 0.5)); + self.kurtosis = self.reservoir.kurtosis(run_mu, run_sigma); + if self.shape_is_normal() { + self.use_empirical = false; + } else { + self.empirical_floors = Self::FLOOR_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); + self.empirical_cascade = Self::CASCADE_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); + self.use_empirical = true; + } + } + + let mu_drift = run_mu.abs_diff(self.mu); + let sigma_drift = run_sigma.abs_diff(self.sigma); + if mu_drift > self.sigma / 2 || sigma_drift > self.sigma / 4 { + Some(FloorShift { + old_mu: self.mu, + new_mu: run_mu, + old_sigma: self.sigma, + new_sigma: run_sigma, + observations: self.moments.n, + }) + } else { + None + } + } + + /// The reference normality window: `|skew| < 2` and `200 < kurt < 500`. + pub fn shape_is_normal(&self) -> bool { + self.skewness.abs() < 2 && self.kurtosis > 200 && self.kurtosis < 500 + } + + /// Calibrated mean. + pub fn mu(&self) -> u32 { + self.mu + } + + /// Calibrated standard deviation. + pub fn sigma(&self) -> u32 { + self.sigma + } + + /// Band floors derived from the calibrated `(μ, σ)`. + pub fn sigma_floors(&self) -> [u32; 4] { + self.sigma_floors + } + + /// Cascade floors derived from the calibrated `(μ, σ)`. + pub fn sigma_cascade(&self) -> [u32; 8] { + self.sigma_cascade + } + + /// Band floors from the reservoir's empirical quantiles. + pub fn empirical_floors(&self) -> [u32; 4] { + self.empirical_floors + } + + /// Cascade floors from the reservoir's empirical quantiles. + pub fn empirical_cascade(&self) -> [u32; 8] { + self.empirical_cascade + } + + /// The band floors in effect: empirical when the shape read as non-normal + /// at the last checkpoint, sigma otherwise. + pub fn active_floors(&self) -> [u32; 4] { + if self.use_empirical { + self.empirical_floors + } else { + self.sigma_floors + } + } + + /// The cascade floors in effect, chosen like [`active_floors`](Self::active_floors). + pub fn active_cascade(&self) -> [u32; 8] { + if self.use_empirical { + self.empirical_cascade + } else { + self.sigma_cascade + } + } + + /// Whether the empirical floors are in effect. + pub fn is_empirical(&self) -> bool { + self.use_empirical + } + + /// Skewness at the last shape evaluation (`0` before any). + pub fn skewness(&self) -> i32 { + self.skewness + } + + /// Kurtosis ×100 at the last shape evaluation (300 before any). + pub fn kurtosis(&self) -> u32 { + self.kurtosis + } + + /// Exact running moments since calibration or the last recalibration. + pub fn moments(&self) -> MomentsU32 { + self.moments + } + + /// Observation count since calibration or the last recalibration. + pub fn observations(&self) -> u64 { + self.moments.n + } + + /// The reservoir behind the empirical floors. + pub fn reservoir(&self) -> &ReservoirU32 { + &self.reservoir + } +} + +fn saturate_u32(x: u128) -> u32 { + u32::try_from(x).unwrap_or(u32::MAX) +} + +/// `Σ(x − q)²` with `q = ⌊Σx / n⌋`, exactly. With `r = Σx − n·q` this is +/// `Σx² − q·Σx − q·r`, and no intermediate goes negative. Requires `n > 0`. +fn centred_on_floor_mean(m: &MomentsU32) -> u128 { + let n = u128::from(m.n); + let q = m.sum / n; + let r = m.sum % n; + m.sum_sq - q * m.sum - q * r +} + +/// `⌊M2 / n⌋` for the exact population variance, with `M2 = Σ(x − mean)²`. +/// Requires `n > 0`. +fn variance_floor(m: &MomentsU32) -> u128 { + let n = u128::from(m.n); + if let (Some(a), Some(b)) = (n.checked_mul(m.sum_sq), m.sum.checked_mul(m.sum)) { + return (a - b) / (n * n); + } + // n·M2 = C·n − r², with C centred on the floor mean. Write C = a·n + b; + // then ⌊(C·n − r²)/n²⌋ = a − [b·n < r²], and b·n, r² < 2^128. + let c = centred_on_floor_mean(m); + let r = m.sum % n; + let (a, b) = (c / n, c % n); + if b * n >= r * r { + a + } else { + a - 1 + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn stream(n: usize, base: u32, spread: u32, mut s: u64) -> Vec { + (0..n) + .map(|_| { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + base + (s % u64::from(spread)) as u32 + }) + .collect() + } + + /// Roughly normal around `mu` with the given spread (sum of 12 uniforms). + fn normalish(n: usize, mu: u32, sigma: u32, mut s: u64) -> Vec { + (0..n) + .map(|_| { + let mut acc = 0i64; + for _ in 0..12 { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + acc += (s % 1000) as i64; + } + // Sum of 12 U(0,1000) has sd ≈ 1000. + (i64::from(mu) + (acc - 6000) * i64::from(sigma) / 1000).max(0) as u32 + }) + .collect() + } + + /// Scalar reference loop: observe, recalibrate on every shift. + fn run_scalar(mut f: RollingFloor, xs: &[u32]) -> (RollingFloor, Vec) { + let mut shifts = Vec::new(); + for &d in xs { + if let Some(s) = f.observe(d) { + f.recalibrate(&s); + shifts.push(s); + } + } + (f, shifts) + } + + fn run_batched(mut f: RollingFloor, xs: &[u32], chunk: usize) -> (RollingFloor, Vec) { + let mut shifts = Vec::new(); + for c in xs.chunks(chunk) { + let mut rest = c; + while !rest.is_empty() { + let (used, s) = f.observe_batch(rest); + rest = &rest[used..]; + if let Some(s) = s { + f.recalibrate(&s); + shifts.push(s); + } + } + } + (f, shifts) + } + + fn same_state(a: &RollingFloor, b: &RollingFloor) { + assert_eq!(a.moments(), b.moments()); + assert_eq!(a.reservoir(), b.reservoir()); + assert_eq!((a.mu(), a.sigma()), (b.mu(), b.sigma())); + assert_eq!(a.active_floors(), b.active_floors()); + assert_eq!(a.active_cascade(), b.active_cascade()); + assert_eq!((a.is_empirical(), a.skewness(), a.kurtosis()), (b.is_empirical(), b.skewness(), b.kurtosis())); + } + + #[test] + fn isqrt_is_floor_sqrt() { + for n in (0..200_000u32).chain([u32::MAX, u32::MAX - 1, 65535 * 65535, 65536 * 65535]) { + let r = isqrt_u32(n); + assert!(u64::from(r) * u64::from(r) <= u64::from(n)); + assert!(u64::from(r + 1) * u64::from(r + 1) > u64::from(n)); + } + } + + #[test] + fn scalar_equals_singleton_batch() { + let xs = stream(5000, 8000, 300, 11); + let mut a = RollingFloor::for_width(16384); + let mut b = RollingFloor::for_width(16384); + for &d in &xs { + let sa = a.observe(d); + let (used, sb) = b.observe_batch(&[d]); + assert_eq!(used, 1); + assert_eq!(sa, sb); + if let Some(s) = sa { + a.recalibrate(&s); + b.recalibrate(&s); + } + } + same_state(&a, &b); + } + + /// Any batching reproduces the scalar loop, including shifts raised in + /// the middle of a batch. + #[test] + fn arbitrary_batching_equals_scalar_stream() { + let mut xs = normalish(4000, 8192, 64, 3); + xs.extend(normalish(6000, 8900, 64, 4)); // forces mid-stream shifts + let (s, s_shifts) = run_scalar(RollingFloor::for_width(16384), &xs); + assert!(!s_shifts.is_empty(), "fixture must raise a shift"); + for chunk in [1, 7, 999, 1000, 1001, 2500, 10_000] { + let (b, b_shifts) = run_batched(RollingFloor::for_width(16384), &xs, chunk); + assert_eq!(s_shifts, b_shifts, "chunk {chunk}"); + same_state(&s, &b); + } + } + + #[test] + fn checkpoint_cadence_is_every_1000_after_the_first() { + // A floor far from the data: every checkpoint must alert. + let mut f = RollingFloor::from_params(100, 1); + let mut at = Vec::new(); + for i in 1..=4000u64 { + if f.observe(9000).is_some() { + at.push(i); + } + } + assert_eq!(at, [2000, 3000, 4000]); + } + + #[test] + fn drift_boundary_is_strict() { + // Calibrated σ = 8: shift iff |Δμ| > 4 or |Δσ| > 2. + for (value, expect) in [(104u32, false), (105, true)] { + let mut f = RollingFloor::from_params(100, 8); + // Alternate value ± 8 so the running σ is exactly 8. + let mut hit = None; + for i in 0..2000u32 { + hit = f.observe(if i % 2 == 0 { value - 8 } else { value + 8 }); + } + assert_eq!(hit.is_some(), expect, "mean {value}"); + } + for (half, expect) in [(10u32, false), (11, true)] { + let mut f = RollingFloor::from_params(100, 8); + let mut hit = None; + for i in 0..2000u32 { + hit = f.observe(if i % 2 == 0 { 100 - half } else { 100 + half }); + } + assert_eq!(hit.is_some(), expect, "σ {half}"); + } + } + + #[test] + fn recalibration_resets_the_running_state() { + let mut f = RollingFloor::calibrate(&normalish(2000, 5000, 50, 9)); + let shift = f + .observe_batch(&normalish(3000, 5400, 50, 10)) + .1 + .expect("shift"); + f.recalibrate(&shift); + assert_eq!((f.mu(), f.sigma()), (shift.new_mu, shift.new_sigma)); + assert_eq!(f.observations(), 0); + assert!(f.reservoir().is_empty()); + assert_eq!(f.reservoir().capacity(), RollingFloor::RESERVOIR_CAP); + assert!(!f.is_empirical()); + assert_eq!((f.skewness(), f.kurtosis()), (0, 300)); + assert_eq!(f.active_floors(), RollingFloor::sigma_floors_of(f.mu(), f.sigma())); + assert_eq!(f.empirical_cascade(), f.sigma_cascade()); + } + + #[test] + fn reservoir_is_deterministic_and_bounded() { + let xs = stream(20_000, 0, 1 << 20, 5); + let mut a = ReservoirU32::new(1000); + let mut b = ReservoirU32::new(1000); + xs.iter().for_each(|&d| a.observe(d)); + xs.iter().for_each(|&d| b.observe(d)); + assert_eq!(a, b); + assert_eq!(a.len(), 1000); + assert_eq!(a.seen(), 20_000); + // Replacement actually happens: the held set is not the first 1000. + assert_ne!(a.samples(), &xs[..1000]); + } + + #[test] + fn quantile_rule_is_floor_index_clamped() { + let s: Vec = (0..1000).collect(); + assert_eq!(quantile_of_sorted(&s, 0.0013), 1); + assert_eq!(quantile_of_sorted(&s, 0.159), 159); + assert_eq!(quantile_of_sorted(&s, 0.5), 500); + assert_eq!(quantile_of_sorted(&s, 1.0), 999); + assert_eq!(quantile_of_sorted(&[], 0.5), 0); + } + + #[test] + fn normal_stream_stays_sigma_and_bimodal_goes_empirical() { + let mut f = RollingFloor::calibrate(&normalish(1000, 8192, 64, 1)); + f.observe_batch(&normalish(1000, 8192, 64, 2)); + assert!(f.shape_is_normal(), "skew {} kurt {}", f.skewness(), f.kurtosis()); + assert!(!f.is_empirical()); + + // Symmetric two-mode mixture: kurtosis ≈ 100, far below the window. + let (lo, hi) = (normalish(1000, 7800, 20, 3), normalish(1000, 8600, 20, 4)); + let bimodal: Vec = lo.iter().zip(&hi).flat_map(|(&a, &b)| [a, b]).collect(); + let mut g = RollingFloor::calibrate(&bimodal[..1000]); + g.observe_batch(&bimodal[1000..]); + assert!(g.is_empirical(), "skew {} kurt {}", g.skewness(), g.kurtosis()); + let sorted = g.reservoir().sorted(); + assert_eq!(g.active_floors()[2], quantile_of_sorted(&sorted, 0.159)); + } + + /// Parameter drift and shape are separate readings. A pure location and + /// spread change of a normal stream raises a parameter shift. At that + /// checkpoint the reservoir still mixes the old and the new population, + /// so the shape reads non-normal: this is the reference behaviour, kept + /// as is. After recalibration the shape layer restarts and the same + /// shifted stream reads normal again. + #[test] + fn parameter_drift_then_recalibration_restores_the_normal_shape() { + let mut f = RollingFloor::calibrate(&normalish(1000, 5000, 40, 7)); + let shifted = normalish(5000, 5600, 90, 8); + let (used, shift) = f.observe_batch(&shifted); + let shift = shift.expect("parameters moved"); + assert!(!f.shape_is_normal(), "mixed reservoir at the drift checkpoint"); + assert_eq!(shift.observations, 2000, "first checkpoint still mixes both"); + f.recalibrate(&shift); + // The first shift adopts the mixture; the next one settles on the + // shifted stream's own parameters. + let (f, later) = run_scalar(f, &shifted[used..]); + assert_eq!(later.len(), 1, "{later:?}"); + assert!(f.mu().abs_diff(5600) <= 5 && f.sigma().abs_diff(90) <= 5, "{} {}", f.mu(), f.sigma()); + assert!(f.shape_is_normal(), "skew {} kurt {}", f.skewness(), f.kurtosis()); + assert!(!f.is_empirical()); + } + + /// Usable from the first observation, and refined by more of them without + /// any reset. + #[test] + fn anytime_use_refines_with_population() { + let mut f = RollingFloor::for_width(16384); + assert_eq!(f.active_floors(), [8000, 8064, 8128, 8192]); + let xs = normalish(50_000, 8192, 64, 12); + let mut errs = Vec::new(); + for (i, &d) in xs.iter().enumerate() { + assert!(f.observe(d).is_none(), "on-prior data must not drift"); + if [100, 1000, 50_000].contains(&(i + 1)) { + let m = f.moments(); + errs.push((m.variance().sqrt() - 64.0).abs()); + } + } + assert_eq!(f.observations(), 50_000); + assert!(errs[2] < errs[0], "{errs:?}"); + } + + /// Large-population running spread is still meaningful: the checkpoint + /// variance comes from exact moments past the u128 product range. + #[test] + fn large_population_variance_floor_is_exact() { + let half = 1u128 << 32; + let hi = u128::from(u32::MAX); + let m = MomentsU32 { + n: 1 << 33, + sum: half * (hi + hi - 2), + sum_sq: half * (hi * hi + (hi - 2) * (hi - 2)), + }; + assert!(u128::from(m.n).checked_mul(m.sum_sq).is_none()); + assert_eq!(variance_floor(&m), 1); // values ±1 around the mean + let m = MomentsU32 { + n: 4, + sum: 10, + sum_sq: 30, + }; // 1,2,3,4: var 1.25 + assert_eq!(variance_floor(&m), 1); + } + + #[test] + fn calibrate_uses_spread_around_the_integer_mean() { + // 0,0,0,1: mean 0.25, floor mean 0, Σ(x−0)² = 1, 1/4 = 0 -> σ 1 (floored). + let f = RollingFloor::calibrate(&[0, 0, 0, 1]); + assert_eq!((f.mu(), f.sigma()), (0, 1)); + let f = RollingFloor::calibrate(&[100, 120, 100, 120]); + assert_eq!((f.mu(), f.sigma()), (110, 10)); + assert_eq!(f.observations(), 4); + } + + // ── The lance-graph reference, kept verbatim as an oracle ─────────── + // + // Old `hdr.rs` arithmetic: integer Welford with truncated means. Used + // only to measure where exact moments change a floor decision. + struct LegacyWelford { + n: u64, + sum: u64, + m2: u64, + } + impl LegacyWelford { + fn observe(&mut self, d: u32) -> Option<(u32, u32)> { + let d = u64::from(d); + self.n += 1; + self.sum += d; + let old = if self.n > 1 { (self.sum - d) / (self.n - 1) } else { d }; + let new = self.sum / self.n; + self.m2 = self + .m2 + .wrapping_add(((d as i64 - old as i64) * (d as i64 - new as i64)) as u64); + if self.n.is_multiple_of(1000) && self.n > 1000 { + Some((new as u32, isqrt_u32((self.m2 / self.n) as u32).max(1))) + } else { + None + } + } + } + + /// At every checkpoint, the running (μ, σ) the drift rule sees match the + /// legacy integer Welford on these streams. + #[test] + fn checkpoint_parameters_match_legacy_welford() { + let mut checked = 0; + for (mu, sigma, seed) in [(8192, 64, 1), (5000, 50, 2), (100, 3, 3), (8192, 1, 4), (16384, 91, 5)] { + let xs = normalish(20_000, mu, sigma, seed); + let mut legacy = LegacyWelford { n: 0, sum: 0, m2: 0 }; + let mut m = MomentsU32::default(); + for &d in &xs { + m.observe(d); + if let Some((lmu, lsig)) = legacy.observe(d) { + let emu = saturate_u32(m.sum / u128::from(m.n)); + let esig = isqrt_u32(saturate_u32(variance_floor(&m))).max(1); + assert_eq!((emu, esig), (lmu, lsig), "mu {mu} sigma {sigma} n {}", m.n); + checked += 1; + } + } + } + assert_eq!(checked, 5 * 19); + } +} diff --git a/src/hpc/statistics.rs b/src/hpc/statistics.rs index de8bdb79..e9c40fc1 100644 --- a/src/hpc/statistics.rs +++ b/src/hpc/statistics.rs @@ -385,6 +385,16 @@ pub struct MomentsU32 { } impl MomentsU32 { + /// Fold one value in. Exactly equal to `self.merge(moments_u32(&[x]))`, + /// so a scalar stream and any batching of it reach the same state. + #[inline] + pub fn observe(&mut self, x: u32) { + let x = u128::from(x); + self.n += 1; + self.sum += x; + self.sum_sq += x * x; + } + /// Moments of the union of two samples — exact integer addition. #[inline] #[must_use] From 6ad34b2919bc5c8c55968de51ed7202a7b30e2b4 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 18:21:59 +0000 Subject: [PATCH 2/9] hdr: drop a redundant reservoir rebuild in RollingFloor::recalibrate from_params already constructs a fresh reservoir at RESERVOIR_CAP. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- src/hpc/rolling_floor.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index ca94599f..7caea304 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -421,9 +421,7 @@ impl RollingFloor { /// observation from scratch: running moments, reservoir, empirical mode /// and shape diagnostics are all reset. `σ` is floored at 1. pub fn recalibrate(&mut self, shift: &FloorShift) { - let capacity = self.reservoir.capacity(); *self = Self::from_params(shift.new_mu, shift.new_sigma.max(1)); - self.reservoir = ReservoirU32::new(capacity); } fn at_checkpoint(&self) -> bool { From 7e4a2e71376a8f98fadec6986c31c046c787b6b1 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 18:28:01 +0000 Subject: [PATCH 3/9] hdr: test the variance fallback branch and each kurtosis bound in isolation Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- src/hpc/rolling_floor.rs | 52 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index 7caea304..16ee1218 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -567,8 +567,15 @@ fn variance_floor(m: &MomentsU32) -> u128 { if let (Some(a), Some(b)) = (n.checked_mul(m.sum_sq), m.sum.checked_mul(m.sum)) { return (a - b) / (n * n); } - // n·M2 = C·n − r², with C centred on the floor mean. Write C = a·n + b; - // then ⌊(C·n − r²)/n²⌋ = a − [b·n < r²], and b·n, r² < 2^128. + variance_floor_centred(m) +} + +/// The overflow-free form of [`variance_floor`], valid for every `n > 0`. +/// `n·M2 = C·n − r²`, with `C` centred on the floor mean. Write +/// `C = a·n + b`; then `⌊(C·n − r²)/n²⌋ = a − [b·n < r²]`, and both `b·n` +/// and `r²` are below `2^128`. +fn variance_floor_centred(m: &MomentsU32) -> u128 { + let n = u128::from(m.n); let c = centred_on_floor_mean(m); let r = m.sum % n; let (a, b) = (c / n, c % n); @@ -848,6 +855,47 @@ mod tests { assert_eq!(variance_floor(&m), 1); } + /// The overflow-free variance form agrees with the direct one on every + /// small sample, including the fractional means that take the `a − 1` + /// correction. + #[test] + fn centred_variance_floor_matches_the_direct_form() { + let mut corrected = 0; + for seed in 1..400u64 { + let len = 2 + (seed % 37) as usize; + let xs = stream(len, (seed * 97 % 5000) as u32, 1 + (seed % 60) as u32, seed); + let m = moments_u32(&xs); + let n = u128::from(m.n); + let direct = (n * m.sum_sq - m.sum * m.sum) / (n * n); + assert_eq!(variance_floor_centred(&m), direct, "{xs:?}"); + let c = centred_on_floor_mean(&m); + corrected += usize::from((c % n) * n < (m.sum % n).pow(2)); + } + assert!(corrected > 0, "fixture must exercise the correction branch"); + } + + /// Each kurtosis bound switches to empirical floors on its own, with the + /// skew inside the window: a uniform stream is too light-tailed, a + /// narrow-core wide-tail mixture too heavy-tailed. + #[test] + fn kurtosis_alone_switches_to_empirical() { + let uniform = stream(2000, 8000, 400, 21); + let mut u = RollingFloor::calibrate(&uniform[..1000]); + u.observe_batch(&uniform[1000..]); + assert!(u.skewness().abs() < 2 && u.kurtosis() <= 200, "skew {} kurt {}", u.skewness(), u.kurtosis()); + assert!(u.is_empirical()); + + let (core, tail) = (normalish(2000, 8192, 10, 22), normalish(200, 8192, 120, 23)); + let mut mix = core; + for (i, t) in tail.into_iter().enumerate() { + mix[i * 9] = t; + } + let mut h = RollingFloor::calibrate(&mix[..1000]); + h.observe_batch(&mix[1000..]); + assert!(h.skewness().abs() < 2 && h.kurtosis() >= 500, "skew {} kurt {}", h.skewness(), h.kurtosis()); + assert!(h.is_empirical()); + } + #[test] fn calibrate_uses_spread_around_the_integer_mean() { // 0,0,0,1: mean 0.25, floor mean 0, Σ(x−0)² = 1, 1/4 = 0 -> σ 1 (floored). From 7269a86f786723cf421a3f557d16102e5ba669d6 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 18:32:40 +0000 Subject: [PATCH 4/9] hdr: pin the normality window at its boundaries Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- src/hpc/rolling_floor.rs | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index 16ee1218..cd40d534 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -896,6 +896,27 @@ mod tests { assert!(h.is_empirical()); } + /// The normality window at each of its boundaries. + #[test] + fn normality_window_boundaries() { + let mut f = RollingFloor::for_width(16384); + for (skew, kurt, normal) in [ + (0, 300, true), + (1, 300, true), + (-1, 300, true), + (2, 300, false), + (-2, 300, false), + (0, 200, false), + (0, 201, true), + (0, 499, true), + (0, 500, false), + ] { + f.skewness = skew; + f.kurtosis = kurt; + assert_eq!(f.shape_is_normal(), normal, "skew {skew} kurt {kurt}"); + } + } + #[test] fn calibrate_uses_spread_around_the_integer_mean() { // 0,0,0,1: mean 0.25, floor mean 0, Σ(x−0)² = 1, 1/4 = 0 -> σ 1 (floored). From fbfc63e1c9c02905ebf2247728e65d9e7a5b047e Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 18:41:30 +0000 Subject: [PATCH 5/9] hdr: gate hdr_rolling_floor_bench on std like the other hpc examples The no-default-features test build compiles every example; hpc exists only under std. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- Cargo.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index 652c2d0c..58a2c9af 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -207,6 +207,10 @@ required-features = ["std"] name = "rolling_floor_probe" required-features = ["std"] +[[example]] +name = "hdr_rolling_floor_bench" +required-features = ["std"] + [[example]] name = "codec_overlap_probe" required-features = ["std"] From c658012378c03323e03424ca93723b119692d0b1 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 19:23:57 +0000 Subject: [PATCH 6/9] =?UTF-8?q?hdr:=20thresholds=20on=20demand=20from=20an?= =?UTF-8?q?=20integer=20=CF=83-lattice?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RollingFloor no longer stores floor or cascade arrays. The accumulator keeps the facts (exact moments, reservoir, shape belief, drift anchor); every threshold is derived when asked. A SigmaLevel(k) is k quarter-σ below the noise floor. The Gaussian shape locates it as μ_t − k·σ_t/4 from the running moments; an empirical shape locates it at its Gaussian-equivalent tail rank Φ(−k/4)·len and moves the sample value into the current frame, μ_t + (x − μ_s)·σ_t/σ_s, in i128. shade(x, levels) counts how many lattice points a response survives. Coordinates are valid from the first observation (μ = x, σ = 0); only n = 0 falls back to the anchor. Recalibration moves the anchor and forgets the evidence but keeps the shape; a drift checkpoint does not read the shape from mixed-regime evidence. All integer; the f32 quantile path is gone. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- examples/hdr_rolling_floor_bench.rs | 38 +- src/hpc/rolling_floor.rs | 910 ++++++++++++++++++---------- 2 files changed, 621 insertions(+), 327 deletions(-) diff --git a/examples/hdr_rolling_floor_bench.rs b/examples/hdr_rolling_floor_bench.rs index c44da012..c48e086e 100644 --- a/examples/hdr_rolling_floor_bench.rs +++ b/examples/hdr_rolling_floor_bench.rs @@ -13,10 +13,14 @@ //! //! Periodic path, once per 1000 observations: //! 5. shape evaluation (sort 1000 samples, median, kurtosis) -//! 6. empirical quantiles (12 lookups into the sorted reservoir) +//! 6. empirical shape: locate 8 σ-lattice levels +//! +//! Query path, on demand (nothing is stored): +//! 7. Gaussian thresholds of 8 levels +//! 8. shade of a response over 8 levels use ndarray::hpc::bitwise::hamming_distance_raw; -use ndarray::hpc::rolling_floor::{quantile_of_sorted, ReservoirU32, RollingFloor}; +use ndarray::hpc::rolling_floor::{quantile_of_sorted, EmpiricalShape, ReservoirU32, RollingFloor, SigmaLevel}; use ndarray::hpc::statistics::{moments_u32, MomentsU32}; use std::hint::black_box; use std::time::Instant; @@ -110,15 +114,35 @@ fn main() { ns_per("5. shape: sort + median + kurtosis", K, || { for _ in 0..K { let sorted = black_box(&r).sorted(); - black_box(quantile_of_sorted(&sorted, 0.5)); + black_box(quantile_of_sorted(&sorted, 5000)); black_box(r.kurtosis(8192, 64)); } }); - let sorted = r.sorted(); - ns_per("6. empirical floors: 12 quantiles (sorted)", K, || { + let lattice = [4u8, 6, 7, 8, 9, 10, 11, 12].map(SigmaLevel); + let shape = EmpiricalShape::from_sample(r.samples()).unwrap(); + ns_per("6. empirical: locate 8 lattice levels", K, || { for _ in 0..K { - black_box(RollingFloor::FLOOR_PERCENTILES.map(|p| quantile_of_sorted(black_box(&sorted), p))); - black_box(RollingFloor::CASCADE_PERCENTILES.map(|p| quantile_of_sorted(black_box(&sorted), p))); + black_box(lattice.map(|l| black_box(&shape).locate(l, 8200, 70))); + } + }); + + println!("-- query path, on demand --"); + let mut g = RollingFloor::for_width(16384); + dists[..3000].iter().for_each(|&d| { + if let Some(s) = g.observe(d) { + g.recalibrate(&s); + } + }); + ns_per("7. Gaussian: thresholds of 8 levels", N, || { + for _ in 0..N { + black_box(black_box(&g).thresholds(&lattice)); } }); + ns_per("8. Gaussian: shade over 8 levels", N, || { + let mut acc = 0usize; + for &d in &dists { + acc += black_box(&g).shade(d, &lattice); + } + black_box(acc); + }); } diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index cd40d534..4c626cd6 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -1,46 +1,62 @@ -//! HDR rolling floor: an online distribution floor for popcount / Hamming -//! observations. +//! HDR rolling distribution: live statistics of popcount / Hamming +//! observations, and σ-lattice thresholds derived from them on demand. //! -//! This is the adaptive half of the HDR exposure meter, harvested from -//! lance-graph's `graph/blasgraph/hdr.rs`. That file stays the behavioural -//! reference: constants, cadence, drift rule and reset semantics are carried -//! over unchanged. What changed is the arithmetic underneath. The old -//! approximate integer Welford is replaced by the exact, mergeable -//! [`MomentsU32`]. +//! Harvested from lance-graph's `graph/blasgraph/hdr.rs`. Its constants +//! (reservoir capacity, cadence, minimum population, normality window, drift +//! rule) carry over unchanged; the arithmetic underneath is #327's exact, +//! mergeable [`MomentsU32`]. //! -//! # Two layers, two cadences +//! # The accumulator owns facts, the query owns the view //! -//! * **Parameters (continuous, cheap).** Every observation folds into exact -//! `(n, Σx, Σx²)`. Location and spread are available at any time, so the -//! floor can be used from the first observation and only gets better as the -//! population grows. There is no training barrier. -//! * **Shape (periodic, amortised).** A deterministic Algorithm-R reservoir -//! feeds median, empirical quantiles, skewness and kurtosis. They are -//! evaluated once every [`RollingFloor::EVAL_CADENCE`] observations, never on -//! the per-observation path. +//! What is stored: //! -//! # Floors +//! * **Moments** — exact `(n, Σx, Σx²)`, folded on every observation. The +//! current coordinates `(μ_t, σ_t)` follow from them at any time, from the +//! first observation on: observation and use are concurrent, and more +//! observations only sharpen the estimate. +//! * **Reservoir** — a deterministic Algorithm-R sample, the *evidence* for +//! the periodic shape check. +//! * **Shape** — the *belief* about the distribution family: [`Shape::Gaussian`] +//! (no data) or [`Shape::Empirical`] (a sorted sample and the `(μ_s, σ_s)` +//! frame it was measured in). +//! * **Anchor** — the calibrated `(μ, σ)`, used only as the reference the drift +//! rule measures against. //! -//! While the shape reads as normal, the floors are the analytical quantiles of -//! the calibrated `(μ, σ)`: `μ − kσ`, i.e. `μ + σ·Φ⁻¹(p)` at the percentiles -//! the empirical tables name. When the shape does not read as normal, the -//! floors are the reservoir's empirical quantiles at those same percentiles. +//! What is derived, never stored: every threshold. A detector asks for a +//! [`SigmaLevel`] — a point on the integer σ-lattice, `k` quarter-σ below the +//! noise floor — and the active shape locates it in the current coordinates: //! -//! # What is and is not mergeable +//! * Gaussian: `μ_t − k·σ_t/4`. +//! * Empirical: the sample value at `k`'s Gaussian-equivalent tail rank, +//! moved into the current frame, `μ_t + (x − μ_s)·σ_t/σ_s`. +//! +//! The percentile is not a second coordinate: it is fixed by `k` through +//! [`SigmaLevel::gaussian_tail_per_10000`]. A detector's shade of a response is +//! how far along its chosen lattice points the response survives +//! ([`RollingFloor::shade`]). Everything is integer. +//! +//! # Two cadences +//! +//! Moments and reservoir update on every observation. Every +//! [`RollingFloor::EVAL_CADENCE`] observations (after the first cadence) a +//! checkpoint checks drift and, when the parameters have not drifted, the +//! shape. //! -//! The moments are exact and partition-independent: shards can be merged in -//! any order. The reservoir is not. It is deterministic for a given -//! observation *order* and has no merge law. [`RollingFloor::observe_batch`] -//! therefore reproduces the scalar stream exactly rather than merging shards. +//! # Parameter drift is not shape drift //! -//! # Reference behaviour this preserves, including its limits +//! A drift alert means the running `(μ, σ)` left the anchor: the evidence now +//! spans two parameter regimes, so the shape is not judged from it at that +//! checkpoint. [`RollingFloor::recalibrate`] moves the anchor and forgets the +//! moments and the reservoir, but keeps the shape: a Gaussian whose `μ` and +//! `σ` moved is still Gaussian. The shape changes only at a drift-free +//! checkpoint. //! -//! * Floors move only when a drift alert is acted on -//! ([`RollingFloor::recalibrate`]), not continuously with the running -//! parameters. -//! * A parameter drift resets the shape layer too (reservoir, empirical mode, -//! skewness, kurtosis). The reference does not separate parameter drift from -//! shape drift; neither does this port. +//! # What is and is not mergeable +//! +//! The moments are exact and partition-independent. The reservoir is +//! deterministic for a given observation *order* and has no merge law; +//! [`RollingFloor::observe_batch`] therefore reproduces the scalar stream +//! exactly rather than merging shards. use super::statistics::{moments_u32, MomentsU32}; @@ -69,6 +85,43 @@ pub fn isqrt_u32(n: u32) -> u32 { } } +/// Rank `⌊per_10000 · len / 10000⌋`, clamped to the last index. `0` for an +/// empty sample. The integer rank rule for every empirical lookup. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::rank_per_10000; +/// assert_eq!(rank_per_10000(1000, 1587), 158); +/// assert_eq!(rank_per_10000(1000, 10_000), 999); +/// assert_eq!(rank_per_10000(0, 5000), 0); +/// ``` +pub fn rank_per_10000(len: usize, per_10000: u32) -> usize { + if len == 0 { + return 0; + } + ((u64::from(per_10000) * len as u64 / 10_000) as usize).min(len - 1) +} + +/// Empirical quantile of an ascending slice at `per_10000 / 10000`, by +/// [`rank_per_10000`]. `0` for an empty slice. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::quantile_of_sorted; +/// let s = [1, 2, 3, 4]; +/// assert_eq!(quantile_of_sorted(&s, 0), 1); +/// assert_eq!(quantile_of_sorted(&s, 5000), 3); +/// assert_eq!(quantile_of_sorted(&s, 10_000), 4); +/// ``` +pub fn quantile_of_sorted(sorted: &[u32], per_10000: u32) -> u32 { + if sorted.is_empty() { + return 0; + } + sorted[rank_per_10000(sorted.len(), per_10000)] +} + /// Deterministic reservoir sample of a `u32` stream (Vitter's Algorithm R). /// /// Every element seen so far has the same chance of being held. The @@ -85,7 +138,7 @@ pub fn isqrt_u32(n: u32) -> u32 { /// r.observe(d); /// } /// assert_eq!(r.len(), 3); -/// assert_eq!(r.quantile(0.5), 20); +/// assert_eq!(r.quantile(5000), 20); /// ``` #[derive(Debug, Clone, PartialEq, Eq)] pub struct ReservoirU32 { @@ -144,21 +197,17 @@ impl ReservoirU32 { &self.samples } - /// The held samples, sorted ascending. Sort once and use - /// [`quantile_of_sorted`] when several quantiles are needed. + /// The held samples, sorted ascending. pub fn sorted(&self) -> Vec { let mut s = self.samples.clone(); s.sort_unstable(); s } - /// Empirical quantile: the sorted sample at index `⌊q·len⌋`, clamped to - /// the last element. `0` for an empty reservoir. - /// - /// O(len · log len); meant for the periodic shape path, not per - /// observation. - pub fn quantile(&self, q: f32) -> u32 { - quantile_of_sorted(&self.sorted(), q) + /// Empirical quantile at `per_10000 / 10000` (see [`rank_per_10000`]). + /// O(len · log len); meant for the periodic shape path. + pub fn quantile(&self, per_10000: u32) -> u32 { + quantile_of_sorted(&self.sorted(), per_10000) } /// Pearson's second skewness `3(μ − median) / σ`, in integer division. @@ -168,7 +217,7 @@ impl ReservoirU32 { if sigma == 0 || self.samples.is_empty() { return 0; } - skewness_from_median(mu, sigma, self.quantile(0.5)) + skewness_from_median(mu, sigma, self.quantile(5000)) } /// Kurtosis ×100: `100 · E[(X − μ)⁴] / σ⁴` over the reservoir, with the @@ -202,40 +251,141 @@ impl ReservoirU32 { } } -/// Empirical quantile of an ascending slice: the element at `⌊q·len⌋`, -/// clamped to the last element. `0` for an empty slice. +fn skewness_from_median(mu: u32, sigma: u32, median: u32) -> i32 { + // i64 so no difference of two u32 values can overflow. + let s = 3 * (i64::from(mu) - i64::from(median)) / i64::from(sigma); + s.clamp(i64::from(i32::MIN), i64::from(i32::MAX)) as i32 +} + +/// `Φ(−k/4)`, the Gaussian lower-tail mass `k` quarter-σ below the mean, in +/// parts per 10 000 (rounded to nearest), for `k = 0..=16`. +const GAUSSIAN_TAIL_PER_10000: [u32; 17] = + [5000, 4013, 3085, 2266, 1587, 1056, 668, 401, 228, 122, 62, 30, 13, 6, 2, 1, 0]; + +/// A point on the integer σ-lattice: `k` quarter-σ below the noise floor. +/// +/// This is the identity of a sensitivity cut. A detector chooses which points +/// it asks for; the active [`Shape`] decides how each point is located in the +/// current distribution. `SigmaLevel(12)` is 3σ, `SigmaLevel(6)` is 1.5σ, +/// `SigmaLevel(0)` is the mean itself. /// /// # Example /// /// ``` -/// use ndarray::hpc::rolling_floor::quantile_of_sorted; -/// let s = [1, 2, 3, 4]; -/// assert_eq!(quantile_of_sorted(&s, 0.0), 1); -/// assert_eq!(quantile_of_sorted(&s, 0.5), 3); -/// assert_eq!(quantile_of_sorted(&s, 1.0), 4); +/// use ndarray::hpc::rolling_floor::SigmaLevel; +/// assert_eq!(SigmaLevel(4).gaussian_tail_per_10000(), 1587); // 1σ +/// assert_eq!(SigmaLevel(12).gaussian_tail_per_10000(), 13); // 3σ /// ``` -pub fn quantile_of_sorted(sorted: &[u32], q: f32) -> u32 { - if sorted.is_empty() { - return 0; +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub struct SigmaLevel(pub u8); + +impl SigmaLevel { + /// The lattice coordinate `k`, in quarter-σ. + pub const fn quarters(self) -> u32 { + self.0 as u32 + } + + /// The Gaussian-equivalent lower-tail mass of this level, `Φ(−k/4)`, in + /// parts per 10 000. This is how an empirical shape locates the same cut: + /// the level fixes the rank, no caller supplies a percentile. Levels + /// beyond 4σ (`k > 16`) have tail `0`, i.e. the sample minimum. + pub const fn gaussian_tail_per_10000(self) -> u32 { + let k = self.0 as usize; + if k < GAUSSIAN_TAIL_PER_10000.len() { + GAUSSIAN_TAIL_PER_10000[k] + } else { + 0 + } } - let idx = ((q * sorted.len() as f32) as usize).min(sorted.len() - 1); - sorted[idx] } -fn skewness_from_median(mu: u32, sigma: u32, median: u32) -> i32 { - // i64 so no difference of two u32 values can overflow. - let s = 3 * (i64::from(mu) - i64::from(median)) / i64::from(sigma); - s.clamp(i64::from(i32::MIN), i64::from(i32::MAX)) as i32 +/// The learned geometry of a non-Gaussian distribution: a sorted sample and +/// the `(μ_s, σ_s)` frame it was measured in. +/// +/// It answers a [`SigmaLevel`] by the sample value at the level's tail rank, +/// moved from its own frame into the current coordinates. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::{EmpiricalShape, SigmaLevel}; +/// let shape = EmpiricalShape::from_sample(&[10, 20, 30, 40]).unwrap(); +/// // In its own frame a level answers with the raw sample value. +/// let (mu, sigma) = (shape.mu(), shape.sigma()); +/// assert_eq!(shape.locate(SigmaLevel(0), mu, sigma), 30); +/// ``` +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct EmpiricalShape { + sorted: Vec, + mu: u32, + sigma: u32, +} + +impl EmpiricalShape { + /// Learn the shape from a sample: sort it and record its own floor mean + /// and `⌊√⌊M2/n⌋⌋` spread. `None` for an empty sample. + pub fn from_sample(sample: &[u32]) -> Option { + if sample.is_empty() { + return None; + } + let mut sorted = sample.to_vec(); + sorted.sort_unstable(); + let m = moments_u32(&sorted); + Some(Self { + mu: saturate_u32(m.sum / u128::from(m.n)), + sigma: isqrt_u32(saturate_u32(variance_floor(&m))), + sorted, + }) + } + + /// The sorted sample. + pub fn sorted(&self) -> &[u32] { + &self.sorted + } + + /// The sample's floor mean. + pub fn mu(&self) -> u32 { + self.mu + } + + /// The sample's spread. + pub fn sigma(&self) -> u32 { + self.sigma + } + + /// Locate `level` in the frame `(mu, sigma)`: + /// `mu + ⌊(x − μ_s)·sigma / σ_s⌋` with `x` the sample value at the level's + /// tail rank, clamped to `u32`. When the sample has no spread (`σ_s = 0`) + /// the offset `x − μ_s` is used unscaled. + pub fn locate(&self, level: SigmaLevel, mu: u32, sigma: u32) -> u32 { + let x = quantile_of_sorted(&self.sorted, level.gaussian_tail_per_10000()); + let delta = i128::from(x) - i128::from(self.mu); + let offset = if self.sigma == 0 { + delta + } else { + (delta * i128::from(sigma)).div_euclid(i128::from(self.sigma)) + }; + (i128::from(mu) + offset).clamp(0, i128::from(u32::MAX)) as u32 + } +} + +/// The distribution family a [`RollingFloor`] currently believes in. +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum Shape { + /// Normal: levels are located analytically as `μ − k·σ/4`. + Gaussian, + /// Not normal: levels are located through a learned sample. + Empirical(EmpiricalShape), } -/// A detected drift of the running parameters away from the calibrated ones. +/// A detected drift of the running parameters away from the anchor. #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub struct FloorShift { - /// Calibrated mean before the shift. + /// Anchor mean before the shift. pub old_mu: u32, /// Running mean at the checkpoint that raised the shift. pub new_mu: u32, - /// Calibrated standard deviation before the shift. + /// Anchor standard deviation before the shift. pub old_sigma: u32, /// Running standard deviation at that checkpoint. pub new_sigma: u32, @@ -243,98 +393,67 @@ pub struct FloorShift { pub observations: u64, } -/// The HDR rolling floor for a stream of `u32` distances. -/// -/// It holds the calibrated `(μ, σ)`, the sigma floors derived from them, and -/// empirical floors from a reservoir. Running parameters come from exact -/// [`MomentsU32`]. Every [`EVAL_CADENCE`](Self::EVAL_CADENCE) observations -/// (after the first cadence) it re-reads the distribution shape, selects sigma -/// or empirical floors, and checks for drift. -/// -/// The floors are lower-is-better thresholds, from strictest to loosest: -/// four band floors at `[μ−3σ, μ−2σ, μ−σ, μ]` and eight cascade floors at -/// quarter-σ steps from `μ−σ` to `μ−3σ`. How a caller turns them into bands -/// is its own policy. +/// Live distribution of a stream of `u32` distances, answering σ-lattice +/// thresholds on demand. /// /// # Example /// /// ``` -/// use ndarray::hpc::rolling_floor::RollingFloor; +/// use ndarray::hpc::rolling_floor::{RollingFloor, SigmaLevel}; /// let mut floor = RollingFloor::for_width(16384); -/// assert_eq!(floor.mu(), 8192); -/// assert_eq!(floor.sigma(), 64); -/// assert_eq!(floor.active_floors(), [8000, 8064, 8128, 8192]); -/// // Usable immediately; observations refine it as they arrive. +/// // No observations yet: the binomial prior μ = 8192, σ = 64 answers. +/// assert_eq!(floor.threshold(SigmaLevel(12)), 8000); +/// // The first observation is already a valid current state. +/// floor.observe(8100); +/// assert_eq!(floor.coordinates(), Some((8100, 0))); +/// // Three sensitivities of one detector; the response's shade is how many +/// // of them it survives. /// for d in 0..3000u32 { /// if let Some(shift) = floor.observe(8192 + (d % 7)) { /// floor.recalibrate(&shift); /// } /// } +/// let levels = [SigmaLevel(6), SigmaLevel(8), SigmaLevel(12)]; +/// assert_eq!(floor.shade(0, &levels), 3); +/// assert_eq!(floor.shade(u32::MAX, &levels), 0); /// ``` #[derive(Debug, Clone)] pub struct RollingFloor { - mu: u32, - sigma: u32, - sigma_floors: [u32; 4], - sigma_cascade: [u32; 8], + anchor_mu: u32, + anchor_sigma: u32, + moments: MomentsU32, reservoir: ReservoirU32, - empirical_floors: [u32; 4], - empirical_cascade: [u32; 8], - use_empirical: bool, + shape: Shape, skewness: i32, kurtosis: u32, - moments: MomentsU32, } impl RollingFloor { /// Reservoir capacity. pub const RESERVOIR_CAP: usize = 1000; - /// Shape and drift are evaluated when the observation count is a multiple - /// of this, and larger than it. + /// Checkpoints fall where the observation count is a multiple of this and + /// larger than it. pub const EVAL_CADENCE: u64 = 1000; - /// Minimum reservoir population before shape is evaluated. + /// Minimum reservoir population before the shape is judged. pub const MIN_SHAPE_SAMPLES: usize = 100; /// Kurtosis ×100 of the normal distribution. pub const NORMAL_KURTOSIS: u32 = 300; - /// Percentiles of the empirical band floors: ≈3σ, 2σ, 1σ, median. - pub const FLOOR_PERCENTILES: [f32; 4] = [0.001, 0.023, 0.159, 0.500]; - /// Percentiles of the empirical cascade floors: 1σ, 1.5σ, 1.75σ, 2σ, - /// 2.25σ, 2.5σ, 2.75σ, 3σ below the mean. - pub const CASCADE_PERCENTILES: [f32; 8] = [0.1587, 0.0668, 0.0401, 0.0228, 0.0122, 0.0062, 0.0030, 0.0013]; - - /// Band floors `[μ−3σ, μ−2σ, μ−σ, μ]`, saturating at zero. - pub fn sigma_floors_of(mu: u32, sigma: u32) -> [u32; 4] { - [mu.saturating_sub(3 * sigma), mu.saturating_sub(2 * sigma), mu.saturating_sub(sigma), mu] - } - - /// Cascade floors at `μ − kσ/4` for `k = 4, 6, 7, 8, 9, 10, 11, 12`, - /// saturating at zero. - pub fn sigma_cascade_of(mu: u32, sigma: u32) -> [u32; 8] { - [4, 6, 7, 8, 9, 10, 11, 12].map(|k: u32| mu.saturating_sub(k * sigma / 4)) - } - /// Floors that assume only a prior `(μ, σ)`, with no observations yet. + /// A floor with only a prior `(μ, σ)`: Gaussian shape, no observations. pub fn from_params(mu: u32, sigma: u32) -> Self { - let sigma_floors = Self::sigma_floors_of(mu, sigma); - let sigma_cascade = Self::sigma_cascade_of(mu, sigma); Self { - mu, - sigma, - sigma_floors, - sigma_cascade, + anchor_mu: mu, + anchor_sigma: sigma, + moments: MomentsU32::default(), reservoir: ReservoirU32::new(Self::RESERVOIR_CAP), - empirical_floors: sigma_floors, - empirical_cascade: sigma_cascade, - use_empirical: false, + shape: Shape::Gaussian, skewness: 0, kurtosis: Self::NORMAL_KURTOSIS, - moments: MomentsU32::default(), } } - /// Resume from calibrated `(μ, σ)` and already-accumulated running - /// moments, with an empty reservoir. The next checkpoint follows from - /// `moments.n`. + /// Resume from an anchor `(μ, σ)` and already-accumulated moments, with an + /// empty reservoir and Gaussian shape. pub fn from_params_and_moments(mu: u32, sigma: u32, moments: MomentsU32) -> Self { let mut floor = Self::from_params(mu, sigma); floor.moments = moments; @@ -349,10 +468,9 @@ impl RollingFloor { /// Calibrate from a warm-up sample of at least two distances. /// - /// `μ` is the floor of the sample mean and `σ` is - /// `max(1, ⌊√⌊Σ(x − μ)² / n⌋⌋)`, spread measured around that integer - /// mean, exactly as the reference does. The sample seeds the reservoir and - /// the running moments; the floors start in sigma mode. + /// The anchor is `μ = ⌊Σx/n⌋` and `σ = max(1, ⌊√⌊Σ(x − μ)²/n⌋⌋)`, spread + /// around that integer mean exactly as the reference does. The sample + /// seeds the moments and the reservoir; the shape starts Gaussian. /// /// # Panics /// @@ -362,41 +480,33 @@ impl RollingFloor { let moments = moments_u32(sample); let mu = saturate_u32(moments.sum / u128::from(moments.n)); let sigma = isqrt_u32(saturate_u32(centred_on_floor_mean(&moments) / u128::from(moments.n))).max(1); - let mut floor = Self::from_params(mu, sigma); + let mut floor = Self::from_params_and_moments(mu, sigma, moments); for &d in sample { floor.reservoir.observe(d); } - let sorted = floor.reservoir.sorted(); - floor.empirical_floors = Self::FLOOR_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); - floor.empirical_cascade = Self::CASCADE_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); - floor.moments = moments; floor } - /// Fold one observation in. Returns a shift when this observation lands - /// on a checkpoint and the running parameters have drifted: - /// `|μ_run − μ| > σ/2` or `|σ_run − σ| > σ/4`, against the calibrated - /// `(μ, σ)`. + /// Fold one observation in. At a checkpoint, returns a shift when the + /// running parameters have drifted from the anchor: + /// `|μ_run − μ| > σ/2` or `|σ_run − σ| > σ/4`. #[inline] pub fn observe(&mut self, distance: u32) -> Option { self.moments.observe(distance); self.reservoir.observe(distance); if self.at_checkpoint() { - self.evaluate() + self.checkpoint() } else { None } } /// Fold a batch in, stopping right after the first checkpoint that - /// raises a shift. + /// raises a shift. Returns how many values were consumed and the shift. /// - /// Returns how many values were consumed and the shift, if any. Feeding - /// the unconsumed rest in after acting on the shift reproduces the scalar - /// loop `if let Some(s) = observe(d) { recalibrate(&s) }` exactly, for any - /// batching: the moments between checkpoints go through - /// [`moments_u32`], the reservoir sees every value in stream order, and - /// every checkpoint is evaluated on the same state the scalar loop sees. + /// Feeding the unconsumed rest in after acting on the shift reproduces the + /// scalar loop `if let Some(s) = observe(d) { recalibrate(&s) }` exactly, + /// for any batching. pub fn observe_batch(&mut self, distances: &[u32]) -> (usize, Option) { let mut consumed = 0; while consumed < distances.len() { @@ -409,7 +519,7 @@ impl RollingFloor { } consumed += take; if take == to_checkpoint && self.at_checkpoint() { - if let Some(shift) = self.evaluate() { + if let Some(shift) = self.checkpoint() { return (consumed, Some(shift)); } } @@ -417,11 +527,14 @@ impl RollingFloor { (consumed, None) } - /// Adopt the shifted parameters as the new calibration and restart - /// observation from scratch: running moments, reservoir, empirical mode - /// and shape diagnostics are all reset. `σ` is floored at 1. + /// Adopt the shifted parameters as the new anchor and forget the + /// accumulated evidence (moments and reservoir). The shape is kept: + /// parameter drift is not shape drift. `σ` is floored at 1. pub fn recalibrate(&mut self, shift: &FloorShift) { - *self = Self::from_params(shift.new_mu, shift.new_sigma.max(1)); + self.anchor_mu = shift.new_mu; + self.anchor_sigma = shift.new_sigma.max(1); + self.moments = MomentsU32::default(); + self.reservoir = ReservoirU32::new(Self::RESERVOIR_CAP); } fn at_checkpoint(&self) -> bool { @@ -429,37 +542,37 @@ impl RollingFloor { n.is_multiple_of(Self::EVAL_CADENCE) && n > Self::EVAL_CADENCE } - /// The periodic path: shape evaluation, floor selection, drift check. - fn evaluate(&mut self) -> Option { + /// The periodic path: drift first; the shape only when there is none. + fn checkpoint(&mut self) -> Option { let run_mu = saturate_u32(self.moments.sum / u128::from(self.moments.n)); let run_sigma = isqrt_u32(saturate_u32(variance_floor(&self.moments))).max(1); + let mu_drift = run_mu.abs_diff(self.anchor_mu); + let sigma_drift = run_sigma.abs_diff(self.anchor_sigma); + if mu_drift > self.anchor_sigma / 2 || sigma_drift > self.anchor_sigma / 4 { + // The evidence spans two parameter regimes; do not read the shape + // from it. + return Some(FloorShift { + old_mu: self.anchor_mu, + new_mu: run_mu, + old_sigma: self.anchor_sigma, + new_sigma: run_sigma, + observations: self.moments.n, + }); + } + if self.reservoir.len() >= Self::MIN_SHAPE_SAMPLES { let sorted = self.reservoir.sorted(); - self.skewness = skewness_from_median(run_mu, run_sigma, quantile_of_sorted(&sorted, 0.5)); + self.skewness = skewness_from_median(run_mu, run_sigma, quantile_of_sorted(&sorted, 5000)); self.kurtosis = self.reservoir.kurtosis(run_mu, run_sigma); - if self.shape_is_normal() { - self.use_empirical = false; + self.shape = if self.shape_is_normal() { + Shape::Gaussian } else { - self.empirical_floors = Self::FLOOR_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); - self.empirical_cascade = Self::CASCADE_PERCENTILES.map(|p| quantile_of_sorted(&sorted, p)); - self.use_empirical = true; - } - } - - let mu_drift = run_mu.abs_diff(self.mu); - let sigma_drift = run_sigma.abs_diff(self.sigma); - if mu_drift > self.sigma / 2 || sigma_drift > self.sigma / 4 { - Some(FloorShift { - old_mu: self.mu, - new_mu: run_mu, - old_sigma: self.sigma, - new_sigma: run_sigma, - observations: self.moments.n, - }) - } else { - None + // `sorted` is non-empty here. + EmpiricalShape::from_sample(&sorted).map_or(Shape::Gaussian, Shape::Empirical) + }; } + None } /// The reference normality window: `|skew| < 2` and `200 < kurt < 500`. @@ -467,66 +580,79 @@ impl RollingFloor { self.skewness.abs() < 2 && self.kurtosis > 200 && self.kurtosis < 500 } - /// Calibrated mean. - pub fn mu(&self) -> u32 { - self.mu + /// Current coordinates `(μ_t, σ_t)` from the running moments: + /// `⌊Σx/n⌋` and `⌊√⌊M2/n⌋⌋`. `None` only before the first observation; + /// after one observation they are `(x, 0)`. + pub fn coordinates(&self) -> Option<(u32, u32)> { + if self.moments.n == 0 { + return None; + } + Some(( + saturate_u32(self.moments.sum / u128::from(self.moments.n)), + isqrt_u32(saturate_u32(variance_floor(&self.moments))), + )) } - /// Calibrated standard deviation. - pub fn sigma(&self) -> u32 { - self.sigma + /// The coordinates a threshold is located in: the running ones, or the + /// anchor when nothing has been observed yet. + fn frame(&self) -> (u32, u32) { + self.coordinates() + .unwrap_or((self.anchor_mu, self.anchor_sigma)) } - /// Band floors derived from the calibrated `(μ, σ)`. - pub fn sigma_floors(&self) -> [u32; 4] { - self.sigma_floors + fn locate(&self, level: SigmaLevel, mu: u32, sigma: u32) -> u32 { + match &self.shape { + Shape::Gaussian => mu.saturating_sub(level.quarters() * sigma / 4), + Shape::Empirical(e) => e.locate(level, mu, sigma), + } } - /// Cascade floors derived from the calibrated `(μ, σ)`. - pub fn sigma_cascade(&self) -> [u32; 8] { - self.sigma_cascade + /// The threshold of one σ-lattice point in the current distribution. + /// With `σ_t = 0` every Gaussian level sits at `μ_t`. + pub fn threshold(&self, level: SigmaLevel) -> u32 { + let (mu, sigma) = self.frame(); + self.locate(level, mu, sigma) } - /// Band floors from the reservoir's empirical quantiles. - pub fn empirical_floors(&self) -> [u32; 4] { - self.empirical_floors + /// Thresholds of several lattice points, reading the coordinates once. + pub fn thresholds(&self, levels: &[SigmaLevel; N]) -> [u32; N] { + let (mu, sigma) = self.frame(); + levels.map(|l| self.locate(l, mu, sigma)) } - /// Cascade floors from the reservoir's empirical quantiles. - pub fn empirical_cascade(&self) -> [u32; 8] { - self.empirical_cascade + /// How many of `levels` the response `x` survives: the number whose + /// threshold lies strictly above `x` (lower distance is stronger). This is + /// the response's shade on the detector's own lattice. + pub fn shade(&self, x: u32, levels: &[SigmaLevel; N]) -> usize { + self.thresholds(levels).iter().filter(|&&t| x < t).count() } - /// The band floors in effect: empirical when the shape read as non-normal - /// at the last checkpoint, sigma otherwise. - pub fn active_floors(&self) -> [u32; 4] { - if self.use_empirical { - self.empirical_floors - } else { - self.sigma_floors - } + /// The current shape belief. + pub fn shape(&self) -> &Shape { + &self.shape } - /// The cascade floors in effect, chosen like [`active_floors`](Self::active_floors). - pub fn active_cascade(&self) -> [u32; 8] { - if self.use_empirical { - self.empirical_cascade - } else { - self.sigma_cascade - } + /// Whether the shape is empirical. + pub fn is_empirical(&self) -> bool { + matches!(self.shape, Shape::Empirical(_)) } - /// Whether the empirical floors are in effect. - pub fn is_empirical(&self) -> bool { - self.use_empirical + /// Anchor mean, the drift reference. + pub fn mu(&self) -> u32 { + self.anchor_mu } - /// Skewness at the last shape evaluation (`0` before any). + /// Anchor standard deviation, the drift reference. + pub fn sigma(&self) -> u32 { + self.anchor_sigma + } + + /// Skewness at the last shape check (`0` before any). pub fn skewness(&self) -> i32 { self.skewness } - /// Kurtosis ×100 at the last shape evaluation (300 before any). + /// Kurtosis ×100 at the last shape check (300 before any). pub fn kurtosis(&self) -> u32 { self.kurtosis } @@ -541,7 +667,7 @@ impl RollingFloor { self.moments.n } - /// The reservoir behind the empirical floors. + /// The reservoir, the shape check's evidence. pub fn reservoir(&self) -> &ReservoirU32 { &self.reservoir } @@ -646,13 +772,25 @@ mod tests { (f, shifts) } + const LATTICE: [SigmaLevel; 9] = [ + SigmaLevel(0), + SigmaLevel(4), + SigmaLevel(6), + SigmaLevel(7), + SigmaLevel(8), + SigmaLevel(9), + SigmaLevel(10), + SigmaLevel(11), + SigmaLevel(12), + ]; + fn same_state(a: &RollingFloor, b: &RollingFloor) { assert_eq!(a.moments(), b.moments()); assert_eq!(a.reservoir(), b.reservoir()); assert_eq!((a.mu(), a.sigma()), (b.mu(), b.sigma())); - assert_eq!(a.active_floors(), b.active_floors()); - assert_eq!(a.active_cascade(), b.active_cascade()); - assert_eq!((a.is_empirical(), a.skewness(), a.kurtosis()), (b.is_empirical(), b.skewness(), b.kurtosis())); + assert_eq!(a.shape(), b.shape()); + assert_eq!((a.skewness(), a.kurtosis()), (b.skewness(), b.kurtosis())); + assert_eq!(a.thresholds(&LATTICE), b.thresholds(&LATTICE)); } #[test] @@ -664,6 +802,75 @@ mod tests { } } + /// `Φ(−k/4)` at the quarter-σ lattice, per 10 000. + #[test] + fn gaussian_tail_table_is_phi() { + let expect = [ + (0, 5000), + (4, 1587), + (6, 668), + (7, 401), + (8, 228), + (9, 122), + (10, 62), + (11, 30), + (12, 13), + (16, 0), + (40, 0), + ]; + for (k, v) in expect { + assert_eq!(SigmaLevel(k).gaussian_tail_per_10000(), v, "k {k}"); + } + // Monotone: a deeper cut is rarer. + for k in 0..16u8 { + assert!(SigmaLevel(k).gaussian_tail_per_10000() >= SigmaLevel(k + 1).gaussian_tail_per_10000()); + } + } + + /// The integer rank rule against the reference `f32` rule + /// `⌊(q as f32)·len⌋`, for every reservoir length 1..=1000. + /// + /// The eight cascade levels match the reference exactly. The three old + /// band percentiles (0.159, 0.023, 0.001) were coarse approximations of + /// 1σ, 2σ and 3σ; unifying them onto the lattice ranks changes the sample + /// index at exactly these many lengths, pinned so the delta is explicit. + #[test] + fn integer_rank_rule_against_the_reference() { + let f32_rank = |q: f32, len: usize| ((q * len as f32) as usize).min(len - 1); + let cascade = [ + (4u8, 0.1587f32), + (6, 0.0668), + (7, 0.0401), + (8, 0.0228), + (9, 0.0122), + (10, 0.0062), + (11, 0.0030), + (12, 0.0013), + ]; + for len in 1..=1000 { + for (k, q) in cascade { + assert_eq!( + rank_per_10000(len, SigmaLevel(k).gaussian_tail_per_10000()), + f32_rank(q, len), + "k {k} len {len}" + ); + } + assert_eq!(rank_per_10000(len, SigmaLevel(0).gaussian_tail_per_10000()), f32_rank(0.5, len)); + } + let changed = |k: u8, q: f32| { + (1..=1000) + .filter(|&len| rank_per_10000(len, SigmaLevel(k).gaussian_tail_per_10000()) != f32_rank(q, len)) + .count() + }; + assert_eq!(changed(4, 0.159), 150); + assert_eq!(changed(8, 0.023), 98); + assert_eq!(changed(12, 0.001), 230); + // At the full reservoir: 1σ 159 → 158, 2σ 23 → 22, 3σ 1 → 1. + assert_eq!((f32_rank(0.159, 1000), rank_per_10000(1000, 1587)), (159, 158)); + assert_eq!((f32_rank(0.023, 1000), rank_per_10000(1000, 228)), (23, 22)); + assert_eq!((f32_rank(0.001, 1000), rank_per_10000(1000, 13)), (1, 1)); + } + #[test] fn scalar_equals_singleton_batch() { let xs = stream(5000, 8000, 300, 11); @@ -699,7 +906,7 @@ mod tests { #[test] fn checkpoint_cadence_is_every_1000_after_the_first() { - // A floor far from the data: every checkpoint must alert. + // An anchor far from the data: every checkpoint must alert. let mut f = RollingFloor::from_params(100, 1); let mut at = Vec::new(); for i in 1..=4000u64 { @@ -712,10 +919,9 @@ mod tests { #[test] fn drift_boundary_is_strict() { - // Calibrated σ = 8: shift iff |Δμ| > 4 or |Δσ| > 2. + // Anchor σ = 8: shift iff |Δμ| > 4 or |Δσ| > 2. for (value, expect) in [(104u32, false), (105, true)] { let mut f = RollingFloor::from_params(100, 8); - // Alternate value ± 8 so the running σ is exactly 8. let mut hit = None; for i in 0..2000u32 { hit = f.observe(if i % 2 == 0 { value - 8 } else { value + 8 }); @@ -732,22 +938,111 @@ mod tests { } } + /// The ballot box: no running coordinates before the first observation, + /// a valid (noisy) state after it, and σ = 0 collapses every Gaussian cut + /// onto μ rather than borrowing the anchor. + #[test] + fn coordinates_are_valid_from_the_first_observation() { + let mut f = RollingFloor::from_params(5000, 50); + assert_eq!(f.coordinates(), None); + assert_eq!(f.threshold(SigmaLevel(8)), 4900, "prior answers before any observation"); + f.observe(7000); + assert_eq!(f.coordinates(), Some((7000, 0))); + assert_eq!(f.thresholds(&[SigmaLevel(0), SigmaLevel(4), SigmaLevel(12)]), [7000, 7000, 7000]); + let levels = [SigmaLevel(4), SigmaLevel(8), SigmaLevel(12)]; + assert_eq!(f.shade(6999, &levels), 3); + assert_eq!(f.shade(7000, &levels), 0); + f.observe(7010); + assert_eq!(f.coordinates(), Some((7005, 5))); + } + + /// Gaussian thresholds follow the running coordinates on every + /// observation, with no recalibration: exactly `μ_t − k·σ_t/4`. + #[test] + fn gaussian_thresholds_follow_current_moments() { + let mut f = RollingFloor::calibrate(&normalish(1000, 5000, 40, 5)); + let anchored = f.thresholds(&LATTICE); + let mut moved = false; + // A slow drift, small enough never to raise a shift. + for (i, d) in normalish(900, 5015, 40, 6).into_iter().enumerate() { + assert!(f.observe(d).is_none()); + let (mu, s) = f.coordinates().unwrap(); + let want = LATTICE.map(|l| mu.saturating_sub(l.quarters() * s / 4)); + assert_eq!(f.thresholds(&LATTICE), want, "observation {i}"); + moved |= want != anchored; + } + assert!(moved, "the floor must move without recalibration"); + } + + /// Half-σ lattice points reproduce the integer `SigmaGate::custom` tiers + /// bit for bit. #[test] - fn recalibration_resets_the_running_state() { - let mut f = RollingFloor::calibrate(&normalish(2000, 5000, 50, 9)); - let shift = f - .observe_batch(&normalish(3000, 5400, 50, 10)) - .1 - .expect("shift"); + fn half_sigma_levels_match_the_sigma_gate() { + use crate::hpc::kernels::SigmaGate; + for (mu, s) in [(8192u32, 64u32), (5000, 51), (1000, 7), (100, 1)] { + let f = RollingFloor::from_params(mu, s); + let g = SigmaGate::custom(mu, s); + let t = f.thresholds(&[SigmaLevel(12), SigmaLevel(10), SigmaLevel(8), SigmaLevel(6)]); + assert_eq!(t, [g.discovery, g.strong, g.evidence, g.hint], "mu {mu} sigma {s}"); + } + } + + #[test] + fn shade_counts_the_levels_a_response_survives() { + let f = RollingFloor::from_params(8192, 64); + let levels = [SigmaLevel(12), SigmaLevel(6), SigmaLevel(8)]; // order is free + // Cuts at 8000 (3σ), 8096 (1.5σ), 8064 (2σ). + assert_eq!(f.shade(7999, &levels), 3); + assert_eq!(f.shade(8000, &levels), 2); + assert_eq!(f.shade(8063, &levels), 2); + assert_eq!(f.shade(8064, &levels), 1); + assert_eq!(f.shade(8095, &levels), 1); + assert_eq!(f.shade(8096, &levels), 0); + } + + /// Recalibration forgets the evidence and moves the anchor, but keeps the + /// shape. + #[test] + fn recalibration_forgets_evidence_and_keeps_shape() { + let (lo, hi) = (normalish(1000, 7800, 20, 3), normalish(1000, 8600, 20, 4)); + let bimodal: Vec = lo.iter().zip(&hi).flat_map(|(&a, &b)| [a, b]).collect(); + let mut f = RollingFloor::calibrate(&bimodal[..1000]); + f.observe_batch(&bimodal[1000..]); + assert!(f.is_empirical()); + let shape = f.shape().clone(); + let shift = FloorShift { + old_mu: f.mu(), + new_mu: 9000, + old_sigma: f.sigma(), + new_sigma: 300, + observations: f.observations(), + }; f.recalibrate(&shift); - assert_eq!((f.mu(), f.sigma()), (shift.new_mu, shift.new_sigma)); + assert_eq!((f.mu(), f.sigma()), (9000, 300)); assert_eq!(f.observations(), 0); assert!(f.reservoir().is_empty()); assert_eq!(f.reservoir().capacity(), RollingFloor::RESERVOIR_CAP); - assert!(!f.is_empirical()); - assert_eq!((f.skewness(), f.kurtosis()), (0, 300)); - assert_eq!(f.active_floors(), RollingFloor::sigma_floors_of(f.mu(), f.sigma())); - assert_eq!(f.empirical_cascade(), f.sigma_cascade()); + assert_eq!(f.shape(), &shape); + } + + /// A pure location and spread change of a normal stream is parameter + /// drift only: the drift checkpoint does not read the mixed evidence as a + /// new shape, and the shape stays Gaussian throughout. + #[test] + fn pure_parameter_shift_keeps_the_gaussian_shape() { + let mut f = RollingFloor::calibrate(&normalish(1000, 5000, 40, 7)); + let shifted = normalish(8000, 5600, 90, 8); + let (used, shift) = f.observe_batch(&shifted); + let shift = shift.expect("parameters moved"); + assert_eq!(shift.observations, 2000); + assert_eq!(f.shape(), &Shape::Gaussian, "mixed evidence must not relearn the shape"); + f.recalibrate(&shift); + let (f, later) = run_scalar(f, &shifted[used..]); + assert!(later.len() <= 1, "{later:?}"); + assert_eq!(f.shape(), &Shape::Gaussian); + assert!(f.shape_is_normal(), "skew {} kurt {}", f.skewness(), f.kurtosis()); + let (mu, s) = f.coordinates().unwrap(); + assert!(mu.abs_diff(5600) <= 5 && s.abs_diff(90) <= 5, "{mu} {s}"); } #[test] @@ -760,82 +1055,104 @@ mod tests { assert_eq!(a, b); assert_eq!(a.len(), 1000); assert_eq!(a.seen(), 20_000); - // Replacement actually happens: the held set is not the first 1000. assert_ne!(a.samples(), &xs[..1000]); } #[test] - fn quantile_rule_is_floor_index_clamped() { - let s: Vec = (0..1000).collect(); - assert_eq!(quantile_of_sorted(&s, 0.0013), 1); - assert_eq!(quantile_of_sorted(&s, 0.159), 159); - assert_eq!(quantile_of_sorted(&s, 0.5), 500); - assert_eq!(quantile_of_sorted(&s, 1.0), 999); - assert_eq!(quantile_of_sorted(&[], 0.5), 0); - } - - #[test] - fn normal_stream_stays_sigma_and_bimodal_goes_empirical() { + fn normal_stream_stays_gaussian_and_bimodal_goes_empirical() { let mut f = RollingFloor::calibrate(&normalish(1000, 8192, 64, 1)); f.observe_batch(&normalish(1000, 8192, 64, 2)); assert!(f.shape_is_normal(), "skew {} kurt {}", f.skewness(), f.kurtosis()); - assert!(!f.is_empirical()); + assert_eq!(f.shape(), &Shape::Gaussian); - // Symmetric two-mode mixture: kurtosis ≈ 100, far below the window. let (lo, hi) = (normalish(1000, 7800, 20, 3), normalish(1000, 8600, 20, 4)); let bimodal: Vec = lo.iter().zip(&hi).flat_map(|(&a, &b)| [a, b]).collect(); let mut g = RollingFloor::calibrate(&bimodal[..1000]); g.observe_batch(&bimodal[1000..]); assert!(g.is_empirical(), "skew {} kurt {}", g.skewness(), g.kurtosis()); - let sorted = g.reservoir().sorted(); - assert_eq!(g.active_floors()[2], quantile_of_sorted(&sorted, 0.159)); } - /// Parameter drift and shape are separate readings. A pure location and - /// spread change of a normal stream raises a parameter shift. At that - /// checkpoint the reservoir still mixes the old and the new population, - /// so the shape reads non-normal: this is the reference behaviour, kept - /// as is. After recalibration the shape layer restarts and the same - /// shifted stream reads normal again. + /// An empirical shape answers with the raw sample quantile in its own + /// frame, translates with μ and scales with σ. #[test] - fn parameter_drift_then_recalibration_restores_the_normal_shape() { - let mut f = RollingFloor::calibrate(&normalish(1000, 5000, 40, 7)); - let shifted = normalish(5000, 5600, 90, 8); - let (used, shift) = f.observe_batch(&shifted); - let shift = shift.expect("parameters moved"); - assert!(!f.shape_is_normal(), "mixed reservoir at the drift checkpoint"); - assert_eq!(shift.observations, 2000, "first checkpoint still mixes both"); - f.recalibrate(&shift); - // The first shift adopts the mixture; the next one settles on the - // shifted stream's own parameters. - let (f, later) = run_scalar(f, &shifted[used..]); - assert_eq!(later.len(), 1, "{later:?}"); - assert!(f.mu().abs_diff(5600) <= 5 && f.sigma().abs_diff(90) <= 5, "{} {}", f.mu(), f.sigma()); - assert!(f.shape_is_normal(), "skew {} kurt {}", f.skewness(), f.kurtosis()); - assert!(!f.is_empirical()); + fn empirical_shape_projects_into_current_coordinates() { + let (lo, hi) = (normalish(500, 7800, 20, 3), normalish(500, 8600, 20, 4)); + let sample: Vec = lo.into_iter().chain(hi).collect(); + let e = EmpiricalShape::from_sample(&sample).unwrap(); + let (mu, s) = (e.mu(), e.sigma()); + for l in LATTICE { + let x = quantile_of_sorted(e.sorted(), l.gaussian_tail_per_10000()); + assert_eq!(e.locate(l, mu, s), x, "own frame, k {}", l.0); + assert_eq!(e.locate(l, mu + 100, s), x + 100, "translation, k {}", l.0); + let d = i64::from(x) - i64::from(mu); + let doubled = i64::from(mu) + (d * 2 * i64::from(s)).div_euclid(i64::from(s)); + assert_eq!(i64::from(e.locate(l, mu, 2 * s)), doubled, "scale, k {}", l.0); + } + assert_eq!(EmpiricalShape::from_sample(&[]), None); + // No spread: the offset is used unscaled. + let flat = EmpiricalShape::from_sample(&[7, 7, 7]).unwrap(); + assert_eq!((flat.sigma(), flat.locate(SigmaLevel(12), 100, 50)), (0, 100)); } - /// Usable from the first observation, and refined by more of them without + /// Each kurtosis bound switches to the empirical shape on its own, with + /// the skew inside the window. + #[test] + fn kurtosis_alone_switches_to_empirical() { + let uniform = stream(2000, 8000, 400, 21); + let mut u = RollingFloor::calibrate(&uniform[..1000]); + u.observe_batch(&uniform[1000..]); + assert!(u.skewness().abs() < 2 && u.kurtosis() <= 200, "skew {} kurt {}", u.skewness(), u.kurtosis()); + assert!(u.is_empirical()); + + let (core, tail) = (normalish(2000, 8192, 10, 22), normalish(200, 8192, 120, 23)); + let mut mix = core; + for (i, t) in tail.into_iter().enumerate() { + mix[i * 9] = t; + } + let mut h = RollingFloor::calibrate(&mix[..1000]); + h.observe_batch(&mix[1000..]); + assert!(h.skewness().abs() < 2 && h.kurtosis() >= 500, "skew {} kurt {}", h.skewness(), h.kurtosis()); + assert!(h.is_empirical()); + } + + /// The normality window at each of its boundaries. + #[test] + fn normality_window_boundaries() { + let mut f = RollingFloor::for_width(16384); + for (skew, kurt, normal) in [ + (0, 300, true), + (1, 300, true), + (-1, 300, true), + (2, 300, false), + (-2, 300, false), + (0, 200, false), + (0, 201, true), + (0, 499, true), + (0, 500, false), + ] { + f.skewness = skew; + f.kurtosis = kurt; + assert_eq!(f.shape_is_normal(), normal, "skew {skew} kurt {kurt}"); + } + } + + /// Usable from the first observation, refined by more of them without /// any reset. #[test] fn anytime_use_refines_with_population() { let mut f = RollingFloor::for_width(16384); - assert_eq!(f.active_floors(), [8000, 8064, 8128, 8192]); let xs = normalish(50_000, 8192, 64, 12); let mut errs = Vec::new(); for (i, &d) in xs.iter().enumerate() { assert!(f.observe(d).is_none(), "on-prior data must not drift"); if [100, 1000, 50_000].contains(&(i + 1)) { - let m = f.moments(); - errs.push((m.variance().sqrt() - 64.0).abs()); + errs.push(f.coordinates().unwrap().1.abs_diff(64)); } } assert_eq!(f.observations(), 50_000); - assert!(errs[2] < errs[0], "{errs:?}"); + assert!(errs[2] <= errs[0], "{errs:?}"); } - /// Large-population running spread is still meaningful: the checkpoint - /// variance comes from exact moments past the u128 product range. #[test] fn large_population_variance_floor_is_exact() { let half = 1u128 << 32; @@ -846,12 +1163,12 @@ mod tests { sum_sq: half * (hi * hi + (hi - 2) * (hi - 2)), }; assert!(u128::from(m.n).checked_mul(m.sum_sq).is_none()); - assert_eq!(variance_floor(&m), 1); // values ±1 around the mean + assert_eq!(variance_floor(&m), 1); let m = MomentsU32 { n: 4, sum: 10, sum_sq: 30, - }; // 1,2,3,4: var 1.25 + }; assert_eq!(variance_floor(&m), 1); } @@ -874,52 +1191,8 @@ mod tests { assert!(corrected > 0, "fixture must exercise the correction branch"); } - /// Each kurtosis bound switches to empirical floors on its own, with the - /// skew inside the window: a uniform stream is too light-tailed, a - /// narrow-core wide-tail mixture too heavy-tailed. - #[test] - fn kurtosis_alone_switches_to_empirical() { - let uniform = stream(2000, 8000, 400, 21); - let mut u = RollingFloor::calibrate(&uniform[..1000]); - u.observe_batch(&uniform[1000..]); - assert!(u.skewness().abs() < 2 && u.kurtosis() <= 200, "skew {} kurt {}", u.skewness(), u.kurtosis()); - assert!(u.is_empirical()); - - let (core, tail) = (normalish(2000, 8192, 10, 22), normalish(200, 8192, 120, 23)); - let mut mix = core; - for (i, t) in tail.into_iter().enumerate() { - mix[i * 9] = t; - } - let mut h = RollingFloor::calibrate(&mix[..1000]); - h.observe_batch(&mix[1000..]); - assert!(h.skewness().abs() < 2 && h.kurtosis() >= 500, "skew {} kurt {}", h.skewness(), h.kurtosis()); - assert!(h.is_empirical()); - } - - /// The normality window at each of its boundaries. - #[test] - fn normality_window_boundaries() { - let mut f = RollingFloor::for_width(16384); - for (skew, kurt, normal) in [ - (0, 300, true), - (1, 300, true), - (-1, 300, true), - (2, 300, false), - (-2, 300, false), - (0, 200, false), - (0, 201, true), - (0, 499, true), - (0, 500, false), - ] { - f.skewness = skew; - f.kurtosis = kurt; - assert_eq!(f.shape_is_normal(), normal, "skew {skew} kurt {kurt}"); - } - } - #[test] fn calibrate_uses_spread_around_the_integer_mean() { - // 0,0,0,1: mean 0.25, floor mean 0, Σ(x−0)² = 1, 1/4 = 0 -> σ 1 (floored). let f = RollingFloor::calibrate(&[0, 0, 0, 1]); assert_eq!((f.mu(), f.sigma()), (0, 1)); let f = RollingFloor::calibrate(&[100, 120, 100, 120]); @@ -928,9 +1201,6 @@ mod tests { } // ── The lance-graph reference, kept verbatim as an oracle ─────────── - // - // Old `hdr.rs` arithmetic: integer Welford with truncated means. Used - // only to measure where exact moments change a floor decision. struct LegacyWelford { n: u64, sum: u64, From 59276d6c3a5cf0cf6aa30bb8f2e2f55b892e71df Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 19:25:27 +0000 Subject: [PATCH 7/9] hdr: pin the empirical rounding rule and the unscaled-offset branch Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- src/hpc/rolling_floor.rs | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index 4c626cd6..d5b09e19 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -1088,10 +1088,23 @@ mod tests { let doubled = i64::from(mu) + (d * 2 * i64::from(s)).div_euclid(i64::from(s)); assert_eq!(i64::from(e.locate(l, mu, 2 * s)), doubled, "scale, k {}", l.0); } + // Rounding is floor (toward −∞), also for offsets below the mean. + let below = SigmaLevel(12); + let x = quantile_of_sorted(e.sorted(), below.gaussian_tail_per_10000()); + let d = i64::from(x) - i64::from(mu); + assert!( + d < 0 && (d * i64::from(s + 1)) % i64::from(s) != 0, + "fixture must exercise a fractional negative offset" + ); + assert_eq!( + i64::from(e.locate(below, mu, s + 1)), + i64::from(mu) + (d * i64::from(s + 1)).div_euclid(i64::from(s)) + ); assert_eq!(EmpiricalShape::from_sample(&[]), None); - // No spread: the offset is used unscaled. - let flat = EmpiricalShape::from_sample(&[7, 7, 7]).unwrap(); - assert_eq!((flat.sigma(), flat.locate(SigmaLevel(12), 100, 50)), (0, 100)); + // Floor spread 0 with a non-zero offset: the offset is used unscaled. + let flat = EmpiricalShape::from_sample(&[7, 8, 8, 8]).unwrap(); + assert_eq!((flat.mu(), flat.sigma()), (7, 0)); + assert_eq!(flat.locate(SigmaLevel(0), 100, 50), 101); } /// Each kurtosis bound switches to the empirical shape on its own, with From 525055e2c579355de626ac06879d8534c32a03f8 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 19:57:38 +0000 Subject: [PATCH 8/9] hdr: state the SigmaLevel domain explicitly MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Gaussian answers every k in 0..=255. Empirical answers k in 0..=16 from the fixed integer Φ(−k/4) table, one rank per level; k > 16 resolves to the sample minimum. No interpolation, no float. The eight cascade cuts are ordinary lattice points, pinned by a test over off-cascade k. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- src/hpc/rolling_floor.rs | 52 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index d5b09e19..0c16a78e 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -269,6 +269,24 @@ const GAUSSIAN_TAIL_PER_10000: [u32; 17] = /// current distribution. `SigmaLevel(12)` is 3σ, `SigmaLevel(6)` is 1.5σ, /// `SigmaLevel(0)` is the mean itself. /// +/// # Domain +/// +/// The lattice is every `k` in `0..=255`. No fixed set of cuts is built in; +/// the eight cascade cuts `4, 6, 7, …, 12` a consumer uses today are ordinary +/// points on it. +/// +/// * **Gaussian shape**: every `k` is answered exactly as `μ − k·σ/4` +/// (saturating at 0). +/// * **Empirical shape**: every `k` in `0..=16` (0 to 4σ) has its own +/// Gaussian-equivalent rank, from a fixed integer table of `Φ(−k/4)` +/// rounded to parts per 10 000. Every `k > 16` has tail mass below +/// 0.5 / 10 000 and resolves to rank 0, the sample minimum. Nothing is +/// interpolated and nothing is evaluated in floating point. +/// +/// Rank resolution is also bounded by the sample. With `len` samples, two +/// levels land on distinct ranks only if their `⌊tail·len/10000⌋` differ. At +/// the 1000-sample reservoir, every `k ≥ 13` already reads the minimum. +/// /// # Example /// /// ``` @@ -287,8 +305,9 @@ impl SigmaLevel { /// The Gaussian-equivalent lower-tail mass of this level, `Φ(−k/4)`, in /// parts per 10 000. This is how an empirical shape locates the same cut: - /// the level fixes the rank, no caller supplies a percentile. Levels - /// beyond 4σ (`k > 16`) have tail `0`, i.e. the sample minimum. + /// the level fixes the rank, no caller supplies a percentile. Defined from + /// the table for `k = 0..=16`; `k > 16` returns `0` (the sample minimum), + /// see the type docs. pub const fn gaussian_tail_per_10000(self) -> u32 { let k = self.0 as usize; if k < GAUSSIAN_TAIL_PER_10000.len() { @@ -871,6 +890,35 @@ mod tests { assert_eq!((f32_rank(0.001, 1000), rank_per_10000(1000, 13)), (1, 1)); } + /// The lattice is every `k`, not the eight historical cascade cuts: + /// off-cascade points answer in both shapes, the empirical table covers + /// `0..=16` with its own ranks, and beyond it saturates to the minimum. + #[test] + fn sigma_level_domain_is_the_whole_lattice() { + let g = RollingFloor::from_params(8192, 64); + for k in [1u8, 2, 3, 5, 13, 14, 15, 16, 17, 40, 128] { + assert_eq!(g.threshold(SigmaLevel(k)), 8192 - u32::from(k) * 16, "k {k}"); + } + assert_eq!(g.threshold(SigmaLevel(255)), 8192 - 255 * 16); + assert_eq!(RollingFloor::from_params(100, 64).threshold(SigmaLevel(255)), 0, "saturates"); + + // A 10 000-sample empirical shape resolves every table entry to its + // own rank, including the off-cascade ones. + let sample: Vec = (0..10_000).collect(); + let e = EmpiricalShape::from_sample(&sample).unwrap(); + let located: Vec = (0..=16u8) + .map(|k| e.locate(SigmaLevel(k), e.mu(), e.sigma())) + .collect(); + let want: Vec = (0..=16u8) + .map(|k| SigmaLevel(k).gaussian_tail_per_10000()) + .collect(); + assert_eq!(located, want, "rank = tail · len / 10000 on 0..10000"); + for k in [17u8, 40, 255] { + assert_eq!(SigmaLevel(k).gaussian_tail_per_10000(), 0); + assert_eq!(e.locate(SigmaLevel(k), e.mu(), e.sigma()), 0, "k {k}: sample minimum"); + } + } + #[test] fn scalar_equals_singleton_batch() { let xs = stream(5000, 8000, 300, 11); From aeda1f46174192b5e11f409b9e2c5c0ffa201211 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 20:09:24 +0000 Subject: [PATCH 9/9] =?UTF-8?q?hdr:=20full=20u32=20range=20=E2=80=94=20wid?= =?UTF-8?q?e=20sqrt,=20wide=20k=C2=B7=CF=83,=20overflow-safe=20kurtosis?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Codex P2 on #328: - variance of u32 values reaches 2^62; take its square root in u128 (isqrt_u128) instead of narrowing to u32 first, which clamped σ at 65 535 (calibrate(&[0, u32::MAX]) now gives σ = 2 147 483 647). - form k·σ in u128 before the saturating subtraction. - kurtosis: exact when Σd⁴ and ·100 fit u128; otherwise a 16.16 fixed-point (d/σ)² path within one unit of the exact form, saturating only where the true value is astronomically large. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- src/hpc/rolling_floor.rs | 139 ++++++++++++++++++++++++++++++++------- 1 file changed, 117 insertions(+), 22 deletions(-) diff --git a/src/hpc/rolling_floor.rs b/src/hpc/rolling_floor.rs index 0c16a78e..98dae99f 100644 --- a/src/hpc/rolling_floor.rs +++ b/src/hpc/rolling_floor.rs @@ -85,6 +85,31 @@ pub fn isqrt_u32(n: u32) -> u32 { } } +/// Floor of `√n` for any `u128`, integer Newton iteration. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::rolling_floor::isqrt_u128; +/// assert_eq!(isqrt_u128(1 << 64), 1 << 32); +/// assert_eq!(isqrt_u128(u128::MAX), u128::from(u64::MAX)); +/// ``` +pub fn isqrt_u128(n: u128) -> u128 { + if n == 0 { + return 0; + } + // Start >= floor(√n) so Newton descends monotonically; x ≤ 2^64 and + // n / x ≤ 2^64, so x + n / x cannot overflow. + let mut x = 1u128 << ((129 - n.leading_zeros()) / 2); + loop { + let x1 = (x + n / x) / 2; + if x1 >= x { + return x; + } + x = x1; + } +} + /// Rank `⌊per_10000 · len / 10000⌋`, clamped to the last index. `0` for an /// empty sample. The integer rank rule for every empirical lookup. /// @@ -227,19 +252,7 @@ impl ReservoirU32 { if sigma == 0 || self.samples.len() < 4 { return 300; } - let n = self.samples.len() as u128; - // u128: a fourth power of a u32 difference is below 2^128. - let m4: u128 = self - .samples - .iter() - .map(|&d| { - let diff = u128::from(d.abs_diff(mu)); - diff * diff * diff * diff - }) - .sum::() - / n; - let s4 = u128::from(sigma).pow(4); - u32::try_from(m4 * 100 / s4).unwrap_or(u32::MAX) + kurtosis_x100_exact(&self.samples, mu, sigma).unwrap_or_else(|| kurtosis_x100_scaled(&self.samples, mu, sigma)) } /// Deterministic splitmix64-style hash that drives replacement. @@ -352,7 +365,7 @@ impl EmpiricalShape { let m = moments_u32(&sorted); Some(Self { mu: saturate_u32(m.sum / u128::from(m.n)), - sigma: isqrt_u32(saturate_u32(variance_floor(&m))), + sigma: sqrt_u32(variance_floor(&m)), sorted, }) } @@ -498,7 +511,7 @@ impl RollingFloor { assert!(sample.len() > 1, "need at least 2 samples to calibrate"); let moments = moments_u32(sample); let mu = saturate_u32(moments.sum / u128::from(moments.n)); - let sigma = isqrt_u32(saturate_u32(centred_on_floor_mean(&moments) / u128::from(moments.n))).max(1); + let sigma = sqrt_u32(centred_on_floor_mean(&moments) / u128::from(moments.n)).max(1); let mut floor = Self::from_params_and_moments(mu, sigma, moments); for &d in sample { floor.reservoir.observe(d); @@ -564,7 +577,7 @@ impl RollingFloor { /// The periodic path: drift first; the shape only when there is none. fn checkpoint(&mut self) -> Option { let run_mu = saturate_u32(self.moments.sum / u128::from(self.moments.n)); - let run_sigma = isqrt_u32(saturate_u32(variance_floor(&self.moments))).max(1); + let run_sigma = sqrt_u32(variance_floor(&self.moments)).max(1); let mu_drift = run_mu.abs_diff(self.anchor_mu); let sigma_drift = run_sigma.abs_diff(self.anchor_sigma); @@ -606,10 +619,7 @@ impl RollingFloor { if self.moments.n == 0 { return None; } - Some(( - saturate_u32(self.moments.sum / u128::from(self.moments.n)), - isqrt_u32(saturate_u32(variance_floor(&self.moments))), - )) + Some((saturate_u32(self.moments.sum / u128::from(self.moments.n)), sqrt_u32(variance_floor(&self.moments)))) } /// The coordinates a threshold is located in: the running ones, or the @@ -621,7 +631,9 @@ impl RollingFloor { fn locate(&self, level: SigmaLevel, mu: u32, sigma: u32) -> u32 { match &self.shape { - Shape::Gaussian => mu.saturating_sub(level.quarters() * sigma / 4), + Shape::Gaussian => { + saturate_u32(u128::from(mu).saturating_sub(u128::from(level.quarters()) * u128::from(sigma) / 4)) + } Shape::Empirical(e) => e.locate(level, mu, sigma), } } @@ -692,6 +704,43 @@ impl RollingFloor { } } +/// `⌊100 · E[(X − μ)⁴] / σ⁴⌋` exactly, or `None` when an intermediate +/// leaves `u128` (only for spreads far beyond any popcount width). +fn kurtosis_x100_exact(samples: &[u32], mu: u32, sigma: u32) -> Option { + let mut sum: u128 = 0; + for &d in samples { + let diff = u128::from(d.abs_diff(mu)); + sum = sum.checked_add((diff * diff).checked_mul(diff * diff)?)?; + } + let m4 = sum / samples.len() as u128; + Some(u32::try_from(m4.checked_mul(100)? / u128::from(sigma).pow(4)).unwrap_or(u32::MAX)) +} + +/// The same ratio from per-sample `(d/σ)²` in 16.16 fixed point, used only +/// when the exact form overflows. Where it would overflow too the kurtosis is +/// astronomically large and saturates at `u32::MAX`. Agrees with the exact +/// form to within one unit where both fit. +fn kurtosis_x100_scaled(samples: &[u32], mu: u32, sigma: u32) -> u32 { + let s2 = u128::from(sigma) * u128::from(sigma); + let mut sum: u128 = 0; + for &d in samples { + let diff = u128::from(d.abs_diff(mu)); + let r = ((diff * diff) << 16) / s2; // (d/σ)², 16 fractional bits + match r.checked_mul(r).and_then(|t| sum.checked_add(t)) { + Some(v) => sum = v, + None => return u32::MAX, + } + } + let m4 = sum / samples.len() as u128; // 32 fractional bits + u32::try_from(m4.saturating_mul(100) >> 32).unwrap_or(u32::MAX) +} + +/// `⌊√n⌋` for a variance of `u32` values. That variance is below `2^62`, so +/// the root fits `u32`. +fn sqrt_u32(n: u128) -> u32 { + saturate_u32(isqrt_u128(n)) +} + fn saturate_u32(x: u128) -> u32 { u32::try_from(x).unwrap_or(u32::MAX) } @@ -919,6 +968,52 @@ mod tests { } } + /// The full `u32` range: σ above 65 535, `k·σ` above `u32`, fourth-power + /// sums above `u128` — none may clamp, wrap or panic. + #[test] + fn full_u32_range_does_not_clamp_or_overflow() { + const Q: u32 = u32::MAX / 2; // 2 147 483 647 + // Variance of {0, MAX} is Q² + Q (floor); its root is Q, not 65 535. + let f = RollingFloor::calibrate(&[0, u32::MAX]); + assert_eq!((f.mu(), f.sigma()), (Q, Q)); + assert_eq!(f.coordinates(), Some((Q, Q))); + let e = EmpiricalShape::from_sample(&[0, u32::MAX]).unwrap(); + assert_eq!((e.mu(), e.sigma()), (Q, Q)); + + // k·σ is formed wide: 4 · 1.5e9 does not fit u32. + let g = RollingFloor::from_params(u32::MAX, 1_500_000_000); + assert_eq!(g.threshold(SigmaLevel(4)), u32::MAX - 1_500_000_000); + assert_eq!(g.threshold(SigmaLevel(12)), 0, "saturates, never wraps"); + + // A two-point distribution has kurtosis exactly 1, i.e. 100. + let mut r = ReservoirU32::new(1000); + (0..1000u32).for_each(|i| r.observe(if i % 2 == 0 { 0 } else { u32::MAX })); + assert!((99..=101).contains(&r.kurtosis(Q, Q)), "{}", r.kurtosis(Q, Q)); + } + + /// Where the exact kurtosis fits, the overflow-safe path agrees with it + /// to within one unit of the ×100 scale. + #[test] + fn kurtosis_fallback_matches_the_exact_path() { + for (mu, sigma, seed) in [(8192, 64, 1), (5000, 3, 2), (100_000, 900, 3)] { + let xs = normalish(1000, mu, sigma, seed); + let mut r = ReservoirU32::new(1000); + xs.iter().for_each(|&d| r.observe(d)); + let exact = kurtosis_x100_exact(r.samples(), mu, sigma).expect("fits"); + let scaled = kurtosis_x100_scaled(r.samples(), mu, sigma); + assert!(exact.abs_diff(scaled) <= 1, "exact {exact} scaled {scaled}"); + } + } + + #[test] + fn isqrt_u128_is_floor_sqrt() { + for n in (0..100_000u128).chain([u128::MAX, u128::MAX - 1, (1 << 64) - 1, 1 << 64, (1 << 126) + 12345]) { + let r = isqrt_u128(n); + assert!(r.checked_mul(r).is_some_and(|sq| sq <= n), "n {n}"); + assert!((r + 1).checked_mul(r + 1).is_none_or(|sq| sq > n), "n {n}"); + } + } + #[test] fn scalar_equals_singleton_batch() { let xs = stream(5000, 8000, 300, 11); @@ -1298,7 +1393,7 @@ mod tests { m.observe(d); if let Some((lmu, lsig)) = legacy.observe(d) { let emu = saturate_u32(m.sum / u128::from(m.n)); - let esig = isqrt_u32(saturate_u32(variance_floor(&m))).max(1); + let esig = sqrt_u32(variance_floor(&m)).max(1); assert_eq!((emu, esig), (lmu, lsig), "mu {mu} sigma {sigma} n {}", m.n); checked += 1; }