"""Cassini/Mercury compatibility attempt and a metric-extension viability gate.

Published marginal constraints are NOT independent likelihoods. The metric
sector is a new phenomenological assumption, not a quantum derivation.
"""
from __future__ import annotations

import csv
import hashlib
import json
from pathlib import Path

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq, minimize_scalar

from dark_transition_tweak import Tweak, REFERENCE
import dark_mercury_check as mercury

BASE = Path(__file__).resolve().parent
OUT = BASE / "outputs/dark_joint_solar"
ALPHA_K = 1e-9
LAMBDA_K = .01
B_K = 0.0
PREFERRED_ALPHA1_BOUND = 1e-4  # Quoted Solar-System bound in Bonetti & Barausse.


def digest(path):
    return hashlib.sha256(Path(path).read_bytes()).hexdigest()


def response_offsets(law, x):
    """mu-1 and (mu+x*mu')-1 without subtracting rounded ones at high x."""
    x = np.asarray(x, dtype=float)
    m, r = law.evaluate(x)
    safe = np.maximum(x, 100.)
    inv_t = law.position / safe
    correction = (27/4)*law.amplitude * inv_t**4 / (1 + inv_t**2)**3
    radial_correction = correction * 3 * (inv_t**2 - 1) / (1 + inv_t**2)
    tail_m = -REFERENCE.c / safe**3 + 2*REFERENCE.normalization*REFERENCE.d / safe**4 + correction
    tail_r = 2*REFERENCE.c / safe**3 - 6*REFERENCE.normalization*REFERENCE.d / safe**4 + radial_correction
    return np.where(x > 100, tail_m, m-1), np.where(x > 100, tail_r, r-1)


def lapse_coefficients(law, x, alpha=ALPHA_K):
    dm, dr = response_offsets(law, x)
    # F'(a)/(2a) and F''(a)/2 for the exact static matching function.
    return alpha - (2-alpha)*dm, alpha - (2-alpha)*dr


def ppn(alpha, lam, b=0.):
    a1 = 4*(alpha-2*b)/(b-1)
    a2 = (alpha-2*b)*(-b*(3+b+3*lam)-lam+alpha*(1+b+2*lam)) / ((alpha-2)*(b-1)*(b+lam))
    return {"beta_PPN": 1., "gamma_PPN": 1., "alpha1_PPN": a1, "alpha2_PPN": a2,
            "tensor_speed_squared_over_c_squared": 1/(1-b),
            "scalar_high_acceleration_speed_squared_over_c_squared": (alpha-2)*(b+lam)/(alpha*(b-1)*(2+b+3*lam))}


def principal_reduction(lapse_coefficient, lam):
    """Numerical Schur elimination of shift and lapse, b_K=0.

    This is the leading local quadratic action on a slowly varying
    stationary background, not a full background-dependent mode spectrum.
    Spatial coefficient is evaluated for one wave-vector direction.
    """
    lambda_h = 1 + lam
    kinetic = np.array([[3*(1-3*lambda_h), 3*lambda_h-1],
                        [3*lambda_h-1, 1-lambda_h]])
    spatial = np.array([[2., 2.], [2., lapse_coefficient]])
    reduced_kinetic = kinetic[0, 0] - kinetic[0, 1]**2 / kinetic[1, 1]
    reduced_spatial = spatial[0, 0] - spatial[0, 1]**2 / spatial[1, 1]
    speed_squared = -reduced_spatial / reduced_kinetic
    formula = lam/(2+3*lam) * (2-lapse_coefficient)/lapse_coefficient
    assert abs(speed_squared/formula-1) < 2e-11
    return {"reduced_kinetic": float(reduced_kinetic),
            "reduced_spatial": float(reduced_spatial),
            "speed_squared_over_c_squared": float(speed_squared)}


def geodesic_control():
    """General beta/gamma geodesic acceleration against secular formula."""
    sample = mercury.orbit(256)
    r, v, radius = sample["r"], sample["v"], sample["radius"]
    speed_squared = np.einsum("ij,ij->i", v, v)
    rv = np.einsum("ij,ij->i", r, v)
    gr = 3*mercury.MEAN_MOTION*mercury.GM / (
        mercury.LIGHT_SPEED**2*mercury.SEMIMAJOR*(1-mercury.ECCENTRICITY**2))
    rows = []
    for beta, gamma in [(1., 1.), (.9, 1.1), (1.2, .8), (0., 0.)]:
        acceleration = mercury.GM/(mercury.LIGHT_SPEED**2*radius[:, None]**3) * (
            (2*(beta+gamma)*mercury.GM/radius-gamma*speed_squared)[:, None]*r
            + 2*(1+gamma)*rv[:, None]*v)
        measured = mercury.secular_rates(sample, acceleration)["in_plane_apsidal_rad_s"]
        predicted = gr*(2+2*gamma-beta)/3
        error = abs(measured/predicted-1)
        assert error < 1e-10
        rows.append({"beta_PPN": beta, "gamma_PPN": gamma, "relative_error": float(error)})
    return rows


def shapiro_control(gamma):
    """Illustrative straight-ray first-PN delay, not Cassini tracking data."""
    r_earth, r_probe, impact = mercury.AU, 8*mercury.AU, 2*695700000.
    left = -np.sqrt(r_earth**2-impact**2)/impact
    right = np.sqrt(r_probe**2-impact**2)/impact
    numerical, quadrature_error = quad(lambda u: 1/np.sqrt(1+u*u), left, right,
                                       points=[0.], epsabs=1e-12, epsrel=1e-12)
    exact = np.arcsinh(right)-np.arcsinh(left)
    prefactor = 2*(1+gamma)*mercury.GM/mercury.LIGHT_SPEED**3
    assert abs(numerical/exact-1) < 1e-12
    return {"scope": "Synthetic geometry: Earth 1 au, probe 8 au, impact 2 nominal solar radii",
            "round_trip_delay_seconds": float(prefactor*numerical),
            "analytic_relative_error": float(abs(numerical/exact-1)),
            "quadrature_error_seconds": float(prefactor*quadrature_error),
            "largest_endpoint_radius_m": r_probe,
            "upper_path_length_m": r_earth+r_probe}


def response_peak(law, component):
    x = np.geomspace(1e-5, 1e5, 20001)
    delta = response_offsets(law, x)[component]
    i = int(np.argmax(delta))
    assert 0 < i < len(x)-1
    fit = minimize_scalar(lambda z: -float(response_offsets(law, np.exp(z))[component]),
                          bounds=(np.log(x[i-1]), np.log(x[i+1])), method="bounded",
                          options={"xatol": 1e-12})
    assert fit.success
    return float(np.exp(fit.x)), float(1-fit.fun)


def negative_intervals(law, component):
    x = np.geomspace(1e-5, 1e7, 48001)
    values = lapse_coefficients(law, x)[component]
    crossings = np.flatnonzero((values[:-1] < 0) != (values[1:] < 0))
    roots = [brentq(lambda z: float(lapse_coefficients(law, np.exp(z))[component]),
                    np.log(x[i]), np.log(x[i+1]), xtol=1e-12) for i in crossings]
    edges = [np.log(x[0]), *roots, np.log(x[-1])]
    return [[float(np.exp(a)), float(np.exp(b))] for a, b in zip(edges[:-1], edges[1:])
            if lapse_coefficients(law, np.exp((a+b)/2))[component] < 0]


def main():
    pub_path = BASE / "outputs/dark_transition_tweak/publication_summary.json"
    mercury_path = BASE / "outputs/dark_mercury/summary.json"
    pub = json.loads(pub_path.read_text())
    previous = json.loads(mercury_path.read_text())
    assert digest(BASE / "dark_mercury_check.py") == previous["source_sha256"]
    for name, expected in previous["input_sha256"].items():
        assert digest(BASE / name) == expected, name
    candidate = pub["candidate"]
    law = Tweak(candidate["amplitude"], candidate["position"])
    metric = ppn(ALPHA_K, LAMBDA_K, B_K)
    assert abs(metric["alpha1_PPN"]) < PREFERRED_ALPHA1_BOUND
    assert abs(metric["alpha2_PPN"]) < 1e-7
    light = shapiro_control(metric["gamma_PPN"])

    observations = {
        "cassini_Q2": {"mean": 1.6e-27, "sigma": 1.8e-27, "units": "s^-2",
                       "source": "https://arxiv.org/abs/2602.17884v2"},
        "cassini_gamma_minus_one": {"mean": 2.1e-5, "sigma": 2.3e-5,
                                   "source": "https://www.nature.com/articles/nature01997"},
        "messenger_beta_minus_one": {"mean": -2.7e-5, "sigma": 3.9e-5,
            "source": "https://doi.org/10.3847/1538-3881/aa5be2",
            "scope": "Published abstract, already conditional on the Cassini gamma prior; Table 1 prints -2.6e-5 with the same error",
            "independent_of_cassini_gamma": False}}
    gamma_obs = observations["cassini_gamma_minus_one"]
    beta_obs = observations["messenger_beta_minus_one"]
    gamma_offset = (metric["gamma_PPN"]-1-gamma_obs["mean"])/gamma_obs["sigma"]
    beta_offset = (metric["beta_PPN"]-1-beta_obs["mean"])/beta_obs["sigma"]
    cases = []
    gr_rate = previous["GR_reference_arcsec_century"]
    radial = abs(previous["radial_tail"]["total_arcsec_century"])
    for row in previous["quadrupole_cases"]:
        q2 = row["Q2_s_minus2"]
        peri_bound = row["max_abs_longitude_precession_arcsec_century"] + radial
        beta_bias_bound = 3*peri_bound/gr_rate
        q_offset = (q2-1.6e-27)/1.8e-27
        # Separate marginal screens; no joint likelihood or p-value.
        marginal_compatible = (abs(q_offset) <= 2 and abs(gamma_offset) <= 2
                               and abs(beta_offset)+beta_bias_bound/beta_obs["sigma"] <= 2)
        # Conservative straight-ray bound if the added potential has gamma=1.
        delay_bound = (2*(1+abs(metric["gamma_PPN"]))*abs(q2)
                       * light["largest_endpoint_radius_m"]**2/3
                       * light["upper_path_length_m"]/mercury.LIGHT_SPEED**3)
        cases.append({"external_m_s2": row["external_acceleration_m_s2"], "Q2_s_minus2": q2,
                      "Q2_offset_in_quoted_sigma": q_offset,
                      "mercury_additional_precession_bound_arcsec_century": peri_bound,
                      "equivalent_beta_bias_bound": beta_bias_bound,
                      "illustrative_tidal_roundtrip_delay_bound_seconds": delay_bound,
                      "inside_each_quoted_two_sigma_marginal": bool(marginal_compatible)})

    x_peak, radial_peak = response_peak(law, 1)
    a_t, a_r = map(float, lapse_coefficients(law, x_peak))
    reduced = principal_reduction(a_r, LAMBDA_K)
    reductions = []
    for lam in [1e-4, .001, .01, .03]:
        for coefficient in [a_r, .1, .5, 1.5]:
            reductions.append({"lambda_K": lam, "lapse_coefficient": coefficient,
                               **principal_reduction(coefficient, lam)})
    # Independent derivative of x*mu at the problematic point.
    step = x_peak*1e-5
    stiffness_fd = ((x_peak+step)*law.evaluate(x_peak+step)[0]
                    -(x_peak-step)*law.evaluate(x_peak-step)[0])/(2*step)
    assert abs(stiffness_fd-radial_peak) < 1e-7
    # Static equation matching: (1-chi/2)/(1-alpha_K/2) = mu.
    probe_x = np.geomspace(1e-3, 100, 201)
    mu = law.evaluate(probe_x)[0]
    chi = lapse_coefficients(law, probe_x)[0]
    matching_error = float(np.max(abs((1-chi/2)/(1-ALPHA_K/2)-mu)))
    assert matching_error < 1e-14
    required_alpha = 2*(1-1/radial_peak)

    shapes_path = BASE / "outputs/dark_transition_tweak/all_shape_fits.csv"
    solar_path = BASE / "outputs/dark_transition_tweak/all_solar_results.csv"
    with shapes_path.open() as f:
        shapes = list(csv.DictReader(f))
    with solar_path.open() as f:
        solar = list(csv.DictReader(f))
    screened_shapes = []
    for row in shapes:
        amplitude, position = float(row["amplitude"]), float(row["position"])
        test_law = Tweak(amplitude, position)
        xm, peak_mu = response_peak(test_law, 0)
        xr, peak_radial = response_peak(test_law, 1)
        maximum = max(peak_mu, peak_radial)
        minimum_alpha = max(0., 2*(1-1/maximum))
        any_cassini = any(float(s["amplitude"]) == amplitude and float(s["position"]) == position
                          and -2e-27 <= float(s["Q2_s_minus2"]) <= 5.2e-27 for s in solar)
        screened_shapes.append({**{key: float(value) for key, value in row.items()},
            "max_mu": peak_mu, "max_radial": peak_radial,
            "x_at_max_mu": xm, "x_at_max_radial": xr,
            "minimum_alpha_K_for_nonnegative_principal_lapse_coefficients": minimum_alpha,
            "any_archived_solar_case_inside_cassini_band": bool(any_cassini),
            "compatible_with_quoted_alpha1_bound_in_bK_zero_family": bool(4*minimum_alpha <= PREFERRED_ALPHA1_BOUND)})
    inputs = [pub_path, mercury_path, shapes_path, solar_path, BASE / "dark_mercury_check.py",
              BASE / "dark_transition_tweak.py", BASE / "dark_band_report.py",
              BASE / "outputs/dark_lattice_bands.csv", BASE / "outputs/dark_lattice_controls.csv",
              BASE / "CASSINI_MERCURY_PROTOCOL.md"]
    report = {
        "date": "2026-09-05",
        "status": "Conditional Solar-System marginal compatibility; candidate metric extension fails local principal scalar stability diagnostic",
        "original_fitted_model": {key: candidate[key] for key in ["amplitude", "position", "galaxy_scale_m_s2", "rar_rms_dex"]},
        "metric_extension": {"action_c_equals_one": "S=(16*pi*Gbare)^-1 integral N*sqrt(h) [Kij*Kij-(1+lambda_K)*K^2+R3+F(a)] + Sm[g,matter]",
            "F_matching": "F(a)=2*a^2-2*(2-alpha_K)*a_star^2*U(a/a_star)+constant; U'(x)=x*mu(x)",
            "a_star_inverse_m": candidate["galaxy_scale_m_s2"]/mercury.LIGHT_SPEED**2,
            "alpha_K": ALPHA_K, "lambda_K": LAMBDA_K, "b_K": B_K,
            "Gbare_over_Gmeasured": 1-ALPHA_K/2,
            "assumptions": ["Known metric/foliation action imported, not derived from LQ8 bands",
                            "Universal matter coupling to this metric postulated",
                            "Kinetic values illustrative; Mercury does not independently determine them",
                            "No nonlocal or higher-spatial-derivative operators in this candidate",
                            "Additive constant chosen to neglect cosmological curvature in this local test"],
            "high_acceleration_PPN": metric,
            "PPN_scope": "Known constant-coupling high-acceleration limit; no full transition-region relativistic boundary-value solve"},
        "observations": observations,
        "gamma_marginal_offset_in_quoted_sigma": float(gamma_offset),
        "beta_marginal_offset_in_quoted_sigma": float(beta_offset),
        "no_joint_chi_squared_computed": True,
        "conditional_cases": cases,
        "conditional_all_marginal_screens_pass": all(row["inside_each_quoted_two_sigma_marginal"] for row in cases),
        "Mercury_GR_benchmark_arcsec_century": gr_rate,
        "geodesic_controls": geodesic_control(),
        "shapiro_control": light,
        "stationary_AQUAL_matching_max_error": matching_error,
        "stability": {"scope": "Leading local quadratic action on slowly varying stationary backgrounds; wavelengths shorter than background variation scale within the two-derivative effective theory",
            "lapse_coefficients": "alpha_transverse=2-(2-alpha_K)*mu; alpha_parallel=2-(2-alpha_K)*(mu+x*mu')",
            "scalar_speed_formula_bK_zero": "cs^2/c^2=lambda_K/(2+3*lambda_K)*(2-alpha_direction)/alpha_direction",
            "x_at_max_radial_stiffness": x_peak, "max_radial_stiffness": radial_peak,
            "radial_stiffness_finite_difference_error": float(abs(stiffness_fd-radial_peak)),
            "lapse_transverse_at_peak": a_t, "lapse_parallel_at_peak": a_r,
            "parallel_mode_at_peak": reduced,
            "negative_parallel_lapse_intervals_x": negative_intervals(law, 1),
            "negative_transverse_lapse_intervals_x": negative_intervals(law, 0),
            "minimum_alpha_K_to_remove_parallel_failure": required_alpha,
            "resulting_abs_alpha1_PPN_at_minimum": 4*required_alpha,
            "quoted_solar_abs_alpha1_bound": PREFERRED_ALPHA1_BOUND,
            "elimination_controls": reductions,
            "passes_local_principal_stability_gate": bool(reduced["speed_squared_over_c_squared"] > 0)},
        "archive_shape_screen": screened_shapes,
        "scope_limits": ["No raw Cassini or MESSENGER tracking data, covariance matrix, ephemeris refit, or independent joint likelihood",
            "The published MESSENGER beta constraint already includes Cassini's gamma information",
            "PPN beta=gamma=1 comes from the added known metric framework, not from quantum emergence",
            "Preferred-frame parameters checked against quoted bounds, not integrated into a full Mercury ephemeris",
            "Higher multipoles, finite-acceleration PN corrections and global metric/foliation response remain uncomputed",
            "The stability failure applies to this explicit minimal metric extension, not every possible nonlocal completion",
            "Passing the marginal observational screens is not presented as a viable simultaneous theory fit"],
        "sources": {"metric_and_PPN": "https://arxiv.org/abs/1502.05554v3",
                    "original_metric_framework": "https://arxiv.org/abs/1107.5264",
                    "related_kinetic_stability_conditions": "https://arxiv.org/abs/2302.14846"},
        "source_sha256": digest(__file__),
        "input_sha256": {str(path.relative_to(BASE)): digest(path) for path in inputs}}
    OUT.mkdir(parents=True, exist_ok=True)
    (OUT / "summary.json").write_text(json.dumps(report, indent=2, allow_nan=False)+"\n")
    with (OUT / "archived_shapes_metric_screen.csv").open("w") as f:
        writer = csv.DictWriter(f, fieldnames=list(screened_shapes[0]))
        writer.writeheader()
        writer.writerows(screened_shapes)
    x = np.geomspace(.001, 1e5, 2401)
    m, r = law.evaluate(x)
    at, ar = lapse_coefficients(law, x)
    np.savetxt(OUT / "selected_metric_response.csv", np.column_stack((x, m, r, at, ar)),
               delimiter=",", header="x,mu,radial,alpha_transverse,alpha_parallel", comments="")
    print(json.dumps({"conditional_marginal_compatibility": report["conditional_all_marginal_screens_pass"],
        "gamma_offset_sigma": gamma_offset, "beta_offset_sigma": beta_offset,
        "max_radial": radial_peak, "x_at_peak": x_peak, "parallel_lapse": a_r,
        "parallel_cs_squared_over_c_squared": reduced["speed_squared_over_c_squared"],
        "minimum_alpha_K_for_repair": required_alpha,
        "shapes_examined": len(screened_shapes),
        "shapes_passing_both_stability_and_quoted_alpha1_bound": sum(row["compatible_with_quoted_alpha1_bound_in_bK_zero_family"] for row in screened_shapes)}, indent=2))


if __name__ == "__main__":
    main()
