zspace: Fisher-Z entry tax without materializing z (ZGamma codes) + binomial Hamming null; batch Welford recovered from #323

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<f32> = (-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(())
+}

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<ShiftAlert> {
+        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<u32> {
+        (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;

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.

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];