"""Gaussian one-sided threshold teaching model; synthetic analytical probabilities."""
from pathlib import Path
import math,json,csv
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
OUT=Path(__file__).resolve().parent/'signal-study';OUT.mkdir(exist_ok=True)
def tail(x):return .5*math.erfc(x/math.sqrt(2))
def family(p,n):return -math.expm1(n*math.log1p(-p)) if p<1 else 1.0
rows=[]
for n in [1,100,10000]:
 for i in range(81):
  threshold=i/10;p=tail(threshold);rows.append({'threshold_sigma':threshold,'independent_tests':n,'per_test_false_alarm':p,'at_least_one_false_alarm':family(p,n),'detection_at_signal_3sigma':tail(threshold-3),'detection_at_signal_5sigma':tail(threshold-5)})
with (OUT/'study.csv').open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
(OUT/'study.json').write_text(json.dumps({'model':'One-sided Gaussian score: noise N(0,1), signal N(amplitude,1).','assumptions':['Known noise mean and variance','Independent tests for family probability','A known scalar score, not unknown-template search','No correlated interference, drift or variance estimation'],'equations':{'false_alarm':'Q(threshold)','detection':'Q(threshold-amplitude)','family':'1-(1-false_alarm)^N'},'cases':rows},indent=2))
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','pdf.fonttype':42})
fig,ax=plt.subplots(1,2,figsize=(10,4.3),layout='constrained')
x=[i/100 for i in range(801)]
ax[0].plot(x,[tail(t) for t in x],label='Noise only: false alarm',ls='--',color='#a31b78');ax[0].plot(x,[tail(t-3) for t in x],label='Signal mean = 3σ',color='#146b8c');ax[0].plot(x,[tail(t-5) for t in x],label='Signal mean = 5σ',color='#8a49a3');ax[0].set(xlabel='Threshold (noise standard deviations)',ylabel='Probability per test',title='(a) Detection versus threshold');ax[0].legend(frameon=False,fontsize=8)
for n,col in zip([1,100,10000],['#146b8c','#8a49a3','#b76a20']):ax[1].semilogy(x,[family(tail(t),n) for t in x],label=f'{n:,} independent tests',color=col)
ax[1].set(xlabel='Threshold (noise standard deviations)',ylabel='Probability of at least one false alarm',title='(b) Searching more places');ax[1].legend(frameon=False,fontsize=8)
for a in ax:a.grid(alpha=.2)
for ext in ['svg','png','pdf']:fig.savefig(OUT/f'threshold-study.{ext}',dpi=200,bbox_inches='tight')
assert tail(0)==.5 and abs(family(.01,1)-.01)<1e-14
assert all(rows[i]['per_test_false_alarm']>=rows[i+1]['per_test_false_alarm'] for i in range(80))
assert abs(family(tail(5),10000)-(1-(1-tail(5))**10000))<1e-11
print(json.dumps({'rows':len(rows),'per_test_5sigma':tail(5),'family_5sigma_10000':family(tail(5),10000),'detection_3sigma_at_threshold5':tail(2)}))