#!/usr/bin/env python3
"""COMPLEXA C321-C330 — reproducible exact scientific figures.

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

C321-C329: original deterministic mathematical examples. C330: published
retrospective summary statistics reproduced from the existing chapter only;
ECG raw signals are NOT archived here, so no raw-data reproduction is claimed.
"""
from __future__ import annotations
import argparse
import math
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]
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.ordinal=1; 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}.{self.ordinal}. {X(title)}</title>')
        self.a.append(f'<desc id="des{n}">{X(desc)}. {"Métricas publicadas en el capítulo sobre ECG histórico; datos crudos no incluidos en este suplemento." if n==330 else "Modelo matemático didáctico: no son observaciones empíricas."}</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,"RESUMEN DE RESULTADOS EMPÍRICOS" if n==330 else "MODELO MATEMÁTICO / DATOS DIDÁCTICOS",12,MUTED,700)
        self.text(1138,718,f"complexa.science  ·  {n}.{self.ordinal}",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")







NAMES={
 321:("Observación, filtrado y suavizado", "Un estado térmico reconstruido a partir de dos lecturas ruidosas"),
 322:("Dos sensores no equivalen a dos pruebas independientes", "Efecto de la correlación de errores en una actualización de Kalman"),
 323:("Una observación cuadrática conserva dos posibilidades", "Posterior bimodal exacta frente a una previa normal unimodal"),
 324:("Las partículas pueden concentrar sus pesos", "Tres estados: actualización exacta, cuatro partículas y remuestreo"),
 325:("Filtrar y suavizar una cadena oculta", "Cuatro lecturas: alta, baja, baja, alta; estados A y B"),
 326:("La separación temporal modifica el pronóstico", "Proceso Ornstein–Uhlenbeck: transición exacta y paso de Euler"),
 327:("Observar selectivamente puede desplazar una media", "Mil unidades, dos grupos: casos completos y ponderación inversa"),
 328:("La población modifica la estimación individual", "Tres trayectorias y cuatro hipótesis jerárquicas de persistencia"),
 329:("Mejoran ambos grupos y cae el promedio", "Dos semanas y distinta composición: paradoja de Simpson"),
 330:("Estados latentes: distintos criterios, distintas ventajas", "Métricas retrospectivas publicadas para ECG real, sin nueva señal cruda"),
}

# Deterministic calculations, no random draws, explicit model and units.
def logbeta(a,b):
    return math.lgamma(a)+math.lgamma(b)-math.lgamma(a+b)

def make_data():
    D={}
    # C321. Random-walk state with process variance 0.25, reading variance 1.
    state_m=20.; state_P=1.; rows=[]; filtered=[]
    for t,y in enumerate((22.,19.),1):
        mp=state_m; Pp=state_P+.25
        gain=Pp/(Pp+1.)
        state_m=mp+gain*(y-mp);state_P=(1-gain)*Pp
        rows.append([t,y,mp,Pp,gain,state_m,state_P,None,None]);filtered.append((state_m,state_P,Pp))
    j=filtered[0][1]/filtered[1][2]
    rows[0][-2]=filtered[0][0]+j*(filtered[1][0]-rows[1][2]);rows[0][-1]=filtered[0][1]+j*j*(filtered[1][1]-rows[1][3])
    rows[1][-2]=rows[1][5];rows[1][-1]=rows[1][6]
    assert abs(rows[0][5]-21.1111111111)<1e-9
    assert abs(rows[0][-2]-20.4615384615)<1e-9
    assert abs(rows[1][5]-20.1692307692)<1e-9
    D[321]=rows

    # C322. Two correlated equal-variance sensors, both report 2, prior N(0,1).
    rows=[]
    for k in range(61):
        rho=-.5+k*.025
        if rho>1.:rho=1.
        posterior_var=(1+rho)/(3+rho)
        posterior_mean=4/(3+rho)
        gain_each=1/(3+rho)
        rows.append([rho,posterior_mean,posterior_var,gain_each])
    s=min(rows,key=lambda r:abs(r[0]-.8))
    assert abs(s[1]-20/19)<1e-12 and abs(s[2]-9/19)<1e-12
    assert abs(rows[20][1]-4/3)<1e-12
    D[322]=rows

    # C323. Posterior p(x|y=4) ∝ N(x;0,1) N(4;x²,0.25).
    # Integrate using composite Simpson over [-4,4], 600 subintervals.
    rows=[]
    for k in range(601):
        x=-4+8*k/600
        prior=math.exp(-x*x/2)/math.sqrt(2*math.pi)
        likelihood=math.exp(-.5*(4-x*x)**2/.25)/math.sqrt(2*math.pi*.25)
        rows.append([x,prior,likelihood,prior*likelihood,0.])
    dx=8/600
    def simpson(cols):
        return dx/3*(cols[0]+cols[-1]+sum((4 if i%2 else 2)*v for i,v in enumerate(cols[1:-1],1)))
    z=simpson([r[3] for r in rows]);
    for r in rows:r[-1]=r[3]/z
    mean=simpson([r[0]*r[4] for r in rows]);
    var=simpson([(r[0]-mean)**2*r[4] for r in rows])
    assert abs(mean)<1e-9 and abs(var-3.841876)<2e-5,(mean,var)
    D[323]=rows

    # C324. Prior [-1,0,1] = [1,2,1]/4, likelihood [0.1,0.2,0.8].
    vals=(-1,0,1);prior=(.25,.5,.25);likelihood=(.1,.2,.8)
    post=tuple(a*b/sum(x*y for x,y in zip(prior,likelihood)) for a,b in zip(prior,likelihood))
    D[324]=[[vals[i],prior[i],likelihood[i],post[i],[0,2,2][i]] for i in range(3)]
    weights=[1/13,2/13,2/13,8/13]
    ess=1/sum(x*x for x in weights)
    assert abs(ess-169/73)<1e-14
    assert all(abs(post[i]-n/13)<1e-12 for i,n in enumerate((1,4,8)))
    assert abs(sum(vals[i]*post[i] for i in range(3))-7/13)<1e-14

    # C325. HMM, exact forward and backward recursions.
    A=((.9,.1),(.2,.8)); emission=(.2,.8);observations=(1,0,0,1)
    raw=[];alpha=(.6,.4)
    for i,o in enumerate(observations):
        predicted=alpha if i==0 else (alpha[0]*A[0][0]+alpha[1]*A[1][0],alpha[0]*A[0][1]+alpha[1]*A[1][1]);
        w=[predicted[j]*(emission[j] if o else 1-emission[j]) for j in range(2)]
        scale=sum(w);posterior=(w[0]/scale,w[1]/scale)
        raw.append([i+1,'alta' if o else 'baja',predicted[1],posterior[1],scale])
        alpha=posterior
    beta=[(1.,1.)]*4
    for i in range(2,-1,-1):
        o=observations[i+1];nxt=beta[i+1]
        b=[sum(A[j][k]*(emission[k] if o else 1-emission[k])*nxt[k] for k in range(2)) for j in range(2)]
        scale=max(b);beta[i]=(b[0]/scale,b[1]/scale)
    D[325]=[]
    for i,r in enumerate(raw):
        a=r[3];b=beta[i]
        sm=a*b[1]/((1-a)*b[0]+a*b[1])
        D[325].append(r+[sm]);
    assert abs(math.prod(row[4] for row in raw)-.03569152)<1e-10
    assert abs(D[325][-1][-1]-D[325][-1][3])<1e-14

    # C326. Continuous-time OU E[x_t]=10+4*e^(-0.5 t), Var=4(1-e^-t).
    D[326]=[]
    for i in range(91):
        h=i/30
        exact_mu=10+4*math.exp(-.5*h)
        exact_var=4*(1-math.exp(-h))
        euler_mu=14-2*h
        euler_var=4*h
        D[326].append([h,exact_mu,exact_var,euler_mu,euler_var])
    assert abs(D[326][90][1]-10.8925206405)<1e-9
    assert abs(D[326][90][2]-3.8008517265)<1e-9
    assert D[326][90][3]==8 and D[326][90][4]==12

    # C327. 500 per group, response probabilities .8/.2, means 2/8.
    D[327]=[[1,500,400,2.,.8,1.25],[2,500,100,8.,.2,5.]]
    naive=sum(r[2]*r[3] for r in D[327])/sum(r[2] for r in D[327])
    corrected=sum(r[2]*r[-1]*r[3] for r in D[327])/sum(r[2]*r[-1] for r in D[327])
    assert abs(naive-3.2)<1e-14 and abs(corrected-5.)<1e-14

    # C328. Four equiprobable population mean/concentration configurations.
    configs=((.55,4.),(.55,20.),(.8,4.),(.8,20.))
    ntr=((10,2),(4,0),(6,6))
    logweights=[]
    for m,c in configs:
        lw=sum(logbeta(m*c+s,(1-m)*c+f)-logbeta(m*c,(1-m)*c) for s,f in ntr)
        logweights.append(lw)
    ref=max(logweights);W=[math.exp(v-ref) for v in logweights];W=[w/sum(W) for w in W]
    D[328]=[]
    for (m,c),w in zip(configs,W):
        means=[(m*c+s)/(c+s+f) for s,f in ntr]
        D[328].append([m,c,w,*means])
    posterior=[sum(r[2]*r[j+3] for r in D[328]) for j in range(3)]
    assert all(abs(a-b)<.001 for a,b in zip(posterior,(.789,.822,.597))),posterior

    # C329. Two groups: rates improve but mix shifts to harder cases.
    D[329]=[['Sencillas',80,72,20,19,.9,.95],['Complejas',20,10,80,48,.5,.6]]
    total0=sum(r[2] for r in D[329])/sum(r[1] for r in D[329]);total1=sum(r[4] for r in D[329])/sum(r[3] for r in D[329])
    standard0=sum(r[5] for r in D[329])/2;standard1=sum(r[6] for r in D[329])/2
    assert all(abs(a-b)<1e-12 for a,b in zip((total0,total1,standard0,standard1),(.82,.67,.7,.775)))
    D[329].append(['Total',100,82,100,67,.82,.67])

    # C330. These four summary metrics are transcribed from the EXISTING,
    # published empirical case (ECG MIT-BIH record 208, 60 held-out seconds).
    # Raw 108,000-sample ECG and full model-fitting code are not in this repo;
    # therefore the figure is a transparent reinterpretation of reported metrics,
    # not independent reanalysis of the ECG or proof of physiological states.
    D[330]=[
      ['Gaussiana',60,-.07335,.21353,.16245],
      ['AR(1)',60,-.09261,.25523,.19798],
      ['Mezcla-3',60,-.06386,.21353,.16245],
      ['HMM-3',60,-.04421,.23386,.17814],
    ]
    assert D[330][-1][2]>max(r[2] for r in D[330][:-1])
    assert abs(60*D[330][-1][2]+2.6526)<.001
    return D

def grouped(s,area,heading,labels,sets,ymax,unit='Proporción',digits=2):
    """Accessible grouped comparison of values [0,ymax]; groups are categories."""
    x,y,w,h=area
    s.rect(x,y,w,h,'#FFFFFF',12,GRID);s.text(x+20,y+34,heading,17,NAVY,700)
    left=x+74;right=x+w-23;top=y+108;bottom=y+h-65
    def yy(v):return bottom-v/ymax*(bottom-top)
    for q in range(5):
        val=ymax*q/4;cy=yy(val)
        s.line(left,cy,right,cy,GRID);s.text(left-10,cy+4,f'{val:.2f}' if ymax<5 else f'{val:g}',12,MUTED,400,'end')
    s.line(left,bottom,right,bottom,NAVY,2)
    lx=x+20
    for name,color,vals in sets:
        s.line(lx,y+68,lx+17,y+68,color,5);s.text(lx+22,y+73,name,13,INK,500)
        lx+=len(name)*8+65
    n=len(labels);slot=(right-left)/n;groupw=min(slot*.79,120);gap=5;barw=(groupw-gap*(len(sets)-1))/len(sets)
    for i,lab in enumerate(labels):
        middle=left+slot*(i+.5)
        for j,(name,color,vals) in enumerate(sets):
            val=vals[i];bx=middle-groupw/2+j*(barw+gap)
            s.rect(bx,yy(val),barw,bottom-yy(val),color,3)
            if groupw>=110 and val/ymax>.09:
                fmt=f'{val:.{digits}f}'
                s.text(bx+barw/2,yy(val)-9,fmt,12,color,700,'middle')
        s.text(middle,bottom+24,lab,12,INK,600,'middle')
    s.text((left+right)/2,y+h-18,unit,13,MUTED,500,'middle')


def charts(n,D):
    title,desc=NAMES[n];s=SVG(n,title,desc);rows=D[n]
    if n==321:
        plot(s,(55,175,530,452),1,2,18,22.5,[1,2],[18,19,20,21,22],
             'Minuto','Temperatura',[
             ('Lectura',MUTED,[(r[0],r[1]) for r in rows]),
             ('Filtrado',TEAL,[(r[0],r[5]) for r in rows]),
             ('Suavizado',OCHRE,[(r[0],r[7]) for r in rows])],
             'DOS OBSERVACIONES',legend=True)
        plot(s,(609,175,536,452),1,2,0,1.4,[1,2],[0,.25,.5,.75,1,1.25],
             'Minuto','Varianza',[
             ('Predicha',MUTED,[(r[0],r[3]) for r in rows]),
             ('Filtrada',TEAL,[(r[0],r[6]) for r in rows]),
             ('Suavizada',OCHRE,[(r[0],r[8]) for r in rows])],
             'DIFERENTES CONJUNTOS DE DATOS')
        s.note('La lectura del minuto 2 revisa la temperatura del minuto 1: 21,111 → 20,462.')
    elif n==322:
        plot(s,(55,175,531,453),-.5,1,0,1.8,[-.5,0,.5,1],[0,.4,.8,1.2,1.6],
             'Correlación de errores','Media posterior',[('Media',TEAL,[(r[0],r[1]) for r in rows])],
             'AMBOS SENSORES INFORMAN 2',highlight=[(.8,20/19,'rho=0,8',TEAL)])
        plot(s,(604,175,540,453),-.5,1,0,.65,[-.5,0,.5,1],[0,.15,.3,.45,.6],
             'Correlación de errores','Varianza posterior',[('Varianza',OCHRE,[(r[0],r[2]) for r in rows])],
             'INFORMACIÓN COMPARTIDA',highlight=[(.8,9/19,'9/19',OCHRE)])
        s.note('Con correlación 0,8, media = 20/19 y varianza = 9/19; con independencia: 4/3 y 1/3.')
    elif n==323:
        plot(s,(55,175,1090,451),-4,4,0,1.65,[-4,-3,-2,-1,0,1,2,3,4],[0,.4,.8,1.2,1.6],
             'Estado latente x','Densidad normalizada',[
             ('Previa N(0,1)',MUTED,[(r[0],r[1]) for r in rows]),
             ('Posterior | observación x²≈4',PURPLE,[(r[0],r[4]) for r in rows])],
             'DOS MODOS: LA MEDICIÓN NO IDENTIFICA EL SIGNO')
        s.note('Integración numérica: media posterior ≈ 0; varianza ≈ 3,841876 (superior a 1).')
    elif n==324:
        labels=['Estado −1','Estado 0','Estado +1']
        grouped(s,(55,175,530,451),'ESTADOS ANTES Y DESPUÉS DE OBSERVAR',labels,[
            ('Previa',MUTED,[r[1] for r in rows]),('Posterior',TEAL,[r[3] for r in rows])],.75,'Probabilidad',3)
        grouped(s,(610,175,535,451),'PARTÍCULAS ANTES Y DESPUÉS',labels,[
            ('Inicial',PURPLE,[1,2,1]),('Tras remuestreo',OCHRE,[r[4] for r in rows])],3.4,'Recuento de representantes',0)
        s.note('N = 4; ESS ponderado = 169/73 ≈ 2,315. El remuestreo elimina el estado −1.')
    elif n==325:
        plot(s,(55,175,1090,451),1,4,0,1,[1,2,3,4],[0,.2,.4,.6,.8,1],
             'Instante: alta, baja, baja, alta','Probabilidad del estado B',[
             ('Predicha',MUTED,[(r[0],r[2]) for r in rows]),
             ('Filtrada',TEAL,[(r[0],r[3]) for r in rows]),
             ('Suavizada',PURPLE,[(r[0],r[5]) for r in rows])],
             'USAR EL FUTURO CAMBIA LA RECONSTRUCCIÓN DEL PASADO')
        s.note('Secuencia observada alta-baja-baja-alta: verosimilitud exacta 0,03569152.')
    elif n==326:
        plot(s,(55,175,530,451),0,3,7.5,14.5,[0,.5,1,1.5,2,2.5,3],[8,10,12,14],
             'Intervalo (horas)','Media del estado',[
             ('Exacta',TEAL,[(r[0],r[1]) for r in rows]),
             ('Euler: un paso',OCHRE,[(r[0],r[3]) for r in rows])],
             'RETORNO HACIA EL NIVEL 10')
        plot(s,(609,175,536,451),0,3,0,12.8,[0,.5,1,1.5,2,2.5,3],[0,3,6,9,12],
             'Intervalo (horas)','Varianza predictiva',[
             ('Exacta',TEAL,[(r[0],r[2]) for r in rows]),
             ('Euler: un paso',OCHRE,[(r[0],r[4]) for r in rows])],
             'INCERTIDUMBRE NO LINEAL EN EL TIEMPO')
        s.note('A las 3 h: media exacta = 10,8925 y varianza = 3,8009; Euler = 8 y 12.')
    elif n==327:
        labels=['Grupo 1','Grupo 2']
        grouped(s,(55,175,530,451),'COMPOSICIÓN DE LA MUESTRA',labels,[
            ('Población',MUTED,[r[1] for r in rows]),
            ('Respuestas',TEAL,[r[2] for r in rows])],590,'Unidades esperadas',0)
        grouped(s,(610,175,535,451),'MEDIA SEGÚN EL PROCEDIMIENTO',
             ['Población','Observada','Pond. inversa'],[
             ('Media',OCHRE,[5.,3.2,5.])],6.2,'Media de la variable',2)
        s.note('Las respuestas dependen del grupo: 80 % y 20 %. La ponderación recupera la media 5.')
    elif n==328:
        labels=['0,55/4','0,55/20','0,80/4','0,80/20']
        grouped(s,(55,175,531,453),'CUATRO HIPÓTESIS POBLACIONALES',labels,
             [('Peso posterior',PURPLE,[r[2] for r in rows])],.43,'Media / concentración',3)
        observations=[10/12,1.,6/12]
        posterior=[sum(r[2]*r[j+3] for r in rows) for j in range(3)]
        grouped(s,(609,175,535,453),'EFECTO DE LA COMBINACIÓN PARCIAL',
             ['Persona A','Persona B','Persona C'],[
             ('Frecuencia observada',MUTED,observations),
             ('Media posterior',TEAL,posterior)],1.08,'Probabilidad de permanecer',3)
        s.note('Las medias posteriores personales son 0,789; 0,822 y 0,597: integración de 4 modelos.')
    elif n==329:
        labels=['Sencillas','Complejas']
        grouped(s,(55,175,529,453),'MEJORA EN AMBAS CATEGORÍAS',labels,[
             ('Semana 1',MUTED,[r[5] for r in rows[:2]]),
             ('Semana 2',TEAL,[r[6] for r in rows[:2]])],1.06,'Proporción resuelta',2)
        grouped(s,(611,175,534,453),'CAÍDA DEL TOTAL POR COMPOSICIÓN',
             ['Total real','Estand. 50/50'],[
             ('Semana 1',MUTED,[.82,.7]),
             ('Semana 2',OCHRE,[.67,.775])],1.05,'Proporción resuelta',3)
        s.note('Descomposición simétrica exacta: tasas +7,5 puntos; composición −22,5 puntos.')
    elif n==330:
        labs=[r[0] for r in rows]
        grouped(s,(55,175,531,452),'PÉRDIDA LOGARÍTMICA MEDIA',labs,[
            ('Menor = mejor',PURPLE,[-r[2] for r in rows])],.112,'Negativo de la logdensidad media',4)
        grouped(s,(610,175,535,452),'ERROR CUADRÁTICO DE PREDICCIÓN',labs,[
            ('RMSE; menor = mejor',OCHRE,[r[3] for r in rows])],.32,'Escala log(RMS)',4)
        s.note('Resumen de C330 (60 s de ECG): HMM mejora logdensidad; gaussiana y mezcla empatan en RMSE.')
    else: raise ValueError(n)
    return s.svg()

HEADERS={
 321:['minute','observed_temperature','predicted_state_mean','predicted_state_variance','kalman_gain','filtered_mean','filtered_variance','smoothed_mean','smoothed_variance'],
 322:['sensor_error_correlation','posterior_mean','posterior_variance','gain_per_sensor'],
 323:['latent_state','prior_pdf','likelihood_given_reading','unnormalized_posterior','normalized_posterior_pdf'],
 324:['state','prior_probability','observation_likelihood','exact_posterior_probability','systematic_resampled_count'],
 325:['time','observation','predicted_state_B_probability','filtered_state_B_probability','predictive_observation_probability','smoothed_state_B_probability'],
 326:['hours','exact_conditional_mean','exact_conditional_variance','one_step_euler_mean','one_step_euler_variance'],
 327:['group','original_population_count','respondents_count_expected','group_mean','response_probability','inverse_probability_weight'],
 328:['population_mean','concentration','model_posterior_weight','individual_A_posterior_mean','individual_B_posterior_mean','individual_C_posterior_mean'],
 329:['category','received_week1','resolved_week1','received_week2','resolved_week2','rate_week1','rate_week2'],
 330:['model','held_out_ecg_seconds','reported_mean_log_density','reported_RMSE_log_RMS','reported_MAE_log_RMS']
}


def models():
    d=make_data();files={}
    for n in d:
        files[f'figures/complexa-c{n}-fig01.svg']=charts(n,d)
        files[f'repro/c321-330/c{n}.csv']=csv_bytes(HEADERS[n],d[n])
    for n in d:
        fig=ElementTree.fromstring(files[f'figures/complexa-c{n}-fig01.svg'])
        assert fig.tag.endswith('svg')
        assert fig.find('{http://www.w3.org/2000/svg}title') is not None
        assert fig.find('{http://www.w3.org/2000/svg}desc') is not None
        assert len(d[n])>1
    return files,d


def archive(files,readme,script):
    entries=dict(files)
    entries['repro/c321-330/README.md']=readme
    entries['repro/c321-330/make_figures.py']=script
    buf=BytesIO()
    with zipfile.ZipFile(buf,'w',compression=zipfile.ZIP_STORED,allowZip64=False) as z:
        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
            z.writestr(info,entries[name])
    return buf.getvalue()


def outputs():
    files,D=models()
    path=Path(__file__)
    zipdata=archive(files,path.with_name('README.md').read_bytes(),path.read_bytes())
    files['zip-data/complexa-c321-c330-reproducible.zip.b64']=base64.b64encode(zipdata)+b'\n'
    return files,D,zipdata


def main():
    ap=argparse.ArgumentParser(description=__doc__)
    gp=ap.add_mutually_exclusive_group()
    gp.add_argument('--out',type=Path,help='Write files in content/supplementary')
    gp.add_argument('--check',type=Path,help='Recreate and compare byte for byte')
    args=ap.parse_args();dest=args.out or args.check or ROOT
    files,data,zipdata=outputs()
    for rel,blob in sorted(files.items()):
        target=dest/rel
        if args.check:
            if not target.exists() or target.read_bytes()!=blob:
                raise SystemExit(f'FAIL reproducibility: {rel}')
        else:
            target.parent.mkdir(parents=True,exist_ok=True)
            target.write_bytes(blob)
    with zipfile.ZipFile(BytesIO(zipdata)) as z:
        assert len(z.namelist())==22 and z.testzip() is None
    assert len(files)==21 and len(data)==10
    print(f'{"VERIFIED" if args.check else "GENERATED"}: 10 SVG, 10 CSV, 22 ZIP entries; rows {sum(len(x) for x in data.values())}; SHA-256 {sha256(zipdata).hexdigest()}')

if __name__=='__main__':main()