1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2518
  6np.random.seed(SEED); random.seed(SEED)
  7
  8
  9def cocycle_run(l1, l2, K=200, noise=0.20, theta0=None):
 10    """Random diagonal positive cocycle; log singular increments have known means."""
 11    if theta0 is None:
 12        theta0 = np.random.uniform(0.25, 1.25, size=2)
 13    # tangent vectors are represented by slopes y/x; projective distance is angle.
 14    slopes = np.tan(theta0)
 15    logs = np.zeros(2)
 16    angles = [theta0.copy()]
 17    q = np.eye(2)
 18    qr_diag = []
 19    for t in range(K):
 20        # independent temporal randomness, but the same scalar derivative applies to both vectors
 21        inc = np.array([l1, l2]) + noise*np.random.randn(2)
 22        A = np.diag(np.exp(inc))
 23        # QR iteration estimates both exponents
 24        z = A @ q
 25        q, R = np.linalg.qr(z)
 26        # enforce positive diagonal, without changing Q directions
 27        signs = np.sign(np.diag(R)); signs[signs == 0] = 1
 28        q = q @ np.diag(signs); R = np.diag(signs) @ R
 29        qr_diag.append(np.log(np.abs(np.diag(R))))
 30        logs += inc
 31        slopes *= np.exp(inc[1]-inc[0])
 32        # angles to the top axis (modulo orientation)
 33        angles.append(np.arctan(np.abs(slopes)))
 34    return np.asarray(angles), np.asarray(qr_diag), logs/K
 35
 36
 37def fit_slope(y, start=20, end=None):
 38    if end is None: end = len(y)
 39    x = np.arange(start, end)
 40    return float(np.polyfit(x, np.log(np.maximum(y[start:end], 1e-300)), 1)[0])
 41
 42
 43def toy_verification():
 44    # Prediction 1: projective log-angle slope is -gamma.
 45    gaps = [0.02, 0.05, 0.10, 0.20, 0.40]
 46    rows = []
 47    for g in gaps:
 48        # average many paths to suppress finite-sample noise
 49        slopes = []; estimates = []
 50        for _ in range(80):
 51            a, qr, means = cocycle_run(0.0, -g, K=180, noise=0.20)
 52            slopes.append(fit_slope(a[:,0], 30, 170))
 53            estimates.append(np.mean(qr[:,0]-qr[:,1]))
 54        obs_s = float(np.mean(slopes)); obs_g = float(np.mean(estimates))
 55        rows.append({'gamma_true':g, 'predicted_slope':-g, 'observed_slope':obs_s,
 56                     'predicted_gap':g, 'observed_qr_gap':obs_g,
 57                     'slope_rel_error':abs(obs_s+g)/g})
 58
 59    # Prediction 2: time to reach angle epsilon scales log(theta0/epsilon)/gamma.
 60    threshold_rows = []
 61    eps = 1e-3
 62    for g in [0.05, 0.10, 0.20, 0.40]:
 63        measured=[]; predicted=[]
 64        for _ in range(100):
 65            a, _, _ = cocycle_run(0.0, -g, K=500, noise=0.20,
 66                                  theta0=np.array([0.8, 0.8]))
 67            hit = np.flatnonzero(a[:,0] <= eps)
 68            measured.append(int(hit[0]) if len(hit) else 500)
 69            predicted.append(math.log(math.tan(0.8)/eps)/g)
 70        threshold_rows.append({'gamma':g, 'predicted_steps':float(np.mean(predicted)),
 71                              'observed_steps':float(np.mean(measured)),
 72                              'ratio_observed_predicted':float(np.mean(measured)/np.mean(predicted))})
 73
 74    # Prediction 3: no gap means no projective contraction; with exactly equal rates,
 75    # the angle remains constant even though the common exponent can be negative.
 76    zero=[]
 77    for l in [0.0, -0.2]:
 78        a, qr, means = cocycle_run(l, l, K=150, noise=0.20)
 79        zero.append({'common_lambda':l, 'angle_log_slope':fit_slope(a[:,0], 10, 140),
 80                     'qr_gap':float(np.mean(qr[:,0]-qr[:,1])),
 81                     'angle_ratio_final_initial':float(a[-1,0]/a[0,0])})
 82    return {'gap_sweep':rows, 'threshold_sweep':threshold_rows, 'zero_gap':zero}
 83
 84
 85# Optional practical check: train a tiny tanh RNN, regularizing a finite-window
 86# Jacobian gap estimated by explicit 2x2 Jacobians (small enough to be transparent).
 87def rnn_check(steps=180, seq_len=40, batch=32):
 88    try:
 89        import torch
 90        torch.manual_seed(SEED)
 91        device = 'cuda' if torch.cuda.is_available() else 'cpu'
 92        class RNN(torch.nn.Module):
 93            def __init__(self):
 94                super().__init__(); self.W=torch.nn.Parameter(torch.randn(12,12)*.18)
 95                self.U=torch.nn.Parameter(torch.randn(12,2)*.18); self.b=torch.nn.Parameter(torch.zeros(12))
 96                self.out=torch.nn.Linear(12,1)
 97            def forward(self,x, regularize=False):
 98                z=torch.zeros(x.shape[0],12,device=x.device); gs=[]; l1s=[]
 99                for t in range(x.shape[1]):
100                    z=torch.tanh(z@self.W.T+x[:,t]@self.U.T+self.b)
101                    # Exact Jacobian singular values for the hidden map at this state.
102                    if regularize:
103                        D=1-z*z
104                        J=D.unsqueeze(2)*self.W.unsqueeze(0)
105                        s=torch.linalg.svdvals(J)
106                        l1s.append(torch.log(s[:,0].clamp_min(1e-8)))
107                        gs.append(torch.log((s[:,0]/s[:,1]).clamp_min(1e-8)))
108                return self.out(z), (torch.stack(l1s,1).mean() if regularize else None), (torch.stack(gs,1).mean() if regularize else None)
109        def run(reg):
110            torch.manual_seed(SEED+int(reg)); m=RNN().to(device); opt=torch.optim.Adam(m.parameters(),lr=2e-3)
111            losses=[]; gaps=[]; exps=[]
112            for k in range(steps):
113                x=torch.rand(batch,seq_len,2,device=device); target=x[:,:,0].sum(1)+x[:,:,1].sum(1)
114                y,l1,g=m(x,reg); loss=((y[:,0]-target)**2).mean()
115                if reg: loss=loss+0.5*torch.relu(0.05-g)**2+0.5*l1*l1
116                opt.zero_grad(); loss.backward(); torch.nn.utils.clip_grad_norm_(m.parameters(),1.0); opt.step()
117                losses.append(float(((y[:,0]-target)**2).mean().detach().cpu()))
118                if reg: gaps.append(float(g.detach().cpu())); exps.append(float(l1.detach().cpu()))
119            return {'final_loss':float(np.mean(losses[-20:])), 'gap':float(np.mean(gaps[-20:])) if reg else None,
120                    'lambda1':float(np.mean(exps[-20:])) if reg else None, 'device':device}
121        return {'baseline':run(False),'regularized':run(True)}
122    except Exception as e:
123        return {'error':repr(e), 'fallback':'toy verification still valid'}
124
125
126if __name__ == '__main__':
127    result={'seed':SEED, 'toy':toy_verification(), 'rnn':rnn_check()}
128    Path('results.json').write_text(json.dumps(result, indent=2))
129    print(json.dumps(result, indent=2))