#!/usr/bin/env python3
"""COMPLEXA C301-C310: source-derived, deterministic, standard-library figures.

Usage:
    python make_figures.py --out content/supplementary
    python make_figures.py --check content/supplementary

Data are exact evaluations of explicit didactic models, not empirical samples.
SVG is generated directly so graphical elements remain vector and searchable.
"""
from __future__ import annotations
import argparse
import base64
import csv
from hashlib import sha256
from html import escape
from io import BytesIO, StringIO
from math import exp, log, sqrt, isclose, sin, cos, pi, comb
from xml.etree import ElementTree
from pathlib import Path
import zipfile

ROOT = Path(__file__).resolve().parents[2]
NAMES = {
 301:("Correlación perfecta no garantiza acuerdo", "Las mismas cinco lecturas, sesgo constante de tres unidades"),
 302:("Una verosimilitud no es una posterior", "Ocho éxitos en diez frente a ochenta éxitos en cien"),
 303:("La función de pérdida cambia la calibración", "Dos ajustes de una misma decaída exponencial"),
 304:("Identificación práctica según la región medida", "Dos curvas de saturación casi idénticas a baja entrada"),
 305:("Una previa cambia tras observar veinte intentos", "Distribuciones beta previa y posterior: catorce éxitos"),
 306:("Ridge contrae; lasso puede anular", "Tres coeficientes con columnas ortonormales"),
 307:("Una media adecuada puede ocultar demasiados ceros", "Comprobación predictiva posterior de diez recuentos"),
 308:("AIC y BIC pueden elegir modelos distintos", "Recta y curva cúbica evaluadas con cien observaciones"),
 309:("La información depende de la coordenada", "Información esperada Bernoulli en p y en logit(p)"),
 310:("Un ajuste excelente y parámetros sin identificar", "Dos resortes: solo se observa la suma de compliancias"),
}

NAVY="#183448"; TEAL="#087F87"; OCHRE="#C56D35"; PURPLE="#7465AD"
INK="#263B4A"; MUTED="#66808B"; GRID="#E2EAEF"; BG="#F5F8FA"; GREEN="#4A9870"


def v(x):
    if isinstance(x,str): return x
    if isinstance(x,int): return str(x)
    return f"{x:.10f}"


def csv_bytes(header, rows):
    f=StringIO(newline="")
    w=csv.writer(f,lineterminator="\n")
    w.writerow(header)
    for row in rows: w.writerow([v(x) for x in row])
    return f.getvalue().encode("utf-8")


def X(text): return escape(str(text),quote=True)


class SVG:
    def __init__(self,n,title,desc):
        self.n=n; self.a=[]
        self.a.append(f'<svg xmlns="http://www.w3.org/2000/svg" width="1200" height="760" viewBox="0 0 1200 760" role="img" aria-labelledby="ttl{n} des{n}">')
        self.a.append(f'<title id="ttl{n}">Figura {n}.1. {X(title)}</title>')
        self.a.append(f'<desc id="des{n}">{X(desc)}. Valores de modelos definidos en el manuscrito, no resultados observacionales.</desc>')
        self.rect(0,0,1200,760,"#FFFFFF")
        self.rect(0,0,13,760,TEAL)
        self.rect(43,28,1114,99,BG,15)
        self.text(65,56,"C O M P L E X A   /   LABORATORIO REPRODUCIBLE",13,TEAL,700)
        self.text(65,93,title,29,NAVY,700)
        self.text(1133,62,f"C{n}",24,NAVY,700,anchor="end")
        self.text(68,150,desc,17,MUTED,400)
        self.line(55,691,1153,691,GRID,2)
        self.text(64,718,"MODELO MATEMÁTICO / DATOS DIDÁCTICOS",12,MUTED,700)
        self.text(1138,718,f"complexa.science  ·  {n}.1",13,MUTED,500,anchor="end")
    def rect(self,x,y,w,h,color,r=0,stroke=None,sw=1):
        self.a.append(f'<rect x="{x}" y="{y}" width="{w}" height="{h}" rx="{r}" fill="{color}"'+(f' stroke="{stroke}" stroke-width="{sw}"' if stroke else "")+'/>')
    def text(self,x,y,s,size=16,color=INK,weight=400,anchor="start",rotate=None):
        rot=f' transform="rotate({rotate} {x} {y})"' if rotate is not None else ""
        self.a.append(f'<text x="{x}" y="{y}" text-anchor="{anchor}" fill="{color}" font-family="DejaVu Sans, Arial, sans-serif" font-size="{size}" font-weight="{weight}"{rot}>{X(s)}</text>')
    def line(self,x1,y1,x2,y2,color,sw=1,dash=None):
        self.a.append(f'<line x1="{x1:.2f}" y1="{y1:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" stroke="{color}" stroke-width="{sw}" stroke-linecap="round"'+(f' stroke-dasharray="{dash}"' if dash else "")+'/>')
    def circle(self,x,y,r,color,stroke=None):
        self.a.append(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="{r}" fill="{color}"'+(f' stroke="{stroke}"' if stroke else "")+'/>')
    def path(self,p,color,sw=3,fill="none",dash=None):
        self.a.append(f'<path d="{p}" stroke="{color}" stroke-width="{sw}" stroke-linejoin="round" stroke-linecap="round" fill="{fill}"'+(f' stroke-dasharray="{dash}"' if dash else "")+'/>')
    def note(self,s):
        self.rect(65,641,1070,33,BG,8)
        self.text(80,663,s,14,INK,500)
    def svg(self):return ("\n".join(self.a)+"\n</svg>\n").encode("utf-8")


def plot(svg,area,xmin,xmax,ymin,ymax,xticks,yticks,xlabel,ylabel,series,tag=None,highlight=None,legend=True):
    """area=(x,y,w,h); series list (label,color,[(x,y)...])."""
    x,y,w,h=area
    svg.rect(x,y,w,h,"#FFFFFF",12,GRID)
    svg.text(x+20,y+33,tag or "MODELO Y COMPARACIÓN",17,NAVY,700)
    left=x+72; right=x+w-25; top=y+(91 if legend else 65); bottom=y+h-64
    def xx(t):return left+(t-xmin)*(right-left)/(xmax-xmin)
    def yy(t):return bottom-(t-ymin)*(bottom-top)/(ymax-ymin)
    for val in yticks:
        cy=yy(val);svg.line(left,cy,right,cy,GRID,1)
        label=f"{val:.2f}" if abs(val)<2 and val!=int(val) else f"{val:g}"
        svg.text(left-12,cy+5,label,12,MUTED,400,anchor="end")
    for val in xticks:
        cx=xx(val);svg.line(cx,bottom,cx,bottom+5,MUTED,1.1)
        label=f"{val:.2f}" if abs(val)<2 and val!=int(val) else f"{val:g}"
        svg.text(cx,bottom+23,label,12,MUTED,400,anchor="middle")
    svg.line(left,bottom,right,bottom,NAVY,1.7);svg.line(left,top,left,bottom,NAVY,1.7)
    svg.text((left+right)/2,y+h-17,xlabel,14,MUTED,500,anchor="middle")
    svg.text(x+23,(top+bottom)/2,ylabel,13,MUTED,500,anchor="middle",rotate=-90)
    for lab,color,vals in series:
        points=[(xx(a),yy(b)) for a,b in vals]
        d=" ".join(("M" if j==0 else "L")+f"{a:.2f},{b:.2f}" for j,(a,b) in enumerate(points))
        svg.path(d,color,3.8)
    if highlight:
        for item in highlight:
            hx,hy,itemtext,color=item
            cx=xx(hx);cy=yy(hy)
            svg.circle(cx,cy,5.3,"#FFFFFF",stroke=color)
            svg.circle(cx,cy,2.4,color)
            svg.text(cx+8,cy-10,itemtext,12,color,700)
    if legend:
        # Draw compact legend below panel heading, above plot.
        lx=x+w-22
        for label,color,vals in reversed(series):
            length=min(280,len(label)*8.3+28);lx-=length
            svg.line(lx,y+60,lx+18,y+60,color,3)
            svg.text(lx+24,y+65,label,13,INK,600)
            lx-=14
    return xx,yy


def bar_panel(svg,area,title,bars,ymax,unit="puntos",baseline=0):
    x,y,w,h=area
    svg.rect(x,y,w,h,"#FFFFFF",12,GRID)
    svg.text(x+20,y+33,title,17,NAVY,700)
    x0=x+70; x1=x+w-25; y0=y+66; y1=y+h-64
    def scale(z):return y1-(z-baseline)*(y1-y0)/(ymax-baseline)
    for q in range(5):
        z=baseline+(ymax-baseline)*q/4
        cy=scale(z);svg.line(x0,cy,x1,cy,GRID,1)
        svg.text(x0-9,cy+4,f"{z:.2f}" if ymax<3 else f"{z:.1f}",12,MUTED,400,"end")
    count=len(bars);block=(x1-x0)/count
    for i,(label,val,color) in enumerate(bars):
        bw=min(67,block*.61);bx=x0+i*block+block/2-bw/2
        by=scale(val);svg.rect(bx,by,bw,y1-by,color,4)
        svg.text(bx+bw/2,by-10,f"{val:.3f}" if ymax<=3 else f"{val:.2f}",13,color,700,anchor="middle")
        svg.text(bx+bw/2,y1+26,label,12,INK,600,anchor="middle")
    svg.text((x0+x1)/2,y+h-17,unit,13,MUTED,500,anchor="middle")





# Each record is derived from equations and explicit didactic data in the chapter.
# No row represents a new real-world or clinical observation.
def beta_density(a,b,p):
    if p <= 0 or p >= 1:return 0.0
    return exp((a-1)*log(p)+(b-1)*log(1-p)+__import__('math').lgamma(a+b)-__import__('math').lgamma(a)-__import__('math').lgamma(b))


def gamma_poisson_zeros(n,alpha,rate):
    """Exact analytic predictive P(K=k), gamma-rate prior updated by 10 counts."""
    from decimal import Decimal, localcontext
    with localcontext() as ctx:
        ctx.prec=75
        b=Decimal(rate)
        out=[]
        for k in range(n+1):
            value=sum((Decimal(-1)**j)*Decimal(comb(n-k,j))*(b/(b+Decimal(k+j)))**alpha for j in range(n-k+1))
            out.append(float(Decimal(comb(n,k))*value))
    assert abs(sum(out)-1)<1e-12 and min(out)>-1e-12
    return out


def make_data():
    D={}
    # C301: five records compared using both code versions, exactly as in the manuscript.
    D[301]=(['record','method_A_omissions','method_B_omissions','B_minus_A','average_methods'],
      [(i,a,a+3,3,a+1.5) for i,a in enumerate((2,4,6,8,10),1)])
    # C302: normalized LIKELIHOOD, not probabilities over p. Count data differ in n.
    rows=[]
    for j in range(101):
        p=j/100
        l10=comb(10,8)*p**8*(1-p)**2
        rel10=(p/.8)**8*((1-p)/.2)**2 if 0<p<1 else 0
        rel100=(p/.8)**80*((1-p)/.2)**20 if 0<p<1 else 0
        rows.append((p,l10,rel10,rel100,8,10,80,100))
    D[302]=(['candidate_p','likelihood_eight_of_ten','relative_likelihood_n10','relative_likelihood_n100','successes_n10','trials_n10','successes_n100','trials_n100'],rows)
    # C303: the two source-reported fits optimize DIFFERENT loss functions.
    obs=[11.0,7.1,4.6,3.4,2.3]
    rows=[]
    A,k=10.8907,.4072;AL,kL=10.6058,.3866
    for t,y in enumerate(obs):
        a=A*exp(-k*t);b=AL*exp(-kL*t)
        rows.append((t,y,a,b,y-a,y-b,log(y),log(b)))
    D[303]=(['minute','observed_didactic_signal','fit_original_scale','fit_log_scale','residual_original','residual_log_fit_original_units','log_observed','log_fitted_log_scale'],rows)
    # C304: practical identifiability depends on concentration range; curves not observed samples.
    ss=[1,2,5,10]+[10+50*i for i in range(1,101)]
    D[304]=(['input_S','response_V100_K1000','response_V200_K2000','absolute_difference','relative_difference'],[(s,100*s/(1000+s),200*s/(2000+s),200*s/(2000+s)-100*s/(1000+s),(200*s/(2000+s))/(100*s/(1000+s))-1) for s in ss])
    # C305: prior Beta(2,2), posterior Beta(16,8) after 14 successes and six failures.
    D[305]=(['probability_theta','prior_beta_2_2_density','posterior_beta_16_8_density','prior_mean','posterior_mean'],[(i/100,beta_density(2,2,i/100),beta_density(16,8,i/100),.5,2/3) for i in range(101)])
    # C306: exact orthonormal-design shrinkage with lambda=0.5.
    ordinary=[.2,-1.3,3]; lam=.5
    D[306]=(['coefficient_index','ordinary_z','ridge_lambda_0_5','lasso_lambda_0_5','lambda'],[(i,z,z/(1+lam),(1 if z>=0 else -1)*max(abs(z)-lam,0),lam) for i,z in enumerate(ordinary,1)])
    # C307: exact posterior predictive of zeros under Poisson-gamma posterior Gamma(21,11).
    pz=gamma_poisson_zeros(10,21,11)
    fixed=exp(-21/11)
    D[307]=(['number_zeros','posterior_predictive_probability','plug_in_Poisson_probability','observed_number_zeros','posterior_tail_at_least_5'],[(j,pz[j],comb(10,j)*fixed**j*(1-fixed)**(10-j),5,sum(pz[5:])) for j in range(11)])
    # C308: two model selection scores; both compared on exactly the same n=100 data.
    D[308]=(['number_parameters','log_likelihood','AIC','BIC','n_observations','is_cubic'],[(k,ell,-2*ell+2*k,-2*ell+k*log(100),100,int(k==4)) for k,ell in [(2,-140),(4,-136)]])
    # C309: local Fisher metric changes with parameter scale; no information disappears.
    D[309]=(['probability_p','information_p_for_n10','information_logit_for_n10','inverse_information_p','inverse_information_logit','sample_size'],[(i/100,10/((i/100)*(1-i/100)),10*(i/100)*(1-i/100),(i/100)*(1-i/100)/10,1/(10*(i/100)*(1-i/100)),10) for i in range(5,96)])
    # C310: exact continuum of compliances. Constant profile vs false fixed-other slice.
    # Use c1 0.04...2.94; all profiled c2 values strictly positive.
    D[310]=(['compliance_1','profiled_compliance_2','profile_sum','profile_SSE','fixed_c2_1_99_SSE','observed_total_compliance','fixed_compliance_2'],[(i*0.03,2.99-i*.03,2.99,.097,.097+30*(i*.03+1.99-2.99)**2,2.99,1.99) for i in range(1,99)])
    return D


def point_legend(s,points,xstart=120,ystart=245,step=260):
    for i,(label,col) in enumerate(points):
        x=xstart+(i%3)*step;y=ystart+(i//3)*28
        s.line(x,y-6,x+23,y-6,col,3.8)
        s.text(x+30,y,label,13,INK,600)


def charts(n,data):
    title,desc=NAMES[n];s=SVG(n,title,desc)
    if n==301:
        plot(s,(55,174,537,453),1,14,1,15,[2,4,6,8,10,12,14],[2,4,6,8,10,12,14],
             'Omisiones con método A','Omisiones con método B',[
             ('Coincidencia perfecta',MUTED,[(2,2),(13,13)]),
             ('Método B = A + 3',TEAL,[(r[1],r[2]) for r in data])],
             'CORRELACIÓN DE PEARSON = 1')
        plot(s,(608,174,537,453),3,12,0,4.2,[3,5,7,9,11],[0,1,2,3,4],
             'Promedio de ambos métodos','Diferencia B − A',[
             ('Sesgo observado: +3',OCHRE,[(r[4],r[3]) for r in data])],
             'ACUERDO: SESGO CONSTANTE',legend=False)
        s.note('La correlación mide asociación; la diferencia entre versiones revela un error de acuerdo.')
    elif n==302:
        plot(s,(55,174,1090,453),0,1,0,1.05,[0,.2,.4,.6,.8,1],[0,.25,.5,.75,1],
             'Parámetro hipotético p','Verosimilitud relativa',[
             ('8 aciertos / 10',TEAL,[(r[0],r[2]) for r in data]),
             ('80 aciertos / 100',OCHRE,[(r[0],r[3]) for r in data])],
             'MISMA PROPORCIÓN · DISTINTA INFORMACIÓN',[(.8,1,'máximo p = 0,80',TEAL)])
        s.note('Normalizar por el máximo NO convierte la verosimilitud en probabilidad posterior.')
    elif n==303:
        plot(s,(55,174,537,453),0,4,0,12,[0,1,2,3,4],[0,2,4,6,8,10,12],
             'Tiempo (min)','Señal',[
             ('Observado didáctico',NAVY,[(r[0],r[1]) for r in data]),
             ('Mínimos cuadrados',TEAL,[(r[0],r[2]) for r in data]),
             ('Ajuste en log',OCHRE,[(r[0],r[3]) for r in data])],
             'MISMAS LECTURAS · DISTINTO CRITERIO',legend=False)
        plot(s,(608,174,537,453),0,4,-.55,.55,[0,1,2,3,4],[-.5,-.25,0,.25,.5],
             'Tiempo (min)','Residuo en escala original',[
             ('Cuadrados originales',TEAL,[(r[0],r[4]) for r in data]),
             ('Cuadrados logarítmicos',OCHRE,[(r[0],r[5]) for r in data])],
             'PONDERACIÓN DISTINTA DEL ERROR',legend=False)
        s.note('Original: A=10,8907; k=0,4072 · log: A=10,6058; k=0,3866. Son objetivos distintos.')
    elif n==304:
        x=[r for r in data if r[0]<=10]
        plot(s,(55,174,537,453),0,10,0,1.04,[0,2,4,6,8,10],[0,.2,.4,.6,.8,1],
             'Entrada S (rango observado)','Respuesta v',[
             ('V=100, K=1000',TEAL,[(r[0],r[1]) for r in x]),
             ('V=200, K=2000',OCHRE,[(r[0],r[2]) for r in x])],
             'SEPARACIÓN PRÁCTICA MUY PEQUEÑA',legend=False)
        plot(s,(608,174,537,453),0,5010,0,205,[0,1000,2000,3000,4000,5000],[0,50,100,150,200],
             'Entrada S (ampliación propuesta)','Respuesta v',[
             ('V=100, K=1000',TEAL,[(0,0)]+[(r[0],r[1]) for r in data]),
             ('V=200, K=2000',OCHRE,[(0,0)]+[(r[0],r[2]) for r in data])],
             'OTRA REGIÓN DEL DISEÑO',legend=False)
        s.note('Dos curvas diferentes pueden resultar casi indistinguibles si solo se mide S entre 1 y 10.')
    elif n==305:
        plot(s,(55,174,1090,453),0,1,0,5.1,[0,.2,.4,.6,.8,1],[0,1,2,3,4,5],
             'Probabilidad desconocida θ','Densidad sobre θ',[
             ('Previa Beta(2,2)',OCHRE,[(r[0],r[1]) for r in data]),
             ('Posterior Beta(16,8)',TEAL,[(r[0],r[2]) for r in data])],
             '14 ÉXITOS Y 6 FRACASOS · PREVIA EXPLÍCITA',[(2/3,beta_density(16,8,2/3),'media posterior = 2/3',TEAL)])
        s.note('Una distribución posterior se obtiene combinando previa y verosimilitud; no es la curva L.')
    elif n==306:
        plot(s,(55,174,1090,453),1,3,-1.6,3.3,[1,2,3],[-1,0,1,2,3],
             'Coeficiente j','Estimación de coeficiente',[
             ('Mínimos cuadrados',NAVY,[(r[0],r[1]) for r in data]),
             ('Ridge (λ=0,5)',TEAL,[(r[0],r[2]) for r in data]),
             ('Lasso (λ=0,5)',OCHRE,[(r[0],r[3]) for r in data])],
             'DISEÑO ORTONORMAL · CONTRACCIÓN ANALÍTICA',legend=False)
        point_legend(s,[('Sin penalización',NAVY),('Ridge',TEAL),('Lasso',OCHRE)],95,235,220)
        s.note('Ridge contrae todos los valores; lasso lleva el coeficiente 0,2 exactamente a cero.')
    elif n==307:
        plot(s,(55,174,1090,453),0,10,0,.44,[0,2,4,5,6,8,10],[0,.1,.2,.3,.4],
             'Número de recuentos iguales a cero','Probabilidad predictiva',[
             ('Gamma–Poisson posterior',TEAL,[(r[0],r[1]) for r in data]),
             ('Poisson tasa fija 21/11',OCHRE,[(r[0],r[2]) for r in data])],
             '10 RÉPLICAS DE RECUENTO · 5 CEROS OBSERVADOS',legend=False)
        point_legend(s,[('Predictiva posterior',TEAL),('Tasa fija',OCHRE)],92,235,245)
        s.note(f'La media total se ajusta; P(5 o más ceros | datos, Poisson–Gamma) = {data[0][4]:.4f}.')
    elif n==308:
        plot(s,(55,174,537,453),1.65,4.35,275,295,[2,3,4],[275,280,285,290,295],
             'Parámetros ajustados','AIC (menor es preferible)',[
             ('AIC',TEAL,[(r[0],r[2]) for r in data])],
             'RECTA k=2 / CÚBICA k=4',legend=False)
        plot(s,(608,174,537,453),1.65,4.35,285,296,[2,3,4],[285,288,291,294,296],
             'Parámetros ajustados','BIC (menor es preferible)',[
             ('BIC',OCHRE,[(r[0],r[3]) for r in data])],
             'PENALIZA MÁS CON n=100',legend=False)
        s.note('AIC elige la curva cúbica (280); BIC elige la recta (289,210 frente a 290,421).')
    elif n==309:
        plot(s,(55,174,537,453),.05,.95,0,220,[.05,.2,.4,.6,.8,.95],[0,50,100,150,200],
             'Probabilidad p','Información I(p)',[
             ('Escala p',TEAL,[(r[0],r[1]) for r in data])],
             'UNIDADES: PROBABILIDAD⁻²',legend=False)
        plot(s,(608,174,537,453),.05,.95,0,2.7,[.05,.2,.4,.6,.8,.95],[0,.5,1,1.5,2,2.5],
             'Probabilidad p','Información I(logit p)',[
             ('Escala logit(p)',OCHRE,[(r[0],r[2]) for r in data])],
             'UNIDADES: LOG-ODDS⁻²',legend=False)
        s.note('I_p = 10/[p(1−p)]; I_logit = 10p(1−p). La escala del parámetro cambia.')
    elif n==310:
        plot(s,(55,174,537,453),0,3,0,3,[0,.5,1,1.5,2,2.5,3],[0,.5,1,1.5,2,2.5,3],
             'Compliancia del resorte 1','Compliancia del resorte 2',[
             ('Suma ajustada 2,99',TEAL,[(r[0],r[1]) for r in data])],
             'UN CONTINUO DE SOLUCIONES',legend=False)
        plot(s,(608,174,537,453),0,3,0,125,[0,.5,1,1.5,2,2.5,3],[0,30,60,90,120],
             'Compliancia del resorte 1','Suma de residuos al cuadrado',[
             ('Perfil: reajustar c2',TEAL,[(r[0],r[3]) for r in data]),
             ('Corte: fijar c2=1,99',OCHRE,[(r[0],r[4]) for r in data])],
             'PERFIL ≠ CORTE CON PARÁMETRO FIJO',legend=False)
        s.note('Todas las parejas positivas c1+c2=2,99 reproducen los mismos datos: R² = 0,99776.')
    return s.svg()


def assert_near(x,y,comment,tol=1e-6):
    if not isclose(x,y,abs_tol=tol,rel_tol=tol):raise AssertionError(f'{comment}: {x} != {y}')


def audit_model(D):
    assert len(D)==10
    rows={301:5,302:101,303:5,304:104,305:101,306:3,307:11,308:2,309:91,310:98}
    assert sum(rows.values())==521
    for n,(headers,values) in D.items():
        assert len(values)==rows[n] and len(headers)>=4
        assert len(set(headers))==len(headers)
        assert all(len(row)==len(headers) for row in values)
    assert [r[3] for r in D[301][1]]==[3]*5
    assert_near(D[302][1][80][1],comb(10,8)*.8**8*.2**2,'C302 binomial observed')
    assert_near(D[302][1][80][2],1,'C302 n10 likelihood maximum')
    assert_near(D[302][1][80][3],1,'C302 n100 likelihood maximum')
    assert_near(sum(r[4]**2 for r in D[303][1]),.1466,'C303 original SSE',tol=.003)
    assert_near(sum(r[5]**2 for r in D[303][1]),.2605,'C303 log model original SSE',tol=.003)
    low=D[304][1][:4]
    assert_near(low[-1][1],.990099,'C304 reference low signal',tol=.0001)
    assert_near(low[-1][3],.0049256,'C304 near-equivalent curves',tol=.00002)
    assert_near(sum(.01*r[2] for r in D[305][1]),1,'C305 beta posterior integral',tol=.002)
    assert_near(D[305][1][50][1],1.5,'C305 prior Beta(2,2)')
    assert_near(D[306][1][0][3],0,'C306 lasso exact zero')
    assert_near(D[306][1][2][2],2,'C306 ridge strongest')
    assert_near(sum(r[1] for r in D[307][1]),1,'C307 predictive normalized')
    assert_near(D[307][1][0][4],.0280,'C307 count >=5 posterior predictive',tol=.002)
    assert_near(D[308][1][0][2],284,'C308 linear AIC')
    assert_near(D[308][1][1][2],280,'C308 cubic AIC')
    assert_near(D[308][1][0][3],289.21034,'C308 linear BIC',tol=.00003)
    assert_near(D[309][1][15][1],62.5,'C309 Fisher p=.20')
    assert_near(D[309][1][15][2],1.6,'C309 Fisher logit p=.20')
    assert all(abs(r[2]-2.99)<1e-10 and abs(r[3]-.097)<1e-10 for r in D[310][1])
    assert_near(D[310][1][0][4],.097+30*(.03+1.99-2.99)**2,'C310 fixed other spring')


def byte_outputs():
    D=make_data();audit_model(D);out={}
    for n in range(301,311):
        head,values=D[n]
        out[f'repro/c301-310/c{n}.csv']=csv_bytes(head,values)
        fig=charts(n,values)
        root=ElementTree.fromstring(fig);ns='{http://www.w3.org/2000/svg}'
        assert root.find(ns+'title') is not None and root.find(ns+'desc') is not None
        out[f'figures/complexa-c{n}-fig01.svg']=fig
    return out


def make_archive(files,readme,script):
    entries=dict(files)
    entries['repro/c301-310/README.md']=readme
    entries['repro/c301-310/make_figures.py']=script
    stream=BytesIO()
    with zipfile.ZipFile(stream,'w',compression=zipfile.ZIP_STORED,allowZip64=False) as archive:
        for name in sorted(entries):
            info=zipfile.ZipInfo(name,date_time=(2026,10,27,12,0,0))
            info.compress_type=zipfile.ZIP_STORED;info.create_system=3;info.external_attr=0o100644<<16
            archive.writestr(info,entries[name])
    return stream.getvalue()


def outputs():
    d=byte_outputs();here=Path(__file__)
    readme=here.with_name('README.md').read_bytes();script=here.read_bytes()
    archive=make_archive(d,readme,script)
    d['zip-data/complexa-c301-c310-reproducible.zip.b64']=base64.b64encode(archive)+b'\n'
    return d,archive


def run():
    ap=argparse.ArgumentParser(description=__doc__)
    gr=ap.add_mutually_exclusive_group()
    gr.add_argument('--out',type=Path,help='Output path content/supplementary')
    gr.add_argument('--check',type=Path,help='Check exact deterministic bytes')
    opts=ap.parse_args();dest=opts.out or opts.check or ROOT
    files,archive=outputs()
    for rel,blob in sorted(files.items()):
        target=dest/rel
        if opts.check:
            if not target.is_file() or target.read_bytes()!=blob:raise SystemExit(f'FAIL reproduction: {rel}')
        else:
            target.parent.mkdir(parents=True,exist_ok=True)
            target.write_bytes(blob)
    assert len(files)==21
    with zipfile.ZipFile(BytesIO(archive)) as z:
        assert len(z.namelist())==22 and z.testzip() is None
    print(f"{'CHECK' if opts.check else 'GENERATED'}: 10 SVG + 10 CSV + reproducible ZIP SHA-256 {sha256(archive).hexdigest()}")


if __name__=='__main__':run()