1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 17
  6np.random.seed(SEED)
  7random.seed(SEED)
  8
  9
 10def regulator(q, alpha, s):
 11    x = np.asarray(q, dtype=float) ** 2 / (s * s)
 12    return alpha * s * s / np.expm1(np.minimum(x, 700.0))
 13
 14
 15def crossover_root(alpha, s):
 16    # Unique positive solution of R_s(q)=q.
 17    lo, hi = 1e-12, max(10.0 * s, 10.0 * alpha * s * s + 1.0)
 18    for _ in range(100):
 19        mid = (lo + hi) / 2.0
 20        if regulator(mid, alpha, s) > mid:
 21            lo = mid
 22        else:
 23            hi = mid
 24    return (lo + hi) / 2.0
 25
 26
 27def stability_boundary(q, alpha, s):
 28    # For x_{t+1}=(1-eta*q/(q+R))x_t, stability is eta < 2/max(rate).
 29    rates = q / (q + regulator(q, alpha, s))
 30    return 2.0 / np.max(rates)
 31
 32
 33def exact_relaxation(q, alpha, s, eta):
 34    rate = eta * q / (q + regulator(q, alpha, s))
 35    multiplier = np.abs(1.0 - rate)
 36    t = -1.0 / np.log(np.maximum(multiplier, 1e-300))
 37    return t
 38
 39
 40def local_exponent(q, t):
 41    return -np.gradient(np.log(t), np.log(q))
 42
 43
 44def run_math_checks():
 45    alpha, s = 0.7, 1.0
 46    q = np.logspace(-4, 2, 1000)
 47    # Prediction A: R ~ alpha*s^4/q^2 at q << s; R/q -> 0 at q >> s.
 48    low = q < 0.01 * s
 49    high = q > 5.0 * s
 50    low_ratio = float(np.median(regulator(q[low], alpha, s) /
 51                               (alpha * s**4 / q[low]**2)))
 52    high_ratio = float(np.median(regulator(q[high], alpha, s) / q[high]))
 53
 54    # Prediction B: crossover solves R=q and has scale q*=s*y(alpha*s).
 55    cases = [(0.2, 0.7), (0.7, 1.0), (1.5, 1.4), (3.0, 0.8)]
 56    cross = []
 57    for a, ss in cases:
 58        pred = crossover_root(a, ss)
 59        grid = np.logspace(-7, 4, 50000) * max(ss, a * ss * ss, 1.0)
 60        obs = float(grid[np.argmin(np.abs(regulator(grid, a, ss) - grid))])
 61        cross.append({'alpha': a, 's': ss, 'predicted_q': pred,
 62                      'observed_q': obs, 'relative_error': abs(obs-pred)/pred})
 63    # Dimensionless scaling check: q*/s is a function only of alpha*s.
 64    scaling = []
 65    for u in [0.15, 0.5, 1.0, 2.0]:
 66        vals = []
 67        for ss in [0.5, 1.0, 2.0]:
 68            aa = u / ss
 69            vals.append(crossover_root(aa, ss) / ss)
 70        scaling.append({'alpha_times_s': u, 'qstar_over_s': vals,
 71                        'spread': max(vals)-min(vals)})
 72
 73    # Prediction C: Euler stability boundary is eta*=2/max_q q/(q+R).
 74    qs = np.logspace(-5, 3, 50000)
 75    eta_star = stability_boundary(qs, alpha, s)
 76    rates = qs / (qs + regulator(qs, alpha, s))
 77    test_etas = [0.99 * eta_star, 1.01 * eta_star]
 78    stability = []
 79    for eta in test_etas:
 80        multiplier = np.max(np.abs(1.0 - eta * rates))
 81        stability.append({'eta': eta, 'predicted_stable': bool(eta < eta_star),
 82                          'observed_max_multiplier': float(multiplier),
 83                          'observed_stable': bool(multiplier < 1.0)})
 84
 85    # Claimed exponents: directly measure t(q) from the proposed update.
 86    # This is a falsification check; no exponent is imposed.
 87    t = exact_relaxation(qs, alpha, s, eta=0.5)
 88    z = local_exponent(qs, t)
 89    bands = [(1e-4, 3e-3), (0.03, 0.2), (5.0, 30.0)]
 90    exponents = [{'q_band': list(b), 'observed_z': float(np.median(z[(qs >= b[0]) & (qs <= b[1])]))}
 91                 for b in bands]
 92    return {'asymptotics': {'low_ratio_R_over_alpha_s4_over_q2': low_ratio,
 93                            'high_ratio_R_over_q': high_ratio},
 94            'crossovers': cross, 'dimensionless_scaling': scaling,
 95            'stability': {'eta_star': eta_star, 'tests': stability},
 96            'relaxation_exponents': exponents}
 97
 98
 99def adam_scalar(H, steps=400, lr=0.08):
100    x = np.ones_like(H) * 2.0
101    m = np.zeros_like(x); v = np.zeros_like(x)
102    for t in range(1, steps + 1):
103        g = H * x
104        m = 0.9*m + 0.1*g
105        v = 0.999*v + 0.001*g*g
106        x -= lr * (m/(1-0.9**t)) / (np.sqrt(v/(1-0.999**t)) + 1e-8)
107    return float(0.5*np.sum(H*x*x)), float(np.linalg.norm(x))
108
109
110def smooth_rg(H, alpha=0.7, s=1.0, eta=0.5, steps=400):
111    # Modewise curvature q=H and exact local damping a=q for this quadratic.
112    x = np.ones_like(H) * 2.0
113    for _ in range(steps):
114        g = H*x
115        x -= eta * g / (np.abs(H) + regulator(H, alpha, s) + 1e-8)
116    return float(0.5*np.sum(H*x*x)), float(np.linalg.norm(x))
117
118
119def mini_experiment():
120    # Same diagonal quadratic, spanning shells; compare fixed AdamW-like and RG rule.
121    H = np.logspace(-2, 2, 64)
122    rows = []
123    for s in [0.3, 1.0, 3.0]:
124        b_loss, b_norm = adam_scalar(H)
125        i_loss, i_norm = smooth_rg(H, s=s)
126        rows.append({'s': s, 'adam_loss': b_loss, 'smooth_rg_loss': i_loss,
127                     'adam_norm': b_norm, 'smooth_rg_norm': i_norm})
128    return rows
129
130
131def main():
132    result = {'seed': SEED, 'math_checks': run_math_checks(),
133              'mini_experiment': mini_experiment()}
134    Path('results.json').write_text(json.dumps(result, indent=2))
135    print(json.dumps(result, indent=2))
136
137if __name__ == '__main__':
138    main()