1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6# Scalar center-manifold normal form:
  7# z_{k+1}=-(1+mu) z_k + a z_k^3.
  8# The fixed point has J=-(1+mu), and a symmetric 2-cycle satisfies
  9# f(z)=-z, giving |z|=sqrt(mu/a) for mu>0.
 10
 11def iterate(mu, a=1.0, z0=0.03, steps=4000):
 12    z = float(z0)
 13    hist = []
 14    for _ in range(steps):
 15        z = -(1.0 + mu) * z + a * z**3
 16        hist.append(z)
 17        if not np.isfinite(z) or abs(z) > 1e6:
 18            return np.asarray(hist), False
 19    return np.asarray(hist), True
 20
 21def cycle_amplitude(hist, tail=500):
 22    if len(hist) < tail or not np.all(np.isfinite(hist)):
 23        return float("nan")
 24    return float(np.mean(np.abs(hist[-tail:])))
 25
 26def a_k(hist, tail=500):
 27    if len(hist) < tail + 2 or not np.all(np.isfinite(hist)):
 28        return float("nan")
 29    return float(np.mean(np.abs(hist[-tail:] - hist[-tail-2:-2])) /
 30                 (np.mean(np.abs(hist[-tail:])) + 1e-12))
 31
 32def jacobian_eigenvalue(mu):
 33    return -(1.0 + mu)
 34
 35def corrected_guard_mu(mu, delta=0.05):
 36    # lambda >= -1+delta; lambda=-1-mu, so mu <= -delta.
 37    return min(mu, -delta)
 38
 39def run():
 40    a = 2.0
 41    delta = 0.05
 42
 43    # Prediction 1: the fixed point loses stability at mu=0 (lambda=-1).
 44    mus_boundary = np.linspace(-0.20, 0.20, 81)
 45    stable = []
 46    for mu in mus_boundary:
 47        h, ok = iterate(mu, a=a, z0=1e-4, steps=1500)
 48        stable.append(ok and abs(h[-1]) < 1e-3 and abs(h[-1] - h[-3]) < 1e-3)
 49    stable_mus = mus_boundary[np.asarray(stable)]
 50    unstable_mus = mus_boundary[~np.asarray(stable)]
 51    observed_boundary = float((stable_mus[-1] + unstable_mus[0]) / 2)
 52
 53    # Prediction 2: period-two amplitude scales as sqrt(mu/a).
 54    mu_amp = np.linspace(0.01, 0.18, 12)
 55    observed_amp, predicted_amp = [], []
 56    for mu in mu_amp:
 57        h, _ = iterate(mu, a=a, z0=0.01, steps=5000)
 58        observed_amp.append(cycle_amplitude(h))
 59        predicted_amp.append(math.sqrt(mu / a))
 60    observed_amp, predicted_amp = np.asarray(observed_amp), np.asarray(predicted_amp)
 61    rel_amp_error = float(np.mean(np.abs(observed_amp - predicted_amp) / predicted_amp))
 62    coeff = float(np.dot(np.sqrt(mu_amp), observed_amp) /
 63                  np.dot(np.sqrt(mu_amp), np.sqrt(mu_amp)))
 64
 65    # Prediction 3: amplitude changes as 1/sqrt(a).
 66    mu_scale = 0.08
 67    a_values = np.asarray([0.5, 1.0, 2.0, 4.0])
 68    scale_obs, scale_pred = [], []
 69    for av in a_values:
 70        h, _ = iterate(mu_scale, a=av, z0=0.01, steps=5000)
 71        scale_obs.append(cycle_amplitude(h))
 72        scale_pred.append(math.sqrt(mu_scale / av))
 73    scale_obs, scale_pred = np.asarray(scale_obs), np.asarray(scale_pred)
 74    scale_error = float(np.mean(np.abs(scale_obs-scale_pred)/scale_pred))
 75
 76    # Small same-task comparison: train scalar recurrent gain toward an
 77    # unstable target mu=+0.15. Baseline follows task loss; corrected guard
 78    # projects onto lambda >= -1+delta after each update.
 79    target, beta, lr = 0.15, 5.0, 0.08
 80    mu_base, mu_guard, mu_literal = -0.10, -0.10, -0.10
 81    for _ in range(100):
 82        mu_base -= lr * 2.0 * (mu_base-target)
 83        corrected_grad = 2.0 * beta * max(0.0, mu_guard + delta)
 84        mu_guard -= lr * (2.0*(mu_guard-target) + corrected_grad)
 85        mu_guard = corrected_guard_mu(mu_guard, delta)
 86        # Literal formula from the prompt: [max(0, lambda+1-delta)]^2.
 87        # Since lambda=-1-mu, this is max(0,-mu-delta)^2 and cannot stop
 88        # positive mu (the flip-unstable region).
 89        literal_grad = -2.0 * beta * max(0.0, -mu_literal-delta)
 90        mu_literal -= lr * (2.0*(mu_literal-target) + literal_grad)
 91
 92    hb, _ = iterate(mu_base, a=a, z0=0.01, steps=3000)
 93    hg, _ = iterate(corrected_guard_mu(mu_guard, delta), a=a, z0=0.01, steps=3000)
 94    hl, _ = iterate(mu_literal, a=a, z0=0.01, steps=3000)
 95    test_mu = 0.10
 96    lam = jacobian_eigenvalue(test_mu)
 97
 98    out = {
 99        "normal_form": {"a": a, "delta": delta},
100        "prediction_1_boundary": {
101            "predicted_mu": 0.0, "observed_mu_midpoint": observed_boundary,
102            "last_stable_grid_mu": float(stable_mus[-1]),
103            "first_unstable_grid_mu": float(unstable_mus[0]),
104            "lambda_at_observed_midpoint": jacobian_eigenvalue(observed_boundary)},
105        "prediction_2_sqrt_scaling": {
106            "mu": mu_amp.tolist(), "observed_amplitude": observed_amp.tolist(),
107            "predicted_amplitude": predicted_amp.tolist(),
108            "mean_relative_error": rel_amp_error,
109            "observed_sqrt_coefficient": coeff,
110            "predicted_sqrt_coefficient": 1.0/math.sqrt(a)},
111        "prediction_3_cubic_scaling": {
112            "a": a_values.tolist(), "observed_amplitude": scale_obs.tolist(),
113            "predicted_amplitude": scale_pred.tolist(), "mean_relative_error": scale_error},
114        "guard_training": {
115            "target_mu": target, "baseline_mu": float(mu_base),
116            "corrected_guard_mu": float(mu_guard), "literal_formula_mu": float(mu_literal),
117            "baseline_lambda": jacobian_eigenvalue(mu_base),
118            "corrected_guard_lambda": jacobian_eigenvalue(mu_guard),
119            "literal_formula_lambda": jacobian_eigenvalue(mu_literal),
120            "baseline_A_K": a_k(hb), "corrected_guard_A_K": a_k(hg),
121            "literal_formula_A_K": a_k(hl),
122            "baseline_cycle_amplitude": cycle_amplitude(hb),
123            "corrected_guard_cycle_amplitude": cycle_amplitude(hg),
124            "literal_formula_cycle_amplitude": cycle_amplitude(hl)},
125        "formula_sign_check": {
126            "mu": test_mu, "lambda": lam,
127            "displayed_penalty": max(0.0, lam + 1.0 - delta)**2,
128            "corrected_lower_bound_penalty": max(0.0, (-1.0 + delta) - lam)**2,
129            "note": "The displayed formula penalizes lambda above -1+delta, rather than lambda below it; corrected results are reported separately."}}
130    Path("results.json").write_text(json.dumps(out, indent=2))
131    print(json.dumps(out, indent=2))
132
133if __name__ == "__main__":
134    run()