1import json
  2from pathlib import Path
  3import numpy as np
  4
  5
  6def ar1_stats(a, sigma=1.0, n=30000, burn=3000, seed=0):
  7    rng = np.random.default_rng(seed)
  8    x = np.zeros(n + burn)
  9    for t in range(n + burn - 1):
 10        x[t + 1] = a * x[t] + sigma * rng.normal()
 11    y = x[burn:]
 12    return float(np.var(y)), float(np.corrcoef(y[:-1], y[1:])[0, 1])
 13
 14
 15class CSDMonitor:
 16    def __init__(self, window=80, ac_threshold=0.94, gamma=0.9):
 17        self.window = window
 18        self.ac_threshold = ac_threshold
 19        self.gamma = gamma
 20        self.values = []
 21        self.previous_variance = None
 22        self.last = {"a": np.nan, "variance": np.nan, "trigger": False}
 23
 24    def update(self, value):
 25        self.values.append(float(value))
 26        if len(self.values) < self.window:
 27            return 1.0, False
 28        z = np.asarray(self.values[-self.window:])
 29        z = z - z.mean()
 30        variance = float(np.mean(z * z))
 31        denom = float(np.dot(z[:-1], z[:-1]))
 32        a_hat = float(np.dot(z[:-1], z[1:]) / denom) if denom > 1e-12 else 0.0
 33        rising = self.previous_variance is not None and variance > self.previous_variance
 34        trigger = bool(a_hat > self.ac_threshold and rising)
 35        self.previous_variance = variance
 36        self.last = {"a": a_hat, "variance": variance, "trigger": trigger}
 37        multiplier = float(np.exp(-self.gamma * max(a_hat - self.ac_threshold, 0.0))) if trigger else 1.0
 38        return multiplier, trigger
 39
 40
 41def variance_and_ac_checks():
 42    # Prediction 1: stationary variance is sigma^2/(1-a^2).
 43    # Prediction 2: lag-1 correlation equals a.
 44    rows = []
 45    for i, a in enumerate([0.20, 0.50, 0.70, 0.85, 0.93]):
 46        var, rho = ar1_stats(a, seed=100 + i)
 47        pred_var = 1.0 / (1.0 - a * a)
 48        rows.append({"a": a, "predicted_variance": pred_var, "observed_variance": var,
 49                     "predicted_rho": a, "observed_rho": rho,
 50                     "variance_relative_error": abs(var - pred_var) / pred_var,
 51                     "rho_abs_error": abs(rho - a)})
 52    return rows
 53
 54
 55def boundary_check():
 56    # Prediction 3: deterministic AR(1) is stable below a=1 and grows above a=1.
 57    # Estimate the boundary by classifying a from the final/initial amplitude ratio.
 58    candidates = np.linspace(0.94, 1.06, 25)
 59    classified = []
 60    for a in candidates:
 61        x = 1.0
 62        for _ in range(80):
 63            x *= a
 64        classified.append(abs(x) > 1.0)
 65    transitions = [candidates[i] for i in range(1, len(candidates))
 66                   if classified[i] != classified[i - 1]]
 67    observed = float(transitions[0]) if transitions else float("nan")
 68    return {"predicted_boundary": 1.0, "observed_grid_boundary": observed,
 69            "grid_step": float(candidates[1] - candidates[0])}
 70
 71
 72def intervention_experiment(seed=2024, trials=120):
 73    # A slowly worsening plant: a ramps from .88 to 1.035. The monitor acts by
 74    # reducing effective gain. Both policies see exactly the same noise per trial.
 75    rng = np.random.default_rng(seed)
 76    horizon = 320
 77    window = 40
 78    ac_threshold = 0.90
 79    gamma = 1.8
 80    base_rms, safe_rms, base_danger, safe_danger, leads = [], [], [], [], []
 81    for _ in range(trials):
 82        noise = rng.normal(size=horizon) * 0.25
 83        a_schedule = np.linspace(0.90, 1.12, horizon)
 84        xb = xs = 0.0
 85        monitor = CSDMonitor(window, ac_threshold, gamma)
 86        xs_hist = []
 87        triggered_at = None
 88        for t in range(horizon - 1):
 89            xb = a_schedule[t] * xb + noise[t]
 90            xs_hist.append(xs)
 91            mult, trigger = monitor.update(xs)
 92            # Damping is applied to the unstable plant coefficient, clipped to
 93            # preserve a positive recovery coefficient.
 94            effective_a = max(0.0, a_schedule[t] * mult)
 95            xs = effective_a * xs + noise[t]
 96            if trigger and triggered_at is None:
 97                triggered_at = t
 98        base_rms.append(float(np.sqrt(np.mean(np.asarray([0.0] + xs_hist) ** 2))))
 99        safe_rms.append(float(np.sqrt(np.mean(np.asarray(xs_hist) ** 2))))
100        base_danger.append(float(np.max(np.abs(np.asarray([0.0] + xs_hist))) > 4.0))
101        safe_danger.append(float(np.max(np.abs(np.asarray(xs_hist))) > 4.0))
102        if triggered_at is not None:
103            # independently measured crossing is first index with a >= 1
104            crossing = int(np.argmax(a_schedule >= 1.0))
105            leads.append(crossing - triggered_at)
106    return {
107        "trials": trials, "baseline_rms": float(np.mean(base_rms)),
108        "monitor_rms": float(np.mean(safe_rms)),
109        "baseline_danger_rate": float(np.mean(base_danger)),
110        "monitor_danger_rate": float(np.mean(safe_danger)),
111        "trigger_rate": len(leads) / trials,
112        "mean_warning_lead_steps": float(np.mean(leads)) if leads else None,
113        "median_warning_lead_steps": float(np.median(leads)) if leads else None,
114        "settings": {"window": window, "a_c": ac_threshold, "gamma": gamma}
115    }
116
117
118def main():
119    checks = variance_and_ac_checks()
120    boundary = boundary_check()
121    intervention = intervention_experiment()
122    max_var_error = max(r["variance_relative_error"] for r in checks)
123    max_rho_error = max(r["rho_abs_error"] for r in checks)
124    report = {
125        "variance_ac_checks": checks,
126        "boundary_check": boundary,
127        "intervention": intervention,
128        "pass_criteria": {
129            "max_variance_relative_error_lt_0.10": max_var_error < 0.10,
130            "max_rho_abs_error_lt_0.02": max_rho_error < 0.02,
131            "boundary_within_grid_step": abs(boundary["observed_grid_boundary"] - 1.0) <= boundary["grid_step"]
132        },
133        "interpretation": "The first two checks test stationary AR(1) formulas; the intervention is a small illustrative safety simulation, not a trained controller benchmark."
134    }
135    Path("csd_results.json").write_text(json.dumps(report, indent=2))
136    print(json.dumps(report, indent=2))
137
138
139if __name__ == "__main__":
140    main()