1import json, math, hashlib
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 12345
  6rng = np.random.default_rng(SEED)
  7
  8# Golden ratio is a Pisot number (its conjugate is -1/phi).
  9PHI = (1.0 + math.sqrt(5.0)) / 2.0
 10
 11
 12def beta_orbit(n, u0, beta=PHI):
 13    u = float(u0)
 14    out = np.empty(n, dtype=np.float64)
 15    for k in range(n):
 16        out[k] = u
 17        u = beta * u
 18        u -= math.floor(u)
 19    return out
 20
 21
 22def orbit_coefficients(n, u0, burn=10000):
 23    # Estimate empirical centering/scaling on a separate long orbit, as proposed.
 24    calib = beta_orbit(burn + n, u0, PHI)[burn:]
 25    mu, sigma = float(calib.mean()), float(calib.std())
 26    if sigma == 0:
 27        raise ValueError('degenerate orbit')
 28    vals = beta_orbit(n + burn, u0, PHI)[burn:]
 29    return (vals - mu) / sigma, mu, sigma
 30
 31
 32def make_projection(d, m, u0, gap=1, offset=0, coeff_n=300000):
 33    # Rows are z[offset + r*gap + i], exactly matching the proposed construction.
 34    need = offset + (m - 1) * gap + d + 1
 35    z, mu, sigma = orbit_coefficients(max(need, coeff_n), u0)
 36    P = np.empty((m, d), dtype=np.float64)
 37    for r in range(m):
 38        start = offset + r * gap
 39        P[r] = z[start:start+d] / math.sqrt(m)
 40    return P, mu, sigma
 41
 42
 43def pair_distortion(P, X):
 44    D = X[:, None, :] - X[None, :, :]
 45    orig = np.sum(D * D, axis=2)
 46    proj = np.einsum('nkd,md->nkm', D, P)
 47    proj = np.sum(proj * proj, axis=2)
 48    mask = orig > 1e-12
 49    ratios = proj[mask] / orig[mask]
 50    return float(np.max(np.abs(ratios - 1))), float(np.mean(np.abs(ratios - 1)))
 51
 52
 53def alpha_calibrate(P, X):
 54    a = np.sum(X * X)
 55    b = np.sum((X @ P.T) ** 2)
 56    return math.sqrt(a / b)
 57
 58
 59def corr_decay(z, maxlag=1000):
 60    z = z - z.mean()
 61    den = np.dot(z, z)
 62    cs = np.array([np.dot(z[:-k], z[k:]) / den for k in range(1, maxlag+1)])
 63    # Fit log envelope to lag 10..maxlag, robustly excluding zeros.
 64    lags = np.arange(1, maxlag+1)
 65    env = np.maximum.accumulate(np.abs(cs)[::-1])[::-1] # conservative not useful for fit
 66    sel = (lags >= 10) & (np.abs(cs) > 1e-5)
 67    slope = np.polyfit(lags[sel], np.log(np.abs(cs[sel])), 1)[0]
 68    return cs, float(math.exp(slope))
 69
 70
 71def variance_sweep(d=32, ms=(8,16,32,64,128), trials=3000):
 72    # Random isotropic vectors make the norm estimator's variance a direct JL diagnostic.
 73    out=[]
 74    for m in ms:
 75        vals=[]
 76        for t in range(trials):
 77            u0 = (0.017 + (t+1)*0.61803398875) % 1.0
 78            P,_,_ = make_projection(d,m,u0,gap=1,coeff_n=5000)
 79            x = rng.normal(size=d); x /= np.linalg.norm(x)
 80            vals.append(np.sum((P @ x)**2))
 81        vals=np.asarray(vals)
 82        out.append((m,float(vals.mean()),float(vals.var()),float(vals.var()*m)))
 83    return out
 84
 85
 86def main():
 87    d,m=24,16
 88    u0=0.3141592653
 89    z,mu,sigma=orbit_coefficients(120000,u0)
 90    cs,rho=corr_decay(z,500)
 91    # Prediction 1: gap-dependent correlation should decrease geometrically.
 92    gaps=[1,2,4,8,16,32]
 93    gap_corr=[]
 94    for g in gaps:
 95        # Correlation of same-column coefficients in adjacent rows.
 96        vals=np.array([np.corrcoef(z[:-g],z[g:])[0,1]])
 97        gap_corr.append((g,float(abs(vals[0])),float(abs(cs[g-1]))))
 98
 99    # Calibration search prediction: optimize finite seed/gap and compare to fixed choice.
100    X=rng.normal(size=(18,d)); X -= X.mean(0); X /= X.std(0)
101    candidates=[]
102    for si in range(30):
103        seed=(0.0314159 + si*0.2718281) % 1.0
104        for g in [1,2,4,8,16,32]:
105            P,_,_=make_projection(d,m,seed,gap=g,coeff_n=3000)
106            a=alpha_calibrate(P,X)
107            md,ad=pair_distortion(a*P,X)
108            candidates.append((md,ad,si,g,a))
109    candidates.sort()
110    best=candidates[0]
111    fixed=[x for x in candidates if x[1] == min(y[1] for y in candidates)][:1]
112    # fixed is not a fair baseline; use first seed gap=1 explicitly.
113    P0,_,_=make_projection(d,m,0.0314159,gap=1,coeff_n=3000)
114    a0=alpha_calibrate(P0,X); base=pair_distortion(a0*P0,X)
115
116    var=variance_sweep()
117    result={
118      'seed':SEED,'beta':PHI,'orbit_mu':mu,'orbit_sigma':sigma,
119      'prediction_correlation':{'claim':'|corr(z_t,z_{t+g})| decreases approximately geometrically with gap','fit_rho_per_step':rho,'observed':gap_corr},
120      'prediction_variance':{'claim':'Var(||Px||^2) approximately proportional to 1/m','observed_m_mean_var_mvar':var,'loglog_slope':float(np.polyfit(np.log([x[0] for x in var]),np.log([x[2] for x in var]),1)[0])},
121      'prediction_search':{'claim':'seed/gap calibration search lowers worst pairwise distortion','fixed_seed_gap1':{'max':base[0],'mean':base[1]},'best':{'max':best[0],'mean':best[1],'seed_index':best[2],'gap':best[3],'alpha':best[4]}},
122      'reproducibility_hash':hashlib.sha256(np.asarray(make_projection(d,m,0.123456,gap=7,coeff_n=3000)[0],dtype=np.float64).tobytes()).hexdigest()
123    }
124    Path('results.json').write_text(json.dumps(result,indent=2))
125    print(json.dumps(result,indent=2))
126
127if __name__=='__main__': main()