1import json, math
 2import numpy as np
 3from scipy.special import erfcx
 4from scipy.optimize import nnls
 5
 6
 7def ml_half(t, tau):
 8    return erfcx(np.sqrt(np.asarray(t, dtype=float) / tau))
 9
10
11def bank(lags, taus, weights):
12    rho = np.exp(-1 / np.asarray(taus))
13    return ((weights * (1-rho))[:, None] * rho[:, None]**np.asarray(lags)[None, :]).sum(0)
14
15
16def fit(K, target_tau=32., horizon=16000):
17    lags = np.arange(horizon + 1)
18    target = ml_half(lags + 1, target_tau)
19    # The causal filter is normalized, and each exponential component has unit
20    # infinite-horizon mass, so fit the same unit-mass target.
21    target = target / target.sum()
22    # Fit with simplex weights to preserve unit DC mass.
23    taus = np.geomspace(1, horizon, K)
24    A = ((1-np.exp(-1/taus))[:, None] * np.exp(-1/taus)[:, None]**lags[None, :]).T
25    # Relative error across logarithmic lags, with mild absolute stabilization.
26    ix = np.unique(np.r_[np.arange(100), np.geomspace(100, horizon, 500).astype(int)])
27    scale = 1 / np.maximum(target[ix], 1e-8)
28    w, _ = nnls(A[ix] * scale[:, None], target[ix] * scale)
29    w /= w.sum()
30    pred = bank(lags, taus, w)
31    rel = np.linalg.norm(pred[ix]-target[ix]) / np.linalg.norm(target[ix])
32    return taus, w, float(rel)
33
34
35def normalized(x, taus, w):
36    rho = np.exp(-1/np.asarray(taus)); s=np.zeros(len(rho)); out=[]; mass=0
37    for v in x:
38        s = rho*s + (1-rho)*v
39        mass += np.dot(w, 1-rho)
40        out.append(np.dot(w,s)/mass)
41    return np.asarray(out)
42
43
44def main():
45    rng=np.random.default_rng(7); out={}
46    # Prediction 1: rho^m decay and half-life ln(2)/(-ln rho).
47    rows=[]
48    for rho in [.8,.95,.99,.999]:
49        impulse=(1-rho)*rho**np.arange(20000)
50        obs=int(np.ceil(math.log(.5)/math.log(rho)))
51        rows.append({'rho':rho,'predicted_half_life':math.log(.5)/math.log(rho),
52                     'observed_half_life':obs,'max_impulse':float(impulse.max()),
53                     'stable_max_state':float(np.max(np.cumsum(impulse)))})
54    out['prediction_stability_half_life']=rows
55    # Prediction 2: Mittag-Leffler alpha=.5 has h(l)*sqrt(l) -> 1/sqrt(pi).
56    rows=[]; constant=1/math.sqrt(math.pi)
57    for l in [100,1000,10000,100000,1000000]:
58        scaled=float(ml_half(l,1)*math.sqrt(l))
59        rows.append({'lag':l,'observed':scaled,'predicted':constant,
60                     'relative_error':abs(scaled-constant)/constant})
61    out['prediction_powerlaw']=rows
62    # Prediction 3: increasing logarithmic poles improves fit, at least in the useful range.
63    fits=[]
64    for k in [1,2,4,8,16,32]:
65        taus,w,e=fit(k); fits.append({'K':k,'relative_fit_error':e})
66    out['prediction_bank_scaling']=fits
67    # Secondary controlled signal: fractional bank versus best single exponential.
68    T=16000; x=rng.normal(size=T); h=ml_half(np.arange(T)+1,32.)
69    exact=np.convolve(x,h)[:T]/np.cumsum(h)
70    taus,w,e=fit(8); frac=normalized(x,taus,w)
71    candidates=np.geomspace(1,T,250)
72    singles=[np.mean((normalized(x,[q],[1.])-exact)**2) for q in candidates]
73    frac_mse=float(np.mean((frac-exact)**2)); best=float(min(singles))
74    out['secondary_filter_comparison']={'fractional_K8_mse':frac_mse,
75      'best_single_exponential_mse':best,'single_over_fractional':best/frac_mse}
76    # Delayed retention: compare tails after a pulse at delay from the start.
77    D=4000; pulse=np.zeros(T); pulse[0]=1
78    fr=normalized(pulse,taus,w); one=normalized(pulse,[32],[1])
79    out['secondary_delayed_impulse']={'delay':D,'fractional_at_delay':float(fr[D]),
80      'single_at_delay':float(one[D]),'ratio':float(fr[D]/one[D])}
81    with open('results.json','w') as f: json.dump(out,f,indent=2)
82    print(json.dumps(out,indent=2))
83
84if __name__=='__main__': main()