1"""Toy verification of excitation-gated neural calibration.
  2
  3The calibration model is the paper's local linearization:
  4    y_j = t + psi * S u_j + noise,
  5where t is an unknown translation nuisance and S is the 90-degree rotation.
  6After eliminating t, Fisher information for psi is S_u / sigma^2.
  7"""
  8import csv
  9import json
 10from pathlib import Path
 11import numpy as np
 12
 13SEED = 2192
 14RNG = np.random.default_rng(SEED)
 15
 16
 17def centered_spread(u):
 18    u = np.asarray(u, dtype=float)
 19    return float(np.sum((u - u.mean(axis=0)) ** 2))
 20
 21
 22def fisher_calibration(u, sigma):
 23    """Jacobian Fisher matrix for [translation_x, translation_y, yaw]."""
 24    u = np.asarray(u, dtype=float)
 25    J = np.zeros((len(u), 2, 3))
 26    J[:, :, :2] = np.eye(2)
 27    J[:, :, 2] = np.stack((-u[:, 1], u[:, 0]), axis=1)
 28    F = np.einsum("nai,nbj->ij", J, J) / (sigma * sigma)
 29    return F, J
 30
 31
 32def estimate_yaw(u, sigma, trials=1000, rng=None):
 33    """OLS estimate and empirical variance; translation is fitted, not known."""
 34    if rng is None:
 35        rng = np.random.default_rng(0)
 36    u = np.asarray(u, dtype=float)
 37    # Design rows for two coordinates: y = tx,ty + psi*S*u.
 38    X = np.zeros((2 * len(u), 3))
 39    for k, (ux, uy) in enumerate(u):
 40        X[2*k:2*k+2] = [[1, 0, -uy], [0, 1, ux]]
 41    true = np.array([0.37, -0.22, 0.15])
 42    pinv = np.linalg.pinv(X)
 43    estimates = []
 44    for _ in range(trials):
 45        estimates.append(pinv @ (X @ true + rng.normal(0, sigma, 2 * len(u))))
 46    estimates = np.asarray(estimates)
 47    return float(np.var(estimates[:, 2], ddof=1)), float(np.mean(estimates[:, 2]))
 48
 49
 50def inverse_information_sweep():
 51    sigma = 0.08
 52    L = 20
 53    rows = []
 54    # Amplitude controls spread while keeping the window shape fixed.
 55    for amp in [0.25, 0.5, 1.0, 2.0, 4.0]:
 56        x = np.linspace(-1, 1, L)[:, None]
 57        u = np.concatenate([amp * x, np.zeros_like(x)], axis=1)
 58        spread = centered_spread(u)
 59        empirical, _ = estimate_yaw(u, sigma, trials=1200, rng=np.random.default_rng(SEED + int(amp*10)))
 60        predicted = sigma * sigma / spread
 61        rows.append({"amplitude": amp, "spread": spread, "fisher": spread/sigma**2,
 62                     "empirical_var": empirical, "predicted_var": predicted,
 63                     "ratio_empirical_to_predicted": empirical/predicted})
 64    return rows
 65
 66
 67def threshold_sweep():
 68    sigma = 0.08
 69    L = 20
 70    epsilon = 0.10
 71    gamma = epsilon ** -2
 72    required_spread = gamma * sigma**2
 73    rows = []
 74    for amp in np.linspace(0.15, 1.50, 10):
 75        x = np.linspace(-1, 1, L)[:, None]
 76        u = np.concatenate([amp*x, np.zeros_like(x)], axis=1)
 77        spread = centered_spread(u)
 78        F, _ = fisher_calibration(u, sigma)
 79        # The 3x3 F contains translation gauge; its calibration Schur complement
 80        # is precisely the yaw information. This is the scalar certificate here.
 81        yaw_info = spread / sigma**2
 82        rows.append({"amplitude": float(amp), "spread": spread, "yaw_fisher": yaw_info,
 83                     "certified": bool(yaw_info >= gamma),
 84                     "predicted_certified": bool(spread >= required_spread),
 85                     "lambda_min_calibration_schur": yaw_info})
 86    return {"epsilon": epsilon, "gamma": gamma, "required_spread": required_spread, "rows": rows}
 87
 88
 89def rolling_window_sweep():
 90    """Verify that certification is forgotten exactly as old excitation leaves."""
 91    sigma = 0.08
 92    L = 20
 93    gamma = 100.0
 94    # First 20 samples have spread 20 (information 3125), then constant input.
 95    x = np.linspace(-1, 1, L)
 96    stream = np.r_[x, np.zeros(40)]
 97    cert = []
 98    infos = []
 99    for k in range(len(stream)):
100        window = stream[max(0, k-L+1):k+1, None]
101        info = centered_spread(window) / sigma**2 if len(window) > 1 else 0.0
102        infos.append(info)
103        cert.append(info >= gamma)
104    first_cert = next((i for i, c in enumerate(cert) if c), None)
105    first_loss = next((i for i in range(first_cert or 0, len(cert)) if not cert[i]), None)
106    # Predicted loss: once the initial spread exits, at k=L (zero-based).
107    return {"window": L, "first_certified_step": first_cert,
108            "observed_loss_step": first_loss, "predicted_loss_step": 2 * L - 1,
109            "information_at_peak": max(infos), "information_tail": infos[-1]}
110
111
112def acquisition_comparison():
113    """Compare decaying exploration with feedback-gated, alternating probes.
114
115    Alternating signs make the centered spread nonzero while e remains
116    orthogonal to the nominal task direction [1, 0]. Certification is checked
117    on every rolling window, so the gated controller re-excites after forgetting.
118    """
119    sigma, L, gamma, steps, trials = 0.08, 20, 100.0, 80, 300
120    outcomes = {"fixed_decay": [], "gated": []}
121    for trial in range(trials):
122        for method in outcomes:
123            us = []
124            first_cert = None
125            certified_count = 0
126            for k in range(steps):
127                recent = np.asarray(us[-L:])
128                info = centered_spread(recent) / sigma**2 if len(recent) > 1 else 0.0
129                certified = info >= gamma
130                if certified:
131                    certified_count += 1
132                    if first_cert is None:
133                        first_cert = k
134                # e alternates sign, preserving task direction while creating spread.
135                if method == "gated" and not certified:
136                    amp = 0.65
137                elif method == "fixed_decay":
138                    amp = 0.65 * np.exp(-k / 3.0)
139                else:
140                    amp = 0.0
141                us.append([1.0, amp * (1.0 if k % 2 == 0 else -1.0)])
142            outcomes[method].append({"first_cert": first_cert,
143                                     "certified_at_40": first_cert is not None and first_cert < 40,
144                                     "certified_fraction": certified_count / steps})
145    return {k: {"certification_rate_by_step_40": float(np.mean([x["certified_at_40"] for x in v])),
146                "median_first_certification": float(np.nanmedian([x["first_cert"] if x["first_cert"] is not None else np.nan for x in v])),
147                "mean_certified_fraction": float(np.mean([x["certified_fraction"] for x in v])),
148                "trials": trials} for k, v in outcomes.items()}
149
150
151def main():
152    result = {"seed": SEED, "inverse_information": inverse_information_sweep(),
153              "threshold": threshold_sweep(), "rolling_window": rolling_window_sweep(),
154              "acquisition": acquisition_comparison()}
155    Path("results.json").write_text(json.dumps(result, indent=2))
156    with open("inverse_information.csv", "w", newline="") as f:
157        rows = result["inverse_information"]
158        writer = csv.DictWriter(f, fieldnames=rows[0].keys())
159        writer.writeheader(); writer.writerows(rows)
160    print(json.dumps(result, indent=2))
161
162
163if __name__ == "__main__":
164    main()