"""Independent analytical, derivative, energy and time-refinement checks."""
from dataclasses import replace
import numpy as np
from scipy.integrate import trapezoid
from model import Skid,frequency_response,stiffness_derivative,simulate,floor_bump
from study import candidate,transient_metrics,local_peak_refinement


def numerical_checks(c,selected,nominal,corners):
    skid=candidate(c,selected['stiffness_N_m'],selected['damping_ratio']);checks=[]
    def add(id,description,expected,observed,passed,method):
        checks.append(dict(id=id,description=description,expected=expected,observed=float(observed),passed=bool(passed),method=method))
    for group in [0,1]:
        errors=[]
        for f in [2.,6.,8.,20.,35.]:
            h=skid.stiffness[group]*1e-5;k=list(skid.stiffness);k[group]+=h;plus=replace(skid,stiffness=tuple(k));k[group]-=2*h;minus=replace(skid,stiffness=tuple(k))
            finite=(frequency_response(plus,[f])['q'][0]-frequency_response(minus,[f])['q'][0])/(2*h)
            analytic=stiffness_derivative(skid,f,group)
            errors.append(np.linalg.norm(analytic-finite)/np.linalg.norm(analytic))
        add('derivative-'+str(group),'Complex FRF derivative for support group '+str(group+1),'relative norm error < 1e-6',max(errors),max(errors)<1e-6,'Implicit differentiation compared with centered finite differences at five frequencies, damping fixed.')
    f=27.;w=2*np.pi*f;F=40.;q=frequency_response(skid,[f],force=F)['q'][0]
    times=np.linspace(0,1,2001);initial=np.r_[q.real,(1j*w*q).real]
    harmonic=simulate(skid,times,force=lambda t:F*np.cos(w*t),initial=initial,rtol=1e-10,max_step=.001)
    exact=np.real(np.exp(1j*w*times[:,None])*q[None,:]);error=np.max(abs(exact-harmonic['q']))/np.max(abs(exact))
    add('harmonic-time','Time-domain solution agrees with complex frequency response','relative displacement error < 1e-7',error,error<1e-7,'Exact steady-state initial condition; DOP853 response versus Re(Q exp(i omega t)).')
    symmetric=Skid(100,16,(-.4,.4),(100000,100000),(600,600),0,.2);freqs=np.linspace(1,40,200);r=frequency_response(symmetric,freqs,force=15);ww=2*np.pi*freqs;exact=15/(200000-100*ww**2+1j*1200*ww)
    err=np.max(abs(r['q'][:,0]-exact))/np.max(abs(exact))
    add('sdof-limit','Symmetric centered forcing agrees with the scalar oscillator','relative response error < 1e-12',err,err<1e-12,'Independent closed-form scalar denominator k-m omega^2+i c omega.')
    decay=simulate(skid,np.linspace(0,3,6001),initial=[.001,.001,0,0],rtol=1e-10,max_step=.002)
    residual=np.max(abs(decay['energy_residual']))/decay['energy'][0]
    add('free-energy','Initial energy becomes viscous heat','relative energy balance < 1e-7',residual,residual<1e-7,'Integrated work and dissipation states versus kinetic plus spring energy.')
    loss=np.einsum('ni,ij,nj->n',decay['velocity'],skid.C,decay['velocity']);quadrature=trapezoid(loss,decay['time']);err=abs(quadrature-decay['dissipation'][-1])/decay['energy'][0]
    add('dissipation-quadrature','Independent sampled integral agrees with ODE dissipation','relative integral error < 1e-4',err,err<1e-4,'Composite trapezoidal integral on 0.5 ms samples, compared with augmented ODE state.')
    pulse=simulate(skid,np.linspace(0,3,3001),base=floor_bump(),rtol=1e-10,max_step=.001)
    residual=np.max(abs(pulse['energy_residual']))/max(np.max(pulse['energy']),1e-12)
    add('moving-base-energy','Foundation input work balances dynamic spring energy and heat','relative energy residual < 1e-7',residual,residual<1e-7,'Relative-coordinate damper heat plus -sum(R)*base velocity input power.')
    refined,_,_=transient_metrics(skid,c,dt=.0005,rtol=1e-10)
    err=abs(refined['transient_relative_motion_mm']-nominal['transient_relative_motion_mm'])/refined['transient_relative_motion_mm']
    add('transient-refinement','Peak travel is stable when sampled four times more finely','relative peak change < 0.005',err,err<.005,'2 ms versus 0.5 ms output sampling; tighter integration tolerances. Acceptance includes 1% peak allowance.')
    limiting={}
    for metric in ['transient_relative_motion_mm','total_compression_mm','group_load_max_N','group_load_min_N']:
        corner=(min if metric=='group_load_min_N' else max)(corners,key=lambda row:row[metric])
        limiting.setdefault(corner['case'],corner)
    for case,corner in limiting.items():
        variant=candidate(c,selected['stiffness_N_m'],selected['damping_ratio'],corner['mass_factor'],corner['cg_shift_m'],(corner['front_stiffness_factor'],corner['rear_stiffness_factor']),corner['damping_factor'])
        fine,_,_=transient_metrics(variant,c,dt=.0005,rtol=1e-10)
        errors=[]
        for key in ['transient_relative_motion_mm','total_compression_mm','group_load_max_N','group_load_min_N']:
            errors.append(abs(fine[key]-corner[key])/max(abs(fine[key]),1e-12))
        add('corner-refinement-'+str(case),'Limiting corner '+str(case)+' peak/travel/load refinement','maximum relative metric change < 0.005',max(errors),max(errors)<.005,'2 ms versus0.5 ms samples with tighter integration, for corners limiting travel, compression, minimum or maximum load. A1% peak allowance is retained.')
    peak=local_peak_refinement(skid,c)
    add('band-refinement','Running-band peak agrees after local scalar maximization','refined peak minus grid peak < 0.01 N',peak['grid_to_refined_N'],peak['grid_to_refined_N']<.01,'601-point frequency grid with bounded refinement around interior maxima.')
    return checks,peak