diff --git a/.claude/board/EPIPHANIES.md b/.claude/board/EPIPHANIES.md index f1340465..8e466264 100644 --- a/.claude/board/EPIPHANIES.md +++ b/.claude/board/EPIPHANIES.md @@ -1,3 +1,92 @@ +## 2026-08-12 — E-THE-DISPLACEMENT-FILTER-ATE-THE-STRANDED-STRATUM-1 + +**Status:** FINDING `[G]` — W6 RUN (`comet_tail_w6.py`/`.json`), audited +(raw predictors committed per storm after a self-caught gap), both a +sample-composition finding and a genuine model-level negative result. + +**The mechanism-level result, stated first because it is the load-bearing +one.** Report §10.2's dipole vector-sum model (`D = c_geo·P_geo + +c_bow·P_bow`, a background-high neighbor predictor plus a relative-motion +bow-wave predictor) was fit on CT-F14's 19 stored storms and **VOIDED by its +own pre-registered anti-vacuity control**: single-geo R²=−0.104 (worse than +predicting the mean); the permuted-P_bow control R²=−0.071 and the +rotated-90° control R²=−0.062 BOTH exceed `single-geo + 0.03 = −0.074` — two +deliberately wrong references score as well as or better than the real +predictor. This is not a marginal miss — it is the anti-vacuity control +doing exactly its job: rejecting a fit that has nothing to identify. + +**Checked before concluding it was a clean negative, per the standing +measurement-skeptic discipline: no implementation bug found in the FIT +itself, but a sign convention AND a units error were found in the +NARRATIVE around it (codex + CodeRabbit P2/Major on PR #940, both real, +fixed before merge).** One storm was independently re-fetched and hand- +audited; extended to the full sample by committing the raw +`D`/`P_geo`/`P_bow`/`A_H`/`d_H`/`v_rel_ms` per storm. Every value is +physically sane (`A_H > 0` always, `d_H` inside the 600–2500 km annulus +always, `v_rel` 2.8–27.9 m/s). + +**Sign:** `spine()`'s raw fit coefficient points toward the storm's HIGH +side (the gradient of increasing residual pressure), while `P_geo`/`P_bow` +both point toward the LOW side by construction — the exact convention +`low_pole_bearing()` makes explicit via its own `(ph + π) % (2π)` flip. `D` +is correctly `−spine(...)`; the first draft used the unflipped `coef`. +Corrected: **`c_geo = +0.407` (the physically predicted positive sign — +CORRECT)**, `c_bow = −0.0006` km⁻¹ (predicted positive — **wrong sign, but +small**). Verified algebraically and numerically that this flip changes +NOTHING about R² or the B0/B1 VOID verdicts (OLS is odd-symmetric in the +fit target) — only the coefficient signs and the sentence describing them. + +**Units:** `D`/`P_geo` are [Pa/km]; `P_bow` is [Pa] — `c_geo` is +dimensionless, `c_bow` carries km⁻¹, and OLS coefficients rescale inversely +under column rescaling while R²/fitted-values stay fixed. **Raw `|c_bow|` +was never valid evidence of "no measurable weight"**, and comparing +`|P_bow|` to `|D|` directly (147×, as first reported) compounded the same +mistake — Pa is not comparable to Pa/km at all. The dimensionally valid +measure is the fitted CONTRIBUTION `|c_bow·P_bow|` against `|D|`, both in +Pa/km: mean `|D|`=0.745, mean `|c_geo·P_geo|`=0.186 (25 % of `|D|`), mean +`|c_bow·P_bow|`=0.068 (9 % of `|D|`) — the geo contribution is ~2.7× the +bow contribution, MODEST rather than "no weight," and both remain +consistent with the R²<0 finding that neither predictor meaningfully +explains `D`'s variance. + +**The finding that generalizes past this one probe: an anti-vacuity control +can be voided by SAMPLE COMPOSITION, and the reason is arithmetic, not +physics.** B3's stranded-vs-moving stratification (`|v_storm| < 8 m/s`, +the report's own named test of the "stranded-rescue" reading) came back +**n=0 for the stranded stratum** — every one of the 19 storms has +`|v_storm| ≥ 12.54 m/s`. This is not a null result about storm motion; it is +a DIRECT ARITHMETIC CONSEQUENCE of CT-F14's own qualifying filter +(`displacement_km ≥ 250` over the 6 h window): `250 km / 6 h = 11.574 m/s`, +a hard floor on `|v_storm|` for ANY storm admitted to the sample. **A filter +built to select clearly-moving storms for a displacement-scoring test +silently and permanently excludes the storms a LATER, differently-motivated +test (stranded-rescue) needs to see.** The stranded-rescue claim is +therefore **UNTESTABLE on this sample, not refuted** — the untestability was +knowable from the filter's own arithmetic before a single storm was fetched, +and wasn't checked until B3 came back empty. + +**Consequence, stated as a reusable rule:** before scoring ANY new +hypothesis against an EXISTING filtered sample, check whether the sample's +own selection criterion is compatible with the new hypothesis's own +discriminating variable — arithmetically, not by running the probe and +discovering an empty stratum after the fact. A filter selected for one +purpose (fast, clearly-displaced storms, easy to center-find and score +against displacement) is not neutral with respect to every future question; +it is a specific cut through the underlying population, and every later +probe inherits that cut whether or not it is the cut that probe needs. + +**Consequence for the report and for CT-F17.** The report's §10.2 vector-sum +model, AS SPECIFIED, is disconfirmed on this sample — not "unproven," not +"needs more data" in the ordinary sense, but VOID by its own control. CT-F17 +(the fresh-sample verdict, gated on W6's result + an independent adversarial +audit) is now moot **for this form of the model** — a fresh-sample test of a +model that already fails identifiability on the stored sample is not the +next useful step. Any REVISED form of the vector-sum model (multiple +neighbors, a nonlinear bow term, per-storm coefficients) would need its own +W6-shaped mechanistic test before earning a CT-F17 slot; a genuine +stranded-rescue test needs a sample built without (or explicitly retaining +slow storms despite) a displacement floor. + ## 2026-08-12 — E-ON-A-GOLDEN-LATTICE-LOCALITY-IS-FIBONACCI-MEMBERSHIP-1 **Status:** FINDING `[G]` — W5 RUN (`spiral_adi_probe.py`/`.json`) + a diff --git a/.claude/board/PR_ARC_INVENTORY.md b/.claude/board/PR_ARC_INVENTORY.md index a72d7e3a..c9198b8d 100644 --- a/.claude/board/PR_ARC_INVENTORY.md +++ b/.claude/board/PR_ARC_INVENTORY.md @@ -1,3 +1,71 @@ +## 2026-08-12 — the W6 dipole-deconvolution RUN lands — VOID by its own anti-vacuity control, plus a reusable sample-composition lesson (PR pending) + +- **Added.** `comet_tail_w6.py`/`.json` — the report §10.2 vector-sum + model (`D = c_geo·P_geo + c_bow·P_bow`) fit on CT-F14/F16's 19 stored + storms, per `weather-w-probes-v1.md` §3's pre-registered bars. New + epiphany `E-THE-DISPLACEMENT-FILTER-ATE-THE-STRANDED-STRATUM-1`. +- **Locked — B0 VOIDS, and the model has nothing to identify.** + Single-geo R²=−0.104 (worse than the mean); both anti-vacuity controls + (permuted P_bow, P_bow rotated +90°) score at −0.071/−0.062, both + clearing the ≤`single-geo + 0.03` bar's ceiling of −0.074. B1/B2 + correctly report VOID (B0 failed), not their own numbers, per the + pre-registered rule. +- **Locked — checked for an implementation bug before calling it a clean + negative, per the standing measurement-skeptic discipline; found none + in the fit — but found a sign convention error AND a units error in + the FIRST DRAFT's narrative (codex + CodeRabbit P2/Major on PR #940, + both fixed before merge).** One storm independently re-fetched and + hand-audited (t0=54358); extended to all 19 by committing the raw + predictors per storm (the first run shipped only derived bearings — a + self-caught audit gap, fixed and re-run before any reviewer needed to + find it). **Sign:** `spine()`'s raw coefficient points toward the + storm's HIGH side (increasing residual), while `P_geo`/`P_bow` both + point toward the LOW side by construction — `D = −spine(...)` matches + the same `(ph+π)` convention `low_pole_bearing()` already uses. + Corrected: `c_geo = +0.407` (physically predicted positive — CORRECT, + the first draft's "−0.41, wrong sign" had the polarity backward); + `c_bow = −0.0006` km⁻¹ (predicted positive — wrong sign, small). + Verified algebraically and numerically that the flip changes NOTHING + about R²/B0/B1 (OLS is odd-symmetric in the fit target) — only the + coefficient signs and their description. **Units:** `D`/`P_geo` are + Pa/km, `P_bow` is Pa, so `c_geo` is dimensionless but `c_bow` carries + km⁻¹ — raw `|c_bow|` (or the first draft's `|P_bow|`/`|D|` ≈ 147× + comparison) was never valid evidence of "no measurable weight," since + OLS coefficients rescale inversely under column rescaling. The + dimensionally valid measure is the fitted contribution + `|c_bow·P_bow|` vs `|D|`, both Pa/km: mean `|D|`=0.745, geo + contribution 25 % of `|D|`, bow contribution 9 % of `|D|` — modest, + not "no weight," and consistent throughout with R²<0. +- **Locked — the reusable lesson: B3's stranded stratum is EMPTY (n=0), + and the reason is arithmetic, not physics.** `min(|v_storm|) = 12.54 + m/s` across all 19 storms — CT-F14's own qualifying filter + (`displacement_km ≥ 250` over 6 h) mathematically implies `|v_storm| ≥ + 250 km / 6 h = 11.574 m/s` for ANY admitted storm. **A filter selected + for one purpose (fast, cleanly-displaced storms) is not neutral for a + LATER, differently-motivated test (stranded-rescue) — the + incompatibility is knowable from the filter's own arithmetic before a + single fetch, and wasn't checked until the stratum came back empty.** + Filed as a reusable pre-flight check for any future probe reusing an + existing filtered sample. +- **Locked — consequence for CT-F17.** The vector-sum model AS + SPECIFIED is disconfirmed on this sample, void by its own control — + not "unproven," not "needs more data." CT-F17's gate (W6's result + + independent adversarial audit) is now moot **for this form of the + model**: a fresh-sample verdict on a model that fails identifiability + on the stored sample is not the useful next step. A REVISED model + form needs its own W6-shaped mechanistic test first; a genuine + stranded-rescue test needs a sample built without the displacement + floor. +- **Docs.** `weather-w-probes-v1.md` §3 gets the full RUN section. + `STATUS_BOARD.md` D-W6 moves from Queued to the RUN verdict. +- **Confidence.** High — every figure re-derived from the committed + JSON (which now carries the raw per-storm predictors, not only the + fit summary) or from a live re-fetch audit, not asserted from the + first run's output alone. + +**Status:** RUN complete; awaiting a PR to land these results. +Probe + docs — zero product code. + ## 2026-08-12 — lance-graph #932 (MERGED) — the golden-ratio index floor + the temperament mechanism (merged before its own review findings could be addressed — see #933) - **Added.** The operator-ruled golden-ratio index floor (convergent index diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index 50b36111..bf0d5c1b 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -27,7 +27,7 @@ Wave 1 = parallel, no operator gate beyond go-ahead; gated rows named. |---|---|---|---|---| | D-W5 | Spiral-ADI anisotropy (v2, full-band control, V-matched iterations, bump 3.35σ from mask) | 1 | **B2 FAIL / B3 VOID CONFIRMED / B4 INCOMPLETE** — B2: real diffusion resolved, aniso 1.5251 vs 1.25 bar, clean baseline 1.0046, operator contributes ~0.52. B3: family A 99.68 % + family B 99.56 % (both link families, QUALIFYING population n=4.78M out of the headline lattice N=7.65M, not the 62k sub-sample) land on a pure Fibonacci offset — dominated by the two discovered strides 2584=F(18)/4181=F(19) respectively. B4: downgraded to a DESCRIPTIVE reading over n=8–17 only — n=19 was dropped without pre-authorization, bar not satisfied, stays open | domino.rs gather-design claim REFUTED at this test point (v1's "unblocked" was the same inert-operator artifact that also drove B2's false PASS); `E-ON-A-GOLDEN-LATTICE-LOCALITY-IS-FIBONACCI-MEMBERSHIP-1` strengthened, third dated update; B4's n=19 remains an open follow-up | | D-W2sA | Golden-vs-grid pairing on real cos-lat metric (zero-ties G1, CV G2) | 1 | **RUN — G1 VOID, G2/G4 FAIL (control degenerate: two identical translated grids are symmetry-uniform, CV ~1e-12 — cannot lose any evenness comparison; diagnosed via smoke test, run as-specified, `E-A-CONTROL-THAT-CANNOT-LOSE-IS-NO-CONTROL-1`)** | honest falsifier for evenness DEFERRED (offset/rotation-varied or spacing-mismatched control); §10.5 properties 1–3 untouched | -| D-W6 | Two-component deconvolution (geo + bow, global lstsq, 38 eqs / 2 params; B3 = stranded stratification via v_rel) | 1 | Queued | dipole vector-sum identification; F17 gate | +| D-W6 | Two-component deconvolution (geo + bow, global lstsq, 38 eqs / 2 params; B3 = stranded stratification via v_rel) | 1 | **RUN COMPLETE — B0 VOID** (single-geo R²=−0.104, worse than the mean; both anti-vacuity controls score ≥ the ceiling — the model has nothing to identify); B1/B2 correctly report VOID per rule; B3 stranded stratum EMPTY (n=0 — a structural consequence of CT-F14's displacement≥250km filter implying `\|v_storm\|≥11.57 m/s`, not a physics finding) | vector-sum model DISCONFIRMED as specified on this sample; CT-F17 gate now moot for this model form — a revised model needs its own W6-shaped test first; `E-THE-DISPLACEMENT-FILTER-ATE-THE-STRANDED-STRATUM-1` | | D-W2sB | α-window sweep β∈[0.85,1.15] | gated (W2s-a) | Queued | corridor α discriminator | | D-W7 | Corridor two-regime α field probe | gated (W6) | Queued | §10.3 physics | | D-CT-F17 | FRESH-sample verdict (1959–1979, N=70 candidates, V-test p<0.05 ∧ R̄≥0.35; independent adversarial spec audit MANDATORY before bars commit) | gated (W6 + audit) | Queued | the directional claim's verdict path | diff --git a/.claude/plans/weather-w-probes-v1.md b/.claude/plans/weather-w-probes-v1.md index d7d67203..16c3764f 100644 --- a/.claude/plans/weather-w-probes-v1.md +++ b/.claude/plans/weather-w-probes-v1.md @@ -550,6 +550,56 @@ Checkpoint after every storm; resume skips completed `t0`. --- +### RUN, 2026-08-12 (`comet_tail_w6.py` / `.json`) — B0 VOID, and the vector-sum model does NOT fit; the stranded stratum is EMPTY by construction + +| bar | verdict | measured | +|---|---|---| +| **B0 CONTROLS** | **VOID** (own pre-registered rule) | single-geo R²=**−0.104**; permuted-P_bow control R²=**−0.071**; rotated-90° control R²=**−0.062** — both controls exceed `single-geo + 0.03 = −0.074` | +| **B1 IDENTIFIABILITY** | VOID (B0 failed) | joint R²=**−0.086** vs best-single R²=−0.104, margin +0.018 (needed +0.10) | +| **B2 SIGN** | VOID (B0 failed) | `c_geo = +0.407` (predicted positive — **correct sign**); `c_bow = −0.0006` km⁻¹ (predicted positive — **wrong sign**) | +| **B3 residual resultant** (descriptive) | overall: R̄=**0.153**, μ=−2.3°, p=0.646 — *below* the n=19 uniform-expectation floor (√π/2√19 ≈ 0.203); **stranded (<8 m/s): n=0**; moving (≥8 m/s): n=19, identical to overall | no clustering at all; stranded stratum **structurally empty** | + +**Every single-predictor model already fails on its own** — geo alone R²=−0.104, bow alone R²=−0.147, both *worse than predicting the mean*. The joint model's negligible +0.018 margin over the geo baseline is exactly what an anti-vacuity control should reject, and B0 correctly rejects it: two deliberately-wrong references (permuted, rotated 90°) score as well as or better than the real geo predictor alone. **There is no signal here for a joint fit to identify.** + +**No implementation bug found — checked, not assumed.** One storm was independently re-fetched and hand-audited (t0=54358); extended to the full 19-storm sample (committed in `comet_tail_w6.json`'s `B4_per_storm`, which now carries the raw `D`/`P_geo`/`P_bow`/`A_H`/`d_H`/`v_rel_ms` per storm, not only derived bearings — the first run shipped without these and the gap was self-caught before a reviewer needed to). Every value is physically sane (`A_H > 0` always, `d_H` inside the 600–2500 km annulus always, `v_rel` 2.8–27.9 m/s). + +> **⚠ TWO CORRECTIONS, 2026-08-12 (codex + CodeRabbit P2/Major on PR #940, +> both real, both fixed same day) — the SIGN and the UNITS were both wrong +> in the first draft above.** (1) **Sign:** `spine()`'s raw fit coefficient +> points toward the storm's HIGH side (the gradient of *increasing* +> residual pressure), while `P_geo`/`P_bow` both point toward the LOW side +> by construction — exactly the convention `low_pole_bearing()` makes +> explicit via its own `(ph + π) % (2π)` flip +> (`comet_tail_f16.py:138-160`). `D` is now `−spine(...)`, matching that +> convention. This flips `c_geo` from −0.407 to **+0.407** (now the +> physically PREDICTED positive sign) and `c_bow` from +0.0006 to +> **−0.0006** (now the WRONG sign, where it had looked merely near-zero +> before) — the table above already carries the corrected values. +> Verified algebraically AND numerically before/after the flip: R² and +> every B0/B1 verdict are provably unchanged by this sign convention +> (OLS is odd-symmetric in the fit target), confirmed bit-identical on +> rerun. (2) **Units:** `D` and `P_geo` are both [Pa/km]; `P_bow` is [Pa] +> — so `c_geo` is dimensionless but `c_bow` carries **km⁻¹**, and OLS +> coefficients rescale inversely under column rescaling while leaving +> R²/fitted values unchanged. **Raw `|c_bow|` was never a valid basis for +> "no measurable weight"** — the original `|P_bow|` averaging 147× `|D|` +> comparison compounded the same error (Pa vs Pa/km, not comparable at +> all). The dimensionally valid measure is the FITTED CONTRIBUTION, +> `|c_bow·P_bow|` against `|D|`, both in Pa/km: +> **mean `|D|` = 0.745, mean `|c_geo·P_geo|` = 0.186 (25 % of `|D|`), mean +> `|c_bow·P_bow|` = 0.068 (9 % of `|D|`)**. The geo contribution is ~2.7× +> the bow contribution in fitted terms — modest, not "no weight" — and +> both are consistent with the R²<0 finding that neither predictor +> meaningfully explains `D`'s variance. **B0/B1's VOID verdicts are +> untouched by either correction** — this changes only the supporting +> narrative around B2, not the headline finding. + +**The stranded stratum is empty for a structural reason, not a physics finding.** `min(v_storm) = 12.54 m/s` across all 19 storms — comfortably above `250 km / 6 h = 11.574 m/s`, the speed CT-F14's own `displacement_km ≥ 250` qualifying filter mathematically implies as a floor. **No storm in this displacement-filtered 19-storm set can ever be "stranded" (<8 m/s) — the report §10.2 stranded-rescue reading is UNTESTABLE on this sample by construction, not refuted.** Testing it needs a sample built WITHOUT the fast-motion-selecting displacement filter (or with a filter that explicitly retains slow storms) — a design note for any future stranded-rescue probe, not a task for this one to retrofit. + +**What this means for the report's vector-sum model.** As specified — a global 2-parameter linear combination of a single background-high neighbor predictor and a single relative-motion bow-wave predictor, fit across 19 storms by ordinary least squares — **the model does not fit this data, and the fit is not merely weak, it is void by its own anti-vacuity control.** This does not rule out a richer version of the model (multiple neighbors, a nonlinear bow term, storm-specific coefficients) — but the specific, pre-registered, mechanistically-motivated form named in §10.2 is disconfirmed on this sample as tested. Consequence for CT-F17: its gate ("W6's result AND an independent adversarial spec audit") is now moot for the vector-sum model's CURRENT form — a fresh-sample verdict on a model that already fails its identifiability control on the STORED sample would not be a meaningful next step; the audit gate stands for any REVISED form of the model instead. + +--- + ## §4 BRIEF W2s-b — the α-field on a real H–T pair (GATED on W2s-a G2 pass) **File:** `corridor_alpha_probe.py`. Outline — finalize bars at spawn time diff --git a/probes/weather-p1/comet_tail_w6.json b/probes/weather-p1/comet_tail_w6.json new file mode 100644 index 00000000..cf012ce4 --- /dev/null +++ b/probes/weather-p1/comet_tail_w6.json @@ -0,0 +1,504 @@ +{ + "n": 19, + "seed": 20260812, + "B0_controls": { + "r2_single_geo": -0.10381057418374096, + "r2_permuted": -0.07139097664756067, + "r2_rotated90": -0.062400718434540536, + "bar": "<= r2_single_geo + 0.03", + "verdict": "VOID" + }, + "single_models": { + "geo": { + "c": 0.4361308247679559, + "R2": -0.10381057418374096 + }, + "bow": { + "c": -0.0007731317854461328, + "R2": -0.14654469146937688 + } + }, + "joint_model": { + "c_geo": 0.40656942616866554, + "c_bow": -0.0006195037332984975, + "R2": -0.08608330921952212 + }, + "B1_identifiability": { + "joint_R2": -0.08608330921952212, + "best_single_R2": -0.10381057418374096, + "margin": 0.017727264964218836, + "bar": ">= best_single_R2 + 0.10", + "verdict": "VOID (B0 failed)" + }, + "B2_sign": { + "c_geo": 0.40656942616866554, + "c_bow": -0.0006195037332984975, + "units_note": "c_geo dimensionless (both D and P_geo are Pa/km); c_bow has units km^-1 (D is Pa/km, P_bow is Pa) -- see fitted_contribution_Pa_per_km for a unit-consistent magnitude comparison, not raw |c_bow|", + "verdict": "VOID (B0 failed)" + }, + "fitted_contribution_Pa_per_km": { + "mean_|D|": 0.7454893648206984, + "mean_|c_geo*P_geo|": 0.18568652512701836, + "mean_|c_bow*P_bow|": 0.06793644728208614, + "geo_contrib_frac_of_D": 0.24908004579204, + "bow_contrib_frac_of_D": 0.09113000196646119 + }, + "B3_residual_resultant": { + "overall": { + "R_bar": 0.15311408050017475, + "mu_deg": -2.3016039763941936, + "rayleigh_p": 0.6464383782088862 + }, + "stranded_lt8ms": { + "n": 0, + "R_bar": null, + "mu_deg": null, + "rayleigh_p": null + }, + "moving_ge8ms": { + "n": 19, + "R_bar": 0.15311408050017475, + "mu_deg": -2.3016039763941936, + "rayleigh_p": 0.6464383782088862 + } + }, + "B4_per_storm": [ + { + "t0": 54358, + "bearing_D_deg": 39.85050154139167, + "bearing_Dhat_deg": 135.70001638914562, + "resid_deg": -95.84951484775394, + "v_storm_ms": 18.28531366893943, + "D": [ + 0.5614133833389219, + 0.46859095586174915 + ], + "P_geo": [ + -0.33964208969670856, + 0.23720740369675417 + ], + "P_bow": [ + -43.63373173859991, + -19.264582561146636 + ], + "v_rel_ms": 8.916019274557227, + "A_H_Pa": 905.9376302083401, + "d_H_km": 2186.8004833107357, + "theta_H_rad": -0.6096540392824519 + }, + { + "t0": 55578, + "bearing_D_deg": 95.62746477365677, + "bearing_Dhat_deg": 177.7683544951107, + "resid_deg": -82.14088972145395, + "v_storm_ms": 12.783884621110882, + "D": [ + -0.07928538759450642, + 0.8046431920973721 + ], + "P_geo": [ + -0.5521970281293828, + -0.004784604144652332 + ], + "P_bow": [ + -13.186772014911245, + -16.74852788286718 + ], + "v_rel_ms": 5.960530920810273, + "A_H_Pa": 1372.707899305562, + "d_H_km": 2485.809055907479, + "theta_H_rad": 0.008664451580220747 + }, + { + "t0": 55822, + "bearing_D_deg": 161.09781086412593, + "bearing_Dhat_deg": -151.05468982925726, + "resid_deg": -47.84749930661678, + "v_storm_ms": 15.22434719417646, + "D": [ + -0.8696922749351309, + 0.2977993440558817 + ], + "P_geo": [ + -0.44079578715172685, + -0.26674607644830617 + ], + "P_bow": [ + -10.104461592956396, + -20.65579555862832 + ], + "v_rel_ms": 6.190695471090927, + "A_H_Pa": 1286.806157769097, + "d_H_km": 2497.5728813709798, + "theta_H_rad": 0.5441951816750085 + }, + { + "t0": 59726, + "bearing_D_deg": 50.99573280407391, + "bearing_Dhat_deg": 64.01295378469953, + "resid_deg": -13.01722098062561, + "v_storm_ms": 13.706197408730237, + "D": [ + 0.7443119031528372, + 0.9190086968281168 + ], + "P_geo": [ + 0.1425898473700703, + 0.3399739707305417 + ], + "P_bow": [ + -15.275801356096554, + -0.1950334026701901 + ], + "v_rel_ms": 5.0459631302903825, + "A_H_Pa": 778.6802951388818, + "d_H_km": 2112.160048025693, + "theta_H_rad": -1.967926095263975 + }, + { + "t0": 62410, + "bearing_D_deg": -31.620883758388153, + "bearing_Dhat_deg": -163.50778459972932, + "resid_deg": 131.88690084134117, + "v_storm_ms": 12.576298180205079, + "D": [ + 0.13429403865980433, + -0.08268573381842366 + ], + "P_geo": [ + -0.675225328363115, + -0.23915668565414538 + ], + "P_bow": [ + -42.683103073968994, + -38.39321327239433 + ], + "v_rel_ms": 9.7817694813837, + "A_H_Pa": 1108.7235405816027, + "d_H_km": 1547.7884868103886, + "theta_H_rad": 0.34040086396685526 + }, + { + "t0": 62898, + "bearing_D_deg": 95.54528377998422, + "bearing_Dhat_deg": -97.99125028597646, + "resid_deg": -166.46346593403933, + "v_storm_ms": 14.943193389122174, + "D": [ + -0.16679390712012496, + 1.717987696028101 + ], + "P_geo": [ + -0.09217767395646605, + -0.6102376442820587 + ], + "P_bow": [ + -4.433273636366096, + -1.1487830163504862 + ], + "v_rel_ms": 2.7627571462209923, + "A_H_Pa": 1215.8246744791686, + "d_H_km": 1970.0309330972063, + "theta_H_rad": 1.4208775956251016 + }, + { + "t0": 63142, + "bearing_D_deg": 54.26196494243309, + "bearing_Dhat_deg": 177.05826003046167, + "resid_deg": -122.79629508802859, + "v_storm_ms": 12.717183891649741, + "D": [ + 0.28622693783421815, + 0.39776945504954614 + ], + "P_geo": [ + -0.6827556516666856, + -0.04846924662396028 + ], + "P_bow": [ + 2.385182678874679, + -54.95814395759871 + ], + "v_rel_ms": 9.575130813332041, + "A_H_Pa": 1610.3618489583314, + "d_H_km": 2352.700082338784, + "theta_H_rad": 0.07087172234858281 + }, + { + "t0": 63630, + "bearing_D_deg": -78.45913853725781, + "bearing_Dhat_deg": 160.28153426999518, + "resid_deg": 121.25932719274704, + "v_storm_ms": 21.18286602908681, + "D": [ + 0.05306805970771929, + -0.25988900401614007 + ], + "P_geo": [ + 0.12551788402042788, + -0.15527935773518797 + ], + "P_bow": [ + 344.9480162814228, + -196.01728584720016 + ], + "v_rel_ms": 25.71483678511945, + "A_H_Pa": 496.8787651909661, + "d_H_km": 2488.5526193372934, + "theta_H_rad": 2.2505995249627646 + }, + { + "t0": 65094, + "bearing_D_deg": 43.39450535929249, + "bearing_Dhat_deg": -13.340135138234208, + "resid_deg": 56.7346404975267, + "v_storm_ms": 27.39969823490195, + "D": [ + 0.22978891893623596, + 0.21725883912867677 + ], + "P_geo": [ + 0.21108343052667497, + -0.13400480563377204 + ], + "P_bow": [ + -351.86889705178135, + 28.34324399352165 + ], + "v_rel_ms": 24.25587845374545, + "A_H_Pa": 624.8040364583285, + "d_H_km": 2498.9462227237236, + "theta_H_rad": 2.575946627171843 + }, + { + "t0": 66558, + "bearing_D_deg": -179.03522695431383, + "bearing_Dhat_deg": -11.105709069986547, + "resid_deg": -167.92951788432728, + "v_storm_ms": 19.88968241564589, + "D": [ + -0.14592651753919264, + -0.002457410988132297 + ], + "P_geo": [ + 0.05500058563147917, + -0.007527304372467833 + ], + "P_bow": [ + -133.50704523448445, + 28.35229134318727 + ], + "v_rel_ms": 15.08223961537412, + "A_H_Pa": 138.4839518229128, + "d_H_km": 2494.6092387290937, + "theta_H_rad": 3.0055790222094045 + }, + { + "t0": 67046, + "bearing_D_deg": 162.2070290972292, + "bearing_Dhat_deg": 176.5434613345046, + "resid_deg": -14.336432237275403, + "v_storm_ms": 14.194808034891913, + "D": [ + -1.7286903849258044, + 0.5547878337044447 + ], + "P_geo": [ + -0.776482785822735, + 0.04060590531375081 + ], + "P_bow": [ + 12.683636198354918, + -4.897172452210021 + ], + "v_rel_ms": 4.7602888495379085, + "A_H_Pa": 1241.671636284722, + "d_H_km": 1596.9153649016405, + "theta_H_rad": -0.052247070061416344 + }, + { + "t0": 67778, + "bearing_D_deg": 119.5090522278074, + "bearing_Dhat_deg": -8.976294463196243, + "resid_deg": 128.48534669100366, + "v_storm_ms": 14.226193347487772, + "D": [ + -0.23261607945264737, + 0.4109960072201622 + ], + "P_geo": [ + 0.6017113247733159, + -0.15495481990831922 + ], + "P_bow": [ + -49.88844047158462, + -31.436362856799295 + ], + "v_rel_ms": 9.913538872090644, + "A_H_Pa": 1547.5926052517316, + "d_H_km": 2490.7206192405865, + "theta_H_rad": 2.889545658644649 + }, + { + "t0": 68754, + "bearing_D_deg": -138.16172056958374, + "bearing_Dhat_deg": -43.07068062555434, + "resid_deg": -95.0910399440294, + "v_storm_ms": 13.49216615367673, + "D": [ + -0.43903176313967673, + -0.3930678184962927 + ], + "P_geo": [ + 0.33367620197461706, + -0.3329477688658759 + ], + "P_bow": [ + 0.3719453021337218, + -14.142248613849121 + ], + "v_rel_ms": 4.855776440455581, + "A_H_Pa": 1177.317822265628, + "d_H_km": 2497.626713528983, + "theta_H_rad": 2.357287209402393 + }, + { + "t0": 68998, + "bearing_D_deg": 38.238459044191984, + "bearing_Dhat_deg": -28.744066021246514, + "resid_deg": 66.9825250654385, + "v_storm_ms": 13.07843632631205, + "D": [ + 0.745057352396505, + 0.5871125523414068 + ], + "P_geo": [ + 0.3540219374953975, + -0.24266491111920693 + ], + "P_bow": [ + -25.016829546769788, + -18.101490941298547 + ], + "v_rel_ms": 7.173898077303566, + "A_H_Pa": 1006.0542046440969, + "d_H_km": 2343.989138641649, + "theta_H_rad": 2.540697579284103 + }, + { + "t0": 70218, + "bearing_D_deg": -158.90278612887818, + "bearing_Dhat_deg": -102.15095322073547, + "resid_deg": -56.7518329081427, + "v_storm_ms": 14.70736816777759, + "D": [ + -0.33346941327299046, + -0.12865650778683624 + ], + "P_geo": [ + -0.07728779615182074, + -0.12649048272525307 + ], + "P_bow": [ + -30.98217740937703, + 8.669595937522367 + ], + "v_rel_ms": 7.3226025686378495, + "A_H_Pa": 360.6112575954903, + "d_H_km": 2432.720276365064, + "theta_H_rad": 1.0223156738251673 + }, + { + "t0": 73634, + "bearing_D_deg": -103.83661463453105, + "bearing_Dhat_deg": -119.15707329062008, + "resid_deg": 15.320458656089045, + "v_storm_ms": 27.084474843184424, + "D": [ + -0.17484039563323409, + -0.7098639772091326 + ], + "P_geo": [ + 0.03450959494694185, + -0.4976758586015638 + ], + "P_bow": [ + 366.92949631831254, + 290.4884037071427 + ], + "v_rel_ms": 27.928377054816018, + "A_H_Pa": 1051.943663194441, + "d_H_km": 2108.6490882789735, + "theta_H_rad": 1.6400270181744452 + }, + { + "t0": 73878, + "bearing_D_deg": -111.93148785397005, + "bearing_Dhat_deg": -124.06327095763608, + "resid_deg": 12.13178310366601, + "v_storm_ms": 16.81493066406619, + "D": [ + -0.18760054245727764, + -0.4659309395059687 + ], + "P_geo": [ + -0.22448508252419622, + -0.25271388797470856 + ], + "P_bow": [ + 57.73026340585698, + 137.43318896729545 + ], + "v_rel_ms": 15.762083960308457, + "A_H_Pa": 394.9323947482626, + "d_H_km": 1168.3681596531612, + "theta_H_rad": 0.8444845168265731 + }, + { + "t0": 74122, + "bearing_D_deg": 119.35921843918658, + "bearing_Dhat_deg": 68.46824507880542, + "resid_deg": 50.89097336038117, + "v_storm_ms": 26.05788702271268, + "D": [ + -0.6080720956923845, + 1.080952379830265 + ], + "P_geo": [ + 0.03121803878898213, + 0.5574998389691659 + ], + "P_bow": [ + -152.55053729021833, + -72.69321437555732 + ], + "v_rel_ms": 16.782188833487737, + "A_H_Pa": 1166.9645724826405, + "d_H_km": 2089.9365539158202, + "theta_H_rad": -1.6267344072817411 + }, + { + "t0": 74366, + "bearing_D_deg": 70.31914944701053, + "bearing_Dhat_deg": 16.221232295960224, + "resid_deg": 54.0979171510503, + "v_storm_ms": 12.537805236919697, + "D": [ + 0.25752297964948623, + 0.7199919298029541 + ], + "P_geo": [ + 0.44272351166787993, + 0.12879728832936585 + ], + "P_bow": [ + -35.76846403339691, + -10.40850969356574 + ], + "v_rel_ms": 7.87952148366395, + "A_H_Pa": 1151.9230197482684, + "d_H_km": 2498.3261336095843, + "theta_H_rad": -2.858486468225491 + } + ] +} \ No newline at end of file diff --git a/probes/weather-p1/comet_tail_w6.py b/probes/weather-p1/comet_tail_w6.py new file mode 100644 index 00000000..4b131950 --- /dev/null +++ b/probes/weather-p1/comet_tail_w6.py @@ -0,0 +1,445 @@ +"""W6 -- the dipole deconvolution: is D = c_geo*neighbor + c_bow*bow-wave? + +Per weather-w-probes-v1.md §3 (the report §10.2 vector-sum model). This is a +MECHANISTIC test on CT-F14's 19 stored storms, NOT a verdict -- a fresh- +sample verdict is CT-F17, gated on this result plus an independent +adversarial spec audit. + +Every geometry/statistics primitive below is copied VERBATIM from the +existing arc scripts, named at each site, so this probe carries no new +scoring convention: `wrap_deg`/`err_deg`/`geom_ll`/`disk_mean_uv` from +comet_tail_f16.py; `spine()` -- the CONSTRAINED 2-parameter dipole fit +(ring-profile means + a single global lstsq on [r*cos(theta), r*sin(theta)]) +-- from l4_rail_probe.py (the brief names this file as comet_tail_f16.py; +verified against the tree it actually lives in l4_rail_probe.py, ported from +there instead); `circular()` from the report §10.1 statistics standard. + +UNITS -- corrected 2026-08-12 (codex + CodeRabbit P2/Major on PR #940, both +real, both fixed here). `spine()` regresses pressure [Pa] onto +`[r*cos(theta), r*sin(theta)]` in KM, so `D = (a1, b1)` has units [Pa/km], +NOT dimensionless Pa as the first draft's docstring claimed -- and P_geo +(`A_H/d_H`) is ALSO [Pa/km], so `D = c_geo*P_geo + c_bow*P_bow` forces +`c_geo` DIMENSIONLESS but `c_bow` to carry [Pa/km]/[Pa] = **km^-1**. +Rescaling P_bow's raw units would rescale c_bow's NUMERIC VALUE inversely +while leaving R^2 and the fitted Dhat UNCHANGED (ordinary least squares is +scale-covariant per column) -- so a raw `|c_bow|` close to zero does NOT by +itself prove "no measurable weight"; the dimensionally valid comparison is +the FITTED CONTRIBUTION magnitude `|c_bow * P_bow|` against `|D|`, computed +below and reported instead of the raw coefficient's magnitude. + +SIGN CONVENTION -- also corrected. `spine()`'s raw `coef` is the gradient of +INCREASING residual (`resid = v - prof[rings]`), i.e. it points toward the +storm's HIGH side, not the low pole -- exactly the convention +`low_pole_bearing()` in comet_tail_f16.py:138-160 makes explicit by +returning `(ph + pi) % (2*pi)`, an explicit 180-degree flip AFTER computing +the raw regression-direction angle `ph`. `P_geo`/`P_bow` are both +constructed pointing toward the LOW side (away from the neighbor high; +behind relative motion) -- so comparing them meaningfully to `spine()`'s +raw `coef` needs the SAME flip. `D = -spine(...)` below, not the raw coef. +""" +import datetime +import json +import pathlib +import urllib.request + +import numcodecs +import numpy as np + +SEED = 20260812 +B = ("https://storage.googleapis.com/weatherbench2/datasets/era5/" + "1959-2022-6h-1440x721.zarr") +R_E, R_DISK, RING = 6371.0, 1200.0, 100.0 +ANNULUS_LO, ANNULUS_HI = 600.0, 2500.0 +RHO_AIR = 1.2 # kg/m^3, sea-level air density -- the bow-wave dynamic-pressure constant + +op = urllib.request.build_opener(urllib.request.ProxyHandler({})) +meta = json.loads(op.open(B + "/.zmetadata", timeout=90).read())["metadata"] + +EPOCH = datetime.datetime(1959, 1, 1) + + +def t_index(dt): + """WB2 time index: 6-hourly steps since 1959-01-01. Anchor guard, run once + at import time so a broken store fails loudly before any storm fetch.""" + return int(round((dt - EPOCH).total_seconds() / 3600 / 6)) + + +assert t_index(datetime.datetime(2021, 6, 15, 12)) == 91246 + + +def fetch(var, key): + """Fetch and decode one zarr chunk from the WB2 store.""" + za = meta[f"{var}/.zarray"] + raw = op.open(f"{B}/{var}/{key}", timeout=900).read() + dec = numcodecs.get_codec(za["compressor"]).decode(raw) + return np.frombuffer(dec, dtype=np.dtype(za["dtype"])).reshape(za["chunks"]) + + +def wrap_deg(d): + """Wrap degrees into [-180, 180) -- verbatim from comet_tail_f16.py.""" + return (d + 180.0) % 360.0 - 180.0 + + +def err_deg(low_pole_rad, motion_rad): + """Signed alignment error vs the left-of-motion prediction. Verbatim from + comet_tail_f16.py's err_deg -- kept only for the algebraic-inversion step + below, not used to re-score anything here.""" + return float(wrap_deg(np.rad2deg(low_pole_rad - (motion_rad + np.pi / 2)))) + + +lat = fetch("latitude", "0").astype(np.float64).ravel() +levels = fetch("level", "0").astype(int).ravel() +NY, NX = lat.size, 1440 +phi = np.deg2rad(lat) +lon_deg = np.arange(NX) * 0.25 +LEV_IDX = {int(v): i for i, v in enumerate(levels)} +T_MAX = meta["mean_sea_level_pressure/.zarray"]["shape"][0] - 1 + + +def geom_ll(latc, lonc): + """dx, dy, r (km), azimuth (rad CCW from east) about a continuous centre. + Verbatim from comet_tail_f16.py -- the SAME basis spine()'s coef and + disk_mean_uv's (u,v) both live in, so D/P_geo/P_bow are directly + combinable without a basis-change step.""" + phic = np.deg2rad(latc) + dlon = np.deg2rad((lon_deg[None, :] - lonc + 180) % 360 - 180) + dx = R_E * np.cos(phic) * dlon * np.ones((NY, 1)) + dy = R_E * (phi[:, None] - phic) * np.ones((1, NX)) + return dx, dy, np.hypot(dx, dy), np.arctan2(dy, dx) + + +def spine(la, lo, p0): + """The f64 constrained spine: ring-profile means + the CONSTRAINED + 2-parameter dipole (a single global lstsq on [r*cos(theta), r*sin(theta)] + against the ring-demeaned residual). Ported verbatim from + l4_rail_probe.py:109-122 (the brief names comet_tail_f16.py; the actual + function lives in l4_rail_probe.py -- verified against the tree, ported + from the real location rather than guessed). Returns D = coef = (a1, b1), + the dipole vector this whole probe deconvolves.""" + _, _, r, th = geom_ll(la, lo) + disk = r <= R_DISK + v, rr, tt = p0[disk], r[disk], th[disk] + nb = int(R_DISK / RING) + rings = np.clip((rr / RING).astype(int), 0, nb - 1) + prof = np.array([v[rings == b].mean() if (rings == b).any() else 0.0 + for b in range(nb)]) + resid = v - prof[rings] + X = np.column_stack([rr * np.cos(tt), rr * np.sin(tt)]) + coef, *_ = np.linalg.lstsq(X, resid, rcond=None) + return coef + + +def disk_mean_uv(u3, v3, latc, lonc, lev_list): + """Disk-mean (u, v) over the 1200 km disk, averaged across `lev_list`. + Verbatim from comet_tail_f16.py.""" + _, _, r, _ = geom_ll(latc, lonc) + disk = r <= R_DISK + us, vs = [], [] + for lev in lev_list: + li = LEV_IDX[lev] + us.append(u3[li][disk].mean()) + vs.append(v3[li][disk].mean()) + return float(np.mean(us)), float(np.mean(vs)) + + +def neighbor_predictor(p0, latc, lonc): + """Strongest POSITIVE zonal-anomaly cell in the 600-2500 km annulus, + lat 20-80N, about the storm centre -- the background-high the report's + §10.2 vector-sum model treats as a far-field neighbour. Returns + (A_H [Pa], d_H [km], theta_H [rad]); None if the annulus admits no + candidate (should not occur inside the storm's climatological band, but + checked rather than assumed -- see the NO-VERDICT path in the main loop).""" + fa = p0 - p0.mean(axis=1, keepdims=True) + _, _, r, th = geom_ll(latc, lonc) + lat_ok = (lat >= 20) & (lat <= 80) + mask = (r >= ANNULUS_LO) & (r <= ANNULUS_HI) & lat_ok[:, None] + cand = np.where(mask, fa, -np.inf) + if not np.isfinite(cand).any() or np.nanmax(cand) <= 0: + return None + i, j = np.unravel_index(np.argmax(cand), cand.shape) + return float(fa[i, j]), float(r[i, j]), float(th[i, j]) + + +def circular(errs_deg, n): + """Resultant length R_bar, mean direction mu (deg), Rayleigh p + (Zar/Mardia small-n correction) -- report §10.1's binding statistics + standard, verbatim.""" + th = np.deg2rad(np.asarray(errs_deg)) + c, s = np.cos(th).mean(), np.sin(th).mean() + r = float(np.hypot(c, s)) + mu = float(np.rad2deg(np.arctan2(s, c))) + z = n * r * r + p = float(np.exp(-z) * (1 + (2 * z - z * z) / (4 * n) + - (24 * z - 132 * z**2 + 76 * z**3 - 9 * z**4) + / (288 * n * n))) + return r, mu, max(min(p, 1.0), 0.0) + + +def r2_vec(D, Dhat): + """Vector R^2: 1 - sum|D_i-Dhat_i|^2 / sum|D_i-mean(D)|^2, summed over + BOTH x/y components and all storms -- the identifiability metric B0/B1 + are scored on.""" + D, Dhat = np.asarray(D), np.asarray(Dhat) + sse = float(np.sum((D - Dhat) ** 2)) + sst = float(np.sum((D - D.mean(axis=0)) ** 2)) + return 1.0 - sse / sst if sst > 0 else float("nan") + + +def fit_joint(D, Pg, Pb): + """Solve D = c_geo*Pg + c_bow*Pb by lstsq over the 2N-stacked scalar + equations (38 for 19 storms); returns (c_geo, c_bow, Dhat, R2).""" + D, Pg, Pb = np.asarray(D), np.asarray(Pg), np.asarray(Pb) + y = D.ravel() + X = np.column_stack([Pg.ravel(), Pb.ravel()]) + coef, *_ = np.linalg.lstsq(X, y, rcond=None) + Dhat = (coef[0] * Pg + coef[1] * Pb) + return float(coef[0]), float(coef[1]), Dhat, r2_vec(D, Dhat) + + +def fit_single(D, P): + """Solve D = c*P by lstsq; returns (c, Dhat, R2) for the single-predictor + comparison B1 needs.""" + D, P = np.asarray(D), np.asarray(P) + y = D.ravel() + X = P.ravel()[:, None] + coef, *_ = np.linalg.lstsq(X, y, rcond=None) + Dhat = coef[0] * P + return float(coef[0]), Dhat, r2_vec(D, Dhat) + + +def bearing_deg(vec2): + """arctan2(y,x) in degrees, the same east-CCW convention geom_ll/ + disk_mean_uv use throughout this probe.""" + return float(np.rad2deg(np.arctan2(vec2[1], vec2[0]))) + + +# ---- CT-F14's 19 qualifying storms, read directly from comet_tail_f16.json +# (the brief's own input spec) -- reuses ALL of CT-F16's already-stored +# fields, so nothing about centre-finding, displacement filtering, or the +# steering-level scoring can move between this probe and the prior ones. +src = json.loads( + pathlib.Path(__file__).with_name("comet_tail_f16.json").read_text()) +storms = src["rows"] +assert len(storms) == 19, f"expected CT-F14/F16's 19 qualifying storms, got {len(storms)}" + +out_dir = pathlib.Path(__file__).parent +partial_path = out_dir.with_name("weather-p1") / "comet_tail_w6.partial.jsonl" +tag_path = out_dir / "exec-runs" / "comet_tail_w6.txt" +tag_path.parent.mkdir(exist_ok=True) + + +def load_completed(): + """Resume-skip: read the partial checkpoint file if it exists, return + {t0: row} for every already-fetched storm so a stranded run resumes + instead of re-fetching.""" + done = {} + if partial_path.exists(): + for line in partial_path.read_text().splitlines(): + if line.strip(): + row = json.loads(line) + done[row["t0"]] = row + return done + + +def run(): + """Per-storm fetch + spine/neighbor/bow computation (checkpointed), then + the global 38-equation joint fit and all four pre-registered bars, + B0 (controls, reported FIRST) through B4.""" + # Load the checkpoint BEFORE writing any run metadata (codex/CodeRabbit + # P2 on #940): if load_completed() itself fails on a malformed partial + # file, the tag file must not already claim a run started successfully. + done = load_completed() + with open(tag_path, "a") as tf: + tf.write(f"START seed={SEED} n_storms={len(storms)} " + f"resumed={len(done)}\n") + + with open(partial_path, "a") as pf, open(tag_path, "a") as tf: + for i, s in enumerate(storms): + t0 = s["t0"] + if t0 in done: + continue + if t0 > T_MAX: + tf.write(f"SKIP t0={t0} beyond store coverage (T_MAX={T_MAX})\n") + continue + try: + la, lo = s["center_lat"], s["center_lon"] + p0 = fetch("mean_sea_level_pressure", f"{t0}.0.0")[0].astype(np.float64) + u3 = fetch("u_component_of_wind", f"{t0}.0.0.0")[0].astype(np.float64) + v3 = fetch("v_component_of_wind", f"{t0}.0.0.0")[0].astype(np.float64) + + # Sign-flipped per the module docstring's SIGN CONVENTION + # note: spine()'s raw coef points toward the HIGH side + # (increasing residual); P_geo/P_bow both point toward the + # LOW side, the SAME convention low_pole_bearing()'s "+pi" + # flip encodes. + D = -spine(la, lo, p0) + + # Step 2: motion bearing recovered by ALGEBRA (exact + # inversion of err_deg = wrap(lp - (mth+pi/2))), no tracking. + mth_deg = wrap_deg(np.rad2deg(s["low_pole_rad"]) - 90.0 + - s["err_surface_deg"]) + mth_rad = np.deg2rad(mth_deg) + v_storm_ms = s["displacement_km"] * 1000.0 / (6 * 3600.0) + v_storm = v_storm_ms * np.array( + [np.cos(mth_rad), np.sin(mth_rad)]) + + # Step 3: v_rel = v_storm - v_env850; bow predictor. + u850, v850 = disk_mean_uv(u3, v3, la, lo, (850,)) + v_env850 = np.array([u850, v850]) + v_rel = v_storm - v_env850 + speed_rel = float(np.hypot(*v_rel)) + bear_rel = float(np.arctan2(v_rel[1], v_rel[0])) + P_bow = (0.5 * RHO_AIR * speed_rel ** 2) * np.array( + [np.cos(bear_rel + np.pi), np.sin(bear_rel + np.pi)]) + + # Step 4: neighbor predictor. + nb = neighbor_predictor(p0, la, lo) + if nb is None: + tf.write(f"NO-VERDICT t0={t0}: no positive annulus anomaly\n") + tf.flush() + continue + A_H, d_H, theta_H = nb + P_geo = (A_H / d_H) * np.array( + [np.cos(theta_H + np.pi), np.sin(theta_H + np.pi)]) + except Exception as exc: + # §0's iron rule: a dead fetch/step is recorded, then the + # run stops -- it does not improvise or silently skip. + tf.write(f"ERROR t0={t0}: {type(exc).__name__}: {exc}\n") + tf.flush() + raise + + row = {"t0": t0, "date": s["date"], "D": D.tolist(), + "P_geo": P_geo.tolist(), "P_bow": P_bow.tolist(), + "v_storm_ms": v_storm_ms, "v_rel_ms": speed_rel, + "A_H_Pa": A_H, "d_H_km": d_H, "theta_H_rad": theta_H} + pf.write(json.dumps(row) + "\n") + pf.flush() + done[t0] = row + if (i + 1) % 5 == 0 or (i + 1) == len(storms): + tf.write(f"progress {i + 1}/{len(storms)} t0={t0}\n") + tf.flush() + + rows = [done[s["t0"]] for s in storms if s["t0"] in done] + n = len(rows) + D = np.array([r["D"] for r in rows]) + Pg = np.array([r["P_geo"] for r in rows]) + Pb = np.array([r["P_bow"] for r in rows]) + vstorm = np.array([r["v_storm_ms"] for r in rows]) + + rng = np.random.default_rng(SEED) + + # B0 -- controls FIRST, reported before any real-model number. + Pb_perm = Pb[np.array([(i + 7) % n for i in range(n)])] + _, _, _, r2_perm = fit_joint(D, Pg, Pb_perm) + rot = np.array([[0.0, -1.0], [1.0, 0.0]]) # +90deg rotation matrix + Pb_rot = Pb @ rot.T + _, _, _, r2_rot = fit_joint(D, Pg, Pb_rot) + + c_geo_single, Dhat_geo, r2_geo = fit_single(D, Pg) + c_bow_single, Dhat_bow, r2_bow = fit_single(D, Pb) + best_single = max(r2_geo, r2_bow) + + b0_pass = (r2_perm <= r2_geo + 0.03) and (r2_rot <= r2_geo + 0.03) + + c_geo, c_bow, Dhat_joint, r2_joint = fit_joint(D, Pg, Pb) + + b1_pass = r2_joint >= best_single + 0.10 + b2_pass = (c_bow > 0) and (c_geo > 0) + + # DIMENSIONALLY VALID contribution comparison (codex/CodeRabbit P2 on + # #940): D and c_geo*P_geo are both [Pa/km] and directly comparable; + # c_bow*P_bow is ALSO [Pa/km] once the coefficient's implicit km^-1 is + # applied, so |c_bow*P_bow| vs |D| is the right ratio -- NOT raw + # |c_bow| (dimensionless-looking but actually km^-1) vs anything, and + # NOT raw |P_bow| vs |D| (different units, Pa vs Pa/km, incomparable). + geo_contrib = c_geo * Pg + bow_contrib = c_bow * Pb + mean_D_mag = float(np.mean(np.hypot(D[:, 0], D[:, 1]))) + mean_geo_contrib_mag = float(np.mean(np.hypot(geo_contrib[:, 0], geo_contrib[:, 1]))) + mean_bow_contrib_mag = float(np.mean(np.hypot(bow_contrib[:, 0], bow_contrib[:, 1]))) + + resid_bear = np.array([wrap_deg(bearing_deg(D[i]) - bearing_deg(Dhat_joint[i])) + for i in range(n)]) + rbar_all, mu_all, p_all = circular(resid_bear, n) + stranded = vstorm < 8.0 + moving = ~stranded + rbar_s, mu_s, p_s = (circular(resid_bear[stranded], int(stranded.sum())) + if stranded.sum() >= 2 else (None, None, None)) + rbar_m, mu_m, p_m = (circular(resid_bear[moving], int(moving.sum())) + if moving.sum() >= 2 else (None, None, None)) + + # per_storm carries the RAW predictors too, not only derived bearings -- + # a magnitude/units audit (are P_geo/P_bow physically sane, is there a + # scale mismatch driving a near-zero coefficient) needs the actual + # vectors committed, not just the fit's summary numbers. (Fixed after + # the first run shipped only bearings; re-run rather than leaving the + # audit gap, since a rerun costs the same ~3 min the brief itself + # estimates.) + per_storm = [{"t0": rows[i]["t0"], "bearing_D_deg": bearing_deg(D[i]), + "bearing_Dhat_deg": bearing_deg(Dhat_joint[i]), + "resid_deg": float(resid_bear[i]), + "v_storm_ms": float(vstorm[i]), + "D": D[i].tolist(), "P_geo": Pg[i].tolist(), + "P_bow": Pb[i].tolist(), "v_rel_ms": rows[i]["v_rel_ms"], + "A_H_Pa": rows[i]["A_H_Pa"], "d_H_km": rows[i]["d_H_km"], + "theta_H_rad": rows[i]["theta_H_rad"]} for i in range(n)] + + out = { + "n": n, "seed": SEED, + "B0_controls": {"r2_single_geo": r2_geo, "r2_permuted": r2_perm, + "r2_rotated90": r2_rot, "bar": "<= r2_single_geo + 0.03", + "verdict": "PASS" if b0_pass else "VOID"}, + "single_models": {"geo": {"c": c_geo_single, "R2": r2_geo}, + "bow": {"c": c_bow_single, "R2": r2_bow}}, + "joint_model": {"c_geo": c_geo, "c_bow": c_bow, "R2": r2_joint}, + "B1_identifiability": {"joint_R2": r2_joint, "best_single_R2": best_single, + "margin": r2_joint - best_single, + "bar": ">= best_single_R2 + 0.10", + "verdict": ("PASS" if b1_pass else "FAIL") + if b0_pass else "VOID (B0 failed)"}, + "B2_sign": {"c_geo": c_geo, "c_bow": c_bow, + "units_note": ("c_geo dimensionless (both D and P_geo " + "are Pa/km); c_bow has units km^-1 (D is " + "Pa/km, P_bow is Pa) -- see " + "fitted_contribution_Pa_per_km for a " + "unit-consistent magnitude comparison, " + "not raw |c_bow|"), + "verdict": ("PASS" if b2_pass else "FAIL") + if b0_pass else "VOID (B0 failed)"}, + "fitted_contribution_Pa_per_km": { + "mean_|D|": mean_D_mag, + "mean_|c_geo*P_geo|": mean_geo_contrib_mag, + "mean_|c_bow*P_bow|": mean_bow_contrib_mag, + "geo_contrib_frac_of_D": mean_geo_contrib_mag / mean_D_mag, + "bow_contrib_frac_of_D": mean_bow_contrib_mag / mean_D_mag}, + "B3_residual_resultant": { + "overall": {"R_bar": rbar_all, "mu_deg": mu_all, "rayleigh_p": p_all}, + "stranded_lt8ms": {"n": int(stranded.sum()), "R_bar": rbar_s, + "mu_deg": mu_s, "rayleigh_p": p_s}, + "moving_ge8ms": {"n": int(moving.sum()), "R_bar": rbar_m, + "mu_deg": mu_m, "rayleigh_p": p_m}}, + "B4_per_storm": per_storm, + } + with open(out_dir / "comet_tail_w6.json", "w") as fh: + json.dump(out, fh, indent=2) + with open(tag_path, "a") as tf: + tf.write(f"DONE B0={out['B0_controls']['verdict']} " + f"B1={out['B1_identifiability']['verdict']} " + f"B2={out['B2_sign']['verdict']}\n") + if partial_path.exists(): + partial_path.unlink() + return out + + +if __name__ == "__main__": + res = run() + print(f"n={res['n']}") + print("B0 controls:", res["B0_controls"]) + print("single models:", res["single_models"]) + print("joint model:", res["joint_model"]) + print("fitted contribution (Pa/km):", res["fitted_contribution_Pa_per_km"]) + print("B1:", res["B1_identifiability"]) + print("B2:", res["B2_sign"]) + print("B3 overall:", res["B3_residual_resultant"]["overall"]) + print("B3 stranded (<8 m/s):", res["B3_residual_resultant"]["stranded_lt8ms"]) + print("B3 moving (>=8 m/s):", res["B3_residual_resultant"]["moving_ge8ms"]) diff --git a/probes/weather-p1/exec-runs/comet_tail_w6.txt b/probes/weather-p1/exec-runs/comet_tail_w6.txt new file mode 100644 index 00000000..8c44daaf --- /dev/null +++ b/probes/weather-p1/exec-runs/comet_tail_w6.txt @@ -0,0 +1,6 @@ +START seed=20260812 n_storms=19 resumed=0 +progress 5/19 t0=62410 +progress 10/19 t0=66558 +progress 15/19 t0=70218 +progress 19/19 t0=74366 +DONE B0=VOID B1=VOID (B0 failed) B2=VOID (B0 failed)