"""V3 classical tests + actual DESI DR2 BAO-only Gaussian fit.
No CMB/SNe likelihood, laboratory likelihood, or UV completion is claimed.
"""
from pathlib import Path
import json,hashlib,platform
import numpy as np
import scipy
from scipy.integrate import quad
from scipy.optimize import least_squares
from scipy.linalg import solve_triangular
from scipy.stats import chi2
import sympy as s
R=Path(__file__).resolve().parent;D=R/'data'
x,B,Delta,rho=s.symbols('x B Delta rho',real=True)
A=(1-x*x)**2*(1+3*x*x)
V=B*x*x*(1-x*x)**2+Delta*x*x*(2-x*x)
W=V+rho*A
checks={
'A0':A.subs(x,0)==1,'A1':A.subs(x,1)==0,
'Aprime0':s.diff(A,x).subs(x,0)==0,'Asecond0':s.diff(A,x,2).subs(x,0)==2,
'Wprime_factor':s.simplify(s.diff(W,x)-2*x*(1-x*x)*(B*(1-3*x*x)+2*Delta+rho*(1-9*x*x)))==0,
'curvature_ordinary':s.simplify(s.diff(W,x,2).subs(x,0)-(2*B+4*Delta+2*rho))==0,
'curvature_zero':s.simplify(s.diff(W,x,2).subs(x,1)-8*(B-Delta+4*rho))==0,
'energy_gap':s.simplify(W.subs(x,1)-W.subs(x,0)-(Delta-rho))==0,
'origin_even':s.simplify(A.subs(x,-x)-A)==0 and s.simplify(V.subs(x,-x)-V)==0}
assert all(checks.values())
rows=[]
for density in [0,.001,.1,1,10,100]:
 saddle=np.sqrt((10+2+density)/(30+9*density))
 for state,a in [('ordinary',0),('barrier',saddle),('zero',1)]:
  vv=float(W.subs({x:a,B:10,Delta:1,rho:density}));kk=float(s.diff(W,x,2).subs({x:a,B:10,Delta:1,rho:density}))
  rows.append({'rho_over_Delta':density,'state':state,'chi':float(a),'mass_ratio':float(A.subs(x,a)),'energy_over_Delta':vv,'curvature_over_Delta':kk})
# Natural-unit conversion: 1 eV^4 = eV_J/(hbar*c in eV*m)^3 J/m^3.
c=299792458.;evJ=1.602176634e-19;hbar_ev_s=6.582119569e-16;hc=1.973269804e-7
conv=evJ/hc**3;Mpl=2.435e27;f=Mpl;dv=1e36;bv=10*dv
m=np.sqrt(2*bv+4*dv)/f
H0=67.4*1000/3.085677581491367e22;H0eV=H0*hbar_ev_s
crit=3*Mpl*Mpl*H0eV*H0eV;Om=.315;Or=9.2e-5;Ol=1-Om-Or
Omega_spec=1e-6;spec=Omega_spec*crit
xamp=np.sqrt(spec/(bv+2*dv));aamp=2*Mpl/f*xamp
benchmark={'f_eV':f,'Delta_eV4':dv,'B_eV4':bv,'rhoLambda_eV4':Ol*crit,'scalar_mass_eV':m,'range_m':hc/m,'Delta_J_m3':dv*conv,'global_stability_density_kg_m3':dv*conv/c**2,'vacuum_barrier_J_m3':2.08*dv*conv,'critical_density_eV4':crit,'illustrative_Omega_scalar':Omega_spec,'homogeneous_scalar_amplitude_today':xamp,'max_linear_alpha_amplitude':aamp,'max_weak_source_fifth_force_ratio_no_Yukawa':2*aamp*aamp,'long_range_PPN_gamma_deviation_envelope':4*aamp*aamp,'mass_shift_amplitude_approx':xamp*xamp,'zero_phase_deSitter_radius_m':hc*np.sqrt(3)*Mpl/np.sqrt(dv)}
# Actual Gaussian BAO likelihood using published covariance, preserving row order.
raw=np.loadtxt(D/'desi_gaussian_bao_ALL_GCcomb_mean.txt',dtype=str)
z=raw[:,0].astype(float);data=raw[:,1].astype(float);kinds=raw[:,2]
C=np.loadtxt(D/'desi_gaussian_bao_ALL_GCcomb_cov.txt');assert C.shape==(13,13)
assert np.allclose(C,C.T);Ch=np.linalg.cholesky(C)
def predict(params):
 om,Q=params
 def E(zv):return np.sqrt(om*(1+zv)**3+Or*(1+zv)**4+1-om-Or)
 res=[]
 for zz,kind in zip(z,kinds):
  DM=Q*quad(lambda zz1:1/E(zz1),0,float(zz),epsabs=1e-11,epsrel=1e-11)[0];DH=Q/E(zz)
  if kind=='DV_over_rs':val=(zz*DM*DM*DH)**(1/3)
  elif kind=='DM_over_rs':val=DM
  elif kind=='DH_over_rs':val=DH
  else:raise ValueError(kind)
  res.append(val)
 return np.array(res)
def resid(params):return solve_triangular(Ch,predict(params)-data,lower=True)
sol=least_squares(resid,[.3,30],bounds=([.05,15],[.6,50]),xtol=1e-13,ftol=1e-13,gtol=1e-10)
assert sol.success
sol2=least_squares(resid,[.4,28],bounds=([.05,15],[.6,50]),xtol=1e-13,ftol=1e-13,gtol=1e-10)
assert np.max(abs(sol.x-sol2.x))<1e-6
stat=float(sol.fun@sol.fun);covpar=np.linalg.inv(sol.jac.T@sol.jac)
reference=np.array([.315,c/1000/(67.4*147.09)])
fit={'scope':'flat LambdaCDM branch chi=0; BAO only, not full cosmology likelihood','N':13,'parameters':['Omega_m','Q=c/(H0*r_d)'],'Omega_r_fixed':Or,'bestfit':sol.x.tolist(),'linearized_errors':np.sqrt(np.diag(covpar)).tolist(),'chi2_min':stat,'dof':11,'p_value_approx':float(chi2.sf(stat,11)),'H0_times_rd_km_s':float(c/1000/sol.x[1]),'illustrative_Planck_like_input':{'Omega_m':.315,'H0_km_s_Mpc':67.4,'rd_Mpc':147.09},'chi2_at_reference':float(resid(reference)@resid(reference)),'covariance_parameter':covpar.tolist(),'reference_not_joint_CMB_fit':True}
np.savetxt(D/'bao_predictions.csv',np.c_[z,data,predict(sol.x),np.sqrt(np.diag(C)),predict(reference)],delimiter=',',header='z,observed,bestfit,sigma_diagonal,reference',comments='')
# Spectator-envelope background check, adjusted Lambda so H0 remains fixed.
a=np.geomspace(1e-9,1,1000);E2=Or/a**4+Om/a**3+Ol
E2s=Or/a**4+(Om+Omega_spec)/a**3+(Ol-Omega_spec)
bg={'assumption':'small-amplitude fast oscillation averaged as dust; no direct stiff-field integration','max_fractional_H_change_a_1e-9_to_1':float(np.max(abs(np.sqrt(E2s/E2)-1))),'note':'reference cosmology uses standard CDM; scalar is not required to explain all DM'}
results={'symbolic':checks,'stationary_scan':rows,'benchmark':benchmark,'BAO_fit':fit,'spectator_background':bg,'software':{'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__,'sympy':s.__version__},'script_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest()}
(D/'results.json').write_text(json.dumps(results,ensure_ascii=False,indent=2))
print(json.dumps({'symbolic':checks,'benchmark':benchmark,'BAO_fit':fit,'spectator_background':bg},ensure_ascii=False,indent=2))