"""Actual-header CPU algebra, not CUDA arithmetic/performance evidence."""
import numpy as np
from .test_numerical_operators import compile_run


def test_retained_density_and_shifted_planar_headers(tmp_path):
    rows=compile_run(tmp_path,'shifted_block_operators.cpp',{})
    cx=np.array([0,1,-1,0,0,1,-1,1,-1]);cy=np.array([0,0,0,1,-1,1,-1,-1,1])
    w=np.array([4/9]+[1/9]*4+[1/36]*4)
    it=iter(rows)
    for header in it:
        kind,*numbers=header.split()
        a=np.fromstring(next(it),sep=' ',dtype=np.float32).astype(float)
        if kind=='E':
            delta,u,v=np.array(numbers,dtype=np.float32).astype(float)
            rho=float(np.float32(1+delta));cu=3*(cx*u+cy*v)
            expected=w*(delta+rho*(cu+.5*cu*cu-1.5*(u*u+v*v)))
            np.testing.assert_allclose(a,expected,rtol=3e-6,atol=1e-10)
            if u==0:
                np.testing.assert_allclose(a.sum(),delta,rtol=5e-8,atol=0)
                if delta: assert np.any(a!=0)
        elif kind=='J':
            delta,jx,jy,ux,uy=np.array(numbers,dtype=np.float32).astype(float)
            jd=cx*jx+cy*jy;ud=cx*ux+cy*uy
            expected=w*(delta+3*jd+4.5*jd*ud-1.5*(jx*ux+jy*uy))
            np.testing.assert_allclose(a,expected,rtol=3e-6,atol=2e-9)
            actual=np.array([a.sum(),a@cx,a@cy,a@(cx*cx),a@(cx*cy),a@(cy*cy)])
            target=np.array([delta,jx,jy,delta/3+jx*ux,jx*uy,delta/3+jy*uy])
            np.testing.assert_allclose(actual,target,atol=4e-9,rtol=3e-6)
            if jx==0 and delta: assert np.any(a!=0)
        else:
            mode=int(numbers[0]);rho,u,v=np.array(numbers[1:],dtype=np.float32).astype(float)
            b=np.fromstring(next(it),sep=' ',dtype=np.float32).astype(float)
            known=cy!=(1 if mode==0 else -1)
            np.testing.assert_array_equal(a[known],b[known])
            physical_rho=1+b.sum()
            if mode==2:
                assert abs(physical_rho-rho)<2e-8
                assert abs(b@cx)<2e-9
            else:
                np.testing.assert_allclose([b@cx,b@cy],physical_rho*np.array([u,v]),atol=3e-9,rtol=0)