"""Independent spectrum, quantum geometry, static response and descriptive RAR checks."""
from __future__ import annotations
import hashlib
import json
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.interpolate import CubicHermiteSpline, PchipInterpolator
from scipy.optimize import minimize_scalar
from dark_band_probe import matrices, spectral, lattice_polar
from research import load_rar, empirical_rar

BASE=Path(__file__).resolve().parent
OUT=BASE/"outputs"


class Response:
    def __init__(self, rows, control):
        self.rows=rows
        self.a2=float(control["a2"])
        self.slope=float(control["deep_slope"])
        self.d=float(control["tail_D"])
        self.normalization=3/(2*self.a2)
        self.derivative=CubicHermiteSpline(rows["x"],rows["first"],rows["second"])
        self.integral=self.derivative.antiderivative()
        self.double=self.derivative.antiderivative(2)
        self.c=self.tail_coefficient(100.)

    def moment(self,x):
        return x*self.integral(x)-self.double(x)+self.double(0.)

    def tail_coefficient(self,limit):
        return float(self.normalization*(-self.moment(limit)+2*self.d/limit))

    def evaluate(self,x):
        x=np.asarray(x,dtype=float)
        core=np.clip(x,1e-4,100.)
        moment=self.moment(core)
        mu=1+self.normalization*moment/core**3
        radial=1+self.normalization*(self.derivative(core)/core-2*moment/core**3)
        high_mu=1-self.c/np.maximum(x,1e-4)**3+2*self.normalization*self.d/np.maximum(x,1e-4)**4
        high_radial=1+2*self.c/np.maximum(x,1e-4)**3-6*self.normalization*self.d/np.maximum(x,1e-4)**4
        mu=np.where(x<1e-4,self.slope*x,np.where(x>100,high_mu,mu))
        radial=np.where(x<1e-4,2*self.slope*x,np.where(x>100,high_radial,radial))
        return mu,radial

    def inverse(self):
        x=np.geomspace(1e-10,1e10,12001)
        mu,radial=self.evaluate(x)
        if np.any(mu<=0) or np.any(radial<=0) or np.any(np.diff(x*mu)<=0):
            return None
        return PchipInterpolator(np.log10(x*mu),np.log10(x))


def physical_ground(k,x,w,m=1.,omega=20.):
    h,*_=matrices(k,x,m,omega,w)
    _,v=np.linalg.eigh(h)
    state=v[:,0].copy()
    r=1/np.sqrt(1+x*x)
    f=x*r
    state[2:]=np.array([[f,r],[r,-f]])@state[2:]
    return state


def main():
    spectra=np.genfromtxt(OUT/"dark_spectra.csv",delimiter=",",names=True)
    bands=np.genfromtxt(OUT/"dark_lattice_bands.csv",delimiter=",",names=True)
    controls=np.genfromtxt(OUT/"dark_lattice_controls.csv",delimiter=",",names=True)
    checks=np.genfromtxt(OUT/"dark_band_checks.csv",delimiter=",",names=True)
    assert len(spectra)==108 and len(bands)==1456 and len(checks)==64
    spectral_checks=[]
    for row in spectra:
        result=spectral(row["k"],row["x"],row["mass"],row["omega"],row["wilson_mass"])
        errors={"eigenvalues":float(np.max(np.abs(result[0]-[row[f"e{i}"] for i in range(4)])))}
        for key,value in zip(("first","second","contact","response","fidelity"),result[1:]):
            errors[key]=float(abs(value-row[key]))
        assert max(errors.values())<1e-8,(tuple(row),errors)
        spectral_checks.append({"omega":float(row["omega"]),"k":float(row["k"]),"x":float(row["x"]),
                                "wilson_mass":float(row["wilson_mass"]),"errors":errors})
    fidelity_checks=[]
    for k in [-.7,.3]:
        for x in [.1,1.,3.]:
            delta=1e-4
            a=physical_ground(k,x-delta/2,.3)
            b=physical_ground(k,x+delta/2,.3)
            measured=(1-abs(np.vdot(a,b))**2)/delta**2
            expected=float(spectral(k,x,1.,20.,.3)[-1])
            assert abs(measured-expected)<3e-6
            fidelity_checks.append({"k":k,"x":x,"metric":expected,"overlap_estimate":measured,
                                    "absolute_error":abs(measured-expected)})
    print("Full complex spectra, physical source derivatives and fidelity overlap checks passed.",flush=True)

    integral_checks=[]
    for row in checks:
        # Independent raw polar quadrature; Rust integrates the IR cone analytically.
        expected=lattice_polar(row["x"],1.,20.,row["wilson"],128,48)
        error=np.abs(expected[:3]-np.array([row["energy_refined"],row["first_refined"],row["second_refined"]]))
        assert max(error)<2e-7,(tuple(row),error)
        refinement=np.abs(np.array([row["energy"],row["first"],row["second"]])-
                          np.array([row["energy_refined"],row["first_refined"],row["second_refined"]]))
        assert max(refinement)<2e-7,(tuple(row),refinement)
        integral_checks.append({"wilson":float(row["wilson"]),"x":float(row["x"]),
                                "independent_errors":error.tolist(),"Rust_refinement_errors":refinement.tolist()})
    print("Independent lattice-zone integrals passed.",flush=True)

    # A sufficient analytic filling bound, uniform in source x and in the zone.
    omega,m=20.,1.
    mu_d=-2*m*m/(omega+np.sqrt(omega*omega+4*m*m))
    kmax,wmax=np.sqrt(2.),2.
    a=omega-kmax-mu_d
    bound=a*(-omega*mu_d)-wmax**4/(4*a)
    assert bound>0 and kmax*kmax+wmax*wmax<omega*omega
    assert np.all(bands["occupied_margin"]>0) and np.all(bands["empty_margin"]>0)
    assert controls[-1]["a2"]<0

    response_data=[]
    summaries=[]
    responses={}
    xgrid=np.geomspace(1e-4,1e4,1201)
    for control in controls[:-1]:
        w=float(control["wilson"])
        rows=bands[np.isclose(bands["wilson"],w)]
        response=Response(rows,control)
        responses[w]=response
        mu,radial=response.evaluate(xgrid)
        # Coarsen source sampling without changing the microscopic quadrature.
        indices=np.r_[0,np.arange(1,len(rows),2),len(rows)-1]
        indices=np.unique(indices)
        coarse=Response(rows[indices],control)
        cmu,cradial=coarse.evaluate(xgrid)
        source_errors={"mu":float(np.max(np.abs(cmu-mu))),"radial":float(np.max(np.abs(cradial-radial)))}
        assert max(source_errors.values())<2e-4,(w,source_errors)
        energy_integral=response.integral(rows["x"])-response.integral(0.)
        energy_direct=rows["energy"]-rows["energy"][0]
        summary={"wilson":w,"a2":response.a2,"deep_mu_slope":response.slope,
                 "tail_cubed_coefficient":response.c,"tail_coefficient_limit30_difference":abs(response.c-response.tail_coefficient(30.)),
                 "min_mu_on_scan":float(mu.min()),"min_radial_on_scan":float(radial.min()),
                 "max_mu_on_scan":float(mu.max()),"source_grid_check":source_errors,
                 "integrated_derivative_energy_max_error":float(np.max(np.abs(energy_integral-energy_direct))),
                 "static_positive_on_scan":bool(np.all(mu>0)&np.all(radial>0))}
        summaries.append(summary)
        response_data.extend(zip(np.full_like(xgrid,w),xgrid,mu,radial))
    np.savetxt(OUT/"dark_response.csv",response_data,delimiter=",",header="wilson,x,mu,radial",comments="")

    data,checksum=load_rar()
    logbar,_,logobs,_=data.T
    fits=[]
    for w,response in responses.items():
        inverse=response.inverse()
        if inverse is None:
            fits.append({"wilson":w,"status":"not fitted: static positivity or invertibility failed"})
            continue
        def objective(logscale):
            prediction=logscale+inverse(logbar-logscale)
            return float(np.mean((prediction-logobs)**2))
        fit=minimize_scalar(objective,bounds=(-12,-8),method="bounded",options={"xatol":1e-10})
        assert fit.success
        scale=10**float(fit.x)
        fits.append({"wilson":w,"g_dark_m_s2":scale,"deep_a0_m_s2":scale/response.slope,
                     "rms_dex":float(np.sqrt(fit.fun)),"fitted_scales_per_fixed_wilson":1,
                     "scope":"Descriptive full-sample spherical RAR proxy; Wilson choice is additional family selection"})
    empirical=minimize_scalar(lambda z:float(np.mean((np.log10(empirical_rar(10**logbar,10**z))-logobs)**2)),
                              bounds=(-12,-8),method="bounded",options={"xatol":1e-10})
    assert empirical.success
    comparison={"rows":len(data),"source_sha256":checksum,"fits":fits,
                "empirical_RAR_one_scale_rms_dex":float(np.sqrt(empirical.fun)),
                "limitations":"No held-out galaxy validation, no disc PDE, correlated radii, no dark-halo fit, and no solar prediction in this calculation"}
    report={"status":"LQ8 finite-band nonlinear mechanism: lattice spectrum, quantum geometry and static response checked; not a complete gravity theory",
            "rust_source_sha256":hashlib.sha256((BASE/"rust/src/bin/dark-bands.rs").read_bytes()).hexdigest(),
            "python_sources_sha256":{name:hashlib.sha256((BASE/name).read_bytes()).hexdigest()
                                     for name in ("dark_band_report.py","dark_band_probe.py")},
            "spectral_checks":spectral_checks,"fidelity_overlap_checks":fidelity_checks,
            "lattice_integral_checks":integral_checks,"response_summaries":summaries,"RAR":comparison,
            "fixed_filling":{"chemical_potential":mu_d,"analytic_positive_bound":bound,
                              "occupied_margin_min_on_quadrature":float(bands["occupied_margin"].min()),
                              "empty_margin_min_on_quadrature":float(bands["empty_margin"].min())},
            "negative_stiffness_control":{"mass":.5,"omega":20.,"wilson":.5,"a2":float(controls[-1]["a2"])},
            "limitations":["Quadratic source dependence and critical gravitational compensation are chosen microscopic inputs",
                           "The free-band projection is not a symmetry protection against all interactions",
                           "The lattice is a physical model choice, not a derived relativistic spacetime",
                           "Band filling carries matter and stress; the vacuum energy problem is not solved",
                           "At source zero the integrated fidelity has a logarithmic IR divergence; its finite quadrature value is not a physical prediction",
                           "No full coupling to a dynamical metric, wall network or original microscopic bulk theory is derived",
                           "No RAR score establishes novelty, quantum gravity or agreement with solar and relativistic tests"]}
    (OUT/"dark_summary.json").write_text(json.dumps(report,indent=2,allow_nan=False)+"\n")

    fig,axes=plt.subplots(1,2,figsize=(11,4.2),layout="constrained")
    for summary in summaries:
        w=summary["wilson"]
        if w not in [.25,.4,.45,.5]:continue
        mu,radial=responses[w].evaluate(xgrid)
        axes[0].semilogx(xgrid,mu,label=f"r_W={w:g}")
        axes[1].semilogx(xgrid,radial,label=f"r_W={w:g}")
    axes[0].axhline(1,color="gray",ls="--",lw=.8)
    axes[1].axhline(0,color="gray",lw=.8)
    axes[0].set(xlabel="Accélération / g_dark",ylabel="Réponse μ",xlim=(1e-3,100),title="Même Hamiltonien, différents paramètres de réseau")
    axes[1].set(xlabel="Accélération / g_dark",ylabel="Raideur radiale μ + x μ′",xlim=(1e-3,100),title="Contrôle de la convexité statique")
    for ax in axes:ax.grid(alpha=.2);ax.legend(fontsize=8)
    fig.supxlabel("M=1, Ω=20 ; compensation critique réglée ; ni fond gravitationnel ni modèle cosmologique validé.",fontsize=9)
    for suffix in ("png","pdf"):fig.savefig(OUT/f"dark_response.{suffix}",dpi=170)
    plt.close(fig)
    print(json.dumps({"responses":summaries,"RAR":comparison},indent=2))


if __name__=="__main__":main()
