#!/usr/bin/env python3
"""COMPLEXA C311-C320: 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.ordinal=3 if n==320 else 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)}. {"Comparación retrospectiva de datos reales históricos del caso C320." if n==320 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,"REGISTRO HISTÓRICO / RETROSPECTIVO" if n==320 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")






# Numerical companions reproduce the parameter choices printed in each chapter.
# The only real-world data here are the original C320 historical-model summary
# statistics, already published alongside a versioned original data archive.
NAMES = {
 311:("Dos frecuencias distintas, lecturas idénticas", "Aliasing: ondas de 4 y 6 ciclos/min al muestrear 10 veces/min"),
 312:("Una tendencia no es un paseo aleatorio", "Varianzas esperadas: ruido estacionario, raíz unitaria y AR(1)"),
 313:("ACF y PACF responden preguntas diferentes", "AR(1) teórico con phi = 0,6; sin observaciones simuladas"),
 314:("Una perturbación conserva efectos diferentes", "Pesos exactos de un AR(1) y de un ARMA(1,1)"),
 315:("Dos canales conservan memoria cruzada", "VAR(1) estable desde x = 2 e y = −1: evolución del pronóstico"),
 316:("Granger mide mejora de pronóstico", "Tres procesos explícitos: información adicional y horizonte temporal"),
 317:("La independencia no se ve en una sola ACF", "Residuos AR(1) frente a AR(2); memoria ARCH en cuadrados"),
 318:("Pronosticar solo con datos ya confirmados", "Ocho orígenes, horizonte de dos días y demora de registro"),
 319:("Un episodio breve necesita dos fronteras", "Doce observaciones, mínimo de tres y penalización 10 por corte"),
 320:("La frecuencia del registro cambia el pronóstico", "MAE del consumo anual 1990–2008; misma fecha objetivo"),
}

def matmul(A,x): return [sum(a*b for a,b in zip(row,x)) for row in A]


def make_data():
    D={}
    # C311. Aliasing: at t=n/10, cos(2π4t)=cos(2π6t).
    rows=[]
    for n in range(61):
        t=n/100
        a=10+2*cos(2*pi*4*t)
        b=10+2*cos(2*pi*6*t)
        sampled=int(n%10==0)
        rows.append((t,a,b,sampled))
        if sampled: assert isclose(a,b,abs_tol=1e-12)
    assert rows[0][1]==rows[0][2]==12
    D[311]=rows

    # C312. True variance of linear trend with iid additive noise is 1.
    # Random walk variance increases linearly. AR0.98 relaxes to 25.2525.
    rows=[]
    phi=.98
    for t in range(101):
        ar=(1-phi**(2*t))/(1-phi**2)
        rows.append((t,1.,float(t),ar))
    assert isclose(rows[50][3],(1-.98**100)/(1-.98**2))
    assert isclose(1/(1-.98**2),25.252525252525235)
    D[312]=rows

    # C313. Theoretical PACF of a stationary AR(1) is zero beyond lag 1.
    D[313]=[(h,.6**h,1. if h==0 else .6 if h==1 else 0.,.6**h-(.6 if h==1 else 0.)) for h in range(9)]
    assert isclose(D[313][2][1],.36) and D[313][2][2]==0

    # C314. ARMA with phi .6 and theta .4: ψ0=1, ψ1=1, ψh=.6^(h-1).
    rows=[]
    for h in range(13):
        arma=1 if h==0 else .6**(h-1)
        ar=.8**h
        ma=1 if h==0 else .4 if h==1 else 0
        rows.append((h,arma,ar,ma))
    assert isclose(sum((1 if h==0 else .6**(h-1))**2 for h in range(150)),2.5625,abs_tol=1e-9)
    D[314]=rows

    # C315. The one-step VAR forecast without new observations is A^h x0.
    A=((.5,.4),(.2,.3));state=[2.,-1.]
    rows=[]
    for h in range(16):
        rows.append((h,state[0],state[1],state[0]+state[1]))
        state=matmul(A,state)
    assert isclose(rows[1][1],.6) and isclose(rows[1][2],.1)
    assert isclose(rows[2][1],.34) and isclose(rows[2][2],.15)
    D[315]=rows

    # C316. The first experiment is a direct predictive link. The second
    # is a common cause, so these ratios are NOT direct causal effects.
    rows=[
      ('X anticipa Y',1,1.,1.64,log(1.64)),
      ('Causa común no observada',1,1.5,2.,log(2/1.5)),
      ('Vía mediada por Z',2,1.25,1.41,log(1.41/1.25))]
    assert isclose(rows[0][4],.494696241836107)
    D[316]=rows

    # C317. Sample residual ACF: source chapter's realized data estimates
    # are used only at lags 1 and 2. ARCH ACFs are population results.
    rows=[]
    for h in range(1,9):
        ar1={1:-.1743,2:.237}.get(h,"")
        ar2={1:.0058,2:.0197}.get(h,"")
        rows.append((h,ar1,ar2,0.,.5**h))
    assert rows[0][4]==.5
    D[317]=rows

    # C318. Each forecast uses only Y_{t-1} and Y_{t-4}.
    y=[20,22,21,24,23,25,27,26,29,28,30,33,31,34,35,34,36]
    rows=[]
    for t in range(8,16):
        actual=y[t+2-1] # Python index 0-based for day t+2
        fa=y[t-2]      # Y_{t-1}
        fb=2*y[t-2]-y[t-5] # Y_{t-4}
        rows.append((t,t+2,actual,fa,fb,actual-fa,actual-fb))
    assert [(r[3],r[4]) for r in rows]==[(27,30),(26,29),(29,33),(28,29),(30,34),(33,37),(31,34),(34,38)]
    assert isclose(sum(abs(r[5]) for r in rows)/8,2.875)
    assert isclose(sum(abs(r[6]) for r in rows)/8,1.125)
    D[318]=rows

    # C319. All admissible partitions with minimum length 3 enumerated.
    y=[0]*4+[3]*4+[0]*4
    def cost(a,b):
        xx=y[a:b];mu=sum(xx)/len(xx)
        return sum((z-mu)**2 for z in xx)
    solutions={}
    from itertools import combinations
    for m in range(4):
        candidates=[]
        for cuts in combinations(range(3,10),m):
            edges=(0,)+cuts+(12,)
            if min(edges[k+1]-edges[k] for k in range(len(edges)-1))<3:continue
            loss=sum(cost(edges[i],edges[i+1]) for i in range(len(edges)-1))
            candidates.append((loss+10*m,loss,cuts))
        assert candidates
        solutions[m]=min(candidates)
    assert [round(solutions[m][1]) for m in range(4)]==[24,18,0,12]
    assert [round(solutions[m][0]) for m in range(4)]==[24,28,20,42]
    assert solutions[2][2]==(4,8)
    D[319]=[(t+1,ys,0 if t<4 or t>=8 else 3,int(t in (3,7)),
             solutions[t][0] if t<4 else "",solutions[t][1] if t<4 else "")
            for t,ys in enumerate(y)]

    # C320. These are real-data *retrospective summary metrics* reproduced
    # by the already versioned C320/macrodata.csv and C320/reproduce.py.
    # No fabricated observation, no claim of real-time vintage data.
    D[320]=[
        ('Anual','Persistencia',19,217.22,237.34),
        ('Anual','Crecimiento',19,76.66,112.06),
        ('Trimestral','Persistencia',19,137.76,148.61),
        ('Trimestral','Crecimiento',19,57.13,82.77),
        ('Trimestral','AR(1)',19,53.93,78.68)
    ]
    assert D[320][-1][3]==53.93
    return D


def charts(n,data):
    title,desc=NAMES[n]
    s=SVG(n,title,desc)
    if n==311:
        plot(s,(52,176,545,449),0,.6,7.5,12.5,[0,.1,.2,.3,.4,.5,.6],[8,9,10,11,12],
             'Tiempo (min)','Magnitud',[('4 ciclos/min',TEAL,[(r[0],r[1]) for r in data]),('6 ciclos/min',OCHRE,[(r[0],r[2]) for r in data])],
             'DOS ONDAS DISTINTAS')
        plot(s,(613,176,533,449),0,.6,7.5,12.5,[0,.1,.2,.3,.4,.5,.6],[8,9,10,11,12],
             'Intervalos de 6 segundos','Lecturas', [('Ambas en el registro',PURPLE,[(r[0],r[1]) for r in data if r[3]])],
             'UNA SOLA SECUENCIA OBSERVADA',legend=False)
        # Sample dots mark actual discrete observations, not continuous paths.
        s.note('Los dos procesos coinciden en cada lectura, aunque entre ellas son diferentes.')
    elif n==312:
        plot(s,(55,175,1090,451),0,100,0,105,[0,20,40,60,80,100],[0,25,50,75,100],
             'Tiempo observado','Varianza esperada',[('Tendencia determinista + ruido',MUTED,[(r[0],r[1]) for r in data]),('Paseo aleatorio',OCHRE,[(r[0],r[2]) for r in data]),('AR(1) phi = 0,98',TEAL,[(r[0],r[3]) for r in data])],
             'VARIANZAS: MISMA PENDIENTE MEDIA, DISTINTA MEMORIA')
        s.note('La deriva lineal de la media no distingue un paseo de una tendencia determinista.')
    elif n==313:
        plot(s,(55,174,1090,455),0,8,0,1.05,[0,1,2,3,4,5,6,7,8],[0,.2,.4,.6,.8,1],
             'Retardo h','Coeficiente poblacional',[('ACF: rho(h) = 0,6^h',TEAL,[(r[0],r[1]) for r in data]),('PACF: solo lag 1',OCHRE,[(r[0],r[2]) for r in data])],
             'UNA ASOCIACIÓN LEJANA PUEDE ESTAR MEDIADA POR EL RETARDO 1')
        s.note('ACF(2) = 0,36, pero PACF(2) = 0 al controlar el antecedente inmediato.')
    elif n==314:
        plot(s,(55,174,1090,455),0,12,0,1.1,[0,2,4,6,8,10,12],[0,.2,.4,.6,.8,1],
             'Pasos tras la perturbación','Efecto relativo',[('ARMA(0,6; 0,4)',TEAL,[(r[0],r[1]) for r in data]),('AR(0,8)',OCHRE,[(r[0],r[2]) for r in data]),('MA(0,4)',PURPLE,[(r[0],r[3]) for r in data])],
             'PESOS DE RESPUESTA A UNA INNOVACIÓN UNITARIA')
        s.note('ARMA tiene peso 1 al paso siguiente y luego una cola geométrica con factor 0,6.')
    elif n==315:
        plot(s,(55,174,1090,456),0,15,-1.1,2.1,[0,3,6,9,12,15],[-1,0,1,2],
             'Horizonte (pasos)','Valor esperado',[('Canal x',TEAL,[(r[0],r[1]) for r in data]),('Canal y',OCHRE,[(r[0],r[2]) for r in data])],
             'PRONÓSTICO VAR(1) SIN CONOCER FUTURAS INNOVACIONES')
        s.note('Desde (2, −1): horizonte 1 = (0,6; 0,1); horizonte 2 = (0,34; 0,15).')
    elif n==316:
        # A grouped bar display: three experiments on exactly their own scale.
        s.rect(57,174,1087,452,'#FFFFFF',12,GRID)
        s.text(79,210,'VARIANZA DEL ERROR DE PRONÓSTICO',18,NAVY,700)
        s.line(140,574,1090,574,NAVY,1.8)
        for val in [0,.5,1,1.5,2]:
            cy=574-val*177
            s.line(135,cy,1090,cy,GRID,1)
            s.text(123,cy+4,f'{val:.1f}',12,MUTED,400,'end')
        for i,(lab,h,full,reduced,g) in enumerate(data):
            mid=320+i*320;left=mid-63;right=mid+18
            s.rect(left,574-full*177,61,full*177,TEAL,5)
            s.rect(right,574-reduced*177,61,reduced*177,OCHRE,5)
            s.text(left+30,564-full*177,f'{full:.2f}',15,TEAL,700,'middle')
            s.text(right+30,564-reduced*177,f'{reduced:.2f}',15,OCHRE,700,'middle')
            s.text(mid+8,610,['Enlace directo','Causa compartida','Ruta indirecta'][i],13,INK,600,'middle')
        s.line(790,207,815,207,TEAL,5);s.text(824,212,'Modelo completo',14,INK,500)
        s.line(969,207,993,207,OCHRE,5);s.text(1003,212,'Modelo reducido',14,INK,500)
        s.note('La mejora predictiva también aparece por una causa común; no identifica intervención.')
    elif n==317:
        plot(s,(55,175,530,453),1,2,-.25,.55,[1,2],[-.2,0,.2,.4],
             'Retardo','Correlación residual',[('AR(1)',OCHRE,[(r[0],r[1]) for r in data[:2]]),('AR(2)',TEAL,[(r[0],r[2]) for r in data[:2]])],
             'EFECTO DE UN RETARDO OMITIDO')
        plot(s,(605,175,540,453),1,8,0,.55,[1,2,3,4,5,6,7,8],[0,.1,.2,.3,.4,.5],
             'Retardo','ACF poblacional',[('Errores ARCH',PURPLE,[(r[0],r[3]) for r in data]),('Errores al cuadrado',TEAL,[(r[0],r[4]) for r in data])],
             'DEPENDENCIA EN LA ESCALA')
        s.note('ARCH: ACF de errores = 0; ACF de sus cuadrados a un paso = 0,5.')
    elif n==318:
        plot(s,(55,175,1090,452),10,17,24,39,[10,11,12,13,14,15,16,17],[24,27,30,33,36,39],
             'Día que se pronostica','Demanda',[('Valor observado',NAVY,[(r[1],r[2]) for r in data]),('Último dato confirmado',OCHRE,[(r[1],r[3]) for r in data]),('Extrapolación',TEAL,[(r[1],r[4]) for r in data])],
             'OCHO PRONÓSTICOS EMITIDOS SIN INFORMACIÓN FUTURA')
        s.note('Dos días adelante: MAE de la persistencia = 2,875; extrapolación = 1,125.')
    elif n==319:
        plot(s,(55,174,532,452),1,12,-.3,3.5,[1,3,5,7,9,11],[0,1,2,3],
             'Observación','Nivel',[('Registro 0−3−0',OCHRE,[(r[0],r[1]) for r in data]),('Tres tramos',TEAL,[(r[0],r[2]) for r in data])],
             'DOS FRONTERAS EN 4 Y 8')
        ss=[(i,[24,28,20,42][i]) for i in range(4)]
        plot(s,(610,174,535,452),0,3,0,45,[0,1,2,3],[0,10,20,30,40],
             'Número de cortes','Coste penalizado',[('Coste óptimo',PURPLE,ss)],
             'SELECCIÓN GLOBAL DEL CRITERIO',legend=False,highlight=[(2,20,'mínimo',TEAL)])
        s.note('El mejor corte aislado no compensa; los dos cortes juntos sí lo hacen.')
    elif n==320:
        bars=[('Anual/P',217.22,MUTED),('Anual/G',76.66,OCHRE),('Trim/P',137.76,PURPLE),('Trim/G',57.13,TEAL),('Trim/AR',53.93,GREEN)]
        bar_panel(s,(55,174,1090,455),'ERROR ABSOLUTO MEDIO · 19 AÑOS OBJETIVO',bars,260,'1990–2008 · MAE anual / miles de millones USD (2005)')
        s.note('Es el mismo promedio anual objetivo; las predicciones difieren en memoria y regla.')
    else: raise AssertionError(n)
    return s.svg()


def assert_near(got,expected,label,tol=1e-6):
    if abs(got-expected)>tol: raise AssertionError(f'{label}: {got} expected {expected}')


def audit_model(D):
    assert len(D)==10
    for n in range(311,321):assert len(D[n])>2 and len(D[n][0])>=4
    assert_near(D[311][10][1],D[311][10][2],'aliasing')
    assert_near(D[312][100][2],100,'random walk variance')
    assert_near(D[313][2][1],.36,'ACF lag2')
    assert_near(D[313][2][2],0,'PACF lag2')
    assert_near(D[314][1][1],1,'ARMA psi1')
    assert_near(D[315][2][1],.34,'VAR x2')
    assert_near(D[316][1][3]-D[316][1][2],.5,'common cause variance')
    assert_near(D[317][0][4],.5,'ARCH squares autocorrelation')
    assert_near(sum(abs(r[5]) for r in D[318])/8,2.875,'forecast MAE A')
    assert_near(sum(abs(r[6]) for r in D[318])/8,1.125,'forecast MAE B')
    assert_near(D[319][2][4],20,'best 2 cut objective')
    assert_near(D[320][-1][3],53.93,'historical annual MAE')
    return True


def csv_companions(D):
    headers={
      311:['time_min','signal_4_cycles_min','signal_6_cycles_min','is_sampled'],
      312:['time','trend_noise_variance','random_walk_variance','ar098_variance'],
      313:['lag','acf_AR06','pacf_AR06','acf_minus_pacf'],
      314:['lag','arma06_04_weight','ar08_weight','ma04_weight'],
      315:['horizon','forecast_x','forecast_y','forecast_sum'],
      316:['experiment','horizon','full_error_variance','reduced_error_variance','geweke_ln_ratio'],
      317:['lag','observed_residual_AR1','observed_residual_AR2','population_ARCH_ACF','population_ARCH_squared_ACF'],
      318:['issue_day','target_day','observed','forecast_persistence','forecast_extrapolation','error_persistence','error_extrapolation'],
      319:['time','value','segment_mean','is_boundary_after','optimum_penalized_first4','optimum_residual_first4'],
      320:['source_calendar','forecast_rule','annual_targets','MAE','RMSE']
    }
    return {f'repro/c311-320/c{n}.csv':csv_bytes(headers[n],D[n]) for n in sorted(D)}


def byte_outputs():
    D=make_data();audit_model(D)
    d={f'figures/complexa-c{n}-fig{3 if n==320 else 1:02d}.svg':charts(n,D[n]) for n in sorted(D)}
    d.update(csv_companions(D))
    for n in D:
        xml=d[f'figures/complexa-c{n}-fig{3 if n==320 else 1:02d}.svg']
        root=ElementTree.fromstring(xml)
        assert root.tag.endswith('svg') and root.find('{http://www.w3.org/2000/svg}title') is not None
        assert root.find('{http://www.w3.org/2000/svg}desc') is not None
    return d


def make_archive(files,readme,script):
    entries=dict(files)
    entries['repro/c311-320/README.md']=readme
    entries['repro/c311-320/make_figures.py']=script
    stream=BytesIO()
    with zipfile.ZipFile(stream,'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 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-c311-c320-reproducible.zip.b64']=base64.b64encode(archive)+b'\n'
    return d,archive


def run():
    ap=argparse.ArgumentParser(description=__doc__)
    group=ap.add_mutually_exclusive_group()
    group.add_argument('--out',type=Path,help='Output in content/supplementary')
    group.add_argument('--check',type=Path,help='Compare exact reproducible bytes')
    opt=ap.parse_args()
    dest=opt.out or opt.check or ROOT
    files,archive=outputs()
    for rel,blob in sorted(files.items()):
        target=dest/rel
        if opt.check:
            if not target.is_file() or target.read_bytes()!=blob:
                raise SystemExit(f'FAIL source 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
    msg='CHECK' if opt.check else 'GENERATED'
    print(f'{msg}: ten SVG, ten CSV, 22 reproducible ZIP entries; SHA-256 {sha256(archive).hexdigest()}')

if __name__=='__main__':run()
