#!/usr/bin/env python3
"""COMPLEXA C231–C240: source-derived computational experiments, Monte Carlo, variance and sensitivity SVG/CSV illustrations.

Pure Python 3.12+ standard library; deterministic charts, numerical audits and ZIP.
Run: python make_figures.py --out /path/to/content/supplementary
     python make_figures.py --check /path/to/content/supplementary
Every panel derives from the manuscript-defined examples, never unreported observations.
"""
from __future__ import annotations
import argparse
import csv
import hashlib
import io
import math
from pathlib import Path
from xml.sax.saxutils import escape
from zipfile import ZipFile, ZipInfo, ZIP_STORED

BASE = Path(__file__).resolve().parent
TEAL = '#087f8c'; BLUE = '#326d98'; ORANGE = '#c26c38'; PLUM = '#754c96'
INK = '#183348'; MUTED = '#566c78'; GRID = '#dae4eb'; BG = '#ffffff'; LT = '#f3f7fa'; GREEN='#2d7865'
COLORS=[BLUE,TEAL,ORANGE,PLUM,GREEN]
W,H=1120,680

def h(p):
    if p<=0 or p>=1:
        assert 0<=p<=1
        return 0.
    return -p*math.log2(p)-(1-p)*math.log2(1-p)

def entropy(p):
    return -sum(v*math.log2(v) for v in p if v>0)

def kl(p,q):
    return sum(a*math.log2(a/b) for a,b in zip(p,q) if a>0)

def fmt(x,n=6):
    return f'{x:.{n}f}'.rstrip('0').rstrip('.') if isinstance(x,float) else str(x)

def number(x):
    return f'{x:.10f}' if isinstance(x,float) else str(x)

def csv_content(cols,rows):
    out=io.StringIO(newline='')
    writer=csv.writer(out,lineterminator='\n');writer.writerow(cols)
    for row in rows: writer.writerow([number(x) for x in row])
    return out.getvalue().encode('utf-8')

class Figure:
    def __init__(self,n,title,subtitle,detail):
        self.n=n;self.a=[]
        self.add(f'<svg xmlns="http://www.w3.org/2000/svg" width="{W}" height="{H}" viewBox="0 0 {W} {H}" role="img" aria-labelledby="title desc">')
        self.add(f'<title id="title">{escape(title)}</title><desc id="desc">{escape(detail)}</desc>')
        self.rect(0,0,W,H,BG);self.rect(0,0,W,9,TEAL)
        self.txt(65,47,f'COMPLEXA · FIGURA {n}.1',14,TEAL,True)
        self.txt(65,87,title,26,INK,True)
        self.txt(65,119,subtitle,14,MUTED)
        self.line(65,139,W-65,139,GRID,1.5)
    def add(self,s): self.a.append(s)
    def rect(self,x,y,w,h,color,r=0,stroke=None):
        b=f' stroke="{stroke}"' if stroke else ''
        self.add(f'<rect x="{x:.2f}" y="{y:.2f}" width="{w:.2f}" height="{h:.2f}" rx="{r:.1f}" fill="{color}"{b}/>')
    def line(self,x,y,x2,y2,color=GRID,width=1,dash=None):
        dd=f' stroke-dasharray="{dash}"' if dash else ''
        self.add(f'<line x1="{x:.2f}" y1="{y:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" stroke="{color}" stroke-width="{width:.2f}"{dd}/>')
    def txt(self,x,y,s,size=14,color=INK,bold=False,anchor='start'):
        self.add(f'<text x="{x:.2f}" y="{y:.2f}" text-anchor="{anchor}" fill="{color}" font-family="Arial,Helvetica,sans-serif" font-size="{size:.1f}" font-weight="{700 if bold else 400}">{escape(str(s))}</text>')
    def circle(self,x,y,r,fill):self.add(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="{r:.2f}" fill="{fill}" stroke="#fff" stroke-width="1.6"/>')
    def poly(self,xy,col,w=3,dash=None):
        pts=' '.join(f'{x:.2f},{y:.2f}' for x,y in xy)
        d=f' stroke-dasharray="{dash}"' if dash else ''
        self.add(f'<polyline points="{pts}" fill="none" stroke="{col}" stroke-width="{w:.2f}" stroke-linejoin="round" stroke-linecap="round"{d}/>')
    def legend(self,parts,x=735,y=166):
        for i,(name,color) in enumerate(parts):
            yy=y+i*23
            self.rect(x,yy-12,24,9,color,2);self.txt(x+34,yy-3,name,13,MUTED)
    def note(self,lead,txt):
        self.rect(65,583,W-130,68,LT,8,GRID)
        self.txt(82,610,lead,14,TEAL,True)
        self.txt(82,635,txt,13,MUTED)
    def output(self):return ('\n'.join(self.a+['</svg>'])+'\n').encode('utf-8')

class Axes:
    def __init__(self,f,x0=119,y0=188,w=847,h=322,xmin=0,xmax=1,ymin=0,ymax=1,xlabel='',ylabel=''):
        assert xmax>xmin and ymax>ymin
        self.f=f;self.x0=x0;self.y0=y0;self.w=w;self.h=h
        self.xmin=xmin;self.xmax=xmax;self.ymin=ymin;self.ymax=ymax
        f.rect(x0,y0,w,h,BG,4,GRID)
        f.line(x0,y0+h,x0+w,y0+h,INK,1.8);f.line(x0,y0,x0,y0+h,INK,1.8)
        f.txt(x0+w/2,y0+h+53,xlabel,15,INK,anchor='middle')
        f.txt(x0-13,y0-14,ylabel,13,MUTED)
    def xy(self,x,y):return (self.x0+(x-self.xmin)/(self.xmax-self.xmin)*self.w,
                                  self.y0+self.h-(y-self.ymin)/(self.ymax-self.ymin)*self.h)
    def ticks(self,xs,ys):
        for x,s in xs:
            px,_=self.xy(x,self.ymin);self.f.line(px,self.y0,px,self.y0+self.h,GRID,1)
            self.f.txt(px,self.y0+self.h+24,s,12,MUTED,anchor='middle')
        for y,s in ys:
            _,py=self.xy(self.xmin,y);self.f.line(self.x0,py,self.x0+self.w,py,GRID,1)
            self.f.txt(self.x0-12,py+4,s,12,MUTED,anchor='end')
    def curve(self,pts,col,width=3):self.f.poly([self.xy(*xy) for xy in pts],col,width)
    def dot(self,x,y,color=ORANGE,r=5):self.f.circle(*self.xy(x,y),r,color)
    def dotted(self,x,y,col=ORANGE):
        p1=self.xy(x,self.ymin);p2=self.xy(x,y);self.f.line(*p1,*p2,col,1.5,'4 6')



# The examples below reproduce manuscript equations and parameters; no measured data.
# Their mathematical checks are independent of the drawing routines.


def bar(f, axes, x, y, width, color):
    x1, y0 = axes.xy(x-width/2, 0)
    x2, y1 = axes.xy(x+width/2, y)
    if y > 0: f.rect(x1, y1, x2-x1, y0-y1, color, 3)


def backlog(arrivals,capacity,initial=0):
    q=initial; out=[]
    for a in arrivals:
        available=q+a
        served=min(available,capacity)
        q=available-served
        assert q>=0 and 0<=served<=capacity and q+served==available
        out.append(q)
    return out


def build231():
    plain=[8]*10
    burst=[16 if t%2==0 else 0 for t in range(10)]
    values={('Constante',10):backlog(plain,10),('Constante',12):backlog(plain,12),
            ('Alternante',10):backlog(burst,10),('Alternante',12):backlog(burst,12)}
    assert values['Alternante',10]==[6,0]*5
    assert values['Alternante',12]==[4,0]*5
    assert all(v==0 for v in values['Constante',10]+values['Constante',12])
    assert sum(plain)==sum(burst)==80
    f=Figure(231,'La media de llegadas no predice la cola','Diez días · ochenta pedidos · dos patrones · capacidades de 10 y 12 pedidos/día',
      'Saldo al cierre diario con llegadas constantes de ocho pedidos o alternantes de dieciséis y cero. La cola es cero en patrón constante y alterna entre seis y cero con capacidad diez o cuatro y cero con capacidad doce.')
    a=Axes(f,xmin=.8,xmax=10.2,ymin=0,ymax=7.1,xlabel='Día del experimento',ylabel='Pedidos pendientes al cierre')
    a.ticks([(i,str(i)) for i in range(1,11)],[(v,str(v)) for v in (0,2,4,6)])
    for cap,color in [(10,ORANGE),(12,TEAL)]:
        points=list(zip(range(1,11),values['Alternante',cap]));a.curve(points,color,3.5)
        for x,y in points:a.dot(x,y,color,4.7)
    a.curve([(i,0) for i in range(1,11)],MUTED,2)
    f.legend([('Alternante, C=10',ORANGE),('Alternante, C=12',TEAL),('Constante: cola cero',MUTED)],775,165)
    f.note('Misma demanda total; distinta congestión','Alternante: promedio 3 / máximo 6 (C=10); promedio 2 / máximo 4 (C=12).')
    rows=[(i+1,plain[i],burst[i],values['Constante',10][i],values['Alternante',10][i],values['Constante',12][i],values['Alternante',12][i]) for i in range(10)]
    return f.output(),csv_content(['dia','llegada_constante','llegada_alternante','saldo_constante_C10','saldo_alternante_C10','saldo_constante_C12','saldo_alternante_C12'],rows)


def valid_no_adjacent(n):
    a,b=1,2
    if n==0:return a
    if n==1:return b
    for _ in range(2,n+1):a,b=b,a+b
    return b


def build232():
    rows=[]
    for n in range(1,31):
        valid=valid_no_adjacent(n)
        p10=1-2**(-n)
        p12=1-valid/2**n
        assert 0<=p12<=p10<=1
        rows.append((n,p10,p12,valid,2**n))
    assert valid_no_adjacent(30)==2178309
    assert math.isclose(rows[-1][2],.997971,abs_tol=6e-7)
    f=Figure(232,'Un umbral extremo exige otra verificación','Probabilidad exacta de superar cuatro pedidos pendientes; lotes de 16 con probabilidad 1/2',
      'Probabilidad de excedencia según número de días para capacidades diez y doce, derivada de ausencia total de lotes o de evitar lotes consecutivos. Valores exactos del modelo, no frecuencias simuladas.')
    a=Axes(f,xmin=0,xmax=30,ymin=0,ymax=1.06,xlabel='Horizonte de simulación (días)',ylabel='P(máximo de cola > 4)')
    a.ticks([(v,str(v)) for v in (0,5,10,15,20,25,30)],[(v,fmt(v,1)) for v in (0,.2,.4,.6,.8,1)])
    a.curve([(0,0)]+[(n,p10) for n,p10,*_ in rows],ORANGE,3.8)
    a.curve([(0,0)]+[(n,p12) for n,_,p12,*_ in rows],TEAL,3.8)
    f.legend([('Capacidad 10: algún lote',ORANGE),('Capacidad 12: lotes contiguos',TEAL)],733,167)
    f.note('Comprobación combinatoria exacta','En 30 días: 2.178.309 secuencias sin dos lotes consecutivos entre 2³⁰.')
    return f.output(),csv_content(['dias','p_exceso_C10','p_exceso_C12','secuencias_sin_consecutivos','secuencias_totales'],rows)


def build233():
    rows=[]
    N=10000
    for i in range(100):
        rho=i/100
        tau=(1+rho*rho)/(1-rho*rho)
        ess=N/tau
        assert tau>=1 and 0<ess<=N
        rows.append((rho,tau,ess))
    assert math.isclose(rows[-1][1],(1+.99**2)/(1-.99**2),rel_tol=1e-12)
    assert 99<rows[-1][2]<102
    f=Figure(233,'Gibbs: aceptación perfecta, muestras redundantes','Objetivo normal bivariado · barrido sistemático · 10.000 barridos completos',
      'Número efectivo asintótico de barridos para estimar la media de una coordenada frente a correlación rho. A rho 0,99, diez mil barridos proporcionan una precisión equivalente a cerca de cien observaciones independientes.')
    a=Axes(f,xmin=0,xmax=1,ymin=0,ymax=4.15,xlabel='Correlación ρ del objetivo',ylabel='log₁₀ del tamaño efectivo Nef')
    a.ticks([(v,fmt(v,2)) for v in (0,.2,.4,.6,.8,1)],[(i,'10'+('⁰' if i==0 else '¹' if i==1 else '²' if i==2 else '³' if i==3 else '⁴')) for i in range(5)])
    a.curve([(rho,math.log10(ess)) for rho,_,ess in rows],TEAL,3.8)
    a.dot(.99,math.log10(rows[-1][2]),ORANGE,7)
    f.note('La métrica depende del estimando','ρ=0,99: τ≈99,5; Nef≈100,5 para la media de una coordenada.')
    return f.output(),csv_content(['rho','factor_autocorrelacion_integrada','N_efectivo_10000_barridos'],rows)


def build234():
    N=10000
    rows=[]
    for i in range(141):
        phi=-.5+i/100
        tau=(1+phi)/(1-phi)
        ess=N/tau
        rows.append((phi,tau,ess,1/math.sqrt(ess)))
    assert math.isclose((1+.9)/(1-.9),19)
    assert math.isclose(rows[140][2],N/19)
    assert math.isclose(rows[0][2],30000)
    f=Figure(234,'Dependencia y precisión de una media MCMC','Proceso AR(1) estacionario · varianza unitaria · 10.000 observaciones',
      'Tamaño efectivo de la media frente a coeficiente AR1 desde menos 0,5 hasta 0,9. Con correlación negativa, el tamaño efectivo supera diez mil, pero con phi 0,9 desciende a 526.')
    a=Axes(f,xmin=-.5,xmax=.9,ymin=0,ymax=31500,xlabel='Coeficiente autorregresivo φ',ylabel='N efectivo de la media (observaciones)')
    a.ticks([(v,fmt(v,2)) for v in (-.5,-.25,0,.25,.5,.75,.9)],[(v,f'{v:,}'.replace(',','.')) for v in (0,5000,10000,20000,30000)])
    a.curve([(phi,ess) for phi,_,ess,_ in rows],BLUE,3.5)
    a.curve([(-.5,10000),(.9,10000)],MUTED,2)
    a.dot(.9,N/19,ORANGE,6)
    f.legend([('AR(1): N efectivo',BLUE),('Independencia: 10.000',MUTED)],745,165)
    f.note('Precisión, no número literal de observaciones','φ=0,9: τ=19, Nef≈526,32; φ=−0,5: Nef=30.000 para esta media.')
    return f.output(),csv_content(['phi','tau_integrada','Nef_media_N10000','error_estandar_media_var1'],rows)


def normal_tail(z):
    return .5*math.erfc(z/math.sqrt(2))


def build235():
    a0=5;mu=5
    p=normal_tail(a0)
    second=math.exp(mu**2)*normal_tail(a0+mu)
    vi=second-p*p
    vd=p*(1-p)
    assert p>0 and vi>0 and vi<vd
    assert math.isclose(p,2.866515718791933e-7,rel_tol=1e-9)
    assert math.isclose(second,5.487e-13,rel_tol=.001)
    rows=[]
    for exponent in [2+i*.1 for i in range(71)]:
        n=10**exponent
        rd=math.sqrt(vd/n)/p
        ri=math.sqrt(vi/n)/p
        rows.append((n,rd,ri,p,vi))
    ri_10k=math.sqrt(vi/10000)/p
    assert .023<ri_10k<.025
    f=Figure(235,'Una cola rara cambia el presupuesto de simulación','Normal estándar: P(X>5) ≈ 2,867 × 10⁻⁷; importancia con propuesta N(5,1)',
      'Error estándar relativo teórico, en escala logarítmica, para Monte Carlo directo e importancia normal desplazada, frente al número de simulaciones. La estimación de importancia reduce drásticamente la varianza de este evento.')
    a=Axes(f,xmin=2,xmax=9,ymin=-4,ymax=2.8,xlabel='Número de sorteos N (escala log₁₀)',ylabel='Error estándar relativo (escala log₁₀)')
    a.ticks([(k,f'10^{k}') for k in range(2,10)],[(k,'10'+('⁻⁴' if k==-4 else '⁻³' if k==-3 else '⁻²' if k==-2 else '⁻¹' if k==-1 else '⁰' if k==0 else '¹' if k==1 else '²')) for k in range(-4,3)])
    a.curve([(math.log10(n),math.log10(rd)) for n,rd,*_ in rows],ORANGE,3.8)
    a.curve([(math.log10(n),math.log10(ri)) for n,_,ri,*_ in rows],TEAL,3.8)
    f.legend([('Directo',ORANGE),('Importancia: q=N(5,1)',TEAL)],737,163)
    f.note('Cálculo de segundos momentos, no una simulación','Con N=10.000, error relativo teórico ≈2,38 % usando importancia.')
    return f.output(),csv_content(['N_sorteos','error_relativo_directo','error_relativo_importancia','p_objetivo','varianza_importancia_una_observacion'],rows)


def build236():
    rows=[('Directo (1 evaluación)',1,4/45,1),('Directo (2 evaluaciones)',2,2/45,2),
          ('Antitético (2 evaluaciones)',2,1/180,3),('Control conocido (1 evaluación)',1,1/180,4)]
    assert math.isclose((2/45)/(1/180),8)
    assert math.isclose((4/45)/(1/180),16)
    f=Figure(236,'Dos técnicas, distinto costo de evaluación','Integral ∫₀¹x²dx = 1/3 · varianzas analíticas de estimadores insesgados',
      'Cuatro barras comparan la varianza del estimador directo con una o dos evaluaciones, promedio antitético con dos y variable de control con una. Antitético reduce ocho veces la varianza a igual número de evaluaciones; control, dieciséis veces frente a una.')
    a=Axes(f,xmin=.45,xmax=4.55,ymin=0,ymax=.10,xlabel='Diseño (enumeración 1–4, ver tabla y nota)',ylabel='Varianza del estimador')
    a.ticks([(i,str(i)) for i in range(1,5)],[(v,fmt(v,3)) for v in (0,.02,.04,.06,.08,.10)])
    for _,_,v,i in rows:
        bar(f,a,i,v,.53,COLORS[i-1]);xp,yp=a.xy(i,v);f.txt(xp,yp-12,fmt(v,5),14,INK,True,anchor='middle')
    f.note('Comparar siempre con el mismo presupuesto','1 directo→control: varianza /16; 2 directos→antitético: /8.')
    return f.output(),csv_content(['metodo','evaluaciones','varianza_estadistico','indice'],rows)


def build237():
    rows=[]
    for b in [-1,1]:
        for a0 in [-1,1]:
            val=10+2*a0-b+3*a0*b
            rows.append((a0,b,val,a0*b))
    got={(a0,b):y for a0,b,y,_ in rows}
    assert [got[-1,-1],got[1,-1],got[-1,1],got[1,1]]==[12,10,4,14]
    main_A=(got[1,1]+got[1,-1]-got[-1,1]-got[-1,-1])/2
    main_B=(got[1,1]+got[-1,1]-got[1,-1]-got[-1,-1])/2
    intAB=(got[1,1]+got[-1,-1]-got[1,-1]-got[-1,1])/2
    assert (main_A,main_B,intAB)==(4,-2,6)
    f=Figure(237,'Un efecto puede invertir su dirección','Factorial completo 2×2 · Y(A,B)=10+2A−B+3AB · niveles codificados −1/+1',
      'Diagrama de interacción de un factorial con dos niveles. Cuando B es bajo, aumentar A reduce Y de 12 a 10; cuando B es alto, aumenta Y de 4 a 14. Efectos principales A=4 B=-2 e interacción=6.')
    a=Axes(f,xmin=-1.25,xmax=1.25,ymin=0,ymax=16,xlabel='Factor A (nivel codificado)',ylabel='Respuesta Y (unidades de modelo)')
    a.ticks([(-1,'−1'),(1,'+1')],[(v,str(v)) for v in (0,4,8,12,16)])
    for b,color in [(-1,BLUE),(1,ORANGE)]:
        pts=[(i,got[i,b]) for i in [-1,1]];a.curve(pts,color,4)
        for x,y in pts:a.dot(x,y,color,7);xx,yy=a.xy(x,y);f.txt(xx+(19 if x<0 else -19),yy-12,str(y),15,color,True,'start' if x<0 else 'end')
    f.legend([('B=−1',BLUE),('B=+1',ORANGE)],765,166)
    f.note('Interacción matemática explícita','Efecto simple de A: −2 con B bajo y +10 con B alto. Contraste AB = 6.')
    return f.output(),csv_content(['A_codificado','B_codificado','respuesta_Y','producto_AB'],rows)


def build238():
    components=[('X',1/3,1/8),('2Z',4/3,1/2),('3XZ',1,3/8)]
    assert math.isclose(sum(r[1] for r in components),8/3)
    assert math.isclose(sum(r[2] for r in components),1)
    sx=1/8;sz=1/2;stx=1/2;stz=7/8
    rows=[('X',1/3,sx,stx),('Z',4/3,sz,stz),('Interacción XZ',1,3/8,3/8)]
    f=Figure(238,'Sobol: una entrada puede importar por interacción','X y Z independientes U(−1,1) · Y=X+2Z+3XZ · varianza total 8/3',
      'Tres barras muestran contribución de varianza de X, Z e interacción XZ, que suman uno. Dos pares de barras comparan índices Sobol de primer orden y totales, donde X pasa de 0,125 a 0,5 al contar interacción.')
    a=Axes(f,xmin=.5,xmax=5.5,ymin=0,ymax=1.06,xlabel='1: X, 2: Z, 3: interacción, 4: ST X, 5: ST Z',ylabel='Fracción de varianza (adimensional)')
    a.ticks([(i,str(i)) for i in range(1,6)],[(v,fmt(v,3)) for v in (0,.25,.5,.75,1)])
    for i,(_,_,fraction) in enumerate(components,1):bar(f,a,i,fraction,.49,COLORS[i-1])
    bar(f,a,4,stx,.49,BLUE);bar(f,a,5,stz,.49,TEAL)
    f.note('Índices de primer orden ≠ índices totales','S_X=0,125 y S_Z=0,5; ST_X=0,5 y ST_Z=0,875; interacción=0,375.')
    return f.output(),csv_content(['componente','varianza_absoluta','fraccion_varianza_o_Sprimero','indice_total_si_aplicable'],rows)


def build239():
    def f0(x):return x*x
    def g(x):return x*x+8*x*(x-.5)*(x-1)
    def interp(x):
        if x<=.5:return .5*x
        return 1.5*x-.5
    rows=[]
    for j in range(101):
        x=j/100;rows.append((x,f0(x),g(x),interp(x),interp(x)-f0(x)))
    for x in (0,.5,1):assert math.isclose(f0(x),g(x),abs_tol=1e-14)
    assert math.isclose(g(.25),.4375) and math.isclose(g(.75),.1875)
    assert math.isclose(max(z[4] for z in rows),1/16,abs_tol=1e-12)
    assert math.isclose(sum((1/3,1/24)),3/8)
    f=Figure(239,'Tres nodos no identifican un simulador','f(x)=x², g(x)=x²+8x(x−½)(x−1); interpolación lineal en 0, ½ y 1',
      'Tres curvas sobre el intervalo unidad coinciden en cero, un medio y uno, pero difieren en puntos no entrenados. Una función cúbica rival ajusta exactamente los tres datos; la interpolación lineal tiene máximo error 0,0625 respecto de x cuadrado.')
    a=Axes(f,xmin=0,xmax=1,ymin=0,ymax=1.07,xlabel='Entrada x (sin unidades)',ylabel='Respuesta del modelo (adimensional)')
    a.ticks([(v,fmt(v,2)) for v in (0,.25,.5,.75,1)],[(v,fmt(v,1)) for v in (0,.2,.4,.6,.8,1)])
    for index,col in [(1,BLUE),(2,ORANGE),(3,TEAL)]:a.curve([(row[0],row[index]) for row in rows],col,3.2)
    for x in (0,.5,1):a.dot(x,f0(x),INK,5)
    f.legend([('f: x²',BLUE),('g: cúbica rival',ORANGE),('Interpolación lineal',TEAL)],756,161)
    f.note('La validación necesita entradas nuevas','En x=0,25: f=0,0625, g=0,4375; error máximo lineal vs. f=1/16.')
    return f.output(),csv_content(['x','f_cuadratica','g_cubica','interpolacion_lineal','error_interpolacion_menos_f'],rows)


def exact_queue_mean(capacity,block):
    # Deterministic count aggregation: (count, sum over all prior areas) per queue state.
    states={0:(1,0)}
    steps=15 if block else 30
    for _ in range(steps):
        new={}
        for Q,(count,area) in states.items():
            for incoming in (0,16):
                q=max(0,Q+incoming-capacity)
                increment=q
                if block:
                    q2=max(0,q+incoming-capacity)
                    increment+=q2
                    q=q2
                oldcount,oldarea=new.get(q,(0,0))
                new[q]=(oldcount+count,oldarea+area+count*increment)
        states=new
    ways=2**steps
    assert sum(c for c,_ in states.values())==ways
    total_area=sum(a for _,a in states.values())
    from fractions import Fraction
    return Fraction(total_area,ways*30)


def build240():
    i10=exact_queue_mean(10,False);i12=exact_queue_mean(12,False)
    b10=exact_queue_mean(10,True);b12=exact_queue_mean(12,True)
    # Chapter C240.6 publishes rounded decimals, not measured data.
    assert abs(float(i10)-9.495447)<.000002
    assert abs(float(i12)-4.412897)<.000002
    assert abs(float(b10)-15.054012)<.000002
    assert abs(float(b12)-7.641032)<.000002
    diff_i=float(i12-i10);diff_b=float(b12-b10)
    assert math.isclose(diff_i+6,.917450,abs_tol=2e-6)
    assert math.isclose(diff_b+6,-1.412980,abs_tol=2e-6)
    threshold_i=-diff_i/2;threshold_b=-diff_b/2
    assert threshold_i<3<threshold_b
    rows=[(j*.05,diff_i+2*j*.05,diff_b+2*j*.05,threshold_i,threshold_b) for j in range(101)]
    f=Figure(240,'La preferencia cambia cuando cambia la dependencia','Dos capacidades: 10 y 12 · horizonte 30 días · diferencia exacta de costos frente a λ',
      'Dos líneas de diferencia esperada del costo capacidad doce menos diez frente a penalización de capacidad lambda. A lambda tres, la independencia favorece capacidad diez y los lotes emparejados favorecen capacidad doce.')
    a=Axes(f,xmin=0,xmax=5,ymin=-8,ymax=6,xlabel='Costo unitario adicional λ (por capacidad y día)',ylabel='J(12) − J(10) (unidades de costo)')
    a.ticks([(v,str(v)) for v in range(6)],[(v,str(v)) for v in (-8,-6,-4,-2,0,2,4,6)])
    a.curve([(0,0),(5,0)],MUTED,2.2)
    a.curve([(l,d) for l,d,*_ in rows],ORANGE,3.5)
    a.curve([(l,d) for l,_,d,*_ in rows],TEAL,3.5)
    a.dot(3,diff_i+6,ORANGE,7);a.dot(3,diff_b+6,TEAL,7)
    f.legend([('Llegadas independientes',ORANGE),('Lotes en pares',TEAL),('Empate de costos',MUTED)],748,165)
    f.note('Referencia por conteo exacto de trayectorias','λ=3: ΔJ≈+0,91745 (elige C=10) y −1,41298 (elige C=12).')
    return f.output(),csv_content(['lambda','J12_menos_J10_independiente','J12_menos_J10_pares','lambda_critico_independiente','lambda_critico_pares'],rows)


CASES={n:globals()[f'build{n}'] for n in range(231,241)}

def generated():
    result={}
    for n,fn in CASES.items():
        svg,data=fn()
        result[f'figures/complexa-c{n}-fig01.svg']=svg
        result[f'repro/c231-240/c{n}.csv']=data
    return result


def bundle(files):
    buf=io.BytesIO()
    with ZipFile(buf,'w',compression=ZIP_STORED) as z:
        for path,content in sorted(files.items()):
            name=path.removeprefix('figures/').removeprefix('repro/c231-240/')
            inf=ZipInfo(name,(2026,1,1,0,0,0));inf.compress_type=ZIP_STORED
            inf.create_system=3;inf.external_attr=0o644<<16
            z.writestr(inf,content)
    return buf.getvalue()


def run(check:bool,root:Path):
    files=generated()
    import xml.etree.ElementTree as ET
    for n in CASES:
        svg=files[f'figures/complexa-c{n}-fig01.svg']
        element=ET.fromstring(svg)
        assert element.attrib.get('role')=='img'
        assert any(child.tag.endswith('title') for child in element)
        assert any(child.tag.endswith('desc') for child in element)
        data=list(csv.reader(io.StringIO(files[f'repro/c231-240/c{n}.csv'].decode())))
        assert len(data)>=4 and len(data[0])>=3,n
        assert all(len(x)==len(data[0]) for x in data[1:]),n
    extra={'repro/c231-240/make_figures.py':BASE.joinpath('make_figures.py').read_bytes(),
           'repro/c231-240/README.md':BASE.joinpath('README.md').read_bytes()}
    zip_bytes=bundle({**files,**extra})
    if check:
        for path,expected in files.items():
            target=root/path
            assert target.is_file() and target.read_bytes()==expected,f'CHECK FAIL: {target}'
        zp=root/'zip-data/complexa-c231-c240-reproducible.zip.b64'
        import base64
        assert zp.is_file() and base64.b64decode(zp.read_text().strip())==zip_bytes,'CHECK FAIL: zip'
        print('PASS C231–C240: 10 SVG, 10 CSV, ZIP bitwise reproducible; numerical checks.')
        return
    for path,data in files.items():
        target=root/path;target.parent.mkdir(parents=True,exist_ok=True);target.write_bytes(data)
    import base64
    zp=root/'zip-data/complexa-c231-c240-reproducible.zip.b64'
    zp.parent.mkdir(parents=True,exist_ok=True)
    zp.write_text(base64.b64encode(zip_bytes).decode()+'\n')
    direct=root/'complexa-c231-c240-reproducible.zip'
    direct.write_bytes(zip_bytes)
    print(f'PASS C231–C240: {len(files)} SVG+CSV; ZIP {len(zip_bytes)} bytes; SHA256 {hashlib.sha256(zip_bytes).hexdigest()}')

if __name__=='__main__':
    ap=argparse.ArgumentParser();g=ap.add_mutually_exclusive_group(required=True)
    g.add_argument('--out',type=Path);g.add_argument('--check',type=Path)
    args=ap.parse_args()
    run(args.check is not None,args.check or args.out)