From eb78265a7016c99c75f5ecc4a162973c448bf4ee Mon Sep 17 00:00:00 2001 From: LightChainr Date: Sat, 12 Sep 2026 14:54:26 +0800 Subject: [PATCH 1/3] Repair covariance support inference; recover N580 covariance from committed deletions --- .github/workflows/n580-recovery.yml | 25 ++ notes/projective-support-repair-20260912.md | 63 +++ .../projective-support.json | 19 + scripts/projective_inference.py | 361 +++++++++--------- scripts/recover_n580_covariance.py | 106 +++++ tests/test_n580_covariance_recovery.py | 22 ++ tests/test_projective_inference.py | 259 +++++-------- 7 files changed, 502 insertions(+), 353 deletions(-) create mode 100644 .github/workflows/n580-recovery.yml create mode 100644 notes/projective-support-repair-20260912.md create mode 100644 results/research-control-20260912/projective-support.json create mode 100644 scripts/recover_n580_covariance.py create mode 100644 tests/test_n580_covariance_recovery.py diff --git a/.github/workflows/n580-recovery.yml b/.github/workflows/n580-recovery.yml new file mode 100644 index 00000000..34d99df2 --- /dev/null +++ b/.github/workflows/n580-recovery.yml @@ -0,0 +1,25 @@ +name: N580 existing-block recovery +on: + pull_request: + branches: [main] + paths: + - 'scripts/recover_n580_covariance.py' + - '.github/workflows/n580-recovery.yml' +permissions: + contents: read +jobs: + recover: + runs-on: ubuntu-latest + timeout-minutes: 3 + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: '3.11' + - run: python -m pip install mpmath + - name: Targeted mathematical tests + run: | + python -m unittest discover -s tests -p 'test_projective_inference.py' + python -m unittest discover -s tests -p 'test_n580_covariance_recovery.py' + - name: Recover and report without rerunning the experiment + run: python scripts/recover_n580_covariance.py diff --git a/notes/projective-support-repair-20260912.md b/notes/projective-support-repair-20260912.md new file mode 100644 index 00000000..f3d967e4 --- /dev/null +++ b/notes/projective-support-repair-20260912.md @@ -0,0 +1,63 @@ +# P3: covariance support repair and N580 recovery without production + +Date: 2026-09-12. Companion to the exact-foundations correction (#702). + +## Mathematical correction + +For nonsingular S and a fixed full-column-rank signed-real design V, +D=min_a (y-Va)^T S^{-1}(y-Va) has the ordinary Gaussian-reference residual +interpretation. Two coordinates/one line give exactly Fieller's contrast squared. +Neither this algebra nor correct Fieller inference creates new information. +The three-versus-two-rung comparison changes the data used, not the two-rung test. + +For singular S whose nullspace is KNOWN as an exact deterministic constraint, +feasibility requires U0^T(y-Va)=0, where U0 spans ker S. Find a particular +solution a0 and a null basis Z for U0^T V. On the stochastic support whiten +with W=Lambda_+^{-1/2} U_+^T and fit W V Z against W(y-Va0). The correct df is +rank(S)-rank(W V Z), not blindly rank(S)-number_of_original_amplitudes. +If the support constraint is infeasible, there is no compatible model mean. + +Concrete regressions: S=diag(1,1,0), y=(0,0,1), v=(1,0,0) was D=0,p=1 in +blob `1b6270238db586c5c48849bd607376bf3c7d14fa`; it now raises an explicit +incompatibility in exact_support mode. For v=(1,0,1), the deterministic +coordinate forces a=1, giving D=1,df=2. A matrix with eigenvalue -0.1 is +rejected rather than clipped into an apparently valid covariance. + +Default `covariance_mode='strict'` refuses rank-deficient/truncated inputs. +`exact_support` is an explicit scientific assertion, NOT something inferred +from too few jackknife batches. Small positive variances cut off numerically +are not silently made deterministic. Exact-support computation is numerical +at the recorded precision/tolerance; it is not an exact rational certificate. +For empirical low rank, explicitly choose/justify a reduced observation space +or acquire adequate covariance information before using a reference p-value. + +A signed line {a*v:a in R} is not a positive ray. Cone/bounded-amplitude tests +need their own reference law. Estimated covariance likewise does not supply +exact finite-sample chi-square coverage. These distinctions also correct the +broader claims in the P3 manuscript; this patch does not promote its readiness. + +## N580: recover before replaying + +The main-branch files `results/aspect-ladder-n580/shards/rung_r{1,2,4}.json` +retain 100 aligned `_deleted` rows, plus seed, replica offset and sample count. +The missing cov(r2,r4) therefore need not be set to zero or commissioned as a +new run. `scripts/recover_n580_covariance.py` reads these arrays, checks the +#577 alignment metadata and previously saved diagonal/cross entries, then +reports the complete covariance, raw curvature and the existing eight ray tests. +It writes only stdout: historical results and freezes remain untouched. + +An additional labelled sensitivity fits a common-modulus H8 nuisance plane and +bounded common ratios |rho|<=0.2 and <=1. The bounded results are MINIMUM +RESIDUALS ONLY, not chi-square-calibrated tests. This does not identify H8, +validate a common ratio across moduli, or turn the reanalysis into new evidence. + +## Validation and limits + +Locally: 14 projective checks and 2 covariance-reader controls passed. The +full-rank implementation exactly reproduced all returned fields in 80 +fixed-seed comparisons against the saved original module. No production data +were changed. The real-shard end-to-end command is also put in a bounded, +read-only, path-triggered GitHub Actions job (three-minute cap), so it can run +where a complete repository checkout is available. Its actual result must be +read before any numerical N580 verdict is claimed. Full repository CI is a +separate check; it was not run in this local connector-only checkout. diff --git a/results/research-control-20260912/projective-support.json b/results/research-control-20260912/projective-support.json new file mode 100644 index 00000000..bfe073c5 --- /dev/null +++ b/results/research-control-20260912/projective-support.json @@ -0,0 +1,19 @@ +{ + "local_tests_passed": 14, + "full_rank_exact_output_matches": 80, + "historical_source_blob": "1b6270238db586c5c48849bd607376bf3c7d14fa", + "exact_support_constraint_example": { + "statistic": 1.0, + "degrees_of_freedom": 2, + "amplitudes": [1.0], + "covariance_rank": 2, + "covariance_condition_number": 1.0, + "p_value": 0.6065306597126334, + "equivalent_sigma": 0.515031998812217, + "covariance_mode": "exact_support", + "deterministic_constraint_rank": 1, + "fitted_stochastic_rank": 0, + "support_numerical_tolerance": "1.0e-25" + }, + "full_repository_suite": "not run locally" +} diff --git a/scripts/projective_inference.py b/scripts/projective_inference.py index 1b627023..ccca1e62 100644 --- a/scripts/projective_inference.py +++ b/scripts/projective_inference.py @@ -1,214 +1,209 @@ #!/usr/bin/env python3 -"""Test a measured response vector against a model's ray or subspace, without -ever forming a ratio. - -Why this module exists. The N=580 modulus ladder measured a three-entry -response `(A4(i), A4(2i), A4(4i))` and then scored competitors on the ratio -`A4(4i)/A4(i)`. `A4(i)` came in 3.6 sigma from zero, so the ratio's sampling -distribution was nothing like normal and its standard error at the observed -point understated the spread out where the large predictions sit. Every -verdict computed that way was wrong; three of eight flipped when the test was -moved onto the linear contrast `Y - R0 X` instead (Fieller). - -Fieller is the two-entry, one-ray case of a general fact: **a model that -predicts only proportions predicts a ray, and a model that predicts a -low-dimensional family predicts a subspace.** Neither is a statement about any -single coordinate, so neither should be tested by dividing one coordinate by -another. Dividing chooses a denominator, and the choice is arbitrary; when the -chosen denominator is poorly resolved the test degrades for a reason that has -nothing to do with the physics. - -What to compute instead. For a model ray `v` and a measured `y` with covariance -`S`, the statistic is the covariance-weighted distance from `y` to the line -through `v`: - - D(v) = min_a (y - a v)^T S^+ (y - a v), - -referred to chi-square with `rank(S) - 1` degrees of freedom. For a model whose -image is an `r`-dimensional subspace `V` the same expression minimises over all -of `V` and carries `rank(S) - r` degrees of freedom. There is no denominator, -so there is nothing to be badly conditioned except `S` itself -- and `S`'s -condition number is reported rather than hidden, because that *is* the honest -statement of how much the geometry can support. - -This also buys something the ratio test could not have. A systematic that -enters along a known direction -- the spin-8 leakage in the ladder's r=2 rung, -say -- can be carried as an extra basis vector of the model subspace instead of -by discarding the contaminated entry. A nuisance direction costs one degree of -freedom; discarding an entry costs the whole entry. - -Precision. The linear algebra runs in mpmath at a caller-chosen precision, not -in float64. The subject of the module is conditioning, and a routine that -silently lost three digits inverting an ill-conditioned covariance would be -answering a question about itself rather than about the data. -""" - -from __future__ import annotations +"""Covariance-weighted tests of fixed signed-real model subspaces. -from typing import Any, Mapping, Sequence +Full-rank inference retains the original GLS/Fieller calculation. Rank deficiency +is not automatically missing information or deterministic information: the +caller must say which. Default mode rejects singular/truncated covariance. +``covariance_mode='exact_support'`` imposes its nullspace as a known affine +support, profiles the remaining amplitudes, and computes the constrained df. +Empirical low rank from too few batches does NOT license exact_support. +These are Gaussian-reference/asymptotic tests, not exact finite-sample coverage +claims for an estimated covariance. A signed line is not a nonnegative ray. +""" +from __future__ import annotations +from typing import Any, Sequence from mpmath import mp - DEFAULT_DPS = 50 -# Eigenvalues below this fraction of the largest are treated as null directions -# of the covariance rather than as very small variances. A measured covariance -# with a genuine null direction is a covariance whose data cannot speak in that -# direction at all; inverting it there manufactures certainty. -DEFAULT_RELATIVE_TOLERANCE = "1e-12" +DEFAULT_RELATIVE_TOLERANCE = '1e-12' -def _as_matrix(rows: Sequence[Sequence[Any]]): - size = len(rows) - if any(len(row) != size for row in rows): - raise ValueError("covariance must be square") - return mp.matrix([[mp.mpf(str(value)) for value in row] for row in rows]) +def _as_matrix(rows): + n = len(rows) + if not n or any(len(row) != n for row in rows): + raise ValueError('covariance must be nonempty and square') + result = mp.matrix([[mp.mpf(str(x)) for x in row] for row in rows]) + if any(not mp.isfinite(x) for x in result): + raise ValueError('covariance entries must be finite') + return result -def _as_vector(values: Sequence[Any]): - return mp.matrix([mp.mpf(str(value)) for value in values]) +def _as_vector(values): + result = mp.matrix([mp.mpf(str(x)) for x in values]) + if any(not mp.isfinite(x) for x in result): + raise ValueError('vector entries must be finite') + return result -def spectral_pseudo_inverse(covariance: Sequence[Sequence[Any]], - relative_tolerance: str = DEFAULT_RELATIVE_TOLERANCE): - """Moore-Penrose inverse of a symmetric covariance, with its rank and condition. - - Returns ``(pinv, rank, condition_number, eigenvalues)``. Directions whose - eigenvalue falls below ``relative_tolerance`` times the largest are dropped: - a covariance estimated from a finite jackknife can be singular by - construction (100 batches cannot support more than 99 directions), and - inverting such a direction invents information. - """ +def _spectral(covariance, relative_tolerance): + tol = mp.mpf(relative_tolerance) + if not 0 < tol < 1: + raise ValueError('relative_tolerance must lie in (0,1)') matrix = _as_matrix(covariance) - size = matrix.rows - for i in range(size): - for j in range(i + 1, size): - if abs(matrix[i, j] - matrix[j, i]) > mp.mpf("1e-30") * ( - abs(matrix[i, j]) + abs(matrix[j, i]) + mp.mpf(1)): - raise ValueError(f"covariance is not symmetric at ({i}, {j})") - matrix[j, i] = matrix[i, j] + n = matrix.rows + for i in range(n): + for j in range(i+1,n): + if abs(matrix[i,j]-matrix[j,i]) > mp.mpf('1e-30')*( + abs(matrix[i,j])+abs(matrix[j,i])+1): + raise ValueError('covariance is not symmetric') + matrix[j,i] = matrix[i,j] values, vectors = mp.eigsy(matrix) - eigenvalues = [values[i] for i in range(size)] + eigenvalues = [values[i] for i in range(n)] + scale = max(abs(v) for v in eigenvalues) + # PSD validation is NOT the user-selected statistical rank cutoff. + roundoff = 100*mp.eps*scale + if min(eigenvalues) < -roundoff: + raise ValueError('covariance is not positive semidefinite') largest = max(eigenvalues) - if largest <= 0: - raise ValueError("covariance has no positive direction") - floor = mp.mpf(relative_tolerance) * largest - kept = [i for i in range(size) if eigenvalues[i] > floor] - if not kept: - raise ValueError("covariance has no direction above the tolerance") - pinv = mp.zeros(size, size) - for index in kept: - inverse = 1 / eigenvalues[index] - for i in range(size): - for j in range(size): - pinv[i, j] += inverse * vectors[i, index] * vectors[j, index] - smallest_kept = min(eigenvalues[i] for i in kept) - return pinv, len(kept), largest / smallest_kept, eigenvalues - - -def subspace_residual(observed: Sequence[Any], - covariance: Sequence[Sequence[Any]], + kept = [i for i,v in enumerate(eigenvalues) if v > tol*largest] + pinv = mp.zeros(n,n) + for k in kept: + inverse = 1/eigenvalues[k] + for i in range(n): + for j in range(n): + pinv[i,j] += inverse*vectors[i,k]*vectors[j,k] + condition = largest/min(eigenvalues[k] for k in kept) if kept else None + return pinv, kept, condition, eigenvalues, vectors, roundoff + + +def spectral_pseudo_inverse(covariance, relative_tolerance=DEFAULT_RELATIVE_TOLERANCE): + """Algebraic pseudoinverse and diagnostics; not by itself a test.""" + pinv,kept,condition,values,_,_ = _spectral(covariance,relative_tolerance) + return pinv,len(kept),condition,values + + +def _least_squares_and_kernel(matrix, target, tol): + """Minimum-norm fit, null basis and identifiable rank via SVD.""" + width = matrix.cols + if matrix.rows == 0: + return mp.zeros(width,1), mp.eye(width), 0 + if width == 0: + return mp.zeros(0,1), mp.zeros(0,0), 0 + U,s,Vt = mp.svd(matrix,full_matrices=True) + largest = max(s) if len(s) else mp.mpf(0) + retained = [i for i in range(len(s)) if s[i] > tol*largest] + solution = mp.zeros(width,1) + for i in retained: + coefficient = mp.fsum(U[j,i]*target[j] for j in range(matrix.rows))/s[i] + for j in range(width): + solution[j] += Vt[i,j]*coefficient + free = [i for i in range(width) if i not in retained] + null = mp.matrix(width,len(free)) + for col,i in enumerate(free): + for j in range(width): + null[j,col] = Vt[i,j] + return solution,null,len(retained) + + +def _exact_support_fit(y,directions,kept,values,vectors,tol,roundoff): + n,width = y.rows,len(directions) + dropped = [i for i in range(n) if i not in kept] + if any(abs(values[i]) > roundoff for i in dropped): + raise ValueError('a numerical cutoff is not an exact support: lower relative_tolerance') + V = mp.matrix([[directions[j][i] for j in range(width)] for i in range(n)]) + U0 = mp.matrix([[vectors[i,j] for i in range(n)] for j in dropped]) + constraints, rhs = U0*V,U0*y + a0,Z,constraint_rank = _least_squares_and_kernel(constraints,rhs,tol) + support_tol = mp.power(10,-mp.dps//2) + if mp.norm(constraints*a0-rhs) > support_tol*max(1,mp.norm(y),mp.norm(V*a0)): + raise ValueError('model is incompatible with exact covariance support') + W = mp.matrix(len(kept),n) + for row,k in enumerate(kept): + for i in range(n): + W[row,i] = vectors[i,k]/mp.sqrt(values[k]) + if Z.cols and kept: + fit,_,fitted_rank = _least_squares_and_kernel(W*V*Z,W*(y-V*a0),tol) + amplitudes = a0+Z*fit + else: + fitted_rank, amplitudes = 0,a0 + residual = W*(y-V*amplitudes) + statistic = mp.fsum(x*x for x in residual) + degrees = len(kept)-fitted_rank + if degrees == 0 and statistic < support_tol**2: + statistic = mp.mpf(0) + return statistic,degrees,amplitudes,constraint_rank,fitted_rank,support_tol + + +def subspace_residual(observed: Sequence[Any], covariance: Sequence[Sequence[Any]], basis: Sequence[Sequence[Any]], - relative_tolerance: str = DEFAULT_RELATIVE_TOLERANCE, - dps: int = DEFAULT_DPS) -> dict[str, Any]: - """Covariance-weighted distance from ``observed`` to the span of ``basis``. - - ``basis`` is a sequence of direction vectors. One direction is a ray -- a - model that predicts proportions only. Two or more span the image of a model - with that many free amplitudes, or a model plus a known nuisance direction. - """ + relative_tolerance=DEFAULT_RELATIVE_TOLERANCE, + dps=DEFAULT_DPS, *, covariance_mode='strict') -> dict[str,Any]: + if covariance_mode not in ('strict','exact_support'): + raise ValueError('covariance_mode must be strict or exact_support') + if dps < 30: + raise ValueError('use at least 30 decimal digits') with mp.workdps(dps): y = _as_vector(observed) - size = y.rows - directions = [_as_vector(row) for row in basis] - if not directions: - raise ValueError("need at least one model direction") - for direction in directions: - if direction.rows != size: - raise ValueError("model direction has the wrong length") - pinv, rank, condition, _ = spectral_pseudo_inverse(covariance, relative_tolerance) - if pinv.rows != size: - raise ValueError("covariance and observation disagree on dimension") - - # Normal equations for min_a (y - Va)^T S^+ (y - Va): (V^T S^+ V) a = V^T S^+ y. - width = len(directions) - gram = mp.zeros(width, width) - rhs = mp.zeros(width, 1) - weighted = [pinv * direction for direction in directions] - for i in range(width): - for j in range(width): - gram[i, j] = (directions[i].T * weighted[j])[0] - rhs[i] = (directions[i].T * (pinv * y))[0] - model_rank = _rank_of(gram, relative_tolerance) - if model_rank < width: - raise ValueError( - "model directions are linearly dependent after weighting; drop one" - ) - amplitudes = mp.lu_solve(gram, rhs) - residual = y - sum((amplitudes[i] * directions[i] for i in range(width)), - mp.zeros(size, 1)) - statistic = (residual.T * (pinv * residual))[0] - degrees = rank - width - if degrees < 0: - raise ValueError("model has more directions than the covariance has rank") - return { - "statistic": float(statistic), - "degrees_of_freedom": degrees, - "amplitudes": [float(amplitudes[i]) for i in range(width)], - "covariance_rank": rank, - "covariance_condition_number": float(condition), - "p_value": chi_square_upper_tail(float(statistic), degrees), - "equivalent_sigma": _equivalent_sigma(float(statistic), degrees), - } - - -def ray_residual(observed: Sequence[Any], - covariance: Sequence[Sequence[Any]], - direction: Sequence[Any], - relative_tolerance: str = DEFAULT_RELATIVE_TOLERANCE, - dps: int = DEFAULT_DPS) -> dict[str, Any]: - """The one-direction case: a model that predicts proportions and nothing else. - - For a two-entry observation this is exactly Fieller's test on the ratio of - the two entries, squared -- see the tests. For three or more entries it is - the thing Fieller cannot do: it uses every entry at once, and never has to - nominate one of them as a denominator. - """ - return subspace_residual(observed, covariance, [direction], relative_tolerance, dps) - - -def _rank_of(matrix, relative_tolerance: str) -> int: - values, _ = mp.eigsy(matrix) + directions = [_as_vector(v) for v in basis] + if not directions or any(v.rows != y.rows for v in directions): + raise ValueError('need model directions of the observed dimension') + pinv,kept,condition,values,vectors,roundoff = _spectral(covariance,relative_tolerance) + if pinv.rows != y.rows: + raise ValueError('covariance and observation disagree on dimension') + rank,width = len(kept),len(directions) + extra = {} + if rank < y.rows: + if covariance_mode == 'strict': + raise ValueError('singular or truncated covariance: exact support must be declared; ' + 'empirical low rank does not establish deterministic constraints') + statistic,degrees,amplitudes,cr,mr,st = _exact_support_fit( + y,directions,kept,values,vectors,mp.mpf(relative_tolerance),roundoff) + extra = {'covariance_mode':'exact_support','deterministic_constraint_rank':cr, + 'fitted_stochastic_rank':mr,'support_numerical_tolerance':str(st)} + else: + # Preserve the original full-rank arithmetic and return fields. + gram = mp.zeros(width,width) + rhs = mp.zeros(width,1) + weighted = [pinv*v for v in directions] + for i in range(width): + for j in range(width): + gram[i,j] = (directions[i].T*weighted[j])[0] + rhs[i] = (directions[i].T*(pinv*y))[0] + if _rank_of(gram,relative_tolerance) < width: + raise ValueError('model directions are linearly dependent after weighting; drop one') + amplitudes = mp.lu_solve(gram,rhs) + residual = y-sum((amplitudes[i]*directions[i] for i in range(width)),mp.zeros(y.rows,1)) + statistic = (residual.T*(pinv*residual))[0] + degrees = rank-width + if degrees < 0: + raise ValueError('model has more directions than covariance rank') + return {'statistic':float(statistic),'degrees_of_freedom':degrees, + 'amplitudes':[float(x) for x in amplitudes], 'covariance_rank':rank, + 'covariance_condition_number':float(condition) if condition is not None else None, + 'p_value':chi_square_upper_tail(float(statistic),degrees), + 'equivalent_sigma':_equivalent_sigma(float(statistic),degrees), **extra} + + +def ray_residual(observed,covariance,direction, + relative_tolerance=DEFAULT_RELATIVE_TOLERANCE,dps=DEFAULT_DPS, + *,covariance_mode='strict'): + """Signed-real line test; in two nonsingular dimensions this is Fieller z^2.""" + return subspace_residual(observed,covariance,[direction],relative_tolerance,dps, + covariance_mode=covariance_mode) + + +def _rank_of(matrix,relative_tolerance): + values,_ = mp.eigsy(matrix) magnitudes = [abs(values[i]) for i in range(matrix.rows)] largest = max(magnitudes) - if largest == 0: - return 0 - floor = mp.mpf(relative_tolerance) * largest - return sum(1 for value in magnitudes if value > floor) + return sum(x > mp.mpf(relative_tolerance)*largest for x in magnitudes) if largest else 0 -def chi_square_upper_tail(statistic: float, degrees: int) -> float: - """P(chi^2_degrees >= statistic), via the regularised incomplete gamma.""" +def chi_square_upper_tail(statistic,degrees): if degrees <= 0: return 1.0 if statistic <= 0 else 0.0 if statistic <= 0: return 1.0 with mp.workdps(40): - return float(mp.gammainc(mp.mpf(degrees) / 2, mp.mpf(statistic) / 2, - mp.inf, regularized=True)) - + return float(mp.gammainc(mp.mpf(degrees)/2,mp.mpf(statistic)/2,mp.inf,regularized=True)) -def _equivalent_sigma(statistic: float, degrees: int) -> float: - """The two-sided normal deviate with the same tail probability. - Reported so that a many-degree-of-freedom result can be compared with the - one-degree-of-freedom numbers the project already quotes, without anyone - having to take a square root that is only correct when ``degrees == 1``. - """ - tail = chi_square_upper_tail(statistic, degrees) - if tail <= 0.0: - return float("inf") - if tail >= 1.0: +def _equivalent_sigma(statistic,degrees): + tail = chi_square_upper_tail(statistic,degrees) + if tail <= 0: + return float('inf') + if tail >= 1: return 0.0 with mp.workdps(40): - return float(-mp.sqrt(2) * mp.erfinv(mp.mpf(tail) - 1)) + return float(-mp.sqrt(2)*mp.erfinv(mp.mpf(tail)-1)) diff --git a/scripts/recover_n580_covariance.py b/scripts/recover_n580_covariance.py new file mode 100644 index 00000000..f8666096 --- /dev/null +++ b/scripts/recover_n580_covariance.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python3 +"""Reanalysis of the existing N580 shards. No simulation or data overwrite. + +Uses the aligned delete-one ordering recorded by the #577 replay, checking +its seed/offset/batch metadata and previously published covariance entries. +Outputs a new result to stdout. All model scores are conditional on the +stated angular/readout contract; this is not new independent evidence. +""" +from __future__ import annotations +import argparse +import json +import math +from pathlib import Path +from projective_inference import ray_residual,subspace_residual + +CHANNEL='P4_S_prime' +RUNGS=(1,2,4) + + +def jackknife_covariance(columns): + n=len(columns[0]) + if n<2 or any(len(c)!=n for c in columns): + raise ValueError('aligned columns with at least two deletions required') + if any(not math.isfinite(x) for c in columns for x in c): + raise ValueError('delete-one estimates must be finite') + means=[math.fsum(c)/n for c in columns] + return [[(n-1)/n*math.fsum((columns[i][k]-means[i])*(columns[j][k]-means[j]) + for k in range(n)) for j in range(len(columns))] + for i in range(len(columns))] + + +def recover(root): + directory=root/'results'/'aspect-ladder-n580' + shards=[json.loads((directory/'shards'/f'rung_r{r}.json').read_text()) for r in RUNGS] + for key in ('seed','replica_offset','samples','batches'): + if len({str(s[key]) for s in shards})!=1: + raise ValueError(f'unaligned replay metadata: {key}') + n=shards[0]['batches'] + columns=[[row[CHANNEL] for row in s['measured']['_deleted']] for s in shards] + if any(len(c)!=n for c in columns): + raise ValueError('incomplete delete-one array') + y=[s['measured']['channels'][CHANNEL]['value'] for s in shards] + S=jackknife_covariance(columns) + for i,s in enumerate(shards): + expected=s['measured']['channels'][CHANNEL]['standard_error']**2 + if not math.isclose(S[i][i],expected,rel_tol=1e-8,abs_tol=1e-20): + raise ValueError('jackknife diagonal does not reproduce saved standard error') + old=json.loads((directory/'latest.json').read_text()) + for entry,j in (('r2_over_r1',1),('r4_over_r1',2)): + item=old['ratios'][entry][CHANNEL] + if not math.isclose(S[0][j],item['covariance'],rel_tol=1e-7,abs_tol=1e-20): + raise ValueError('paired covariance does not reproduce #577') + reference=json.loads((root/'results'/'aspect-ladder-n580-projective'/'latest.json').read_text()) + if any(not math.isclose(a,b,rel_tol=1e-12,abs_tol=1e-15) + for a,b in zip(y,reference['response_vector'])): + raise ValueError('response vector changed') + scores={} + leak=1148/21025 + for name,item in reference['competitors'].items(): + v=item['ray'] + score=ray_residual(y,S,v) + score['excluded_at_nominal_3_sigma']=score['p_value'] None: - """The anchor. Fieller is the 2x1 case, so this must reproduce it. - - The wrong number this stops us believing is a generalisation that has - drifted into being a different test. Fieller's z on the ratio of two - entries is already established; if the covariance-weighted distance to - the model ray is not its square, one of the two is wrong, and this test - does not say which -- it says stop. - """ - covariance = self._two_entry_covariance() - observed = [self.Y[0], self.Y[2]] - for predicted in (1.0, 4.0, 10.9908008589, 16.0, 120.79770352, 2080.30719731): - fieller = ( - (self.Y[2] - predicted * self.Y[0]) - / math.sqrt(self.VAR[2] + predicted ** 2 * self.VAR[0] - - 2 * predicted * covariance[0][1]) - ) - result = projective.ray_residual(observed, covariance, [1.0, predicted]) - self.assertEqual(result["degrees_of_freedom"], 1) - self.assertAlmostEqual( - result["statistic"], fieller ** 2, - delta=abs(fieller ** 2) * 1e-9 + 1e-12, - ) - - def test_no_denominator_is_ever_chosen(self) -> None: - """Scaling the model direction must not move the statistic. - - The wrong number here is one that depends on which entry was written - first. A ray is a ray: (1, R) and (c, cR) are the same model, and a - test that distinguished them would be testing our bookkeeping. - """ - covariance = self._two_entry_covariance() - observed = [self.Y[0], self.Y[2]] - base = projective.ray_residual(observed, covariance, [1.0, 4.0])["statistic"] - for scale in (1e-6, 0.5, 7.0, 1e6, -3.0): - scaled = projective.ray_residual( - observed, covariance, [scale, 4.0 * scale] - )["statistic"] - self.assertAlmostEqual(scaled, base, delta=abs(base) * 1e-9 + 1e-12) - - def test_a_model_through_the_data_has_no_residual(self) -> None: - """A ray that passes through the observation must score zero. - - The wrong number this catches is a statistic that has picked up an - offset -- from a mis-signed cross term, say -- which would inflate every - model equally and be invisible in a comparison between models. - """ - covariance = self._two_entry_covariance() - observed = [self.Y[0], self.Y[2]] - exact = projective.ray_residual( - observed, covariance, [self.Y[0], self.Y[2]] - ) - self.assertLess(exact["statistic"], 1e-20) - self.assertGreater(exact["p_value"], 0.999) - - def test_a_nuisance_direction_costs_one_degree_of_freedom(self) -> None: - """Carrying a known systematic must be cheaper than discarding an entry. - - The frozen ladder dropped the r=2 rung from its decision because that - rung carries spin-8 leakage. The wrong belief this test guards is that - dropping was the only option: a systematic with a known direction is a - second basis vector, which costs one degree of freedom rather than a - whole measurement. - """ - s = [math.sqrt(v) for v in self.VAR] - c14 = self.RHO_14 * s[0] * s[2] - c12 = 0.0209 * s[0] * s[1] - covariance = [[self.VAR[0], c12, c14], - [c12, self.VAR[1], 0.0], - [c14, 0.0, self.VAR[2]]] - ray = projective.subspace_residual(self.Y, covariance, [[1.0, 2.0, 4.0]]) - with_nuisance = projective.subspace_residual( - self.Y, covariance, [[1.0, 2.0, 4.0], [-1.0, 1.0, -1.0]] - ) - self.assertEqual(ray["degrees_of_freedom"], 2) - self.assertEqual(with_nuisance["degrees_of_freedom"], 1) - # Absorbing a direction can only reduce the distance to the model set. - self.assertLessEqual(with_nuisance["statistic"], ray["statistic"] + 1e-15) - - def test_dependent_model_directions_are_refused_rather_than_fitted(self) -> None: - """Stops a model set that is secretly smaller than it looks. - - Two directions that span one line give a singular normal-equation - matrix. Solving it anyway returns an arbitrary split of one amplitude - between two coefficients, and a degrees-of-freedom count that is one too - low -- which makes every p-value too small. - """ - covariance = self._two_entry_covariance() - with self.assertRaisesRegex(ValueError, "linearly dependent"): - projective.subspace_residual( - [self.Y[0], self.Y[2]], covariance, [[1.0, 4.0], [2.0, 8.0]] - ) - - def test_a_null_direction_of_the_covariance_is_dropped_not_inverted(self) -> None: - """Stops manufactured certainty along a direction with no information. - - A jackknife over 100 batches cannot support more than 99 directions, so - a measured covariance can be singular by construction. Inverting a zero - eigenvalue would make the residual in that direction infinitely - significant. The rank must fall instead, and the degrees of freedom with - it. - """ - covariance = [[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 1.0]] - result = projective.ray_residual([1.0, 1.0, 2.0], covariance, [1.0, 1.0, 1.0]) - self.assertEqual(result["covariance_rank"], 2) - self.assertEqual(result["degrees_of_freedom"], 1) - self.assertTrue(math.isfinite(result["statistic"])) - - def test_the_chi_square_tail_matches_values_that_can_be_checked_by_hand(self) -> None: - """The p-values are computed, not tabulated, so check them once. - - At one degree of freedom the tail is the two-sided normal tail, which is - known in closed form; at two it is a plain exponential. - """ - for sigma in (1.0, 2.0, 3.0, 5.0): - expected = math.erfc(sigma / math.sqrt(2.0)) - self.assertAlmostEqual( - projective.chi_square_upper_tail(sigma ** 2, 1) / expected, 1.0, places=9 - ) - for statistic in (0.5, 2.0, 9.0): - self.assertAlmostEqual( - projective.chi_square_upper_tail(statistic, 2), - math.exp(-statistic / 2.0), places=12, - ) - - def test_equivalent_sigma_is_a_tail_not_a_square_root(self) -> None: - """Stops sqrt(chi-square) being quoted as sigma at more than one df. - - The wrong number here is the one everybody writes: sqrt(D) reported as - a sigma when D has two or more degrees of freedom, which overstates the - significance because a larger D is expected under the null. - """ - self.assertAlmostEqual(projective._equivalent_sigma(9.0, 1), 3.0, places=9) - two = projective._equivalent_sigma(9.0, 2) - self.assertLess(two, 3.0) - self.assertGreater(two, 2.0) - - -if __name__ == "__main__": + Y=[9.01643304e-4,2.91097703e-3,4.13180763e-3] + VAR=[6.207021e-8,6.713835e-8,3.807515e-8] + S=[[1,0,0],[0,1,0],[0,0,0]] + + def pair(self): + c=-.1648*math.sqrt(self.VAR[0]*self.VAR[2]) + return [self.Y[0],self.Y[2]],[[self.VAR[0],c],[c,self.VAR[2]]] + + def test_fieller_identity(self): + y,s=self.pair() + for r in (1,4,10.9908008589,16,120.79770352,2080.30719731): + z=(y[1]-r*y[0])/math.sqrt(s[1][1]-2*r*s[0][1]+r*r*s[0][0]) + got=p.ray_residual(y,s,[1,r]) + self.assertEqual(got['degrees_of_freedom'],1) + self.assertAlmostEqual(got['statistic'],z*z,delta=1e-9*z*z+1e-12) + + def test_signed_line_scaling(self): + y,s=self.pair() + reference=p.ray_residual(y,s,[1,4])['statistic'] + for c in (1e-6,.5,7,1e6,-3): + self.assertAlmostEqual(p.ray_residual(y,s,[c,4*c])['statistic'],reference, + delta=1e-9*reference+1e-12) + + def test_exact_fit(self): + y,s=self.pair() + self.assertLess(p.ray_residual(y,s,y)['statistic'],1e-20) + + def test_nuisance_lowers_df_and_distance(self): + y=[1,3,4];s=[[1,0,0],[0,1,0],[0,0,1]] + one=p.subspace_residual(y,s,[[1,2,4]]) + two=p.subspace_residual(y,s,[[1,2,4],[-1,1,-1]]) + self.assertEqual((one['degrees_of_freedom'],two['degrees_of_freedom']),(2,1)) + self.assertLessEqual(two['statistic'],one['statistic']) + + def test_dependent_full_rank_design_refused(self): + y,s=self.pair() + with self.assertRaisesRegex(ValueError,'linearly dependent'): + p.subspace_residual(y,s,[[1,4],[2,8]]) + + def test_singular_covariance_requires_semantics(self): + with self.assertRaisesRegex(ValueError,'must be declared'): + p.ray_residual([0,0,1],self.S,[1,0,0]) + + def test_exact_support_contradiction_is_not_p_one(self): + with self.assertRaisesRegex(ValueError,'incompatible'): + p.ray_residual([0,0,1],self.S,[1,0,0],covariance_mode='exact_support') + + def test_exact_support_constrains_amplitude_and_df(self): + got=p.ray_residual([0,0,1],self.S,[1,0,1],covariance_mode='exact_support') + self.assertEqual(got['amplitudes'],[1]) + self.assertEqual(got['statistic'],1) + self.assertEqual(got['degrees_of_freedom'],2) + + def test_known_duplicate_coordinate_support(self): + s=[[1,1,0],[1,1,0],[0,0,1]] + got=p.ray_residual([1,1,2],s,[1,1,1],covariance_mode='exact_support') + self.assertEqual(got['covariance_rank'],2) + self.assertEqual(got['degrees_of_freedom'],1) + self.assertAlmostEqual(got['statistic'],.5) + + def test_negative_covariance_refused(self): + with self.assertRaisesRegex(ValueError,'positive semidefinite'): + p.ray_residual([0,0,1],[[1,0,0],[0,1,0],[0,0,-.1]],[1,0,0]) + + def test_cutoff_is_not_structural_null(self): + with self.assertRaisesRegex(ValueError,'numerical cutoff'): + p.ray_residual([0,0,1],[[1,0,0],[0,1,0],[0,0,1e-15]], + [1,0,1],covariance_mode='exact_support') + + def test_all_deterministic_compatible(self): + got=p.ray_residual([2,4],[[0,0],[0,0]],[1,2],covariance_mode='exact_support') + self.assertEqual(got['degrees_of_freedom'],0) + self.assertEqual(got['statistic'],0) + + def test_chi_square_tail(self): + for z in (1,2,3,5): + self.assertAlmostEqual(p.chi_square_upper_tail(z*z,1)/math.erfc(z/math.sqrt(2)),1,places=9) + for x in (.5,2,9): + self.assertAlmostEqual(p.chi_square_upper_tail(x,2),math.exp(-x/2),places=12) + + def test_sigma_is_not_sqrt_for_multiple_df(self): + self.assertAlmostEqual(p._equivalent_sigma(9,1),3,places=9) + self.assertTrue(2 Date: Sat, 12 Sep 2026 15:03:10 +0800 Subject: [PATCH 2/3] Record executed N580 covariance recovery and resolve the 3-sigma bare-aspect decision --- notes/n580-complete-covariance-20260912.md | 90 +++++++++++++++++++ .../n580-complete-covariance-summary.json | 48 ++++++++++ 2 files changed, 138 insertions(+) create mode 100644 notes/n580-complete-covariance-20260912.md create mode 100644 results/research-control-20260912/n580-complete-covariance-summary.json diff --git a/notes/n580-complete-covariance-20260912.md b/notes/n580-complete-covariance-20260912.md new file mode 100644 index 00000000..ab58c5fe --- /dev/null +++ b/notes/n580-complete-covariance-20260912.md @@ -0,0 +1,90 @@ +# N580: the missing covariance recovered, without another experiment + +Date: 2026-09-12. Analysis implemented in PR #703. + +## Executed, not proposed + +GitHub Actions run 34679359201, job 103514945963, completed successfully. +Command: `python scripts/recover_n580_covariance.py`. +Head: `eb78265a7016c99c75f5ecc4a162973c448bf4ee`; +actual tested merge checkout: `dda2a8e2dcb9e5a4a9215cb1f1877a4aa016b189`. +Runner: Python 3.11.16, mpmath 1.4.1. The 14 inference checks and 2 recovery +checks passed in the same job. Full-repository CI is separate from this job. + +The reader consumed the EXISTING 100 aligned delete-one rows in +`results/aspect-ladder-n580/shards/rung_r{1,2,4}.json`. It checked replay +seed/offset/sample/batch metadata and recovered the previously stored standard +errors and cov(r1,r2), cov(r1,r4). No Monte Carlo or histogram reconstruction +was run. All historical results and freezes are unchanged. + +The result summary is in +`results/research-control-20260912/n580-complete-covariance-summary.json`. +Fields there were extracted from the successful job stdout, then all eight +pure-ray statistics were independently recomputed locally from the recovered +vector/matrix; the largest D difference was zero. That local calculation was +NOT a second shard replay or independent experiment. + +## Recovered matrix and decision + +The missing cov(r2,r4) is 5.5621485097801135e-9, correlation +0.11001108422637242. The full matrix has rank 3 and condition number 1.96789. +It is not singular; the exact-support bug fixed in #703 does not invalidate +this well-conditioned calculation. + +| Pure model, retrospective three-rung test | D, df=2 | nominal Gaussian-reference p | equivalent sigma | declared 3-sigma decision | +|---|---:|---:|---:|---| +| bare aspect ratio | 10.8644599 | 0.00437333 | 2.84990 | not rejected | +| no modulus dependence | 89.8585102 | 3.07236e-20 | 9.21640 | rejected | +| area scaling | 57.5124933 | 3.24579e-13 | 7.28379 | rejected | +| Q4 weight-4 shape | 55.5877105 | 8.49732e-13 | 7.15288 | rejected | +| weight-12 E12 | 132.267324 | 1.89896e-29 | 11.26749 | rejected | +| weight-12 E4 cubed | 132.226117 | 1.93849e-29 | 11.26568 | rejected | +| weight-12 delta | 513.659218 | 2.88614e-112 | not resolved by inverse-erf conversion | rejected | +| weight-8 E8 | 115.804378 | 7.13508e-26 | 10.51803 | rejected | + +The previous missing-covariance ambiguity for bare_aspect_ratio is resolved: +it is **not rejected at the declared nominal 3-sigma cutoff**. This does not +confirm that law, make it unique physically, or license a post-hoc 2-sigma +cutoff. The original frozen two-rung comparison remains a different test with +its recorded underpowered verdict; this is the retrospective three-rung test +on the same block, conditional on its pure-amplitude/readout assumptions. + +The raw divided difference is + + (m4 - 3*m2 + 2*m1)/6 = -0.00046630613986373587 + SE = 0.0001499079748299409 + nominal z = -3.110615965513005. + +This 1-df contrast and the 2-df line test need not cross the same cutoff. +The sign alone cannot reject an arbitrary signed multiple of a convex shape; +any class-level sign claim must additionally declare the amplitude sign. +No continuum operator or percolation threshold is identified. + +## What the nuisance sensitivity says, and does not say + +Allowing a second column lambda*(-v1,+v2,-v4), lambda=1148/21025, assumes +that spin-8 has the SAME modulus shape and a common ratio to spin-4. Under +that enlarged, unbounded plane the nominal residuals are bare 0.25187/1, +area 6.54806/1, Q4 4.32890/1. The small-ratio bounds matter: at |rho|<=1, +minimum D is 5.72524, 50.81440, 48.88288 respectively. These bounded-cone +minima are not assigned an ordinary chi-square calibration. + +This is a sensitivity to a changed model, not a rescue of a rejected pure +ray and not a measurement identifying H8. A bound on |rho| across rungs does +not imply a COMMON rho; the latter is a separate stronger assumption. In +particular the older `required_spin8_ratio` wording that called common rho +an implication of a common bound is not logically valid. + +## Decision for the queue and manuscript + +Cancel any proposed N580 replay whose sole purpose is this covariance entry. +Use the recovered matrix and clearly separated pure-line/nuisance readings in +P3. Replace 'undetermined because covariance missing' with the narrow decision +above. Keep prospective and retrospective evidence separate. The manuscript +still needs the support/df, amplitude-sign, shared-ratio and prior-art language +corrected; a numerical recovery is not publication readiness. + +The scientific next step is not an unbudgeted larger amplitude ladder. First +resolve the angular/Smith-class identifiability in #589 and the within-model +shape question in #622. All p-values here are Gaussian-reference diagnostics +using estimated jackknife covariance, not exact finite-sample coverage. diff --git a/results/research-control-20260912/n580-complete-covariance-summary.json b/results/research-control-20260912/n580-complete-covariance-summary.json new file mode 100644 index 00000000..2bcea397 --- /dev/null +++ b/results/research-control-20260912/n580-complete-covariance-summary.json @@ -0,0 +1,48 @@ +{ + "schema": "matching-one.n580-complete-covariance-summary.20260912", + "provenance": { + "origin": "Selected fields extracted from successful GitHub Actions stdout; eight D scores independently recomputed locally from the recovered vector and covariance, not from a local shard checkout.", + "run_id": 34679359201, + "job_id": 103514945963, + "head_sha": "eb78265a7016c99c75f5ecc4a162973c448bf4ee", + "checkout_sha": "dda2a8e2dcb9e5a4a9215cb1f1877a4aa016b189", + "command": "python scripts/recover_n580_covariance.py", + "runner_python": "3.11.16", + "runner_mpmath": "1.4.1", + "local_max_statistic_difference": 0.0 + }, + "standing": "Retrospective three-rung reanalysis of the same existing block. Original frozen two-rung verdict unchanged.", + "batches": 100, + "rungs": [1, 2, 4], + "response_vector": [0.0009016433036753112, 0.0029109770266212776, 0.004131807633330795], + "covariance": [ + [6.207020770273794e-8, 1.3480872826970252e-9, -8.011196220988015e-9], + [1.3480872826970252e-9, 6.713835319545046e-8, 5.5621485097801135e-9], + [-8.011196220988015e-9, 5.5621485097801135e-9, 3.807514679910218e-8] + ], + "covariance_rank": 3, + "covariance_condition_number": 1.9678931965099733, + "correlation_r2_r4": 0.11001108422637242, + "curvature": { + "value": -0.00046630613986373587, + "standard_error": 0.0001499079748299409, + "z_nominal": -3.110615965513005 + }, + "competitors": { + "bare_aspect_ratio": {"ray": [1,2,4], "statistic": 10.86445991488253, "degrees_of_freedom": 2, "p_value_nominal": 0.004373332575404287, "equivalent_sigma_diagnostic": 2.849897460536986, "excluded_at_nominal_3_sigma": false}, + "no_modulus_dependence": {"ray": [1,1,1], "statistic": 89.85851023621422, "degrees_of_freedom": 2, "p_value_nominal": 3.072362291533769e-20, "equivalent_sigma_diagnostic": 9.216402139538845, "excluded_at_nominal_3_sigma": true}, + "plain_area_scaling": {"ray": [1,4,16], "statistic": 57.512493327374884, "degrees_of_freedom": 2, "p_value_nominal": 3.245792509791233e-13, "equivalent_sigma_diagnostic": 7.2837862071923, "excluded_at_nominal_3_sigma": true}, + "q4_jordan_weight4": {"ray": [1,2.75,10.9908008589], "statistic": 55.58771053243131, "degrees_of_freedom": 2, "p_value_nominal": 8.497321045543934e-13, "equivalent_sigma_diagnostic": 7.1528843204156916, "excluded_at_nominal_3_sigma": true}, + "weight12_E12": {"ray": [1,32.515625,2080.30719731], "statistic": 132.26732359383956, "degrees_of_freedom": 2, "p_value_nominal": 1.898959111319179e-29, "equivalent_sigma_diagnostic": 11.267494938720372, "excluded_at_nominal_3_sigma": true}, + "weight12_E4_cubed": {"ray": [1,20.796875,1327.6635036], "statistic": 132.22611687853123, "degrees_of_freedom": 2, "p_value_nominal": 1.938489879957022e-29, "equivalent_sigma_diagnostic": 11.265680304319943, "excluded_at_nominal_3_sigma": true}, + "weight12_delta": {"ray": [1,0.125,2.79010739704e-5], "statistic": 513.6592181764852, "degrees_of_freedom": 2, "p_value_nominal": 2.8861439189966654e-112, "equivalent_sigma_diagnostic": null, "excluded_at_nominal_3_sigma": true}, + "weight8_E8": {"ray": [1,7.5625,120.79770352], "statistic": 115.80437790454678, "degrees_of_freedom": 2, "p_value_nominal": 7.135080016926461e-26, "equivalent_sigma_diagnostic": 10.518028644032594, "excluded_at_nominal_3_sigma": true} + }, + "boundaries": [ + "Not rejected at the declared nominal 3-sigma threshold does not establish a law.", + "Pure-line decisions are conditional on the declared amplitude/readout contract.", + "A common H8/H4 ratio is stronger than a uniform bound on rung-specific ratios.", + "Estimated jackknife covariance; Gaussian-reference p-values, not exact finite-sample coverage.", + "Null extreme-tail sigma is an unresolved numerical conversion, not infinite evidence." + ] +} From 909d4b3f97f5fdc889630e869f7836a686d84037 Mon Sep 17 00:00:00 2001 From: LightChainr Date: Sat, 12 Sep 2026 15:34:02 +0800 Subject: [PATCH 3/3] Finish P3 saturated-design support repair and install existing numerical-test dependencies Prevent df=0 roundoff from producing a false rejection; preserve positive-df historical scores. Add six mathematical regressions. CI failure 103516017657 was missing NumPy, so install NumPy/SciPy for the existing complete suite without skipping tests. --- .github/workflows/ci.yml | 22 ++------ .github/workflows/n580-recovery.yml | 5 +- requirements-test.txt | 5 ++ scripts/projective_inference.py | 13 ++++- tests/test_projective_saturated_support.py | 58 ++++++++++++++++++++++ 5 files changed, 83 insertions(+), 20 deletions(-) create mode 100644 requirements-test.txt create mode 100644 tests/test_projective_saturated_support.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 48eed7a9..acce126d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,16 +1,8 @@ name: ci -# What CI is for here: catching the breakage you would not notice yourself. -# It is not a gate, not a claim, and not evidence. See GOVERNANCE.md sections 0 and 2. -# -# The tests are deterministic standard-library mathematics, so running the whole -# suite on three interpreters bought almost nothing and cost three times the budget. -# It now runs once, on 3.11; 3.9 and 3.13 are compile-checked, which is what actually -# catches a version incompatibility in this codebase. -# -# If branch protection is ever turned on for main, the checks to require are the ones -# this file defines: "Compile on 3.9", "Compile on 3.13", "Tests", and -# "C++17 build and self-tests". That is a note, not a machine-checked contract. +# CI catches computational breakage, not scientific truth. +# Run the complete suite once on 3.11; compile-check 3.9 and 3.13. +# The full suite includes the P598 NumPy/SciPy balanced-realization tests. on: push: @@ -60,18 +52,15 @@ jobs: - uses: actions/setup-python@v5 with: python-version: "3.11" - - name: Install dependencies shell: bash run: | python -m pip install --upgrade pip - [[ -f requirements.txt ]] && python -m pip install -r requirements.txt - + python -m pip install -r requirements-test.txt - name: Run tests shell: bash run: | [[ -d tests ]] && python -m unittest discover -s tests -p 'test_*.py' - - name: Smoke check shell: bash run: | @@ -83,7 +72,6 @@ jobs: timeout-minutes: 20 steps: - uses: actions/checkout@v4 - - name: Compile standalone C++ sources shell: bash run: | @@ -93,7 +81,6 @@ jobs: echo "No C++ sources on this branch." exit 0 fi - mkdir -p build/ci for source in "${sources[@]}"; do name="$(basename "${source%.cpp}")" @@ -101,7 +88,6 @@ jobs: echo "Compiling $source" g++ -std=c++17 -O2 -Wall -Wextra -Wpedantic -Werror=return-type \ -fopenmp "$source" -o "$binary" - if grep -q -- '--self-test' "$source"; then echo "Running $name --self-test" timeout 120 "$binary" --self-test diff --git a/.github/workflows/n580-recovery.yml b/.github/workflows/n580-recovery.yml index 34d99df2..0800b29c 100644 --- a/.github/workflows/n580-recovery.yml +++ b/.github/workflows/n580-recovery.yml @@ -4,6 +4,9 @@ on: branches: [main] paths: - 'scripts/recover_n580_covariance.py' + - 'scripts/projective_inference.py' + - 'tests/test_projective*.py' + - 'tests/test_n580_covariance_recovery.py' - '.github/workflows/n580-recovery.yml' permissions: contents: read @@ -19,7 +22,7 @@ jobs: - run: python -m pip install mpmath - name: Targeted mathematical tests run: | - python -m unittest discover -s tests -p 'test_projective_inference.py' + python -m unittest discover -s tests -p 'test_projective*.py' python -m unittest discover -s tests -p 'test_n580_covariance_recovery.py' - name: Recover and report without rerunning the experiment run: python scripts/recover_n580_covariance.py diff --git a/requirements-test.txt b/requirements-test.txt new file mode 100644 index 00000000..0080dc33 --- /dev/null +++ b/requirements-test.txt @@ -0,0 +1,5 @@ +-r requirements.txt +# Existing P598 balanced-realization tests import the scientific Python stack. +# Keep these test dependencies separate from the lightweight core runtime. +numpy>=1.24,<3 +scipy>=1.10,<2 diff --git a/scripts/projective_inference.py b/scripts/projective_inference.py index ccca1e62..7ee800d7 100644 --- a/scripts/projective_inference.py +++ b/scripts/projective_inference.py @@ -120,7 +120,11 @@ def _exact_support_fit(y,directions,kept,values,vectors,tol,roundoff): residual = W*(y-V*amplitudes) statistic = mp.fsum(x*x for x in residual) degrees = len(kept)-fitted_rank - if degrees == 0 and statistic < support_tol**2: + if degrees == 0: + # Saturation is an algebraic zero, not a chi-square fluctuation. + scale = max(1, mp.norm(W*y), mp.norm(W*V*amplitudes)) + if mp.norm(residual) > support_tol*scale: + raise ArithmeticError('saturated support solve lost numerical accuracy') statistic = mp.mpf(0) return statistic,degrees,amplitudes,constraint_rank,fitted_rank,support_tol @@ -168,6 +172,13 @@ def subspace_residual(observed: Sequence[Any], covariance: Sequence[Sequence[Any degrees = rank-width if degrees < 0: raise ValueError('model has more directions than covariance rank') + if degrees == 0: + # The design spans all observations. Do not turn roundoff into p=0. + tolerance = mp.power(10, -mp.dps//2) + scale = max(1, mp.norm(y), mp.norm(y-residual)) + if mp.norm(residual) > tolerance*scale: + raise ArithmeticError('saturated solve lost numerical accuracy') + statistic = mp.mpf(0) return {'statistic':float(statistic),'degrees_of_freedom':degrees, 'amplitudes':[float(x) for x in amplitudes], 'covariance_rank':rank, 'covariance_condition_number':float(condition) if condition is not None else None, diff --git a/tests/test_projective_saturated_support.py b/tests/test_projective_saturated_support.py new file mode 100644 index 00000000..b05f23fc --- /dev/null +++ b/tests/test_projective_saturated_support.py @@ -0,0 +1,58 @@ +"""Prevent saturated p=0 artefacts without clipping positive-df residuals.""" +import math +from pathlib import Path +import sys +import unittest +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / 'scripts')) +import projective_inference as p + + +class SaturatedSupportTests(unittest.TestCase): + def test_full_rank_saturation_is_an_algebraic_zero(self): + for scale in (1e-20, 1., 1e20): + with self.subTest(scale=scale): + result = p.subspace_residual( + [.3*scale, .7*scale], [[1, .23], [.23, 2]], + [[1, .123], [.47, .9]]) + self.assertEqual(result['degrees_of_freedom'], 0) + self.assertEqual(result['statistic'], 0.) + self.assertEqual(result['p_value'], 1.) + self.assertEqual(result['equivalent_sigma'], 0.) + + def test_exact_support_saturation(self): + result = p.subspace_residual([.3, .3, .7], + [[1, 1, 0], [1, 1, 0], [0, 0, 2]], + [[1, 1, .123], [.47, .47, .9]], covariance_mode='exact_support') + self.assertEqual(result['degrees_of_freedom'], 0) + self.assertEqual(result['statistic'], 0.) + self.assertEqual(result['p_value'], 1.) + + def test_zero_df_does_not_hide_support_incompatibility(self): + with self.assertRaisesRegex(ValueError, 'incompatible'): + p.subspace_residual([.3, .4, .7], + [[1, 1, 0], [1, 1, 0], [0, 0, 2]], + [[1, 1, .123], [.47, .47, .9]], covariance_mode='exact_support') + + def test_positive_df_residual_is_not_clipped(self): + result = p.ray_residual([1, 1e-15], [[1, 0], [0, 1]], [1, 0]) + self.assertEqual(result['degrees_of_freedom'], 1) + self.assertGreater(result['statistic'], 0.) + self.assertAlmostEqual(result['statistic']/1e-30, 1.) + + def test_covariance_support_and_psd_regressions(self): + s = [[1,0,0],[0,1,0],[0,0,0]] + with self.assertRaisesRegex(ValueError, 'incompatible'): + p.ray_residual([0,0,1],s,[1,0,0],covariance_mode='exact_support') + result=p.ray_residual([0,0,1],s,[1,0,1],covariance_mode='exact_support') + self.assertEqual((result['statistic'],result['degrees_of_freedom']), (1.,2)) + with self.assertRaisesRegex(ValueError, 'positive semidefinite'): + p.ray_residual([0,0,1],[[1,0,0],[0,1,0],[0,0,-.1]],[1,0,0]) + + def test_fieller_identity_unchanged(self): + y = [.3, .7]; s = [[1, .23], [.23, 2]] + for r in (1., 4., 16., -2.): + expected = (y[1]-r*y[0])**2/(s[1][1]-2*r*s[0][1]+r*r*s[0][0]) + self.assertAlmostEqual(p.ray_residual(y,s,[1,r])['statistic'], expected, places=14) + +if __name__ == '__main__': + unittest.main()