1import json
  2import numpy as np
  3from pathlib import Path
  4
  5
  6def controller(gamma, H, S, H0, S0, rho=0.10, gmin=1e-7, gmax=0.05, eps0=1e-6):
  7    score = (H-H0)/max(H0, eps0) + (S-S0)/max(S0, eps0)
  8    return float(np.clip(gamma*np.exp(rho*score), gmin, gmax)), float(score)
  9
 10
 11def ridge_fit(X, y, gamma):
 12    n, d = X.shape
 13    A = X.T @ X / n + gamma*np.eye(d)
 14    return np.linalg.solve(A, X.T @ y / n)
 15
 16
 17def diagnostics(X, eta=0.30, eps=1e-10):
 18    n = X.shape[0]
 19    lam = np.linalg.eigvalsh(X.T @ X / n + eps*np.eye(X.shape[1]))
 20    inv = 1.0/lam
 21    H = float(np.mean(inv))
 22    M = float(np.mean(np.where(lam >= eta, inv, 0.0)))
 23    S = float(np.mean(np.where(lam < eta, inv, 0.0)))
 24    return H, M, S, float(lam.min()), lam
 25
 26
 27def make_data(rng, n, d, hard_edge=False):
 28    # A controlled spectrum: only the final coordinates become nearly unobserved.
 29    eig = np.ones(d)
 30    if hard_edge:
 31        eig[-d//4:] = np.geomspace(1e-4, 2e-2, d//4)
 32    else:
 33        eig[-d//4:] = np.geomspace(0.35, 1.0, d//4)
 34    X = rng.normal(size=(n,d)) * np.sqrt(eig)[None,:]
 35    beta = rng.normal(size=d) / np.sqrt(d)
 36    y = X @ beta + 0.10*rng.normal(size=n)
 37    return X, y, beta, eig
 38
 39
 40def verify_math():
 41    rng = np.random.default_rng(123)
 42    X, _, _, _ = make_data(rng, 128, 24, hard_edge=True)
 43    H, M, S, _, lam = diagnostics(X, eta=0.30)
 44    split_error = abs(H-(M+S))
 45    # The closed-form ridge solution is also checked against its normal equations.
 46    w = ridge_fit(X, rng.normal(size=128), 0.07)
 47    normal_eq = np.linalg.norm((X.T@X/128 + .07*np.eye(24))@w - X.T@rng.normal(size=128))
 48    # use a fresh explicit RHS for the actual residual check
 49    yy = rng.normal(size=128); ww = ridge_fit(X, yy, .07)
 50    residual = np.linalg.norm((X.T@X/128 + .07*np.eye(24))@ww - X.T@yy/128)
 51    return {"H":H, "M":M, "S":S, "decomposition_abs_error":split_error,
 52            "min_eigenvalue":float(lam.min()), "normal_equation_residual":float(residual)}
 53
 54
 55def run():
 56    rng = np.random.default_rng(2025)
 57    n, d, ntest, reps = 192, 48, 512, 80
 58    # Calibration represents the healthy-spectrum regime, as prescribed by the controller.
 59    cal_H, cal_S = [], []
 60    for _ in range(20):
 61        X, _, _, _ = make_data(rng, n, d, hard_edge=False)
 62        h, _, s, _, _ = diagnostics(X, eta=0.30)
 63        cal_H.append(h); cal_S.append(s)
 64    H0, S0 = float(np.median(cal_H)), float(np.median(cal_S))
 65    rows = {k: [] for k in ["ridgeless","fixed_ridge","adaptive"]}
 66    gammas, scores, hvals = [], [], []
 67    for hard in [False, True]:
 68        for _ in range(reps):
 69            X, y, beta, _ = make_data(rng, n, d, hard_edge=hard)
 70            H, M, S, mineig, _ = diagnostics(X, eta=0.30)
 71            # Initialize weakly; update once from this feature batch.
 72            ga, score = controller(1e-4, H, S, H0, S0)
 73            for name, gamma in [("ridgeless",0.0),("fixed_ridge",0.01),("adaptive",ga)]:
 74                w = ridge_fit(X, y, gamma)
 75                # Test corruption is additive feature noise, exposing unstable directions.
 76                Xt, _, _, _ = make_data(rng, ntest, d, hard_edge=hard)
 77                yt = Xt @ beta + 0.10*rng.normal(size=ntest)
 78                pred = Xt @ w
 79                losses = (pred-yt)**2
 80                rows[name].append({"hard":hard, "mse":float(losses.mean()),
 81                                   "p99":float(np.quantile(losses,.99)),
 82                                   "max_output":float(np.max(np.abs(pred))),
 83                                   "mineig":mineig})
 84            if hard: gammas.append(ga); scores.append(score); hvals.append(H)
 85    summary = {}
 86    for name, vals in rows.items():
 87        summary[name] = {}
 88        for regime, label in [(False,"healthy"),(True,"hard_edge")]:
 89            a = [v for v in vals if v["hard"] == regime]
 90            summary[name][label] = {k:float(np.mean([x[k] for x in a])) for k in ["mse","p99","max_output","mineig"]}
 91    # Robustness signal is the hard-edge tail loss; report paired regime counts and controller behavior.
 92    return {"math":verify_math(), "calibration":{"H0":H0,"S0":S0},
 93            "summary":summary, "adaptive_hard_gamma":{"median":float(np.median(gammas)),
 94            "mean":float(np.mean(gammas)),"min":float(np.min(gammas)),"max":float(np.max(gammas)),
 95            "median_score":float(np.median(scores)),"median_H":float(np.median(hvals))}}
 96
 97
 98if __name__ == "__main__":
 99    out = run()
100    print(json.dumps(out, indent=2))
101    Path("results.json").write_text(json.dumps(out, indent=2))