"""Binomial, Clopper-Pearson and sign-test values, megabar baseline errors, and the ZrH3 ramp sensitivity quoted in the section on tests.

Run with the project environment: .venv/bin/python numerics/derived/test_statistics.py
The printed output is stored beside this file as test_statistics.out.txt.
"""
import pathlib
import numpy as np, json
from scipy import stats, optimize, integrate
kB=8.617333e-5
# 1. binomial for 12 of 15
for p in (2/3,0.7,0.8,0.85,0.9):
    print('P(>=12 of 15 | p=%.3f)=%.3f'%(p, stats.binom.sf(11,15,p)))
lo=stats.beta.ppf(0.025,12,4); hi=stats.beta.ppf(0.975,13,3)
print('CP 12/15: %.3f %.3f'%(lo,hi))
# zero false positives in n rows -> one-sided 95% upper bound
for n in (4,5,6,7): print('0 of',n,'upper95 = %.3f'%(1-0.05**(1/n)))
# sign test
for k in (3,4,5,6): print('sign test k wins of k: p=%.4f'%(0.5**k))
# 2. megabar baseline
d=json.load(open(pathlib.Path(__file__).resolve().parents[2] / 'data' / 'h_census.json'))
ref={}
for r in d['rows']:
    if r['id'] in ('quan-SH3-220','quan-LaH10-250','quan-CaH6-150','quan-MgH6-300','quan-YH10-300'):
        ref[r['compound']]=(r['pressure_GPa'],r['volume_A3'],r['rho_H'],r['eta_H_eV_A2'],r['h_eV_A'])
for base in (36.6,28.8,22.0,31.0):
    errs=[]
    for k,(P,V,rho,eta,h) in ref.items():
        errs.append(np.log(base/h))
        if base==36.6: print(k,P,V,round(rho,3),eta,round(h,1),'%+.0f%%'%(100*(base/h-1)))
    errs=np.array(errs)
    print('base',base,'MAE ln',round(np.mean(np.abs(errs)),3),'within30%:',int(np.sum(np.abs(base/np.array([v[4] for v in ref.values()])-1)<=0.30)))
hs=np.array([v[4] for v in ref.values()])
print('geo mean of five',np.exp(np.mean(np.log(hs))),'sd ln',np.std(np.log(hs),ddof=1), 'best-constant MAE ln', np.mean(np.abs(np.log(hs)-np.median(np.log(hs)))))
# 3. ramp kinetics: first-order, T_half on a linear ramp
def thalf(E,nu,beta):  # beta K/s
    f=lambda T: integrate.quad(lambda x: nu*np.exp(-E/(kB*x)), 50, T)[0]/beta-np.log(2)
    return optimize.brentq(f,60,1500)
beta=10/60
for nu in (1e13,):
    E0=optimize.brentq(lambda E: thalf(E,nu,beta)-235,0.3,1.5)
    print('E for T_half 235 K at 1e13:',E0)
    T0=thalf(E0,nu,beta)
    print('dE=+0.1:',thalf(E0+0.1,nu,beta)/T0-1,' dE=-0.1:',thalf(E0-0.1,nu,beta)/T0-1)
    for fac in (30,1/30,1000,1/1000):
        print('prefactor x',fac, thalf(E0,nu*fac,beta)/T0-1)
    # combined 0.1 eV and x30
    print('combined +0.1 & /30:',thalf(E0+0.1,nu/30,beta)/T0-1,' -0.1 & x30:',thalf(E0-0.1,nu*30,beta)/T0-1)
# 4. lifetimes single event
for (T,t,nu) in ((300,3.156e7,1e13),(300,3.156e7,1e10),(200,7*86400,1e13)):
    print(T,t,nu,'E=%.3f'%(kB*T*np.log(nu*t)))
# LRs
print('LR matbench',(0.927/0.073)/(0.153/0.847),'stanev',(0.74/0.26)/(0.034/0.966))
# needed h for 300K at 0.10: 8.7..12.0 /0.10
print(8.7/0.10,12.0/0.10, 87/44.6, 120/44.6, 87/64,120/64)
# floor
print((300/136.54)**2)