"""Rebuild original scientific figures from published aggregate results (Matplotlib)."""
import json,math
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.ticker import PercentFormatter
D=Path('public/lab/data');O=Path('public/science/figures');O.mkdir(parents=True,exist_ok=True)
x=json.loads((D/'results.json').read_text());ex={e['slug']:e for e in x['experiments']}
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':12,'axes.spines.top':False,'axes.spines.right':False,'axes.labelcolor':'#152d3d','text.color':'#152d3d','axes.edgecolor':'#9aaebb','xtick.color':'#304957','ytick.color':'#304957','figure.facecolor':'white','axes.facecolor':'white','svg.fonttype':'none','axes.titleweight':'bold'})
green='#147d65';blue='#315c9b';orange='#bf592e'
def chart(slug,draw):
    fig,ax=plt.subplots(figsize=(9,5.5),layout='constrained');draw(ax);ax.grid(axis='y',alpha=.18);ax.set_axisbelow(True)
    fig.savefig(O/(slug+'.svg'));fig.savefig(O/(slug+'.png'),dpi=160);plt.close(fig)
def bars(ax,slug,ix=1,expected=2):
    rows=ex[slug]['rows'];xs=np.arange(len(rows));obs=np.array([r[ix] for r in rows],float)
    ax.bar(xs-.18,obs,.36,label='Observed simulation',color=green)
    if expected is not None:ax.bar(xs+.18,[r[expected] for r in rows],.36,label='Exact model expectation',color=blue)
    ax.set_xticks(xs,[str(r[0]) for r in rows]);ax.set_ylim(bottom=0);ax.legend(fontsize=10)

def freq(ax):
    rows=ex['ten-million-649-draws']['rows'];mu=rows[0][2];sd=math.sqrt(10000000*6/49*43/49)
    ax.bar([r[0] for r in rows],[(r[1]-mu)/sd for r in rows],color=green)
    ax.axhline(0,color=blue);ax.set(xlabel='Number label',ylabel='Deviation from expectation (single-number SDs)',title='Ten million fair 6/49 draws still have uneven counts')
chart('ten-million-649-draws',freq)
def hot(ax):
    s=ex['hot-cold-holdout']['summary'];ys=[s[k] for k in ['hot','cold','random']]
    ax.errorbar([0,1,2],ys,yerr=np.array(s['standard_errors'])*1.96,fmt='o',color=green,capsize=8,markersize=8,label='Mean ± 1.96 trial-level SE')
    ax.axhline(s['theory'],color=blue,ls='--',label='Exact expectation: 36/49');ax.set_xticks([0,1,2],['Hot','Cold','Random'])
    ax.set(ylabel='Matches per test draw',title='Past rankings, fresh test draws');ax.legend(fontsize=10)
chart('hot-cold-holdout',hot)
def ordered(ax):
    rows=ex['ordered-versus-irregular']['rows'];ax.plot(range(1,121),[r[1] for r in rows],'.',color=green,label='Observed count, each of 120 sets')
    ax.axhline(1e6/120,color=blue,label='Equal model expectation');ax.set(xlabel='Combination index (lexicographic order)',ylabel='Count in one million 3/10 draws',title='Equal probabilities do not force equal counts');ax.legend(fontsize=10)
chart('ordered-versus-irregular',ordered)
def gaps(ax):
    t=np.arange(0,111);ax.plot(t,(43/49)**t,color=blue,label='Exact: (43/49)^d');ax.set_yscale('log');ax.set(xlabel='Consecutive missed draws, d',ylabel='Probability of at least d misses (log scale)',title='A fixed number can wait a very long time');ax.legend()
chart('overdue-gaps',gaps)
def consecutive(ax):
    rows=ex['consecutive-pairs']['rows'];ax.bar([r[0] for r in rows],[r[1]/1e6 for r in rows],color=green);ax.yaxis.set_major_formatter(PercentFormatter(1));ax.set(xlabel='Adjacent pairs in a sorted draw',ylabel='Share of one million simulated draws',title='Consecutive labels are an ordinary category')
chart('consecutive-pairs',consecutive)
def overlap(ax):
    bars(ax,'overlap-between-draws');ax.set(xlabel='Shared labels between two complete draws',ylabel='Number of pairs',title='Independent draws can overlap')
chart('overlap-between-draws',overlap)
def parity(ax):
    bars(ax,'parity-categories');ax.set(xlabel='Odd labels in a six-number draw',ylabel='Number of draws',title='Balanced categories contain more combinations')
chart('parity-categories',parity)
def birthday(ax):
    r=ex['birthday-collisions']['rows'];t=np.arange(1,251);p=[-math.expm1(sum(math.log1p(-j/10000) for j in range(int(n)))) for n in t]
    ax.plot(t,p,color=blue,label='Exact collision probability');ax.scatter([v[0] for v in r],[v[1] for v in r],color=green,label='10,000 simulated sequences');ax.yaxis.set_major_formatter(PercentFormatter(1));ax.set(xlabel='Draws from 10,000 possible outcomes',ylabel='Probability of any repeated outcome',title='Any repeat is easier than one specified repeat');ax.legend(fontsize=10)
chart('birthday-collisions',birthday)
def multiple(ax):
    r=ex['multiple-testing']['rows'];ax.bar([v[0] for v in r],[v[1]/1e5 for v in r],color=green);ax.yaxis.set_major_formatter(PercentFormatter(1));ax.set(xlabel='False alarms among 20 independent tests',ylabel='Share of 100,000 simulated experiments',title='Searching more tests creates more opportunities');ax.set_xlim(-.6,9.6)
chart('multiple-testing',multiple)
def lln(ax):
    r=ex['law-large-numbers']['rows'];ax.loglog([v[0] for v in r],[v[2] for v in r],'o-',color=green,label='Observed mean absolute error');ax.loglog([v[0] for v in r],[v[4] for v in r],'s--',color=blue,label='Theoretical standard deviation');ax.set(xlabel='Tosses per independent sample (log scale)',ylabel='Error in proportion (log scale)',title='Relative error shrinks; the past is not repaid');ax.legend(fontsize=10)
chart('law-large-numbers',lln)
def modulo(ax):
    r=ex['modulo-bias']['rows'];t=np.arange(6);ax.bar(t-.18,[v[1]/1e6 for v in r],.36,color=orange,label='Naive modulo');ax.bar(t+.18,[v[2]/1e6 for v in r],.36,color=green,label='Rejection sampling');ax.axhline(1/6,color=blue,ls='--',label='Fair probability');ax.set_xticks(t,range(1,7));ax.set_ylim(0,.19);ax.yaxis.set_major_formatter(PercentFormatter(1));ax.set(xlabel='Die face',ylabel='Observed proportion',title='A uniform byte does not divide evenly into six');ax.legend(fontsize=10,loc='lower right')
chart('modulo-bias',modulo)
# Original exact mathematical figures supplement empirical figures.
def combin(ax):
    ns=np.arange(6,91);ax.semilogy(ns,[math.comb(int(n),6) for n in ns],color=blue);ax.scatter([49,59,90],[math.comb(n,6) for n in [49,59,90]],color=green);ax.set(xlabel='Pool size n (six labels selected)',ylabel='Combinations C(n,6), logarithmic scale',title='The outcome space grows combinatorially')
chart('combinations',combin)
def coverage(ax):
    t=np.linspace(0,4,250);ax.plot(t,1-np.exp(-t),color=blue,label='Independent trials: 1 − exp(−Np), approximation');ax.axhline(1,color='#777',ls=':');ax.yaxis.set_major_formatter(PercentFormatter(1));ax.set(xlabel='Expected success count Np',ylabel='Probability of at least one success',title='One expected success is not a guarantee');ax.legend(fontsize=9)
chart('coverage',coverage)
def partition(ax):
    from matplotlib.patches import Rectangle
    for j in range(64):
        heads=bin(j).count('1');ax.add_patch(Rectangle((j%8,7-j//8),.85,.85,facecolor=green if heads==3 else '#d9e3eb',edgecolor='white'))
    ax.set(xlim=(-.1,8),ylim=(-.1,8),aspect='equal',title='64 equally likely six-toss strings');ax.set_axis_off();ax.text(8.2,5.5,'20 strings contain\nthree heads and three tails.\n\nEach square still has\nprobability 1/64.',fontsize=13);ax.set_xlim(-.1,14)
chart('category-sketch',partition)
print('Created 14 original SVG + PNG figure pairs')
def sums_plot(ax):
    d=np.zeros((7,280),dtype=np.int64);d[0,0]=1
    for v in range(1,50):
        for k in range(6,0,-1):d[k,v:]+=d[k-1,:280-v]
    ax.plot(range(21,280),d[6,21:],color=blue);ax.axvline(150,color=green,ls='--',label='Expected sum: 150');ax.set(xlabel='Sum of six selected labels',ylabel='Exact number of combinations',title='Middle sums have more representations');ax.set_ylim(bottom=0);ax.legend()
chart('lottery-sums',sums_plot)