"""Generate paired median/IQR curves for the Quick overview.

Every rule sees the same 1,000 Gaussian racks at each sensor count. Labels
remain fixed within a rack across counts; response/noise are freshly drawn.
Historical lecture, dense-refinement, and ridge fixtures are never modified.
Optional --raw-output saves all scored draws outside the published payload.
"""
from pathlib import Path
import argparse
import hashlib
import json
import math
import platform
import time
import numpy as np

HERE = Path(__file__).resolve().parent
MASTER_SEED = 2026102901
RACKS = 1000
PRESETS = {"n5s4": (5, 4), "n20s1": (20, 1), "n5s25": (5, 25), "n20s40": (20, 40)}


def seed_for(*parts):
    name = ":".join(["sensor-overview-paired-v1", str(MASTER_SEED), *map(str, parts)])
    return int.from_bytes(hashlib.sha256(name.encode("ascii")).digest()[:8], "little")


def rng_for(*parts):
    return np.random.Generator(np.random.PCG64(seed_for(*parts)))


def raw_rack(key, rack, p):
    n, variance = PRESETS[key]
    y = rng_for(key, rack, "labels").standard_normal(n)
    a = rng_for(key, rack, p, "response").standard_normal(p)
    length = np.linalg.norm(a)
    if length == 0 or y @ y == 0:
        raise ArithmeticError("Degenerate draw: abort without replacing it")
    a *= math.sqrt(p) / length
    X = y[:, None] * a + math.sqrt(variance) * rng_for(key, rack, p, "noise").standard_normal((n, p))
    return X, y, a, variance


def vector_risk(w, a, variance):
    return float((w @ a - 1) ** 2 + variance * (w @ w))


def score_pair(X, y, a, variance):
    n, p = X.shape
    if p <= n:
        wi = np.linalg.lstsq(X, y, rcond=None)[0]
        interpolation = vector_risk(wi, a, variance)
    else:
        gram = X @ X.T
        alpha = np.linalg.solve(gram, y)
        interpolation = float((alpha @ (X @ a) - 1) ** 2 + variance * (alpha @ gram @ alpha))
    wc = X.T @ y / (p * (y @ y))
    correlator = vector_risk(wc, a, variance)
    if not np.isfinite([interpolation, correlator]).all() or min(interpolation, correlator) < 0:
        raise ArithmeticError("Invalid score: abort without filtering or replacing the rack")
    return interpolation, correlator


def summarize(values):
    lo, median, hi = np.quantile(values, [.25, .5, .75], axis=1, method="linear")
    return dict(lo=lo.tolist(), median=median.tolist(), hi=hi.tolist())


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--raw-output", type=Path)
    args = parser.parse_args()
    started = time.perf_counter()
    base = json.loads((HERE / "base.json").read_text())
    fixtures, raw, validation = {}, {}, []
    for key, (n, variance) in PRESETS.items():
        dense = range(4, 16) if n == 5 else range(10, 36)
        widths = sorted(set(base[key]["p"]) | set(dense))
        risks = np.empty((len(widths), RACKS, 2))
        for ip, p in enumerate(widths):
            for rack in range(RACKS):
                X, y, a, v = raw_rack(key, rack, p)
                risks[ip, rack] = score_pair(X, y, a, v)
                # Independent full-vector checks: direct square solve versus
                # production SVD at p=n; SVD versus Gram solve at wider counts.
                if rack < 3 and p in [n, n + 1, widths[-1]]:
                    weights = np.linalg.solve(X, y) if p == n else np.linalg.lstsq(X, y, rcond=None)[0]
                    direct = vector_risk(weights, a, v)
                    relative = abs(direct - risks[ip, rack, 0]) / max(1., abs(direct))
                    if relative > 2e-8:
                        raise ArithmeticError(f"Independent full-vector cross-check failed: {key}, p={p}, rack={rack}")
                    validation.append(dict(preset=key, p=p, rack=rack,
                        method="direct square solve versus SVD" if p == n else "SVD versus Gram solve", relative_error=relative))
        interp, corr = summarize(risks[:, :, 0]), summarize(risks[:, :, 1])
        fixtures[key] = dict(n=n, variance=variance, racks=RACKS, p=widths,
                             lo=interp["lo"], interp=interp["median"], hi=interp["hi"],
                             corrLo=corr["lo"], correlator=corr["median"], corrHi=corr["hi"])
        raw[key] = risks
        print(f"{key}: {len(widths)} counts × {RACKS} paired racks; {time.perf_counter()-started:.1f}s", flush=True)
    data = dict(schema="sensor-overview-paired-v1", masterSeed=MASTER_SEED, racks=RACKS,
                quantiles=[.25, .5, .75], presets=fixtures)
    payload = json.dumps(data, separators=(",", ":"), allow_nan=False)
    js = "// Generated by data/sensor-fixtures/generate-overview.py; all paired draws retained.\n" + \
        "(function(root,data){if(typeof module==='object'&&module.exports)module.exports=data;else root.SENSOR_OVERVIEW_DATA=data;})(typeof window!=='undefined'?window:globalThis," + payload + ");\n"
    output = HERE.parent.parent / "assets" / "sensor-overview-data.js"
    output.write_text(js)
    manifest = dict(schema=data["schema"], master_seed=MASTER_SEED, racks_per_count=RACKS,
        model="X=y a^T+sqrt(v) Z; iid standard Gaussian y,Z; normalized Gaussian a with norm²=p",
        seed_policy="First 8 SHA256 bytes, little-endian: sensor-overview-paired-v1:master:preset:rack[:p]:component",
        pairing="Both rules use exactly the same X,y,a at every rack/count; Gaussian y is shared across counts, a/noise refresh by count",
        selection_policy="All indices 0..999 retained, without conditioning label energy, filtering, replacement, clipping, winsorizing, or seed search; invalid arithmetic aborts",
        quantile_policy="NumPy quantile(method='linear'), q25/q50/q75, full float64 output; intervals describe racks, not confidence intervals",
        score="Population fresh MSE=(a·w−1)²+v||w||²; min-norm least squares and correlator w=X^T y/(p y^T y)",
        computation="Direct SVD at p<=n; equivalent Gram solve at p>n; 12 independent direct square-solve checks versus SVD at p=n, and 24 full-vector SVD checks versus Gram solve at n+1,last count for first three racks",
        runtime=dict(python=platform.python_version(), numpy=np.__version__),
        generator_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
        data_sha256=hashlib.sha256(output.read_bytes()).hexdigest(),
        presets={key: dict(n=c["n"], variance=c["variance"], p=c["p"], racks=c["racks"], paired_draws=len(c["p"])*RACKS) for key,c in fixtures.items()},
        validation=validation, elapsed_seconds=time.perf_counter()-started)
    if args.raw_output:
        np.savez_compressed(args.raw_output, **raw)
        manifest["raw_archive_sha256"] = hashlib.sha256(args.raw_output.read_bytes()).hexdigest()
    (HERE / "overview-manifest.json").write_text(json.dumps(manifest, indent=2) + "\n")
    print(f"Wrote {output.name}; {sum(len(c['p'])*RACKS for c in fixtures.values())} paired draws, {len(validation)} independent full-vector checks", flush=True)


if __name__ == "__main__":
    main()