import json, math
from pathlib import Path
import numpy as np

SEED = 12345
ALPHA = 1.5

def stable_symmetric(rng, n, alpha=ALPHA, scale=1.0):
    v = rng.uniform(-np.pi/2, np.pi/2, size=n)
    w = rng.exponential(1.0, size=n)
    x = (np.sin(alpha*v) / np.cos(v)**(1.0/alpha) *
         (np.cos((1-alpha)*v) / w)**((1-alpha)/alpha))
    return scale*x

def robust_mad(x):
    med = np.median(x)
    return float(np.median(np.abs(x-med)))

def boundary_sweep():
    qs = np.array([.25, .5, 1., 1.5, 1.9, 2.0, 2.1, 2.5])
    rows = []
    for q in qs:
        a = 1-q
        x = 1.0
        for _ in range(80):
            x *= a
        logmag = math.log10(abs(x)) if x != 0 else -300.0
        rows.append({'q': float(q), 'abs_multiplier': abs(a),
                     'log10_abs_x80': float(logmag),
                     'predicted_stable': bool(0 < q < 2),
                     'observed_stable': bool(logmag < -3)})
    return rows

def stationary_sweep():
    rng = np.random.default_rng(SEED)
    q, a, n, burn = .8, .2, 220000, 3000
    sigmas = np.array([.25, .5, 1., 2.])
    mads, coeffs = [], []
    us = np.array([.25, .5, .75, 1.0])
    for sigma in sigmas:
        x = 0.0
        samples = np.empty(n)
        for t in range(n + burn):
            x = a*x + float(stable_symmetric(rng, 1, ALPHA, sigma)[0])
            if t >= burn: samples[t-burn] = x
        mads.append(robust_mad(samples))
        phi = np.array([np.mean(np.cos(u*samples)) for u in us])
        # -log|phi(u)| = C |u|^alpha; regress through zero.
        y = -np.log(np.maximum(phi, 1e-12))
        coeffs.append(float(np.dot(us**ALPHA, y) / np.dot(us**ALPHA, us**ALPHA)))
    slope = float(np.polyfit(np.log(sigmas), np.log(mads), 1)[0])
    cf_slope = float(np.polyfit(np.log(sigmas), np.log(coeffs), 1)[0])
    # Discrete AR prediction; continuous formula is approached as q is small.
    predicted_cf_coeff_at_sigma1 = 1.0/(1.0-abs(a)**ALPHA)
    return {'q': q, 'alpha': ALPHA, 'sigmas': sigmas.tolist(),
            'mads': [float(x) for x in mads], 'mad_loglog_slope': slope,
            'cf_coefficients': coeffs, 'cf_coefficient_loglog_slope': cf_slope,
            'predicted_mad_slope': 1.0, 'predicted_cf_slope': ALPHA,
            'predicted_discrete_cf_coeff_sigma1': predicted_cf_coeff_at_sigma1}

def cf_alpha_sweep():
    rng = np.random.default_rng(SEED + 9)
    results = []
    us = np.array([.35, .45, .55, .70, .85, 1.0, 1.2])
    for alpha in [1.2, 1.5, 1.8]:
        x = stable_symmetric(rng, 1000000, alpha, 1.0)
        phi = np.array([np.mean(np.cos(u*x)) for u in us])
        y = -np.log(np.maximum(phi, 1e-12))
        fit = np.polyfit(np.log(us), np.log(y), 1)
        results.append({'true_alpha': alpha, 'fitted_cf_power': float(fit[0]), 'fit_intercept': float(fit[1])})
    return results

def optimizer_comparison():
    # Same 1-D quadratic and same jump scale; report median final |x| over seeds.
    # This is illustrative, not a claim of universal optimization improvement.
    out = []
    for method in ['gaussian', 'levy']:
        finals = []
        for seed in range(30):
            rng = np.random.default_rng(9000 + seed)
            x, m = 5.0, 5.0
            eta, lam, beta, noise = .04, .8, .95, .035
            for _ in range(500):
                g = x
                m = beta*m + (1-beta)*x
                z = (rng.normal() if method == 'gaussian' else stable_symmetric(rng, 1, ALPHA)[0])
                x = x - eta*g - eta*lam*(x-m) + noise*z
            finals.append(abs(x))
        out.append({'method': method, 'median_final_abs_x': float(np.median(finals)),
                    'iqr_final_abs_x': [float(np.percentile(finals,25)), float(np.percentile(finals,75))]})
    return out

def main():
    result = {'seed': SEED, 'boundary': boundary_sweep(),
              'stationary_scaling': stationary_sweep(),
              'cf_power': cf_alpha_sweep(),
              'quadratic_optimizer': optimizer_comparison()}
    Path('results.json').write_text(json.dumps(result, indent=2))
    print(json.dumps(result, indent=2))

if __name__ == '__main__':
    main()