diff --git a/crates/wasm-simd-parity/src/lib.rs b/crates/wasm-simd-parity/src/lib.rs index b84c4a6c..247e74d5 100644 --- a/crates/wasm-simd-parity/src/lib.rs +++ b/crates/wasm-simd-parity/src/lib.rs @@ -35,6 +35,9 @@ pub extern "C" fn selfcheck() -> u32 { if let Err(code) = check_i32x16_compare() { return code; } + if let Err(code) = check_zgamma_golden() { + return code; + } 0 } @@ -424,3 +427,38 @@ fn check_i32x16_compare() -> Result<(), u32> { } Ok(()) } + +/// `ZGamma` codes are bit-exact on every target: the same grid and pinned +/// digest as `ndarray::hpc::zspace::golden::GOLDEN_CODES` (keep in sync). +/// Before the deterministic `ln`, 624 of these 8 539 z values differed between +/// x86-64 glibc and wasm32, so this is the target where drift would show. +fn check_zgamma_golden() -> Result<(), u32> { + use ndarray::hpc::zspace::ZGamma; + const GOLDEN_CODES: u64 = 0xfb0b_294d_2a65_3dbb; + let mut grid: Vec = (-4096..=4096).map(|i| i as f32 / 4096.0).collect(); + let mut r = 0.999f32; + while r < 1.0 { + grid.push(r); + grid.push(-r); + r = f32::from_bits(r.to_bits() + 97); + } + if grid.len() != 8539 { + return Err(0x500); + } + let env = ZGamma::fit(&grid); + let mut codes = vec![0i8; grid.len()]; + env.encode_batch(&grid, &mut codes); + let mut h = 0xcbf2_9ce4_8422_2325u64; + for &c in &codes { + h = (h ^ u64::from(c as u8)).wrapping_mul(0x0100_0000_01b3); + } + if h != GOLDEN_CODES { + return Err(0x501); + } + for (&c, &k) in grid.iter().zip(&codes) { + if env.encode(c) != k { + return Err(0x502); + } + } + Ok(()) +} diff --git a/src/hpc/cascade.rs b/src/hpc/cascade.rs index 4673c10d..4fb1d99b 100644 --- a/src/hpc/cascade.rs +++ b/src/hpc/cascade.rs @@ -208,6 +208,56 @@ impl Cascade { } } + /// Fold a whole batch of distances into the rolling floor at once. + /// + /// The batch is reduced to exact integer moments with + /// [`moments_u32`](crate::hpc::statistics::moments_u32) and merged into + /// the running `(n, μ, σ)` with the parallel-variance merge + /// (Chan, Golub & LeVeque), so the resulting `mu`/`sigma`/`observations` + /// match calling [`observe`](Self::observe) once per distance up to f64 + /// rounding. Shards computed on different threads can each be folded in + /// this way, in any order. + /// + /// Drift is judged once per batch, not per element: an alert fires when + /// the batch moves μ by more than 2σ of the pre-batch state (and the state + /// had > 10 observations and σ > 0), the same test `observe` applies to a + /// single step. + pub fn observe_batch(&mut self, distances: &[u32]) -> Option { + let b = crate::hpc::statistics::moments_u32(distances); + if b.n == 0 { + return None; + } + let old_mu = self.mu; + let old_sigma = self.sigma; + let old_n = self.observations; + let n_a = old_n as f64; + let n_b = b.n as f64; + let n = n_a + n_b; + let (mean_b, m2_b) = (b.mean(), b.variance() * n_b); + if old_n == 0 { + self.mu = mean_b; + self.sigma = (m2_b / n_b).sqrt(); + } else { + let delta = mean_b - old_mu; + self.mu = old_mu + delta * n_b / n; + let m2 = old_sigma * old_sigma * n_a + m2_b + delta * delta * n_a * n_b / n; + self.sigma = (m2 / n).sqrt(); + } + self.observations = old_n + b.n as usize; + + if old_n > 10 && old_sigma > 0.0 && (self.mu - old_mu).abs() > 2.0 * old_sigma { + Some(ShiftAlert { + old_mu, + new_mu: self.mu, + old_sigma, + new_sigma: self.sigma, + observations: self.observations, + }) + } else { + None + } + } + pub fn recalibrate(&mut self, alert: &ShiftAlert) { self.mu = alert.new_mu; self.sigma = alert.new_sigma; @@ -768,6 +818,76 @@ mod tests { assert_eq!(got, exact_hits(&expected, threshold)); } + fn close(a: f64, b: f64) -> bool { + (a - b).abs() <= 1e-9 * a.abs().max(b.abs()).max(1.0) + } + + fn noisy(n: usize, base: u32, spread: u32, mut s: u64) -> Vec { + (0..n) + .map(|_| { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + base + (s as u32) % spread + }) + .collect() + } + + /// One `observe_batch` lands on the same rolling floor as `observe` + /// called per distance, from an empty state and from a warm one. + #[test] + fn observe_batch_matches_sequential_observe() { + let warm = noisy(300, 8000, 400, 1); + let batch = noisy(2000, 8100, 500, 2); + for start in [&[][..], &warm[..]] { + let mut seq = Cascade::from_threshold(8000, 2048); + let mut bat = Cascade::from_threshold(8000, 2048); + for &d in start { + seq.observe(d); + bat.observe(d); + } + for &d in &batch { + seq.observe(d); + } + bat.observe_batch(&batch); + assert_eq!(seq.observations(), bat.observations()); + assert!(close(seq.mu(), bat.mu()), "mu {} vs {}", seq.mu(), bat.mu()); + assert!(close(seq.sigma(), bat.sigma()), "sigma {} vs {}", seq.sigma(), bat.sigma()); + } + } + + /// Shard-parallel use: folding shards in any order gives the same floor + /// as folding the whole batch. + #[test] + fn observe_batch_shards_merge_in_any_order() { + let x = noisy(3001, 8000, 700, 3); + let mut whole = Cascade::from_threshold(8000, 2048); + whole.observe_batch(&x); + let shards: Vec<&[u32]> = x.chunks(640).collect(); + for order in [[0usize, 1, 2, 3, 4], [4, 2, 0, 3, 1]] { + let mut c = Cascade::from_threshold(8000, 2048); + for i in order { + c.observe_batch(shards[i]); + } + assert!(close(c.mu(), whole.mu())); + assert!(close(c.sigma(), whole.sigma())); + assert_eq!(c.observations(), whole.observations()); + } + } + + /// The drift alert can fire (a batch from a distribution shifted far + /// past 2σ) and stays silent on a batch from the same distribution. + #[test] + fn observe_batch_alerts_on_a_shift_only() { + let mut c = Cascade::from_threshold(8000, 2048); + assert!(c.observe_batch(&noisy(500, 8000, 100, 4)).is_none(), "first batch has no prior"); + assert!(c.observe_batch(&noisy(500, 8000, 100, 5)).is_none(), "same distribution"); + let alert = c + .observe_batch(&noisy(5000, 9000, 100, 6)) + .expect("shifted batch must alert"); + assert!(alert.new_mu > alert.old_mu + 2.0 * alert.old_sigma); + } + #[test] fn packed_database_roundtrip() { let vec_bytes = 256; diff --git a/src/hpc/mod.rs b/src/hpc/mod.rs index 074921b6..8067cf8f 100644 --- a/src/hpc/mod.rs +++ b/src/hpc/mod.rs @@ -25,6 +25,9 @@ pub mod blas_level2; pub mod blas_level3; pub mod reductions; pub mod statistics; +/// z-space entry points: Fisher-Z for cosine-shaped values, the binomial +/// null for bitpacked Hamming distances. +pub mod zspace; /// 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/reliability.rs b/src/hpc/reliability.rs index 56adf785..a529ea32 100644 --- a/src/hpc/reliability.rs +++ b/src/hpc/reliability.rs @@ -215,6 +215,12 @@ pub fn icc_a1(ratings: &[&[f64]]) -> f64 { /// One call computes all four coefficients plus the relative-L2 error and /// cosine similarity, so a harness can print a row per codec/flavor without /// recomputing means four times. +/// +/// `pearson`, `spearman` and `cosine` are cosine-shaped and stay raw for +/// display; anything that thresholds, averages or compares them goes through +/// [`pearson_z`](Self::pearson_z) / [`spearman_z`](Self::spearman_z) / +/// [`cosine_z`](Self::cosine_z) (the Fisher-Z entry point in +/// [`zspace`](crate::hpc::zspace)). #[derive(Debug, Clone, Copy, PartialEq)] pub struct FidelityReport { /// Pearson product-moment correlation (linear association). @@ -246,7 +252,9 @@ impl FidelityReport { /// let truth = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]; /// let est = [0.1, 0.9, 2.1, 2.9, 4.2, 4.8]; /// let r = FidelityReport::compute(&truth, &est); - /// assert!(r.pearson > 0.99 && r.spearman > 0.99); + /// // Cosine-shaped coefficients are consumed in Fisher-Z space. + /// use ndarray::hpc::zspace::fisher_z; + /// assert!(r.pearson_z() > fisher_z(0.99) && r.spearman_z() > fisher_z(0.99)); /// assert!(r.rel_l2 < 0.1); /// // Mismatched lengths → degenerate (does NOT truncate-then-score): /// let bad = FidelityReport::compute(&[1.0, 2.0, 100.0], &[1.0, 2.0]); @@ -295,12 +303,42 @@ impl FidelityReport { cosine, } } + + /// Fisher-Z of [`pearson`](Self::pearson) — the form every threshold, + /// average or confidence interval over it must use. + pub fn pearson_z(&self) -> f64 { + crate::hpc::zspace::fisher_z(self.pearson) + } + + /// Fisher-Z of [`spearman`](Self::spearman). + pub fn spearman_z(&self) -> f64 { + crate::hpc::zspace::fisher_z(self.spearman) + } + + /// Fisher-Z of [`cosine`](Self::cosine). + pub fn cosine_z(&self) -> f64 { + crate::hpc::zspace::fisher_z(self.cosine) + } } #[cfg(test)] mod tests { use super::*; + /// The z accessors are exactly Fisher-Z of the raw coefficients, and a + /// perfect score stays finite (the rim clamp, not `atanh(1) = inf`). + #[test] + fn fidelity_z_accessors_are_fisher_z() { + use crate::hpc::zspace::fisher_z; + let truth = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]; + let r = FidelityReport::compute(&truth, &[0.1, 0.9, 2.1, 2.9, 4.2, 4.8]); + assert_eq!(r.pearson_z(), fisher_z(r.pearson)); + assert_eq!(r.spearman_z(), fisher_z(r.spearman)); + assert_eq!(r.cosine_z(), fisher_z(r.cosine)); + let perfect = FidelityReport::compute(&truth, &truth); + assert!(perfect.pearson_z().is_finite() && perfect.pearson_z() > 10.0); + } + #[test] fn pearson_perfect_and_anti() { let x = [1.0, 2.0, 3.0, 4.0, 5.0]; diff --git a/src/hpc/statistics.rs b/src/hpc/statistics.rs index 14d97010..b08b0124 100644 --- a/src/hpc/statistics.rs +++ b/src/hpc/statistics.rs @@ -360,3 +360,203 @@ mod tests { } } } + +// ── Batch moments for shard-parallel Welford ────────────────────────────── + +/// Exact first and second moments of a `u32` sample: count, `Σx` and `Σx²`. +/// +/// These three integers are the sufficient statistics of a Welford rolling +/// floor. Unlike a running `(mean, M2)` pair they merge by plain addition, so +/// [`MomentsU32::merge`] is exact, associative and commutative: shards of a +/// sample can be reduced in parallel, in any order, and combined into the +/// same result as one sequential pass. Floats appear only in [`mean`] and +/// [`variance`], at the end. +/// +/// [`mean`]: MomentsU32::mean +/// [`variance`]: MomentsU32::variance +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)] +pub struct MomentsU32 { + /// Number of values. + pub n: u64, + /// `Σx`. + pub sum: u128, + /// `Σx²`. + pub sum_sq: u128, +} + +impl MomentsU32 { + /// Moments of the union of two samples — exact integer addition. + #[inline] + #[must_use] + pub fn merge(self, other: Self) -> Self { + Self { + n: self.n + other.n, + sum: self.sum + other.sum, + sum_sq: self.sum_sq + other.sum_sq, + } + } + + /// Sample mean; `0.0` for an empty sample. + pub fn mean(&self) -> f64 { + if self.n == 0 { + 0.0 + } else { + self.sum as f64 / self.n as f64 + } + } + + /// Population variance `E[(X - μ)²]`; `0.0` for an empty sample. + /// + /// Computed as `(n·Σx² − (Σx)²) / n²`. The numerator is formed exactly in + /// `u128` whenever it fits (always, for Hamming-scale data: n ≤ 2³², x ≤ + /// 2¹⁷), so there is no cancellation between two large floats; only the + /// final division rounds. Past that range it falls back to `f64`. + pub fn variance(&self) -> f64 { + if self.n == 0 { + return 0.0; + } + let n = u128::from(self.n); + match (n.checked_mul(self.sum_sq), self.sum.checked_mul(self.sum)) { + (Some(a), Some(b)) => (a - b) as f64 / (self.n as f64 * self.n as f64), + _ => { + let mean = self.mean(); + (self.sum_sq as f64 / self.n as f64 - mean * mean).max(0.0) + } + } + } +} + +/// Exact [`MomentsU32`] of `values`, eight lanes at a time through `U64x8`. +/// +/// Each value is widened to `u64` and its square split into 32-bit halves, so +/// every lane add is below 2³², and the lanes are drained into `u128` totals +/// every 2²⁸ chunks, before an 8-lane reduction could overflow `u64`. The +/// widening multiply is plain lane-wise Rust that the compiler vectorizes +/// (`vpmuludq` on x86); the accumulation and the final reduction go through +/// the polyfill, so every backend runs the same code. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::statistics::moments_u32; +/// +/// let m = moments_u32(&[1, 2, 3, 4]); +/// assert_eq!((m.n, m.sum, m.sum_sq), (4, 10, 30)); +/// assert_eq!(m.mean(), 2.5); +/// assert_eq!(m.variance(), 1.25); +/// ``` +pub fn moments_u32(values: &[u32]) -> MomentsU32 { + use crate::simd::U64x8; + + /// Chunks per drain. Every lane add is below 2^32, so after 2^28 chunks + /// a lane is below 2^60 and the 8-lane `reduce_sum` below 2^63 — the + /// reduction, not the lane, is the binding limit. + const DRAIN: usize = 1 << 28; + + let (chunks, tail) = values.as_chunks::<8>(); + let mut out = MomentsU32 { + n: values.len() as u64, + ..MomentsU32::default() + }; + for block in chunks.chunks(DRAIN) { + let (mut s, mut lo, mut hi) = (U64x8::splat(0), U64x8::splat(0), U64x8::splat(0)); + for c in block { + let x: [u64; 8] = core::array::from_fn(|i| u64::from(c[i])); + let sq: [u64; 8] = core::array::from_fn(|i| x[i] * x[i]); + s += U64x8::from_array(x); + lo += U64x8::from_array(core::array::from_fn(|i| sq[i] & 0xFFFF_FFFF)); + hi += U64x8::from_array(core::array::from_fn(|i| sq[i] >> 32)); + } + out.sum += u128::from(s.reduce_sum()); + out.sum_sq += u128::from(lo.reduce_sum()) + (u128::from(hi.reduce_sum()) << 32); + } + for &v in tail { + out.sum += u128::from(v); + out.sum_sq += u128::from(v) * u128::from(v); + } + out +} + +#[cfg(test)] +mod moments_tests { + use super::*; + + fn xorshift(n: usize, mut s: u64, mask: u32) -> Vec { + (0..n) + .map(|_| { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + (s as u32) & mask + }) + .collect() + } + + fn reference(x: &[u32]) -> MomentsU32 { + MomentsU32 { + n: x.len() as u64, + sum: x.iter().map(|&v| u128::from(v)).sum(), + sum_sq: x.iter().map(|&v| u128::from(v) * u128::from(v)).sum(), + } + } + + /// Exact against a u128 reference at every length across the 8-lane + /// chunk boundary, for small (Hamming-scale) and full-range values. + #[test] + fn moments_u32_is_exact() { + for mask in [0x3FFF, u32::MAX] { + let x = xorshift(1000, 0x9E37_79B9_7F4A_7C15, mask); + for n in (0..=40).chain([63, 64, 65, 999, 1000]) { + assert_eq!(moments_u32(&x[..n]), reference(&x[..n]), "n={n} mask={mask:#x}"); + } + } + } + + /// Worst case for the square accumulator: every square is (2^32-1)^2, + /// which alone nearly fills a u64. A lane that summed whole squares + /// would overflow on the second element. + #[test] + fn moments_u32_does_not_overflow_at_u32_max() { + let x = vec![u32::MAX; 100_003]; + let m = moments_u32(&x); + let v = u128::from(u32::MAX); + assert_eq!(m.n, 100_003); + assert_eq!(m.sum, 100_003 * v); + assert_eq!(m.sum_sq, 100_003 * v * v); + } + + /// Merging is exact integer addition: the moments of a concatenation + /// equal the merge of the parts' moments at any split point and in + /// either order — which is what makes shard-parallel statistics exact. + #[test] + fn merge_is_exact_and_order_independent() { + let x = xorshift(777, 42, u32::MAX); + let whole = moments_u32(&x); + for split in [0, 1, 7, 8, 9, 400, 776, 777] { + let (a, b) = x.split_at(split); + assert_eq!(moments_u32(a).merge(moments_u32(b)), whole, "split={split}"); + assert_eq!(moments_u32(b).merge(moments_u32(a)), whole, "reversed split={split}"); + } + } + + /// Mean and population variance against a two-pass f64 reference. + #[test] + fn mean_and_variance_match_two_pass() { + for mask in [0x3FFF, u32::MAX] { + let x = xorshift(5000, 7, mask); + let n = x.len() as f64; + let mean = x.iter().map(|&v| f64::from(v)).sum::() / n; + let var = x + .iter() + .map(|&v| (f64::from(v) - mean).powi(2)) + .sum::() + / n; + let m = moments_u32(&x); + assert!((m.mean() - mean).abs() <= 1e-9 * mean.abs(), "mean mask={mask:#x}"); + assert!((m.variance() - var).abs() <= 1e-9 * var, "variance mask={mask:#x}"); + } + assert_eq!(MomentsU32::default().mean(), 0.0); + assert_eq!(MomentsU32::default().variance(), 0.0); + assert_eq!(moments_u32(&[5, 5, 5]).variance(), 0.0, "constant input has zero variance"); + } +} diff --git a/src/hpc/zspace.rs b/src/hpc/zspace.rs new file mode 100644 index 00000000..df9d1d2e --- /dev/null +++ b/src/hpc/zspace.rs @@ -0,0 +1,558 @@ +//! The z-space entry points: every similarity or distance enters the +//! substrate as a z-score before anything thresholds, averages, bands or +//! ranks it by margin. +//! +//! Two doors, no third: +//! +//! - **Cosine-shaped** values (cosine, Pearson r, normalized dot products — +//! anything bounded in `[-1, 1]`) pay the Fisher-Z entry tax: +//! [`ZGamma`] (codes) / [`fisher_z`] (one scalar). No exceptions, including lab and +//! calibration code. A raw cosine is not a z-score: its variance shrinks +//! towards the rim, so equal steps do not mean equal evidence and a fixed +//! `cos > t` threshold means a different confidence at every `t`. +//! - **Bitpacked** values go through HDR popcount stacking on the Hamming +//! distance itself, never through a cosine reinterpretation such as +//! `1 - 2d/N`. Their null distribution is known exactly — `Binomial(N, ½)` +//! for two random N-bit vectors — so [`hamming_null_z`] needs no samples. +//! +//! The Fisher transform here is the helix convention +//! (`helix::fisher_z::Similarity::fisher_z`, mirrored by `jc::stats::fisher_2z`): +//! `z = ½·(ln(1+r) − ln(1−r))` in `ln` form, with `r` clamped to +//! `[−1+ε, 1−ε]`, `ε = `[`FISHER_CLAMP_EPS`] `= 1e-9`, so every finite input +//! gives a finite z and NaN propagates. Hyperbolic depth `2z` is +//! [`hyperbolic_depth`]. +//! +//! **z is never materialized.** Populations are encoded through a +//! [`ZGamma`] envelope straight to `i8` codes and stay there; there is no +//! z buffer and no inverse back to cosine. [`fisher_z`] exists for single +//! scalars — a report statistic, a threshold — not for arrays. + +use crate::simd::F32x16; + +/// Rim clamp for the `f64` Fisher transform — identical to +/// `helix::fisher_z::Similarity::CLAMP_EPS`, so this and helix agree bit for +/// bit. At `r = 1 − ε` the transform is `≈ 10.708`. +pub const FISHER_CLAMP_EPS: f64 = 1e-9; + +/// Largest `f32` below `1.0` (`1 − 2⁻²⁴`): the `f32` rim clamp. +/// +/// `1 − 1e-9` is not representable in `f32` — it rounds to exactly `1.0`, +/// where `ln(1 − r)` is `−∞` — so the batch path clamps to the nearest +/// representable value instead. At that rim the transform is `≈ 8.664`. +pub const FISHER_CLAMP_F32: f32 = 1.0 - f32::EPSILON / 2.0; + +/// Fisher-Z transform `z = atanh(r) = ½·(ln(1+r) − ln(1−r))`. +/// +/// `r` is clamped to `[−1+ε, 1−ε]` ([`FISHER_CLAMP_EPS`]), so the result is +/// finite for every finite input, including exact `±1` and out-of-range +/// values; NaN propagates. Under the usual bivariate-normal model `z` is +/// approximately normal with variance `1/(n−3)`, which is what makes +/// confidence bands on it mean the same thing everywhere on the scale. +/// +/// **Not bit-exact across targets.** This is the `f64` scalar for single +/// values (a report statistic, a threshold) and uses the platform `ln`; its +/// last ulp can differ between targets (measured: x86-64 glibc vs wasm32). +/// The substrate path — [`ZGamma`] codes — does not use it and is bit-exact. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::zspace::fisher_z; +/// +/// assert_eq!(fisher_z(0.0), 0.0); +/// assert!((fisher_z(0.5) - 0.5_f64.atanh()).abs() < 1e-15); +/// assert!(fisher_z(1.0).is_finite()); +/// ``` +#[inline] +pub fn fisher_z(r: f64) -> f64 { + if r.is_nan() { + return f64::NAN; + } + let s = r.clamp(-1.0 + FISHER_CLAMP_EPS, 1.0 - FISHER_CLAMP_EPS); + 0.5 * ((1.0 + s).ln() - (1.0 - s).ln()) +} + +/// Hyperbolic (Poincaré-disk) depth `2·atanh(r)` — exactly `2 ×` [`fisher_z`], +/// the "Fisher 2z" of helix and `jc`. +#[inline] +pub fn hyperbolic_depth(r: f64) -> f64 { + 2.0 * fisher_z(r) +} + +/// The z-space envelope: the metadata that makes an `i8` code mean a z-score. +/// +/// **Fisher-Z is never materialized.** There is no buffer of z values and no +/// way back to cosine. A cosine enters, is transformed in registers, is +/// normalized against this envelope, and leaves as an `i8` code. Everything +/// after that — thresholds, bands, ranking — works on codes. A threshold is +/// converted into code space once with [`ZGamma::threshold_code`]; it is never +/// converted back. +/// +/// The code is `((z − z_min) / z_range)·254 − 127`, **rounded** to the +/// nearest integer (ties to even) and saturated to the **symmetric** range +/// `[−127, 127]`. `z_min`/`z_range` travel with the codes as 8 bytes of +/// little-endian `f32` ([`ZGamma::to_le_bytes`]). +/// +/// Two deliberate departures from `bgz_tensor::fisher_z::FamilyGamma`, whose +/// byte layout this shares but whose code does not: +/// +/// - **Rounding, not truncation.** `as i8` truncates toward zero, which makes +/// the code-0 bucket twice as wide as every other and pulls every code half +/// a step toward the centre (measured: mean error `+0.49` below the centre, +/// `−0.49` above it). Rounding removes that median bias. +/// - **Symmetric range with a sentinel.** Two's-complement `i8` is asymmetric +/// (`−128..=127`); using `−128` for data would give one side an extra level. +/// Data codes stay in `−127..=127` and `−128` is reserved as +/// [`ZGamma::NAN_CODE`], so NaN can never be mistaken for a real value. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct ZGamma { + /// z of the smallest cosine the envelope covers. + pub z_min: f32, + /// `z_max − z_min`; always `> 0`. + pub z_range: f32, +} + +impl ZGamma { + /// Serialized size: two little-endian `f32`s. + pub const BYTES: usize = 8; + + /// The code for a NaN input. It is outside the data range `−127..=127`, + /// so no finite cosine ever encodes to it. + pub const NAN_CODE: i8 = i8::MIN; + + /// Fit the envelope to a population of cosines. + /// + /// Fisher-Z is monotone, so the z range is the transform of the cosine + /// range: only the cosine minimum and maximum are scanned, and no z value + /// is stored. NaN elements are ignored. An empty, all-NaN, or constant + /// population gets `z_range = 1.0` so that codes stay finite. + pub fn fit(cosines: &[f32]) -> Self { + let (mut lo, mut hi) = (f32::INFINITY, f32::NEG_INFINITY); + for &c in cosines { + if !c.is_nan() { + lo = lo.min(c); + hi = hi.max(c); + } + } + if lo > hi { + return Self { + z_min: 0.0, + z_range: 1.0, + }; + } + let z_min = fisher_z_f32(lo); + let range = fisher_z_f32(hi) - z_min; + Self { + z_min, + z_range: if range > 0.0 { range } else { 1.0 }, + } + } + + /// Encode one cosine as a z-space code. + /// + /// NaN encodes as [`ZGamma::NAN_CODE`]; values outside the envelope + /// saturate to `±127`. + #[inline] + pub fn encode(&self, cosine: f32) -> i8 { + let n = (fisher_z_f32(cosine) - self.z_min) / self.z_range; + quantize(n * 254.0 - 127.0) + } + + /// Encode a slice of cosines, sixteen lanes at a time. + /// + /// Clamp, Fisher transform and normalization are fused in `F32x16` + /// registers; the only thing written is the `i8` code. The result is + /// bit-identical to [`ZGamma::encode`] per element. + /// + /// # Panics + /// + /// Panics if `cosines.len() != codes.len()`. + pub fn encode_batch(&self, cosines: &[f32], codes: &mut [i8]) { + assert_eq!(cosines.len(), codes.len(), "ZGamma::encode_batch: length mismatch"); + let (cs, ts) = cosines.as_chunks::<16>(); + let (cd, td) = codes.as_chunks_mut::<16>(); + let one = F32x16::splat(1.0); + let half = F32x16::splat(0.5); + let z_min = F32x16::splat(self.z_min); + let z_range = F32x16::splat(self.z_range); + let scale = F32x16::splat(254.0); + let shift = F32x16::splat(127.0); + for (s, d) in cs.iter().zip(cd) { + let x = F32x16::from_array(core::array::from_fn(|i| clamp_f32(s[i]))); + let (p, m) = ((one + x).to_array(), (one - x).to_array()); + let lp = F32x16::from_array(core::array::from_fn(|i| ln_det(p[i]))); + let lm = F32x16::from_array(core::array::from_fn(|i| ln_det(m[i]))); + let z = (lp - lm) * half; + let v = (((z - z_min) / z_range) * scale - shift).to_array(); + *d = core::array::from_fn(|i| quantize(v[i])); + } + for (s, d) in ts.iter().zip(td) { + *d = self.encode(*s); + } + } + + /// Convert a cosine threshold into code space, once. + /// + /// Because the code is monotone in the cosine, `encode(c) >= t` with + /// `t = threshold_code(c₀)` holds for every `c ≥ c₀`. A cosine slightly + /// below `c₀` can share its code, so the test is inclusive at the + /// resolution of one code step. + #[inline] + pub fn threshold_code(&self, cosine: f32) -> i8 { + self.encode(cosine) + } + + /// The envelope as 8 little-endian bytes: `z_min` then `z_range`. + pub fn to_le_bytes(&self) -> [u8; Self::BYTES] { + let mut b = [0u8; Self::BYTES]; + b[..4].copy_from_slice(&self.z_min.to_le_bytes()); + b[4..].copy_from_slice(&self.z_range.to_le_bytes()); + b + } + + /// Read an envelope written by [`ZGamma::to_le_bytes`]. + pub fn from_le_bytes(b: [u8; Self::BYTES]) -> Self { + Self { + z_min: f32::from_le_bytes([b[0], b[1], b[2], b[3]]), + z_range: f32::from_le_bytes([b[4], b[5], b[6], b[7]]), + } + } +} + +/// Round to the nearest code (ties to even), saturate to the symmetric data +/// range, and map NaN to the sentinel. Shared by the scalar and batch paths so +/// they cannot drift. +#[inline] +fn quantize(v: f32) -> i8 { + if v.is_nan() { + return ZGamma::NAN_CODE; + } + v.round_ties_even().clamp(-127.0, 127.0) as i8 +} + +/// Scalar `f32` Fisher-Z with the `f32` rim clamp. Private on purpose: its +/// result only ever lives in a register on its way to a code. +/// +/// Uses [`ln_det`], not `f32::ln`: the platform libm differs in the last one +/// or two ulp between targets (measured: 624 of 8 539 grid values differ +/// between x86-64 glibc and wasm32), and a one-ulp difference at a rounding +/// boundary flips a code. +#[inline] +fn fisher_z_f32(r: f32) -> f32 { + let s = clamp_f32(r); + (ln_det(1.0 + s) - ln_det(1.0 - s)) * 0.5 +} + +/// Deterministic natural log for positive normal `f32`, bit-identical on +/// every target. +/// +/// Built only from IEEE add, sub, mul, div and bit operations, evaluated in a +/// fixed order (Rust never contracts into FMA), so every conforming target +/// produces the same bits. Algorithm: musl / FreeBSD `logf` — split off the +/// exponent `k`, reduce the mantissa to `[√½, √2)`, and evaluate a degree-8 +/// polynomial in `s = f/(2+f)`. Error is below one ulp. +/// +/// Domain: the Fisher path only ever passes `1 ± s` with `|s| ≤ 1 − 2⁻²⁴`, +/// i.e. `[2⁻²⁴, 2)` — positive and normal. NaN propagates. Zero, negative, +/// subnormal and infinite inputs are outside the domain and not handled. +#[inline] +fn ln_det(x: f32) -> f32 { + const LN2_HI: f32 = f32::from_bits(0x3f31_7180); // 6.9313812256e-01 + const LN2_LO: f32 = f32::from_bits(0x3717_f7d1); // 9.0580006145e-06 + const LG1: f32 = f32::from_bits(0x3f2a_aaaa); // 0.66666662693 + const LG2: f32 = f32::from_bits(0x3ecc_ce13); // 0.40000972152 + const LG3: f32 = f32::from_bits(0x3e91_e9ee); // 0.28498786688 + const LG4: f32 = f32::from_bits(0x3e78_9e26); // 0.24279078841 + if x.is_nan() { + return x; + } + let mut ix = x.to_bits(); + ix = ix.wrapping_add(0x3f80_0000 - 0x3f35_04f3); + let k = (ix >> 23) as i32 - 0x7f; + ix = (ix & 0x007f_ffff) + 0x3f35_04f3; + let f = f32::from_bits(ix) - 1.0; + let s = f / (2.0 + f); + let z = s * s; + let w = z * z; + let t1 = w * (LG2 + w * LG4); + let t2 = z * (LG1 + w * LG3); + let r = t2 + t1; + let hfsq = 0.5 * f * f; + let dk = k as f32; + s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI +} + +#[inline] +fn clamp_f32(r: f32) -> f32 { + // `f32::clamp` returns NaN for a NaN input, so NaN propagates. + r.clamp(-FISHER_CLAMP_F32, FISHER_CLAMP_F32) +} + +/// z-score of a Hamming distance against the random-vector null. +/// +/// Two independent uniformly random `n_bits`-bit vectors differ in +/// `Binomial(n_bits, ½)` positions: `μ = n/2`, `σ = √n/2`. The result is +/// `(d − μ)/σ`, so similar vectors score **negative** (`−3` is three sigma +/// closer than chance) with no calibration samples at all. Returns `0.0` for +/// `n_bits == 0`. +/// +/// # Example +/// +/// ``` +/// use ndarray::hpc::zspace::hamming_null_z; +/// +/// // 16 384 bits: μ = 8192, σ = 64. +/// assert_eq!(hamming_null_z(8192, 16_384), 0.0); +/// assert_eq!(hamming_null_z(8192 - 3 * 64, 16_384), -3.0); +/// ``` +#[inline] +pub fn hamming_null_z(d: u64, n_bits: u64) -> f64 { + if n_bits == 0 { + return 0.0; + } + let n = n_bits as f64; + (d as f64 - n / 2.0) / (n.sqrt() / 2.0) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn fisher_z_pins_known_values() { + assert_eq!(fisher_z(0.0), 0.0); + assert!((fisher_z(0.5) - 0.549_306_144_334_054_8).abs() < 1e-15); + // Rim value with ε = 1e-9, computed independently in f64. + assert!((fisher_z(1.0) - 10.708_206_522_644_144).abs() < 1e-9); + assert_eq!(hyperbolic_depth(0.5), 2.0 * fisher_z(0.5)); + } + + #[test] + fn fisher_z_is_odd_and_monotone() { + let grid: Vec = (-99..=99).map(|i| f64::from(i) / 100.0).collect(); + for w in grid.windows(2) { + assert!(fisher_z(w[0]) < fisher_z(w[1]), "not increasing at {}", w[0]); + } + for &r in &grid { + assert_eq!(fisher_z(-r), -fisher_z(r), "not odd at {r}"); + } + } + + #[test] + fn fisher_z_is_finite_everywhere_and_propagates_nan() { + for r in [1.0, -1.0, 2.0, -2.0, f64::INFINITY, f64::NEG_INFINITY] { + assert!(fisher_z(r).is_finite(), "r = {r}"); + } + assert_eq!(fisher_z(1.0), fisher_z(5.0), "everything past the rim clamps to it"); + assert!(fisher_z(f64::NAN).is_nan()); + } + + fn cosines() -> Vec { + let mut v: Vec = (0..200) + .map(|i| ((i * 37) % 199) as f32 / 99.0 - 1.0) + .collect(); + v[3] = 1.0; + v[20] = -1.0; + v[33] = 7.5; + v[50] = f32::NAN; + v + } + + /// The fused batch path must equal the scalar encoder bit for bit, at + /// every length across the 16-lane boundary, at the rim, and on NaN. + #[test] + fn encode_batch_is_bit_identical_to_encode() { + let src = cosines(); + let g = ZGamma::fit(&src[..150]); + for n in (0..=40).chain([199, 200]) { + let mut dst = vec![0i8; n]; + g.encode_batch(&src[..n], &mut dst); + for (i, (&c, &k)) in src[..n].iter().zip(&dst).enumerate() { + assert_eq!(k, g.encode(c), "n={n} i={i} c={c}"); + } + } + } + + /// The fitted envelope spans the full code range: the population minimum + /// lands at −127, the maximum at 127. + #[test] + fn fit_spans_the_code_range() { + let pop = [-0.2f32, 0.1, 0.4, 0.9]; + let g = ZGamma::fit(&pop); + assert_eq!(g.encode(-0.2), -127); + assert_eq!(g.encode(0.9), 127); + assert!(g.encode(0.1) > -127 && g.encode(0.4) < 127); + // Outside the envelope saturates, symmetrically, rather than wrapping. + assert_eq!(g.encode(-0.9), -127); + assert_eq!(g.encode(0.99), 127); + } + + /// Codes are spaced in z, not in cosine: two cosine steps of the same + /// size get more codes near the rim, where the evidence is stronger. + #[test] + fn codes_are_spaced_in_z_not_in_cosine() { + let g = ZGamma::fit(&[0.0, 0.99]); + let mid = i32::from(g.encode(0.10)) - i32::from(g.encode(0.00)); + let rim = i32::from(g.encode(0.99)) - i32::from(g.encode(0.89)); + assert!(rim > 2 * mid, "rim step {rim} vs centre step {mid}"); + } + + #[test] + fn fit_is_finite_on_degenerate_populations() { + for pop in [&[][..], &[f32::NAN][..], &[0.5, 0.5][..], &[1.0, 1.0][..]] { + let g = ZGamma::fit(pop); + assert!(g.z_min.is_finite() && g.z_range > 0.0, "{pop:?}"); + assert!(g.encode(0.3) >= -127); + } + } + + #[test] + fn threshold_code_is_inclusive_and_monotone() { + let src = cosines(); + let g = ZGamma::fit(&src); + let t = g.threshold_code(0.3); + for &c in src.iter().filter(|c| !c.is_nan()) { + if c >= 0.3 { + assert!(g.encode(c) >= t, "c={c}"); + } + } + } + + #[test] + fn envelope_round_trips_through_8_le_bytes() { + let g = ZGamma { + z_min: -0.75, + z_range: 3.25, + }; + let b = g.to_le_bytes(); + assert_eq!(b.len(), 8); + assert_eq!(&b[..4], &(-0.75f32).to_le_bytes()); + assert_eq!(ZGamma::from_le_bytes(b), g); + } + + /// No median bias: over z uniform on the envelope, code 0 holds the same + /// share as its neighbours and the signed error is centred on both sides. + /// Truncation toward zero fails this: code 0 gets twice its share and + /// every code is pulled half a step inward. + #[test] + fn rounding_has_no_median_bias() { + let g = ZGamma { + z_min: -1.0, + z_range: 2.0, + }; + let mut count = [0u32; 256]; + let (mut err_neg, mut n_neg, mut err_pos, mut n_pos) = (0.0f64, 0u32, 0.0f64, 0u32); + for i in 0..=20_000 { + let z = -1.0 + 2.0 * i as f32 / 20_000.0; + let k = g.encode(z.tanh()); + count[(i32::from(k) + 128) as usize] += 1; + let v = f64::from(((z - g.z_min) / g.z_range) * 254.0 - 127.0); + let e = f64::from(k) - v; + if v < -0.5 { + err_neg += e; + n_neg += 1; + } else if v > 0.5 { + err_pos += e; + n_pos += 1; + } + } + let (c0, c1) = (count[128], count[129]); + assert!(c0 * 10 <= c1 * 12, "code 0 over-full: {c0} vs neighbour {c1}"); + assert!((err_neg / f64::from(n_neg)).abs() < 0.05, "bias below centre {}", err_neg / f64::from(n_neg)); + assert!((err_pos / f64::from(n_pos)).abs() < 0.05, "bias above centre {}", err_pos / f64::from(n_pos)); + } + + /// NaN has its own code, and no finite input — however far outside the + /// envelope — can produce it. + #[test] + fn nan_sentinel_is_reserved() { + let g = ZGamma::fit(&[-0.2, 0.3]); + assert_eq!(g.encode(f32::NAN), ZGamma::NAN_CODE); + let mut out = [0i8; 1]; + g.encode_batch(&[f32::NAN], &mut out); + assert_eq!(out[0], ZGamma::NAN_CODE); + for c in [-1.0f32, -0.99, -0.5, 0.0, 0.5, 1.0, 5.0, -5.0, f32::INFINITY, f32::NEG_INFINITY] { + let k = g.encode(c); + assert!((-127..=127).contains(&k), "c={c} -> {k}"); + } + } + + #[test] + fn hamming_null_z_matches_the_binomial_null() { + assert_eq!(hamming_null_z(8192, 16_384), 0.0); + assert_eq!(hamming_null_z(8192 - 64, 16_384), -1.0); + assert_eq!(hamming_null_z(8192 + 192, 16_384), 3.0); + assert_eq!(hamming_null_z(5, 0), 0.0); + } +} + +#[cfg(test)] +mod golden { + use super::*; + + fn fnv(bytes: impl IntoIterator) -> u64 { + let mut h = 0xcbf2_9ce4_8422_2325u64; + for b in bytes { + h = (h ^ u64::from(b)).wrapping_mul(0x0100_0000_01b3); + } + h + } + + fn grid() -> Vec { + // Dense in the centre, and walking every f32 bit pattern near the rim, + // where ln(1 − r) is most sensitive to the last ulp. + let mut v: Vec = (-4096..=4096).map(|i| i as f32 / 4096.0).collect(); + let mut r = 0.999f32; + while r < 1.0 { + v.push(r); + v.push(-r); + r = f32::from_bits(r.to_bits() + 97); + } + v + } + + /// `ln_det` stays within one ulp of the correctly rounded result over + /// its whole domain `[2⁻²⁴, 2)`, sampled every 257th bit pattern. + #[test] + fn ln_det_is_within_one_ulp() { + let (lo, hi) = (2f32.powi(-24).to_bits(), 2f32.to_bits()); + let mut worst = 0u32; + let mut b = lo; + while b < hi { + let x = f32::from_bits(b); + let exact = (f64::from(x).ln()) as f32; + let got = ln_det(x); + let ulp = got.to_bits().abs_diff(exact.to_bits()); + // The sign can differ only at x = 1, where both are ±0. + if exact != 0.0 { + worst = worst.max(ulp); + } + b += 257; + } + assert!(worst <= 1, "worst error {worst} ulp"); + assert!(ln_det(f32::NAN).is_nan()); + assert_eq!(ln_det(1.0), 0.0); + } + + /// Pinned digests over a fixed grid (dense centre plus every 97th f32 + /// bit pattern near the rim). The same constants must hold on every + /// target — x86-64 (all realizations), aarch64 and wasm32. Measured + /// before `ln_det`: with libm, 624 of the 8 539 z values differed between + /// x86-64 glibc and wasm32. + pub const GOLDEN_CODES: u64 = 0xfb0b_294d_2a65_3dbb; + pub const GOLDEN_Z32: u64 = 0x2e92_a08e_fd1f_cb69; + + #[test] + fn digests_are_pinned() { + let g = grid(); + assert_eq!(g.len(), 8539); + let env = ZGamma::fit(&g); + let mut codes = vec![0i8; g.len()]; + env.encode_batch(&g, &mut codes); + assert_eq!(fnv(codes.iter().map(|&c| c as u8)), GOLDEN_CODES, "codes drifted"); + let z = fnv(g + .iter() + .flat_map(|&r| fisher_z_f32(r).to_bits().to_le_bytes())); + assert_eq!(z, GOLDEN_Z32, "z bits drifted"); + } +} diff --git a/src/simd.rs b/src/simd.rs index 172e5797..58e25f7b 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -617,6 +617,8 @@ pub use crate::hpc::fft::{wht_f32, wht_f32_new}; pub use crate::hpc::fingerprint::{ vector_config, Fingerprint, Fingerprint1K, Fingerprint2K, Fingerprint64K, VectorConfig, VectorWidth, }; +pub use crate::hpc::statistics::{moments_u32, MomentsU32}; +pub use crate::hpc::zspace::{fisher_z, hamming_null_z, hyperbolic_depth, ZGamma}; // PR-X1 — SoA carrier + const-size slice helpers, dispatched from their // respective `simd_{type}.rs` modules. The W1a consumer contract forbids