1import math, json
 2import numpy as np
 3
 4# E_{1/2}(-x): exact series for small x and erfcx asymptotic for large x.
 5def mlhalf_x(x):
 6    x=float(x)
 7    if x<3.5:
 8        term=total=1.0
 9        for n in range(1,500):
10            term *= -x*math.gamma(.5*(n-1)+1)/math.gamma(.5*n+1)
11            total += term
12            if abs(term)<2e-15*max(1,abs(total)): break
13        return max(0.,total)
14    z=1/x**2
15    return (1/math.sqrt(math.pi)/x)*(1-z/2+3*z*z/4-15*z**3/8+105*z**4/16)
16
17def target(l,tau,alpha=.5): return mlhalf_x(((l+1)/tau)**alpha) if alpha==.5 else None
18
19def simplex(v):
20    u=np.sort(v)[::-1]; css=np.cumsum(u)-1; q=np.arange(1,len(v)+1)
21    ok=np.where(u-css/q>0)[0]; j=int(ok[-1]) if len(ok) else 0
22    return np.maximum(v-css[j]/(j+1),0)
23
24def fit(tau=16, taus=(1,4,16,64,256,1024), L=4096):
25    l=np.unique(np.r_[0,np.logspace(0,math.log10(L),320).astype(int)])
26    r=np.exp(-1/np.asarray(taus,float)); A=np.array([[(1-x)*x**z for x in r] for z in l])
27    y=np.array([target(z,tau) for z in l]); w=np.ones(len(r))/len(r)
28    step=1/(np.linalg.norm(A,2)**2+1e-12)
29    for _ in range(10000):
30        nw=simplex(w-step*1.8*A.T@(A@w-y))
31        if np.max(abs(nw-w))<1e-11: w=nw; break
32        w=nw
33    return l,r,w,y,A@w
34
35def sparse_retention(gap, taus, w):
36    r=np.exp(-1/np.asarray(taus,float)); v=np.zeros(len(r))
37    v=(1-r)*1 # impulse at t=0
38    for _ in range(gap): v=r*v
39    return float(np.dot(w,v))
40
41def run_opt(kind, seed=7, steps=500):
42    rng=np.random.default_rng(seed); d=20; x=rng.normal(size=d); Q=np.linspace(1,30,d)
43    m=np.zeros(d); taus=np.array([1,4,16,64,256,1024.]); r=np.exp(-1/taus); w=np.ones(6)/6
44    vv=np.zeros(d) if kind=='adam' else np.zeros((6,d))
45    if kind=='idea':
46        _,r,w,_,_=fit(16,tuple(taus.astype(int)))
47    for t in range(1,steps+1):
48        g=Q*x + .03*rng.normal(size=d)
49        m=.9*m+.1*g
50        if kind=='adam':
51            vv=.999*vv+.001*g*g; den=np.sqrt(vv/(1-.999**t))+1e-8
52        else:
53            vv=r[:,None]*vv+(1-r[:,None])*g[None,:]*g[None,:]
54            den=np.sqrt(sum(w[k]*vv[k]/(1-r[k]**t) for k in range(len(r))))+1e-8
55        x-=.03*m/(1-.9**t)/den
56    return float(.5*np.sum(Q*x*x))
57
58def main():
59    taus=np.array([1,4,16,64,256,1024.]); l,r,w,y,hat=fit()
60    # Prediction 1: asymptotic ML slope is -alpha=-.5.
61    ll=np.arange(512,4096); slope=np.polyfit(np.log(ll),np.log([target(z,16) for z in ll]),1)[0]
62    # Prediction 2: exponential bank has a less negative late slope than Adam beta=.999.
63    late=np.arange(1000,4000); bank=np.array([sum(w[k]*(1-r[k])*r[k]**z for k in range(6)) for z in late])
64    ok=bank>1e-300; bs=np.polyfit(np.log(late[ok]),np.log(bank[ok]),1)[0]
65    adam_slope=math.log(.999)
66    # Prediction 3: old impulse retention grows with longest available memory.
67    gaps=[1,16,128,512]; ret=[sparse_retention(g,taus,w) for g in gaps]
68    # Compare equal toy steps.
69    base=[run_opt('adam',s) for s in range(3)]; idea=[run_opt('idea',s) for s in range(3)]
70    out={'alpha':.5,'fit_weights':w.tolist(),'kernel_rel_rmse':float(np.linalg.norm(hat-y)/np.linalg.norm(y)),
71      'predictions':{'tail_slope_pred':-.5,'tail_slope_obs':float(slope),'tol':.08,
72      'late_bank_log_slope_obs':float(bs),'single_ema_log_slope_obs':adam_slope,
73      'sparse_gaps':gaps,'sparse_retention':ret,'retention_monotone':bool(np.all(np.diff(ret)<0))},
74      'toy_final_losses':{'adamw':base,'mittag_leffler':idea},
75      'math_normalization_note':'The stated h_hat uses (1-rho)rho^l while sum_l h_hat=1, whereas h_0 is generally <1; fitting with sum(w)=1 cannot exactly match h_0.'}
76    print(json.dumps(out,indent=2))
77if __name__=='__main__': main()