"""Reckoning 001, follow-up: sensitivity to the drift model (criticism by Tern, 2026-10-24).
For each row (and the session means) fit JOINTLY: polynomial drift of degree p in the mark index k (k=0..16, 17 readings,
the 17th is the 1st mark again) + the harmonic-2 terms cos(2th), sin(2th), th = 22.5deg*k.  Joint fitting means drift and
signal cannot leak into each other except through the real degeneracy of the sampling.  Report h2 and its 1-sigma from the
fit residual, for p = 1..4.  Then INJECT a synthetic signal of the expected amplitude (0.184 fringe) with random phase
into the same rows and report how much of it each model recovers, so over-aggressive drift models can be seen to eat signal.
Run: python3 mm1887_sensitivity.py
"""
import json, math, numpy as np
d = json.load(open('data/mm1887.json')); A_exp = 0.184
k = np.arange(17); th = np.radians(22.5) * k
def design(p, harmonics=(2,)):
    cols = [k.astype(float) ** j for j in range(p + 1)]
    for h in harmonics: cols += [np.cos(h * th), np.sin(h * th)]
    return np.array(cols).T
def fit(y, p, harmonics=(2,)):
    X = design(p, harmonics); beta, res, rank, _ = np.linalg.lstsq(X, y, rcond=None)
    r = y - X @ beta; dof = len(y) - X.shape[1]; s2 = (r @ r) / dof
    cov = s2 * np.linalg.inv(X.T @ X); i = p + 1       # index of the first harmonic column pair
    a, b = beta[i], beta[i + 1]; amp = math.hypot(a, b)
    # sigma of amplitude ~ sqrt of the mean of the variances of a and b (adequate for a summary)
    sig = math.sqrt((cov[i, i] + cov[i + 1, i + 1]) / 2)
    return amp, sig, math.sqrt(s2)
rows = {}
for sess in ('noon', 'evening'):
    for day, x in d[sess].items(): rows[sess + ' ' + day] = np.array(x) * 0.02
    rows[sess + ' MEAN'] = np.mean([np.array(x) * 0.02 for x in d[sess].values()], axis=0)
print('h2 amplitude (fringe) [+/- 1 sigma] under drift polynomial of degree p (joint fit):')
print('%-18s' % 'row' + ''.join('   p=%d            ' % p for p in range(1, 5)))
for name, y in rows.items():
    cells = []
    for p in range(1, 5):
        a, s, r = fit(y, p); cells.append('%.3f +/- %.3f   ' % (a, s))
    print('%-18s' % name + ' '.join(cells))
print('\nExpected for an orbital-speed ether wind: %.3f. Fraction expected, by model, session means:' % A_exp)
for name in ('noon MEAN', 'evening MEAN'):
    print('  %-14s' % name + '  '.join('p=%d: %.2f' % (p, fit(rows[name], p)[0] / A_exp) for p in range(1, 5)))
print('\nInjection test: add %.3f cos(2th+phi) (random phi, 400 draws) to each session MEAN row, recover amplitude by model:' % A_exp)
rng = np.random.default_rng(7)
for name in ('noon MEAN', 'evening MEAN'):
    for p in range(1, 5):
        rec = []
        for _ in range(400):
            ph = rng.uniform(0, 2 * math.pi); y = rows[name] + A_exp * np.cos(2 * th + ph); rec.append(fit(y, p)[0])
        base = fit(rows[name], p)[0]
        print('  %-14s p=%d  recovered mean %.3f (sd %.3f); no-injection %.3f' % (name, p, np.mean(rec), np.std(rec), base))