"""Assemble checked solar comparisons for the constrained metric candidate.

Published marginal summaries are correlated and do not constitute a joint
likelihood. This report concerns an assumed effective action, not emergence.
"""
import csv
import json
from pathlib import Path
import numpy as np
from dark_cmc_model import BASE, OUT, CMCResponse, candidate, digest
import dark_mercury_check as mercury


def read(name,source):
    report=json.loads((OUT/name).read_text())
    assert report["source_sha256"]==digest(BASE/source),source
    for path,expected in report.get("input_sha256",{}).items():
        assert digest(BASE/path)==expected,path
    return report


def mercury_tail(law,record,nodes):
    sample=mercury.orbit(nodes)
    gN=mercury.GM/sample["radius"]**2
    x=gN/record["galaxy_scale_m_s2"]
    offset=law.offsets(x)[0]
    acceleration=(gN*offset)[:,None]*sample["radial"]
    value=mercury.secular_rates(sample,acceleration)["longitude_perihelion_rad_s"]*mercury.RAD_S_TO_ARCSEC_CENTURY
    return {"nodes":nodes,"arcsec_century":value,"max_abs_mu_minus_one":float(abs(offset).max()),
            "max_extra_acceleration_m_s2":float(abs(gN*offset).max()),"x_range":[float(x.min()),float(x.max())]}


def multipole_envelopes(green,external_integrals,scale):
    # A norm bound on the complete 3D Gauss longitude-rate functional.
    sample=mercury.orbit(256)
    e=mercury.ECCENTRICITY
    s=np.sqrt(1-e*e)
    r=sample["radius"]
    cf,sf=sample["radial"][:,0],sample["radial"][:,1]
    prefactor=s/(mercury.MEAN_MOTION*mercury.SEMIMAJOR*e)
    kr=-prefactor*cf
    kt=prefactor*(1+r/(mercury.SEMIMAJOR*s*s))*sf
    kn=np.tan(mercury.INCLINATION/2)*r*(np.sin(mercury.ARG_PERIHELION)*cf
         +np.cos(mercury.ARG_PERIHELION)*sf)/(mercury.MEAN_MOTION*mercury.SEMIMAJOR**2*s)
    kernel=np.sqrt(kr*kr+kt*kt+kn*kn)
    length=np.sqrt(6.67430e-11*1.98847e30/scale)
    records=[]
    for degree in [3,4]:
        values=[item["dimensionless_potential_multipole_coefficients"][str(degree)]
                for grid in green for item in grid["integrals"]]
        external_values=[row["integral"]["dimensionless_potential_multipole_coefficients"][str(degree)]
                         for row in external_integrals]
        # Include every sampled magnitude and the central refinement spread.
        bound=max(map(abs,values+external_values))+max(values)-min(values)
        # |P_l|<=1 and |P_l'|<=l(l+1)/2 give a conservative gradient bound.
        angular_bound=np.sqrt(degree**2+(degree*(degree+1)/2)**2)
        force_bound=bound*scale*(r/length)**(degree-1)*angular_bound
        precession=float(np.dot(sample["weights"],kernel*force_bound)*mercury.RAD_S_TO_ARCSEC_CENTURY)
        records.append({"degree":degree,"dimensionless_coefficients_384_512_orders_2_4":values,
            "coefficient_absolute_spread":float(max(values)-min(values)),
            "coefficients_at_all_sampled_external_fields":external_values,
            "precession_orientation_envelope_arcsec_century":precession})
    return records


def main():
    reports={key:read(key+"_summary.json",source) for key,source in [
        ("fields","dark_cmc_model.py"),("controls","dark_cmc_model.py"),("ir","dark_cmc_model.py"),
        ("theory","dark_cmc_theory.py"),("3d","dark_cmc_3d.py"),
        ("variation","dark_cmc_variation.py"),("solar_checks","dark_cmc_solar_checks.py"),
        ("inner","dark_cmc_inner_check.py")]}
    law=CMCResponse()
    record=candidate(law)
    for stage in ["fields","controls","ir"]:
        for item in reports[stage]["cases"]:
            field=item["solar"]
            assert digest(BASE/field["field"])==field["field_sha256"]
            assert digest(BASE/field["table"])==field["table_sha256"]
            assert field["result"]["converged"]
            assert field["spec"]["source_sha256"]==digest(BASE/"dark_cmc_model.py")
    for item in reports["3d"]["fields"]:
        assert digest(OUT/f"field_3d_n{item['n']}.npz")==item["field_sha256"]
    tail=[mercury_tail(law,record,n) for n in [128,256]]
    assert abs(tail[0]["arcsec_century"]-tail[1]["arcsec_century"])<1e-25
    eigenvalues=np.linalg.eigvalsh(mercury.orientation_matrix())
    gr=3*mercury.MEAN_MOTION*mercury.GM/(mercury.LIGHT_SPEED**2*mercury.SEMIMAJOR
           *(1-mercury.ECCENTRICITY**2))*mercury.RAD_S_TO_ARCSEC_CENTURY
    observations={
        "Q2":{"mean":1.6e-27,"sigma":1.8e-27,"unit":"s^-2","source":"https://arxiv.org/abs/2602.17884v2"},
        "gamma_minus_one":{"mean":2.1e-5,"sigma":2.3e-5,"source":"https://www.nature.com/articles/nature01997"},
        "beta_minus_one":{"mean":-2.7e-5,"sigma":3.9e-5,"source":"https://doi.org/10.3847/1538-3881/aa5be2",
            "scope":"Abstract value; Table 1 gives -2.6e-5. This result already uses the Cassini gamma prior.",
            "independent_of_gamma":False}}
    gamma_offset=-observations["gamma_minus_one"]["mean"]/observations["gamma_minus_one"]["sigma"]
    beta_offset=-observations["beta_minus_one"]["mean"]/observations["beta_minus_one"]["sigma"]
    multipoles=multipole_envelopes(reports["solar_checks"]["green_checks"],
                reports["solar_checks"]["external_integrals"],record["galaxy_scale_m_s2"])
    higher_bound=sum(row["precession_orientation_envelope_arcsec_century"] for row in multipoles)
    rows=[]
    for item in reports["fields"]["cases"]:
        field=item["solar"]
        q=field["result"]["Q2_s_minus2"][1]
        tide=np.sort(q/mercury.MEAN_MOTION*eigenvalues*mercury.RAD_S_TO_ARCSEC_CENTURY)
        additional=float(max(abs(tide))+abs(tail[-1]["arcsec_century"])+higher_bound)
        bias=3*additional/gr
        offset=(q-observations["Q2"]["mean"])/observations["Q2"]["sigma"]
        inside=abs(offset)<=2 and abs(gamma_offset)<=2 and abs(beta_offset)+bias/observations["beta_minus_one"]["sigma"]<=2
        rows.append({"external_acceleration_m_s2":field["spec"]["external"],"nr":field["spec"]["nr"],
            "nu":field["spec"]["nu"],"Q2_s_minus2":q,"Q2_offset_in_quoted_sigma":offset,
            "longitude_tide_orientation_interval_arcsec_century":[float(tide[0]),float(tide[-1])],
            "additional_precession_envelope_arcsec_century":additional,"equivalent_beta_bias_envelope":bias,
            "inside_each_quoted_two_sigma_marginal":bool(inside)})
    central=min(rows,key=lambda r:abs(r["external_acceleration_m_s2"]-2.32e-10))
    qc=central["Q2_s_minus2"]
    control_rows=[{**{key:f["solar"]["spec"][key] for key in ["nr","nu","rmin","rmax","table_nodes"]},
                   "Q2_s_minus2":f["solar"]["result"]["Q2_s_minus2"][1]} for f in reports["controls"]["cases"]]
    max_control_change=max(abs(r["Q2_s_minus2"]-qc) for r in control_rows)
    assert max_control_change<.06e-27
    ir_changes=[abs(item["solar"]["result"]["Q2_s_minus2"][1]-control_rows[0]["Q2_s_minus2"])
                for item in reports["ir"]["cases"]]
    minimum_margin=min(min(r["Q2_s_minus2"]+2e-27,5.2e-27-r["Q2_s_minus2"]) for r in rows)
    assert all(r["inside_each_quoted_two_sigma_marginal"] for r in rows)
    assert minimum_margin>10*max_control_change
    report={"date":"2026-09-05","status":"Conditional Cassini and Mercury compatibility in a constrained effective metric model",
        "candidate":record,"action":reports["theory"]["action"],"response_checks":reports["theory"]["response"],
        "observations":observations,"high_acceleration_metric":reports["theory"]["metric"],
        "gamma_offset_in_quoted_sigma":gamma_offset,"beta_offset_in_quoted_sigma":beta_offset,
        "no_joint_likelihood_computed":True,"all_sampled_marginal_screens_pass":True,
        "central_solar_case":central,"external_field_cases":rows,
        "Mercury_first_PN_benchmark_arcsec_century":float(gr),"mercury_radial_tail":tail[-1],
        "mercury_multipoles":multipoles,
        "largest_extra_precession_envelope_arcsec_century":max(r["additional_precession_envelope_arcsec_century"] for r in rows),
        "largest_equivalent_beta_bias":max(r["equivalent_beta_bias_envelope"] for r in rows),
        "numerical_controls":{"solar":control_rows,"maximum_solar_control_shift_s_minus2":max_control_change,
            "minimum_sampled_margin_to_Q2_band_s_minus2":minimum_margin,
            "maximum_IR_sensitivity_shift_s_minus2":max(ir_changes),
            "Green_extractions":reports["solar_checks"]["green_checks"],
            "inner_grid_diagnostic":reports["solar_checks"]["inner_extended_grid"],
            "inner_precision_diagnostic":reports["inner"],
            "three_dimensional_final":{k:v for k,v in reports["3d"]["fields"][-1].items() if k!="history"},
            "maximum_nonlinear_variation_error_at_smallest_step":max(
                max(r["finite_variations"][-1]["lapse_absolute_error"],r["finite_variations"][-1]["conformal_absolute_error"])
                for r in reports["variation"]["checks"])},
        "assumptions":["Constrained metric action and universal matter coupling are postulated",
            "The galaxy transition A=-0.45,T=4.25 was selected using known Cassini and Gaia constraints",
            "The momentum constraint, q(t)=-2/t clock and ultraweak Newtonian return are new assumptions",
            "High-acceleration GR branch and quasistatic Galactic external field approximation",
            "Additive cosmological vacuum term remains freely adjustable"],
        "limitations":["No raw Cassini/MESSENGER tracking fit or joint covariance likelihood",
            "No global relativistic matching of the Solar System, moving Galactic environment and cosmological foliation",
            "Positivity, constraint-counting and principal 3D checks are not global nonlinear stability or quantum consistency proofs",
            "The higher-multipole diagnostic evaluates degrees 3 and 4 at all 17 sampled fields, with central refinement controls; higher degrees and a rigorous numerical error bound are not supplied",
            "The galaxy result is the same 2693-point full-sample spherical proxy, not held-out or resolved disk fits",
            "17 external magnitudes sample the Gaia two-sigma interval; they do not prove compatibility at every intermediate value",
            "The homogeneous equations and formal infinitesimal GR limit do not fit nonlinear cosmological perturbations, CMB or dark matter",
            "No quantum derivation or novelty claim for the known constrained-gravity ingredients"],
        "source_sha256":digest(__file__),"input_sha256":{str(p.relative_to(BASE)):digest(p) for p in
            [OUT/(key+"_summary.json") for key in reports]+[BASE/"OVERNIGHT_SOLAR_PROTOCOL.md",BASE/"CMC_METRIC_DERIVATION.md"]}}
    (OUT/"publication_summary.json").write_text(json.dumps(report,indent=2,allow_nan=False)+"\n")
    for name,values in [("external_fields.csv",rows),("solar_controls.csv",control_rows)]:
        with (OUT/name).open("w",newline="") as stream:
            writer=csv.DictWriter(stream,fieldnames=list(values[0]));writer.writeheader();writer.writerows(values)
    print(json.dumps({k:report[k] for k in ["status","central_solar_case","largest_extra_precession_envelope_arcsec_century","largest_equivalent_beta_bias"]},indent=2),flush=True)


if __name__=="__main__":
    main()
