Spectral-gap adaptive polynomial filtering / analyze.py

Failed on benchmark

Raw ⬇ ZIP
 1import json, math
 2import numpy as np
 3from spectral_gap_filter import *
 4
 5def plain_coeffs(K):
 6    c=np.zeros(K+1); c[K]=1.; return c
 7rows=[]
 8for K in [7,15,31,63,127]:
 9  for a in np.linspace(.25,12,48):
10    s=min(.95,a/K); lam,_=rotation_spectrum(s,n=8192)
11    f=fejer_coeffs(K); j=jackson_coeffs(K); b=plain_coeffs(K)
12    rf=residual_ratio(lam,f); rj=residual_ratio(lam,j); rb=residual_ratio(lam,b)
13    accept=(rj <= 1.1*rf)
14    rows.append((K,a,rf,rj,rb,accept))
15# Find first a where Jackson is accepted and remains useful (sampled)
16for K in [15,31,63,127]:
17 q=[x for x in rows if x[0]==K]
18 first=next((x for x in q if x[5]),None)
19 print('K',K,'first_accept_sK',None if first is None else round(first[1],3), 'ratio_at_2', round(next(x[3]/x[2] for x in q if abs(x[1]-2)<.2),3))
20# scaling at fixed a=8 and fixed s for K
21print('scaling_fixed_sK8')
22for K,a,rf,rj,rb,ac in rows:
23 if abs(a-8)<.13: print(K, 's',a/K,'jackson*K2s',rj*K*K*(a/K),'plain',rb,'fejer',rf)
24print('scaling_fixed_s=.25')
25for K in [7,15,31,63,127]:
26 s=.25; lam,_=rotation_spectrum(s,n=8192); r=residual_ratio(lam,jackson_coeffs(K)); print(K,r,r*K*K*s)
27# save compact report
28out={'math_checks':math_checks(),'critical_prediction':{'predicted': 'Fejer residual*d0^{-1} = 1/(K+1)', 'observed_Kplus1_residual':[float(residual_ratio(rotation_spectrum(1.0/K)[0],fejer_coeffs(K))*(K+1)) for K in [7,15,31,63,127]]},
29'gap_scaling_prediction':{'predicted':'Jackson residual proportional to 1/(K^2*s) at fixed sK large','observed_K2s_at_sK8':[(K,float(rj*K*K*(8/K))) for K,a,rf,rj,rb,ac in rows if abs(a-8)<.13]},
30'switch_prediction':{'predicted':'safe below sK=2; safeguard prevents >10% increase','observed':[]}, 'rows':[]}
31for K in [15,31,63,127]:
32 q=[x for x in rows if x[0]==K]
33 for a in [1,2,4,8]:
34  x=min(q,key=lambda z:abs(z[1]-a)); out['switch_prediction']['observed'].append({'K':K,'sK':a,'jackson_over_fejer':float(x[3]/x[2]),'accepted':bool(x[5])})
35json.dump(out,open('report.json','w'),indent=2)