1import json
  2import math
  3import random
  4from pathlib import Path
  5
  6import numpy as np
  7import torch
  8
  9SEED = 2466
 10random.seed(SEED)
 11np.random.seed(SEED)
 12torch.manual_seed(SEED)
 13A = 1.15
 14T = 60
 15
 16
 17def rollout(k, e0=1.0, a=A, T=T):
 18    rng = np.random.default_rng(SEED)
 19    v = rng.normal(0.0, 0.35, size=T)
 20    q, z = 0.0, e0
 21    errors = [z - q]
 22    for t in range(T):
 23        u = k * (q - z)
 24        q = a * q + v[t]
 25        z = a * z + u + v[t]
 26        errors.append(z - q)
 27    return np.asarray(errors)
 28
 29
 30def fit_slope(errors, burn=3):
 31    # Fit only the numerically resolved exponential segment. Exact synchronization
 32    # can otherwise turn roundoff into a spurious positive slope.
 33    vals = np.abs(errors[burn:])
 34    peak = float(np.max(vals))
 35    keep = vals > max(peak * 1e-12, 1e-14)
 36    vals = vals[keep]
 37    x = np.arange(len(errors))[burn:][keep].astype(float)
 38    if len(vals) < 2:
 39        return 0.0
 40    return float(np.polyfit(x, np.log(vals), 1)[0])
 41
 42
 43def analytic_lambda(k, a=A):
 44    m = abs(a - k)
 45    return -math.inf if m == 0 else math.log(m)
 46
 47
 48def numerical_checks():
 49    # Prediction 1: transverse stability boundary is k*=a-1 for positive gains.
 50    gains = np.linspace(-0.8, 2.8, 37)
 51    rows = []
 52    for k in gains:
 53        errors = rollout(float(k))
 54        measured = fit_slope(errors)
 55        predicted = analytic_lambda(float(k))
 56        rows.append({"k": float(k), "predicted_lambda": predicted,
 57                     "measured_slope": measured,
 58                     "stable_predicted": bool(abs(A-k) < 1),
 59                     "stable_observed": bool(measured < 0)})
 60
 61    # Prediction 2: measured log-residual slope equals transverse exponent.
 62    slope_errors = [abs(r["measured_slope"] - r["predicted_lambda"])
 63                    for r in rows if np.isfinite(r["predicted_lambda"]) and abs(A-r["k"]) > 1e-6]
 64    max_slope_error = max(slope_errors)
 65
 66    # Prediction 3: close to the boundary, decay time scales as 1/|lambda|.
 67    near_gains = [0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25]
 68    scaling = []
 69    for k in near_gains:
 70        lam = analytic_lambda(k)
 71        slope = fit_slope(rollout(k))
 72        scaling.append({"k": k, "abs_lambda": abs(lam),
 73                        "measured_abs_slope": abs(slope),
 74                        "time_constant_pred": 1.0/abs(lam),
 75                        "time_constant_obs": (1.0/abs(slope) if abs(slope) > 1e-15 else float("inf"))})
 76
 77    # Prediction 3: the convergence margin is linear in coupling gain,
 78    # lambda(k)=log|a-k|, with slope -1/(a-k) locally.
 79    scaling_gains = [0.20, 0.30, 0.40, 0.50]
 80    gain_scaling = [{"k": k, "predicted_multiplier": A-k,
 81                     "observed_multiplier": float(abs(rollout(k)[1] / rollout(k)[0]))}
 82                    for k in scaling_gains]
 83
 84    # Locate the observed sign transition by linear interpolation between sweep points.
 85    crossing = None
 86    for x, y in zip(rows[:-1], rows[1:]):
 87        if x["measured_slope"] * y["measured_slope"] <= 0:
 88            crossing = x["k"] + (0.0-x["measured_slope"]) * (y["k"]-x["k"]) / (y["measured_slope"]-x["measured_slope"])
 89            break
 90    stable_agreement = sum(r["stable_predicted"] == r["stable_observed"] for r in rows) / len(rows)
 91    return {"analytic_boundary": A - 1.0, "observed_boundary": crossing,
 92            "boundary_error": abs(crossing-(A-1.0)) if crossing is not None else None,
 93            "classification_agreement": stable_agreement,
 94            "max_abs_slope_error": max_slope_error,
 95            "slope_rows": rows,
 96            "scaling_rows": scaling, "gain_scaling": gain_scaling}
 97
 98
 99def train_controller(sync_weight, steps=500, batch=64, horizon=25):
100    # Same task for both methods: make plant track reference under shared forcing.
101    torch.manual_seed(SEED)
102    k = torch.nn.Parameter(torch.tensor(0.0))
103    opt = torch.optim.Adam([k], lr=0.03)
104    rng = torch.Generator().manual_seed(SEED)
105    losses = []
106    for _ in range(steps):
107        v = torch.randn(batch, horizon, generator=rng) * 0.35
108        q = torch.zeros(batch)
109        z = torch.randn(batch, generator=rng) * 1.5
110        total = torch.zeros(())
111        for t in range(horizon):
112            u = k * (q-z)
113            q = A*q + v[:, t]
114            z = A*z + u + v[:, t]
115            total = total + (z-q).square().mean()
116        # task and synchrony coincide here; sync_weight tests explicit residual shaping.
117        loss = total / horizon * (1.0 + sync_weight)
118        opt.zero_grad(); loss.backward(); opt.step()
119        losses.append(float(loss.detach()))
120    with torch.no_grad():
121        test = []
122        for e0 in [-2.0, -1.0, 1.0, 2.0]:
123            test.append(float(np.mean(rollout(float(k), e0=e0)[-10:]**2)))
124    return {"final_k": float(k.detach()), "train_loss_last": losses[-1],
125            "test_tail_mse": float(np.mean(test)), "loss_curve": losses}
126
127
128def main():
129    checks = numerical_checks()
130    baseline = train_controller(0.0)
131    idea = train_controller(2.0)
132    out = {"seed": SEED, "a": A, "horizon": T, "checks": checks,
133           "training": {"baseline_task_only": baseline,
134                         "idea_task_plus_synchrony": idea}}
135    Path("results.json").write_text(json.dumps(out, indent=2))
136    print(json.dumps({"analytic_boundary": checks["analytic_boundary"],
137                      "observed_boundary": checks["observed_boundary"],
138                      "boundary_error": checks["boundary_error"],
139                      "classification_agreement": checks["classification_agreement"],
140                      "max_abs_slope_error": checks["max_abs_slope_error"],
141                      "baseline": {k:v for k,v in baseline.items() if k != "loss_curve"},
142                      "idea": {k:v for k,v in idea.items() if k != "loss_curve"}}, indent=2))
143
144
145if __name__ == "__main__":
146    main()