diff --git a/.claude/board/LATEST_STATE.md b/.claude/board/LATEST_STATE.md index 15518c5dc..73fdeaa20 100644 --- a/.claude/board/LATEST_STATE.md +++ b/.claude/board/LATEST_STATE.md @@ -1,3 +1,13 @@ +## 2026-08-11 — lance-graph #924 (MERGED) — evaluation plan ACTIVE after a 0-of-11 spec audit + +### Current Contract Inventory — no new types (plan + audit + corrections) + +- **`.claude/plans/weather-substrate-evaluation-v1.md` is now ACTIVE** (audited 2026-08-11, §8). §3 carries **v2, post-audit** specs — every v1 spec had failed. Read §8 before citing any EV. +- **The audit's shape is the headline:** **22 of 24** KNOWN claims CONFIRMED with `file:line`, **0 of 11** EV specs SOUND. Reliable about *what is*, unreliable about *what would falsify* — banked as `E-ZERO-FOR-ELEVEN-THE-AUTHOR-CANNOT-AUDIT-HIS-OWN-FALSIFIERS-1`. **Pre-registration review by independent adversarial readers is load-bearing, not ceremony.** +- **Two `[G]` ledger corrections:** the bitboard primitives are at **`ndarray/src/bitwise.rs`**, not `src/hpc/`; and **a SPRITE moves for 2 bytes, not the field** — `morton_shift_motion_probe` legA rigid-translates a 24×24 sprite within a 256×256 field, toroidally (§12.18 corrects §12.16; whole-field advection is one `(dx,dy)` **per tile**, which is what EV-1 actually tests). +- **Status of the queue:** EV-1..EV-10 all **Queued, none RUN** — a v2 spec is *audited*, not *validated*. **EV-9 (Wave 0) needs no data** and closes the only two `[H]` rows (K-12/K-13, whose scratch reproducers were deleted). Everything else in Wave 1 waits on one `fetch.py` re-fetch; Wave 2 is the 16k×16k scale run. +- **Still open, operator calls:** D-1 noise floor · D-2 saturation window · D-3 `from_bearing` API · D-4 dormant-lane fix shape · D-5 helix CI wiring · D-6 harness-of-record. + ## 2026-08-11 — lance-graph #923 (MERGED) — the evaluation plan + `DistanceLut::circular()` ### Current Contract Inventory — one new public constructor; one new plan diff --git a/.claude/board/PR_ARC_INVENTORY.md b/.claude/board/PR_ARC_INVENTORY.md index 73d97d6e2..81027cb10 100644 --- a/.claude/board/PR_ARC_INVENTORY.md +++ b/.claude/board/PR_ARC_INVENTORY.md @@ -1,3 +1,15 @@ +## 2026-08-11 — lance-graph #924 (MERGED) — the audit that failed 11 of 11 test specs, and the plan it made ACTIVE + +- **Added.** `weather-substrate-evaluation-v1.md` **§8 audit record** + a full **§3 rewrite to v2 specs**; the plan header flips **DRAFT → ACTIVE**. Knowledge doc **§12.18**. EPIPHANIES `E-ZERO-FOR-ELEVEN-THE-AUTHOR-CANNOT-AUDIT-HIS-OWN-FALSIFIERS-1`. `AGENT_LOG` run entry. Board hygiene for #923. +- **The audit** (13 agents — 7 Sonnet source-verifiers + 6 Opus falsifier-auditors; 13/13 completed, 0 errors, ~4.42 M subagent tokens, 165 tool calls, ~20 min; run `wf_99d677e6-b45`). +- **Locked — the meta-finding.** **Pre-registration review by independent adversarial readers is load-bearing, not ceremony.** The pass split sharply by KIND of statement: **22 of 24** KNOWN claims CONFIRMED against source with `file:line` (2 PARTIAL, corrected), but **0 of 11** EV test specs SOUND (10 VACUOUS, 1 UNDERSPECIFIED). The author was reliable about **what IS** and near-uniformly unreliable about **what would FALSIFY** — one document after being corrected through §12.10–§12.17 for exactly these modes, and citing the falsifiability rule *inside the specs themselves*. **A spec's author is structurally the wrong person to find its vacuous pass routes**; same-author review had already failed measurably once this arc (§12.10 announced the right frame and shipped the wrong one). +- **Corrections to the KNOWN ledger `[G]`.** **K-23** — the bitboard primitives live at `ndarray/src/bitwise.rs` (`lib.rs:290`; fns at `:280/:286/:201`), **not** `src/hpc/bitwise.rs`. **K-21** — *"whole-field shift" OVERSTATES*: `morton_shift_motion_probe.rs` legA rigid-translates a **24×24 `Sprite` within a 256×256 field** via one `(dx,dy)` address-delta (bit-exact vs ground-truth re-render, genuinely O(1) in sprite pixels), and the shift is **toroidal** (`wrapping_add` per lane) with a separate clipped helper. §12.16's *"2 bytes moves the whole field"* is corrected by **§12.18**: 2 bytes moves one sprite/tile; whole-field advection is one `(dx,dy)` **per tile**. legB/legC verified exact. Independently re-proven in passing: K-25's absence (zero `value_offset(ValueTenant::HelixResidue)` callers anywhere; `mailbox_scan.rs:237-240` documents the Signed360 tier as *"named; wired as they land"*). +- **The 11 spec defects, by class** (each a house classic reborn): identity-comparison silence half (EV-1 — the shift rounds to 0 below 3.86 m/s, so calm tiles compared the field to ITSELF); arithmetic-tautology silence half (EV-2 — `circular == linear` for `|a−b| ≤ 128`, and a 90° sector spans 64 indices); `count > 0` implied by sampling design (EV-2 — the `elimination_rate()` defect verbatim); guard pre-satisfied by committed data (EV-3); monotone-but-dead-knob (EV-4 — a constant curve is monotone, which is exactly what a parameter that never reaches the code produces); fixture that **404s at the pinned timestep** feeding a pipeline whose anomaly step **destroys the bound-mass under test**, with Kelvin floors on a dimensionless field (EV-5); an *"exact to 1e-4"* gate over **3-sig-fig prose constants AND Monte-Carlo statistics** — the rounding incident rebuilt inside the gate (EV-6); a claim naming top-k with **no criterion touching top-k** (EV-7a); author-satisfiable anti-vacuity (EV-7a); an instrument its own citation bans, plus a 2-arm design that could only prove divisor-vs-coprime arithmetic (EV-7b); a pass criterion that **cannot fail** plus the wrong `n` (EV-8 — 200 000 sampled pairs, not the 1 038 240 grid); assertions implied by their own subject (EV-9); a **stale-fixture route where a "winter" run re-measures summer bytes** (EV-10 — `fetch_bg.py` hardcodes `t` and skips-if-exists); and a physical-constant error creating a false-fail route (EV-1 — missing `cos(lat)`, wrong by up to **2.6×** exactly where winds are strong). +- **Deferred.** EV-1..EV-10 remain **Queued** — the v2 specs are audited, none has RUN. D-1..D-6 open. The `[H]` rows (K-12/K-13) still lack a committed reproducer; **EV-9 (Wave 0) needs no data and is the unblocked next move**. Dormant-lane defect filed, not fixed. `helix` still in no CI (`[G-absence]`, D-5). +- **Docs.** Plan §3 (v2) + §8; knowledge §12.18; EPIPHANIES `E-ZERO-FOR-ELEVEN-…-1`; `AGENT_LOG` (sole-writer rule observed — agents wrote no board file; the orchestrator consolidated). + +**Confidence (2026-08-11):** merged. The audit's verify half is `[G]` (source quotes); its attack half is a **design review, not a measurement** — a v2 spec is *audited*, not *validated*, and validation is the probe run itself. The plan is ACTIVE with **zero EVs executed**. + ## 2026-08-11 — lance-graph #923 (MERGED) — the evaluation plan (KNOWN vs TO TEST), `DistanceLut::circular()`, and the two verdict inversions that reshaped both - **Added.** `.claude/plans/weather-substrate-evaluation-v1.md` — the arc's known-vs-test ledger: **§1 KNOWN K-1..K-25** each with `file:line` and an honest grade (`[G]` / `[G-session]` probe-source-committed-fixture-fetched / `[H]` measured-but-no-committed-reproducer / `[G-absence]`), **§2** the honesty split naming the only two `[H]` rows, **§3 EV-1..EV-10** in three waves, **§4** wave gating, **§5** the D-1..D-6 operator decision register. `DistanceLut::circular()` + 3 tests (`crates/helix/src/distance.rs`). Knowledge-doc §12.14–§12.17. Three EPIPHANIES entries. diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index 62af04002..b31247df6 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -1,3 +1,11 @@ +## 2026-08-11 — EV queue re-graded to v2 specs (post-audit; rows below unchanged in STATUS) + +All eleven EV rows in the block below still read **Queued** — correctly, none has +run. But their **specifications are now v2**: the 13-agent pre-registration audit +(plan §8) found **0 of 11 v1 specs sound** and every one was rewritten. Cite the +plan's §3 v2 text, never the v1 shape summarized in the row descriptions. +**EV-9 is Wave 0 and unblocked** (no data needed). + ## weather-substrate-evaluation-v1 — EV queue (PRE-REGISTERED 2026-08-11) Plan: `.claude/plans/weather-substrate-evaluation-v1.md` (known-vs-test ledger; diff --git a/.gitignore b/.gitignore index 1ce49ada0..1c0b3cb5d 100644 --- a/.gitignore +++ b/.gitignore @@ -105,3 +105,10 @@ crates/lance-graph-planner/examples/data/rosetta/* # Main was left holding the measurement receipts while the instruments that # produced them lived on another branch. Reproducibility is not a per-PR # size concern; the generators are permanent. Do not re-add such a block. + +# Weather-probe fixtures — FETCHED, never committed (evaluation plan §0). +# 36 MB of ERA5 slabs; re-fetch with probes/weather-p1/fetch.py. The manifest +# IS committed as provenance (time_index / shape / nonfinite / min / max), so a +# fresh clone can verify which bytes a recorded result was produced from. +probes/weather-p1/fixture/**/*.npy +probes/weather-p1/fixture/*.npy diff --git a/crates/helix/Cargo.toml b/crates/helix/Cargo.toml index 72146f613..ec0621b7d 100644 --- a/crates/helix/Cargo.toml +++ b/crates/helix/Cargo.toml @@ -32,3 +32,12 @@ name = "prove_residue" # Cargo.toml `exclude` so workspace-wide commands never pull it into the # deterministic main compile graph. [workspace] + +# EV-7b carries its pre-registered bars as `#[cfg(test)] mod tests`. Cargo's +# default for an [[example]] target is `test = false`, so without this stanza +# `cargo test` BUILDS the example and silently runs NONE of its 11 tests — +# measured: 0 `test tests::` lines under a plain `cargo test`. The bars would +# be decoration. (PP-13/PP-15 P0, wave wf_daae6e63-d39.) +[[example]] +name = "ev7b_comma_antimoire" +test = true diff --git a/crates/helix/examples/ev7b_comma_antimoire.rs b/crates/helix/examples/ev7b_comma_antimoire.rs new file mode 100644 index 000000000..56d38566a --- /dev/null +++ b/crates/helix/examples/ev7b_comma_antimoire.rs @@ -0,0 +1,605 @@ +//! PROBE-EV7B-COMMA-ANTIMOIRE — the v2, three-arm falsifier for D-QUANTGATE's +//! coprime requirement (`.claude/plans/weather-substrate-evaluation-v1.md` +//! §3 EV-7b v2; the §8 audit found the v1 spec NOT SOUND on two counts: a +//! self-contradicting instrument — its "or lag-k autocorrelation" fallback is +//! exactly the tool `mu_hydration_probe.rs` (see its `three_gap_distinct_count` +//! doc comment) rules out for a deterministic equidistributed walk, since +//! Pearson autocorrelation mis-measures a near-constant step as periodicity — +//! and a two-arm design that could only prove divisor-vs-coprime arithmetic, +//! never "coprime generically avoids aliasing" vs "4-over-17 is uniquely +//! special"). +//! +//! **The claim under test** is D-QUANTGATE's *coprime requirement*: a +//! quantized-layer perturbation phase generated by a stride-`s`-over-`m` +//! modular walk must use a modulus `m` INCOMMENSURATE with the surrounding +//! tile pitch (never a divisor of it, never equal to it) — or the walk +//! aliases against the tile grid, regardless of whether `gcd(s, m) = 1`. +//! This is explicitly **NOT** a claim that stride 4 / modulus 17 (K-22's +//! specific comma, `17 = 4² + 1`) is uniquely special; arm (c) below uses a +//! DIFFERENT coprime stride over the SAME incommensurate modulus, which is +//! what separates "any coprime stride over an incommensurate modulus avoids +//! aliasing" from "4 is magic". +//! +//! **Wording fix carried from the v2 spec (loosens K-22's phrasing):** the +//! comma walk (stride 4 over 17) is PERIODIC with period 17 — it is a +//! *permutation* of `Z_17`, not aperiodic in the mathematical sense. K-22's +//! "aperiodic" means *incommensurate with the 16-pitch* (17 neither divides +//! nor is a multiple of the tile pitch 16), not "has no period at all". +//! +//! **Instrument, pinned to ONE:** a 2-D DFT of the perturbation field. +//! Statistic = energy at the FIXED tile-pitch bin `(U0, 0)`, +//! `U0 = FIELD_SIZE / TILE_PITCH`, divided by the MEDIAN energy over every +//! OTHER non-DC bin — a unitless peak-to-median ratio. `U0` is derived ONLY +//! from `FIELD_SIZE` and `TILE_PITCH`, never from any arm's own +//! `(stride, modulus)` — the same fixed bin is checked identically across +//! all three arms, so nothing is cherry-picked per arm. +//! +//! **Why `FIELD_SIZE = 272` (not a smaller, cheaper power of two):** +//! `272 = 16 × 17 = lcm(TILE_PITCH, COMMA_MODULUS)`. This is the SMALLEST +//! field size that both: +//! - lets the REGULAR arm's comb land EXACTLY on `U0` (harmonic 1 of a +//! period-16 signal inside a 272-wide field is bin `272/16 = 17 = U0`); +//! - lets the COMMA/COPRIME2 arms' comb land on bins that are PROVABLY +//! DISJOINT from `U0` (harmonics of a period-17 signal inside a 272-wide +//! field sit at multiples of `272/17 = 16`; since `gcd(16,17)=1`, `17` is +//! never a nonzero integer multiple of `16` — checked directly in +//! `tests::u0_is_never_a_comma_harmonic`, not merely asserted here). +//! A smaller field (e.g. 64, not a multiple of 17) would leave the comma +//! arm's true fundamental at the NON-integer bin `64/17 ≈ 3.76`, which is +//! close enough to a plausible `U0` choice that ordinary DFT spectral +//! leakage could smear a non-trivial, hard-to-bound fraction of its energy +//! onto the neighbouring integer bin — turning "does the comma arm alias" +//! into a leakage-magnitude question this file cannot answer without +//! actually running the numbers. `FIELD_SIZE = 272` removes that ambiguity: +//! the separation below is EXACT (up to floating-point rounding), not a +//! probabilistic estimate of leakage. +//! +//! **Field construction (identical formula for all three arms — only +//! `(stride, modulus)` differ, which is what makes the comparison +//! non-cherry-picked):** +//! `addr(x,y) = y*FIELD_SIZE + x` — the same row-major address convention +//! `fire_forget_replay_probe.rs`'s `Region::addr` uses; +//! `phase_walk(addr) = (stride * (addr % modulus)) % modulus` — the exact +//! shape of `fire_forget_replay_probe.rs::comma_phase` and the comma-family +//! walks in `morton_shift_motion_probe.rs`, generalized to take +//! `(stride, modulus)` as parameters instead of hardcoding K-22's pair; +//! `field(x,y) = phase_walk(addr(x,y)) + noise(x,y)`, where `noise(x,y)` is +//! a FIXED broadband component — same salt, same address, IDENTICAL across +//! all three arms — added on top of the (now exactly commensurate, for +//! EVERY arm, given `FIELD_SIZE`'s construction) coherent comb. +//! +//! **Why the noise floor is still needed even with an exact comb +//! separation:** without it, EVERY arm's phase-only field is an exactly +//! periodic, row-invariant signal (see the LCM reasoning above), so its DFT +//! is EXACTLY zero (up to floating-point rounding noise) at every +//! non-harmonic bin — including, for the comma/coprime2 arms, at `U0` +//! itself. Both the numerator (`U0`'s energy) and the "median over non-DC +//! bins" denominator would then be dominated by floating-point ROUNDING +//! ARTIFACTS rather than a real physical quantity, making the ratio for the +//! comma/coprime2 arms an unpredictable function of numerical noise instead +//! of a meaningful measurement. Adding a REAL, fixed broadband floor +//! (`noise(x,y)`, arm-independent) makes both `U0`'s energy and the median +//! dominated by genuine, well-defined noise power instead — turning the +//! comma/coprime2 ratio into an ordinary "is this one bin an outlier +//! relative to ~74k other same-distribution bins" statistical question, +//! which is what the instrument is actually meant to measure. The regular +//! arm's peak still dominates regardless of the noise floor's presence, +//! since its coherent comb concentrates comparable total power into a +//! handful of bins instead of spreading it across the whole field (see the +//! bars' reasoning below). +//! +//! **Three arms, pre-registered:** +//! (a) REGULAR — `stride=1, modulus=TILE_PITCH(16)`. The modulus EQUALS the +//! tile pitch, so it commensurately divides `FIELD_SIZE` +//! (`272 = 17×16`): `addr % 16 = (y*272 + x) % 16 = x % 16` for every +//! row (since `272 % 16 == 0`), collapsing the coherent component to a +//! pure period-16 ramp in `x`, IDENTICAL across all 272 rows. This is +//! the worst-case "regular stride that divides the tile pitch" +//! perturbation: every tile-row's phase pattern repeats exactly, and +//! its fundamental harmonic sits EXACTLY at `(U0, 0)`. +//! (b) COMMA — `stride=4, modulus=17` (K-22's discrete Pythagorean comma). +//! `17` divides `FIELD_SIZE` too (`272 = 16×17`), so this arm is ALSO +//! exactly row-invariant (`addr % 17 = x % 17`) — but its comb sits on +//! harmonics of `272/17 = 16`, a set that PROVABLY EXCLUDES `U0 = 17` +//! (since `17` is never a multiple of `16`). Row-invariance here is a +//! feature of this specific `FIELD_SIZE`, not evidence of "genuine 2-D +//! coupling" — what matters for the claim under test is which BINS the +//! energy concentrates on, and those bins are disjoint from `U0` by +//! construction, exactly, not approximately. +//! (c) A DIFFERENT COPRIME — `stride=5, modulus=17` (same incommensurate +//! modulus as (b), a different stride; `gcd(5,17)=1` since 17 is +//! prime). Shares the same disjoint-comb guarantee as (b) — proves the +//! coprime-with-incommensurate-modulus property generalizes past +//! stride 4 specifically, which is what separates "any coprime stride +//! over 17 avoids the tile-pitch bin" from "4-over-17 is magic". +//! +//! **Bars, pre-registered BEFORE the run:** `R_HI = 150.0`, `R_LO = 15.0`, +//! with `R_HI >= 10.0 * R_LO` (150 >= 150, held as an explicit runtime +//! assertion so the pre-registration cannot silently drift). Reasoning +//! (order-of-magnitude, not a measured number — this file has NOT been +//! compiled or run; edit-only worker, orchestrator compiles/runs centrally): +//! the regular arm concentrates essentially its entire coherent power into +//! ~8 harmonic bins (out of `272² - 1 ≈ 74 000` non-DC bins total), so even +//! a crude "comparable total power, ~8 bins vs ~74 000 bins" concentration +//! argument alone predicts a peak-to-median ratio in the tens of thousands — +//! `R_HI = 150` clears with very wide margin. For the comma/coprime2 arms, +//! `U0`'s energy is (by the exact-disjoint-comb argument above) PURELY the +//! noise floor's contribution at that one bin — an ordinary single draw from +//! the same distribution the ~74 000-bin median is computed over, so its +//! ratio to that median should sit near 1 with only ordinary statistical +//! spread; `R_LO = 15` gives generous headroom against that spread being +//! larger than typically expected. If the measured ratios land outside these +//! bars, the correct response is a recorded amendment to the bars (per the +//! plan's own discipline), not a silent re-pin. +//! +//! **Silence-side fixture guard (the v1 defect this closes — it guarded only +//! the can-fire half):** the comma arm's PHASE COMPONENT ALONE (no noise — +//! see `build_phase_only`) must carry NONZERO total energy +//! (`Σ phase(x,y)² > 0`). Checked in isolation from the (arm-independent) +//! noise floor so the guard actually tests the perturbation generator, not +//! the noise floor that would trivially keep the combined field nonzero even +//! if `phase_walk` had a bug that made it identically zero. "No peak" must +//! not silently mean "no signal". +//! +//! std-only, fully deterministic — no RNG. `noise(x,y)` is a stateless hash +//! of `(x,y)`, not a seeded stream, so there is nothing to seed and no run- +//! to-run variance. Prints `EV7B …` greppable lines to stderr, one PASS/FAIL +//! line per bar. + +use std::f64::consts::PI; + +/// Field axis: FIELD_SIZE × FIELD_SIZE cells. `272 = lcm(TILE_PITCH, 17)` — +/// see the module doc's "Why FIELD_SIZE = 272" section for why this exact +/// value (not a smaller, cheaper power of two) is what makes the +/// regular-vs-coprime bin separation EXACT rather than a leakage estimate. +const FIELD_SIZE: usize = 272; +/// The surrounding tile grid's pitch — K-22's 16-pitch (matches the v2 +/// spec's own wording-fix phrase "incommensurate with the 16-pitch"). +const TILE_PITCH: usize = 16; +/// The fixed target bin: horizontal spatial frequency corresponding to +/// period = TILE_PITCH, zero vertical frequency. Derived ONLY from +/// FIELD_SIZE and TILE_PITCH — never from any arm's own generative +/// parameters — so the same bin is checked identically for all three arms. +const U0: usize = FIELD_SIZE / TILE_PITCH; // = 17 + +// Arm (a) — regular: modulus EQUALS the tile pitch (commensurate). +const REGULAR_STRIDE: u32 = 1; +const REGULAR_MODULUS: u32 = TILE_PITCH as u32; // 16 + +// Arm (b) — the comma walk (K-22): stride 4 over 17. +const COMMA_STRIDE: u32 = 4; +const COMMA_MODULUS: u32 = 17; + +// Arm (c) — a different coprime stride over the SAME incommensurate +// modulus. `gcd(5, 17) = 1` since 17 is prime. +const COPRIME2_STRIDE: u32 = 5; +const COPRIME2_MODULUS: u32 = 17; + +/// Broadband noise-floor amplitude scale — roughly matches the phase walk's +/// own value range (0..17) so neither term structurally dominates the other +/// by construction. See the module doc's noise-floor reasoning. +const NOISE_AMPLITUDE: f64 = 17.0; +/// Salt for the noise hash ("EV7BNOIS" as big-endian bytes) — arm- +/// independent, so the SAME noise realization backs every arm. +const NOISE_SALT: u64 = 0x_4556_3742_4E4F_4953; + +/// Pre-registered bars (stated BEFORE the run, per the v2 spec's discipline; +/// see the module doc for the reasoning behind these specific values). +const R_HI: f64 = 150.0; // regular arm's ratio must clear this +const R_LO: f64 = 15.0; // comma / coprime arms' ratios must stay under this + +/// SplitMix64 finalizer (stateless), same generator family as +/// `fire_forget_replay_probe.rs::mix` — used only to derive the fixed, +/// arm-independent noise floor from `(x, y)`, never as a stateful stream. +fn mix64(mut z: u64) -> u64 { + z = z.wrapping_add(0x9E37_79B9_7F4A_7C15); + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) +} + +/// Row-major address, identical convention to +/// `fire_forget_replay_probe.rs::Region::addr`. +fn addr(x: usize, y: usize) -> u64 { + (y as u64) * (FIELD_SIZE as u64) + (x as u64) +} + +/// The shared phase-walk formula: `(stride * (addr % modulus)) % modulus`. +/// Identical shape to `fire_forget_replay_probe.rs::comma_phase` and the +/// comma-family walks in `morton_shift_motion_probe.rs`, generalized to take +/// `(stride, modulus)` so the SAME formula produces all three arms — the +/// arms differ only in their two integer parameters, never in the formula +/// itself. +fn phase_walk(addr: u64, stride: u32, modulus: u32) -> u32 { + let r = (addr % modulus as u64) as u32; + stride.wrapping_mul(r) % modulus +} + +/// Fixed, arm-independent broadband noise floor at `(x, y)`, in +/// `[0, NOISE_AMPLITUDE)`. Deterministic hash, not a seeded RNG stream. +fn noise_at(x: usize, y: usize) -> f64 { + let h = mix64(addr(x, y) ^ NOISE_SALT); + ((h & 0xFFFF) as f64 / 65536.0) * NOISE_AMPLITUDE +} + +/// The coherent phase component ALONE, no noise — used by the silence-side +/// fixture guard, which must test the PERTURBATION signal in isolation from +/// the noise floor (see the module doc). +fn build_phase_only(stride: u32, modulus: u32) -> Vec { + let mut field = vec![0.0f64; FIELD_SIZE * FIELD_SIZE]; + for y in 0..FIELD_SIZE { + for x in 0..FIELD_SIZE { + field[y * FIELD_SIZE + x] = phase_walk(addr(x, y), stride, modulus) as f64; + } + } + field +} + +/// Full field = coherent phase + the shared broadband noise floor. This is +/// the input to the 2-D DFT ratio measurement. +fn build_field(stride: u32, modulus: u32) -> Vec { + let mut field = build_phase_only(stride, modulus); + for y in 0..FIELD_SIZE { + for x in 0..FIELD_SIZE { + field[y * FIELD_SIZE + x] += noise_at(x, y); + } + } + field +} + +/// Separable 2-D DFT (row pass then column pass) — mathematically identical +/// to a direct 2-D DFT, `O(FIELD_SIZE^3)` instead of `O(FIELD_SIZE^4)`. +/// Returns `energy[v * FIELD_SIZE + u] = |X(u,v)|^2` for every bin. +fn dft2d_energy(field: &[f64]) -> Vec { + let n = FIELD_SIZE; + // Row pass: row_{re,im}[y][u] = 1-D DFT of row y over x. + let mut row_re = vec![0.0f64; n * n]; + let mut row_im = vec![0.0f64; n * n]; + for y in 0..n { + for u in 0..n { + let mut re = 0.0f64; + let mut im = 0.0f64; + for x in 0..n { + let angle = -2.0 * PI * (u * x) as f64 / n as f64; + let v = field[y * n + x]; + re += v * angle.cos(); + im += v * angle.sin(); + } + row_re[y * n + u] = re; + row_im[y * n + u] = im; + } + } + // Column pass: energy[v][u] = |1-D DFT over y of row_dft[y][u]|^2. + let mut energy = vec![0.0f64; n * n]; + for u in 0..n { + for v in 0..n { + let mut re = 0.0f64; + let mut im = 0.0f64; + for y in 0..n { + let angle = -2.0 * PI * (v * y) as f64 / n as f64; + let c = angle.cos(); + let s = angle.sin(); + let rre = row_re[y * n + u]; + let rim = row_im[y * n + u]; + // Complex multiply: (rre + i*rim) * (c + i*s). + re += rre * c - rim * s; + im += rre * s + rim * c; + } + energy[v * n + u] = re * re + im * im; + } + } + energy +} + +/// The pinned statistic: energy at the fixed tile-pitch bin `(U0, 0)` +/// divided by the median energy over every OTHER non-DC bin. `(U0, 0)` +/// itself stays IN the pool the median is drawn from (it is a non-DC bin +/// like any other) — only the true DC bin `(0, 0)` is excluded. +fn tile_pitch_ratio(energy: &[f64]) -> f64 { + let n = FIELD_SIZE; + let target = energy[U0]; // v=0, u=U0 + let mut others: Vec = Vec::with_capacity(n * n - 1); + for v in 0..n { + for u in 0..n { + if u == 0 && v == 0 { + continue; // exclude DC only + } + others.push(energy[v * n + u]); + } + } + others.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let mid = others.len() / 2; + let median = if others.len() % 2 == 0 { + (others[mid - 1] + others[mid]) / 2.0 + } else { + others[mid] + }; + // Defensive floor: avoid a literal 0.0/0.0 -> NaN if every "should be + // near-zero" bin happened to land on exact fp zero. With the noise floor + // added this should never trigger in practice; kept as a guard, not + // relied on. + target / median.max(f64::MIN_POSITIVE) +} + +/// Total energy of a field: `Σ field(x,y)²`. Used by the silence-side +/// fixture guard on the coherent-phase-only field (see `build_phase_only`). +fn total_energy(field: &[f64]) -> f64 { + field.iter().map(|v| v * v).sum() +} + +fn main() { + // Pre-registration self-check: the ratio invariant must hold BEFORE any + // measurement runs, or the bars themselves are inconsistent. + assert!( + R_HI >= 10.0 * R_LO, + "pre-registration violated: R_HI ({R_HI}) must be >= 10 * R_LO ({R_LO})" + ); + + let regular_field = build_field(REGULAR_STRIDE, REGULAR_MODULUS); + let comma_field = build_field(COMMA_STRIDE, COMMA_MODULUS); + let coprime2_field = build_field(COPRIME2_STRIDE, COPRIME2_MODULUS); + + let regular_ratio = tile_pitch_ratio(&dft2d_energy(®ular_field)); + let comma_ratio = tile_pitch_ratio(&dft2d_energy(&comma_field)); + let coprime2_ratio = tile_pitch_ratio(&dft2d_energy(&coprime2_field)); + + eprintln!( + "EV7B field_size={FIELD_SIZE} tile_pitch={TILE_PITCH} u0={U0} r_hi={R_HI} r_lo={R_LO}" + ); + eprintln!( + "EV7B arm=regular stride={REGULAR_STRIDE} modulus={REGULAR_MODULUS} ratio={regular_ratio:.6e}" + ); + eprintln!( + "EV7B arm=comma stride={COMMA_STRIDE} modulus={COMMA_MODULUS} ratio={comma_ratio:.6e}" + ); + eprintln!( + "EV7B arm=coprime2 stride={COPRIME2_STRIDE} modulus={COPRIME2_MODULUS} ratio={coprime2_ratio:.6e}" + ); + + // Bar 1 — can-fire: the regular (divisor, commensurate-modulus) arm must + // alias: ratio >= R_HI. + let bar1 = regular_ratio >= R_HI; + eprintln!( + "EV7B BAR1 regular_aliases: ratio={regular_ratio:.6e} >= R_HI={R_HI} -> {}", + if bar1 { "PASS" } else { "FAIL" } + ); + + // Bar 2 — can-stay-silent: the comma arm must NOT alias: ratio <= R_LO. + let bar2 = comma_ratio <= R_LO; + eprintln!( + "EV7B BAR2 comma_silent: ratio={comma_ratio:.6e} <= R_LO={R_LO} -> {}", + if bar2 { "PASS" } else { "FAIL" } + ); + + // Bar 3 — the SAME silence under a DIFFERENT coprime stride over the + // SAME incommensurate modulus — separates "coprime + incommensurate + // modulus generically avoids aliasing" from "stride 4 specifically". + let bar3 = coprime2_ratio <= R_LO; + eprintln!( + "EV7B BAR3 coprime2_silent: ratio={coprime2_ratio:.6e} <= R_LO={R_LO} -> {}", + if bar3 { "PASS" } else { "FAIL" } + ); + + // Bar 4 — silence-side fixture guard: the comma arm's PHASE COMPONENT + // ALONE (no noise) must carry nonzero total energy. "No peak at the + // tile-pitch bin" must not silently mean "no perturbation signal at + // all" — a degenerate all-zero phase generator would trivially clear + // bar 2 for the wrong reason. + let comma_phase_total = total_energy(&build_phase_only(COMMA_STRIDE, COMMA_MODULUS)); + let bar4 = comma_phase_total > 0.0; + eprintln!( + "EV7B BAR4 comma_fixture_nonzero: phase_total_energy={comma_phase_total:.6e} -> {}", + if bar4 { "PASS" } else { "FAIL" } + ); + // Same guard for the second coprime arm too, cheaply — not the v2 + // spec's literal wording ("the comma arm's total perturbation energy"), + // but free extra rigor since it shares the guard helper. + let coprime2_phase_total = total_energy(&build_phase_only(COPRIME2_STRIDE, COPRIME2_MODULUS)); + eprintln!( + "EV7B BAR4b coprime2_fixture_nonzero (extra, not spec-required): phase_total_energy={coprime2_phase_total:.6e} -> {}", + if coprime2_phase_total > 0.0 { "PASS" } else { "FAIL" } + ); + + let overall = bar1 && bar2 && bar3 && bar4; + eprintln!("EV7B OVERALL {}", if overall { "PASS" } else { "FAIL" }); +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Pre-registration invariant: R_HI >= 10 * R_LO, checked independently + /// of `main`'s runtime assertion (so a test run alone still catches a + /// drifted bar pair without needing to execute `main`). + #[test] + fn bars_satisfy_the_ten_x_separation() { + assert!( + R_HI >= 10.0 * R_LO, + "R_HI ({R_HI}) must be >= 10 * R_LO ({R_LO})" + ); + } + + /// FIELD_SIZE must be an EXACT common multiple of TILE_PITCH and the + /// comma modulus — this is what the whole exact-disjoint-comb argument + /// rests on, checked directly rather than merely asserted in prose. + #[test] + fn field_size_is_a_common_multiple_of_tile_pitch_and_comma_modulus() { + assert_eq!(FIELD_SIZE % TILE_PITCH, 0); + assert_eq!(FIELD_SIZE % (COMMA_MODULUS as usize), 0); + assert_eq!(FIELD_SIZE % (COPRIME2_MODULUS as usize), 0); + } + + /// The mathematical core of the whole design: `U0` (the regular arm's + /// harmonic-1 bin) is NEVER an integer multiple of `FIELD_SIZE / + /// COMMA_MODULUS` (the comma arm's harmonic spacing) — i.e. the comma + /// arm's comb structurally cannot land on `U0`. Two-sided against a real + /// colliding modulus (the regular arm's own 16, see below) rather than + /// merely asserted in prose, since a hand-derived counter-example is + /// exactly the kind of claim this file's own falsifiability discipline + /// says must be exercised, not asserted. + #[test] + fn u0_is_never_a_comma_harmonic() { + let comma_harmonic_spacing = FIELD_SIZE / (COMMA_MODULUS as usize); // 16 + assert_ne!(comma_harmonic_spacing, 0); + assert_ne!( + U0 % comma_harmonic_spacing, + 0, + "U0={U0} must NOT be a multiple of the comma harmonic spacing {comma_harmonic_spacing}" + ); + // Silence-side twin: the REGULAR arm's own modulus (16) DOES + // collide — its harmonic spacing is `272/16 = 17 = U0` exactly + // (harmonic 1 IS the target bin). This proves the check above can + // distinguish the colliding case (modulus=16) from the + // non-colliding one (modulus=17), not just always report + // "disjoint" regardless of input. (17 is prime, so among small + // divisors of FIELD_SIZE only spacing=17 — i.e. modulus=16 — or the + // degenerate spacing=1 collide with U0=17; modulus=16 is the + // non-degenerate one and is also the file's own REGULAR_MODULUS.) + let regular_harmonic_spacing = FIELD_SIZE / (REGULAR_MODULUS as usize); // 17 + assert_eq!( + U0 % regular_harmonic_spacing, + 0, + "sanity: the regular arm's own modulus (16) must collide with U0" + ); + } + + /// The comma walk (K-22: stride 4 over 17) is a full permutation of + /// Z_17 — the discrete Pythagorean comma property. Two-sided by + /// construction: `seen.iter().all(|&s| s)` fails on any non-bijective + /// walk (e.g. a non-coprime stride), so this is not implied by the code + /// under test — a stride sharing a factor with 17 would break it. + #[test] + fn comma_walk_is_a_full_permutation_of_z17() { + let mut seen = [false; 17]; + for k in 0..17u64 { + let w = phase_walk(k, COMMA_STRIDE, COMMA_MODULUS) as usize; + assert!(!seen[w], "comma walk repeated residue {w} at k={k}"); + seen[w] = true; + } + assert!(seen.iter().all(|&s| s)); + } + + /// The second coprime arm (stride 5 over 17) is ALSO a full permutation + /// of Z_17 — proving `gcd(5,17)=1` holds structurally, not just + /// asserted in the doc comment. + #[test] + fn coprime2_walk_is_a_full_permutation_of_z17() { + let mut seen = [false; 17]; + for k in 0..17u64 { + let w = phase_walk(k, COPRIME2_STRIDE, COPRIME2_MODULUS) as usize; + assert!(!seen[w], "coprime2 walk repeated residue {w} at k={k}"); + seen[w] = true; + } + assert!(seen.iter().all(|&s| s)); + } + + /// The regular arm's collapse to pure-x dependency is a consequence of + /// `FIELD_SIZE % TILE_PITCH == 0`, not an assumption — checked directly, + /// which is what makes arm (a)'s expected aliasing signature + /// (constant-across-rows) a proven property of this file's own + /// constants rather than an unverified claim in prose. + #[test] + fn regular_arm_phase_depends_only_on_x() { + assert_eq!( + FIELD_SIZE % TILE_PITCH, + 0, + "field size must be a tile-pitch multiple" + ); + for y in [0usize, 1, 100, 271] { + for x in [0usize, 5, 15, 130, 271] { + let a = phase_walk(addr(x, y), REGULAR_STRIDE, REGULAR_MODULUS); + let b = phase_walk(addr(x, 0), REGULAR_STRIDE, REGULAR_MODULUS); + assert_eq!( + a, b, + "regular arm must be row-invariant at x={x} (y={y} vs y=0)" + ); + } + } + } + + /// The comma arm is ALSO row-invariant at this FIELD_SIZE (a consequence + /// of `FIELD_SIZE % COMMA_MODULUS == 0`) — the module doc's explicit + /// point that row-invariance here is a property of the chosen field + /// size, not evidence of "genuine 2-D coupling"; what actually separates + /// the arms is which BINS their (row-invariant) combs land on. + #[test] + fn comma_arm_phase_also_depends_only_on_x_at_this_field_size() { + assert_eq!(FIELD_SIZE % (COMMA_MODULUS as usize), 0); + for y in [0usize, 1, 100, 271] { + for x in [0usize, 5, 15, 130, 271] { + let a = phase_walk(addr(x, y), COMMA_STRIDE, COMMA_MODULUS); + let b = phase_walk(addr(x, 0), COMMA_STRIDE, COMMA_MODULUS); + assert_eq!( + a, b, + "comma arm must be row-invariant at x={x} (y={y} vs y=0) for FIELD_SIZE=272" + ); + } + } + } + + /// Silence-side guard, unit-level: the comma arm's phase-only field is + /// genuinely nonzero (proves `build_phase_only` + `phase_walk` are wired + /// correctly, independent of the DFT machinery). + #[test] + fn comma_phase_only_field_is_nonzero() { + let field = build_phase_only(COMMA_STRIDE, COMMA_MODULUS); + assert!(total_energy(&field) > 0.0); + // Anti-vacuity: a MAJORITY of cells must be nonzero too — a field + // that is all-zero except one stray cell would clear `> 0.0` for + // the wrong reason. + let nonzero = field.iter().filter(|&&v| v != 0.0).count(); + assert!( + nonzero * 2 > field.len(), + "comma phase field must be nonzero on a majority of cells, got {nonzero}/{}", + field.len() + ); + } + + /// Can-it-STAY-SILENT twin for the same guard: with an intentionally + /// degenerate (stride=0) walk, the phase field IS all-zero — proving the + /// guard's `> 0.0` check is not itself vacuously true for any input. + #[test] + fn zero_stride_walk_is_the_genuine_silent_case() { + let field = build_phase_only(0, COMMA_MODULUS); + assert_eq!(total_energy(&field), 0.0); + } + + /// DFT sanity: the DC bin (0,0) equals the (squared) sum of the field — + /// a basic Parseval-adjacent correctness check on `dft2d_energy` itself, + /// independent of any arm-specific reasoning. + #[test] + fn dc_bin_equals_squared_field_sum() { + let field = build_field(COMMA_STRIDE, COMMA_MODULUS); + let sum: f64 = field.iter().sum(); + let energy = dft2d_energy(&field); + let dc = energy[0]; + assert!( + (dc - sum * sum).abs() < 1e-6 * (sum * sum).abs().max(1.0), + "DC bin ({dc}) must equal (sum of field)^2 ({}), field sum={sum}", + sum * sum + ); + } + + /// The regular arm's tile-pitch-bin ratio must clear R_HI, and the + /// comma arm's must stay under R_LO — a can-fire vs can-be-silent + /// comparison computed directly, independent of the printed PASS/FAIL + /// lines in `main`. This is the most expensive test in the file + /// (O(FIELD_SIZE^3) per arm, FIELD_SIZE=272) — it duplicates `main`'s + /// own computation deliberately, so a `cargo test` run alone (without + /// ever running the example binary) still gates the pre-registered bars. + #[test] + fn regular_ratio_clears_comma_ratio_by_the_pre_registered_margin() { + let regular_ratio = + tile_pitch_ratio(&dft2d_energy(&build_field(REGULAR_STRIDE, REGULAR_MODULUS))); + let comma_ratio = + tile_pitch_ratio(&dft2d_energy(&build_field(COMMA_STRIDE, COMMA_MODULUS))); + assert!( + regular_ratio >= R_HI, + "regular arm ratio {regular_ratio:.3e} did not clear R_HI={R_HI}" + ); + assert!( + comma_ratio <= R_LO, + "comma arm ratio {comma_ratio:.3e} did not stay under R_LO={R_LO}" + ); + } +} diff --git a/crates/helix/tests/ev9_latitude_bands.rs b/crates/helix/tests/ev9_latitude_bands.rs new file mode 100644 index 000000000..ec9f662af --- /dev/null +++ b/crates/helix/tests/ev9_latitude_bands.rs @@ -0,0 +1,474 @@ +//! EV-9 (v2) — commits the orphaned K-12/K-13 measurements from +//! `.claude/plans/weather-substrate-evaluation-v1.md` §1a/§3 as falsifiable +//! tests. +//! +//! K-12/K-13 were measured in scratch tests that were deliberately DELETED — +//! graded `[H]` in the plan's KNOWN ledger precisely because no committed +//! reproducer existed. This file is that reproducer, built to the **v2** +//! spec (§3 EV-9) that survived the 2026-08-11 audit (§8): the v1 spec's two +//! bare-point assertions (the `0.99-1.02x` equal-budget ratio, the bare +//! `0.352°` mean) were struck as near-tautologies — both follow from +//! arithmetic alone (the two arms' step sizes; `step/4` for a uniform +//! quantizer's mean absolute error) — and are NOT asserted here in that bare +//! form. What IS asserted: (i) the two byte-encoding arms are genuinely +//! different code paths, not one aliased twice; (ii) the equal-budget ratio +//! against a tolerance DERIVED from the two arms' analytic step sizes; (iii) +//! the pole/equator error SPREAD (the actual K-12 finding, not two point +//! values that could drift together); (iv) the palette-`circular()` vs +//! nearest-`(n,sign)` comparison (K-13 vs K-7) on the SAME sample; (v) the +//! sampling scheme (`N`, latitude stratification, azimuth sweep) pinned +//! explicitly in code rather than left to the mechanism's reputation. +//! +//! Decode convention (matches `bearing_encode_paths.rs`'s `decode` — the +//! crate's own house convention for measuring `Signed360` angular error): +//! `y = polar>=128 ? (polar-128)/127 : -((127-polar)/127)`, +//! `r = sqrt(max(0, 1-y^2))`, `az = azimuth/65536 * TAU`, +//! `v = (r*sin(az), r*cos(az), y)`. +//! +//! NOTE: `helix` is excluded from the root workspace and appears in no CI +//! workflow (K-8), so these run only when invoked by hand in this crate. +use helix::placement::{HemispherePoint, Sign}; +use helix::residue::{ResidueEncoder, Signed360}; +use helix::PALETTE_SIZE; +use std::f64::consts::TAU; +use std::ops::RangeInclusive; + +// ── (v) sampling scheme, pinned explicitly ────────────────────────────── +// +// N, the latitude stratification, and the azimuth-sweep width are all named +// constants (not magic numbers buried in a loop), and +// `azimuth_actually_sweeps_within_each_latitude_band` below asserts the +// sweep is real rather than trusting golden-angle equidistribution by +// reputation. + +/// Total residue count, matching K-12/K-13's own provenance ("measured this +/// session"/"N=65536" throughout `bearing_encode_paths.rs` and +/// `signed360_claims.rs`). +const N: usize = 65_536; + +/// Equator band, matching K-12's own band definition ("equator (0-5°)"). +const EQUATOR_BAND_DEG: RangeInclusive = 0.0..=5.0; + +/// Pole band, matching K-12's own band definition ("pole (85-90°)"). +const POLE_BAND_DEG: RangeInclusive = 85.0..=90.0; + +/// How many azimuth targets the comparative test (iv) sweeps at elevation 0 +/// (pure-azimuth, the axis K-13 measures) — evenly spaced over the full +/// circle, every 5°. +const AZ_SAMPLES: usize = 72; + +/// A fixed, arbitrary PLACE anchor. `place` feeds only `rim` (the metric +/// carrier); it never reaches `polar`/`azimuth` (`residue.rs:182-204`), so +/// any fixed value is representative. +const PLACE: u64 = 0x1234; + +// ── shared helpers (match the crate's existing test-file conventions) ── + +/// `y` treated as `sin(latitude)`, matching `placement.rs`'s own doc +/// ("the pole (`y = 1`)"; "the equator (`y -> 0`)") — `y` IS the height +/// coordinate on the unit sphere, so `lat = asin(y)`. +fn latitude_deg(y: f64) -> f64 { + y.clamp(-1.0, 1.0).asin().to_degrees() +} + +/// Decode a `Signed360` to a unit-sphere Cartesian point. Identical to +/// `bearing_encode_paths.rs::decode` (the house convention for measuring +/// angular error against this codec). +fn decode(s: &Signed360) -> (f64, f64, f64) { + let y = if s.polar >= 128 { + (s.polar as f64 - 128.0) / 127.0 + } else { + -((127.0 - s.polar as f64) / 127.0) + }; + let r = (1.0 - y * y).max(0.0).sqrt(); + let a = s.azimuth as f64 / 65536.0 * TAU; + (r * a.sin(), r * a.cos(), y) +} + +/// Ground truth for residue `n`, `sign` — `signed_lift` per the brief: +/// `(r*sin(raw_azimuth), r*cos(raw_azimuth), y)`, raw azimuth `n*phi` +/// unbounded (matches `bearing_encode_paths.rs`'s `nearest` closure). +fn true_vec(n: usize, total: usize, sign: Sign) -> (f64, f64, f64) { + let p = HemispherePoint::signed_lift(n, total, sign); + (p.r * p.azimuth.sin(), p.r * p.azimuth.cos(), p.y) +} + +/// Angle in degrees between two (unit-ish) vectors. Identical to +/// `bearing_encode_paths.rs::ang`. +fn ang(a: (f64, f64, f64), b: (f64, f64, f64)) -> f64 { + (a.0 * b.0 + a.1 * b.1 + a.2 * b.2) + .clamp(-1.0, 1.0) + .acos() + .to_degrees() +} + +/// The naive 8-bit-over-`[-1,1]` linear arm — never implemented in the +/// crate; built inline here purely as the comparison reference (i)/(ii) +/// need. `y=-1 -> 0`, `y=0 -> 128` (rounds up), `y=1 -> 255`. +fn naive_linear8_polar(y: f64) -> u8 { + let y = y.clamp(-1.0, 1.0); + (((y + 1.0) / 2.0) * 255.0).round().clamp(0.0, 255.0) as u8 +} + +/// Inverse of `naive_linear8_polar`, for the equal-budget error +/// measurement in (ii). +fn naive_linear8_decode_y(polar: u8) -> f64 { + (polar as f64 / 255.0) * 2.0 - 1.0 +} + +/// The house arm's own decode (`polar` -> `y`), matching `residue.rs`'s +/// `encode_signed` encode direction exactly (inverse of `128 + round(y*127)` +/// / `127 - round(|y|*127)`). +fn house_decode_y(polar: u8) -> f64 { + if polar >= 128 { + (polar as f64 - 128.0) / 127.0 + } else { + -((127.0 - polar as f64) / 127.0) + } +} + +// ═════════════════════════════════════════════════════════════════════ +// (i) ANTI-VACUITY — the two byte-encoding arms are genuinely different +// code paths, not one aliased twice. +// ═════════════════════════════════════════════════════════════════════ + +/// The house 7-bit-magnitude+sign arm (the REAL `encode_signed` polar byte) +/// vs the naive 8-bit-over-`[-1,1]` linear arm must produce DIFFERENT bytes +/// on a real, non-trivial fraction of samples. +/// +/// Analytic derivation (full algebra in +/// `equal_budget_ratio_matches_the_analytic_step_sizes`'s doc comment): +/// the two arms' pre-rounding difference is `0.5*(1-y)`, ranging over +/// `(0, 0.5]` for `y in (0,1)` and concentrated near the equator (`y` +/// small) — so disagreement is EXPECTED on a real, but not dominant, +/// fraction, not "almost always" or "almost never". Hand-computed check at +/// `N=8`: 1/8 (12.5%) differ, at the smallest-`y` (equator-most) sample. +/// Bar set to 5% — comfortably below that spot-check and below the +/// continuous-approximation estimate (~16.7%, via +/// `E[0.5*(1-Y)] = 0.5*(1 - 2/3)` for the disk-equal-area `Y` density), +/// while remaining a real, non-zero claim. +/// +/// Falsifier: an implementation where the "naive" arm is accidentally +/// computed with the SAME formula as the house arm — two names for one +/// code path — produces a 0% differing fraction and fails this test +/// outright. That is exactly the `[H]`-not-`[G]` gap EV-9 exists to close: +/// a measurement that never actually built a second, independent arm. +#[test] +fn house_and_naive_polar_arms_differ_on_a_stated_minimum_fraction() { + let enc = ResidueEncoder::new(N); + let mut differing = 0usize; + for n in 0..N { + let y = HemispherePoint::lift(n, N).y; + let house = enc.encode_signed(PLACE, n, Sign::Pos).polar; + let naive = naive_linear8_polar(y); + if house != naive { + differing += 1; + } + } + let frac = differing as f64 / N as f64; + assert!( + frac >= 0.05, + "house vs naive polar bytes must differ on >=5% of {N} samples \ + (proves two distinct code paths exist); measured {:.4} ({differing}/{N})", + frac + ); +} + +// ═════════════════════════════════════════════════════════════════════ +// (ii) EQUAL-BUDGET RATIO — tolerance from the analytic step sizes, not +// an eyeballed band. +// ═════════════════════════════════════════════════════════════════════ + +/// Analytic quantization step of the house 7-bit-magnitude+sign arm: +/// `|y|` quantized into 127 levels over `[0,1]` +/// (`mag = round(|y|*127).clamp(0,127)`, `residue.rs:191`) — resolution +/// `1/127` per unit `y`. +const STEP_HOUSE: f64 = 1.0 / 127.0; + +/// Analytic quantization step of the naive 8-bit-over-`[-1,1]` arm: `y` +/// quantized into 256 levels (255 steps) over the FULL `[-1,1]` range — +/// resolution `2/255` per unit `y`. +const STEP_NAIVE: f64 = 2.0 / 255.0; + +/// For a round-to-nearest uniform quantizer with step `d`, applied to an +/// input whose density is smooth relative to `d` (N=65536 samples spread +/// over ~127-256 levels — hundreds of samples per level, the regime where +/// the dithering approximation holds), the mean ABSOLUTE quantization error +/// converges to `d/4` (the mean of a value uniformly distributed over +/// `[-d/2, d/2]`), independent of the input distribution's shape. So the +/// two arms' mean-error RATIO converges to their step-size ratio. +const ANALYTIC_RATIO: f64 = STEP_HOUSE / STEP_NAIVE; // = 255/254 ~= 1.003937 + +/// Tolerance band around `ANALYTIC_RATIO`. `N=65536` over ~127-256 +/// levels gives a well-converged mean (hundreds of samples per level), so +/// this stays tight relative to the per-latitude-band `0.99-1.02x` range +/// K-12 measured (those bands hold only ~500 samples each — far noisier +/// than the full-`N` mean measured here) while leaving real headroom +/// above pure float noise. +const RATIO_TOLERANCE: f64 = 0.03; + +/// The equal-budget ratio, WITH the analytic tolerance derived from +/// `STEP_HOUSE`/`STEP_NAIVE` — the v1's bare `"0.99-1.02x"` was struck +/// as a near-tautology (§8 EV-9: "assertion implied by its own subject", +/// both arms have ~equal step size BY CONSTRUCTION). The bar here is the +/// ANALYTIC ratio derived from the two step sizes, pinned BEFORE the run, +/// not an empirically-eyeballed band. +/// +/// Disable-run (orchestrator — the v2 spec's "ASYMMETRIC disable-run" +/// requirement): edit `crates/helix/src/residue.rs`'s `encode_signed` +/// magnitude-scale literal (`127.0` at `residue.rs:191`, in +/// `(p.y.abs() * 127.0).round().clamp(0.0, 127.0) as u8`) to something +/// else, e.g. `63.0` — this changes ONLY the house arm's real +/// behavior; `STEP_HOUSE`/`ANALYTIC_RATIO` are fixed constants computed +/// once in THIS file from the crate's documented design, never read back +/// from the running code, so the measured ratio drifts away from the +/// unmoved analytic bound and this assertion goes red. This is +/// deliberately asymmetric: a SYMMETRIC wiring bug that broke both arms' +/// effective precision by the same factor would leave the measured ratio +/// at ~1.000 and this test alone would not catch it — that is exactly why +/// (i)'s anti-vacuity check (byte-identity, not ratio-based) exists as an +/// independent tripwire. +#[test] +fn equal_budget_ratio_matches_the_analytic_step_sizes() { + let enc = ResidueEncoder::new(N); + let (mut sum_house, mut sum_naive) = (0.0f64, 0.0f64); + for n in 0..N { + let y_true = HemispherePoint::lift(n, N).y; + let house_polar = enc.encode_signed(PLACE, n, Sign::Pos).polar; + let naive_polar = naive_linear8_polar(y_true); + sum_house += (y_true - house_decode_y(house_polar)).abs(); + sum_naive += (y_true - naive_linear8_decode_y(naive_polar)).abs(); + } + let mean_house = sum_house / N as f64; + let mean_naive = sum_naive / N as f64; + let measured_ratio = mean_house / mean_naive; + assert!( + (measured_ratio - ANALYTIC_RATIO).abs() <= RATIO_TOLERANCE, + "measured ratio {measured_ratio:.6} must sit within {RATIO_TOLERANCE} \ + of the analytic ratio {ANALYTIC_RATIO:.6} \ + (= STEP_HOUSE/STEP_NAIVE = {STEP_HOUSE:.6}/{STEP_NAIVE:.6}); \ + mean_house={mean_house:.6}, mean_naive={mean_naive:.6}" + ); +} + +// ═════════════════════════════════════════════════════════════════════ +// (iii) THE SPREAD — pole/equator error RATIO, not two point values. +// ═════════════════════════════════════════════════════════════════════ + +/// Pole-band mean angular error / equator-band mean angular error, per +/// K-12 ("equator (0-5°) 0.112° mean ... pole (85-90°) 3.332° ... ~30x +/// spread"). A v1 that pinned the two means as independent point +/// assertions (`equator_mean ~= 0.112`, `pole_mean ~= 3.332`) could drift +/// TOGETHER under a shared systematic bug (e.g. a constant offset baked +/// into the decode formula) and stay green — the RATIO is the actual +/// finding under test, per the v2 spec's own framing ("independent point +/// assertions can drift together"). Bar loosened to 20x (K-12 measured +/// ~29.75x, i.e. 3.332/0.112) to leave headroom for run-to-run float +/// noise while staying far above any value a symmetric bug could produce. +/// +/// Falsifier: swapping which band is the numerator/denominator (see the +/// disable-run below), or a decode bug that makes error latitude- +/// INDEPENDENT (e.g. a constant error regardless of `y`), collapses this +/// ratio toward 1.0 and fails the `>=20x` bar. +/// +/// Disable-run (orchestrator): in this test, swap the ratio's numerator +/// and denominator (`eq_mean / pole_mean` instead of `pole_mean / +/// eq_mean`) — the ratio drops well below 1.0 (nowhere near `>=20`), +/// proving the `>=20x` bar is genuinely direction-sensitive rather than +/// something any two positive means would satisfy regardless of which +/// band is which. +#[test] +fn pole_band_error_dominates_equator_band_error_by_the_spread() { + let enc = ResidueEncoder::new(N); + let (mut eq_sum, mut eq_n) = (0.0f64, 0usize); + let (mut pole_sum, mut pole_n) = (0.0f64, 0usize); + for n in 0..N { + let p = HemispherePoint::lift(n, N); + let lat = latitude_deg(p.y); + let in_equator = EQUATOR_BAND_DEG.contains(&lat); + let in_pole = POLE_BAND_DEG.contains(&lat); + if !(in_equator || in_pole) { + continue; + } + let s = enc.encode_signed(PLACE, n, Sign::Pos); + let truth = true_vec(n, N, Sign::Pos); + let err = ang(truth, decode(&s)); + if in_equator { + eq_sum += err; + eq_n += 1; + } else { + pole_sum += err; + pole_n += 1; + } + } + assert!( + eq_n > 0 && pole_n > 0, + "the sampling scheme (N={N}) must populate both bands: eq_n={eq_n}, pole_n={pole_n}" + ); + let eq_mean = eq_sum / eq_n as f64; + let pole_mean = pole_sum / pole_n as f64; + let spread = pole_mean / eq_mean; + assert!( + spread >= 20.0, + "pole/equator mean-error spread must be >=20x (K-12 measured ~29.75x); \ + got {spread:.2}x ({pole_n} pole samples mean {pole_mean:.4} deg, \ + {eq_n} equator samples mean {eq_mean:.4} deg)" + ); +} + +// ═════════════════════════════════════════════════════════════════════ +// (iv) COMPARATIVE — palette-`circular()` beats nearest-`(n,sign)` on the +// SAME sample. +// ═════════════════════════════════════════════════════════════════════ + +/// The nearest-`(n,sign)` lattice search, matching +/// `bearing_encode_paths.rs`'s "Path A" exactly (the same O(N) search over +/// the golden-spiral lattice, picking the `(n, sign)` whose `true_vec` is +/// angularly closest to `target`, then decoding via the real +/// `encode_signed`). +fn nearest_encode(enc: &ResidueEncoder, target: (f64, f64, f64)) -> Signed360 { + let mut best = (f64::MAX, 0usize, Sign::Pos); + for n in 0..N { + for sg in [Sign::Pos, Sign::Neg] { + let v = true_vec(n, N, sg); + let e = ang(target, v); + if e < best.0 { + best = (e, n, sg); + } + } + } + enc.encode_signed(PLACE, best.1, best.2) +} + +/// The palette-`circular()` azimuth arm: quantize the true azimuth into a +/// `PALETTE_SIZE` (256, i.e. 8-bit — the domain `DistanceLut::circular()` +/// operates over) circular index, decode at the bucket's own angle — a +/// PURE 1-D azimuth quantization, elevation held at 0 to isolate the axis +/// K-13 measures. +fn palette_circular_encode(az_true: f64) -> (f64, f64, f64) { + let step = TAU / PALETTE_SIZE as f64; + let raw_idx = (az_true.rem_euclid(TAU) / step).round() as i64; + let idx = raw_idx.rem_euclid(PALETTE_SIZE as i64); + let az_decoded = idx as f64 * step; + (az_decoded.sin(), az_decoded.cos(), 0.0) +} + +/// Comparative K-13: on the SAME azimuth targets, the u8-palette +/// `circular()` arm's mean angular error must be LOWER than the +/// nearest-`(n,sign)` lattice arm's — the mechanism K-7/K-13 both name: +/// the golden-spiral lattice couples latitude and azimuth through ONE +/// index, so near the horizon (elevation 0, the axis this test sweeps) +/// its azimuth resolution is starved by the sparse near-equator lattice +/// density, while the palette arm writes azimuth as an independent field. +/// (Recorded, K-7/K-13: nearest-n ~0.972° mean, palette-circular ~0.352° +/// mean.) +/// +/// Falsifier: a palette arm that silently ignores the azimuth input (a +/// constant decode) would fail this in the OBVIOUS direction — its error +/// would be a large, target-independent average, losing badly to +/// nearest-n rather than winning; the strict `<` catches that as readily +/// as a genuine regression in the palette encode. +/// +/// Disable-run (orchestrator): in `palette_circular_encode`, zero the +/// `az_true` input before quantizing (`let raw_idx = (0.0_f64.rem_euclid(TAU) +/// / step).round() as i64;`) — the palette arm now ignores azimuth +/// entirely, its mean error rises to roughly the average angular distance +/// from a fixed point over a swept circle, and the `<` assertion goes +/// red. +#[test] +fn palette_circular_azimuth_beats_nearest_n_lattice_search() { + let enc = ResidueEncoder::new(N); + let (mut sum_palette, mut sum_nearest) = (0.0f64, 0.0f64); + for i in 0..AZ_SAMPLES { + let az_true = (i as f64 / AZ_SAMPLES as f64) * TAU; + let target = (az_true.sin(), az_true.cos(), 0.0); + + let palette_decoded = palette_circular_encode(az_true); + sum_palette += ang(target, palette_decoded); + + let nearest = nearest_encode(&enc, target); + sum_nearest += ang(target, decode(&nearest)); + } + let mean_palette = sum_palette / AZ_SAMPLES as f64; + let mean_nearest = sum_nearest / AZ_SAMPLES as f64; + assert!( + mean_palette < mean_nearest, + "palette-circular azimuth mean {mean_palette:.4} deg must be LOWER \ + than nearest-(n,sign) mean {mean_nearest:.4} deg \ + (K-7 recorded ~0.972 deg, K-13 recorded ~0.352 deg)" + ); +} + +// ═════════════════════════════════════════════════════════════════════ +// (v) SAMPLING SCHEME — pinned explicitly, not trusted by reputation. +// ═════════════════════════════════════════════════════════════════════ + +/// Azimuth is genuinely SWEPT (not degenerate/constant) within each +/// latitude band `pole_band_error_dominates_equator_band_error_by_the_spread` +/// draws from. Golden-angle stepping makes near-full coverage likely (K-3: +/// 256/256 coarse buckets over the full `N=65536` sweep), but "likely" is +/// not "pinned in code" — this asserts it directly for the ~500-sample +/// bands (iii) actually uses, rather than trusting the mechanism by +/// assumption. +/// +/// Bucketing matches `signed360_claims.rs::azimuth_spans_the_full_circle_ +/// not_merely_varies` exactly (`(a >> 8) as usize` on the REAL encoded +/// `u16` azimuth field — no float-division bucketing, no wrap-then-scale +/// edge cases) — reusing the crate's own K-3 coverage technique rather +/// than reinventing it. +/// +/// Bar (200/256 coarse buckets touched per band): a coupon-collector +/// estimate for ~500 i.i.d. uniform draws into 256 bins gives ~220 covered +/// buckets; golden-angle sequences are LOWER discrepancy than i.i.d. +/// (better equidistributed), so coverage should exceed that estimate. 200 +/// leaves margin below it while still requiring broad coverage, not a +/// handful of buckets. +/// +/// Falsifier: an off-by-something bug that fixes `n` (and hence azimuth) +/// near-constant within a band collapses `covered` toward 1 and fails the +/// `>=200` bar. +/// +/// Disable-run (orchestrator): replace `s.azimuth` below with a constant +/// (e.g. `0u16`) before bucketing — `covered` drops to 1 for both bands +/// and the assertion goes red. +#[test] +fn azimuth_actually_sweeps_within_each_latitude_band() { + let enc = ResidueEncoder::new(N); + let mut eq_buckets = [false; 256]; + let mut pole_buckets = [false; 256]; + let (mut eq_n, mut pole_n) = (0usize, 0usize); + for n in 0..N { + let p = HemispherePoint::lift(n, N); + let lat = latitude_deg(p.y); + let in_equator = EQUATOR_BAND_DEG.contains(&lat); + let in_pole = POLE_BAND_DEG.contains(&lat); + if !(in_equator || in_pole) { + continue; + } + let s = enc.encode_signed(PLACE, n, Sign::Pos); + let bucket = (s.azimuth >> 8) as usize; + if in_equator { + eq_buckets[bucket] = true; + eq_n += 1; + } else { + pole_buckets[bucket] = true; + pole_n += 1; + } + } + assert!( + eq_n > 0 && pole_n > 0, + "the sampling scheme (N={N}) must populate both bands: eq_n={eq_n}, pole_n={pole_n}" + ); + let eq_covered = eq_buckets.iter().filter(|b| **b).count(); + let pole_covered = pole_buckets.iter().filter(|b| **b).count(); + assert!( + eq_covered >= 200, + "equator band ({eq_n} samples) must sweep azimuth broadly: \ + only {eq_covered}/256 coarse buckets hit" + ); + assert!( + pole_covered >= 200, + "pole band ({pole_n} samples) must sweep azimuth broadly: \ + only {pole_covered}/256 coarse buckets hit" + ); +} diff --git a/probes/weather-p1/ev10_winter.json b/probes/weather-p1/ev10_winter.json new file mode 100644 index 000000000..acba1567b --- /dev/null +++ b/probes/weather-p1/ev10_winter.json @@ -0,0 +1,106 @@ +{ + "variable": "2m_temperature", + "run": "A (season/timestep only, variable held fixed at 2m_temperature)", + "time_index": { + "summer": 547476, + "winter": 543852 + }, + "stale_fixture_guard": { + "checks": [ + { + "season": "summer", + "key": "t547476/2m_temperature", + "declared_time_index": 547476, + "manifest_time_index": 547476, + "shape": [ + 721, + 1440 + ], + "time_index_matches": true + }, + { + "season": "winter", + "key": "t543852/2m_temperature", + "declared_time_index": 543852, + "manifest_time_index": 543852, + "shape": [ + 721, + 1440 + ], + "time_index_matches": true + } + ], + "max_abs_diff_K": 72.55035400390625 + }, + "flip_points": { + "summer": { + "linear": { + "flip_point_ci": 0.094119530779043, + "n_buckets_at_crossing": 6, + "mass_fraction_at_crossing": 0.01029589902849999, + "percentile_cross_check_ci": 0.094119530779043, + "sat_frac": 0.008476845430728927, + "interior_ci_min": 0.09411953077903945, + "interior_ci_med": 0.09411953077904034, + "interior_ci_max": 0.094119530779043, + "n_interior_points": 1029439, + "n_nonempty_interior_buckets": 254, + "apparatus_note": "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); expected apparatus behaviour (the control), not a defect -- the flip point is a single threshold, not a survival curve." + }, + "fisher_z": { + "flip_point_ci": 0.42898798181318165, + "n_buckets_at_crossing": 1, + "mass_fraction_at_crossing": 0.14222354728063052, + "percentile_cross_check_ci": 0.42898798181318165, + "sat_frac": 0.008202342425643397, + "interior_ci_min": 1.0828864116518844e-06, + "interior_ci_med": 0.005607407743958959, + "interior_ci_max": 0.42898798181318165, + "n_interior_points": 1029724, + "n_nonempty_interior_buckets": 178 + } + }, + "winter": { + "linear": { + "flip_point_ci": 0.09959033389547045, + "n_buckets_at_crossing": 3, + "mass_fraction_at_crossing": 0.012259439476852837, + "percentile_cross_check_ci": 0.09959033389547045, + "sat_frac": 0.00889678687008784, + "interior_ci_min": 0.0995903338954669, + "interior_ci_med": 0.09959033389546779, + "interior_ci_max": 0.09959033389547045, + "n_interior_points": 1029003, + "n_nonempty_interior_buckets": 254, + "apparatus_note": "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); expected apparatus behaviour (the control), not a defect -- the flip point is a single threshold, not a survival curve." + }, + "fisher_z": { + "flip_point_ci": 0.4708932274009494, + "n_buckets_at_crossing": 1, + "mass_fraction_at_crossing": 0.11613493518111118, + "percentile_cross_check_ci": 0.4708932274009494, + "sat_frac": 0.008359338881183542, + "interior_ci_min": 7.322778952101316e-07, + "interior_ci_med": 0.006317516502671872, + "interior_ci_max": 0.4708932274009494, + "n_interior_points": 1029561, + "n_nonempty_interior_buckets": 183 + } + } + }, + "season_comparison": { + "linear": { + "summer_flip_point_K": 0.094119530779043, + "winter_flip_point_K": 0.09959033389547045, + "ratio_max_over_min": 1.0581261197452292, + "within_factor_of_2": true + }, + "fisher_z": { + "summer_flip_point_K": 0.42898798181318165, + "winter_flip_point_K": 0.4708932274009494, + "ratio_max_over_min": 1.0976839617059877, + "within_factor_of_2": true + } + }, + "NO_VERDICT": false +} \ No newline at end of file diff --git a/probes/weather-p1/ev10_winter.py b/probes/weather-p1/ev10_winter.py new file mode 100644 index 000000000..25a89c5c2 --- /dev/null +++ b/probes/weather-p1/ev10_winter.py @@ -0,0 +1,266 @@ +"""EV-10 (v2, plan .claude/plans/weather-substrate-evaluation-v1.md sec.3) -- +second season, ONE factor varied, with a mandatory stale-fixture guard. + +Run A ONLY: SAME variable (2m_temperature). The factor varied is +season/timestep -- summer (t547476, the pinned 2021-06-15 12:00 UTC +timestep) vs winter (t543852, 2021-01-15 12:00Z). The v1 spec varied +season + timestep + variable simultaneously, which made any failure +unattributable (plan sec.8 attack-pass finding for EV-10); this run holds +the variable fixed so a divergence can only be read as a seasonal effect. + +MANDATORY STALE-FIXTURE GUARD: the v1 fetcher hardcoded `t` and skipped +re-fetch if the target file already existed, so a "winter" run could +silently re-measure summer bytes under a winter label. This script reads +fixture/manifest.json and asserts the time_index of every array it loads +against the timestep it declares, and additionally asserts the summer and +winter arrays are NOT byte-identical (a genuine stale-fixture symptom no +time_index check alone would catch, since a manifest entry can be correct +while the file underneath it is not). Either failure produces a loud +NO-VERDICT result -- never a silent pass on stale data. + +Compared quantity: the EV-3 closed-form 1% flip-point (sort interior +buckets by CI descending, cumulate occupancy/n_interior, report the CI at +which the cumulant first reaches 1% -- equivalently the 99th percentile of +ci[idx][interior]), computed independently here (not imported from +ev3_flip_points.py, which is owned by another agent) using the identical +definition so the cross-season comparison is apples-to-apples. Both the +linear and fisher_z arms are compared, per the v2 spec's correction of +v1's "same qualitative ordering" (undefined at floors >= 0.5 K, where both +arms read identically 0.0 per K-9). + +Pre-registered pass bar: flip-points must agree across seasons within a +factor of 2 (max/min <= 2.0). This is reported, not gated on -- the script +reports the actual ratio for both arms and lets the plan's decision +register (D-1) consume it; it does not silently redefine "agree" if the +ratio comes out above 2. +""" +import json + +import numpy as np + +VAR = "2m_temperature" +# season -> (fixture subdir, declared time_index) +SEASONS = { + "summer": ("t547476", 547476), + "winter": ("t543852", 543852), +} +FLIP_MASS = 0.01 # 1% flip-point, pre-registered (plan EV-3, inherited here) + + +def fz(s, eps=1e-9): + """Fisher-Z: arctanh(s), clipped off the +/-1 poles by eps.""" + s = np.clip(s, -1 + eps, 1 - eps) + return 0.5 * (np.log1p(s) - np.log1p(-s)) + + +def build_linear(anom): + """LINEAR path: uniform buckets in original units (K). CI is CONSTANT + across every interior bucket -- (hi-lo)/512 -- by construction.""" + lo, hi = np.percentile(anom, [0.4, 99.6]) + idx = np.clip(np.floor((anom - lo) / (hi - lo) * 256), 0, 255).astype(np.uint8) + edges = lo + np.arange(257) / 256 * (hi - lo) + return idx, edges + + +def build_fisher_z(anom): + """FISHER-Z path: uniform buckets in z, NON-uniform in K after tanh -- + CI genuinely varies bucket-to-bucket.""" + lo_q, hi_q = np.percentile(anom, [0.4, 99.6]) + scale = max(abs(lo_q), abs(hi_q)) + z = fz(anom / scale) + zlo, zhi = np.percentile(z, [0.4, 99.6]) + idxz = np.clip(np.floor((z - zlo) / (zhi - zlo) * 256), 0, 255).astype(np.uint8) + edges = np.tanh(zlo + np.arange(257) / 256 * (zhi - zlo)) * scale + return idxz, edges + + +def closed_form_flip_point(idx, edges, mass=FLIP_MASS): + """The EV-3 closed-form 1% flip-point (independently reimplemented here + to the same definition; ev3_flip_points.py is owned by another agent + and is neither imported nor edited by this script). + + Sort interior buckets (1..254) by CI descending, cumulate occupancy / + n_interior, report the CI at which the cumulant first reaches `mass`, + plus the number of (non-empty) buckets and the exact mass fraction + that determine the crossing. No grid, no step size. + + Cross-checked in-run against np.percentile(ci[idx][interior], + 100*(1-mass)) -- the two framings ("survival-function crossing" vs + "percentile of the interior population") describe the same step + function over the same population and must agree exactly. + """ + ci = np.diff(edges) / 2.0 # 256 half-widths, native units + interior = (idx > 0) & (idx < 255) + sat_frac = float((~interior).mean()) + + occ = np.bincount(idx.ravel(), minlength=256) + n_interior = int(occ[1:255].sum()) + assert n_interior > 0, "no interior points -- fixture is entirely saturated" + + buckets = [(ci[b], int(occ[b])) for b in range(1, 255) if occ[b] > 0] + buckets.sort(key=lambda t: -t[0]) # CI descending + + cum_occ = 0 + flip_point_ci = None + n_buckets_at_crossing = None + mass_fraction_at_crossing = None + for i, (b_ci, b_occ) in enumerate(buckets): + cum_occ += b_occ + frac = cum_occ / n_interior + if frac >= mass: + flip_point_ci = float(b_ci) + n_buckets_at_crossing = i + 1 + mass_fraction_at_crossing = float(frac) + break + + per_point_ci = ci[idx][interior] + percentile_cross_check = float(np.percentile(per_point_ci, 100.0 * (1.0 - mass))) + + interior_vals = ci[1:255] + return { + "flip_point_ci": flip_point_ci, + "n_buckets_at_crossing": n_buckets_at_crossing, + "mass_fraction_at_crossing": mass_fraction_at_crossing, + "percentile_cross_check_ci": percentile_cross_check, + "sat_frac": sat_frac, + "interior_ci_min": float(interior_vals.min()), + "interior_ci_med": float(np.median(interior_vals)), + "interior_ci_max": float(interior_vals.max()), + "n_interior_points": n_interior, + "n_nonempty_interior_buckets": len(buckets), + } + + +def main(): + result = { + "variable": VAR, + "run": "A (season/timestep only, variable held fixed at 2m_temperature)", + } + + # ---- MANDATORY STALE-FIXTURE GUARD ---- + manifest = json.load(open("fixture/manifest.json")) + guard = {"checks": []} + no_verdict_reason = None + arrays = {} + declared_time_index = {} + + for season, (tdir, declared_t) in SEASONS.items(): + declared_time_index[season] = declared_t + key = f"{tdir}/{VAR}" + entry = manifest.get(key) + if entry is None: + no_verdict_reason = f"manifest missing key {key!r}" + break + if entry["time_index"] != declared_t: + no_verdict_reason = ( + f"manifest time_index mismatch for {key}: declared {declared_t}, " + f"manifest says {entry['time_index']}" + ) + break + arr = np.load(f"fixture/{tdir}/{VAR}.npy").astype(np.float64) + if list(arr.shape) != entry["shape"]: + no_verdict_reason = ( + f"shape mismatch for {key}: array {list(arr.shape)} vs " + f"manifest {entry['shape']}" + ) + break + if not np.isfinite(arr).all(): + no_verdict_reason = f"{key}: nonfinite values in fixture" + break + arrays[season] = arr + guard["checks"].append( + { + "season": season, + "key": key, + "declared_time_index": declared_t, + "manifest_time_index": entry["time_index"], + "shape": entry["shape"], + "time_index_matches": True, + } + ) + print(f"guard: {key} time_index={entry['time_index']} shape={entry['shape']} OK") + + if no_verdict_reason is None: + max_abs_diff = float(np.max(np.abs(arrays["summer"] - arrays["winter"]))) + guard["max_abs_diff_K"] = max_abs_diff + print(f"guard: max_abs_diff(summer, winter) = {max_abs_diff:.6f} K") + if max_abs_diff == 0.0: + no_verdict_reason = ( + "summer and winter arrays are byte-identical (max_abs_diff == 0.0) " + "-- stale fixture (winter re-measured summer bytes); refusing to " + "report a verdict" + ) + + result["time_index"] = declared_time_index + result["stale_fixture_guard"] = guard + + if no_verdict_reason is not None: + result["NO_VERDICT"] = True + result["no_verdict_reason"] = no_verdict_reason + print(f"\nNO-VERDICT: {no_verdict_reason}") + json.dump(result, open("ev10_winter.json", "w"), indent=2) + raise SystemExit(1) + + # ---- per-season flip-points (both arms) ---- + flip = {} + for season, arr in arrays.items(): + print(f"\n{season} ({SEASONS[season][0]}, time_index={SEASONS[season][1]})") + clim = arr.mean(axis=1, keepdims=True) + anom = arr - clim + flip[season] = {} + + idx_lin, edges_lin = build_linear(anom) + r_lin = closed_form_flip_point(idx_lin, edges_lin) + r_lin["apparatus_note"] = ( + "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); " + "expected apparatus behaviour (the control), not a defect -- the " + "flip point is a single threshold, not a survival curve." + ) + flip[season]["linear"] = r_lin + print( + f" linear : flip_ci={r_lin['flip_point_ci']:.6f} " + f"n_buckets={r_lin['n_buckets_at_crossing']} " + f"mass={r_lin['mass_fraction_at_crossing']*100:.4f}% " + f"(cross-check {r_lin['percentile_cross_check_ci']:.6f}) " + f"sat={r_lin['sat_frac']*100:.3f}%" + ) + + idx_fz, edges_fz = build_fisher_z(anom) + r_fz = closed_form_flip_point(idx_fz, edges_fz) + flip[season]["fisher_z"] = r_fz + print( + f" fisher_z : flip_ci={r_fz['flip_point_ci']:.6f} " + f"n_buckets={r_fz['n_buckets_at_crossing']} " + f"mass={r_fz['mass_fraction_at_crossing']*100:.4f}% " + f"(cross-check {r_fz['percentile_cross_check_ci']:.6f}) " + f"sat={r_fz['sat_frac']*100:.3f}%" + ) + + result["flip_points"] = flip + + # ---- cross-season comparison, pre-registered factor-of-2 agreement ---- + comparisons = {} + print(f"\n{'arm':<10s} {'summer_K':>12s} {'winter_K':>12s} {'ratio':>8s} {'<=2x':>6s}") + for arm_name in ["linear", "fisher_z"]: + fp_s = flip["summer"][arm_name]["flip_point_ci"] + fp_w = flip["winter"][arm_name]["flip_point_ci"] + lo_fp, hi_fp = (fp_s, fp_w) if fp_s <= fp_w else (fp_w, fp_s) + ratio = (hi_fp / lo_fp) if lo_fp > 0 else float("inf") + within_factor_2 = bool(ratio <= 2.0) + comparisons[arm_name] = { + "summer_flip_point_K": fp_s, + "winter_flip_point_K": fp_w, + "ratio_max_over_min": ratio, + "within_factor_of_2": within_factor_2, + } + print(f"{arm_name:<10s} {fp_s:12.6f} {fp_w:12.6f} {ratio:8.4f} {str(within_factor_2):>6s}") + + result["season_comparison"] = comparisons + result["NO_VERDICT"] = False + + json.dump(result, open("ev10_winter.json", "w"), indent=2) + print("\nwrote ev10_winter.json") + + +if __name__ == "__main__": + main() diff --git a/probes/weather-p1/ev3_flip_points.json b/probes/weather-p1/ev3_flip_points.json new file mode 100644 index 000000000..1fcf17cd9 --- /dev/null +++ b/probes/weather-p1/ev3_flip_points.json @@ -0,0 +1,83 @@ +{ + "2m_temperature": { + "linear": { + "flip_point_ci": 0.094119530779043, + "n_buckets_at_crossing": 6, + "mass_fraction_at_crossing": 0.01029589902849999, + "percentile_cross_check_ci": 0.094119530779043, + "sat_frac": 0.008476845430728927, + "interior_ci_min": 0.09411953077903945, + "interior_ci_med": 0.09411953077904034, + "interior_ci_max": 0.094119530779043, + "n_interior_points": 1029439, + "n_nonempty_interior_buckets": 254, + "apparatus_note": "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); this is expected apparatus behaviour (the control), not a defect -- the flip point is a single threshold, not a survival curve." + }, + "fisher_z": { + "flip_point_ci": 0.42898798181318165, + "n_buckets_at_crossing": 1, + "mass_fraction_at_crossing": 0.14222354728063052, + "percentile_cross_check_ci": 0.42898798181318165, + "sat_frac": 0.008202342425643397, + "interior_ci_min": 1.0828864116518844e-06, + "interior_ci_med": 0.005607407743958959, + "interior_ci_max": 0.42898798181318165, + "n_interior_points": 1029724, + "n_nonempty_interior_buckets": 178 + } + }, + "2m_dewpoint_temperature": { + "linear": { + "flip_point_ci": 0.09714244377853376, + "n_buckets_at_crossing": 16, + "mass_fraction_at_crossing": 0.01021359857147576, + "percentile_cross_check_ci": 0.09714244377853376, + "sat_frac": 0.008597241485591, + "interior_ci_min": 0.09714244377853021, + "interior_ci_med": 0.09714244377853198, + "interior_ci_max": 0.09714244377853376, + "n_interior_points": 1029314, + "n_nonempty_interior_buckets": 254, + "apparatus_note": "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); this is expected apparatus behaviour (the control), not a defect -- the flip point is a single threshold, not a survival curve." + }, + "fisher_z": { + "flip_point_ci": 0.44686068352030467, + "n_buckets_at_crossing": 1, + "mass_fraction_at_crossing": 0.11313452867330083, + "percentile_cross_check_ci": 0.44686068352030467, + "sat_frac": 0.008205231930960087, + "interior_ci_min": 9.57218221131484e-07, + "interior_ci_med": 0.0052651464074298104, + "interior_ci_max": 0.44686068352030467, + "n_interior_points": 1029721, + "n_nonempty_interior_buckets": 182 + } + }, + "10m_u_component_of_wind": { + "linear": { + "flip_point_ci": 0.05318731786305797, + "n_buckets_at_crossing": 25, + "mass_fraction_at_crossing": 0.010778521064504033, + "percentile_cross_check_ci": 0.05318731786305797, + "sat_frac": 0.008551009400523963, + "interior_ci_min": 0.053187317863056194, + "interior_ci_med": 0.05318731786305797, + "interior_ci_max": 0.05318731786305797, + "n_interior_points": 1029362, + "n_nonempty_interior_buckets": 254, + "apparatus_note": "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); this is expected apparatus behaviour (the control), not a defect -- the flip point is a single threshold, not a survival curve." + }, + "fisher_z": { + "flip_point_ci": 0.2275215876242522, + "n_buckets_at_crossing": 1, + "mass_fraction_at_crossing": 0.044552594749411686, + "percentile_cross_check_ci": 0.2275215876242522, + "sat_frac": 0.0082861380798274, + "interior_ci_min": 1.654674086637442e-06, + "interior_ci_med": 0.005041045763830354, + "interior_ci_max": 0.2275215876242522, + "n_interior_points": 1029637, + "n_nonempty_interior_buckets": 181 + } + } +} \ No newline at end of file diff --git a/probes/weather-p1/ev3_flip_points.py b/probes/weather-p1/ev3_flip_points.py new file mode 100644 index 000000000..01a94dba4 --- /dev/null +++ b/probes/weather-p1/ev3_flip_points.py @@ -0,0 +1,179 @@ +"""EV-3 (v2, plan .claude/plans/weather-substrate-evaluation-v1.md sec.3) -- +floor flip-points computed EXACTLY, no grid search. + +The exceedance-vs-floor curve is the occupancy-weighted survival function of +INTERIOR bucket CI half-widths: for a floor F, the fraction of interior points +whose own bucket's CI exceeds F is a step function of F (each step at a +distinct bucket CI value, height = that bucket's occupancy / N_interior). The +1% flip-point is therefore closed-form: sort interior buckets by CI +descending, cumulate occupancy/N_interior, and read off the CI at which the +cumulant first reaches 0.01 -- equivalently the 99th percentile of +ci[idx][interior] (population = interior points only, so the two framings +agree exactly; both are reported below and cross-checked against each other). + +No grid, no step size, no undefined "sharp". This replaces the v1 grid sweep, +whose only guard was pre-satisfied by committed data (plan sec.8 attack-pass +finding for EV-3). + +APPARATUS CONTROL: the LINEAR arm's CI is (hi-lo)/512, a CONSTANT across every +interior bucket -- its "curve" is therefore a SINGLE THRESHOLD (the flip point +equals that constant, and the closed-form crossing degenerates to "however +many top-mass buckets are needed to reach 1%, all tied at the same CI"). That +is expected apparatus behaviour, not a defect, and is reported as the control. +The two-regime expectation (CI varying smoothly across buckets, producing a +genuine survival curve rather than a step at one value) applies to the +FISHER-Z arm only. + +Frame: bucket CI vs noise floor, never decoded round-trip error (plan sec.0). +Reference structure: p1_ci_vs_floor.py (read in full before editing this +file); parameterized here over VARS instead of hardcoded to one variable. +""" +import json + +import numpy as np + +VARS = ["2m_temperature", "2m_dewpoint_temperature", "10m_u_component_of_wind"] +FLIP_MASS = 0.01 # 1% flip-point, pre-registered (plan EV-3) + + +def fz(s, eps=1e-9): + """Fisher-Z: arctanh(s), clipped off the +/-1 poles by eps.""" + s = np.clip(s, -1 + eps, 1 - eps) + return 0.5 * (np.log1p(s) - np.log1p(-s)) + + +def load_anomaly(var): + """404 discipline (plan EV-3): probe availability first, never save a + fill array as a fixture. This raises loudly (FileNotFoundError) rather + than silently substituting -- the fixture set for this wave is limited + to the three variables named in the brief; a missing one is a real + blocker, not a soft skip.""" + a = np.load(f"fixture/{var}.npy").astype(np.float64) + assert np.isfinite(a).all(), f"{var}: nonfinite values in fixture" + clim = a.mean(axis=1, keepdims=True) + return a - clim + + +def build_linear(anom): + """LINEAR path: uniform buckets in original units. CI is CONSTANT + across every interior bucket -- (hi-lo)/512 -- by construction.""" + lo, hi = np.percentile(anom, [0.4, 99.6]) + idx = np.clip(np.floor((anom - lo) / (hi - lo) * 256), 0, 255).astype(np.uint8) + edges = lo + np.arange(257) / 256 * (hi - lo) + return idx, edges + + +def build_fisher_z(anom): + """FISHER-Z path: uniform buckets in z, NON-uniform in original units + after tanh -- CI genuinely varies bucket-to-bucket.""" + lo_q, hi_q = np.percentile(anom, [0.4, 99.6]) + scale = max(abs(lo_q), abs(hi_q)) + z = fz(anom / scale) + zlo, zhi = np.percentile(z, [0.4, 99.6]) + idxz = np.clip(np.floor((z - zlo) / (zhi - zlo) * 256), 0, 255).astype(np.uint8) + edges = np.tanh(zlo + np.arange(257) / 256 * (zhi - zlo)) * scale + return idxz, edges + + +def closed_form_flip_point(idx, edges, mass=FLIP_MASS): + """Closed-form 1% flip-point: sort interior buckets by CI descending, + cumulate occupancy / N_interior, report the CI at which the cumulant + first reaches `mass`, plus the number of (non-empty) buckets and the + exact mass fraction that determine the crossing. + + Cross-checked in-run against np.percentile(ci[idx][interior], 100*(1-mass)) + -- the two framings ("survival-function crossing" vs "percentile of the + interior population") must agree exactly since both describe the same + step function over the same population. + """ + ci = np.diff(edges) / 2.0 # 256 interior-relevant half-widths, in native units + interior = (idx > 0) & (idx < 255) + sat_frac = float((~interior).mean()) + + occ = np.bincount(idx.ravel(), minlength=256) + n_interior = int(occ[1:255].sum()) + assert n_interior > 0, "no interior points -- fixture is entirely saturated" + + # Only buckets that actually hold interior mass "determine" a crossing. + buckets = [(ci[b], int(occ[b])) for b in range(1, 255) if occ[b] > 0] + buckets.sort(key=lambda t: -t[0]) # CI descending + + cum_occ = 0 + flip_point_ci = None + n_buckets_at_crossing = None + mass_fraction_at_crossing = None + for i, (b_ci, b_occ) in enumerate(buckets): + cum_occ += b_occ + frac = cum_occ / n_interior + if frac >= mass: + flip_point_ci = float(b_ci) + n_buckets_at_crossing = i + 1 + mass_fraction_at_crossing = float(frac) + break + + # Cross-check: same population, same threshold, computed the "percentile" + # way instead of the "sorted cumulant" way. Must match to within one + # bucket-width's worth of population granularity. + per_point_ci = ci[idx][interior] + percentile_cross_check = float(np.percentile(per_point_ci, 100.0 * (1.0 - mass))) + + interior_vals = ci[1:255] + return { + "flip_point_ci": flip_point_ci, + "n_buckets_at_crossing": n_buckets_at_crossing, + "mass_fraction_at_crossing": mass_fraction_at_crossing, + "percentile_cross_check_ci": percentile_cross_check, + "sat_frac": sat_frac, + "interior_ci_min": float(interior_vals.min()), + "interior_ci_med": float(np.median(interior_vals)), + "interior_ci_max": float(interior_vals.max()), + "n_interior_points": n_interior, + "n_nonempty_interior_buckets": len(buckets), + } + + +def main(): + out = {} + for var in VARS: + print(f"\n{var}") + anom = load_anomaly(var) + out[var] = {} + + idx_lin, edges_lin = build_linear(anom) + r_lin = closed_form_flip_point(idx_lin, edges_lin) + r_lin["apparatus_note"] = ( + "LINEAR CI is CONSTANT across every interior bucket ((hi-lo)/512); " + "this is expected apparatus behaviour (the control), not a defect -- " + "the flip point is a single threshold, not a survival curve." + ) + out[var]["linear"] = r_lin + print( + f" linear : flip_ci={r_lin['flip_point_ci']:.6f} " + f"n_buckets={r_lin['n_buckets_at_crossing']} " + f"mass={r_lin['mass_fraction_at_crossing']*100:.4f}% " + f"(cross-check {r_lin['percentile_cross_check_ci']:.6f}) " + f"sat={r_lin['sat_frac']*100:.3f}%" + ) + + idx_fz, edges_fz = build_fisher_z(anom) + r_fz = closed_form_flip_point(idx_fz, edges_fz) + out[var]["fisher_z"] = r_fz + print( + f" fisher_z : flip_ci={r_fz['flip_point_ci']:.6f} " + f"n_buckets={r_fz['n_buckets_at_crossing']} " + f"mass={r_fz['mass_fraction_at_crossing']*100:.4f}% " + f"(cross-check {r_fz['percentile_cross_check_ci']:.6f}) " + f"sat={r_fz['sat_frac']*100:.3f}%" + ) + print( + f" fisher_z interior CI: min={r_fz['interior_ci_min']:.6f} " + f"med={r_fz['interior_ci_med']:.6f} max={r_fz['interior_ci_max']:.6f} " + f"(two-regime expectation applies HERE, not to linear)" + ) + + json.dump(out, open("ev3_flip_points.json", "w"), indent=2) + print("\nwrote ev3_flip_points.json") + + +if __name__ == "__main__": + main() diff --git a/probes/weather-p1/ev4_apparatus_sensitivity.json b/probes/weather-p1/ev4_apparatus_sensitivity.json new file mode 100644 index 000000000..ba19b1af8 --- /dev/null +++ b/probes/weather-p1/ev4_apparatus_sensitivity.json @@ -0,0 +1,160 @@ +{ + "variable": "2m_temperature", + "windows": [ + [ + 0.4, + 99.6 + ], + [ + 0.2, + 99.8 + ], + [ + 0.1, + 99.9 + ], + [ + 0.02, + 99.98 + ] + ], + "eps_sweep": { + "1e-06": { + "ci_interior_med": [ + 0.006406100961275918, + 0.012075583562084269, + 0.01814681691043951, + 0.04238015763386027 + ], + "shape": "increasing" + }, + "1e-09": { + "ci_interior_med": [ + 0.005607407743958959, + 0.002447382649405583, + 0.001722375404531462, + 0.015086534238865745 + ], + "shape": "non-monotone" + }, + "1e-12": { + "ci_interior_med": [ + 0.0049060395205593466, + 0.0004758896875554086, + 0.00015271215780732206, + 0.005145889293277861 + ], + "shape": "non-monotone" + }, + "1e-15": { + "ci_interior_med": [ + 0.004290501730232599, + 9.028732159777064e-05, + 1.3063868848917082e-05, + 0.0017167570178413527 + ], + "shape": "non-monotone" + } + }, + "method_sweep": { + "linear": { + "ci_interior_med": [ + 0.005607407743958959, + 0.002447382649405583, + 0.001722375404531462, + 0.015086534238865745 + ], + "shape": "non-monotone" + }, + "lower": { + "ci_interior_med": [ + 0.00023541037606378268, + 0.0003219115829136143, + 0.0005561433400078641, + 0.0004576132750715445 + ], + "shape": "non-monotone" + }, + "higher": { + "ci_interior_med": [ + 0.00023520082821981703, + 0.00032183911893657324, + 0.0005553902360837526, + 0.0004599002213874215 + ], + "shape": "non-monotone" + }, + "nearest": { + "ci_interior_med": [ + 0.00023523546920589666, + 0.00032186212811335935, + 0.0005557104616311648, + 0.0004597472700842431 + ], + "shape": "non-monotone" + }, + "midpoint": { + "ci_interior_med": [ + 0.0021692790139580254, + 0.002657946662043109, + 0.005306876987440923, + 0.00747809731100002 + ], + "shape": "increasing" + } + }, + "scale_controlling_tail": [ + { + "window": [ + 0.4, + 99.6 + ], + "lo_q": -22.823090004815, + "hi_q": 25.36610975405375, + "scale": 25.36610975405375, + "set_by": "HI" + }, + { + "window": [ + 0.2, + 99.8 + ], + "lo_q": -25.26567415601943, + "hi_q": 26.998929019631287, + "scale": 26.998929019631287, + "set_by": "HI" + }, + { + "window": [ + 0.1, + 99.9 + ], + "lo_q": -27.521971438471457, + "hi_q": 28.253163339297046, + "scale": 28.253163339297046, + "set_by": "HI" + }, + { + "window": [ + 0.02, + 99.98 + ], + "lo_q": -31.99861621968799, + "hi_q": 30.526803604950345, + "scale": 31.99861621968799, + "set_by": "LO" + } + ], + "verdict": { + "shape_changes_with_eps": true, + "shapes_across_eps": [ + "increasing", + "non-monotone", + "non-monotone", + "non-monotone" + ], + "method_max_over_min_ratio": 64.14320201613397, + "scale_controlling_tail_flips": true, + "conclusion": "APPARATUS-DOMINATED: the Fisher-Z interior-median CI shape is not stable under eps or percentile-interpolation choice, so the non-monotone curve must NOT be published as a data finding or fed to D-2. The scale-controlling-tail flip IS real and independent." + } +} \ No newline at end of file diff --git a/probes/weather-p1/ev4_apparatus_sensitivity.py b/probes/weather-p1/ev4_apparatus_sensitivity.py new file mode 100644 index 000000000..2b3c0a87e --- /dev/null +++ b/probes/weather-p1/ev4_apparatus_sensitivity.py @@ -0,0 +1,112 @@ +"""EV-4 P0 FOLLOW-UP — is the Fisher-Z non-monotone curve a DATA property or an +APPARATUS artifact? (PP-13 P0, wave wf_daae6e63-d39.) + +`ev4_window_sweep.py` reported `check2_fisher_z_monotone = False` on all three +variables and the orchestrator elevated that to a headline finding. PP-13 +called it a numerical artifact of `np.percentile`'s interpolation against the +`eps` arctanh clip pole. This script settles it by perturbing ONLY the +apparatus knobs and leaving the data untouched: if the curve's SHAPE moves +with `eps` or with the percentile interpolation method, the shape is not a +property of the field. + +Verdict is printed, not assumed. Emits ev4_apparatus_sensitivity.json. +""" +import json + +import numpy as np + +WINDOWS = [(0.4, 99.6), (0.2, 99.8), (0.1, 99.9), (0.02, 99.98)] +EPS_GRID = [1e-6, 1e-9, 1e-12, 1e-15] +METHODS = ["linear", "lower", "higher", "nearest", "midpoint"] +VAR = "2m_temperature" + + +def fisher_z(s, eps): + """arctanh, clipped off the +/-1 poles by eps (the shipped formula).""" + s = np.clip(s, -1 + eps, 1 - eps) + return 0.5 * (np.log1p(s) - np.log1p(-s)) + + +def ci_med(anom, w, eps, method): + """Interior-median bucket CI for the Fisher-Z arm at one window.""" + lo_q, hi_q = np.percentile(anom, [w[0], w[1]], method=method) + scale = max(abs(lo_q), abs(hi_q)) + z = fisher_z(anom / scale, eps) + zlo, zhi = np.percentile(z, [w[0], w[1]], method=method) + edges = np.tanh(zlo + np.arange(257) / 256 * (zhi - zlo)) * scale + ci = np.diff(edges) / 2.0 + return float(np.median(ci[1:255])), float(zlo), float(zhi), float(scale) + + +def shape(vals): + if all(x < y for x, y in zip(vals, vals[1:])): + return "increasing" + if all(x > y for x, y in zip(vals, vals[1:])): + return "decreasing" + return "non-monotone" + + +def main(): + a = np.load(f"fixture/{VAR}.npy").astype(np.float64) + anom = a - a.mean(axis=1, keepdims=True) + out = {"variable": VAR, "windows": [list(w) for w in WINDOWS]} + + out["eps_sweep"] = {} + shapes_eps = [] + for eps in EPS_GRID: + vals = [ci_med(anom, w, eps, "linear")[0] for w in WINDOWS] + s = shape(vals) + shapes_eps.append(s) + out["eps_sweep"][repr(eps)] = {"ci_interior_med": vals, "shape": s} + print(f" eps={eps:<8} {[f'{v:.6f}' for v in vals]} {s}") + + out["method_sweep"] = {} + shapes_m = [] + for m in METHODS: + vals = [ci_med(anom, w, 1e-9, m)[0] for w in WINDOWS] + s = shape(vals) + shapes_m.append(s) + out["method_sweep"][m] = {"ci_interior_med": vals, "shape": s} + print(f" method={m:<9} {[f'{v:.6f}' for v in vals]} {s}") + + # The one structural fact that is NOT apparatus: which tail sets `scale`. + out["scale_controlling_tail"] = [] + for w in WINDOWS: + lo_q, hi_q = np.percentile(anom, [w[0], w[1]]) + out["scale_controlling_tail"].append( + {"window": list(w), "lo_q": float(lo_q), "hi_q": float(hi_q), + "scale": float(max(abs(lo_q), abs(hi_q))), + "set_by": "HI" if abs(hi_q) >= abs(lo_q) else "LO"} + ) + flips = len({r["set_by"] for r in out["scale_controlling_tail"]}) > 1 + + # A shape that survives neither knob is not a property of the field. + eps_unstable = len(set(shapes_eps)) > 1 + method_spread = max( + max(v["ci_interior_med"]) for v in out["method_sweep"].values() + ) / max( + 1e-30, min(min(v["ci_interior_med"]) for v in out["method_sweep"].values()) + ) + out["verdict"] = { + "shape_changes_with_eps": eps_unstable, + "shapes_across_eps": shapes_eps, + "method_max_over_min_ratio": method_spread, + "scale_controlling_tail_flips": flips, + "conclusion": ( + "APPARATUS-DOMINATED: the Fisher-Z interior-median CI shape is not " + "stable under eps or percentile-interpolation choice, so the " + "non-monotone curve must NOT be published as a data finding or fed " + "to D-2. The scale-controlling-tail flip IS real and independent." + if eps_unstable or method_spread > 5 + else "STABLE: the shape survives both apparatus knobs." + ), + } + print(f"\n shape across eps: {shapes_eps}") + print(f" method max/min spread: {method_spread:.1f}x") + print(f" scale-controlling tail flips across windows: {flips}") + print(f"\n {out['verdict']['conclusion']}") + json.dump(out, open("ev4_apparatus_sensitivity.json", "w"), indent=2) + + +if __name__ == "__main__": + main() diff --git a/probes/weather-p1/ev4_window_sweep.json b/probes/weather-p1/ev4_window_sweep.json new file mode 100644 index 000000000..c2ad077c0 --- /dev/null +++ b/probes/weather-p1/ev4_window_sweep.json @@ -0,0 +1,411 @@ +{ + "2m_temperature": { + "windows": [ + { + "w_lo": 0.4, + "w_hi": 99.6 + }, + { + "w_lo": 0.2, + "w_hi": 99.8 + }, + { + "w_lo": 0.1, + "w_hi": 99.9 + }, + { + "w_lo": 0.02, + "w_hi": 99.98 + } + ], + "linear": [ + { + "lo": -22.823090004815, + "hi": 25.36610975405375, + "hi_minus_lo": 48.18919975886875, + "ci_uniform": 0.09411953077904123, + "ci_interior_min": 0.09411953077903945, + "ci_interior_med": 0.09411953077904034, + "ci_interior_max": 0.094119530779043, + "sat_frac": 0.008476845430728927, + "out_of_window_frac": 0.008000077053475111 + }, + { + "lo": -25.26567415601943, + "hi": 26.998929019631287, + "hi_minus_lo": 52.26460317565072, + "ci_uniform": 0.10207930307744206, + "ci_interior_min": 0.10207930307744206, + "ci_interior_med": 0.10207930307744206, + "ci_interior_max": 0.10207930307744562, + "sat_frac": 0.004314031437817846, + "out_of_window_frac": 0.0040010016951764526 + }, + { + "lo": -27.521971438471457, + "hi": 28.253163339297046, + "hi_minus_lo": 55.7751347777685, + "ci_uniform": 0.10893581011282905, + "ci_interior_min": 0.10893581011282727, + "ci_interior_med": 0.10893581011282905, + "ci_interior_max": 0.1089358101128326, + "sat_frac": 0.0022229927569733395, + "out_of_window_frac": 0.002001464016027123 + }, + { + "lo": -31.99861621968799, + "hi": 30.526803604950345, + "hi_minus_lo": 62.525419824638334, + "ci_uniform": 0.12211996059499697, + "ci_interior_min": 0.12211996059499342, + "ci_interior_med": 0.12211996059499697, + "ci_interior_max": 0.12211996059499697, + "sat_frac": 0.0004719525350593312, + "out_of_window_frac": 0.0004006780705809832 + } + ], + "fisher_z": [ + { + "scale": 25.36610975405375, + "zlo": -1.4708913538615156, + "zhi": 7.188807216567313, + "zhi_minus_zlo": 8.659698570428828, + "ci_interior_min": 1.0828864116518844e-06, + "ci_interior_med": 0.005607407743958959, + "ci_interior_max": 0.42898798181318165, + "sat_frac": 0.008202342425643397, + "out_of_window_frac": 0.008000077053475111 + }, + { + "scale": 26.998929019631287, + "zlo": -1.7031591377772433, + "zhi": 8.479073156978998, + "zhi_minus_zlo": 10.182232294756242, + "ci_interior_min": 1.0449270959611567e-07, + "ci_interior_med": 0.002447382649405583, + "ci_interior_max": 0.5367743605941592, + "sat_frac": 0.0040655339805825245, + "out_of_window_frac": 0.0040010016951764526 + }, + { + "scale": 28.253163339297046, + "zlo": -2.1672037968569406, + "zhi": 9.475108865384968, + "zhi_minus_zlo": 11.642312662241908, + "ci_interior_min": 1.7350950898276096e-08, + "ci_interior_med": 0.001722375404531462, + "ci_interior_max": 0.6423033847390104, + "sat_frac": 0.0020236168901217443, + "out_of_window_frac": 0.002001464016027123 + }, + { + "scale": 31.99861621968799, + "zlo": -6.833096834695244, + "zhi": 1.874539259938536, + "zhi_minus_zlo": 8.70763609463378, + "ci_interior_min": 2.7994206757142592e-06, + "ci_interior_med": 0.015086534238865745, + "ci_interior_max": 0.5440557792881052, + "sat_frac": 0.00041512559716443213, + "out_of_window_frac": 0.0004006780705809832 + } + ], + "check1_linear_strictly_increasing": true, + "check1_ci_ratio_widest_over_narrowest": 1.2974986125004242, + "check1_hilo_ratio_widest_over_narrowest": 1.2974986125004315, + "check1_abs_diff": 7.327471962526033e-15, + "check1_ratio_matches_to_1e-12": true, + "check2_fisher_z_ci_interior_med_values": [ + 0.005607407743958959, + 0.002447382649405583, + 0.001722375404531462, + 0.015086534238865745 + ], + "check2_fisher_z_monotone": false, + "check2_fisher_z_increasing": false, + "check2_fisher_z_decreasing": false + }, + "2m_dewpoint_temperature": { + "windows": [ + { + "w_lo": 0.4, + "w_hi": 99.6 + }, + { + "w_lo": 0.2, + "w_hi": 99.8 + }, + { + "w_lo": 0.1, + "w_hi": 99.9 + }, + { + "w_lo": 0.02, + "w_hi": 99.98 + } + ], + "linear": [ + { + "lo": -23.497032365798972, + "hi": 26.239898848809123, + "hi_minus_lo": 49.73693121460809, + "ci_uniform": 0.09714244377853198, + "ci_interior_min": 0.09714244377853021, + "ci_interior_med": 0.09714244377853198, + "ci_interior_max": 0.09714244377853376, + "sat_frac": 0.008597241485591, + "out_of_window_frac": 0.008000077053475111 + }, + { + "lo": -26.333363460392416, + "hi": 27.497747978994585, + "hi_minus_lo": 53.831111439387, + "ci_uniform": 0.10513888953005335, + "ci_interior_min": 0.1051388895300498, + "ci_interior_med": 0.10513888953005335, + "ci_interior_max": 0.10513888953005335, + "sat_frac": 0.004411311450146401, + "out_of_window_frac": 0.0040010016951764526 + }, + { + "lo": -28.659016241221956, + "hi": 28.24588552614851, + "hi_minus_lo": 56.90490176737046, + "ci_uniform": 0.11114238626439565, + "ci_interior_min": 0.11114238626439388, + "ci_interior_med": 0.11114238626439565, + "ci_interior_max": 0.11114238626439743, + "sat_frac": 0.002297156726768377, + "out_of_window_frac": 0.002001464016027123 + }, + { + "lo": -32.167278445362506, + "hi": 29.360957387701635, + "hi_minus_lo": 61.528235833064144, + "ci_uniform": 0.12017233561145346, + "ci_interior_min": 0.12017233561144991, + "ci_interior_med": 0.12017233561145346, + "ci_interior_max": 0.12017233561145346, + "sat_frac": 0.0005518955154877485, + "out_of_window_frac": 0.0004006780705809832 + } + ], + "fisher_z": [ + { + "scale": 26.239898848809123, + "zlo": -1.4488721028878058, + "zhi": 7.271249873419383, + "zhi_minus_zlo": 8.720121976307189, + "ci_interior_min": 9.57218221131484e-07, + "ci_interior_med": 0.0052651464074298104, + "ci_interior_max": 0.44686068352030467, + "sat_frac": 0.008205231930960087, + "out_of_window_frac": 0.008000077053475111 + }, + { + "scale": 27.497747978994585, + "zlo": -1.9168294717433199, + "zhi": 8.52758902782773, + "zhi_minus_zlo": 10.44441849957105, + "ci_interior_min": 9.937516054492335e-08, + "ci_interior_med": 0.003014843649760124, + "ci_interior_max": 0.5606382871133886, + "sat_frac": 0.004053012790876868, + "out_of_window_frac": 0.0040010016951764526 + }, + { + "scale": 28.659016241221956, + "zlo": -9.927965427970003, + "zhi": 2.4626864448846653, + "zhi_minus_zlo": 12.390651872854669, + "ci_interior_min": 7.639318511110105e-09, + "ci_interior_med": 0.0015892583035306274, + "ci_interior_max": 0.6931914108677193, + "sat_frac": 0.002046732932655263, + "out_of_window_frac": 0.002001464016027123 + }, + { + "scale": 32.167278445362506, + "zlo": -7.548888543681487, + "zhi": 1.543810901393665, + "zhi_minus_zlo": 9.092699445075151, + "ci_interior_min": 7.053151733771301e-07, + "ci_interior_med": 0.005612358975672471, + "ci_interior_max": 0.571203519286449, + "sat_frac": 0.00047676837725381415, + "out_of_window_frac": 0.0004006780705809832 + } + ], + "check1_linear_strictly_increasing": true, + "check1_ci_ratio_widest_over_narrowest": 1.237073424726947, + "check1_hilo_ratio_widest_over_narrowest": 1.2370734247269535, + "check1_abs_diff": 6.439293542825908e-15, + "check1_ratio_matches_to_1e-12": true, + "check2_fisher_z_ci_interior_med_values": [ + 0.0052651464074298104, + 0.003014843649760124, + 0.0015892583035306274, + 0.005612358975672471 + ], + "check2_fisher_z_monotone": false, + "check2_fisher_z_increasing": false, + "check2_fisher_z_decreasing": false + }, + "10m_u_component_of_wind": { + "windows": [ + { + "w_lo": 0.4, + "w_hi": 99.6 + }, + { + "w_lo": 0.2, + "w_hi": 99.8 + }, + { + "w_lo": 0.1, + "w_hi": 99.9 + }, + { + "w_lo": 0.02, + "w_hi": 99.98 + } + ], + "linear": [ + { + "lo": -14.334214259695178, + "hi": 12.897692486190296, + "hi_minus_lo": 27.231906745885475, + "ci_uniform": 0.05318731786305797, + "ci_interior_min": 0.053187317863056194, + "ci_interior_med": 0.05318731786305797, + "ci_interior_max": 0.05318731786305797, + "sat_frac": 0.008551009400523963, + "out_of_window_frac": 0.008000077053475111 + }, + { + "lo": -15.849680583599293, + "hi": 13.715110051029175, + "hi_minus_lo": 29.564790634628466, + "ci_uniform": 0.05774373170825875, + "ci_interior_min": 0.05774373170825697, + "ci_interior_med": 0.05774373170825875, + "ci_interior_max": 0.05774373170825875, + "sat_frac": 0.004289952226845431, + "out_of_window_frac": 0.0040010016951764526 + }, + { + "lo": -17.17035068847926, + "hi": 14.542062617564676, + "hi_minus_lo": 31.712413306043935, + "ci_uniform": 0.061938307238367685, + "ci_interior_min": 0.06193830723836591, + "ci_interior_med": 0.061938307238366797, + "ci_interior_max": 0.061938307238367685, + "sat_frac": 0.0021777238403451996, + "out_of_window_frac": 0.002001464016027123 + }, + { + "lo": -18.87960365753722, + "hi": 16.804778804801817, + "hi_minus_lo": 35.68438246233904, + "ci_uniform": 0.06969605949675639, + "ci_interior_min": 0.06969605949675284, + "ci_interior_med": 0.06969605949675639, + "ci_interior_max": 0.06969605949675639, + "sat_frac": 0.000470026198181538, + "out_of_window_frac": 0.0004006780705809832 + } + ], + "fisher_z": [ + { + "scale": 14.334214259695178, + "zlo": -6.656606406370108, + "zhi": 1.4710822881560268, + "zhi_minus_zlo": 8.127688694526135, + "ci_interior_min": 1.654674086637442e-06, + "ci_interior_med": 0.005041045763830354, + "ci_interior_max": 0.2275215876242522, + "sat_frac": 0.0082861380798274, + "out_of_window_frac": 0.008000077053475111 + }, + { + "scale": 15.849680583599293, + "zlo": -8.431134114806294, + "zhi": 1.314159343107991, + "zhi_minus_zlo": 9.745293457914286, + "ci_interior_min": 6.428621368570475e-08, + "ci_interior_med": 0.0009782631941273046, + "ci_interior_max": 0.3016426413237304, + "sat_frac": 0.004245646478656187, + "out_of_window_frac": 0.0040010016951764526 + }, + { + "scale": 17.17035068847926, + "zlo": -9.642741729283811, + "zhi": 1.245187761568175, + "zhi_minus_zlo": 10.887929490851986, + "ci_interior_min": 6.990033440956722e-09, + "ci_interior_med": 0.00032948091470696994, + "ci_interior_max": 0.3650480781147191, + "sat_frac": 0.002209508398828787, + "out_of_window_frac": 0.002001464016027123 + }, + { + "scale": 18.87960365753722, + "zlo": -6.923408108760464, + "zhi": 1.4224181926145996, + "zhi_minus_zlo": 8.345826301375064, + "ci_interior_min": 1.3158534368784558e-06, + "ci_interior_med": 0.004988429979905007, + "ci_interior_max": 0.30771299606339836, + "sat_frac": 0.0004305362921867776, + "out_of_window_frac": 0.0004006780705809832 + } + ], + "check1_linear_strictly_increasing": true, + "check1_ci_ratio_widest_over_narrowest": 1.3103886846899044, + "check1_hilo_ratio_widest_over_narrowest": 1.3103886846899058, + "check1_abs_diff": 1.3322676295501878e-15, + "check1_ratio_matches_to_1e-12": true, + "check2_fisher_z_ci_interior_med_values": [ + 0.005041045763830354, + 0.0009782631941273046, + 0.00032948091470696994, + 0.004988429979905007 + ], + "check2_fisher_z_monotone": false, + "check2_fisher_z_increasing": false, + "check2_fisher_z_decreasing": false + }, + "_summary": { + "vars": [ + "2m_temperature", + "2m_dewpoint_temperature", + "10m_u_component_of_wind" + ], + "windows": [ + { + "w_lo": 0.4, + "w_hi": 99.6 + }, + { + "w_lo": 0.2, + "w_hi": 99.8 + }, + { + "w_lo": 0.1, + "w_hi": 99.9 + }, + { + "w_lo": 0.02, + "w_hi": 99.98 + } + ], + "overall_check1_pass": true, + "per_var_fisher_z_monotone": { + "2m_temperature": false, + "2m_dewpoint_temperature": false, + "10m_u_component_of_wind": false + } + } +} \ No newline at end of file diff --git a/probes/weather-p1/ev4_window_sweep.py b/probes/weather-p1/ev4_window_sweep.py new file mode 100644 index 000000000..9ec32ea7f --- /dev/null +++ b/probes/weather-p1/ev4_window_sweep.py @@ -0,0 +1,191 @@ +"""EV-4 (v2) -- Saturation-window sweep with a knob-reached-the-code proof. + +Plan: .claude/plans/weather-substrate-evaluation-v1.md section 3, EV-4 (v2). +Reference frame consumed, not re-implemented: probes/weather-p1/p1_ci_vs_floor.py +(bucket CI vs floor, never decoded round-trip error). The banned metric is +documented in probes/weather-p1/p1_noise_floor.py's own header (SUPERSEDED). + +Windows swept: (0.4,99.6), (0.2,99.8), (0.1,99.9), (0.02,99.98) [percentile +pairs]. Variable: 2m_temperature (primary, per the brief); also run on +2m_dewpoint_temperature and 10m_u_component_of_wind since it is cheap and the +brief invites it. + +For each window x variable x arm (LINEAR / FISHER_Z), report BOTH quantities +the v1 spec conflated (v1 defect: "sat% counts legitimately in-range points in +the edge buckets"): + (a) bucket-0/255 occupancy ("sat_frac") -- fraction of points landing in the + two edge buckets. This over-counts relative to (b) because np.clip() + forces every truly-out-of-window point into bucket 0 or 255 ALONGSIDE + the legitimately in-range points whose value happens to fall in that + edge bucket's own [lo, lo+bucket_width) / (hi-bucket_width, hi] span. + (b) true out-of-window mass ("out_of_window_frac") -- fraction of points + whose ORIGINAL (pre-quantization) value falls outside [lo, hi]. This is + the nominal mass the percentile window excludes by construction + (0.8% -> 0.4% -> 0.2% -> 0.04% as the four windows widen). + +MANDATORY CHECKS (this is the whole point of v2 -- catching a dead knob; v1's +"monotone" pass route was satisfiable by a constant curve, which is exactly +what a window parameter that never reaches the code would produce): + + 1. LINEAR arm: ci_linear = (hi-lo)/512 (the uniform half-bucket-width, in + the original CI-vs-floor frame's units) must be STRICTLY INCREASING + across the four windows (narrowest -> widest), AND + ci_linear[widest] / ci_linear[narrowest] must equal + (hi[widest]-lo[widest]) / (hi[narrowest]-lo[narrowest]) to 1e-12. + A dead knob (the window never actually reaching np.percentile / + the quantizer) would leave ci_linear constant across windows while + (hi-lo) still visibly varies with the window in the raw percentile + output -- the ratio equality is exactly what a dead knob cannot + satisfy, because a constant ci_ratio of 1.0 would not equal a + hilo_ratio > 1.0. + + 2. FISHER_Z arm: reported, NOT gated on monotonicity. A non-monotone + interior-median-CI curve across windows is a FINDING to investigate + (both `scale` and `(zlo,zhi)` move with the window -- report all three + per window so the finding is diagnosable), never auto-assigned to + "apparatus wrong". The v1 wording ("or the apparatus is wrong") + immunized the expectation against falsification; this version reports + the monotonicity verdict as data, not as a pass/fail gate. +""" +import numpy as np +import json + +WINDOWS = [(0.4, 99.6), (0.2, 99.8), (0.1, 99.9), (0.02, 99.98)] +VARS = ['2m_temperature', '2m_dewpoint_temperature', '10m_u_component_of_wind'] + + +def fz(s, eps=1e-9): + """Fisher-Z: arctanh(s), clipped off the +/-1 poles by eps.""" + s = np.clip(s, -1 + eps, 1 - eps) + return 0.5 * (np.log1p(s) - np.log1p(-s)) + + +def load_anom(var): + a = np.load(f'fixture/{var}.npy').astype(np.float64) + assert np.isfinite(a).all(), f"{var}: nonfinite values present in fixture" + clim = a.mean(axis=1, keepdims=True) + return a - clim + + +def linear_arm(anom, w_lo, w_hi): + """LINEAR: uniform buckets in original units (K or m/s).""" + lo, hi = np.percentile(anom, [w_lo, w_hi]) + idx = np.clip(np.floor((anom - lo) / (hi - lo) * 256), 0, 255).astype(np.uint8) + edges = lo + np.arange(257) / 256 * (hi - lo) + ci = np.diff(edges) / 2.0 # 256 half-widths; uniform arm -> all equal + interior = (idx > 0) & (idx < 255) + sat_frac = float((~interior).mean()) + out_of_window_frac = float(((anom < lo) | (anom > hi)).mean()) + return { + 'lo': float(lo), 'hi': float(hi), 'hi_minus_lo': float(hi - lo), + 'ci_uniform': float(ci[0]), # every bucket has the same CI in this arm + 'ci_interior_min': float(ci[1:255].min()), + 'ci_interior_med': float(np.median(ci[1:255])), + 'ci_interior_max': float(ci[1:255].max()), + 'sat_frac': sat_frac, + 'out_of_window_frac': out_of_window_frac, + } + + +def fisher_z_arm(anom, w_lo, w_hi): + """FISHER_Z: uniform buckets in z=arctanh(anom/scale), non-uniform in K.""" + lo_q, hi_q = np.percentile(anom, [w_lo, w_hi]) + scale = max(abs(lo_q), abs(hi_q)) + z = fz(anom / scale) + zlo, zhi = np.percentile(z, [w_lo, w_hi]) + idxz = np.clip(np.floor((z - zlo) / (zhi - zlo) * 256), 0, 255).astype(np.uint8) + edges_z = zlo + np.arange(257) / 256 * (zhi - zlo) + edges_K = np.tanh(edges_z) * scale + ci = np.diff(edges_K) / 2.0 # 256 half-widths; non-uniform (Fisher-Z arm) + interior = (idxz > 0) & (idxz < 255) + sat_frac = float((~interior).mean()) + out_of_window_frac = float(((z < zlo) | (z > zhi)).mean()) + return { + 'scale': float(scale), 'zlo': float(zlo), 'zhi': float(zhi), + 'zhi_minus_zlo': float(zhi - zlo), + 'ci_interior_min': float(ci[1:255].min()), + 'ci_interior_med': float(np.median(ci[1:255])), + 'ci_interior_max': float(ci[1:255].max()), + 'sat_frac': sat_frac, + 'out_of_window_frac': out_of_window_frac, + } + + +results = {} +for var in VARS: + print(f"\n{'=' * 76}\n{var}\n{'=' * 76}") + anom = load_anom(var) + lin_rows = [] + z_rows = [] + for (w_lo, w_hi) in WINDOWS: + lin = linear_arm(anom, w_lo, w_hi) + zr = fisher_z_arm(anom, w_lo, w_hi) + lin_rows.append(lin) + z_rows.append(zr) + print(f"\nwindow ({w_lo},{w_hi}):") + print(f" LINEAR ci={lin['ci_uniform']:.6f} hi-lo={lin['hi_minus_lo']:.6f} " + f"sat%={lin['sat_frac'] * 100:.4f} true_out%={lin['out_of_window_frac'] * 100:.4f}") + print(f" FISHER_Z ci_med={zr['ci_interior_med']:.6f} " + f"(min={zr['ci_interior_min']:.6f} max={zr['ci_interior_max']:.6f}) " + f"sat%={zr['sat_frac'] * 100:.4f} true_out%={zr['out_of_window_frac'] * 100:.4f} " + f"scale={zr['scale']:.6f} zlo={zr['zlo']:.6f} zhi={zr['zhi']:.6f}") + + # ---- MANDATORY CHECK 1: linear CI strictly increasing + ratio proof ---- + ci_vals = [r['ci_uniform'] for r in lin_rows] + strictly_increasing = all(ci_vals[i] < ci_vals[i + 1] for i in range(len(ci_vals) - 1)) + hilo_ratio = lin_rows[-1]['hi_minus_lo'] / lin_rows[0]['hi_minus_lo'] + ci_ratio = ci_vals[-1] / ci_vals[0] + ratio_matches = abs(ci_ratio - hilo_ratio) < 1e-12 + print(f"\n [CHECK 1] linear ci strictly increasing across windows: {strictly_increasing}") + print(f" [CHECK 1] ci_vals = {ci_vals}") + print(f" [CHECK 1] ci_ratio(widest/narrowest) = {ci_ratio!r}") + print(f" [CHECK 1] hilo_ratio(widest/narrowest) = {hilo_ratio!r}") + print(f" [CHECK 1] |ci_ratio - hilo_ratio| = {abs(ci_ratio - hilo_ratio):.3e} " + f"matches_to_1e-12={ratio_matches}") + + # ---- MANDATORY CHECK 2: fisher-z curve reported, not gated ---- + z_ci_vals = [r['ci_interior_med'] for r in z_rows] + z_increasing = all(z_ci_vals[i] < z_ci_vals[i + 1] for i in range(len(z_ci_vals) - 1)) + z_decreasing = all(z_ci_vals[i] > z_ci_vals[i + 1] for i in range(len(z_ci_vals) - 1)) + z_monotone = z_increasing or z_decreasing + print(f" [CHECK 2] fisher-z interior-median ci across windows: {z_ci_vals}") + print(f" [CHECK 2] fisher-z monotone: {z_monotone} " + f"(increasing={z_increasing}, decreasing={z_decreasing})") + if not z_monotone: + print(" [CHECK 2] NON-MONOTONE -- reported as a FINDING to investigate, " + "not treated as an apparatus failure. See per-window scale/zlo/zhi above.") + + results[var] = { + 'windows': [{'w_lo': w_lo, 'w_hi': w_hi} for (w_lo, w_hi) in WINDOWS], + 'linear': lin_rows, + 'fisher_z': z_rows, + 'check1_linear_strictly_increasing': bool(strictly_increasing), + 'check1_ci_ratio_widest_over_narrowest': float(ci_ratio), + 'check1_hilo_ratio_widest_over_narrowest': float(hilo_ratio), + 'check1_abs_diff': float(abs(ci_ratio - hilo_ratio)), + 'check1_ratio_matches_to_1e-12': bool(ratio_matches), + 'check2_fisher_z_ci_interior_med_values': z_ci_vals, + 'check2_fisher_z_monotone': bool(z_monotone), + 'check2_fisher_z_increasing': bool(z_increasing), + 'check2_fisher_z_decreasing': bool(z_decreasing), + } + +overall_check1_pass = all( + results[v]['check1_linear_strictly_increasing'] and results[v]['check1_ratio_matches_to_1e-12'] + for v in VARS +) +print(f"\n{'=' * 76}") +print(f"OVERALL check1 (linear knob-reached-code proof) pass on all vars: {overall_check1_pass}") +for v in VARS: + print(f" {v}: fisher_z monotone = {results[v]['check2_fisher_z_monotone']}") + +results['_summary'] = { + 'vars': VARS, + 'windows': [{'w_lo': w_lo, 'w_hi': w_hi} for (w_lo, w_hi) in WINDOWS], + 'overall_check1_pass': bool(overall_check1_pass), + 'per_var_fisher_z_monotone': {v: results[v]['check2_fisher_z_monotone'] for v in VARS}, +} + +with open('ev4_window_sweep.json', 'w') as f: + json.dump(results, f, indent=2) +print("\nwrote ev4_window_sweep.json") diff --git a/probes/weather-p1/ev8_estimator_stability.json b/probes/weather-p1/ev8_estimator_stability.json new file mode 100644 index 000000000..d5b892723 --- /dev/null +++ b/probes/weather-p1/ev8_estimator_stability.json @@ -0,0 +1,588 @@ +{ + "n_pairs_per_resample": 200000, + "grid_n_gridpoints_NOT_USED_recorded_for_contrast": 1038240, + "k_seed_resamples": 25, + "seeds": [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25 + ], + "vars": [ + "2m_temperature", + "2m_dewpoint_temperature", + "10m_u_component_of_wind" + ], + "units": { + "2m_temperature": "K", + "2m_dewpoint_temperature": "K", + "10m_u_component_of_wind": "m/s" + }, + "shared_window_z": [ + -3.285941867486222, + 3.611137273885975 + ], + "jirak": { + "dropped_for_this_probe": true, + "reason": "index pairs are drawn independently and uniformly WITH replacement per resample (rng.integers, default replace=True); each of the N=200000 samples is therefore an IID draw from a fixed bivariate population regardless of the spatial correlation structure of the underlying field. Classical Fisher-z delta-method SE is used as a reference scale for n_eff instead of a Jirak weak-dependence rate." + }, + "preregistration": { + "expected_finding": "n_eff ~= n_pairs_per_resample (200000) for these IID index-pair draws. This is the EXPECTED outcome and is a FINDING, not a wiring alarm. n_eff << n_pairs_per_resample would indict the ESTIMATOR (a hidden dependence in the index-pair draw), not the underlying field.", + "correct_n_for_this_estimator": 200000, + "wrong_n_used_by_v1_spec": 1038240, + "wrong_n_reason": "v1 cited the grid size 721*1440=1038240 gridpoints as n; the P2 estimator (p2_probe.py:42,48) draws N=200000 INDEPENDENT UNIFORM INDEX PAIRS, not the grid -- that is the correct n for any CI statement about this estimator." + }, + "statistic_class": { + "is_codec_consistency_statistic": true, + "note": "rho(d_shared, truth) and rho(d_per_variable, truth) measure self-consistency of a DETERMINISTIC MONOTONE QUANTIZER against the SAME z-values it quantizes -- not a relationship between two independently-measured quantities. The deliverable of this probe is CI WIDTH and CROSS-SEED STABILITY of that statistic, not a p-value dressed up as inference." + }, + "no_pvalue_computed": true, + "results": { + "2m_temperature x 2m_dewpoint_temperature": { + "shared": { + "seed_count": 25, + "values_by_seed": [ + 0.9995729423739355, + 0.9996160797652496, + 0.9995898234370383, + 0.9995185184118438, + 0.9995599823322103, + 0.9995664005639123, + 0.9995562984081224, + 0.9995134898904035, + 0.9995418116918264, + 0.999583037153081, + 0.9995258863021326, + 0.9995540002247972, + 0.9995993811783749, + 0.9995526737445098, + 0.9995390300889124, + 0.9994841431629672, + 0.9995445229130183, + 0.9995394733198859, + 0.9995493949526264, + 0.999525903609043, + 0.9995426417364786, + 0.9995453006644321, + 0.9995791617898333, + 0.9995117836049295, + 0.9995288300414553 + ], + "mean": 0.9995496204544407, + "sd": 2.9877361602183034e-05, + "min": 0.9994841431629672, + "max": 0.9996160797652496, + "range": 0.0001319366022823898, + "ci_95_width_normal_approx": 0.00011711710631052214, + "ci_95_lo": 0.9994910619012854, + "ci_95_hi": 0.999608179007596, + "classical_se_theory_at_n_pairs": 2.013720094141001e-06, + "n_eff_estimate": 911.5256687338014, + "n_eff_over_n_pairs": 0.004557628343669007 + }, + "per_variable": { + "seed_count": 25, + "values_by_seed": [ + 0.999447632425782, + 0.9994629668200713, + 0.9994392987649783, + 0.9993906074871051, + 0.9994417167883134, + 0.9994329727358583, + 0.9994263157997967, + 0.9993851992221366, + 0.9994200441240733, + 0.9994512372179365, + 0.9994029203008958, + 0.9994198231245351, + 0.9994631521534922, + 0.9994306006871864, + 0.9994071721945487, + 0.9993668531630137, + 0.9994247331724968, + 0.9994253916457106, + 0.999417939258237, + 0.9993973684132235, + 0.999409253230225, + 0.9994225323425356, + 0.9994284467003558, + 0.9993915646258281, + 0.9994012051107202 + ], + "mean": 0.9994202779003624, + "sd": 2.405647717272568e-05, + "min": 0.9993668531630137, + "max": 0.9994631521534922, + "range": 9.62989904785827e-05, + "ci_95_width_normal_approx": 9.429965845072824e-05, + "ci_95_lo": 0.999373128071137, + "ci_95_hi": 0.9994674277295877, + "classical_se_theory_at_n_pairs": 2.591863992138682e-06, + "n_eff_estimate": 2324.5809870368735, + "n_eff_over_n_pairs": 0.011622904935184367 + } + }, + "2m_temperature x 10m_u_component_of_wind": { + "shared": { + "seed_count": 25, + "values_by_seed": [ + 0.9997611319659437, + 0.9997120303692658, + 0.9997649932307332, + 0.9997477285810774, + 0.9997111578470012, + 0.9997373335213048, + 0.9997359339313134, + 0.9997289269636326, + 0.9997185800319441, + 0.9997456071543322, + 0.9997322557442551, + 0.9997439790544449, + 0.9997438720286735, + 0.9997281443116485, + 0.9997105391659257, + 0.9997198589440766, + 0.9997276936398893, + 0.9997419448275038, + 0.9997259869913002, + 0.999748648447371, + 0.9997509074228678, + 0.9997476030719649, + 0.9997335273559806, + 0.9997507096681071, + 0.9997322936898208 + ], + "mean": 0.9997360555184152, + "sd": 1.4883947693765254e-05, + "min": 0.9997105391659257, + "max": 0.9997649932307332, + "range": 5.445406480752979e-05, + "ci_95_width_normal_approx": 5.834400331532585e-05, + "ci_95_lo": 0.9997068835167575, + "ci_95_hi": 0.9997652275200728, + "classical_se_theory_at_n_pairs": 1.1802486786645316e-06, + "n_eff_estimate": 1260.5760468653418, + "n_eff_over_n_pairs": 0.006302880234326709 + }, + "per_variable": { + "seed_count": 25, + "values_by_seed": [ + 0.8742823160635257, + 0.8725580657940212, + 0.8742700455978657, + 0.8735996919395955, + 0.8734372140171526, + 0.8738751421293475, + 0.874801168764891, + 0.8731443493589824, + 0.8731400690681725, + 0.8738913040037994, + 0.8739284651341969, + 0.8732065534558907, + 0.8733230278678422, + 0.8730539733332696, + 0.8741097391306385, + 0.8741890738412174, + 0.8732766076000185, + 0.8742074775419999, + 0.8735544666573449, + 0.8742943254490758, + 0.8738985861609684, + 0.8735254789202249, + 0.873824996832977, + 0.8754307980795385, + 0.8739315739317696 + ], + "mean": 0.8737901804269731, + "sd": 0.0006108124963128413, + "min": 0.8725580657940212, + "max": 0.8754307980795385, + "range": 0.002872732285517321, + "ci_95_width_normal_approx": 0.0023943410070466032, + "ci_95_lo": 0.8725930099234498, + "ci_95_hi": 0.8749873509304964, + "classical_se_theory_at_n_pairs": 0.0005288132934003761, + "n_eff_estimate": 149906.71025918314, + "n_eff_over_n_pairs": 0.7495335512959157 + } + }, + "2m_dewpoint_temperature x 10m_u_component_of_wind": { + "shared": { + "seed_count": 25, + "values_by_seed": [ + 0.9997294188605014, + 0.999718399207049, + 0.9997245347361846, + 0.9997392166450608, + 0.9997262055339732, + 0.9996964954924907, + 0.9997218339974447, + 0.9996869392549376, + 0.9997380431009277, + 0.9997206883650892, + 0.9996997928325538, + 0.9997392727611731, + 0.9997161978044206, + 0.9997105730891135, + 0.9997423125962694, + 0.9997121612171507, + 0.9997199169805333, + 0.999715685830799, + 0.999712097071861, + 0.999720577452122, + 0.9997346437378679, + 0.9997320375022964, + 0.9997164700434366, + 0.99969905671275, + 0.9996999084551953 + ], + "mean": 0.9997188991712481, + "sd": 1.4722561220648876e-05, + "min": 0.9996869392549376, + "max": 0.9997423125962694, + "range": 5.537334133187688e-05, + "ci_95_width_normal_approx": 5.771137996053571e-05, + "ci_95_lo": 0.9996900434812679, + "ci_95_hi": 0.9997477548612284, + "classical_se_theory_at_n_pairs": 1.2569538615359053e-06, + "n_eff_estimate": 1460.7915092552923, + "n_eff_over_n_pairs": 0.007303957546276461 + }, + "per_variable": { + "seed_count": 25, + "values_by_seed": [ + 0.8569472185981152, + 0.8571158060484332, + 0.8554697564352293, + 0.8564302854418944, + 0.8560805047156866, + 0.8559815352672817, + 0.8569840346681692, + 0.855368915009638, + 0.8554303559744142, + 0.8559549415355489, + 0.8565240982560552, + 0.8562644298649486, + 0.8563180124526995, + 0.8567849406883056, + 0.8564591914435901, + 0.8559186285615478, + 0.8561091713032641, + 0.8561268414372141, + 0.8553312752164172, + 0.8564720265995728, + 0.8564498669879053, + 0.8563259533159077, + 0.8561887531816814, + 0.8560671373800887, + 0.8572382879923308 + ], + "mean": 0.8562536787350377, + "sd": 0.0005227785478376521, + "min": 0.8553312752164172, + "max": 0.8572382879923308, + "range": 0.0019070127759136524, + "ci_95_width_normal_approx": 0.0020492542674681517, + "ci_95_lo": 0.8552290516013036, + "ci_95_hi": 0.8572783058687717, + "classical_se_theory_at_n_pairs": 0.0005966536831225125, + "n_eff_estimate": 260517.9295386108, + "n_eff_over_n_pairs": 1.302589647693054 + } + }, + "2m_temperature (within)": { + "shared": { + "seed_count": 25, + "values_by_seed": [ + 0.9995904732766816, + 0.9995624664517562, + 0.9995872130880127, + 0.9996010555488439, + 0.9995599006725214, + 0.9995962922621393, + 0.9995571417817292, + 0.999536404145236, + 0.9995564065988919, + 0.9995559013071772, + 0.9996146391412978, + 0.9995318706155661, + 0.9995828630794265, + 0.9995287574584798, + 0.9995469936866127, + 0.9995407940980712, + 0.9995611045725258, + 0.9995959812162757, + 0.9995273604578604, + 0.9995172913419434, + 0.9995516967710965, + 0.9995519916972958, + 0.9995533206032533, + 0.9995746255434059, + 0.9995323418389409 + ], + "mean": 0.9995605954902015, + "sd": 2.6155852508828202e-05, + "min": 0.9995172913419434, + "max": 0.9996146391412978, + "range": 9.734779935433036e-05, + "ci_95_width_normal_approx": 0.00010252905861322592, + "ci_95_lo": 0.9995093309608949, + "ci_95_hi": 0.9996118600195081, + "classical_se_theory_at_n_pairs": 1.9646597102783187e-06, + "n_eff_estimate": 1131.3930428926533, + "n_eff_over_n_pairs": 0.0056569652144632665 + }, + "per_variable": { + "seed_count": 25, + "values_by_seed": [ + 0.9996283705670711, + 0.9996188869566083, + 0.9996358964943012, + 0.999652721613433, + 0.9996117569508083, + 0.999640776560886, + 0.999613473800991, + 0.9995949977265098, + 0.9996199233109466, + 0.9996078748872677, + 0.9996554286529802, + 0.9995854199573383, + 0.9996321387947629, + 0.9995936445547785, + 0.9996020664232971, + 0.9996064013865522, + 0.9996190874810104, + 0.9996476324645321, + 0.999593717671712, + 0.9995777439685474, + 0.9996158243157592, + 0.9996100987458988, + 0.9996113856910736, + 0.999630528736039, + 0.9995953670024192 + ], + "mean": 0.999616046588621, + "sd": 2.0681929482516518e-05, + "min": 0.9995777439685474, + "max": 0.9996554286529802, + "range": 7.768468443281051e-05, + "ci_95_width_normal_approx": 8.107167447254201e-05, + "ci_95_lo": 0.9995755107513847, + "ci_95_hi": 0.9996565824258573, + "classical_se_theory_at_n_pairs": 1.716775090296173e-06, + "n_eff_estimate": 1381.060173559357, + "n_eff_over_n_pairs": 0.0069053008677967855 + } + }, + "2m_dewpoint_temperature (within)": { + "shared": { + "seed_count": 25, + "values_by_seed": [ + 0.9995401064396318, + 0.9995459193095025, + 0.9995310732835152, + 0.9995021591761837, + 0.9995777191638366, + 0.9995479638532814, + 0.9995312723074434, + 0.9995210197021076, + 0.999471630808959, + 0.9995835898371563, + 0.9995974337233902, + 0.9995345619299123, + 0.9995108006004857, + 0.9994961227779684, + 0.9995326517930051, + 0.9995423125967201, + 0.9995412556336873, + 0.9995442259747992, + 0.999514791982735, + 0.9995636185914979, + 0.9995633647742278, + 0.999529548204843, + 0.9995204429049701, + 0.9995518340620182, + 0.9995432968059945 + ], + "mean": 0.9995375486495149, + "sd": 2.7733904985825377e-05, + "min": 0.999471630808959, + "max": 0.9995974337233902, + "range": 0.00012580291443120561, + "ci_95_width_normal_approx": 0.0001087149107032765, + "ci_95_lo": 0.9994831911941633, + "ci_95_hi": 0.9995919061048666, + "classical_se_theory_at_n_pairs": 2.0676826113237655e-06, + "n_eff_estimate": 1114.6528491143047, + "n_eff_over_n_pairs": 0.005573264245571524 + }, + "per_variable": { + "seed_count": 25, + "values_by_seed": [ + 0.9996410439375155, + 0.9996247230350557, + 0.999622358634043, + 0.9995981074395485, + 0.9996559805851423, + 0.9996320797708269, + 0.9996178181789303, + 0.9995981480468349, + 0.9995887644249526, + 0.9996521432050296, + 0.999667483047408, + 0.999621531670056, + 0.9996077001977185, + 0.999591532931126, + 0.9996202905639765, + 0.9996261046011138, + 0.9996238481247987, + 0.9996326191144407, + 0.999590397810163, + 0.9996379862309024, + 0.9996336595862365, + 0.9996190080594924, + 0.9996081294178191, + 0.9996350821044071, + 0.9996216872313268 + ], + "mean": 0.9996227291179545, + "sd": 2.029375245691261e-05, + "min": 0.9995887644249526, + "max": 0.999667483047408, + "range": 7.871862245545369e-05, + "ci_95_width_normal_approx": 7.955004848092053e-05, + "ci_95_lo": 0.999582954093714, + "ci_95_hi": 0.9996625041421949, + "classical_se_theory_at_n_pairs": 1.686901061200086e-06, + "n_eff_estimate": 1384.9043597205355, + "n_eff_over_n_pairs": 0.006924521798602677 + } + }, + "10m_u_component_of_wind (within)": { + "shared": { + "seed_count": 25, + "values_by_seed": [ + 0.9998173514221279, + 0.9998222654495103, + 0.9998111177141497, + 0.9998173772311906, + 0.9998243467745702, + 0.999826845930373, + 0.9998171215242767, + 0.9998121030725273, + 0.9998235532848196, + 0.9998096014013144, + 0.9998155368612089, + 0.9998090899270288, + 0.9998209442453634, + 0.9998074893679497, + 0.9998193922097544, + 0.9998241338600086, + 0.9998253370867749, + 0.9998135000854033, + 0.9998203537780489, + 0.9998186942186216, + 0.9998163395897881, + 0.9998225996403607, + 0.9998236445537364, + 0.9998269136611495, + 0.9998201612336429 + ], + "mean": 0.9998186325649481, + "sd": 5.686174776525023e-06, + "min": 0.9998074893679497, + "max": 0.9998269136611495, + "range": 1.942429319978789e-05, + "ci_95_width_normal_approx": 2.2289395719394183e-05, + "ci_95_lo": 0.9998074878670884, + "ci_95_hi": 0.9998297772628078, + "classical_se_theory_at_n_pairs": 8.11032356579582e-07, + "n_eff_estimate": 4071.7352744474542, + "n_eff_over_n_pairs": 0.02035867637223727 + }, + "per_variable": { + "seed_count": 25, + "values_by_seed": [ + 0.9996596335269514, + 0.9996304420295159, + 0.9996149115255794, + 0.9996540913348506, + 0.9996397814729688, + 0.9996425338044466, + 0.9996220613759306, + 0.9996538173774505, + 0.9996340104303487, + 0.9996093777269214, + 0.9996174759308607, + 0.9996115166119767, + 0.9996509479778846, + 0.9995961603225605, + 0.9996131343427581, + 0.9996529831559928, + 0.9996302831797714, + 0.9996162416314941, + 0.9996362967586293, + 0.9996324974972011, + 0.9996391691307069, + 0.9996162825071988, + 0.999636222689724, + 0.9996520821343936, + 0.9996502585249007 + ], + "mean": 0.9996324885200407, + "sd": 1.7314925754709582e-05, + "min": 0.9995961603225605, + "max": 0.9996596335269514, + "range": 6.347320439092918e-05, + "ci_95_width_normal_approx": 6.787326228380723e-05, + "ci_95_lo": 0.9995985518888988, + "ci_95_hi": 0.9996664251511826, + "classical_se_theory_at_n_pairs": 1.6432716141608782e-06, + "n_eff_estimate": 1804.3616827622006, + "n_eff_over_n_pairs": 0.009021808413811004 + } + } + }, + "finding": { + "measured_pattern": "MIXED, not uniformly n_eff ~= n_pairs as pre-registered. 2/12 series land within 0.3x-3x of n_eff/n_pairs == 1 (the pre-registered expectation); 10/12 land BELOW 0.3x (n_eff far smaller than n_pairs); 0/12 land ABOVE 3x.", + "series_near_expected_0p3_to_3x": { + "2m_temperature x 10m_u_component_of_wind / per_variable": 0.7495335512959157, + "2m_dewpoint_temperature x 10m_u_component_of_wind / per_variable": 1.302589647693054 + }, + "series_far_below_0p3x": { + "2m_temperature x 2m_dewpoint_temperature / shared": 0.004557628343669007, + "2m_temperature x 2m_dewpoint_temperature / per_variable": 0.011622904935184367, + "2m_temperature x 10m_u_component_of_wind / shared": 0.006302880234326709, + "2m_dewpoint_temperature x 10m_u_component_of_wind / shared": 0.007303957546276461, + "2m_temperature (within) / shared": 0.0056569652144632665, + "2m_temperature (within) / per_variable": 0.0069053008677967855, + "2m_dewpoint_temperature (within) / shared": 0.005573264245571524, + "2m_dewpoint_temperature (within) / per_variable": 0.006924521798602677, + "10m_u_component_of_wind (within) / shared": 0.02035867637223727, + "10m_u_component_of_wind (within) / per_variable": 0.009021808413811004 + }, + "series_far_above_3x": {}, + "interpretation_caveat": "The near-expected cases are the two cross-variable per_variable-floor pairs whose mean rho sits well below unity (~0.86-0.87: T x u and Td x u, per-variable floor). The far-below cases are every series whose mean rho sits at or above ~0.9994 (all shared-floor series, both within-variable controls under both floors, and the shared-floor cross-variable pairs). This split tracks correlation LEVEL, not pair identity, which is consistent with a KNOWN limitation of the classical Fisher-z delta-method reference formula rather than necessarily a hidden dependence in the index-pair draw: d_shared/d_per_variable are INTEGER-valued (256-bucket quantizations), so Spearman rank correlation on them carries heavy TIES, and the untied delta-method variance formula (1-rho^2)^2/(n-3) is known to understate variance near rho~1 under tie-heavy discretization. This probe does NOT adjudicate that hypothesis (a ties-corrected Spearman variance estimator would be needed, out of scope here) -- it reports the measured n_eff/n_pairs split honestly and flags the caveat rather than asserting either \"the estimator is indicted\" or \"n_eff ~= n confirmed\" uniformly." + } +} \ No newline at end of file diff --git a/probes/weather-p1/ev8_estimator_stability.py b/probes/weather-p1/ev8_estimator_stability.py new file mode 100644 index 000000000..5edd78966 --- /dev/null +++ b/probes/weather-p1/ev8_estimator_stability.py @@ -0,0 +1,274 @@ +"""EV-8 (v2) -- stability of the P2 estimator via seed-resampling. + +Per `.claude/plans/weather-substrate-evaluation-v1.md` Section 3, EV-8 (v2, +post-audit): the v1 spec was a rubber-stamp ("state effective-n and re-grade" +can never fail) and used the WRONG n (the 1,038,240 gridpoints instead of the +200,000 INDEPENDENT UNIFORM INDEX PAIRS p2_probe.py actually draws -- see +p2_probe.py:42,48). + +This probe re-runs p2_probe.py's estimator K>=20 times, varying ONLY the RNG +seed that drives the index-pair draw (nothing else -- the z-values, the +shared/per-variable quantization floors, and the bucket boundaries are all +deterministic and computed once, exactly as p2_probe.py computes them), and +reports per-correlation mean / sd / min / max / CI width across seeds. + +PRE-REGISTERED (before interpreting any number below): for two independent +uniform index draws (ia, ib each drawn fresh per seed via rng.integers with +replace=True -- see p2_probe.py:48), n_eff ~= n=200,000 is the EXPECTED +finding. That outcome is a FINDING (the estimator behaves as a genuinely IID +sampling design), NOT a wiring alarm. n_eff << n would indict the ESTIMATOR +(a hidden dependence smuggled into the pairing/index-draw), not the +underlying field. + +rho(d_shared, truth) and rho(d_per_variable, truth) are CODEC-CONSISTENCY +statistics: d_shared / d_per_variable are DETERMINISTIC MONOTONE +QUANTIZATIONS of the SAME z-values `truth` is computed from (p2_probe.py +:29-33,49-51) -- the correlation measures self-consistency of a lossy +quantizer against its own lossless input, not a relationship between two +independently-measured quantities. Accordingly the deliverable of this probe +is CI WIDTH and CROSS-SEED STABILITY of that self-consistency statistic -- +NOT a p-value dressed up as inference. No p-value is computed anywhere in +this file. + +Jirak 2016 weak-dependence rates (I-NOISE-FLOOR-JIRAK) are the mandatory +citation for POOLED/SPATIAL statistics on this workspace's fields, because +the underlying grid has weakly-dependent (spatially correlated) values. That +does NOT apply to the quantity resampled here: p2_probe.py's index pairs +(ia, ib) are drawn independently and uniformly WITH replacement per resample +(rng.integers, default replace=True), so each of the N=200,000 samples per +seed is an IID draw (with replacement) from a fixed bivariate population -- +the spatial correlation of the underlying field is irrelevant to the +SAMPLING design of THIS estimator (it would matter for a statistic computed +directly over spatially-adjacent points, which this is not). Per the plan's +own escape clause for EV-8 ("Jirak's `p` pinned in the artifact or the Jirak +citation dropped for the resampling CI"): the Jirak citation is DROPPED here +and the classical IID Fisher-z delta-method SE is used instead, strictly as +a REFERENCE SCALE to express n_eff -- never as a significance/p-value claim. +""" +import numpy as np +import json +from scipy import stats as st + +K_SEEDS = list(range(1, 26)) # K=25 independent seed-resamples (>= 20 required by spec) +N = 200_000 # the P2 estimator's own n -- p2_probe.py:42 +GRID_N = 721 * 1440 # 1,038,240 -- the WRONG n the v1 spec used; recorded for contrast, never used as n below + +VARS = ['2m_temperature', '2m_dewpoint_temperature', '10m_u_component_of_wind'] +UNITS = {'2m_temperature': 'K', '2m_dewpoint_temperature': 'K', '10m_u_component_of_wind': 'm/s'} + + +def load(v): + """Identical to p2_probe.py:load -- zonal-mean anomaly of the raw field.""" + a = np.load(f'fixture/{v}.npy').astype(np.float64) + assert np.isfinite(a).all() + return a - a.mean(axis=1, keepdims=True) + + +def quant(x, lo, hi): + """Identical to p2_probe.py:quant -- 256 uniform buckets over [lo, hi]; + only the bucket INDEX is needed here (p2_probe.py's `cen` is unused by + the correlation computation, so it is dropped).""" + idx = np.clip(np.floor((x - lo) / (hi - lo) * 256), 0, 255).astype(np.uint8) + return idx + + +anom = {v: load(v) for v in VARS} +# Identical to p2_probe.py: standardize (identity after scaling; §12.1 +# measured arctanh as actively harmful here -- not re-litigated by this probe). +z = {v: (anom[v].ravel() - anom[v].mean()) / anom[v].std() for v in VARS} + +# Deterministic (no randomness, computed ONCE): the shared canonical floor +# and the per-variable floors, exactly as p2_probe.py computes them. +pool = np.concatenate([z[v] for v in VARS]) +s_lo, s_hi = np.percentile(pool, [0.4, 99.6]) +shared_idx = {v: quant(z[v], s_lo, s_hi) for v in VARS} +per_idx = {v: quant(z[v], *np.percentile(z[v], [0.4, 99.6])) for v in VARS} + +pairs = [(VARS[i], VARS[j]) for i in range(len(VARS)) for j in range(i + 1, len(VARS))] + +# key -> floor_type -> list of rho values, one entry per seed +series = {} +for va, vb in pairs: + series.setdefault(f"{va} x {vb}", {})['shared'] = [] + series[f"{va} x {vb}"]['per_variable'] = [] +for v in VARS: + series.setdefault(f"{v} (within)", {})['shared'] = [] + series[f"{v} (within)"]['per_variable'] = [] + +for seed in K_SEEDS: + rng = np.random.default_rng(seed) + # cross-variable pairs (in the SAME order p2_probe.py draws them, so a + # given seed draws the same index streams for the same key across runs) + for va, vb in pairs: + ia = rng.integers(0, z[va].size, N) + ib = rng.integers(0, z[vb].size, N) + truth = np.abs(z[va][ia] - z[vb][ib]) + d_sh = np.abs(shared_idx[va][ia].astype(np.int32) - shared_idx[vb][ib].astype(np.int32)) + d_pv = np.abs(per_idx[va][ia].astype(np.int32) - per_idx[vb][ib].astype(np.int32)) + r_sh = st.spearmanr(d_sh, truth).statistic + r_pv = st.spearmanr(d_pv, truth).statistic + series[f"{va} x {vb}"]['shared'].append(float(r_sh)) + series[f"{va} x {vb}"]['per_variable'].append(float(r_pv)) + # within-variable controls + for v in VARS: + ia = rng.integers(0, z[v].size, N) + ib = rng.integers(0, z[v].size, N) + truth = np.abs(z[v][ia] - z[v][ib]) + r_sh = st.spearmanr( + np.abs(shared_idx[v][ia].astype(np.int32) - shared_idx[v][ib].astype(np.int32)), truth + ).statistic + r_pv = st.spearmanr( + np.abs(per_idx[v][ia].astype(np.int32) - per_idx[v][ib].astype(np.int32)), truth + ).statistic + series[f"{v} (within)"]['shared'].append(float(r_sh)) + series[f"{v} (within)"]['per_variable'].append(float(r_pv)) + + +def classical_se(rho, n): + """Fisher-z delta-method SE of an IID-sampled correlation coefficient: + Var(rho) ~= (1-rho^2)^2 / (n-3) => SE(rho) ~= (1-rho^2)/sqrt(n-3). + Used strictly as a REFERENCE SCALE for n_eff -- see module docstring + re: dropping the Jirak weak-dependence citation for THIS probe.""" + return (1.0 - rho ** 2) / np.sqrt(max(n - 3, 1)) + + +results = {} +for key, by_floor in series.items(): + results[key] = {} + for floor_type, vals in by_floor.items(): + arr = np.array(vals) + mean = float(arr.mean()) + sd = float(arr.std(ddof=1)) # sample sd across the K seed-resamples + se_theory = classical_se(mean, N) + # Solve se_theory(n_eff) = observed sd for n_eff at the SAME mean rho: + # sd = (1-rho^2)/sqrt(n_eff-3) => n_eff = ((1-rho^2)/sd)^2 + 3 + n_eff = float(((1.0 - mean ** 2) / sd) ** 2 + 3.0) if sd > 0 else None + results[key][floor_type] = { + 'seed_count': len(vals), + 'values_by_seed': vals, + 'mean': mean, + 'sd': sd, + 'min': float(arr.min()), + 'max': float(arr.max()), + 'range': float(arr.max() - arr.min()), + 'ci_95_width_normal_approx': float(2 * 1.959964 * sd), + 'ci_95_lo': float(mean - 1.959964 * sd), + 'ci_95_hi': float(mean + 1.959964 * sd), + 'classical_se_theory_at_n_pairs': float(se_theory), + 'n_eff_estimate': n_eff, + 'n_eff_over_n_pairs': (n_eff / N) if n_eff is not None else None, + } + +out = { + 'n_pairs_per_resample': N, + 'grid_n_gridpoints_NOT_USED_recorded_for_contrast': GRID_N, + 'k_seed_resamples': len(K_SEEDS), + 'seeds': K_SEEDS, + 'vars': VARS, + 'units': UNITS, + 'shared_window_z': [float(s_lo), float(s_hi)], + 'jirak': { + 'dropped_for_this_probe': True, + 'reason': ( + 'index pairs are drawn independently and uniformly WITH replacement ' + 'per resample (rng.integers, default replace=True); each of the ' + 'N=200000 samples is therefore an IID draw from a fixed bivariate ' + 'population regardless of the spatial correlation structure of the ' + 'underlying field. Classical Fisher-z delta-method SE is used as a ' + 'reference scale for n_eff instead of a Jirak weak-dependence rate.' + ), + }, + 'preregistration': { + 'expected_finding': ( + 'n_eff ~= n_pairs_per_resample (200000) for these IID index-pair ' + 'draws. This is the EXPECTED outcome and is a FINDING, not a wiring ' + 'alarm. n_eff << n_pairs_per_resample would indict the ESTIMATOR ' + '(a hidden dependence in the index-pair draw), not the underlying ' + 'field.' + ), + 'correct_n_for_this_estimator': N, + 'wrong_n_used_by_v1_spec': GRID_N, + 'wrong_n_reason': ( + 'v1 cited the grid size 721*1440=1038240 gridpoints as n; the P2 ' + 'estimator (p2_probe.py:42,48) draws N=200000 INDEPENDENT UNIFORM ' + 'INDEX PAIRS, not the grid -- that is the correct n for any CI ' + 'statement about this estimator.' + ), + }, + 'statistic_class': { + 'is_codec_consistency_statistic': True, + 'note': ( + 'rho(d_shared, truth) and rho(d_per_variable, truth) measure ' + 'self-consistency of a DETERMINISTIC MONOTONE QUANTIZER against the ' + 'SAME z-values it quantizes -- not a relationship between two ' + 'independently-measured quantities. The deliverable of this probe is ' + 'CI WIDTH and CROSS-SEED STABILITY of that statistic, not a p-value ' + 'dressed up as inference.' + ), + }, + 'no_pvalue_computed': True, + 'results': results, +} + +# Honest, PROGRAMMATICALLY-derived summary of what the pre-registered +# "n_eff ~= n_pairs" expectation actually found. This is a description of +# the measured pattern (bucketed by ratio band), NOT a p-value and NOT a +# verdict -- see 'interpretation_caveat' for why a single verdict would +# overreach what this probe alone can establish. +n_eff_ratios = { + f"{key} / {floor_type}": r['n_eff_over_n_pairs'] + for key, by_floor in results.items() + for floor_type, r in by_floor.items() + if r['n_eff_over_n_pairs'] is not None +} +near_expected = {k: v for k, v in n_eff_ratios.items() if 0.3 <= v <= 3.0} +far_below = {k: v for k, v in n_eff_ratios.items() if v < 0.3} +far_above = {k: v for k, v in n_eff_ratios.items() if v > 3.0} + +out['finding'] = { + 'measured_pattern': ( + 'MIXED, not uniformly n_eff ~= n_pairs as pre-registered. ' + f'{len(near_expected)}/{len(n_eff_ratios)} series land within 0.3x-3x ' + f'of n_eff/n_pairs == 1 (the pre-registered expectation); ' + f'{len(far_below)}/{len(n_eff_ratios)} land BELOW 0.3x (n_eff far ' + f'smaller than n_pairs); {len(far_above)}/{len(n_eff_ratios)} land ' + f'ABOVE 3x.' + ), + 'series_near_expected_0p3_to_3x': near_expected, + 'series_far_below_0p3x': far_below, + 'series_far_above_3x': far_above, + 'interpretation_caveat': ( + 'The near-expected cases are the two cross-variable per_variable-floor ' + 'pairs whose mean rho sits well below unity (~0.86-0.87: T x u and ' + 'Td x u, per-variable floor). The far-below cases are every series ' + 'whose mean rho sits at or above ~0.9994 (all shared-floor series, ' + 'both within-variable controls under both floors, and the ' + 'shared-floor cross-variable pairs). This split tracks correlation ' + 'LEVEL, not pair identity, which is consistent with a KNOWN ' + 'limitation of the classical Fisher-z delta-method reference formula ' + 'rather than necessarily a hidden dependence in the index-pair draw: ' + 'd_shared/d_per_variable are INTEGER-valued (256-bucket ' + 'quantizations), so Spearman rank correlation on them carries heavy ' + 'TIES, and the untied delta-method variance formula ' + '(1-rho^2)^2/(n-3) is known to understate variance near rho~1 under ' + 'tie-heavy discretization. This probe does NOT adjudicate that ' + 'hypothesis (a ties-corrected Spearman variance estimator would be ' + 'needed, out of scope here) -- it reports the measured n_eff/n_pairs ' + 'split honestly and flags the caveat rather than asserting either ' + '"the estimator is indicted" or "n_eff ~= n confirmed" uniformly.' + ), +} + +with open('ev8_estimator_stability.json', 'w') as f: + json.dump(out, f, indent=2) + +print(f"K={len(K_SEEDS)} seed-resamples, N={N} pairs/resample " + f"(NOT the grid size {GRID_N})\n") +print(f"{'key':<40s} {'floor':<13s} {'mean':>10s} {'sd':>10s} " + f"{'min':>10s} {'max':>10s} {'ci95_width':>12s} {'n_eff/n':>10s}") +for key, by_floor in results.items(): + for floor_type, r in by_floor.items(): + n_eff_ratio = r['n_eff_over_n_pairs'] + ratio_str = f"{n_eff_ratio:10.4f}" if n_eff_ratio is not None else f"{'None':>10s}" + print(f"{key:<40s} {floor_type:<13s} {r['mean']:10.6f} {r['sd']:10.2e} " + f"{r['min']:10.6f} {r['max']:10.6f} {r['ci_95_width_normal_approx']:12.2e} {ratio_str}") diff --git a/probes/weather-p1/fixture/manifest.json b/probes/weather-p1/fixture/manifest.json new file mode 100644 index 000000000..c4544a04e --- /dev/null +++ b/probes/weather-p1/fixture/manifest.json @@ -0,0 +1,62 @@ +{ + "t547476/2m_temperature": { + "time_index": 547476, + "shape": [ + 721, + 1440 + ], + "nonfinite": 0, + "min": 195.8684539794922, + "max": 323.6919250488281 + }, + "t547476/2m_dewpoint_temperature": { + "time_index": 547476, + "shape": [ + 721, + 1440 + ], + "nonfinite": 0, + "min": 192.04747009277344, + "max": 302.839111328125 + }, + "t547476/10m_u_component_of_wind": { + "time_index": 547476, + "shape": [ + 721, + 1440 + ], + "nonfinite": 0, + "min": -24.74393081665039, + "max": 20.80101203918457 + }, + "t543852/2m_temperature": { + "time_index": 543852, + "shape": [ + 721, + 1440 + ], + "nonfinite": 0, + "min": 219.7332763671875, + "max": 315.8507080078125 + }, + "t543852/2m_dewpoint_temperature": { + "time_index": 543852, + "shape": [ + 721, + 1440 + ], + "nonfinite": 0, + "min": 216.30738830566406, + "max": 301.58349609375 + }, + "t543852/10m_u_component_of_wind": { + "time_index": 543852, + "shape": [ + 721, + 1440 + ], + "nonfinite": 0, + "min": -22.079614639282227, + "max": 22.135465621948242 + } +} \ No newline at end of file