1import json, math, os
  2import numpy as np
  3from scipy.linalg import eigvals, expm, svdvals
  4
  5
  6def norm2(a):
  7    return float(svdvals(a)[0])
  8
  9
 10def make_case(p=32, n=128, seed=0):
 11    rng = np.random.default_rng(seed)
 12    X = rng.normal(size=(n, p))
 13    H = X.T @ X / n
 14    K = rng.normal(size=(p, p))
 15    np.fill_diagonal(K, 0.0)
 16    # explicitly non-reciprocal, with a substantial antisymmetric component
 17    K = K - K.T * 0.35 + 0.35 * (K + K.T) / 2
 18    K /= norm2(K)
 19    return H, K
 20
 21
 22def susceptibility(J, gamma=0.02):
 23    p = J.shape[0]
 24    R = np.linalg.inv(gamma * np.eye(p) - J)
 25    return float(np.mean(np.abs(R) ** 2))
 26
 27
 28def transient_gain(J, times):
 29    return max(norm2(expm(J * t)) for t in times)
 30
 31
 32def run_mechanism():
 33    p, n = 32, 128
 34    q = p / n
 35    mp_delta = (1 - math.sqrt(q)) ** 2
 36    H, K = make_case(p, n, 7)
 37    observed_delta = float(np.min(np.linalg.eigvalsh(H)))
 38    k_norm = norm2(K)
 39    pred_ec = mp_delta / k_norm
 40    empirical_ec = observed_delta / k_norm
 41    eps_grid = np.linspace(0, 1.5, 31)
 42    rows = []
 43    for eps in eps_grid:
 44        J = -H + eps * K
 45        ev = eigvals(J)
 46        rows.append({
 47            'eps': float(eps),
 48            'ratio': float(eps * k_norm / observed_delta),
 49            'chi': susceptibility(J),
 50            'max_real_eig': float(np.max(ev.real)),
 51            'gain': float(transient_gain(J, np.linspace(0, 20, 41)))
 52        })
 53    # For the claimed crossover, report the first point whose chi differs
 54    # from passive chi by 25 percent (either suppression or amplification).
 55    chi0 = rows[0]['chi']
 56    changed = [r for r in rows if abs(r['chi'] / chi0 - 1) >= .25]
 57    observed_chi_cross = changed[0]['eps'] if changed else float('nan')
 58    # instability boundary is a separate falsifiable prediction: eigenvalues
 59    # should remain stable for purely dissipative perturbations, but not in general.
 60    unstable = [r for r in rows if r['max_real_eig'] >= 0]
 61    instability_eps = unstable[0]['eps'] if unstable else float('nan')
 62    return {
 63        'q': q, 'mp_delta': mp_delta, 'observed_delta': observed_delta,
 64        'k_norm': k_norm, 'predicted_ec_mp': pred_ec,
 65        'predicted_ec_empirical_gap': empirical_ec,
 66        'observed_chi_25pct_eps': observed_chi_cross,
 67        'observed_instability_eps': instability_eps,
 68        'rows': rows
 69    }
 70
 71
 72def train_linear_demo(seed=11, steps=500, penalty=False):
 73    # Same tiny supervised long-horizon state-space problem for both methods.
 74    rng = np.random.default_rng(seed)
 75    p, T, batch = 12, 35, 64
 76    torch = __import__('torch')
 77    torch.manual_seed(seed)
 78    device = 'cuda' if torch.cuda.is_available() else 'cpu'
 79    try:
 80        dtype = torch.float32
 81        Atrue = torch.tensor(rng.normal(size=(p,p)), dtype=dtype, device=device)
 82        Atrue = Atrue / (1.15 * torch.linalg.matrix_norm(Atrue))
 83        Btrue = torch.tensor(rng.normal(size=(p,1)), dtype=dtype, device=device)
 84        W = torch.randn(p,p, device=device) * .08
 85        U = torch.randn(p,1, device=device) * .08
 86        W.requires_grad_(); U.requires_grad_()
 87        opt = torch.optim.Adam([W,U], lr=.025)
 88        last_loss = None
 89        for step in range(steps):
 90            x = torch.randn(batch,T,1,device=device)
 91            with torch.no_grad():
 92                h = torch.zeros(batch,p,device=device)
 93                ys=[]
 94                for t in range(T):
 95                    h = torch.tanh(h @ Atrue.T + x[:,t] @ Btrue.T)
 96                    ys.append(h)
 97                y = torch.stack(ys,1)
 98            h = torch.zeros(batch,p,device=device)
 99            pred=[]
100            for t in range(T):
101                h = torch.tanh(h @ W.T + x[:,t] @ U.T)
102                pred.append(h)
103            pred=torch.stack(pred,1)
104            loss=((pred-y)**2).mean()
105            if penalty:
106                # Exact small-batch local resolvent proxy at a fixed damping.
107                # Use W's symmetric Gram passive estimate and residual activity.
108                H = (W.T @ W) / p
109                J = W
110                R = torch.linalg.inv(.08*torch.eye(p,device=device)-J)
111                chi=(R*R).mean()
112                loss = loss + .002 * torch.relu(torch.log1p(chi)-torch.tensor(math.log(3.0),device=device))**2
113            opt.zero_grad(); loss.backward(); torch.nn.utils.clip_grad_norm_([W,U], 5.0); opt.step()
114            last_loss=float(loss.detach().cpu())
115        with torch.no_grad():
116            ev=torch.linalg.eigvals(W).real.max().item()
117            # sampled transient gain using torch matrix exponentials
118            gain=max(torch.linalg.matrix_norm(torch.matrix_exp(W*t),2).item() for t in np.linspace(0,8,17))
119            R=torch.linalg.inv(.08*torch.eye(p,device=device)-W)
120            chi=(R*R).mean().item()
121        return {'loss':last_loss,'max_real_eig':ev,'gain':gain,'chi':chi,'device':device}
122    except Exception as exc:
123        # CUDA failures and numerical backend failures fall back to CPU.
124        if device == 'cuda':
125            os.environ['CUDA_VISIBLE_DEVICES']=''
126            return train_linear_demo(seed, steps, penalty)
127        raise exc
128
129
130if __name__ == '__main__':
131    out={'mechanism':run_mechanism(), 'baseline':train_linear_demo(penalty=False), 'idea':train_linear_demo(penalty=True)}
132    with open('results.json','w') as f: json.dump(out,f,indent=2)
133    m=out['mechanism']
134    print(json.dumps({'mechanism_summary':{k:v for k,v in m.items() if k!='rows'}, 'baseline':out['baseline'],'idea':out['idea']},indent=2))