#!/usr/bin/env python3
"""Exact designed checks for I05. All expectations are analytic, not fitted.
A check group may compare many components. Replays are not independent review.
"""
import argparse,json,platform
from pathlib import Path
from itertools import product,combinations
import sympy as s
from core import ETA,IDX,sm,simp,source_connection,geometry,connection,curvature,ricci,source_ricci
from matter import field_quantities,covariant_divergence,frame_divergence
from partd import *
checks=[];fixtures={}
def flatten(v):
    if isinstance(v,s.MatrixBase):return list(v)
    if isinstance(v,(list,tuple)):return [a for b in v for a in flatten(b)]
    return [v]
def zero(name,expr):
    vals=[simp(v) for v in flatten(expr)];bad=[str(v) for v in vals if v!=0]
    checks.append(dict(name=name,kind='exact_identity',passed=not bad,components=len(vals),nonzero=bad[:5]))
    print(('PASS ' if not bad else 'FAIL ')+name,flush=True)
def control(name,ok,witness=None):
    checks.append(dict(name=name,kind='expected_failure_control',passed=bool(ok),witness=str(witness) if witness is not None else None))
    print(('PASS ' if ok else 'FAIL ')+name,flush=True)
def notzero(v):return any(simp(a)!=0 for a in flatten(v))
x,y,z,t=s.symbols('x y z t',real=True);co=[x,y,z,t]
kappa=s.symbols('kappa',positive=True)
# Non-diagonal rational first jet, independent potential and field data.
E=s.Matrix([[2,1,0,0],[0,1,1,0],[0,0,3,1],[1,0,0,2]])
dE=[s.Matrix(4,4,lambda a,b:s.Rational(((m+2)*(a+3)+(b+1)**2)%13-6,7)) for m in IDX]
Fi,w,C=source_connection(E,dE);g=E.T*ETA*E;q=g.inv()
dg=[d.T*ETA*E+E.T*ETA*d for d in dE];dq=[-q*d*q for d in dg]
zero('baseline_connection_recovered_at_dense_first_jet',[C[m]-connection(q,dg)[m] for m in IDX])
phi=s.Matrix([2,-3,1,4]);dphi=[s.Matrix([s.Rational((m+1)*(a+2)-3,5) for a in IDX]) for m in IDX]
alpha=Fi.T*phi
dalpha=[Fi.T*dphi[m]-Fi.T*dE[m].T*alpha for m in IDX]
f=field_from_potential(alpha,dalpha,Fi,w)
Fc=s.Matrix(4,4,lambda a,b:dphi[b][a]-dphi[a][b])
zero('potential_field_operation_intertwines',E.T*f*E-Fc)
zero('field_antisymmetry',f+f.T)
zero('potential_reverse_map',E.T*alpha-phi)
# General raw antisymmetric field: NOT forced to equal a potential derivative.
Fc=s.Matrix(4,4,lambda a,b:(a-b)*(a+b+2))
dFc=[s.Matrix(4,4,lambda a,b:s.Rational((a-b)*((m+2)*(a+b+1)+a*b),3)) for m in IDX]
f=sm(Fi.T*Fc*Fi)
dFi=[-Fi*d*Fi for d in dE]
df=[sm(dFi[m].T*Fc*Fi+Fi.T*dFc[m]*Fi+Fi.T*Fc*dFi[m]) for m in IDX]
J=s.Matrix([2,-1,3,5]);j=E*J
fr=field_ops_frame(f,df,j,Fi,w);cr=field_ops_coord(Fc,dFc,J,g,q,dg,dq,C)
zero('current_reverse_readout',Fi*j-J)
zero('raised_field_operation_intertwines',Fi*ETA*f*ETA*Fi.T-cr['up'])
zero('maxwell_derivative_intertwines',Fi*fr['maxwell']-cr['maxwell'])
zero('homogeneous_derivative_intertwines',[sum(E[a,m]*E[b,n]*E[c,r]*fr['hom'][a][b][c] for a,b,c in product(IDX,repeat=3))-cr['hom'][m][n][r] for m,n,r in product(IDX,repeat=3)])
zero('em_stress_algebra_intertwines',E.T*fr['Sigma']*E-cr['T'])
zero('em_stress_symmetric',fr['Sigma']-fr['Sigma'].T)
zero('em_stress_trace_zero',s.trace(ETA*fr['Sigma']))
zero('em_scalar_invariant_intertwines',fr['invariant']-cr['invariant'])
zero('force_operation_intertwines',E.T*fr['force']-cr['force'])
zero('stress_derivative_intertwines',E.T*fr['div']-cr['div'])
zero('complete_raw_maxwell_stress_identity',fr['div']-fr['force']-fr['defect'])
control('discarding_raw_homogeneous_and_maxwell_defects_fails',notzero(fr['div']-fr['force']))
zero('coordinate_complete_raw_stress_identity',cr['div']-cr['force']-Fc*cr['maxwell']-s.Matrix([sum(cr['up'][a,b]*cr['hom'][mu][a][b]/2 for a,b in product(IDX,repeat=2)) for mu in IDX]))
# Lorentz/frame operations, fixed geometric preparation.
L=s.eye(4);L[0,0]=L[3,3]=s.Rational(5,4);L[0,3]=L[3,0]=s.Rational(3,4)
zero('declared_frame_change_is_lorentz',L.T*ETA*L-ETA)
fp=L.inv().T*f*L.inv();Ep=L*E
zero('frame_change_preserves_em_field',Ep.T*fp*Ep-Fc)
zero('frame_change_preserves_em_stress',Ep.T*em_stress(fp)[0]*Ep-cr['T'])
zero('frame_change_preserves_force',Ep.T*fp*(L*j)-cr['force'])
control('frame_change_of_geometry_alone_changes_field',notzero(Ep.T*f*Ep-Fc))
zero('field_sign_reversal_has_same_stress',em_stress(-f)[0]-em_stress(f)[0])
control('stress_does_not_determine_electromagnetic_field',notzero(f-(-f)))
# Fluid at a generic normalized jet; derivatives tangent to norm one.
v=s.Matrix([s.Rational(3,4),0,0,s.Rational(5,4)]);p=s.Rational(2,7);rho=s.Rational(11,5)
dv=[]
for m in IDX:
    a=s.Matrix([s.Rational((m+1)*(i+2)-4,9) for i in IDX]);dv.append(a-v*(v.T*ETA*a)[0])
dp=[s.Rational(m-2,5) for m in IDX];drho=[s.Rational(2*m+1,7) for m in IDX]
fl=fluid_differential(v,p,rho,dv,dp,drho,Fi,w)
zero('fluid_norm_and_tangent_coverage',[(v.T*ETA*v)[0]-1]+[(v.T*ETA*a)[0] for a in dv])
u=Fi*v;Tf=E.T*fl['Sigma']*E
zero('fluid_original_contravariant_stress',q*Tf*q-(-p*q+rho*u*u.T))
zero('fluid_original_mixed_stress',q*Tf-(-p*s.eye(4)+rho*u*(g*u).T))
zero('fluid_energy_eigenvalue',fl['Sigma']*v-(rho-p)*(ETA*v))
zero('fluid_exact_longitudinal_transverse_split',fl['div']-(ETA*v)*fl['energy']-fl['euler'])
zero('fluid_transverse_euler_has_three_independent_slots',(v.T*fl['euler'])[0])
zero('fluid_energy_projection',(v.T*fl['div'])[0]-fl['energy'])
dTf=[dE[m].T*fl['Sigma']*E+E.T*fl['dSigma'][m]*E+E.T*fl['Sigma']*dE[m] for m in IDX]
zero('fluid_divergence_intertwines',E.T*fl['div']-covariant_divergence(Tf,g,q,C,dTf,dq))
zero('fluid_stress_trace_convention',s.trace(ETA*fl['Sigma'])-(rho-4*p))
control('rho_is_not_rest_energy_for_nonzero_pressure',simp((rho-p)-rho)!=0)
zero('fluid_zero_rho_loses_velocity',fluid_stress(v,p,0)-fluid_stress(s.Matrix([0,0,0,1]),p,0))
control('fluid_velocity_not_recoverable_from_degenerate_stress',v!=s.Matrix([0,0,0,1]))
vbad=2*v;badfl=fluid_differential(vbad,p,rho,[2*a for a in dv],dp,drho,Fi,w)
control('normalization_cannot_be_removed_from_fluid_projection',notzero(badfl['div']-(ETA*vbad)*badfl['energy']-badfl['euler']))
# Changing local boost: a constant coordinate potential has exactly zero field.
Lb=s.eye(4);Lb[0,0]=Lb[3,3]=s.cosh(x);Lb[0,3]=Lb[3,0]=s.sinh(x)
Fb,wb,Cb=source_connection(Lb,[Lb.diff(c) for c in co]);phib=s.Matrix([0,0,0,1]);ab=sm(Fb.T*phib)
proper=field_from_potential(ab,[ab.diff(c) for c in co],Fb,wb)
bare=field_from_potential(ab,[ab.diff(c) for c in co],Fb,[s.zeros(4) for _ in IDX])
zero('varying_frame_potential_field_zero',proper)
control('bare_frame_derivatives_create_spurious_field',notzero(bare))
# Potential contact, gauge and explicit local reverse map on a star-shaped patch.
phi0=s.Matrix([x*t+y*z, x**2*z+t*y, y**2*t+x, t**2*x+y])
F0=coordinate_field(phi0,co);chi=x*y*z+t*x**2
zero('electromagnetic_gauge_invariance',coordinate_field(phi0+s.Matrix([s.diff(chi,c) for c in co]),co)-F0)
zero('potential_implies_homogeneous_equations',[s.diff(F0[b,c],co[a])+s.diff(F0[c,a],co[b])+s.diff(F0[a,b],co[c]) for a,b,c in product(IDX,repeat=3)])
lam=s.symbols('lam',real=True)
scaled={co[a]:lam*co[a] for a in IDX}
phirec=s.Matrix([-s.integrate(sum(lam*co[a]*F0[a,b].subs(scaled,simultaneous=True) for a in IDX),(lam,0,1)) for b in IDX])
zero('radial_homotopy_reconstructs_potential_field',coordinate_field(phirec,co)-F0)
zero('radial_gauge_section_condition',sum(co[a]*phirec[a] for a in IDX))
control('potential_reverse_section_not_identical_potential',notzero(phirec-phi0))
# Raw homogeneous defect cannot be deleted even with J=0 and vanishing Maxwell divergence.
rawF=s.zeros(4);rawF[0,1]=z;rawF[1,0]=-z
raw=field_ops_coord(rawF,[rawF.diff(c) for c in co],s.zeros(4,1),ETA,ETA,[s.zeros(4)]*4,[s.zeros(4)]*4,[s.zeros(4)]*4)
zero('raw_field_has_zero_maxwell_divergence',raw['maxwell'])
control('raw_nonclosed_field_is_not_potential',raw['hom'][2][0][1]!=0)
zero('homogeneous_defect_energy_term',raw['div']-s.Matrix([0,0,z,0]))
# Flat limit: all e,h components, source stress and Poynting readout.
e=s.Matrix(s.symbols('e1:4',real=True));h=s.Matrix(s.symbols('h1:4',real=True));Fe=s.zeros(4)
Fe[1,2]=h[0];Fe[2,0]=h[1];Fe[0,1]=h[2]
for a in range(3):Fe[a,3]=e[a]
Fe=Fe-Fe.T
Se,Ie=em_stress(Fe);mixed=ETA*Se;energy=(e.dot(e)+h.dot(h))/2;Poy=e.cross(h)
zero('flat_electromagnetic_energy',mixed[3,3]-energy)
zero('flat_poynting_flux',s.Matrix([mixed[i,3] for i in range(3)])-Poy)
zero('flat_spatial_stress',Se[:3,:3]-(s.eye(3)*energy-e*e.T-h*h.T))
zero('flat_em_invariant',Ie-2*(h.dot(h)-e.dot(e)))
jvec=s.Matrix(s.symbols('j1:4',real=True));charge=s.symbols('charge',real=True);xi=Fe*s.Matrix([*jvec,charge])
zero('original_force_spatial_components',xi[:3,0]-(charge*e+jvec.cross(h)))
zero('original_force_time_component',xi[3]+e.dot(jvec))
# Metric constitutive map must not be replaced by the flat one.
Ec=s.diag(2,3,1,s.Rational(1,6));gc=Ec.T*ETA*Ec;qc=gc.inv();upc=qc*Fe*qc
zero('nontrivial_unit_volume_metric',Ec.det()-1)
zero('nonflat_coordinate_constitutive_e',[upc[i,3]+[9,4,36][i]*e[i] for i in range(3)])
zero('nonflat_coordinate_constitutive_h',[upc[1,2]-h[0]/9,upc[2,0]-h[1]/4,upc[0,1]-h[2]/36])
control('flat_constitutive_identification_fails_in_general_chart',notzero(upc-ETA*Fe*ETA))
# Explicit electric source with nonzero energy-momentum transfer, fixed flat background.
phis=s.Matrix([0,0,0,-x**2/2]);Fs=coordinate_field(phis,co);Js=s.Matrix([0,0,0,1])
fsops=field_ops_coord(Fs,[Fs.diff(c) for c in co],Js,ETA,ETA,[s.zeros(4)]*4,[s.zeros(4)]*4,[s.zeros(4)]*4)
zero('charged_maxwell_equations_fixed_background',fsops['maxwell'])
zero('charged_example_closed_field',fsops['hom'])
zero('charged_example_nonzero_transfer',fsops['force']-s.Matrix([x,0,0,0]))
zero('charged_example_stress_balance',fsops['div']-fsops['force'])
control('maxwell_does_not_imply_zero_transfer',notzero(fsops['force']))
control('flat_nonzero_field_not_full_gravitational_solution',notzero(kappa*fsops['T']))
# Surface boundary operation F=-dphi, rectangle orientation dx wedge dt.
flux=s.integrate(Fs[0,3],(x,0,1),(t,0,1))
boundary=s.integrate(phis[3].subs(x,1),(t,0,1))-s.integrate(phis[3].subs(x,0),(t,0,1))
zero('electromagnetic_surface_boundary_sign',flux+boundary)
control('finite_surface_flux_not_zero',flux!=0)
# Variable-volume Maxwell residual; time-dependent spatial calibration.
a=1+t;Ev=s.diag(a,a,a,1);Fv,wv,Cv=source_connection(Ev,[Ev.diff(c) for c in co]);gv=Ev.T*ETA*Ev;qv=gv.inv();nuv=Ev.det()
phiv=s.Matrix([t**2/2,0,0,0]);Fcov=coordinate_field(phiv,co);Fup=sm(qv*Fcov*qv)
Jv=sm(s.Matrix([sum(s.diff(nuv*Fup[m,n],co[n]) for n in IDX)/nuv for m in IDX]))
vops=field_ops_coord(Fcov,[Fcov.diff(c) for c in co],Jv,gv,qv,[gv.diff(c) for c in co],[qv.diff(c) for c in co],Cv)
zero('volume_density_form_equals_covariant_maxwell',vops['maxwell'])
zero('charge_continuity_density',sum(s.diff(nuv*Jv[m],co[m]) for m in IDX))
naive=sm(s.Matrix([sum(s.diff(Fup[m,n],co[n]) for n in IDX) for m in IDX]))
control('omitting_volume_changes_maxwell_equation',notzero(naive-Jv))
# A non-solenoidal prescribed current is rejected rather than manufactured by the field.
Jbad=Jv+s.Matrix([x,0,0,0]);mbad=vops['maxwell']-s.Matrix([x,0,0,0])
zero('off_shell_charge_residual',sum(s.diff(nuv*(Jbad[m]+mbad[m]),co[m]) for m in IDX))
control('prescribed_nonconserved_current_not_maxwell_solution',simp(sum(s.diff(nuv*Jbad[m],co[m]) for m in IDX))!=0)
# Pressure-gradient/acceleration fluid on a fixed background; not claimed Einstein-coupled.
N=1+x;El=s.diag(1,1,1,N);Fl,wl,Cl=source_connection(El,[El.diff(c) for c in co]);vl=s.Matrix([0,0,0,1]);pl=2/N-1;rl=2/N
flap=fluid_differential(vl,pl,rl,[vl.diff(c) for c in co],[s.diff(pl,c) for c in co],[s.diff(rl,c) for c in co],Fl,wl)
zero('hydrostatic_fluid_energy_equation',flap['energy'])
zero('hydrostatic_fluid_euler_equation',flap['euler'])
zero('hydrostatic_fluid_full_divergence',flap['div'])
control('hydrostatic_test_has_nonzero_acceleration',notzero(flap['acc']))
control('hydrostatic_test_has_nonzero_pressure_gradient',notzero(flap['gradp']))
# I04 joint fluids now reconstructed from constitutive variables and EOS.
r=s.symbols('r',positive=True);crd=[x,y,z,r];vr=s.Matrix([0,0,0,1])
for name,Ed,eps,pd in [
 ('dust',s.diag(r**s.Rational(2,9),r**s.Rational(2,9),r**s.Rational(2,9),r**s.Rational(-2,3)),4*r**s.Rational(-2,3)/(27*kappa),s.S.Zero),
 ('radiation_fluid',s.diag(r**s.Rational(1,5),r**s.Rational(1,5),r**s.Rational(1,5),r**s.Rational(-3,5)),3*r**s.Rational(-4,5)/(25*kappa),r**s.Rational(-4,5)/(25*kappa))]:
    G=geometry(Ed,crd);rh=eps+pd
    FF=fluid_differential(vr,pd,rh,[vr.diff(c) for c in crd],[s.diff(pd,c) for c in crd],[s.diff(rh,c) for c in crd],G['F'],G['w'])
    MQ=field_quantities(G,FF['Sigma'],crd,kappa)
    zero(name+'_constitutive_stress',FF['Sigma']-s.diag(pd,pd,pd,eps))
    zero(name+'_energy_equation',FF['energy']);zero(name+'_euler_equation',FF['euler'])
    zero(name+'_same_realization_einstein_equation',MQ['Kcov'])
    zero(name+'_total_conservation',MQ['divtotal'])
    fixtures[name]={'E':str(Ed),'pressure':str(pd),'Einstein_rho':str(rh),'rest_energy':str(eps),'full_field_residual':'0','fluid_residual':'0'}
# A genuinely joined Einstein-Maxwell wave, not a flat test field plus separate gravity.
# Derived family: H=kappa A(u)^2 (x^2+y^2)/2, F=A(u) dx wedge du.
# Execute exact polynomial instance A(u)=1+u, kappa=2 to keep algebra small.
kw=s.Integer(2);u=t-z;Aw=1+u;Hw=kw*Aw**2*(x*x+y*y)/2
Ew=s.eye(4);Ew[2,2]=1-Hw/2;Ew[2,3]=Hw/2;Ew[3,2]=-Hw/2;Ew[3,3]=1+Hw/2
# Compute from the source structure-coefficient connection, not preloaded metric curvature.
Gw=geometry(Ew,co);phiw=s.Matrix([u+u*u/2,0,0,0]);Fw=coordinate_field(phiw,co)
alphaw=sm(Gw['F'].T*phiw)
fw=field_from_potential(alphaw,[alphaw.diff(c) for c in co],Gw['F'],Gw['w'])
EMw=field_ops_frame(fw,[fw.diff(c) for c in co],s.zeros(4,1),Gw['F'],Gw['w'])
Mw=field_quantities(Gw,EMw['Sigma'],co,kw)
nullcov=s.Matrix([0,0,-1,1])
zero('joined_em_wave_unit_volume',Ew.det()-1)
zero('joined_em_wave_potential_readout',Ew.T*fw*Ew-Fw)
zero('joined_em_wave_homogeneous_law',EMw['hom'])
zero('joined_em_wave_inhomogeneous_law',EMw['maxwell'])
zero('joined_em_wave_null_invariant',EMw['invariant'])
zero('joined_em_wave_stress_specified_before_curvature',Mw['T']-Aw**2*nullcov*nullcov.T)
zero('joined_em_wave_einstein_equation',Mw['Kcov'])
zero('joined_em_wave_ricci_family_formula',Gw['Ric']-kw*Aw**2*nullcov*nullcov.T)
zero('joined_em_wave_source_stress_conservation',EMw['div'])
zero('joined_em_wave_total_current_conservation',Mw['divtotal'])
control('joined_em_wave_has_nonzero_curvature',notzero(Gw['Ric']))
control('joined_em_wave_has_nonzero_field',notzero(fw))
zero('joined_em_wave_wrong_source_sign_is_detectable_by_field',em_stress(-fw)[0]-EMw['Sigma'])
# Verify curvature again using direct coordinate Christoffels, independent constructor.
Cdirect=connection(Gw['q'],[Gw['g'].diff(c) for c in co])
zero('joined_em_wave_two_connection_constructions',[Cdirect[m]-Gw['C'][m] for m in IDX])
zero('joined_em_wave_two_curvature_routes',ricci(curvature(Cdirect,co))-Gw['Ric'])
# Boundary balance for EM stress on this same joint solution: unit spatial cube.
Q=((1+t)**3-t**3)/3
outflow=simp(Mw['mixed'][2,3].subs(z,1)-Mw['mixed'][2,3].subs(z,0))
primitive=-(1+t-z)**3/3
zero('joined_em_wave_boundary_density_primitive',s.diff(primitive,z)-Mw['mixed'][3,3])
zero('joined_em_wave_boundary_integral',primitive.subs(z,1)-primitive.subs(z,0)-Q)
zero('joined_em_wave_finite_flux_balance',s.diff(Q,t)+outflow)
control('joined_em_wave_boundary_flux_not_zero',outflow!=0)
# Wrong field amplitude on the same geometry preserves Maxwell but not gravity.
wrong=field_quantities(Gw,em_stress(2*fw)[0],co,kw)
control('maxwell_solution_on_same_geometry_not_automatically_joined',notzero(wrong['Kcov']))
fixtures['Einstein_Maxwell_wave']={'kappa':str(kw),'u':'t-z','amplitude':str(Aw),'H':str(Hw),'E':str(Ew),'potential':str(phiw),'field':str(Fw),'stress_cov':str(Mw['T']),'Ricci':str(Gw['Ric']),'Einstein_residual':'0','Maxwell_residual':'0','source_stress_divergence':'0','gravitational_current':str(Mw['t']),'Q4_unit_cube':str(Q),'outward_flux':str(outflow),'note':'Potential, amplitude and geometry supplied explicitly before residual evaluation; not T:=G/kappa.'}
fixtures['fixed_background_charge']={'phi4':'-x^2/2','F14':'x','J4':'1','transfer':'(x,0,0,0)','scope':'Exact Maxwell/stress witness; not claimed to solve coupled Einstein equation; does not model trajectories of free charges.'}
fixtures['raw_homogeneous_defect']={'F12':'z','J':'0','hom_312':'1','stress_divergence':'(0,0,z,0)'}
fixtures['hydrostatic_fixed_background']={'N':'1+x','pressure':'2/N-1','Einstein_rho':'2/N','rest_energy':'1','domain':'0<x<1 for positive pressure','scope':'fluid operation test, not a joint gravitational solution'}
parser=argparse.ArgumentParser();parser.add_argument('--output',type=Path,default=Path(__file__).with_name('results.json'));args=parser.parse_args()
result={'checkpoint':'E1916-I05-v1.0','environment':{'python':platform.python_version(),'sympy':s.__version__},'summary':{'groups':len(checks),'passed':sum(c['passed'] for c in checks),'failed':sum(not c['passed'] for c in checks),'expected_failure_controls':sum(c['kind']=='expected_failure_control' for c in checks)},'checks':checks,'fixtures':fixtures}
args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(json.dumps(result,indent=2,ensure_ascii=False)+'\n')
print(json.dumps(result['summary']),flush=True)
raise SystemExit(0 if all(c['passed'] for c in checks) else 1)
