1import json
  2import math
  3import random
  4from pathlib import Path
  5
  6import numpy as np
  7
  8SEED = 982
  9np.random.seed(SEED)
 10random.seed(SEED)
 11
 12
 13def q_weights(q, alpha, K):
 14    """Finite truncation of (q^alpha;q)_k/(q;q)_k from the paper."""
 15    w = np.empty(K, dtype=float)
 16    w[0] = 1.0
 17    for k in range(1, K):
 18        w[k] = w[k - 1] * (1.0 - q ** (alpha + k - 1)) / (1.0 - q ** k)
 19    return w
 20
 21
 22def companion(g, weights):
 23    K = len(weights)
 24    M = np.zeros((K, K), dtype=float)
 25    M[0, :] = g * weights
 26    M[1:, :-1] = np.eye(K - 1)
 27    return M
 28
 29
 30def rho(g, weights):
 31    return float(np.max(np.abs(np.linalg.eigvals(companion(g, weights)))))
 32
 33
 34def controller(g, rho_hat, delta=0.05, eps=1e-8, g_max=10.0, eta_gain=0.01):
 35    target = 1.0 - delta
 36    # Conservative branch from the proposal; optional slow growth is only used
 37    # when the estimate is comfortably below the target.
 38    if rho_hat > target:
 39        return max(0.0, g * min(1.0, target / (rho_hat + eps)))
 40    return min(g_max, g * (1.0 + eta_gain))
 41
 42
 43def threshold_sweep(weights):
 44    predicted = 1.0 / float(np.sum(weights))
 45    gs = predicted * np.linspace(0.70, 1.30, 25)
 46    rs = np.array([rho(float(g), weights) for g in gs])
 47    # Linear interpolation of the closest sign change in rho-1.
 48    j = int(np.argmin(np.abs(rs - 1.0)))
 49    if j == 0 or j == len(gs) - 1:
 50        observed = float(gs[j])
 51    else:
 52        observed = float(gs[j - 1] + (1-rs[j-1]) * (gs[j]-gs[j-1])/(rs[j]-rs[j-1]))
 53    return predicted, observed, gs, rs
 54
 55
 56def gradient_slope(weights, g, n=160):
 57    M = companion(g, weights)
 58    # A generic perturbation avoids selecting a special eigenvector.
 59    v = np.linspace(1.0, 0.3, len(weights))
 60    vals = []
 61    for _ in range(n):
 62        vals.append(np.linalg.norm(v))
 63        v = M @ v
 64    vals = np.maximum(np.asarray(vals), 1e-300)
 65    # Ignore the transient and fit the asymptotic log slope.
 66    slope = float(np.polyfit(np.arange(40, n), np.log(vals[40:]), 1)[0])
 67    return slope
 68
 69
 70def run_controller(weights, g0, delta=0.05, steps=20):
 71    g = float(g0)
 72    rows = []
 73    for t in range(steps):
 74        r = rho(g, weights)
 75        rows.append((t, g, r))
 76        g = controller(g, r, delta=delta)
 77    return rows
 78
 79
 80def main():
 81    q, alpha, K = 0.7, 0.5, 30
 82    w = q_weights(q, alpha, K)
 83    predicted, observed, gs, rs = threshold_sweep(w)
 84
 85    # Prediction 1: positive-memory boundary is g*S=1 (z=1).
 86    boundary_rel_error = abs(observed - predicted) / predicted
 87    below = rho(0.90 * predicted, w)
 88    above = rho(1.10 * predicted, w)
 89
 90    # Prediction 2: near the boundary, rho crosses one monotonically.
 91    monotone = bool(np.all(np.diff(rs) > -1e-9))
 92    crossing = bool(below < 1.0 and above > 1.0)
 93
 94    # Prediction 3: gradient norm has asymptotic slope log(rho).
 95    g_stable = 0.80 * predicted
 96    g_unstable = 1.20 * predicted
 97    r_stable, r_unstable = rho(g_stable, w), rho(g_unstable, w)
 98    slope_stable = gradient_slope(w, g_stable)
 99    slope_unstable = gradient_slope(w, g_unstable)
100    slope_err = abs(slope_stable - math.log(r_stable))
101
102    # Controller starts 35% above the predicted boundary and should move to margin.
103    rows = run_controller(w, 1.35 * predicted, delta=0.05, steps=50)
104    final_g, final_r = rows[-1][1], rows[-1][2]
105    controller_target = 0.95
106    controller_ok = final_r <= controller_target + 0.005
107
108    # A tiny nonlinear rollout provides an interpretable baseline comparison.
109    # Both systems use the same initial state; controller rescales g when rho is high.
110    x0 = np.zeros(K); x0[0] = 1.0
111    def rollout(controlled):
112        state, g, maxabs = x0.copy(), 1.35 * predicted, 1.0
113        losses = []
114        for t in range(80):
115            r = rho(g, w)
116            if controlled:
117                g = controller(g, r, delta=0.05, eta_gain=0.0)
118            state = companion(g, w) @ state
119            maxabs = max(maxabs, float(np.max(np.abs(state))))
120            losses.append(float(np.sum(state * state)))
121        return maxabs, losses[-1]
122    base_max, base_loss = rollout(False)
123    ctrl_max, ctrl_loss = rollout(True)
124
125    report = {
126        "setup": {"q": q, "alpha": alpha, "K": K, "weight_sum": float(np.sum(w))},
127        "prediction_1_boundary": {
128            "predicted_gcrit_1_over_sum_w": predicted,
129            "observed_interpolated_gcrit": observed,
130            "relative_error": boundary_rel_error,
131            "pass_within_20pct": boundary_rel_error <= 0.20,
132        },
133        "prediction_2_transition": {
134            "rho_at_0.90_gcrit": below,
135            "rho_at_1.10_gcrit": above,
136            "crosses_unit_circle": crossing,
137            "rho_sweep_monotone": monotone,
138        },
139        "prediction_3_gradient_scaling": {
140            "stable_g": g_stable, "stable_rho": r_stable,
141            "measured_log_slope": slope_stable,
142            "predicted_log_slope": math.log(r_stable),
143            "absolute_slope_error": slope_err,
144            "unstable_g": g_unstable, "unstable_rho": r_unstable,
145            "unstable_measured_log_slope": slope_unstable,
146        },
147        "controller": {
148            "initial_g": rows[0][1], "initial_rho": rows[0][2],
149            "final_g": final_g, "final_rho": final_r,
150            "target_rho": controller_target, "pass": controller_ok,
151            "trajectory": [{"step": t, "g": g, "rho": r} for t, g, r in rows],
152        },
153        "mini_comparison": {
154            "uncontrolled_max_state": base_max, "uncontrolled_final_squared_norm": base_loss,
155            "controlled_max_state": ctrl_max, "controlled_final_squared_norm": ctrl_loss,
156        },
157        "mechanism_manifested": bool(boundary_rel_error <= .20 and crossing and monotone and slope_err <= .03 and controller_ok),
158    }
159    Path("results.json").write_text(json.dumps(report, indent=2))
160    print(json.dumps(report, indent=2))
161
162
163if __name__ == "__main__":
164    main()