1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 12345
  6ALPHA = 1.5
  7
  8def stable_symmetric(rng, n, alpha=ALPHA, scale=1.0):
  9    v = rng.uniform(-np.pi/2, np.pi/2, size=n)
 10    w = rng.exponential(1.0, size=n)
 11    x = (np.sin(alpha*v) / np.cos(v)**(1.0/alpha) *
 12         (np.cos((1-alpha)*v) / w)**((1-alpha)/alpha))
 13    return scale*x
 14
 15def robust_mad(x):
 16    med = np.median(x)
 17    return float(np.median(np.abs(x-med)))
 18
 19def boundary_sweep():
 20    qs = np.array([.25, .5, 1., 1.5, 1.9, 2.0, 2.1, 2.5])
 21    rows = []
 22    for q in qs:
 23        a = 1-q
 24        x = 1.0
 25        for _ in range(80):
 26            x *= a
 27        logmag = math.log10(abs(x)) if x != 0 else -300.0
 28        rows.append({'q': float(q), 'abs_multiplier': abs(a),
 29                     'log10_abs_x80': float(logmag),
 30                     'predicted_stable': bool(0 < q < 2),
 31                     'observed_stable': bool(logmag < -3)})
 32    return rows
 33
 34def stationary_sweep():
 35    rng = np.random.default_rng(SEED)
 36    q, a, n, burn = .8, .2, 220000, 3000
 37    sigmas = np.array([.25, .5, 1., 2.])
 38    mads, coeffs = [], []
 39    us = np.array([.25, .5, .75, 1.0])
 40    for sigma in sigmas:
 41        x = 0.0
 42        samples = np.empty(n)
 43        for t in range(n + burn):
 44            x = a*x + float(stable_symmetric(rng, 1, ALPHA, sigma)[0])
 45            if t >= burn: samples[t-burn] = x
 46        mads.append(robust_mad(samples))
 47        phi = np.array([np.mean(np.cos(u*samples)) for u in us])
 48        # -log|phi(u)| = C |u|^alpha; regress through zero.
 49        y = -np.log(np.maximum(phi, 1e-12))
 50        coeffs.append(float(np.dot(us**ALPHA, y) / np.dot(us**ALPHA, us**ALPHA)))
 51    slope = float(np.polyfit(np.log(sigmas), np.log(mads), 1)[0])
 52    cf_slope = float(np.polyfit(np.log(sigmas), np.log(coeffs), 1)[0])
 53    # Discrete AR prediction; continuous formula is approached as q is small.
 54    predicted_cf_coeff_at_sigma1 = 1.0/(1.0-abs(a)**ALPHA)
 55    return {'q': q, 'alpha': ALPHA, 'sigmas': sigmas.tolist(),
 56            'mads': [float(x) for x in mads], 'mad_loglog_slope': slope,
 57            'cf_coefficients': coeffs, 'cf_coefficient_loglog_slope': cf_slope,
 58            'predicted_mad_slope': 1.0, 'predicted_cf_slope': ALPHA,
 59            'predicted_discrete_cf_coeff_sigma1': predicted_cf_coeff_at_sigma1}
 60
 61def cf_alpha_sweep():
 62    rng = np.random.default_rng(SEED + 9)
 63    results = []
 64    us = np.array([.35, .45, .55, .70, .85, 1.0, 1.2])
 65    for alpha in [1.2, 1.5, 1.8]:
 66        x = stable_symmetric(rng, 1000000, alpha, 1.0)
 67        phi = np.array([np.mean(np.cos(u*x)) for u in us])
 68        y = -np.log(np.maximum(phi, 1e-12))
 69        fit = np.polyfit(np.log(us), np.log(y), 1)
 70        results.append({'true_alpha': alpha, 'fitted_cf_power': float(fit[0]), 'fit_intercept': float(fit[1])})
 71    return results
 72
 73def optimizer_comparison():
 74    # Same 1-D quadratic and same jump scale; report median final |x| over seeds.
 75    # This is illustrative, not a claim of universal optimization improvement.
 76    out = []
 77    for method in ['gaussian', 'levy']:
 78        finals = []
 79        for seed in range(30):
 80            rng = np.random.default_rng(9000 + seed)
 81            x, m = 5.0, 5.0
 82            eta, lam, beta, noise = .04, .8, .95, .035
 83            for _ in range(500):
 84                g = x
 85                m = beta*m + (1-beta)*x
 86                z = (rng.normal() if method == 'gaussian' else stable_symmetric(rng, 1, ALPHA)[0])
 87                x = x - eta*g - eta*lam*(x-m) + noise*z
 88            finals.append(abs(x))
 89        out.append({'method': method, 'median_final_abs_x': float(np.median(finals)),
 90                    'iqr_final_abs_x': [float(np.percentile(finals,25)), float(np.percentile(finals,75))]})
 91    return out
 92
 93def main():
 94    result = {'seed': SEED, 'boundary': boundary_sweep(),
 95              'stationary_scaling': stationary_sweep(),
 96              'cf_power': cf_alpha_sweep(),
 97              'quadratic_optimizer': optimizer_comparison()}
 98    Path('results.json').write_text(json.dumps(result, indent=2))
 99    print(json.dumps(result, indent=2))
100
101if __name__ == '__main__':
102    main()