1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 913
  6random.seed(SEED)
  7np.random.seed(SEED)
  8
  9# Memory-retaining profile h(x, eta)=eta^alpha F(x/eta^beta).
 10def profile(x, eta, alpha, beta):
 11    u = x / (eta ** beta)
 12    return eta ** alpha * np.exp(-0.5*u*u) * (1.0 + 0.35*np.cos(3.0*u))
 13
 14def effective_width(x, y):
 15    w = np.maximum(y, 0.0)**2
 16    return float(np.sqrt(np.sum(w*x*x)/(np.sum(w)+1e-30)))
 17
 18def log_slope(a, b):
 19    return float(np.polyfit(np.log(a), np.log(b), 1)[0])
 20
 21def math_sweep():
 22    alpha, beta = 0.70, 0.62
 23    etas = np.geomspace(0.35, 3.0, 13)
 24    x = np.linspace(-10, 10, 20001)
 25    peaks = np.array([np.max(profile(x,e,alpha,beta)) for e in etas])
 26    widths = np.array([effective_width(x, profile(x,e,alpha,beta)) for e in etas])
 27    alpha_hat = log_slope(etas, peaks)
 28    beta_hat = log_slope(etas, widths)
 29
 30    # Prediction 1: amplitude slope is alpha.
 31    # Prediction 2: width slope is beta.
 32    # Prediction 3: with x'=x/L and eta'=eta/L^(1/beta),
 33    # h(x',eta') = L^(-alpha/beta) h(x,eta), so normalized residual is zero.
 34    eta = 1.1
 35    exact_residuals, wrong_eta_residuals = [], []
 36    for L in [1.25, 1.5, 2.0, 3.0]:
 37        ep = eta / (L ** (1.0/beta))
 38        lhs = profile(x/L, ep, alpha, beta)
 39        rhs = (L ** (-alpha/beta))*profile(x, eta, alpha, beta)
 40        exact_residuals.append(np.linalg.norm(lhs-rhs)/(np.linalg.norm(rhs)+1e-12))
 41        ep_wrong = eta/L
 42        lhs_wrong = profile(x/L, ep_wrong, alpha, beta)
 43        wrong_eta_residuals.append(np.linalg.norm(lhs_wrong-rhs)/(np.linalg.norm(rhs)+1e-12))
 44
 45    return {
 46        "true_alpha": alpha, "fitted_alpha": alpha_hat,
 47        "true_beta": beta, "fitted_beta": beta_hat,
 48        "amplitude_abs_error": abs(alpha_hat-alpha),
 49        "width_abs_error": abs(beta_hat-beta),
 50        "rg_L": [1.25,1.5,2.0,3.0],
 51        "rg_exact_residual": exact_residuals,
 52        "rg_wrong_eta_residual": wrong_eta_residuals,
 53        "predictions": [
 54            "log peak amplitude slope equals alpha",
 55            "log profile width slope equals beta",
 56            "matched RG rescaling gives zero residual while eta/L does not"
 57        ]
 58    }
 59
 60def mini_experiment():
 61    # A deliberately tiny signal classification task: class is encoded by a
 62    # narrow/wide Gaussian and amplitude. Inputs are randomly shifted profiles.
 63    try:
 64        import torch
 65        import torch.nn as nn
 66        from torch.utils.data import DataLoader, TensorDataset
 67        torch.manual_seed(SEED)
 68        device = "cpu"
 69        # CPU is the reliable fallback when shared CUDA convolution engines fail.
 70    except Exception as e:
 71        return {"error": "torch unavailable: "+str(e)}
 72
 73    ntrain, nval, length = 1600, 500, 64
 74    def make(n, heldout=False):
 75        X, Y, E = [], [], []
 76        for _ in range(n):
 77            y = np.random.randint(0,2)
 78            # Validation uses a shifted scale range, testing extrapolation.
 79            eta = np.random.uniform(0.55,0.95) if (not heldout and y==0) else None
 80            if eta is None:
 81                eta = np.random.uniform(1.05,1.55) if heldout else np.random.uniform(1.05,1.45)
 82            amp = 1.0 if y==0 else 1.55
 83            xx = np.linspace(-3,3,length)
 84            shift = np.random.uniform(-0.25,0.25)
 85            sig = eta*(0.70 if y==0 else 0.95)
 86            z = amp*np.exp(-0.5*((xx-shift)/sig)**2)
 87            z += np.random.normal(0,0.045,length)
 88            X.append(z.astype(np.float32)); Y.append(y); E.append(eta)
 89        return torch.tensor(np.array(X)), torch.tensor(Y), torch.tensor(np.array(E),dtype=torch.float32)
 90    Xtr,Ytr,Etr = make(ntrain,False); Xv,Yv,Ev = make(nval,True)
 91    train = DataLoader(TensorDataset(Xtr,Ytr,Etr), batch_size=64, shuffle=True)
 92
 93    class Baseline(nn.Module):
 94        def __init__(self):
 95            super().__init__(); self.net=nn.Sequential(nn.Conv1d(1,12,5,stride=2,padding=2),nn.ReLU(),nn.Conv1d(12,20,5,stride=2,padding=2),nn.ReLU(),nn.AdaptiveAvgPool1d(1)); self.fc=nn.Linear(20,2)
 96        def forward(self,x,e): return self.fc(self.net(x[:,None,:]).squeeze(-1))
 97    class Retained(nn.Module):
 98        def __init__(self):
 99            super().__init__(); self.conv=nn.Sequential(nn.Conv1d(1,12,5,stride=2,padding=2),nn.ReLU(),nn.Conv1d(12,20,5,stride=2,padding=2),nn.ReLU()); self.film=nn.Linear(1,40); self.fc=nn.Linear(20,2)
100        def forward(self,x,e):
101            h=self.conv(x[:,None,:]); ab=self.film(torch.log(e[:,None]+1e-5)); a,b=ab[:,:20,None],ab[:,20:,None]; h=h*(1+a)+b; return self.fc(h.mean(-1))
102    def train_eval(model):
103        model.to(device); opt=torch.optim.Adam(model.parameters(),lr=2e-3); lossfn=nn.CrossEntropyLoss()
104        for _ in range(10):
105            model.train()
106            for x,y,e in train:
107                x,y,e=x.to(device),y.to(device),e.to(device); opt.zero_grad(); loss=lossfn(model(x,e),y); loss.backward(); opt.step()
108        model.eval(); correct=0; total=0; losses=[]
109        with torch.no_grad():
110            for x,y,e in DataLoader(TensorDataset(Xv,Yv,Ev),batch_size=128):
111                out=model(x.to(device),e.to(device)); losses.append(float(lossfn(out,y.to(device)))); correct += int((out.argmax(1).cpu()==y).sum()); total += len(y)
112        return {"accuracy":correct/total,"loss":float(np.mean(losses))}
113    try:
114        base=train_eval(Baseline()); idea=train_eval(Retained())
115    except Exception as exc:
116        # Retry from scratch on CPU for shared/unsupported CUDA backends.
117        device = "cpu"
118        base=train_eval(Baseline()); idea=train_eval(Retained())
119    return {"device":device,"baseline":base,"idea":idea,"train_size":ntrain,"validation_size":nval}
120
121def main():
122    result={"math":math_sweep(),"mini_experiment":mini_experiment()}
123    Path("results.json").write_text(json.dumps(result,indent=2))
124    print(json.dumps(result,indent=2))
125
126if __name__ == "__main__": main()