"""Reproducible exploratory calculations; see PROTOCOLE.md for scope.

Run with the scientific Python environment. No network access from this script.
Outputs are calculations, not a certification of a theory of gravity.
"""
from __future__ import annotations

import csv
import hashlib
import json
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.interpolate import PchipInterpolator
from scipy.optimize import minimize_scalar

BASE = Path(__file__).resolve().parent
OUT = BASE / "outputs"
QUAD_X, QUAD_W = np.polynomial.legendre.leggauss(96)
U, UW = (QUAD_X + 1) / 2, QUAD_W / 2
A0_REFERENCE = 1.2e-10


def mu_standard(x):
    x = np.asarray(x, dtype=float)
    return x / np.hypot(1, x)


def mu_interface(x, variant="linear"):
    """Critical stiffness-normalized response; x=g/g_star, x must be >0.

    This is the derivative of the specified fermionic ground-state functional,
    angular averaged over isotropic interfaces. tanh is a postulated coupling.
    """
    raw = np.asarray(x, dtype=float)
    shape = raw.shape
    v = raw.reshape(-1)
    if np.any(v <= 0):
        raise ValueError("g/g_star must be positive")
    result = np.empty_like(v)
    small = v < 1e-3
    z = v[small]
    if variant == "linear":
        result[small] = .75*z - .3*z**2 + 3*z**4/56
        q = v[~small, None] * U
        # Rationalization avoids subtraction of two large terms at large x.
        result[~small] = 1 - 3 * np.sum(UW*U**2/(np.hypot(1, q)+q), axis=1)
    elif variant == "tanh":
        result[small] = .75*z + .5*z**2 - (5/6)*z**3
        w = v[~small]
        upper = np.minimum(w, 24.)
        t = upper[:, None] * U
        f = np.tanh(t)
        fp = 1 / np.cosh(t)**2
        integrand = t*f*fp/(np.hypot(1, f)+f)
        integral = upper * np.sum(integrand*UW, axis=1)
        result[~small] = 1 - 3*integral/w**3
    else:
        raise ValueError(variant)
    return result.reshape(shape)


def solve_standard(y):
    """g/a0 for the spherical equation mu(g/a0)*g/a0 = gN/a0."""
    y = np.asarray(y, dtype=float)
    return np.sqrt(.5*y*(y + np.hypot(y, 2)))


def bisect_response(y, mu):
    lo, hi = 0., max(2., 2*y)
    while float(mu(hi))*hi < y:
        hi *= 2
    for _ in range(100):
        mid = (lo+hi)/2
        if float(mu(mid))*mid < y:
            lo = mid
        else:
            hi = mid
    return (lo+hi)/2


def interface_inverse(variant):
    x = np.logspace(-9, 9, 18001)
    y = x*mu_interface(x, variant)
    if not np.all(np.diff(y) > 0):
        raise ValueError(f"Nonmonotonic static response: {variant}")
    interpolator = PchipInterpolator(np.log(y), np.log(x), extrapolate=False)
    return lambda y: np.exp(interpolator(np.log(y)))


def save_figure(fig, name):
    fig.savefig(OUT / f"{name}.png", dpi=170, bbox_inches="tight")
    fig.savefig(OUT / f"{name}.pdf", bbox_inches="tight")
    plt.close(fig)


def local_scalar():
    """Green function of a local nearest-neighbour action, no input 1/r."""
    records, curves, refinement = [], [], []
    radii = np.arange(4, 13)
    reference = 1/(4*np.pi*radii**2)
    for n in (64, 96, 128):
        wave = 2*np.pi*np.fft.fftfreq(n)
        one = 4*np.sin(wave/2)**2
        operator = one[:, None, None]+one[None, :, None]+one[None, None, :]
        inv = np.zeros_like(operator)
        np.divide(1., operator, out=inv, where=operator > 0)
        green = np.fft.ifftn(inv).real
        force = .5*(green[radii-1, 0, 0]-green[radii+1, 0, 0])
        error = force/reference-1
        record = dict(n=n, rms_relative_error=float(np.sqrt(np.mean(error**2))),
                      max_relative_error=float(np.max(np.abs(error))),
                      force_log_slope=float(np.polyfit(np.log(radii), np.log(force), 1)[0]))
        records.append(record)
        curves.append((n, force/reference))
        fixed_r = n//16
        fixed_force = .5*(green[fixed_r-1, 0, 0]-green[fixed_r+1, 0, 0])
        refinement.append({"n": n, "r_in_cells": fixed_r, "r_over_box": fixed_r/n,
                           "relative_error": float(fixed_force*4*np.pi*fixed_r**2-1)})
    fig, ax = plt.subplots(figsize=(7, 4))
    for n, ratio in curves:
        ax.plot(radii, ratio, "o-", label=f"N={n}")
    ax.axhline(1, c="black", lw=1)
    ax.set(xlabel="Separation / lattice spacing", ylabel="Lattice force / continuum force",
           title="Local scalar operator: approach to inverse-square force")
    ax.legend()
    save_figure(fig, "local_scalar")
    return {"grids": records, "descriptive_10pct_gate": records[-1]["max_relative_error"] < .10,
            "posthoc_fixed_physical_separation_diagnostic": refinement,
            "limitations": "Static scalar only; mass-proportional coupling assumed; periodic neutralizing background."}


def spectral_obstruction():
    rng = np.random.default_rng(20260905)
    r = np.logspace(-3, 3, 400)
    violations = 0
    max_increase = -np.inf
    examples = []
    for trial in range(100):
        masses = 10**rng.uniform(-3, 3, 30)
        masses[0] = 0  # Include a massless pole.
        weights = rng.lognormal(0, 2, 30)
        weights /= weights.sum()
        mr = r[:, None]*masses
        boost = np.sum(weights*(1+mr)*np.exp(-mr), axis=1)
        max_increase = max(max_increase, float(np.max(np.diff(boost))))
        violations += int(np.any(np.diff(boost) > 1e-13))
        if trial < 6:
            examples.append(boost)
    fig, ax = plt.subplots(figsize=(7, 4))
    for curve in examples:
        ax.semilogx(r, curve, alpha=.7)
    ax.set(xlabel="Radius (arbitrary units)", ylabel="r² |F| / total source coupling",
           title="Positive attractive spectral mixtures cannot grow in the infrared")
    save_figure(fig, "spectral_obstruction")
    return {"mixtures": 100, "modes_per_mixture": 30, "monotonicity_violations": violations,
            "maximum_discrete_increase": max_increase,
            "proof": "d/dr sum w(1+mr)exp(-mr) = -r sum w m² exp(-mr) <= 0",
            "scope": "Linear, translation-invariant, universal attractive positive spectral representation; no matter halo or environmental response."}


def static_and_quantum_checks():
    y = np.logspace(-8, 8, 81)
    numerical = np.array([bisect_response(t, mu_standard) for t in y])
    relative = float(np.max(np.abs(numerical/solve_standard(y)-1)))
    mass_checks = []
    for m in (.03, .2, .8, 2.):
        numeric = -quad(lambda k: k*(np.hypot(k, m)-k), 0, 1, epsabs=1e-12)[0]
        analytic = -((1+m*m)**1.5-1-abs(m)**3)/3
        mass_checks.append(dict(m=m, integral=numeric, analytic=analytic,
                                relative_error=abs(numeric/analytic-1)))
    x = np.logspace(-5, 5, 601)
    responses = {"AQUAL standard": mu_standard(x)}
    checks = {}
    for name in ("linear", "tanh"):
        mu = mu_interface(x, name)
        h = 1e-4
        derivative = (mu_interface(x*np.exp(h), name)-mu_interface(x*np.exp(-h), name))/(2*h)
        longitudinal = mu+derivative
        inverse = interface_inverse(name)
        ys = np.logspace(-6, 6, 31)
        exact = np.array([bisect_response(t, lambda z: mu_interface(z, name)) for t in ys])
        checks[name] = {
            "mu_min": float(mu.min()), "longitudinal_stiffness_min": float(longitudinal.min()),
            "static_positive_on_scan": bool(np.all(mu > 0) and np.all(longitudinal > 0)),
            "small_x_mu_over_x": float(mu[0]/x[0]),
            "large_x_deficit_coefficient": float((1-mu[-1])*x[-1]) if name == "linear" else None,
            "scalar_completion_c_parallel_squared_low_x": float(longitudinal[0]/mu[0]),
            "inverse_interpolation_relative_error": float(np.max(np.abs(inverse(ys)/exact-1))),
        }
        energy_errors = []
        def energy(w):
            def band(u):
                m = w*u if name == "linear" else np.tanh(w*u)
                return np.expm1(1.5*np.log1p(m*m))-abs(m)**3
            return .5*w*w-quad(band, 0, 1, epsabs=1e-12, epsrel=1e-12)[0]
        for w in (.01, .1, 1., 3., 10.):
            step = 1e-4*w
            slope = (energy(w-2*step)-8*energy(w-step)+8*energy(w+step)-energy(w+2*step))/(12*step)
            energy_errors.append(abs(slope/w/float(mu_interface(w, name))-1))
        checks[name]["energy_derivative_max_relative_error"] = max(energy_errors)
        responses[f"Quantum interfaces: {name}"] = mu
    # tanh asymptote, evaluated as an integral to avoid rounding 1-mu to zero.
    c_inf = 3*quad(lambda t: t*np.tanh(t)/np.cosh(t)**2 /
                  (np.hypot(1, np.tanh(t))+np.tanh(t)), 0, 24, epsabs=1e-12)[0]
    checks["tanh"]["asymptotic_deficit_coefficient_cubed"] = c_inf
    fig, axes = plt.subplots(1, 2, figsize=(11, 4))
    for label, mu in responses.items():
        axes[0].loglog(x, mu, label=label)
    axes[0].set(xlabel="Gradient / model scale", ylabel="Static response μ",
                title="Static response: known law and declared quantum toy")
    axes[0].legend(fontsize=8)
    for eps in (0., 1e-4, 1e-2):
        mu = (mu_interface(x, "tanh")+eps)/(1+eps)
        axes[1].loglog(x, mu, label=f"Critical detuning ε={eps:g}")
    axes[1].set(xlabel="g / g*", ylabel="μ", title="MOND asymptote requires criticality")
    axes[1].legend(fontsize=8)
    save_figure(fig, "quantum_response")
    return {"standard_inverse_max_relative_error": relative,
            "standard_inverse_gate": relative < 1e-10,
            "fermion_band_integral": mass_checks,
            "fermion_integral_gate": max(x["relative_error"] for x in mass_checks) < 1e-8,
            "interface_response": checks,
            "common_scalar_cone_gate": False,
            "conformal_scalar_lensing": "Phi_scalar + Psi_scalar = 0; no extra photon deflection at leading order."}


def load_rar():
    rows = []
    path = BASE / "data" / "RAR.mrt"
    for line in path.read_text().splitlines():
        fields = line.split()
        if len(fields) != 4:
            continue
        try:
            rows.append([float(v) for v in fields])
        except ValueError:
            pass
    data = np.asarray(rows)
    if data.shape[0] < 1000 or data.shape[1] != 4 or not np.all(np.isfinite(data)):
        raise ValueError("Unexpected RAR dataset; inspect the source file.")
    return data, hashlib.sha256(path.read_bytes()).hexdigest()


def empirical_rar(gbar, a0):
    return gbar / (-np.expm1(-np.sqrt(gbar/a0)))


def fit_rar():
    data, checksum = load_rar()
    logbar, errbar, logobs, errobs = data.T
    gbar = 10**logbar
    inverses = {variant: interface_inverse(variant) for variant in ("linear", "tanh")}
    models = {
        "AQUAL_standard": lambda g, a: a*solve_standard(g/a),
        "empirical_RAR": empirical_rar,
        "interfaces_linear": lambda g, s: s*inverses["linear"](g/s),
        "interfaces_tanh": lambda g, s: s*inverses["tanh"](g/s),
    }
    scores = {"Newton_baryons_only": {
        "rms_dex": float(np.sqrt(np.mean((logobs-logbar)**2))),
        "median_residual_dex": float(np.median(logobs-logbar)),
        "fitted_parameters": 0, "warning": "Not a LambdaCDM halo fit."}}
    grid = np.logspace(logbar.min()-.15, logbar.max()+.15, 400)
    curves = {}
    residuals = {"Newton_baryons_only": logobs-logbar}
    for name, fn in models.items():
        def objective(logscale):
            residual = np.log10(fn(gbar, 10**logscale))-logobs
            return float(np.mean(residual**2))
        fit = minimize_scalar(objective, bounds=(-12, -8), method="bounded",
                              options={"xatol": 1e-10})
        if not fit.success:
            raise RuntimeError(f"Fit did not converge: {name}")
        scale = 10**float(fit.x)
        residual = logobs-np.log10(fn(gbar, scale))
        scores[name] = {
            "scale_m_s2": scale, "scale_name": "g_star" if "interfaces" in name else "a0",
            "deep_mond_a0_m_s2": 4*scale/3 if "interfaces" in name else scale,
            "rms_dex": float(np.sqrt(np.mean(residual**2))),
            "median_residual_dex": float(np.median(residual)),
            "fitted_parameters": 1,
        }
        curves[name] = fn(grid, scale)
        residuals[name] = residual
    fig, axes = plt.subplots(1, 2, figsize=(12, 4.6), gridspec_kw={"width_ratios": [1.4, 1]})
    axes[0].hexbin(logbar, logobs, gridsize=48, mincnt=1, bins="log", cmap="Greys", alpha=.75)
    axes[0].plot(np.log10(grid), np.log10(grid), "k--", lw=1, label="Newton: baryons only")
    colors = ["#0072B2", "#009E73", "#D55E00", "#CC79A7"]
    for (name, curve), color in zip(curves.items(), colors):
        axes[0].plot(np.log10(grid), np.log10(curve), label=name, color=color, lw=1.6)
    axes[0].set(xlabel="log10 baryonic acceleration [m/s²]",
                ylabel="log10 observed acceleration [m/s²]",
                title="Public SPARC RAR: descriptive fits, one global scale")
    axes[0].legend(fontsize=7)
    names = list(scores)
    axes[1].barh(range(len(names)), [scores[n]["rms_dex"] for n in names],
                 color=["#888888"]+colors)
    axes[1].set_yticks(range(len(names)), names, fontsize=8)
    axes[1].set(xlabel="RMS residual [dex]", title="In-sample scores; no covariance model")
    fig.tight_layout()
    save_figure(fig, "rar_comparison")
    with (OUT / "rar_residuals.csv").open("w", newline="") as f:
        writer = csv.writer(f)
        writer.writerow(["log10_gbar", "error_log10_gbar", "log10_gobs", "error_log10_gobs"]+list(residuals))
        for i, row in enumerate(data):
            writer.writerow(list(row)+[float(residuals[n][i]) for n in residuals])
    return {"source": "https://astroweb.case.edu/SPARC/RAR.mrt", "sha256": checksum,
            "points": len(data), "galaxy_ids_present": False,
            "validation": "Descriptive in-sample only: public table has no galaxy IDs. No row-level split.",
            "loss": "Unweighted squared log residual; errors provided but not used to fake independent likelihood.",
            "limitations": "Correlated radii; baryonic mass assumptions retained; spherical algebraic law used as RAR proxy, not full AQUAL disc solutions; not a comparison with dark-matter halo models.",
            "scores": scores}


def physical_scales():
    G, msun, c = 6.67430e-11, 1.98847e30, 299792458.
    kpc, rsun = 3.085677581491367e19, 6.957e8
    galactic_mass = 1e11*msun
    return {
        "a0_reference_m_s2": A0_REFERENCE,
        "galaxy_mass_msun": 1e11,
        "flat_speed_km_s": float((G*galactic_mass*A0_REFERENCE)**.25/1000),
        "transition_radius_kpc": float(np.sqrt(G*galactic_mass/A0_REFERENCE)/kpc),
        "solar_limb_deflection_arcsec": {
            "GR_gamma_1": float(4*G*msun/(rsun*c*c)*206264.806247),
            "tempo_only_gamma_0": float(2*G*msun/(rsun*c*c)*206264.806247),
            "pure_conformal_scalar_gamma_minus1": 0.},
        "interpretation": "Benchmarks and consequences of assumed laws, not emergent relativistic results."}


def main():
    OUT.mkdir(exist_ok=True)
    summary = {
        "status": "Exploratory comparison; no complete relativistic quantum Le Sage theory constructed.",
        "local_scalar": local_scalar(),
        "spectral_obstruction": spectral_obstruction(),
        "static_quantum": static_and_quantum_checks(),
        "rar": fit_rar(),
        "reference_scales": physical_scales(),
    }
    (OUT / "results.json").write_text(json.dumps(summary, indent=2, ensure_ascii=False)+"\n")
    print(json.dumps(summary, indent=2, ensure_ascii=False))


if __name__ == "__main__":
    main()
