1import json
  2from pathlib import Path
  3import numpy as np
  4from sklearn.linear_model import LogisticRegression
  5
  6SEED = 3048
  7
  8
  9def primary_action(x):
 10    # Aggressive controller: fast near the goal, but saturation makes it brittle
 11    # from boundary states under disturbances.
 12    return (-0.85 * x if abs(x) < 0.85 else 0.75 * np.sign(x))
 13
 14
 15def fallback_action(x):
 16    # Conservative controller: slower, with smaller action excursions.
 17    return np.clip(-0.42 * x, -0.55, 0.55)
 18
 19
 20def rollout(x0, policy, noise_std, horizon=12, rng=None):
 21    rng = np.random.default_rng() if rng is None else rng
 22    x = float(x0)
 23    reached = False
 24    failed = False
 25    for _ in range(horizon):
 26        x = x + float(policy(x)) + rng.normal(0.0, noise_std)
 27        if abs(x) >= 2.0:
 28            failed = True
 29            break
 30        if abs(x) <= 0.16:
 31            reached = True
 32            break
 33    return reached and not failed, failed, x
 34
 35
 36def features(xs):
 37    xs = np.asarray(xs)
 38    return np.column_stack([np.abs(xs), xs * xs])
 39
 40
 41def make_critic_data(n, noise, policy, seed):
 42    rng = np.random.default_rng(seed)
 43    xs = rng.uniform(-1.95, 1.95, size=n)
 44    successes = np.array([rollout(x, policy, noise, rng=rng)[0] for x in xs], dtype=int)
 45    return features(xs), successes
 46
 47
 48def train_critics():
 49    # P estimates primary-policy success probability.
 50    xp, yp = make_critic_data(8000, 0.42, primary_action, SEED)
 51    # V is failure risk under the conservative fallback, so lower is safer.
 52    xv, yv_success = make_critic_data(8000, 0.62, fallback_action, SEED + 1)
 53    yv_failure = 1 - yv_success
 54    p = LogisticRegression(C=10.0, solver="lbfgs", random_state=SEED).fit(xp, yp)
 55    v = LogisticRegression(C=10.0, solver="lbfgs", random_state=SEED + 1).fit(xv, yv_failure)
 56    return p, v
 57
 58
 59def gate(p_hat, v_hat, eta_p=0.55, eta_r=0.20):
 60    # Execute primary iff both paper conditions pass.
 61    return (p_hat >= eta_p) & (v_hat <= eta_r)
 62
 63
 64def evaluate(mode, p, v, n=5000, noise=0.48, seed=SEED + 20):
 65    rng = np.random.default_rng(seed)
 66    xs = np.clip(rng.normal(0.0, 1.45, n), -1.95, 1.95)
 67    p_hat = p.predict_proba(features(xs))[:, 1]
 68    v_hat = v.predict_proba(features(xs))[:, 1]
 69    if mode == "primary":
 70        use_primary = np.ones(n, dtype=bool)
 71    elif mode == "p_gate":
 72        use_primary = p_hat >= 0.55
 73    elif mode == "dual_gate":
 74        use_primary = gate(p_hat, v_hat)
 75    else:
 76        raise ValueError(mode)
 77    success = np.zeros(n, dtype=bool)
 78    failed = np.zeros(n, dtype=bool)
 79    for i, x in enumerate(xs):
 80        pol = primary_action if use_primary[i] else fallback_action
 81        success[i], failed[i], _ = rollout(x, pol, noise, rng=rng)
 82    return {
 83        "catastrophic_failure_rate": float(failed.mean()),
 84        "goal_success_rate": float(success.mean()),
 85        "primary_coverage": float(use_primary.mean()),
 86        "mean_p_selected": float(p_hat[use_primary].mean()) if use_primary.any() else 0.0,
 87        "mean_v_selected": float(v_hat[use_primary].mean()) if use_primary.any() else 0.0,
 88    }
 89
 90
 91def verify_math():
 92    rng = np.random.default_rng(SEED)
 93    p, v = rng.random(10000), rng.random(10000)
 94    a = gate(p, v, 0.55, 0.20)
 95    direct = (p >= 0.55) & (v <= 0.20)
 96    tighter = gate(p, v, 0.70, 0.10)
 97    return {"formula_exact": bool(np.array_equal(a, direct)),
 98            "stricter_threshold_subset": bool(np.all(~tighter | a)),
 99            "selected_fraction": float(a.mean())}
100
101
102def main():
103    math = verify_math()
104    p, v = train_critics()
105    results = {m: evaluate(m, p, v) for m in ("primary", "p_gate", "dual_gate")}
106    out = {"math_check": math, "results": results,
107           "config": {"seed": SEED, "train_noise_p": 0.42, "train_noise_v": 0.62,
108                      "eval_noise": 0.48, "eta_p": 0.55, "eta_r": 0.20}}
109    Path("results.json").write_text(json.dumps(out, indent=2))
110    print(json.dumps(out, indent=2))
111
112
113if __name__ == "__main__":
114    main()