1import json, math, random
 2import numpy as np
 3import torch
 4from torch import nn
 5
 6SEED = 184
 7random.seed(SEED); np.random.seed(SEED); torch.manual_seed(SEED)
 8
 9def denominators(J, cutoff=None):
10    # Golden-ratio conjugate: all continued-fraction partial quotients are 1.
11    qs=[]; qm1, q=0, 1
12    for _ in range(J):
13        if cutoff is not None and q > cutoff: break
14        qs.append(q)
15        qn = q + qm1
16        qm1, q = q, qn
17    return np.asarray(qs, dtype=np.int64)
18
19def lacunary_features(x, qs, square_root=True):
20    x = np.asarray(x, dtype=np.float64).reshape(-1, 1)
21    q = np.asarray(qs, dtype=np.float64).reshape(1, -1)
22    scale = 1/np.sqrt(q) if square_root else 1/q
23    z = 2*np.pi*x*q
24    return np.concatenate([np.cos(z)*scale, np.sin(z)*scale], axis=1)
25
26def h_signal(x, qs):
27    q=np.asarray(qs,dtype=float)
28    return (np.cos(2*np.pi*np.asarray(x).reshape(-1,1)*q)/q).sum(1)
29
30def regularity_check():
31    # A direct finite approximation: dyadic modulus, and growth of the
32    # finite-truncation Lipschitz upper bound sum_j 2*pi (each term contributes 2*pi).
33    qs=denominators(25)
34    grid=np.linspace(0,1,200001)
35    vals=h_signal(grid,qs)
36    deltas=2.0**np.arange(-12,-3,dtype=float)
37    mod=[]
38    for d in deltas:
39        k=max(1, int(round(d*len(grid))))
40        mod.append(float(np.max(np.abs(vals[k:]-vals[:-k]))))
41    mod=np.asarray(mod)
42    slope=float(np.polyfit(np.log(deltas), np.log(mod), 1)[0])
43    # For the truncated series, derivative absolute sum is exactly 2*pi*J.
44    growth=[]
45    for J in [4,8,12,16,20,25]: growth.append((J, 2*np.pi*J))
46    # Verify recursion and the claimed absolutely summable coefficient series.
47    recursion_ok=bool(np.all(qs[2:] == qs[1:-1]+qs[:-2]))
48    return {'q_first_12': qs[:12].tolist(), 'recursion_ok':recursion_ok,
49            'sum_inverse_q':float(np.sum(1/qs)), 'dyadic_deltas':deltas.tolist(),
50            'modulus':mod.tolist(), 'loglog_slope':slope,
51            'finite_lipschitz_bound_2piJ':growth}
52
53class MLP(nn.Module):
54    def __init__(self,d):
55        super().__init__(); self.net=nn.Sequential(nn.Linear(d,64),nn.Tanh(),nn.Linear(64,64),nn.Tanh(),nn.Linear(64,1))
56    def forward(self,x): return self.net(x)
57
58def train_model(kind, device):
59    # Target is exactly a longer continued-fraction series; held-out region tests
60    # whether the structured features extrapolate rather than merely interpolate.
61    qtarget=denominators(10)
62    xtr=np.linspace(0, .70, 448, endpoint=False)
63    xte=np.linspace(.70, 1.0, 192, endpoint=False)
64    ytr=h_signal(xtr,qtarget); yte=h_signal(xte,qtarget)
65    qfeat=denominators(8) # [1,1,2,3,5,8,13,21], 16 channels for all encoded models
66    def make(x):
67        if kind=='raw': return x[:,None]
68        if kind=='nerf':
69            # 8 standard dyadic frequencies, same 16 sinusoidal channels.
70            z=2*np.pi*x[:,None]*(2.0**np.arange(8)[None,:])
71            return np.concatenate([np.sin(z),np.cos(z)],1)
72        return np.concatenate([x[:,None],lacunary_features(x,qfeat,True)],1)
73    Xtr=torch.tensor(make(xtr),dtype=torch.float32,device=device); Ytr=torch.tensor(ytr[:,None],dtype=torch.float32,device=device)
74    Xte=torch.tensor(make(xte),dtype=torch.float32,device=device); Yte=torch.tensor(yte[:,None],dtype=torch.float32,device=device)
75    torch.manual_seed(SEED); model=MLP(Xtr.shape[1]).to(device); opt=torch.optim.Adam(model.parameters(),lr=2e-3)
76    for step in range(1800):
77        opt.zero_grad(); loss=((model(Xtr)-Ytr)**2).mean(); loss.backward(); opt.step()
78    with torch.no_grad():
79        tr=float(((model(Xtr)-Ytr)**2).mean().sqrt().cpu()); te=float(((model(Xte)-Yte)**2).mean().sqrt().cpu())
80    return {'train_rmse':tr,'extrapolation_rmse':te,'input_dim':int(Xtr.shape[1]),'parameters':sum(p.numel() for p in model.parameters())}
81
82def main():
83    try:
84        device='cuda' if torch.cuda.is_available() else 'cpu'
85        # CUDA can fail on a shared device; this makes the experiment robust.
86        out={'device':device,'math':regularity_check()}
87        try: out['benchmark']={k:train_model(k,device) for k in ['raw','nerf','lacunary']}
88        except Exception as e:
89            device='cpu'; out['device_fallback']=str(e); out['benchmark']={k:train_model(k,device) for k in ['raw','nerf','lacunary']}
90    except Exception as e:
91        raise
92    with open('results.json','w') as f: json.dump(out,f,indent=2)
93    print(json.dumps(out,indent=2))
94if __name__=='__main__': main()