"""Toy verification of excitation-gated neural calibration. The calibration model is the paper's local linearization: y_j = t + psi * S u_j + noise, where t is an unknown translation nuisance and S is the 90-degree rotation. After eliminating t, Fisher information for psi is S_u / sigma^2. """ import csv import json from pathlib import Path import numpy as np SEED = 2192 RNG = np.random.default_rng(SEED) def centered_spread(u): u = np.asarray(u, dtype=float) return float(np.sum((u - u.mean(axis=0)) ** 2)) def fisher_calibration(u, sigma): """Jacobian Fisher matrix for [translation_x, translation_y, yaw].""" u = np.asarray(u, dtype=float) J = np.zeros((len(u), 2, 3)) J[:, :, :2] = np.eye(2) J[:, :, 2] = np.stack((-u[:, 1], u[:, 0]), axis=1) F = np.einsum("nai,nbj->ij", J, J) / (sigma * sigma) return F, J def estimate_yaw(u, sigma, trials=1000, rng=None): """OLS estimate and empirical variance; translation is fitted, not known.""" if rng is None: rng = np.random.default_rng(0) u = np.asarray(u, dtype=float) # Design rows for two coordinates: y = tx,ty + psi*S*u. X = np.zeros((2 * len(u), 3)) for k, (ux, uy) in enumerate(u): X[2*k:2*k+2] = [[1, 0, -uy], [0, 1, ux]] true = np.array([0.37, -0.22, 0.15]) pinv = np.linalg.pinv(X) estimates = [] for _ in range(trials): estimates.append(pinv @ (X @ true + rng.normal(0, sigma, 2 * len(u)))) estimates = np.asarray(estimates) return float(np.var(estimates[:, 2], ddof=1)), float(np.mean(estimates[:, 2])) def inverse_information_sweep(): sigma = 0.08 L = 20 rows = [] # Amplitude controls spread while keeping the window shape fixed. for amp in [0.25, 0.5, 1.0, 2.0, 4.0]: x = np.linspace(-1, 1, L)[:, None] u = np.concatenate([amp * x, np.zeros_like(x)], axis=1) spread = centered_spread(u) empirical, _ = estimate_yaw(u, sigma, trials=1200, rng=np.random.default_rng(SEED + int(amp*10))) predicted = sigma * sigma / spread rows.append({"amplitude": amp, "spread": spread, "fisher": spread/sigma**2, "empirical_var": empirical, "predicted_var": predicted, "ratio_empirical_to_predicted": empirical/predicted}) return rows def threshold_sweep(): sigma = 0.08 L = 20 epsilon = 0.10 gamma = epsilon ** -2 required_spread = gamma * sigma**2 rows = [] for amp in np.linspace(0.15, 1.50, 10): x = np.linspace(-1, 1, L)[:, None] u = np.concatenate([amp*x, np.zeros_like(x)], axis=1) spread = centered_spread(u) F, _ = fisher_calibration(u, sigma) # The 3x3 F contains translation gauge; its calibration Schur complement # is precisely the yaw information. This is the scalar certificate here. yaw_info = spread / sigma**2 rows.append({"amplitude": float(amp), "spread": spread, "yaw_fisher": yaw_info, "certified": bool(yaw_info >= gamma), "predicted_certified": bool(spread >= required_spread), "lambda_min_calibration_schur": yaw_info}) return {"epsilon": epsilon, "gamma": gamma, "required_spread": required_spread, "rows": rows} def rolling_window_sweep(): """Verify that certification is forgotten exactly as old excitation leaves.""" sigma = 0.08 L = 20 gamma = 100.0 # First 20 samples have spread 20 (information 3125), then constant input. x = np.linspace(-1, 1, L) stream = np.r_[x, np.zeros(40)] cert = [] infos = [] for k in range(len(stream)): window = stream[max(0, k-L+1):k+1, None] info = centered_spread(window) / sigma**2 if len(window) > 1 else 0.0 infos.append(info) cert.append(info >= gamma) first_cert = next((i for i, c in enumerate(cert) if c), None) first_loss = next((i for i in range(first_cert or 0, len(cert)) if not cert[i]), None) # Predicted loss: once the initial spread exits, at k=L (zero-based). return {"window": L, "first_certified_step": first_cert, "observed_loss_step": first_loss, "predicted_loss_step": 2 * L - 1, "information_at_peak": max(infos), "information_tail": infos[-1]} def acquisition_comparison(): """Compare decaying exploration with feedback-gated, alternating probes. Alternating signs make the centered spread nonzero while e remains orthogonal to the nominal task direction [1, 0]. Certification is checked on every rolling window, so the gated controller re-excites after forgetting. """ sigma, L, gamma, steps, trials = 0.08, 20, 100.0, 80, 300 outcomes = {"fixed_decay": [], "gated": []} for trial in range(trials): for method in outcomes: us = [] first_cert = None certified_count = 0 for k in range(steps): recent = np.asarray(us[-L:]) info = centered_spread(recent) / sigma**2 if len(recent) > 1 else 0.0 certified = info >= gamma if certified: certified_count += 1 if first_cert is None: first_cert = k # e alternates sign, preserving task direction while creating spread. if method == "gated" and not certified: amp = 0.65 elif method == "fixed_decay": amp = 0.65 * np.exp(-k / 3.0) else: amp = 0.0 us.append([1.0, amp * (1.0 if k % 2 == 0 else -1.0)]) outcomes[method].append({"first_cert": first_cert, "certified_at_40": first_cert is not None and first_cert < 40, "certified_fraction": certified_count / steps}) return {k: {"certification_rate_by_step_40": float(np.mean([x["certified_at_40"] for x in v])), "median_first_certification": float(np.nanmedian([x["first_cert"] if x["first_cert"] is not None else np.nan for x in v])), "mean_certified_fraction": float(np.mean([x["certified_fraction"] for x in v])), "trials": trials} for k, v in outcomes.items()} def main(): result = {"seed": SEED, "inverse_information": inverse_information_sweep(), "threshold": threshold_sweep(), "rolling_window": rolling_window_sweep(), "acquisition": acquisition_comparison()} Path("results.json").write_text(json.dumps(result, indent=2)) with open("inverse_information.csv", "w", newline="") as f: rows = result["inverse_information"] writer = csv.DictWriter(f, fieldnames=rows[0].keys()) writer.writeheader(); writer.writerows(rows) print(json.dumps(result, indent=2)) if __name__ == "__main__": main()