1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2037
  6np.random.seed(SEED)
  7random.seed(SEED)
  8
  9class BoundedPredictiveGain:
 10    """Euclidean mirror descent gain, one scalar per group."""
 11    def __init__(self, alpha0, eta, amin=None, amax=None, rho=0.0, qmax=None):
 12        self.alpha0 = float(alpha0)
 13        self.eta = float(eta)
 14        self.amin = 0.1 * alpha0 if amin is None else float(amin)
 15        self.amax = 10.0 * alpha0 if amax is None else float(amax)
 16        self.rho = float(rho)
 17        self.qmax = qmax
 18        self.a = float(alpha0)
 19        self.ref = float(alpha0)
 20        self.clip_count = 0
 21        self.steps = 0
 22    def step(self, q):
 23        q0 = float(np.clip(q, -self.qmax, self.qmax)) if self.qmax is not None else float(q)
 24        raw = (self.a - self.eta * q0 + self.rho * self.ref) / (1.0 + self.rho)
 25        new = float(np.clip(raw, self.amin, self.amax))
 26        self.clip_count += int(new != raw)
 27        self.a = new
 28        self.steps += 1
 29        return new
 30
 31def toy_prediction_1_sign_and_slope():
 32    # At a=aref the exact update predicts delta a = eta*r*g^2/(1+rho).
 33    eta, rho, g, a0 = 0.07, 0.4, 1.3, 0.8
 34    rs = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
 35    observed = []
 36    predicted = []
 37    for r in rs:
 38        opt = BoundedPredictiveGain(a0, eta, amin=0.0, amax=100.0, rho=rho)
 39        old = opt.a
 40        opt.step(-r * g * g)
 41        observed.append(opt.a - old)
 42        predicted.append(eta * r * g * g / (1 + rho))
 43    slope = float(np.polyfit(rs, observed, 1)[0])
 44    expected_slope = eta * g * g / (1 + rho)
 45    max_err = float(np.max(np.abs(np.array(observed) - predicted)))
 46    return {
 47        "name": "sign and magnitude of consecutive-direction product",
 48        "prediction": {"delta_a": "eta*r*g^2/(1+rho)", "slope": expected_slope},
 49        "observed": {"r": rs.tolist(), "delta_a": observed, "slope": slope, "max_abs_error": max_err},
 50        "pass": bool(max_err < 1e-12 and slope > 0 and observed[0] < 0 and observed[-1] > 0)
 51    }
 52
 53def toy_prediction_2_persistence():
 54    # With q=0, a'=(a+rho*ref)/(1+rho); error contracts by 1/(1+rho).
 55    a0, ref, rho = 2.0, 0.7, 0.5
 56    opt = BoundedPredictiveGain(ref, eta=0.1, amin=0, amax=100, rho=rho)
 57    opt.a = a0
 58    errors = [abs(opt.a - ref)]
 59    for _ in range(5):
 60        opt.step(0.0)
 61        errors.append(abs(opt.a - ref))
 62    ratios = np.array(errors[1:]) / np.array(errors[:-1])
 63    expected = 1 / (1 + rho)
 64    return {
 65        "name": "Bregman persistence contraction when predictive product is zero",
 66        "prediction": {"error_ratio": expected},
 67        "observed": {"errors": errors, "ratios": ratios.tolist(), "max_abs_ratio_error": float(np.max(abs(ratios-expected)))},
 68        "pass": bool(np.max(abs(ratios-expected)) < 1e-12)
 69    }
 70
 71def toy_prediction_3_quadratic_boundary():
 72    # For f=.5*lam*theta^2 and fixed a: theta sign first reverses at a*lam>1,
 73    # and linear dynamics is stable iff 0<a*lam<2.
 74    lam = 1.0
 75    gains = np.linspace(0.2, 2.4, 12)
 76    sign_reversal = []
 77    stable = []
 78    for a in gains:
 79        theta = 1.0
 80        signs = []
 81        for _ in range(80):
 82            signs.append(np.sign(theta))
 83            theta = (1-a*lam)*theta
 84        # Strict reversal means consecutive nonzero states have opposite signs;
 85        # at a*lam=1 the iterate lands exactly at zero rather than overshooting.
 86        vals = np.array(signs)
 87        sign_reversal.append(bool(np.any(vals[1:] * vals[:-1] < 0)))
 88        stable.append(abs(1-a*lam) < 1)
 89    first_observed = float(gains[np.where(sign_reversal)[0][0]])
 90    first_grid_above = float(gains[np.where(gains > 1/lam)[0][0]])
 91    stable_mismatch = int(np.sum(np.array(stable) != ((gains > 0) & (gains < 2/lam))))
 92    return {
 93        "name": "quadratic first sign reversal and stability boundary",
 94        "prediction": {"first_reversal_threshold": 1/lam, "stability_interval": [0.0, 2/lam]},
 95        "observed": {"first_grid_reversal": first_observed, "first_grid_point_above_threshold": first_grid_above, "stable_boundary_test_mismatches": stable_mismatch, "grid": gains.tolist(), "reversal": sign_reversal},
 96        "pass": bool(first_observed > 1/lam and first_observed <= first_grid_above and stable_mismatch == 0)
 97    }
 98
 99def make_data(n=800):
100    rng = np.random.RandomState(SEED)
101    x = rng.randn(n, 2).astype(np.float32)
102    y = ((x[:,0] * x[:,1] + 0.25*x[:,0] - 0.15*x[:,1]) > 0).astype(np.int64)
103    return x, y
104
105def mini_training():
106    # Small deterministic numpy MLP keeps the comparison transparent and cheap.
107    x, y = make_data()
108    rng = np.random.RandomState(SEED+1)
109    d, h = 2, 24
110    w1 = rng.randn(d,h).astype(np.float64)*0.5; b1=np.zeros(h)
111    w2 = rng.randn(h,2).astype(np.float64)*0.5; b2=np.zeros(2)
112    init = [z.copy() for z in (w1,b1,w2,b2)]
113    def run(adaptive):
114        ww = [z.copy() for z in init]; gains=[0.8,0.8]; prev_u=[None,None,None,None]
115        losses=[]; accs=[]; clips=0; reversals=0
116        for t in range(160):
117            ix = np.arange((t*64)% (len(x)-64), (t*64)% (len(x)-64)+64)
118            X=x[ix].astype(np.float64); Y=y[ix]
119            z=X@ww[0]+ww[1]; ah=np.maximum(z,0); logits=ah@ww[2]+ww[3]
120            logits-=logits.max(1,keepdims=True); p=np.exp(logits); p/=p.sum(1,keepdims=True)
121            loss=-np.log(p[np.arange(len(Y)),Y]+1e-12).mean(); losses.append(float(loss))
122            dl=p; dl[np.arange(len(Y)),Y]-=1; dl/=len(Y)
123            grads=[X.T@(dl@ww[2].T*(z>0)), (dl@ww[2].T*(z>0)).sum(0), ah.T@dl, dl.sum(0)]
124            us=[grads[0],grads[1],grads[2],grads[3]]
125            if adaptive and t > 0:
126                for gi, inds in enumerate(([0,1],[2,3])):
127                    dot=sum(float(np.mean(prev_u[j]*us[j])) for j in inds)/2
128                    q=-dot
129                    old=gains[gi]
130                    gains[gi]=float(np.clip((old-0.12*q+0.15*0.8)/(1.15),0.08,8.0))
131                    clips += int(gains[gi] in (0.08,8.0)); reversals += int(dot < 0)
132            step_g = gains if adaptive else [0.8,0.8]
133            for j in [0,1]: ww[j] -= step_g[0]*us[j]
134            for j in [2,3]: ww[j] -= step_g[1]*us[j]
135            prev_u=[u.copy() for u in us]
136            # full-data accuracy checkpoint
137            zz=x@ww[0]+ww[1]; pp=np.maximum(zz,0)@ww[2]+ww[3]
138            accs.append(float((pp.argmax(1)==y).mean()))
139        return {"final_loss":losses[-1],"best_loss":min(losses),"final_accuracy":accs[-1],"loss_at_40":losses[39],"gain_trajectories":gains,"clip_events":clips,"reversal_events":reversals}
140    return run(False), run(True)
141
142def main():
143    checks=[toy_prediction_1_sign_and_slope(),toy_prediction_2_persistence(),toy_prediction_3_quadratic_boundary()]
144    baseline, idea=mini_training()
145    out={"seed":SEED,"checks":checks,"mini_experiment":{"baseline":baseline,"bounded_predictive_gain":idea}}
146    Path("results.json").write_text(json.dumps(out, indent=2))
147    print(json.dumps(out, indent=2))
148
149if __name__ == "__main__": main()