1import json, math, os
  2import numpy as np
  3from scipy.special import erfcx
  4from scipy.optimize import nnls
  5
  6# For alpha=1/2: E_{1/2}(-z) = exp(z^2) erfc(z) = erfcx(z).
  7def ml_half(t, tau=1.0):
  8    z = np.sqrt(np.asarray(t, dtype=float) / tau)
  9    return erfcx(z)
 10
 11def bank_kernel(lags, taus, weights):
 12    lags = np.asarray(lags, dtype=float)
 13    taus = np.asarray(taus, dtype=float)
 14    rho = np.exp(-1.0 / taus)
 15    return ((weights * (1-rho))[:, None] * rho[:, None] ** lags[None, :]).sum(axis=0)
 16
 17def fit_bank(alpha, tau_target, K, max_lag=20000):
 18    # Fit the normalized finite-horizon causal kernel from the definition.
 19    lags = np.unique(np.r_[0, np.geomspace(1, max_lag, 600).astype(int)])
 20    raw = ml_half(lags + 1, tau_target)
 21    # Approximate the prefix normalizer using all integer lags.
 22    allraw = ml_half(np.arange(max_lag + 1) + 1, tau_target)
 23    target = raw / allraw.sum()
 24    taus = np.geomspace(1, max_lag, K)
 25    rho = np.exp(-1.0 / taus)
 26    A = ((1-rho)[:, None] * rho[:, None] ** lags[None, :]).T
 27    # Relative weighting prevents the long tail from being ignored.
 28    scale = 1.0 / np.maximum(target, 1e-12)
 29    try:
 30        w, _ = nnls(A * scale[:, None], target * scale, maxiter=100000)
 31    except RuntimeError:
 32        from scipy.optimize import lsq_linear
 33        w = lsq_linear(A * scale[:, None], target * scale, bounds=(0, np.inf), max_iter=10000).x
 34    w = w / max(w.sum(), 1e-30)
 35    test = np.unique(np.r_[np.arange(0, 1000), np.geomspace(1000, max_lag, 1000).astype(int)])
 36    yt = ml_half(test + 1, tau_target) / allraw.sum()
 37    yp = bank_kernel(test, taus, w)
 38    return dict(K=K, rel_l2=float(np.linalg.norm(yp-yt)/np.linalg.norm(yt)),
 39                max_rel_tail=float(np.max(np.abs(yp[-500:]-yt[-500:])/(yt[-500:]+1e-12))),
 40                taus=taus.tolist(), weights=w.tolist())
 41
 42def normalized_filter(x, taus, weights, normalize=True):
 43    state=np.zeros(len(taus)); out=[]; mass=0.0
 44    rho=np.exp(-1/np.asarray(taus))
 45    for v in x:
 46        state=rho*state+(1-rho)*v
 47        mass += np.dot(weights, 1-rho) if normalize else 0.0
 48        out.append(np.dot(weights,state) / mass if normalize else np.dot(weights,state))
 49    return np.asarray(out)
 50
 51def main():
 52    rng=np.random.default_rng(7)
 53    results={}
 54    # Prediction 1: stability and exact exponential half-life ln(.5)/ln(rho).
 55    stability=[]
 56    for rho in [0.5,0.8,0.95,0.99,0.999]:
 57        s=0.; peak=0.; vals=[]
 58        for t in range(10000):
 59            s=rho*s+(1-rho)*(1.0 if t==0 else 0.0)
 60            peak=max(peak,abs(s)); vals.append(s)
 61        arr=np.asarray(vals); ix=int(np.argmin(np.abs(arr-rho*0.5))) # after one? use first crossing half peak
 62        half=int(np.where(arr<=0.5*peak)[0][0])
 63        pred=math.log(.5)/math.log(rho)
 64        stability.append(dict(rho=rho, predicted_half_life=pred, observed_half_life=half,
 65                              max_abs=float(peak), final_abs=float(abs(arr[-1]))))
 66    results['stability_half_life']=stability
 67    # Prediction 2: fractional tail l^alpha h(l) -> 1/Gamma(1-alpha), alpha=.5.
 68    tail=[]
 69    for l in [100,1000,10000,100000,1000000]:
 70        h=float(ml_half(l+1,1.0)); tail.append(dict(l=l, scaled=float(h*math.sqrt(l)), predicted=1/math.sqrt(math.pi)))
 71    results['power_law_tail']=tail
 72    # Prediction 3: more logarithmic poles reduce approximation error.
 73    fits=[fit_bank(.5, 32., K) for K in [2,4,8,16,32]]
 74    results['bank_scaling']=fits
 75    # Secondary delayed-memory comparison: exact fractional normalized filter versus
 76    # the best one-pole normalized filter, on random sequences and a long horizon.
 77    T=16000; x=rng.normal(size=T); target=normalized_filter(x,[1], [1]) # placeholder overwritten
 78    frac=fit_bank(.5,32.,8); taus=np.asarray(frac['taus']); w=np.asarray(frac['weights'])
 79    y=normalized_filter(x,taus,w)
 80    # exact normalized causal kernel convolution, computed by direct recurrence-free convolution
 81    h=ml_half(np.arange(T)+1,32.)
 82    exact=np.convolve(x,h)[:T]/np.cumsum(h)
 83    one_tau=np.geomspace(1, T, 300)
 84    errs=[]
 85    for tt in one_tau:
 86        z=normalized_filter(x,[tt],[1.])
 87        errs.append(np.mean((z-exact)**2))
 88    best=float(min(errs)); fracerr=float(np.mean((y-exact)**2))
 89    results['filter_comparison']={'fractional_K8_mse':fracerr,'best_single_exponential_mse':best,
 90                                  'improvement_ratio':best/(fracerr+1e-30)}
 91    # delayed copy signal: correlation with a value 4000 steps back, highlighting retention.
 92    D=4000; probe=np.zeros(T); probe[D]=1.; response=normalized_filter(probe,taus,w)
 93    single=normalized_filter(probe,[32.],[1.])
 94    results['delayed_impulse']={'delay':D,'fractional_response':float(response[D]),
 95                                'single_tau32_response':float(single[D]),
 96                                'fractional_tail_at_delay':float(response[-1]),
 97                                'single_tail_at_delay':float(single[-1])}
 98    with open('results.json','w') as f: json.dump(results,f,indent=2)
 99    print(json.dumps({k:v for k,v in results.items() if k not in ('bank_scaling',)}, indent=2))
100    print('bank_scaling', [(r['K'],r['rel_l2'],r['max_rel_tail']) for r in fits])
101
102if __name__=='__main__': main()