#!/usr/bin/env python3
"""Regenerate original course SVG charts. Optional --qa-dir writes inspection PNGs."""
from pathlib import Path
import argparse, csv, json, os, re
from html import escape
os.environ.setdefault('MPLCONFIGDIR','/tmp/statistics-storytelling-mpl')
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import FancyBboxPatch
ROOT=Path(__file__).resolve().parents[1]
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':11,'axes.titlesize':15,'axes.titleweight':'bold','axes.labelsize':11,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','svg.hashsalt':'aabha-statistics-v1','figure.facecolor':'white','axes.facecolor':'white','text.color':'#183548','axes.labelcolor':'#183548','xtick.color':'#183548','ytick.color':'#183548'})
BLUE='#075f86'; ORANGE='#9d4d00'; GRAY='#52616b'; TEAL='#066653'
def main():
    parser=argparse.ArgumentParser();parser.add_argument('--qa-dir',type=Path);args=parser.parse_args()
    (ROOT/'assets').mkdir(exist_ok=True)
    if args.qa_dir:args.qa_dir.mkdir(parents=True,exist_ok=True)
    r=json.loads((ROOT/'outputs/reference_results.json').read_text())
    with (ROOT/'data/shipments_clean.csv').open() as f:rows=list(csv.DictReader(f))
    durations=np.array([float(x['actual_minutes']) for x in rows]);distance=np.array([float(x['distance_km']) for x in rows]);priority=np.array([x['service']=='priority' for x in rows])
    values={}
    def save(fig,name,title,desc):
        fig.tight_layout(pad=1.4, rect=(0, .055, 1, 1))
        p=ROOT/'assets'/f'{name}.svg';fig.savefig(p,format='svg',metadata={'Date':None,'Title':title,'Description':desc})
        svg=p.read_text();svg=re.sub(r'<svg\b', '<svg role="img" aria-labelledby="chart-title chart-desc"',svg,count=1)
        # Inject accessible text immediately after the root opening tag, not the XML declaration.
        pos=svg.index('>',svg.index('<svg'))+1
        svg=svg[:pos]+f'<title id="chart-title">{escape(title)}</title><desc id="chart-desc">{escape(desc)}</desc>'+svg[pos:];p.write_text(svg)
        if args.qa_dir:fig.savefig(args.qa_dir/f'{name}.png',dpi=150)
        plt.close(fig)
    fig,ax=plt.subplots(figsize=(8.5,3.8));x=[8,9,10,11,12,30]
    ax.scatter(x,[1]*6,s=85,color=BLUE,zorder=3)
    for v in x:ax.text(v,1.08,str(v),ha='center')
    ax.axvline(10.5,color=TEAL,linestyle='--',label='Median 10.50');ax.axvline(80/6,color=ORANGE,linestyle=':',linewidth=2,label='Mean 13.33')
    ax.set(xlim=(5,33),ylim=(0.6,1.65),yticks=[],xlabel='Task duration (minutes)',title='One long task pulls the mean rightward');ax.legend(loc='upper right',frameon=False)
    fig.text(.5,.01,'Original six-task arithmetic illustration; not shipment or trial data.',ha='center',fontsize=9)
    save(fig,'centers','Mean and median of six task durations','Values 8, 9, 10, 11, 12, 30 minutes. Median 10.5, mean 13.333. The longest value pulls the mean above the median.')
    values['centers']={'durations':x,'mean':80/6,'median':10.5}
    fig,ax=plt.subplots(figsize=(9,4.2));ax.set(xlim=(0,10),ylim=(0,5));ax.axis('off')
    def box(x,y,w,h,txt,color=BLUE):
        ax.add_patch(FancyBboxPatch((x,y),w,h,boxstyle='round,pad=.12',facecolor='#f1f6f8',edgecolor=color,linewidth=1.5));ax.text(x+w/2,y+h/2,txt,ha='center',va='center',fontsize=12)
    box(.3,2,2.1,1.1,'Target\npopulation');box(3.9,2,2.1,1.1,'Enrolled\nstudy sample');box(7.8,3.2,1.8,.9,'Control');box(7.8,.8,1.8,.9,'Checklist')
    ax.annotate('',xy=(3.7,2.55),xytext=(2.6,2.55),arrowprops={'arrowstyle':'->','lw':2,'color':BLUE});ax.text(3.15,3.45,'Random\nsampling',ha='center',fontsize=11)
    for y in [3.65,1.25]:ax.annotate('',xy=(7.6,y),xytext=(6.15,2.55),arrowprops={'arrowstyle':'->','lw':2,'color':ORANGE})
    ax.text(6.9,4.5,'Random\nassignment',ha='center',fontsize=11,color=ORANGE);ax.text(2.1,.55,'Population coverage',ha='center',fontsize=11);ax.text(7.1,.15,'Intervention comparison within study',ha='center',fontsize=10)
    ax.set_title('Sampling and assignment answer different questions',pad=18)
    save(fig,'sampling_assignment','Sampling and assignment','Random sampling selects study units from a target population. Random assignment divides enrolled units into intervention groups. These support different kinds of conclusions under their assumptions.')
    fig,ax=plt.subplots(figsize=(8.5,4.4));ax.axis('off');ax.set_title('An 80% catch rate does not mean 80% of flags are correct',pad=20)
    table=ax.table(cellText=[['80','20','100'],['90','810','900'],['170','830','1,000']],rowLabels=['Delayed','Not delayed','Total'],colLabels=['Flagged','Not flagged','Total'],cellLoc='center',bbox=[.14,.26,.84,.62]);table.auto_set_font_size(False);table.set_fontsize(13)
    for (i,j),c in table.get_celld().items():c.set_edgecolor('#cbd8df');c.set_facecolor('#edf4f7' if i==0 or j==-1 else 'white')
    ax.text(.5,.08,'Among 170 flags, 80 identify delayed jobs: 47.06%',transform=ax.transAxes,ha='center',fontsize=13,color=BLUE)
    fig.text(.5,.01,'Original hypothetical alert example; not measured operational performance.',ha='center',fontsize=9)
    save(fig,'base_rates','Alert base-rate example','Hypothetical 1,000 jobs: 100 delayed, 900 not. Flags identify 80 delayed and 90 not delayed. 80 of 170 flagged jobs are delayed, or 47.06%.')
    values['base_rates']={'delayed_flagged':80,'delayed_unflagged':20,'not_delayed_flagged':90,'not_delayed_unflagged':810}
    fig,ax=plt.subplots(figsize=(8.5,5));intervals=[(22.5,25.0),(23.0,25.2),(22.6,24.6),(23.7,26.0),(21.5,24.5),(24.5,26.2),(23.2,25.1),(22.9,24.7),(23.4,26.1),(21.9,24.2)]
    for i,(lo,hi) in enumerate(intervals):
        hit=lo<=24<=hi;ax.plot([lo,hi],[i+1,i+1],color=BLUE if hit else ORANGE,lw=3,marker='|',markersize=10);ax.scatter([(lo+hi)/2],[i+1],color=BLUE if hit else ORANGE,s=20)
    ax.axvline(24,color=GRAY,linestyle='--',label='Fixed model mean = 24');ax.set(xlabel='Mean duration (minutes)',ylabel='Hypothetical study',yticks=range(1,11),title='The interval changes; the model mean stays fixed');ax.legend(loc='upper right',fontsize=10)
    fig.text(.5,.01,'Illustration only: 9 of these 10 constructed intervals cover 24. Not a 95% coverage estimate.',ha='center',fontsize=9)
    save(fig,'intervals','Repeated interval illustration','Ten constructed illustrative intervals, nine covering a fixed mean of 24 and one missing it. This is not a random simulation or an empirical 95% coverage claim.')
    values['intervals']={'illustrative_only':True,'fixed_mean':24,'intervals':intervals}
    fig,ax=plt.subplots(figsize=(8.5,4.5));counts,bins,_=ax.hist(durations,bins=np.arange(0,151,15),color=BLUE,edgecolor='white');ax.set(xlabel='Actual delivery duration (minutes)',ylabel='Shipment count',title='Delivery times vary widely in the eight-week fixture',xlim=(0,150));ax.text(.98,.95,'Synthetic shipments; n = 240',transform=ax.transAxes,ha='right',fontsize=10)
    save(fig,'shipments_histogram','Shipment duration distribution','Histogram of 240 synthetic shipment durations, equal 15-minute bins from 0 to 150. Actual times range 6.7 to 136.8 minutes. Source data/shipments_clean.csv.')
    values['shipments_histogram']={'bin_edges':bins.tolist(),'counts':counts.astype(int).tolist(),'n':240}
    trial=r['trial'];w=trial['welch_checklist_minus_control'] if 'welch_checklist_minus_control' in trial else None
    if w is None:
        w=next(v for v in trial.values() if isinstance(v,dict) and 'ci_low' in v)
    fig,ax=plt.subplots(figsize=(8.5,4.2));ax.errorbar(w['difference'],1,xerr=[[w['difference']-w['ci_low']],[w['ci_high']-w['difference']]],fmt='o',color=BLUE,capsize=8,linewidth=3,markersize=8)
    ax.axvline(0,color=GRAY,linewidth=1.5,label='No difference (0)');ax.axvline(-2,color=ORANGE,linestyle='--',label='Worthwhile threshold (−2)');ax.set(xlim=(-7,1),ylim=(.4,1.8),yticks=[],xlabel='Checklist minus control mean duration (minutes)',title='Estimated benefit spans the practical threshold')
    ax.text(w['difference'],1.2,'−3.20 minutes',ha='center',fontsize=12,color=BLUE);ax.text(w['ci_low'],.82,'−5.85',ha='center');ax.text(w['ci_high'],.82,'−0.55',ha='center');ax.text(-6.8,.5,'Negative = checklist faster',fontsize=10);ax.legend(loc='upper left',frameon=False,fontsize=10)
    fig.text(.5,.01,'95% Welch confidence interval for mean difference. Simulated trial; 20 per group.',ha='center',fontsize=10)
    save(fig,'trial_effect','Checklist mean effect and uncertainty','Checklist minus control is −3.2015 minutes, with a 95% Welch interval of −5.8496 to −0.5534. It excludes zero but crosses the −2-minute hypothetical worthwhile threshold. Simulated trial, 20 per group.')
    values['trial_effect']={'difference':w['difference'],'ci_low':w['ci_low'],'ci_high':w['ci_high'],'practical_threshold':-2,'n_per_group':20}
    slope=np.cov(distance,durations,ddof=1)[0,1]/np.var(distance,ddof=1);intercept=durations.mean()-slope*distance.mean();fitted=intercept+slope*distance
    for name,y,label,title in [('distance_time',durations,'Actual delivery duration (minutes)','Longer recorded routes tend to take longer'),('residuals',durations-fitted,'Residual: actual minus fitted (minutes)','Inspect the variation the fitted line misses')]:
        fig,ax=plt.subplots(figsize=(8.5,4.5));xx=distance if name=='distance_time' else fitted
        for mask,marker,color,service in [(priority,'o',BLUE,'Priority'),(~priority,'x',ORANGE,'Standard')]:ax.scatter(xx[mask],y[mask],marker=marker,color=color,s=25,alpha=.7,label=f'{service} (n={int(mask.sum())})')
        if name=='distance_time':ax.plot([distance.min(),distance.max()],[intercept+slope*distance.min(),intercept+slope*distance.max()],color=GRAY,linestyle='--',label='Descriptive fitted line')
        else:ax.axhline(0,color=GRAY,linestyle='--')
        ax.set(xlabel='Distance (km)' if name=='distance_time' else 'Fitted delivery duration (minutes)',ylabel=label,title=title);ax.legend(fontsize=9,frameon=False)
        fig.text(.5,.01,'Synthetic shipment fixture; n = 240; weeks 1–8. Descriptive association only.',ha='center',fontsize=9)
        save(fig,name,title,'Original chart of 240 synthetic shipment records. '+('Scatterplot distance versus actual minutes, descriptive line 12.6622 + 1.8058 × km. Pearson r is 0.9226. No causal inference.' if name=='distance_time' else 'Residuals actual minus fitted against fitted minutes, distinguished by service with circles and crosses. Shared weekly conditions and omitted service effects limit inference.'))
    fig,ax=plt.subplots(figsize=(8.5,4.2));rates=[20/51*100,21/189*100];ax.barh(['Late shipments','On-time shipments'],rates,color=[ORANGE,BLUE],height=.5);ax.set(xlim=(0,50),xlabel='Missing rating (%)',title='Rating coverage differs by observed late status');ax.invert_yaxis()
    for i,(v,txt) in enumerate(zip(rates,['20 / 51 = 39.22%','21 / 189 = 11.11%'])):ax.text(v+.7,i,txt,va='center',fontsize=11)
    fig.text(.5,.01,'Synthetic shipment fixture; n = 240; weeks 1–8. Missing opinions remain unknown.',ha='center',fontsize=9)
    save(fig,'missingness','Rating missingness by late status','20 of 51 late shipments lack ratings (39.22%), compared with 21 of 189 on-time shipments (11.11%). Observational synthetic fixture; no missing opinions are inferred.')
    values['missingness']={'late':{'missing':20,'n':51},'on_time':{'missing':21,'n':189}}
    (ROOT/'outputs/chart_values.json').write_text(json.dumps(values,indent=2)+'\n')
    print('Built 9 original accessible SVG charts and outputs/chart_values.json')
if __name__=='__main__':main()