"""Independent checks, mandatory before data publication."""
import sys,json
from pathlib import Path
sys.path.insert(0,str(Path(__file__).resolve().parents[1]/'src'))
from rod import Rod,Controls
import numpy as np
from scipy.integrate import solve_bvp
from scipy.optimize import brentq

def run():
    checks={}
    r=Rod(24); c=Controls(chi=.004,fx=.005,fy=.002,h=.001,k1=.3,k2=-.2,k3=.6)
    z=.4*np.sin(np.linspace(0,np.pi/2,r.d)); e,g,K=r.terms(z,c)
    eps=1e-6; eye=np.eye(r.d)
    gn=np.array([(r.terms(z+eps*v,c)[0]-r.terms(z-eps*v,c)[0])/(2*eps) for v in eye])
    Kn=np.column_stack([(r.terms(z+eps*v,c)[1]-r.terms(z-eps*v,c)[1])/(2*eps) for v in eye])
    checks['gradient_fd_max_abs']=float(np.max(np.abs(g-gn)))
    checks['hessian_fd_max_abs']=float(np.max(np.abs(K-Kn)))
    assert checks['gradient_fd_max_abs']<1e-8 and checks['hessian_fd_max_abs']<1e-8
    axial=r.solve(Controls(fx=.03))
    checks['axial_tip_error']=float(np.linalg.norm(axial['tip']-[1.03,0]))
    assert checks['axial_tip_error']<1e-12
    circular=[]
    for n in [24,48,96]:
        rr=Rod(n); cc=Controls(k1=1.2,k2=1.2,k3=1.2)
        ans=rr.solve(cc,np.linspace(0,1.2,n+1)[1:])
        exact=np.array([np.sin(1.2)/1.2,(1-np.cos(1.2))/1.2])
        circular.append(float(np.linalg.norm(ans['tip']-exact)))
    checks['circular_tip_errors_n24_48_96']=circular
    assert circular[2]<circular[1]<circular[0] and circular[-1]<1e-5
    rr=Rod(96); cc=Controls(chi=1e-4,h=.01,delta=.1)
    pc=(-1+np.sqrt(1+16*cc.h*cc.delta-np.pi**2*cc.chi))/2-2*cc.h
    def eig(p):
        cc.fx=p
        return rr.solve(cc)['lambda_min_over_chi']
    pn=brentq(eig,pc-.001,pc+.001,xtol=1e-13)
    checks['critical_load_analytic']=pc; checks['critical_load_discrete_rod96']=pn
    checks['critical_load_error']=abs(pc-pn)
    assert abs(pc-pn)<1e-7
    # Independent solve_bvp of paper Eqs.52-53, negative branch. Not our FE energy.
    errors=[]; chain_errors=[]
    for p in [-.017159,-.016428,-.016316]:
        chi,h,d=1e-4,.01,.1; a=p+2*h
        s=np.linspace(0,1,101)
        def fun(s,y):
            th=y[0]; stretch=1+a*np.cos(th)
            return np.vstack([y[1],(a*stretch-4*h*d)*np.sin(th)/chi,stretch*np.cos(th),stretch*np.sin(th)])
        def bc(l,r): return [l[0],r[1],l[2],l[3]]
        candidates=[]
        for amp in [-.3,-1.,-2.]:
            y=np.vstack([amp*np.sin(np.pi*s/2),amp*np.pi/2*np.cos(np.pi*s/2),s,amp*s*s/2])
            b=solve_bvp(fun,bc,s,y,tol=1e-9,max_nodes=10000)
            if b.success and b.y[3,-1]<-1e-3: candidates.append(b)
        assert candidates
        b=min(candidates,key=lambda b:abs(b.y[0,-1]))
        target=b.sol(np.linspace(0,1,65))[2:4].T
        c=Controls(chi=chi,h=h,delta=d,fx=p)
        arod=Rod(96); out=arod.solve(c,b.sol(np.linspace(0,1,97))[0,1:])
        assert out['accepted'] and out['tip'][1]<0 and out['stability']=='stable'
        errors.append(float(np.sqrt(np.mean(np.sum((out['centerline']-target)**2,axis=1)))))
        one=[]
        for n in [15,30,60,120]:
            rc=Rod(n,'chain'); v=rc.solve(c,b.sol(np.arange(1,n)/n)[0])
            assert v['accepted'] and v['tip'][1]<0
            one.append(float(np.linalg.norm(v['tip']-target[-1])))
        assert one[-1]<one[0]
        chain_errors.append(one)
    checks['fig9_rod96_centerline_rmse_over_L']=errors
    checks['fig9_chain_tip_errors_n15_30_60_120']=chain_errors
    assert max(errors)<.001
    # Exact chain linear critical load verifies the independently assembled hinge boundary.
    rc=Rod(15,'chain'); c=Controls(chi=1e-4,h=.01)
    omega2=4*np.sin(np.pi/(2*(2*15-1)))**2
    expected=(-1+np.sqrt(1-4*c.chi*15**2*omega2+16*c.h*c.delta))/2-2*c.h
    c.fx=expected; v=rc.solve(c)
    checks['chain_critical_eigenvalue']=v['lambda_min_over_chi']
    assert abs(v['lambda_min_over_chi'])<1e-8
    checks['passed']=True
    return checks

if __name__=='__main__':
    checks=run()
    p=Path(__file__).resolve().parents[1]/'reports/validation.json'
    p.write_text(json.dumps(checks,indent=2))
    print(json.dumps(checks,indent=2))