From 8c07a54d6bd4d68eaf1226a12efa3c4f8569fa0a Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 24 Sep 2026 15:02:08 +0000 Subject: [PATCH] Revert "Merge pull request #325 from AdaWorldAPI/claude/llvm-codegen-polyfill-gni3cw" This reverts commit e3789384bf50f516639acd5464d01c990e46611b, reversing changes made to 2c915389e2a4d9b9048b3546f4a8964c90956fbc. Restores the tree of 2c91538 exactly: removes hpc::zspace (ZGamma, ln_det, fisher_z, hamming_null_z), the FidelityReport z accessors, moments_u32 / MomentsU32 + Cascade::observe_batch, their simd.rs re-exports, and the wasm-simd-parity ZGamma golden check. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm --- crates/wasm-simd-parity/src/lib.rs | 38 -- src/hpc/cascade.rs | 120 ------- src/hpc/mod.rs | 3 - src/hpc/reliability.rs | 40 +-- src/hpc/statistics.rs | 200 ----------- src/hpc/zspace.rs | 558 ----------------------------- src/simd.rs | 2 - 7 files changed, 1 insertion(+), 960 deletions(-) delete mode 100644 src/hpc/zspace.rs diff --git a/crates/wasm-simd-parity/src/lib.rs b/crates/wasm-simd-parity/src/lib.rs index 247e74d5..b84c4a6c 100644 --- a/crates/wasm-simd-parity/src/lib.rs +++ b/crates/wasm-simd-parity/src/lib.rs @@ -35,9 +35,6 @@ 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 } @@ -427,38 +424,3 @@ 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 4fb1d99b..4673c10d 100644 --- a/src/hpc/cascade.rs +++ b/src/hpc/cascade.rs @@ -208,56 +208,6 @@ 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; @@ -818,76 +768,6 @@ 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 8067cf8f..074921b6 100644 --- a/src/hpc/mod.rs +++ b/src/hpc/mod.rs @@ -25,9 +25,6 @@ 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 a529ea32..56adf785 100644 --- a/src/hpc/reliability.rs +++ b/src/hpc/reliability.rs @@ -215,12 +215,6 @@ 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). @@ -252,9 +246,7 @@ 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); - /// // 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.pearson > 0.99 && r.spearman > 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]); @@ -303,42 +295,12 @@ 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 b08b0124..14d97010 100644 --- a/src/hpc/statistics.rs +++ b/src/hpc/statistics.rs @@ -360,203 +360,3 @@ 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 deleted file mode 100644 index df9d1d2e..00000000 --- a/src/hpc/zspace.rs +++ /dev/null @@ -1,558 +0,0 @@ -//! 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 58e25f7b..172e5797 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -617,8 +617,6 @@ 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