1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5
  6def f(h, u, lam, q):
  7    return lam * h + u + q * h * h
  8
  9
 10def affine_rollout(lam, q, H, gamma):
 11    # Nominal h*=u*=0, one scalar initial-state coordinate gamma.
 12    hs = np.zeros(H + 1)
 13    Rs = np.zeros(H + 1)
 14    hs[0] = 0.0
 15    Rs[0] = 1.0
 16    for k in range(H):
 17        hs[k + 1] = f(hs[k], 0.0, lam, q)
 18        # A_k=df/dh at nominal = lam; C_k=1 but U_k=0.
 19        Rs[k + 1] = lam * Rs[k]
 20    pred = hs + Rs * gamma
 21    true = np.zeros(H + 1)
 22    true[0] = gamma
 23    for k in range(H):
 24        true[k + 1] = f(true[k], 0.0, lam, q)
 25    return hs, Rs, pred, true
 26
 27
 28def violation(lam, q, H, gamma):
 29    hs, Rs, pred, true = affine_rollout(lam, q, H, gamma)
 30    # v_k is one-step nonlinear defect evaluated on the affine family.
 31    vals = []
 32    for k in range(H):
 33        vals.append(abs(f(pred[k], 0.0, lam, q) - pred[k + 1]))
 34    return max(vals) if vals else 0.0
 35
 36
 37def max_box_violation(lam, q, H, radius):
 38    # In this scalar box, the endpoint is worst for this positive quadratic map.
 39    return max(violation(lam, q, H, -radius), violation(lam, q, H, radius))
 40
 41
 42def measured_horizon(lam, q, radius, eps, Hmax=80):
 43    good = 0
 44    for H in range(1, Hmax + 1):
 45        if max_box_violation(lam, q, H, radius) <= eps * (1 + 1e-12):
 46            good = H
 47        else:
 48            break
 49    return good
 50
 51
 52def predicted_horizon(lam, q, radius, eps):
 53    # For lambda>1: q*r^2*lambda^(2(H-1)) <= eps.
 54    if q * radius * radius <= eps:
 55        return max(1, int(math.floor(1 + math.log(eps/(q*radius*radius)) / (2*math.log(lam)))))
 56    return 0
 57
 58
 59def radius_for_horizon(lam, q, H, eps):
 60    # Exact one-step defect for the affine family at the last step.
 61    return math.sqrt(eps / (q * lam ** (2 * (H - 1))))
 62
 63
 64def main():
 65    np.random.seed(7); random.seed(7)
 66    lam, q, H, gamma = 1.25, 0.3, 6, 0.04
 67    # Core math sanity: autodiff-free finite difference Jacobian and recurrence.
 68    delta = 1e-6
 69    numeric_A = (f(delta, 0, lam, q) - f(-delta, 0, lam, q)) / (2*delta)
 70    _, Rs, _, _ = affine_rollout(lam, q, H, gamma)
 71    jacobian_check = {"analytic_A": lam, "finite_difference_A": numeric_A,
 72                      "max_R_error": float(np.max(np.abs(Rs - lam**np.arange(H+1))))}
 73
 74    # Prediction 1: quadratic scaling with perturbation radius.
 75    radii = np.array([0.01, 0.02, 0.04, 0.08])
 76    vals = np.array([max_box_violation(lam, q, H, r) for r in radii])
 77    log_slope = float(np.polyfit(np.log(radii), np.log(vals), 1)[0])
 78
 79    # Prediction 2: critical horizon follows the exponential boundary.
 80    eps = 1e-3
 81    hs = list(range(1, 13))
 82    measured = [measured_horizon(lam, q, r, eps, 30) for r in [0.01, 0.02, 0.04]]
 83    predicted = [predicted_horizon(lam, q, r, eps) for r in [0.01, 0.02, 0.04]]
 84
 85    # Prediction 3: trusted radius decays as lambda^-(H-1).
 86    horizons = np.array([2, 4, 6, 8])
 87    rr = np.array([radius_for_horizon(lam, q, int(h), eps) for h in horizons])
 88    decay_slope = float(np.polyfit(horizons - 1, np.log(rr), 1)[0])
 89    predicted_decay = -math.log(lam)
 90
 91    # Parameter sweep over lambda: stable dynamics should not show exponential growth.
 92    lambda_sweep = []
 93    for la in [0.8, 1.0, 1.1, 1.25, 1.5]:
 94        v = max_box_violation(la, q, 10, 0.04)
 95        lambda_sweep.append({"lambda": la, "H10_violation": v,
 96                             "critical_H": measured_horizon(la, q, 0.04, eps, 50)})
 97
 98    # Small monitor comparison: unmonitored affine prediction vs accepted trusted radius.
 99    # Baseline uses fixed radius 0.04; monitor shrinks to the largest radius satisfying eps.
100    baseline_v = max_box_violation(lam, q, 10, 0.04)
101    monitored_r = radius_for_horizon(lam, q, 10, eps)
102    monitored_v = max_box_violation(lam, q, 10, monitored_r)
103
104    out = {
105      "jacobian_check": jacobian_check,
106      "prediction_1_quadratic_radius": {"radii": radii.tolist(), "violations": vals.tolist(), "observed_log_slope": log_slope, "predicted_slope": 2.0},
107      "prediction_2_horizon_boundary": {"radii": [0.01,0.02,0.04], "epsilon": eps, "measured_H": measured, "predicted_H": predicted, "formula": "q*r^2*lambda^(2(H-1)) <= epsilon"},
108      "prediction_3_radius_decay": {"H": horizons.tolist(), "radii": rr.tolist(), "observed_log_decay_per_step": decay_slope, "predicted": predicted_decay},
109      "lambda_sweep": lambda_sweep,
110      "monitor_comparison": {"baseline_radius": 0.04, "baseline_H10_violation": baseline_v, "monitored_radius": monitored_r, "monitored_H10_violation": monitored_v},
111    }
112    Path("results.json").write_text(json.dumps(out, indent=2))
113    print(json.dumps(out, indent=2))
114
115if __name__ == "__main__":
116    main()