# How the search window affects the search p-value: the bundle's 1985..min(2015, end-10)
# against a Beaulieu-style window trimmed 10% at each end of 1970..end. Plug-in AR(1) null as in the bundle.
import csv, sys, math
import numpy as np
rows = list(csv.DictReader(open(sys.argv[1])))
years = np.array([int(r["year"]) for r in rows]); y = np.array([float(r["anomaly_c"]) for r in rows])
rng = np.random.Generator(np.random.PCG64(11))
def basis(yrs, knot=None):
    cols = [np.ones(len(yrs)), (yrs - 1970) / 10.0]
    if knot is not None: cols.append(np.maximum(0, yrs - knot) / 10.0)
    return np.column_stack(cols)
def scanner(yrs, knots):
    x0 = basis(yrs); P = x0 @ np.linalg.pinv(x0)
    H = np.maximum(0, yrs[:, None] - knots[None, :]) / 10.0
    Hr = H - P @ H; U = Hr / np.linalg.norm(Hr, axis=0)
    def f(V):
        R = V - V @ P.T; s0 = np.sum(R * R, axis=1); imp = (R @ U) ** 2
        st = imp / (s0[:, None] - imp) * (V.shape[1] - 3)
        return np.max(st, axis=1), knots[np.argmax(st, axis=1)]
    return f
def ar1(n, rho, count):
    e = rng.standard_normal((count, n)); out = np.empty_like(e)
    out[:, 0] = e[:, 0] / math.sqrt(1 - rho * rho)
    for i in range(1, n): out[:, i] = rho * out[:, i - 1] + e[:, i]
    return out
for end in (2025,):
    m = years <= end; yr, v = years[m], y[m]; n = len(v)
    r = v - basis(yr) @ np.linalg.lstsq(basis(yr), v, rcond=None)[0]
    rho = float(r[1:] @ r[:-1] / (r[:-1] @ r[:-1])); rho = {2023: rho, 2024: rho, 2025: 0.3607}[end]
    noise = ar1(n, rho, 20000)
    for name, knots in (("bundle 1985..min(2015,end-10)", np.arange(1985, min(2015, end - 10) + 1)),
                        ("10% trimmed", np.arange(yr[int(math.ceil(0.1 * n))], yr[int(math.floor(0.9 * n))] + 1))):
        sc = scanner(yr, knots); f, k = sc(v[None, :]); s, _ = sc(noise)
        print(f"through {end} [{name}: {knots[0]}-{knots[-1]}]: maxF {f[0]:.3f} at {k[0]}, rho {rho:.3f}, p {(1 + np.count_nonzero(s >= f[0])) / 20001:.4f}")