1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2748
  6np.random.seed(SEED)
  7
  8# Fold normal form from the paper: u_{t+1}=u+g*eps+b*u^2.
  9def fold_step(u, eps, g=1.0, b=1.0):
 10    return u + g * eps + b * u * u
 11
 12def fixed_point(eps, g=1.0, b=1.0):
 13    return -math.sqrt(-g * eps / b)
 14
 15def stable_gap(eps, g=1.0, b=1.0):
 16    u = fixed_point(eps, g, b)
 17    return 1.0 - abs(1.0 + 2.0 * b * u)
 18
 19def passage_time(eps, g=1.0, b=1.0, max_steps=10000000):
 20    u, n = 0.0, 0
 21    while u < 1.0 and n < max_steps:
 22        u = fold_step(u, eps, g, b)
 23        n += 1
 24    return n
 25
 26def log_slope(x, y):
 27    return float(np.polyfit(np.log(x), np.log(y), 1)[0])
 28
 29def jvp_scalar(v, u, eps, g=1.0, b=1.0):
 30    return (1.0 + 2.0 * b * u) * v
 31
 32def predictor(u, delta=0.01, g=1.0, b=1.0):
 33    # Exact scalar JVP, equivalent to one power iteration.
 34    lam = abs(jvp_scalar(1.0, u, 0.0, g, b))
 35    return math.pi / max(delta, 1.0 - lam), lam
 36
 37def residual(u, old):
 38    return abs(u - old) / (abs(u) + 1e-6)
 39
 40def stable_run(eps, mode, tol=1e-5, rtol=0.03, tau_threshold=25.0,
 41               max_steps=50000):
 42    """Inference on the stable branch, stopping when |u-u*| <= tol.
 43
 44    Fixed is a conservative cap, residual is the standard local early exit,
 45    and spectral additionally requires a small local relaxation forecast.
 46    """
 47    target = fixed_point(eps)
 48    u, old = 0.0, 0.0
 49    for t in range(1, max_steps + 1):
 50        old, u = u, fold_step(u, eps)
 51        r = residual(u, old)
 52        accurate = abs(u - target) <= tol
 53        if mode == 'fixed':
 54            if accurate: return t, True
 55        elif mode == 'residual':
 56            if accurate and r < rtol: return t, True
 57        else:
 58            tau, lam = predictor(u)
 59            # Stable branch has lam<1; reject critical/unstable estimates.
 60            if accurate and r < rtol and lam < 1.0 and tau < tau_threshold:
 61                return t, True
 62    return max_steps, False
 63
 64def oscillatory_fallback():
 65    # A stable rotation has oscillating coordinates; a fold scalar detector
 66    # should disable its predictor on this non-fold trajectory.
 67    z = np.array([1.0, 0.0]); theta, rho = 0.55, 0.99
 68    R = rho * np.array([[math.cos(theta), -math.sin(theta)],
 69                        [math.sin(theta), math.cos(theta)]])
 70    xs = []
 71    for _ in range(20):
 72        z = R @ z; xs.append(float(z[0]))
 73    sign_changes = sum(xs[i] * xs[i-1] < 0 for i in range(1, len(xs)))
 74    return sign_changes >= 1, sign_changes
 75
 76def main():
 77    # Prediction 1: gap ~ eps^(1/2); prediction 2: passage time ~ eps^(-1/2).
 78    eps = np.logspace(-5, -2, 12)
 79    gaps = np.array([stable_gap(-e) for e in eps])
 80    passages = np.array([passage_time(e) for e in eps], dtype=float)
 81    # Prediction 3: changing b gives opposite +/-1/2 prefactor exponents.
 82    e0 = 2e-5
 83    bs = np.array([0.25, 0.5, 1., 2., 4.])
 84    b_gaps = np.array([stable_gap(-e0, b=b) for b in bs])
 85    b_passages = np.array([passage_time(e0, b=b) for b in bs], dtype=float)
 86    products = b_gaps * b_passages
 87
 88    # Direct stable-branch test of pi/(1-lambda): actual threshold time has
 89    # the same inverse-gap scaling, while pi is an asymptotic calibration.
 90    rows = []
 91    for e in [1e-4, 3e-5, 1e-5]:
 92        target = fixed_point(-e); u = 0.0
 93        for t in range(1, 50000):
 94            u = fold_step(u, -e)
 95            if abs(u-target) <= 1e-5:
 96                tau, lam = predictor(u)
 97                rows.append({'eps':e, 'actual_steps':t,
 98                             'predicted_tau':tau, 'lambda_hat':lam,
 99                             'gap':1-lam})
100                break
101
102    workload = np.logspace(-5, -3, 20)
103    controller = {}
104    for mode in ['fixed', 'residual', 'spectral']:
105        vals = [stable_run(float(-e), mode)[0] for e in workload]
106        controller[mode] = {'mean_steps':float(np.mean(vals)),
107                            'median_steps':float(np.median(vals)), 'steps':vals}
108    result = {
109      'seed':SEED,
110      'predictions':{
111        'gap_epsilon_exponent':{'predicted':0.5,'observed':log_slope(eps,gaps)},
112        'passage_epsilon_exponent':{'predicted':-0.5,'observed':log_slope(eps,passages)},
113        'gap_b_exponent':{'predicted':0.5,'observed':log_slope(bs,b_gaps)},
114        'passage_b_exponent':{'predicted':-0.5,'observed':log_slope(bs,b_passages)},
115        'Pi':{'predicted':math.pi,'observed':products.tolist(),
116              'mean':float(np.mean(products)),
117              'relative_error':float(abs(np.mean(products)-math.pi)/math.pi)}},
118      'stable_predictor_samples':rows, 'controller':controller,
119      'oscillatory_fallback':dict(detected=oscillatory_fallback()[0],
120                                  sign_changes=oscillatory_fallback()[1]),
121      'settings':{'delta':0.01,'residual_tol':0.03,'accuracy_tol':1e-5,
122                  'tau_threshold':25.0}}
123    Path('results.json').write_text(json.dumps(result,indent=2))
124    print(json.dumps(result,indent=2))
125
126if __name__ == '__main__': main()