1import json
  2import time
  3import numpy as np
  4from scipy.optimize import minimize
  5
  6SEED = 7
  7
  8def dirs(theta):
  9    return np.stack((np.cos(theta), np.sin(theta)), axis=1)
 10
 11def Bmat(x, theta, k=6.0):
 12    return np.exp(1j * k * (x @ dirs(theta).T))
 13
 14def ridge(B, y, lam=0.0, w=None):
 15    if w is None:
 16        w = np.ones(len(y))
 17    A = B.conj().T @ (w[:, None] * B) + lam * np.eye(B.shape[1])
 18    b = B.conj().T @ (w * y)
 19    c = np.linalg.solve(A, b)
 20    r = B @ c - y
 21    phi = float(np.sum(w * np.abs(r)**2) + lam * np.vdot(c, c).real)
 22    return c, phi, r
 23
 24def vp_obj(theta, x, y, k, lam):
 25    return ridge(Bmat(x, theta, k), y, lam)[1]
 26
 27def vp_fit(x, y, k, lam, init, maxiter=220):
 28    result = minimize(lambda t: vp_obj(t, x, y, k, lam), init,
 29                      method='Nelder-Mead',
 30                      options={'maxiter': maxiter, 'xatol': 2e-7, 'fatol': 1e-11})
 31    B = Bmat(x, result.x, k)
 32    c, phi, r = ridge(B, y, lam)
 33    return result, c, phi, r
 34
 35def relative_error(B, c, y):
 36    return float(np.linalg.norm(B @ c-y) / np.linalg.norm(y))
 37
 38def checks():
 39    rng=np.random.default_rng(SEED)
 40    n=180; k=6.; x=rng.uniform(-1,1,(n,2)); theta=np.array([.1,1.,2.2])
 41    ctrue=np.array([1+.2j,-.7+.4j,.4-.3j])
 42    y=Bmat(x,theta,k)@ctrue + .03*(rng.normal(size=n)+1j*rng.normal(size=n))
 43    ls=np.array([0.,1e-4,1e-3,1e-2,.1,1.])
 44    ph=[]; cn=[]
 45    for lam in ls:
 46        c,p,_=ridge(Bmat(x,theta,k),y,lam); ph.append(p); cn.append(np.vdot(c,c).real)
 47    ph=np.array(ph); cn=np.array(cn)
 48    # Envelope prediction: ridge objective nondecreasing; coefficient norm nonincreasing.
 49    ridge_ok=bool(np.all(np.diff(ph)>=-1e-8) and np.all(np.diff(cn)<=1e-8))
 50
 51    # Gramian prediction: for two nearly duplicate directions, smallest eigenvalue
 52    # scales quadratically in angular separation (Taylor expansion of columns).
 53    base=.65; deltas=np.array([.01,.02,.04,.08,.16]); eig=[]
 54    for d in deltas:
 55        G=Bmat(x,np.array([base,base+d]),k).conj().T @ Bmat(x,np.array([base,base+d]),k)
 56        eig.append(np.linalg.eigvalsh(G)[0]/n)
 57    eig=np.array(eig)
 58    slope=float(np.polyfit(np.log(deltas),np.log(np.maximum(eig,1e-16)),1)[0])
 59    gram_ok=bool(1.5 < slope < 2.5)
 60
 61    # Direction identifiability prediction: increasing sample count improves recovery.
 62    recovery=[]
 63    planted=np.array([.35,1.45]); cc=np.array([1.1-.3j,.65+.5j])
 64    for nn in [30,60,120,240]:
 65        xx=rng.uniform(-1,1,(nn,2)); yy=Bmat(xx,planted,k)@cc
 66        init=np.array([2.7,-1.1])
 67        out,_,_,_=vp_fit(xx,yy,k,1e-8,init,180)
 68        # account for 2pi periodicity and column permutation
 69        est=np.sort(np.mod(out.x,2*np.pi)); true=np.sort(np.mod(planted,2*np.pi))
 70        recovery.append(float(np.linalg.norm(np.angle(np.exp(1j*(est-true))))))
 71    # prediction is lower error with more data, allowing optimizer noise.
 72    recover_ok=bool(recovery[-1] < recovery[0] and recovery[-1] < .25)
 73    return {
 74        'ridge_lambdas':ls.tolist(),'ridge_objectives':ph.tolist(),'coefficient_norms':cn.tolist(),
 75        'ridge_monotonic_predicted':True,'ridge_monotonic_observed':ridge_ok,
 76        'duplicate_deltas':deltas.tolist(),'duplicate_small_eigenvalues':eig.tolist(),
 77        'gramian_scaling_predicted_exponent':2.0,'gramian_scaling_observed_exponent':slope,
 78        'gramian_quadratic_predicted':True,'gramian_quadratic_observed':gram_ok,
 79        'sample_counts':[30,60,120,240],'direction_errors':recovery,
 80        'recovery_improves_predicted':True,'recovery_improves_observed':recover_ok,
 81        'all_checks_pass':ridge_ok and gram_ok and recover_ok}
 82
 83def experiment():
 84    rng=np.random.default_rng(SEED+1); k=6.; n=180
 85    x=rng.uniform(-1,1,(n,2)); planted=np.array([.35,1.45,.35+0.65])
 86    # Duplicate-looking third component intentionally tests adaptive fitting under redundancy.
 87    ctrue=np.array([1.0-.25j,.7+.4j,.15-.1j]); y=Bmat(x,planted,k)@ctrue
 88    init=rng.uniform(0,2*np.pi,3)
 89    t=time.perf_counter(); out,c,phi,r=vp_fit(x,y,k,1e-7,init,300); tv=time.perf_counter()-t
 90    idea_err=relative_error(Bmat(x,out.x,k),c,y)
 91    # Standard fixed random Fourier features: same number of columns, coefficients fitted only.
 92    fixed=rng.uniform(0,2*np.pi,3); t=time.perf_counter(); cf,pf,rf=ridge(Bmat(x,fixed,k),y,1e-7); tf=time.perf_counter()-t
 93    base_err=relative_error(Bmat(x,fixed,k),cf,y)
 94    # A fairer fixed baseline averages several random dictionaries.
 95    errs=[]
 96    for _ in range(12):
 97        th=rng.uniform(0,2*np.pi,3); cb,_,_=ridge(Bmat(x,th,k),y,1e-7)
 98        errs.append(relative_error(Bmat(x,th,k),cb,y))
 99    G=Bmat(x,out.x,k).conj().T@Bmat(x,out.x,k)
100    ev=np.linalg.eigvalsh(G); retained=int(np.sum(ev >= 1e-3*ev.max()))
101    return {'n':n,'features':3,'planted_directions':planted.tolist(),
102            'optimized_directions':out.x.tolist(),'idea_relative_error':idea_err,
103            'fixed_random_relative_error_one_draw':base_err,'fixed_random_relative_error_mean':float(np.mean(errs)),
104            'fixed_random_relative_error_std':float(np.std(errs)), 'retained_by_gramian_cutoff':retained,
105            'idea_seconds':tv,'fixed_seconds':tf,'optimizer_success':bool(out.success),'optimizer_iterations':int(out.nit)}
106
107def main():
108    result={'checks':checks(),'experiment':experiment()}
109    with open('results.json','w') as f: json.dump(result,f,indent=2)
110    print(json.dumps(result,indent=2))
111
112if __name__=='__main__': main()