1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6# Toy RG velocity: a translation-invariant field whose real-space influence
  7# decays as exp(-Lambda * distance). This is the minimal kernel implied by
  8# quasi-local RG flow and lets us measure truncation error exactly.
  9SEED = 2708
 10L = 256
 11LAMBDAS = np.array([0.10, 0.16, 0.25, 0.40, 0.64, 1.00])
 12EPSILONS = np.array([0.20, 0.10, 0.05, 0.02, 0.01])
 13
 14
 15def periodic_kernel(lam, n=L):
 16    d = np.minimum(np.arange(n), n - np.arange(n)).astype(float)
 17    k = np.exp(-lam * d)
 18    # Normalize so the zero mode has unit gain; normalization cancels in the
 19    # relative truncation error but makes the operator easy to interpret.
 20    return k / k.sum()
 21
 22
 23def truncate(k, radius):
 24    n = len(k)
 25    d = np.minimum(np.arange(n), n - np.arange(n))
 26    out = np.where(d <= radius, k, 0.0)
 27    # Preserve the zero-mode response, as a local approximation would normally
 28    # absorb the omitted mass into its learned bias/gain.
 29    return out / out.sum()
 30
 31
 32def relative_operator_error(k, radius):
 33    """RMS error over all Fourier modes, normalized by exact operator power."""
 34    exact = np.fft.rfft(k)
 35    approx = np.fft.rfft(truncate(k, radius))
 36    return float(np.linalg.norm(approx - exact) / np.linalg.norm(exact))
 37
 38
 39def required_radius(lam, eps):
 40    k = periodic_kernel(lam)
 41    for r in range(L // 2):
 42        if relative_operator_error(k, r) <= eps:
 43            return r
 44    return L // 2
 45
 46
 47def verify_field_error(lam, radius, trials=64):
 48    rng = np.random.default_rng(SEED + int(1000 * lam) + radius)
 49    k = periodic_kernel(lam)
 50    kt = truncate(k, radius)
 51    errs = []
 52    for _ in range(trials):
 53        x = rng.normal(size=L)
 54        y = np.fft.irfft(np.fft.rfft(k) * np.fft.rfft(x), n=L)
 55        yt = np.fft.irfft(np.fft.rfft(kt) * np.fft.rfft(x), n=L)
 56        errs.append(np.linalg.norm(y - yt) / np.linalg.norm(y))
 57    return float(np.mean(errs))
 58
 59
 60def linear_fit(x, y):
 61    a, b = np.polyfit(x, y, 1)
 62    pred = a * x + b
 63    r2 = 1.0 - np.sum((y - pred) ** 2) / max(np.sum((y - y.mean()) ** 2), 1e-12)
 64    return float(a), float(b), float(r2)
 65
 66
 67def pyramid_radius_check():
 68    # Lambda_s = pi/a_s. The claimed cell-radius scaling is
 69    # R_s/a_s = c [log L + log(1/eps)] / pi, independent of s.
 70    eps = 0.02
 71    c = 1.0
 72    rows = []
 73    for s in range(6):
 74        a = 2 ** s
 75        lam = math.pi / a
 76        R = c * (math.log(L) + math.log(1 / eps)) / lam
 77        r = math.ceil(R / a)
 78        rows.append({'level': s, 'spacing': a, 'Lambda': lam,
 79                     'physical_R': R, 'cell_radius': r})
 80    return rows
 81
 82
 83def main():
 84    table = []
 85    for lam in LAMBDAS:
 86        for eps in EPSILONS:
 87            r = required_radius(lam, eps)
 88            table.append({'Lambda': float(lam), 'epsilon': float(eps),
 89                          'required_R': r,
 90                          'measured_error': relative_operator_error(periodic_kernel(lam), r),
 91                          'field_error': verify_field_error(lam, r)})
 92
 93    # Prediction 1: at fixed epsilon, R*Lambda is approximately constant.
 94    fixed = [q for q in table if q['epsilon'] == 0.02]
 95    inv_lam = np.array([1 / q['Lambda'] for q in fixed])
 96    radii = np.array([q['required_R'] for q in fixed])
 97    slope, intercept, r2 = linear_fit(inv_lam, radii)
 98    products = [q['Lambda'] * q['required_R'] for q in fixed]
 99
100    # Prediction 2: at fixed Lambda, R grows linearly with log(1/epsilon).
101    lam0 = 0.25
102    tolrows = [q for q in table if q['Lambda'] == lam0]
103    xlog = np.array([math.log(1 / q['epsilon']) for q in tolrows])
104    yR = np.array([q['required_R'] for q in tolrows])
105    slope_log, intercept_log, r2_log = linear_fit(xlog, yR)
106
107    # Prediction 3 / baseline: a fixed radius has increasingly bad error as the
108    # RG length grows (Lambda decreases), while radius selected by the bound
109    # stays below tolerance.
110    fixed_radius = 8
111    baseline = []
112    for lam in LAMBDAS:
113        err = relative_operator_error(periodic_kernel(lam), fixed_radius)
114        rg_r = required_radius(lam, 0.02)
115        rg_err = relative_operator_error(periodic_kernel(lam), rg_r)
116        baseline.append({'Lambda': float(lam), 'fixed_radius': fixed_radius,
117                         'fixed_radius_error': err, 'rg_radius': rg_r,
118                         'rg_error': rg_err})
119
120    result = {
121        'setup': {'L': L, 'kernel': 'normalized exp(-Lambda * periodic_distance)',
122                  'seed': SEED},
123        'radius_sweep': table,
124        'predictions': {
125            'inverse_cutoff': {
126                'statement': 'R is affine in 1/Lambda at fixed epsilon',
127                'fit_R_vs_1_over_Lambda': {'slope': slope, 'intercept': intercept, 'R2': r2},
128                'Lambda_times_R': products,
129                'mean': float(np.mean(products)), 'cv': float(np.std(products) / np.mean(products))},
130            'log_tolerance': {
131                'statement': 'R is affine in log(1/epsilon) at fixed Lambda',
132                'Lambda': lam0,
133                'fit_R_vs_log_inverse_epsilon': {'slope': slope_log, 'intercept': intercept_log, 'R2': r2_log}},
134            'rescaled_grid': {
135                'statement': 'R/a is approximately constant when Lambda=pi/a',
136                'levels': pyramid_radius_check()}
137        },
138        'baseline_comparison': baseline
139    }
140    Path('results.json').write_text(json.dumps(result, indent=2))
141    print(json.dumps(result, indent=2))
142
143
144if __name__ == '__main__':
145    main()