1import json
  2import numpy as np
  3
  4
  5def spectral_radius(a):
  6    return float(np.max(np.abs(np.linalg.eigvals(a))))
  7
  8
  9def propagate(A, w, e0, horizon):
 10    e = np.asarray(e0, dtype=float).copy()
 11    out = [e.copy()]
 12    for _ in range(horizon):
 13        e = np.abs(A) @ e + w
 14        out.append(e.copy())
 15    return np.asarray(out)
 16
 17
 18def fit_log_slope(values, start, stop):
 19    y = np.asarray(values[start:stop])
 20    x = np.arange(start, stop)
 21    return float(np.polyfit(x, np.log(np.maximum(y, 1e-300)), 1)[0])
 22
 23
 24def main():
 25    np.set_printoptions(precision=6, suppress=True)
 26    # Positive transition: the absolute-value bound is exact for aligned errors.
 27    B = np.array([[0.72, 0.18], [0.10, 0.48]], dtype=float)
 28    rho0 = spectral_radius(B)
 29    w = np.array([0.01, 0.02])
 30    e0 = np.array([0.03, 0.01])
 31    H = 120
 32    alpha_star = 1.0 / rho0
 33
 34    # Prediction 1: the boundary is alpha*rho(B)=1.
 35    alphas = alpha_star * np.array([0.80, 0.95, 0.99, 1.00, 1.01, 1.05, 1.20])
 36    sweep = []
 37    for alpha in alphas:
 38        A = alpha * B
 39        traj = propagate(A, w, e0, H)
 40        total = traj.sum(axis=1)
 41        r = alpha * rho0
 42        row = {"alpha": float(alpha), "rho_abs": float(r),
 43               "final_total": float(total[-1]),
 44               "growth_ratio_last": float(total[-1] / max(total[-2], 1e-30))}
 45        if r < 1:
 46            fixed = np.linalg.solve(np.eye(2) - A, w)
 47            row["predicted_plateau_total"] = float(fixed.sum())
 48            row["plateau_relative_error"] = float(np.max(np.abs(traj[-1] - fixed)) /
 49                                                    np.max(np.abs(fixed)))
 50        else:
 51            row["predicted_log_growth"] = float(np.log(r))
 52            row["observed_log_growth"] = fit_log_slope(total, 70, 120)
 53        sweep.append(row)
 54
 55    # Prediction 2: below one, finite tube scales as (1-rho)^-1 near boundary.
 56    near = []
 57    for r in [0.50, 0.70, 0.80, 0.90, 0.95, 0.98]:
 58        A = (r / rho0) * B
 59        fixed = np.linalg.solve(np.eye(2) - A, w)
 60        near.append({"rho": r, "plateau_total": float(fixed.sum()),
 61                     "scaled_plateau": float((1-r) * fixed.sum())})
 62
 63    # Prediction 3: output margin is linear in ebar for y = [1, -0.5] h.
 64    G = np.array([1.0, 0.5])
 65    A = 0.75 / rho0 * B
 66    tr = propagate(A, w, e0, 40)
 67    margins = tr @ G
 68    ratios = margins / np.maximum(np.linalg.norm(tr, axis=1, ord=1), 1e-30)
 69    # Exact validity check with random signed disturbances bounded by w.
 70    rng = np.random.default_rng(7)
 71    actual = np.zeros(2)
 72    bound = np.zeros(2)
 73    max_violation = 0.0
 74    for _ in range(1000):
 75        actual[:] = e0
 76        bound[:] = e0
 77        for j in range(30):
 78            disturbance = rng.uniform(-w, w)
 79            actual[:] = A @ actual + disturbance
 80            bound[:] = np.abs(A) @ bound + w
 81            max_violation = max(max_violation, float(np.max(np.abs(actual) - bound)))
 82
 83    # Mini-experiment: horizon-dependent local tube versus a single worst-case tube.
 84    # The latter is the standard conservative replacement A_j -> elementwise max_j A_j.
 85    Htv = 40
 86    scales = 0.72 + 0.18 * np.sin(np.arange(Htv) * 0.55) ** 2
 87    As = [float(sc) / 0.78 * B for sc in scales]
 88    local = np.zeros(2)
 89    local_path = [local.copy()]
 90    for Aj in As:
 91        local = np.abs(Aj) @ local + w
 92        local_path.append(local.copy())
 93    local_path = np.asarray(local_path)
 94    M = np.max(np.asarray(As), axis=0)
 95    shared = propagate(M, w, np.zeros(2), Htv)
 96    local_total = local_path.sum(axis=1)
 97    shared_total = shared.sum(axis=1)
 98    result_extra = {
 99        "time_varying_local_vs_shared": {
100            "rho_local_max": float(max(spectral_radius(Aj) for Aj in As)),
101            "rho_shared_worst_case": spectral_radius(M),
102            "local_final_tube_sum": float(local_total[-1]),
103            "shared_final_tube_sum": float(shared_total[-1]),
104            "shared_over_local_final": float(shared_total[-1] / local_total[-1]),
105            "mean_shared_over_local": float(np.mean(shared_total[1:] / np.maximum(local_total[1:], 1e-30))),
106            "local_tube_sums_first_8": [float(x) for x in local_total[:8]],
107            "shared_tube_sums_first_8": [float(x) for x in shared_total[:8]]
108        }
109    }
110
111    result = {
112        **result_extra,
113        "base_rho": rho0,
114        "predicted_alpha_boundary": alpha_star,
115        "sweep": sweep,
116        "near_boundary_plateau_scaling": near,
117        "output_margin_linear_ratio_min": float(np.min(ratios)),
118        "output_margin_linear_ratio_max": float(np.max(ratios)),
119        "random_disturbance_max_bound_violation": max_violation,
120        "formula": "ebar[j+1]=abs(A)@ebar[j]+w; fixed=(I-A)^(-1)w"
121    }
122    with open("results.json", "w") as f:
123        json.dump(result, f, indent=2)
124    print(json.dumps(result, indent=2))
125
126
127if __name__ == "__main__":
128    main()