1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4from scipy.optimize import linprog
  5
  6SEED = 2189
  7np.random.seed(SEED)
  8random.seed(SEED)
  9
 10
 11def certificate(q, S, epsilon, wmax=100.0):
 12    """Maximize delta subject to w_i >= delta and ||q+S w||inf <= epsilon."""
 13    q = np.asarray(q, dtype=float)
 14    S = np.asarray(S, dtype=float)
 15    d, L = len(q), S.shape[1]
 16    A, b = [], []
 17    for i in range(L):
 18        row = np.zeros(L + 1); row[i] = -1; row[-1] = 1
 19        A.append(row); b.append(0.)
 20    for k in range(d):
 21        row = np.zeros(L + 1); row[:L] = S[k]
 22        A.append(row); b.append(epsilon - q[k])
 23        row = np.zeros(L + 1); row[:L] = -S[k]
 24        A.append(row); b.append(epsilon + q[k])
 25    result = linprog(np.r_[np.zeros(L), -1.], A_ub=np.array(A), b_ub=np.array(b),
 26                     bounds=[(0., wmax)] * L + [(0., wmax)], method="highs")
 27    if not result.success:
 28        return np.zeros(L), 0., float("inf"), False
 29    w = result.x[:L]
 30    return w, float(result.x[-1]), float(np.max(np.abs(q + S @ w))), True
 31
 32
 33def mechanism_checks():
 34    # Prediction 1: for S=lambda I, q=(-3,-1), exact weights are (3/lambda,1/lambda),
 35    # hence delta=1/lambda until the wmax ceiling is active.
 36    lambdas = np.array([.25, .5, 1., 2., 4.])
 37    rows1 = []
 38    for lam in lambdas:
 39        w, delta, residual, ok = certificate([-3., -1.], lam*np.eye(2), 1e-9, 100.)
 40        rows1.append(dict(lambda_=float(lam), predicted_delta=float(1/lam),
 41                          observed_delta=delta, weights=w.tolist(), residual=residual))
 42
 43    # Prediction 2: feasibility transition for scalar q=-3, S=1: delta can be positive
 44    # iff epsilon >= 3-wmax; for wmax=2 this predicts epsilon*=1.
 45    epsilons = np.array([.25, .75, 1., 1.25, 2.])
 46    rows2 = []
 47    for eps in epsilons:
 48        w, delta, residual, ok = certificate([-3.], [[1.]], eps, 2.)
 49        rows2.append(dict(epsilon=float(eps), predicted_positive=bool(eps >= 1.),
 50                          observed_positive=bool(delta > 1e-8), delta=delta, weight=w.tolist()))
 51
 52    # Prediction 3: increasing slack epsilon permits a larger common margin for
 53    # q=(-2,-1), S=I; delta=1+eps (the upper residual bound on the second coordinate) in this unconstrained-wmax sweep.
 54    eps3 = np.array([0., .25, .5, .75, 1.])
 55    rows3 = []
 56    for eps in eps3:
 57        w, delta, residual, ok = certificate([-2., -1.], np.eye(2), eps, 100.)
 58        predicted = 1. + eps
 59        rows3.append(dict(epsilon=float(eps), predicted_delta=predicted,
 60                          observed_delta=delta, residual=residual))
 61    return {"scaling_lambda": rows1, "feasibility_transition": rows2,
 62            "slack_margin": rows3}
 63
 64
 65def hinge(v):
 66    return np.maximum(v, 0.)
 67
 68
 69def objective(z, target, thresholds, weights):
 70    # Convex 2-D last-layer problem: performance J=1/2||z-target||2,
 71    # ordered constraints V_i=[z_i-threshold_i]_+ (each tier acts on one coordinate).
 72    J = .5 * np.sum((z-target)**2)
 73    V = hinge(z-np.asarray(thresholds))
 74    return float(J + np.dot(weights, V)), float(J), V
 75
 76
 77def train_controller(mode, seed=SEED, steps=160, refresh=8):
 78    rng=np.random.default_rng(seed)
 79    target=np.array([1.0, 1.0, 1.0])
 80    thresholds=np.array([.2, .4, .6])
 81    z=np.zeros(3)
 82    weights=np.ones(3)
 83    history=[]
 84    for t in range(steps):
 85        if mode == "fixed": weights=np.array([100.,10.,1.])
 86        elif mode == "equal": weights=np.ones(3)
 87        elif mode == "certificate" and t % refresh == 0:
 88            # q = grad J and columns are subgradients of hinge tiers at current z.
 89            q=z-target
 90            S=np.eye(3) * (z-thresholds > 0.)
 91            neww, delta, residual, ok=certificate(q, S, epsilon=.05, wmax=100.)
 92            if ok and delta > 1e-5: weights=neww
 93            history.append((t, delta, residual, weights.copy()))
 94        # subgradient descent on weighted objective
 95        grad=z-target + weights*(z-thresholds > 0.)
 96        z -= .08*grad
 97        J=.5*np.sum((z-target)**2); V=hinge(z-thresholds)
 98        history.append((t, J, float(np.sum(V)), float(V[0]), float(np.dot(weights,V))))
 99    return z, history
100
101
102def summarize():
103    checks=mechanism_checks()
104    comparison={}
105    for mode in ["equal", "fixed", "certificate"]:
106        z,h=train_controller(mode)
107        comparison[mode]={"final_z":z.tolist(), "final_J":float(.5*np.sum((z-np.array([1.,1.,1.]))**2)),
108                          "final_total_violation":float(np.sum(hinge(z-np.array([.2,.4,.6])))),
109                          "final_high_priority_violation":float(hinge(z[0]-.2))}
110    # Relative numerical errors for the predicted equations.
111    e1=max(abs(r["observed_delta"]-r["predicted_delta"])/max(r["predicted_delta"],1e-12)
112           for r in checks["scaling_lambda"])
113    e3=max(abs(r["observed_delta"]-r["predicted_delta"]) for r in checks["slack_margin"])
114    transition=all(r["observed_positive"] == r["predicted_positive"] for r in checks["feasibility_transition"])
115    out={"seed":SEED, "mechanism_checks":checks,
116         "check_summary":{"lambda_scaling_max_relative_error":e1,
117                           "transition_all_correct":transition,
118                           "slack_margin_max_absolute_error":e3},
119         "training_comparison":comparison}
120    Path("results.json").write_text(json.dumps(out, indent=2))
121    print(json.dumps(out, indent=2))
122
123if __name__ == "__main__":
124    summarize()