1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2068
  6
  7def update_radius(radius, delta, gamma, alpha):
  8    """Adaptive conformal recursion with the optional nonnegative projection."""
  9    return max(0.0, radius + gamma * (float(delta > radius) - alpha))
 10
 11def conservative_radius(ensemble_radii, beta):
 12    """Empirical upper (1-beta) quantile across estimators."""
 13    return float(np.quantile(np.asarray(ensemble_radii), 1.0 - beta, method="linear"))
 14
 15def stationary_sweep():
 16    rng = np.random.default_rng(SEED)
 17    rows = []
 18    for alpha in (0.10, 0.20, 0.30, 0.40):
 19        for gamma in (0.01, 0.03):
 20            r = 0.5
 21            hits = []
 22            for t in range(120000):
 23                d = rng.exponential(1.0)
 24                hit = d > r
 25                r = update_radius(r, d, gamma, alpha)
 26                if t >= 20000:
 27                    hits.append(hit)
 28            observed = float(np.mean(hits))
 29            rows.append({"alpha": alpha, "gamma": gamma,
 30                         "predicted_exceedance": alpha,
 31                         "observed_exceedance": observed,
 32                         "abs_error": abs(observed - alpha)})
 33    return rows
 34
 35def ramp_sweep():
 36    # If every error exceeds r, dr/dt = gamma(1-alpha).
 37    rows = []
 38    alpha, r0, target = 0.2, 0.2, 1.4
 39    for gamma in (0.02, 0.04, 0.08, 0.16):
 40        r, n = r0, 0
 41        while r < target and n < 10000:
 42            r = update_radius(r, 100.0, gamma, alpha)
 43            n += 1
 44        predicted = math.ceil((target - r0) / (gamma * (1 - alpha)))
 45        rows.append({"gamma": gamma, "predicted_steps": predicted,
 46                     "observed_steps": n,
 47                     "relative_error": abs(n - predicted) / predicted})
 48    return rows
 49
 50def recovery_sweep():
 51    # If every error is below r, dr/dt = -gamma*alpha.
 52    rows = []
 53    alpha, r0, target = 0.2, 1.8, 0.4
 54    for gamma in (0.02, 0.04, 0.08, 0.16):
 55        r, n = r0, 0
 56        while r > target and n < 10000:
 57            r = update_radius(r, 0.0, gamma, alpha)
 58            n += 1
 59        predicted = math.ceil((r0 - target) / (gamma * alpha))
 60        rows.append({"gamma": gamma, "predicted_steps": predicted,
 61                     "observed_steps": n,
 62                     "relative_error": abs(n - predicted) / predicted})
 63    return rows
 64
 65def crossing_comparison():
 66    # Predicted clearance is 1.0 m. A shifted predictor error of .75 m makes
 67    # actual clearance .25 m, below the .55 m physical threshold. Benign error
 68    # is .15 m, giving .85 m actual clearance. Strictly use < for collisions.
 69    alpha, gamma = 0.2, 0.08
 70    errors = np.r_[np.full(60, .15), np.full(80, .75), np.full(60, .15)]
 71    actual_clearance = 1.0 - errors
 72    unsafe = actual_clearance < .55
 73    out = {}
 74    for name, kind in (("fixed_0.20", "fixed"), ("fixed_0.80", "fixed"),
 75                       ("adaptive", "adaptive")):
 76        r = .20
 77        accepted, radii = [], []
 78        for e in errors:
 79            if kind == "adaptive":
 80                radius_used = r
 81                r = update_radius(r, e, gamma, alpha)
 82            else:
 83                radius_used = float(name.split("_")[1])
 84            accepted.append(1.0 > .55 + radius_used)
 85            radii.append(radius_used)
 86        accepted = np.asarray(accepted, dtype=bool)
 87        out[name] = {
 88            "unsafe_frames": int(unsafe.sum()),
 89            "false_safe_frames": int(np.sum(accepted & unsafe)),
 90            "safe_rejections": int(np.sum((~accepted) & (~unsafe))),
 91            "mean_radius": float(np.mean(radii)),
 92            "accept_rate": float(np.mean(accepted))}
 93    # A direct check of the ensemble conservatism formula.
 94    out["quantile_check"] = {"radii": [0.2, 0.4, 0.8, 1.0], "beta": 0.25,
 95                             "upper_quantile": conservative_radius([0.2, 0.4, 0.8, 1.0], 0.25)}
 96    return out
 97
 98def main():
 99    stationary, ramp, recovery = stationary_sweep(), ramp_sweep(), recovery_sweep()
100    crossing = crossing_comparison()
101    coverage_ok = max(x["abs_error"] for x in stationary) < .02
102    ramp_ok = max(x["relative_error"] for x in ramp) == 0.0
103    recovery_ok = max(x["relative_error"] for x in recovery) == 0.0
104    result = {"seed": SEED, "stationary_coverage": stationary,
105              "post_shift_ramp": ramp, "recovery_decay": recovery,
106              "crossing_comparison": crossing,
107              "checks": {"stationary_within_0.02": coverage_ok,
108                         "ramp_formula_exact": ramp_ok,
109                         "recovery_formula_exact": recovery_ok,
110                         "mechanism_manifested": coverage_ok and ramp_ok and recovery_ok}}
111    Path("results.json").write_text(json.dumps(result, indent=2))
112    print(json.dumps(result, indent=2))
113
114if __name__ == "__main__":
115    main()