1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2346
  6
  7
  8def projected_update(lam, vbar, alpha, tau, lam_max):
  9    return float(np.clip(lam + alpha * (vbar - tau), 0.0, lam_max))
 10
 11
 12def policy_target(lam, c=0.8, temperature=0.15):
 13    z = np.clip((c - lam) / temperature, -60.0, 60.0)
 14    return 1.0 / (1.0 + np.exp(-z))
 15
 16
 17def run(alpha, steps=400, beta=0.2, tau=0.2, lam_max=5.0,
 18        batch=64, noisy=False, fixed_penalty=None, seed=SEED):
 19    """Toy policy/dual loop. q is terminal violation probability.
 20    The policy relaxes toward the reward-minus-penalty optimum q*(lambda).
 21    """
 22    rng = np.random.default_rng(seed)
 23    q, lam = 0.8, 0.0
 24    qs, vs, ls, drifts = [], [], [], []
 25    for _ in range(steps):
 26        effective = lam if fixed_penalty is None else fixed_penalty
 27        q += beta * (policy_target(effective) - q)
 28        vbar = float(rng.binomial(batch, np.clip(q, 0.0, 1.0)) / batch) if noisy else q
 29        old = lam
 30        if fixed_penalty is None:
 31            lam = projected_update(lam, vbar, alpha, tau, lam_max)
 32        qs.append(q); vs.append(vbar); ls.append(lam); drifts.append(lam - old)
 33    return {"q": np.asarray(qs), "v": np.asarray(vs),
 34            "lam": np.asarray(ls), "drift": np.asarray(drifts)}
 35
 36
 37def math_check():
 38    r = np.random.default_rng(SEED)
 39    max_err = 0.0
 40    for _ in range(10000):
 41        lam, v, a, tau, cap = r.uniform(0, 10), r.uniform(0, 2), r.uniform(0, 3), r.uniform(0, 1), r.uniform(.1, 10)
 42        expected = min(cap, max(0.0, lam + a * (v - tau)))
 43        max_err = max(max_err, abs(projected_update(lam, v, a, tau, cap) - expected))
 44    x = run(alpha=.08, steps=1500, beta=.2, tau=.2, lam_max=5, noisy=False)
 45    interior = (x["lam"] > 1e-5) & (x["lam"] < 4.99999)
 46    drift_error = abs(np.mean(x["drift"][interior]) - .08 * np.mean((x["v"] - .2)[interior]))
 47    return {"max_projection_formula_error": max_err,
 48            "interior_drift_identity_error": float(drift_error),
 49            "final_lambda": float(x["lam"][-1]), "final_violation": float(x["v"][-1])}
 50
 51
 52def sweep():
 53    out = {}
 54    # Prediction 1: with fixed positive violation error, lambda reaches cap in
 55    # approximately (cap-lambda0)/(alpha*(v-tau)) steps.
 56    cap_rows = []
 57    for alpha in [.02, .05, .1, .2]:
 58        x = run(alpha, steps=300, beta=0.0, tau=0.1, lam_max=2.0, noisy=False)
 59        # beta=0 leaves q at .8, so v-tau=.7 exactly.
 60        hit = np.flatnonzero(x["lam"] >= 2.0 - 1e-10)
 61        observed = int(hit[0] + 1) if len(hit) else None
 62        predicted = int(np.ceil(2.0 / (alpha * .7)))
 63        cap_rows.append({"alpha": alpha, "predicted_hit_step": predicted,
 64                         "observed_hit_step": observed})
 65    out["cap_scaling"] = cap_rows
 66
 67    # Prediction 2: away from boundaries, average violation error is zero in
 68    # a bounded drift-balanced regime.
 69    balance_rows = []
 70    for alpha in [.01, .05, .1, .2, .4]:
 71        x = run(alpha, steps=2000, beta=.2, tau=.2, lam_max=5, noisy=False)
 72        tail = slice(1000, None)
 73        balance_rows.append({"alpha": alpha,
 74                             "tail_mean_violation": float(np.mean(x["v"][tail])),
 75                             "tail_abs_mean_error": float(abs(np.mean(x["v"][tail]) - .2)),
 76                             "tail_lambda_mean": float(np.mean(x["lam"][tail])),
 77                             "tail_lambda_std": float(np.std(x["lam"][tail]))})
 78    out["drift_balance"] = balance_rows
 79
 80    # Prediction 3: positive persistent drift saturates; negative drift projects
 81    # to zero. This tests both projection boundaries.
 82    boundary_rows = []
 83    for v in [0.0, .2, .8]:
 84        x = run(.1, steps=100, beta=0.0, tau=.2, lam_max=1.0, noisy=False)
 85        # beta=0 gives v=.8; use direct recurrence for arbitrary constant v.
 86        lam = 0.; hist=[]
 87        for _ in range(100):
 88            lam = projected_update(lam, v, .1, .2, 1.)
 89            hist.append(lam)
 90        boundary_rows.append({"constant_violation": v, "final_lambda": hist[-1],
 91                              "expected": 1.0 if v > .2 else (0.0 if v < .2 else 0.0)})
 92    out["projection_boundaries"] = boundary_rows
 93
 94    # Secondary baseline: best fixed penalty chosen from a small grid, versus
 95    # adaptive dual controller at equal iterations.
 96    adaptive = run(.08, steps=400, beta=.2, tau=.2, lam_max=5, noisy=True, batch=64)
 97    fixed = []
 98    for p in np.linspace(0, 2, 21):
 99        z = run(0, steps=400, beta=.2, tau=.2, lam_max=5, fixed_penalty=float(p), noisy=True, batch=64, seed=SEED)
100        fixed.append((abs(np.mean(z["v"][-100:]) - .2), p, np.mean(z["v"][-100:])))
101    best = min(fixed)
102    out["comparison"] = {
103        "adaptive_tail_violation": float(np.mean(adaptive["v"][-100:])),
104        "adaptive_tail_abs_error": float(abs(np.mean(adaptive["v"][-100:]) - .2)),
105        "adaptive_tail_lambda": float(np.mean(adaptive["lam"][-100:])),
106        "best_fixed_penalty": best[1], "best_fixed_tail_violation": best[2],
107        "best_fixed_tail_abs_error": best[0]
108    }
109    return out
110
111
112def main():
113    result = {"seed": SEED, "math_check": math_check(), "sweeps": sweep()}
114    Path("results.json").write_text(json.dumps(result, indent=2))
115    print(json.dumps(result, indent=2))
116
117
118if __name__ == "__main__":
119    main()