import random,math,statistics,json,csv,sys,hashlib
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
rng=random.Random(21001);half=185.52;repetitions=200;rows=[];samples=[];groups=[]
for n in (100,1000,10000):
    estimates=[statistics.median(rng.expovariate(math.log(2)/half) for _ in range(n)) for _ in range(repetitions)]
    groups.append(estimates)
    samples.extend(dict(n=n,repetition=i+1,median_estimate_s=x,relative_error=(x-half)/half) for i,x in enumerate(estimates))
    rows.append(dict(n=n,repetitions=repetitions,true_half_s=half,mean_estimate_s=statistics.mean(estimates),relative_std=statistics.stdev(estimates)/half,asymptotic_relative_std=1/(math.log(2)*math.sqrt(n))))
curve=[dict(n=n,asymptotic_relative_std=1/(math.log(2)*math.sqrt(n))) for n in sorted(set(round(10**(2+k/50)) for k in range(101)))]
survival=[dict(time_s=k*half/100,time_in_assigned_medians=k/100,expected_surviving_fraction=math.exp(-math.log(2)*k/100)) for k in range(501)]
def save(name,rr):
    with(P/name).open('w',newline='',encoding='utf8')as f:
        w=csv.DictWriter(f,fieldnames=list(rr[0]));w.writeheader();w.writerows(rr)
save('median-estimates.csv',samples);save('summaries.csv',rows);save('asymptotic-curve.csv',curve);save('survival.csv',survival)
# Independent moments from sums, rather than statistics.mean/stdev.
max_moment=0
for r,x in zip(rows,groups):
    mean=math.fsum(x)/len(x);sd=math.sqrt(math.fsum((v-mean)**2 for v in x)/(len(x)-1))/half
    max_moment=max(max_moment,abs(sd-r['relative_std']));assert abs(sd-r['relative_std'])<1e-14
assert survival[0]['expected_surviving_fraction']==1
assert survival[100]['expected_surviving_fraction']==.5
assert abs(survival[500]['expected_surviving_fraction']-1/32)<1e-16
assert all(a['expected_surviving_fraction']>=b['expected_surviving_fraction'] for a,b in zip(survival,survival[1:]))
assert math.isclose(1/(math.log(2)*math.sqrt(400)),.5/(math.log(2)*math.sqrt(100)))
legacy=P.parent.parent/'prototypes/met001_results.json'
legacy_check=None
if legacy.exists():
    old=json.loads(legacy.read_text(encoding='utf8'));legacy_check=max(abs(r[k]-o[k]) for r,o in zip(rows,old['rows']) for k in r);assert legacy_check==0
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(1,2,figsize=(10,4.8),layout='constrained')
ax[0].boxplot([[100*(v-half)/half for v in x]for x in groups],tick_labels=['100','1,000','10,000'],whis=1.5,showfliers=True)
ax[0].axhline(0,color='gray',ls='--');ax[0].set(xlabel='Independent lifetimes per repetition',ylabel='Sample median error (%)',title='A  200 reproduced estimates at each size')
ax[1].loglog([r['n']for r in curve],[r['asymptotic_relative_std']for r in curve],label='Asymptotic approximation',color='#7957a8')
ax[1].loglog([r['n']for r in rows],[r['relative_std']for r in rows],'o',label='Reproduced sample standard deviation',color='#278c89')
ax[1].set(xlabel='Independent lifetimes per repetition',ylabel='Standard deviation / assigned median',title='B  Spread decreases with sample size');ax[1].legend(fontsize=8);ax[1].grid(alpha=.2)
fig.suptitle('Synthetic exponential lifetimes: complete observations, assigned median 185.52 s')
for ext in ['png','svg','pdf']:fig.savefig(P/f'median-spread.{ext}')
plt.close(fig)
fig,ax=plt.subplots(figsize=(7,4.5),layout='constrained');ax.plot([r['time_s']for r in survival],[r['expected_surviving_fraction']for r in survival],color='#3769a2')
ax.set(xlabel='Time (s)',ylabel='Expected fraction not yet ended',title='Analytic exponential survival\nAssigned teaching median 185.52 s');ax.axhline(.5,color='gray',ls='--');ax.axvline(half,color='gray',ls='--');ax.grid(alpha=.2)
for ext in ['png','svg','pdf']:fig.savefig(P/f'expected-survival.{ext}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,legacy_summary_max_difference=legacy_check,independent_std_max_difference=max_moment,survival_limits=True,asymptotic_scaling=True,rows=dict(estimates=600,summaries=3,asymptotic=len(curve),survival=501)),indent=2)+'\n')
(P/'study.json').write_text(json.dumps(dict(runtime=dict(python=sys.version,numpy=np.__version__,matplotlib=matplotlib.__version__),baseline_sha256={name:hashlib.sha256((P.parent.parent/"prototypes"/name).read_bytes()).hexdigest() for name in ["met001_statistics.py","met001_results.json"] if (P.parent.parent/"prototypes"/name).exists()},synthetic=True,seed=21001,python_rng='random.Random; expovariate; one continuous stream',order='N100 then1000 then10000; 200 repetitions at each N',assigned_median_s=half,estimator='statistics.median; even N averages two central values',spread='sample standard deviation ddof1 divided by assigned median',observation='complete independent exponential lifetimes',boxplot='median, quartiles, whiskers to1.5IQR, outliers shown; not confidence intervals',not_demonstrated=['isotope identity','clock precision','detector behavior','SI unit realization']),indent=2)+'\n')
print(rows);print('legacy difference',legacy_check)