1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 3112
  6
  7
  8def cov_long_memory(n, alpha):
  9    i = np.arange(n)
 10    return (1.0 + np.abs(i[:, None] - i[None, :])) ** (-alpha)
 11
 12
 13def correlated_gaussian(n, alpha, rng):
 14    C = cov_long_memory(n, alpha)
 15    return np.linalg.cholesky(C + 1e-10 * np.eye(n)) @ rng.standard_normal(n)
 16
 17
 18def hermite_prob(x, q):
 19    if q == 1: return x
 20    if q == 2: return x*x - 1.0
 21    if q == 3: return x*x*x - 3.0*x
 22    raise ValueError('q must be 1, 2, or 3')
 23
 24
 25def hermite_unit(x, q):
 26    # For standard Gaussian input, E[He_q(G)^2] = q!.
 27    return hermite_prob(x, q) / math.sqrt(math.factorial(q))
 28
 29
 30def slope(x, y):
 31    return float(np.polyfit(np.log(np.asarray(x)), np.log(np.asarray(y)), 1)[0])
 32
 33
 34def scaling_check(alpha=0.4, q=1, max_n=512, trials=160, rng=None):
 35    # No per-path centering: centering would force sum(z)==0 and destroy
 36    # exactly the phenomenon being measured.
 37    rng = np.random.default_rng(SEED + q) if rng is None else rng
 38    ns = np.array([16, 32, 64, 128, 256, max_n])
 39    variances = []
 40    means = []
 41    for n in ns:
 42        sums = []
 43        for _ in range(trials):
 44            g = correlated_gaussian(n, alpha, rng)
 45            sums.append(np.sum(hermite_unit(g, q)))
 46        variances.append(np.var(sums, ddof=1))
 47        means.append(np.mean(sums))
 48    H = 1.0 - alpha*q/2.0
 49    observed_H = slope(ns, variances) / 2.0
 50    rms_scaled = [float(np.sqrt(v) / n**H) for n, v in zip(ns, variances)]
 51    # For Gaussian g, Cov(He_q(g_i)/sqrt(q!), He_q(g_j)/sqrt(q!))=rho(|i-j|)^q.
 52    exact_variances = []
 53    for n in ns:
 54        rho = (1.0 + np.arange(n)) ** (-alpha)
 55        exact_variances.append(float(n + 2.0 * sum((n-j) * rho[j]**q for j in range(1, n))))
 56    exact_H = slope(ns, exact_variances) / 2.0
 57    return {'q': q, 'theory_H': H, 'observed_H': observed_H, 'exact_finite_n_H': exact_H,
 58            'ns': ns.tolist(), 'sum_variance': [float(v) for v in variances],
 59            'exact_sum_variance': exact_variances,
 60            'scaled_sum_rms': rms_scaled,
 61            'rms_ratio_last_first': rms_scaled[-1] / rms_scaled[0],
 62            'sum_means': [float(x) for x in means]}
 63
 64
 65def residual_stability(alpha=0.4, n_trials=300, dims=32, rng=None):
 66    rng = np.random.default_rng(SEED + 99) if rng is None else rng
 67    depths = [32, 64, 128, 256]
 68    out = {}
 69    for kind in ['critical_correlated', 'iid_sqrt', 'correlated_inverse_sqrt', 'correlated_wrong_inverse']:
 70        vals = []
 71        for L in depths:
 72            final = []
 73            for _ in range(n_trials):
 74                if kind == 'critical_correlated':
 75                    z = correlated_gaussian(L, alpha, rng)
 76                    lam = L ** (-(1-alpha/2))
 77                elif kind == 'iid_sqrt':
 78                    z = rng.standard_normal(L)
 79                    lam = L ** -0.5
 80                elif kind == 'correlated_inverse_sqrt':
 81                    z = correlated_gaussian(L, alpha, rng)
 82                    lam = L ** -0.5
 83                else:
 84                    z = correlated_gaussian(L, alpha, rng)
 85                    lam = L ** -1.0
 86                h = rng.standard_normal(dims)
 87                for zl in z:
 88                    h = h + lam * zl * np.tanh(h)
 89                final.append(np.mean(h*h))
 90            vals.append((float(np.mean(final)), float(np.std(final, ddof=1)/math.sqrt(n_trials))))
 91        out[kind] = {'depths': depths, 'mean_final_variance': [x[0] for x in vals],
 92                     'sem': [x[1] for x in vals]}
 93    return out
 94
 95
 96def main():
 97    random.seed(SEED); np.random.seed(SEED)
 98    checks = [scaling_check(q=q) for q in (1, 2)]
 99    stability = residual_stability()
100    report = {'seed': SEED, 'alpha': 0.4, 'checks': checks, 'stability': stability}
101    Path('results.json').write_text(json.dumps(report, indent=2))
102    print(json.dumps(report, indent=2))
103
104if __name__ == '__main__':
105    main()