#!/usr/bin/env python3
"""COMPLEXA C291-C300: 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 = {
 291:("Un regulador amortigua, pero deja error", "Control proporcional frente a caudal sostenido"),
 292:("Tres respuestas para ocho perturbaciones", "Variedad requerida con tolerancia de una unidad"),
 293:("Más proveedores no eliminan el riesgo común", "Diversificación frente a correlación homogénea"),
 294:("El desacuerdo puede desaparecer de dos formas", "Consenso por contagio o fragmentación de contactos"),
 295:("El refuerzo cambia la distribución de resultados", "Urna de Pólya frente a diez extracciones independientes"),
 296:("Adaptarse al final no asegura viabilidad", "Control integral: exceso transitorio sobre el límite"),
 297:("La agregación puede ocultar estado relevante", "Agrupabilidad de una cadena de cuatro estados"),
 298:("Cien opiniones no equivalen a cien fuentes", "Precisión del promedio con correlación y sesgo"),
 299:("Una conexión nueva puede empeorar la red", "Equilibrio de Braess frente a coordinación óptima"),
 300:("Memoria completa y ventana reciente", "Doce sesiones y costos de preparación comparados"),
}
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")




# Every numerical record follows the chapter's mathematical rules and units.
# Results represent synthetic scenarios, analytic expectations, or numerical ODEs, not empirical observations.

def probs_update(prob,transition):
    nxt=[sum(prob[i]*transition[i][j] for i in range(len(prob))) for j in range(len(prob))]
    assert isclose(sum(nxt),1,abs_tol=1e-12) and min(nxt)>-1e-12
    return nxt


def make_data():
    D={}
    # C291: volume 50 initially; baseline extraction five and proportional controller gain .5.
    rows=[]
    for t in range(16):
        p7=54-4*(.5**t)
        p9=58-8*(.5**t)
        uncontrolled=50+2*t
        rows.append((t,p7,p9,uncontrolled,55,45,5+.5*(p7-50)))
    D[291]=(['minute','proportional_inflow_7_liters','proportional_inflow_9_liters','fixed_pump_inflow_7_liters','upper_band_liters','lower_band_liters','pump_flow_inflow_7_liters'],rows)
    # C292: all 8 perturbations with only three integer response magnitudes.
    chosen=[1,1,1,4,4,4,6,6]
    D[292]=(['disturbance_d','compensation_r','residual_d_minus_r','is_acceptable_abs_le_1','maximum_allowed_error'],[(d,chosen[d],d-chosen[d],int(abs(d-chosen[d])<=1),1) for d in range(8)])
    # C293: standard deviations for identically distributed producers (s.d. 0.2), equal weights.
    D[293]=(['provider_count','sd_independent_rho_0','sd_correlated_rho_0_2','sd_correlated_rho_0_75','sd_identical_rho_1','single_provider_sd'],[(n,.2/sqrt(n),.2*sqrt(.2+.8/n),.2*sqrt(.75+.25/n),.2,.2) for n in range(1,41)])
    # C294: two specified paths from the same six-vertex signed ring.
    rows=[]
    for t in range(4):
        active_split=[2,1,0,0][t];active_consensus=[2,2,2,0][t]
        plus_consensus=3+t
        rows.append((t,active_split/6,active_consensus/6,0,(2*plus_consensus-6)/6,active_split,active_consensus,3,plus_consensus,2 if t>=2 else 1,1))
    D[294]=(['event_index','discordant_edge_fraction_fragmentation','discordant_edge_fraction_consensus','magnetization_fragmentation','magnetization_consensus','discordant_edges_fragmentation','discordant_edges_consensus','plus_nodes_fragmentation','plus_nodes_consensus','components_fragmentation','components_consensus'],rows)
    # C295: uniform Polya urn count distribution vs binomial(10, 1/2).
    D[295]=(['number_A_in_ten_draws','probability_polya_urn','probability_fair_coin','distribution_difference','initial_balls_each_color'],[(k,1/11,comb(10,k)/1024,1/11-comb(10,k)/1024,1) for k in range(11)])
    # C296: analytic y and u after permanent step, a=b=k=1, r=2, d=1.
    rows=[]
    for j in range(201):
        t=j*.05;omega=sqrt(3)/2;f=exp(-t/2)
        y=2+2/sqrt(3)*f*sin(omega*t)
        u=1+f*(cos(omega*t)+sin(omega*t)/sqrt(3))
        rows.append((t,y,u,2,2.4,int(y>2.4)))
    D[296]=(['time','regulated_output_y','control_u','target_y','viability_ceiling_y','violates_limit'],rows)
    # C297: strong lumpability if rows for microstates 1,2 have equal totals into A.
    P=[ [.6,.2,.1,.1], [.1,.7,.15,.05], [.05,.25,.4,.3], [.2,.1,.1,.6] ]
    Q=[p[:] for p in P];Q[1]=[.1,.5,.35,.05]
    orig1=[1,0,0,0];orig2=[0,1,0,0];mod1=orig1[:];mod2=orig2[:]
    rows=[]
    for t in range(16):
        rows.append((t,sum(orig1[:2]),sum(orig2[:2]),sum(mod1[:2]),sum(mod2[:2]),.6))
        orig1=probs_update(orig1,P);orig2=probs_update(orig2,P)
        mod1=probs_update(mod1,Q);mod2=probs_update(mod2,Q)
    D[297]=(['step','pA_lumpable_start_micro1','pA_lumpable_start_micro2','pA_nonlumpable_start_micro1','pA_nonlumpable_start_micro2','lumpable_longrun_pA'],rows)
    # C298: each analyst error has bias 3, s.d. 10, pairwise correlation 0.2.
    D[298]=(['number_estimates','variance_uncorrelated','variance_rho_0_2','mse_uncorrelated_bias_3','mse_correlated_bias_3','effective_independent_n_at_rho_0_2','mse_correlated_asymptote'],[(n,100/n,100*(.2+.8/n),9+100/n,9+100*(.2+.8/n),n/(1+(n-1)*.2),29) for n in range(1,151)])
    # C299: exact Braess equilibria/optimum from 4000 travellers.
    D[299]=(['scenario','upper_route_people','lower_route_people','link_route_people','mean_time_minutes','aggregate_travel_minutes','upper_route_time','lower_route_time','link_route_time'],[
      ('without_link_equilibrium',2000,2000,0,65,260000,65,65,0),
      ('with_link_selfish_equilibrium',0,0,4000,80,320000,85,85,80),
      ('with_link_system_optimum',1750,1750,500,64.6875,258750,67.5,67.5,45)])
    # C300: exactly the published 12-session demand sequence and Bayesian rolling estimates.
    demands=[2,2,2,2,6,2,2,6,6,6,6,6]
    rows=[]
    for index,Demand in enumerate(demands):
        previous=demands[:index]
        window=previous[-4:]
        full=(1+sum(x==6 for x in previous))/(len(previous)+2)
        short=(1+sum(x==6 for x in window))/(len(window)+2)
        qf=6 if full>1/3 else 2;qw=6 if short>1/3 else 2
        costf=qf+3*max(0,Demand-qf);costw=qw+3*max(0,Demand-qw)
        rows.append((index+1,Demand,full,short,qf,qw,costf,costw,sum(x==6 for x in previous),len(previous)))
    D[300]=(['session','demand_boxes','probability_full_history_before_session','probability_rolling4_before_session','prepared_full_history','prepared_rolling4','cost_full_history','cost_rolling4','previous_high_demand_count','previous_total_sessions'],rows)
    return D


def charts(n,data):
    title,subtitle=NAMES[n];s=SVG(n,title,subtitle)
    if n==291:
        plot(s,(55,176,1090,450),0,15,45,84,[0,3,6,9,12,15],[45,50,55,60,70,80],
             'Minutos','Volumen del depósito (L)',[
             ('Bomba fija · entrada 7',PURPLE,[(r[0],r[3]) for r in data]),
             ('Control P · entrada 7',TEAL,[(r[0],r[1]) for r in data]),
             ('Control P · entrada 9',OCHRE,[(r[0],r[2]) for r in data]),
             ('Límite superior (55 L)',NAVY,[(0,55),(15,55)])],
             'GANANCIA 0,5 / REFERENCIA 50 L',legend=False)
        for i,(txt,c) in enumerate([('Bomba fija, q=7',PURPLE),('Control P, q=7',TEAL),('Control P, q=9',OCHRE),('Banda 55 L',NAVY)]):
            x=95+(i%2)*310;y=231+(i//2)*28
            s.line(x,y-5,x+22,y-5,c,3.4);s.text(x+28,y,txt,13,INK,600)
        s.note('El control proporcional es estable con q=9, pero su equilibrio de 58 L rebasa la banda.')
    elif n==292:
        plot(s,(55,176,540,450),0,7,0,7.4,[0,1,2,3,4,5,6,7],[0,1,2,3,4,5,6,7],
             'Perturbación d','Acción r', [('Respuesta elegida',TEAL,[(r[0],r[1]) for r in data])],
             'POLÍTICA: 1, 4, 6',legend=False)
        plot(s,(606,176,539,450),0,7,-1.3,1.3,[0,1,2,3,4,5,6,7],[-1,0,1],
             'Perturbación d','Residuo e=d−r', [('Compensación',OCHRE,[(r[0],r[2]) for r in data])],
             'TODOS LOS RESIDUOS ADMISIBLES',legend=False)
        s.note('Tres acciones bien situadas cubren los ocho casos con |e| ≤ 1; el número solo no basta.')
    elif n==293:
        plot(s,(55,176,1090,450),1,40,0,.205,[1,5,10,20,30,40],[0,.05,.10,.15,.20],
             'Número de proveedores n','Desviación típica total',[
             ('ρ=0',TEAL,[(r[0],r[1]) for r in data]),
             ('ρ=0,20',OCHRE,[(r[0],r[2]) for r in data]),
             ('ρ=0,75',PURPLE,[(r[0],r[3]) for r in data])],
             'IGUAL PESO, DESVIACIÓN INDIVIDUAL 0,2')
        s.note('Diversificar reduce el riesgo idiosincrásico, pero no elimina la fluctuación compartida.')
    elif n==294:
        plot(s,(55,176,539,450),0,3,0,.40,[0,1,2,3],[0,.1,.2,.3,.4],
             'Eventos de cambio','Fracción de enlaces discordantes',[
             ('Fragmentación',PURPLE,[(r[0],r[1]) for r in data]),
             ('Consenso',TEAL,[(r[0],r[2]) for r in data])],
             'MISMA MEDIDA FINAL = 0')
        plot(s,(606,176,539,450),0,3,-.06,1.05,[0,1,2,3],[0,.25,.5,.75,1],
             'Eventos de cambio','Magnetización de opiniones',[
             ('Fragmentación',PURPLE,[(r[0],r[3]) for r in data]),
             ('Consenso',TEAL,[(r[0],r[4]) for r in data])],
             'DESENLACE NO EQUIVALENTE')
        s.note('Cero desacuerdo local puede significar dos comunidades separadas o acuerdo de toda la red.')
    elif n==295:
        plot(s,(55,176,1090,450),0,10,0,.29,[0,2,4,6,8,10],[0,.05,.10,.15,.20,.25],
             'Extracciones de A entre diez','Probabilidad exacta',[
             ('Urna de Pólya',TEAL,[(r[0],r[1]) for r in data]),
             ('Moneda independiente',OCHRE,[(r[0],r[2]) for r in data])],
             'CONTEO EXACTO SIN SIMULACIÓN')
        s.note('Bajo refuerzo todos los conteos K=0…10 son equiprobables; una moneda concentra K cerca de 5.')
    elif n==296:
        plot(s,(55,176,1090,450),0,10,1.75,2.65,[0,2,4,6,8,10],[1.8,2,2.2,2.4,2.6],
             'Tiempo después de la perturbación','Salida regulada y',[
             ('Respuesta integral',TEAL,[(r[0],r[1]) for r in data]),
             ('Límite admisible 2,4',OCHRE,[(0,2.4),(10,2.4)]),
             ('Objetivo 2,0',NAVY,[(0,2),(10,2)])],
             'ADAPTACIÓN FINAL ≠ VIABILIDAD TRANSITORIA',[(2*pi/(3*sqrt(3)),2.54629,'pico 2,546',TEAL)],legend=False)
        for i,(txt,c) in enumerate([('Salida',TEAL),('Límite 2,4',OCHRE),('Objetivo 2',NAVY)]):
            x=98+i*264;s.line(x,236,x+20,236,c,3.2);s.text(x+28,241,txt,13,INK,600)
        s.note('El sistema vuelve a y=2, pero supera transitoriamente el límite de 2,4.')
    elif n==297:
        plot(s,(55,176,539,450),0,15,0,1.02,[0,3,6,9,12,15],[0,.2,.4,.6,.8,1],
             'Transiciones','Probabilidad de estar en A',[
             ('Inicio en 1',TEAL,[(r[0],r[1]) for r in data]),
             ('Inicio en 2',OCHRE,[(r[0],r[2]) for r in data])],
             'AGREGACIÓN EXACTA')
        plot(s,(606,176,539,450),0,15,0,1.02,[0,3,6,9,12,15],[0,.2,.4,.6,.8,1],
             'Transiciones','Probabilidad de estar en A',[
             ('Inicio en 1',TEAL,[(r[0],r[3]) for r in data]),
             ('Inicio en 2',OCHRE,[(r[0],r[4]) for r in data])],
             'UNA FILA ALTERADA ROMPE EL CIERRE')
        s.note('Los mismos estados macro A iniciales pueden predecir destinos distintos si importa el microestado.')
    elif n==298:
        plot(s,(55,176,539,450),1,150,0,120,[1,20,50,100,150],[0,20,40,60,80,100,120],
             'Número de personas','Error cuadrático esperado',[
             ('Errores independientes',TEAL,[(r[0],r[3]) for r in data]),
             ('Correlación 0,20',OCHRE,[(r[0],r[4]) for r in data])],
             'MISMO SESGO = 3')
        plot(s,(606,176,539,450),1,150,0,155,[1,20,50,100,150],[0,25,50,75,100,125,150],
             'Número de personas','Número independiente equivalente',[
             ('ρ=0,20',PURPLE,[(r[0],r[5]) for r in data]),
             ('ρ=0',TEAL,[(r[0],r[0]) for r in data])],
             'TAMAÑO INFORMATIVO')
        s.note('Con 100 respuestas y correlación 0,2, el tamaño efectivo equivale a unas 4,81 respuestas.')
    elif n==299:
        labels=['Sin enlace','Enlace: egoísta','Enlace: óptimo']
        bars=[(labels[i],r[4],[TEAL,OCHRE,PURPLE][i]) for i,r in enumerate(data)]
        bar_panel(s,(55,176,539,450),'TIEMPO MEDIO DE VIAJE',bars,90,'Minutos por usuario')
        x,y,w,h=(606,176,539,450);s.rect(x,y,w,h,'#FFFFFF',12,GRID)
        s.text(x+20,y+33,'REPARTO DE 4.000 VIAJES',17,NAVY,700)
        left=x+54;right=x+w-33;top=y+92;bottom=y+h-66
        for i,row in enumerate(data):
            yy=top+38+i*92
            s.text(left,yy-10,labels[i],14,INK,600)
            cx=left
            for val,col in zip(row[1:4],(TEAL,PURPLE,OCHRE)):
                width=(right-left)*val/4000
                if width>0:s.rect(cx,yy+2,width,30,col,0)
                cx+=width
        for j,(lab,col) in enumerate([('Superior',TEAL),('Inferior',PURPLE),('Enlace U–V',OCHRE)]):
            lx=left+j*150;s.rect(lx,bottom+14,12,12,col,2);s.text(lx+20,bottom+25,lab,12,INK,500)
        s.note('La conexión permite un óptimo mejor que antes, pero el equilibrio individual tarda 80 min.')
    elif n==300:
        plot(s,(55,176,539,450),1,12,0,1,[1,3,5,7,9,12],[0,.2,.4,.6,.8,1],
             'Sesión','Probabilidad predictiva alta',[
             ('Historia completa',TEAL,[(r[0],r[2]) for r in data]),
             ('Ventana de cuatro',OCHRE,[(r[0],r[3]) for r in data]),
             ('Umbral de decisión',NAVY,[(1,1/3),(12,1/3)])],
             'PRONÓSTICO PREVIO A CADA SESIÓN',legend=False)
        for i,(lab,col) in enumerate([('Historia',TEAL),('Ventana 4',OCHRE),('Umbral',NAVY)]):
            x=118+i*135;s.line(x,244,x+20,244,col,3);s.text(x+26,248,lab,12,INK,500)
        bar_panel(s,(606,176,539,450),'COSTO EN 12 SESIONES',[
           ('Todo historial',sum(r[6] for r in data),TEAL),
           ('Ventana 4',sum(r[7] for r in data),OCHRE),
           ('Siempre 6',72,PURPLE),
           ('Siempre 2',sum(2+3*max(0,r[1]-2) for r in data),NAVY)],110,'Unidades de costo didáctico')
        s.note('Las previsiones son anteriores al resultado; la ventana reacciona antes al cambio de demanda.')
    return s.svg()


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


def audit_model(D):
    assert_near(D[291][1][2][1],53,'C291 volume at minute 2')
    assert_near(D[291][1][5][1],53.875,'C291 volume after 5')
    assert_near(D[291][1][15][2],58-8/(2**15),'C291 under high inflow')
    assert [r[1] for r in D[292][1]]==[1,1,1,4,4,4,6,6]
    assert all(abs(r[2])<=1 for r in D[292][1]); assert {r[1] for r in D[292][1]}=={1,4,6}
    for n,sd in [(4,.1),(16,.05)]:assert_near(D[293][1][n-1][1],sd,'C293 independent producers')
    assert_near(D[293][1][15][3],.175,'C293 correlated suppliers sd',tol=.001)
    assert [r[5] for r in D[294][1]]==[2,1,0,0]
    assert [r[6] for r in D[294][1]]==[2,2,2,0]
    assert [r[4] for r in D[294][1]]==[0,1/3,2/3,1]
    assert_near(sum(r[1] for r in D[295][1]),1,'C295 Polya normalized')
    assert_near(sum(r[2] for r in D[295][1]),1,'C295 coin normalized')
    assert_near(D[295][1][10][1],1/11,'C295 tail count')
    assert_near(D[295][1][10][2],1/1024,'C295 binomial tail')
    ma=max(r[1] for r in D[296][1]);assert_near(ma,2.5463,'C296 transient max',tol=.001)
    assert D[296][1][0][5]==0 and any(r[5]==1 for r in D[296][1])
    assert_near(D[297][1][1][1],.8,'C297 lumpable from micro1')
    assert_near(D[297][1][1][2],.8,'C297 lumpable from micro2')
    assert_near(D[297][1][1][4],.6,'C297 broken lumpability')
    assert_near(D[297][1][5][1],.6125,'C297 from A=1 fifth iteration',tol=.0000001)
    assert_near(D[298][1][4][2],36,'C298 correlated variance n=5')
    assert_near(D[298][1][19][4],33,'C298 correlated mse n20')
    assert_near(D[298][1][99][5],100/20.8,'C298 n_eff n100')
    assert_near(D[299][1][2][4],64.6875,'C299 optimal time')
    assert_near(D[299][1][1][5],320000,'C299 selfish aggregate time')
    assert_near(sum(r[6] for r in D[300][1]),76,'C300 full history losses')
    assert_near(sum(r[7] for r in D[300][1]),68,'C300 window 4 losses')
    assert_near(D[300][1][8][2],.3,'C300 n9 full predictive')
    assert_near(D[300][1][8][3],.5,'C300 n9 window predictive')
    for n,(header,rows) in D.items():
        assert len(header)>=4 and rows and all(len(r)==len(header) for r in rows),f'C{n} CSV malformed'
    assert len(D)==10


def byte_outputs():
    D=make_data();audit_model(D);out={}
    for n in range(291,301):
        header,rows=D[n]
        out[f'repro/c291-300/c{n}.csv']=csv_bytes(header,rows)
        fig=charts(n,rows)
        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/c291-300/README.md']=readme
    entries['repro/c291-300/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,26,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-c291-c300-reproducible.zip.b64']=base64.b64encode(archive)+b'\n'
    return d,archive


def run():
    ap=argparse.ArgumentParser(description=__doc__)
    args=ap.add_mutually_exclusive_group()
    args.add_argument('--out',type=Path,help='Directory content/supplementary to write')
    args.add_argument('--check',type=Path,help='Directory to verify generated 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()