#!/usr/bin/env python3
"""COMPLEXA C191–C200. Figures derived from published analytical examples.

Python 3.12 standard library only.  --generate updates local supplement assets,
--check regenerates all numerical SVG/CSV bytes and the canonical ZIP in memory.
No observation, inference, random sample or undocumented model parameter is used.
"""
import argparse
import base64
import csv
import hashlib
import html
import io
import json
import math
import zipfile
from pathlib import Path

ROOT = Path(__file__).resolve().parents[2]
TEAL = '#087e83'
BLUE = '#285e91'
AMBER = '#b66b28'
INK = '#1a3140'
MUTED = '#506579'
GRID = '#e2eaf0'
BG = '#f3f8fa'


def esc(s):
    return html.escape(str(s), quote=True)


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


def csv_bytes(cols, rows):
    f = io.StringIO(newline='')
    writer = csv.writer(f, lineterminator='\n')
    writer.writerow(cols)
    writer.writerows([[num(x) for x in row] for row in rows])
    return f.getvalue().encode('utf-8')


class Figure:
    def __init__(self, n, title, subtitle, description):
        self.parts = [
            '<svg xmlns="http://www.w3.org/2000/svg" width="1080" height="670" viewBox="0 0 1080 670" role="img" aria-labelledby="title desc">',
            f'<title id="title">{esc(title)}</title><desc id="desc">{esc(description)}</desc>',
            '<rect width="1080" height="670" fill="white"/>',
            f'<rect width="1080" height="8" fill="{TEAL}"/>',
        ]
        self.text(72, 51, f'COMPLEXA · FIGURA {n}.1', 15, TEAL, bold=True)
        self.text(72, 94, title, 27, INK, bold=True)
        self.text(72, 125, subtitle, 15, MUTED)
        self.line(72, 145, 1008, 145, GRID)

    def add(self, s): self.parts.append(s)

    def line(self, x1, y1, x2, y2, color=INK, width=1.5, dash=''):
        d = f' stroke-dasharray="{dash}"' if dash else ''
        self.add(f'<line x1="{x1:.2f}" y1="{y1:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" stroke="{color}" stroke-width="{width}"{d}/>')

    def text(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}" font-family="Arial,Helvetica,sans-serif" font-size="{size}" font-weight="{700 if bold else 400}" fill="{color}">{esc(s)}</text>')

    def rect(self, x, y, w, h, fill, radius=0):
        self.add(f'<rect x="{x:.2f}" y="{y:.2f}" width="{w:.2f}" height="{h:.2f}" rx="{radius}" fill="{fill}"/>')

    def circle(self, x, y, r=5, color=TEAL, hollow=False):
        self.add(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="{r}" fill="{ "white" if hollow else color}" stroke="{color}" stroke-width="2"/>')

    def poly(self, pts, color, width=3.5):
        coords = ' '.join(f'{x:.2f},{y:.2f}' for x, y in pts)
        self.add(f'<polyline points="{coords}" fill="none" stroke="{color}" stroke-width="{width}" stroke-linejoin="round" stroke-linecap="round"/>')

    def legend(self, entries):
        x = 125
        for label, color in entries:
            self.line(x, 572, x + 32, 572, color, 5)
            self.text(x + 41, 577, label, 15, INK)
            x += max(185, 53 + len(label) * 8.7)

    def note(self, caption):
        self.rect(72, 602, 936, 46, BG, 8)
        self.text(87, 630, caption, 13, MUTED)

    def bytes(self):
        return ('\n'.join(self.parts + ['</svg>']) + '\n').encode('utf-8')


class Plot:
    def __init__(self, f, xmin, xmax, ymin, ymax, xlabel, ylabel, xticks, yticks, logy=False):
        self.f = f
        self.x0, self.y0, self.w, self.h = 126, 188, 837, 327
        self.xmin, self.xmax, self.ymin, self.ymax = xmin, xmax, ymin, ymax
        self.logy = logy
        assert xmin < xmax and ymin < ymax and (not logy or ymin > 0)
        f.rect(self.x0, self.y0, self.w, self.h, '#ffffff')
        for x, label in xticks:
            px, _ = self.xy(x, ymin)
            f.line(px, self.y0, px, self.y0 + self.h, GRID, 1)
            f.text(px, self.y0 + self.h + 23, label, 12, MUTED, anchor='middle')
        for y, label in yticks:
            _, py = self.xy(xmin, y)
            f.line(self.x0, py, self.x0 + self.w, py, GRID, 1)
            f.text(self.x0 - 11, py + 4, label, 12, MUTED, anchor='end')
        f.line(self.x0, self.y0 + self.h, self.x0 + self.w, self.y0 + self.h, INK)
        f.line(self.x0, self.y0, self.x0, self.y0 + self.h, INK)
        f.text(self.x0 + self.w / 2, 557, xlabel, 15, INK, anchor='middle')
        f.text(self.x0, 176, ylabel, 14, MUTED)

    def xy(self, x, y):
        y0, y1 = (math.log10(self.ymin), math.log10(self.ymax)) if self.logy else (self.ymin, self.ymax)
        yy = math.log10(y) if self.logy else y
        return (self.x0 + self.w * (x-self.xmin)/(self.xmax-self.xmin),
                self.y0 + self.h * (y1-yy)/(y1-y0))

    def curve(self, xy, color, width=3.5):
        self.f.poly([self.xy(x, y) for x,y in xy], color, width)

    def dots(self, xy, color, size=5):
        for x, y in xy: self.f.circle(*self.xy(x, y), size, color)

    def hline(self, y, color=MUTED, dash='5,7'):
        self.f.line(*self.xy(self.xmin, y), *self.xy(self.xmax, y), color, 1.7, dash)


def build191():
    f = Figure(191, 'Tres modos y una relajación colectiva',
        'Cadena de tres nodos · Laplaciano con autovalores 0, 1 y 3',
        'Curvas de los tres estados, suma conservada y convergencia a un tercio.')
    p = Plot(f, 0, 5, 0, 1, 'Tiempo t (unidades del modelo)', 'Estado xᵢ(t) (unidades normalizadas)',
             [(i,str(i)) for i in range(6)], [(i/5, f'{i/5:.1f}') for i in range(6)])
    rows = []
    for j in range(101):
        t = j/20
        a = 1/3 + .5*math.exp(-t) + math.exp(-3*t)/6
        b = (1-math.exp(-3*t))/3
        c = 1/3 - .5*math.exp(-t) + math.exp(-3*t)/6
        rows.append((t,a,b,c,a+b+c,a-c))
        assert abs(a+b+c-1)<1e-12 and abs(a-c-math.exp(-t))<1e-12
    assert all(abs(a-b)<1e-13 for a,b in [(rows[0][1],1),(rows[0][2],0),(rows[0][3],0)])
    for i, col in [(1,BLUE),(2,TEAL),(3,AMBER)]:p.curve([(r[0],r[i]) for r in rows],col)
    p.hline(1/3)
    f.legend([('Nodo 1',BLUE),('Nodo 2',TEAL),('Nodo 3',AMBER),('Límite 1/3',MUTED)])
    f.note('Solución exacta del C191.3; suma x₁+x₂+x₃=1 para todo t.')
    return f.bytes(), csv_bytes(['t','x1','x2','x3','suma','x1_menos_x3'],rows)


def build192():
    f = Figure(192, 'Una recuperación compartida cambia el brote',
        'Estrella de dos hojas · tasas de transmisión y recuperación = 1',
        'Probabilidad de cero, una o dos hojas infectadas bajo duración compartida e independencia artificial.')
    p=Plot(f,-.4,2.4,0,.6,'Número de hojas infectadas','Probabilidad (adimensional)',
        [(j,str(j)) for j in range(3)],[(j/10,f'{j/10:.1f}') for j in range(7)])
    shared=[1/3,1/3,1/3]; independent=[1/4,1/2,1/4]
    for j in range(3):
        for val,col,dx in [(shared[j],TEAL,-.18),(independent[j],AMBER,.18)]:
            px,py=p.xy(j+dx,val); _,base=p.xy(0,0)
            f.rect(px-36,py,72,base-py,col,4)
            f.text(px,py-9,f'{val:.3f}',13,INK,anchor='middle')
    assert sum(shared)==1 and sum(independent)==1 and abs(sum(j*shared[j] for j in range(3))-1)<1e-12
    f.legend([('Duración compartida',TEAL),('Enlaces independientes',AMBER)])
    f.note('La media coincide (=1 hoja); P(2) es 1/3 frente a 1/4 (C192.6).')
    return f.bytes(),csv_bytes(['hojas_infectadas','duracion_compartida','enlaces_independientes'],[(j,shared[j],independent[j]) for j in range(3)])


def build193():
    f=Figure(193,'Avalanchas críticas: soporte impar',
        'Descendencia 0 o 2, cada una con probabilidad 1/2 · enumeración exacta',
        'Probabilidades Catalan de tamaños totales impares desde uno a sesenta y uno.')
    p=Plot(f,1,61,1e-4,1,'Tamaño total T (individuos)','P(T) · escala logarítmica',
           [(i,str(i)) for i in [1,11,21,31,41,51,61]],[(10**i,f'10^{i}') for i in [-4,-3,-2,-1,0]],logy=True)
    rows=[]
    for k in range(31):
        t=2*k+1;prob=math.comb(2*k,k)/(k+1)/2**(2*k+1)
        rows.append((k,t,prob))
    assert [rows[i][2] for i in range(4)]==[.5,.125,.0625,5/128]
    p.curve([(r[1],r[2]) for r in rows],BLUE)
    p.dots([(r[1],r[2]) for r in rows],BLUE,3)
    f.legend([('P(T=2k+1), exacta',BLUE)])
    f.note('T pares imposibles; cola mostrada truncada en T=61, no renormalizada.')
    return f.bytes(),csv_bytes(['k','tamano_total_T','probabilidad_exacta'],rows)


def build194():
    f=Figure(194,'Percolación de enlaces en red 3-regular',
        'Límite localmente arbóreo · u=(1−p)/p para p>1/2',
        'Tamaño relativo de componente gigante S=1−u³ frente a retención de enlaces.')
    p=Plot(f,0,1,0,1,'Probabilidad de retener enlace p','Fracción gigante S (población original)',
           [(i/10,f'{i/10:.1f}') for i in range(0,11,2)],[(i/5,f'{i/5:.1f}') for i in range(6)])
    rows=[]
    for j in range(101):
        prob=j/100
        u=1 if prob<=.5 else (1-prob)/prob
        size=1-u**3
        rows.append((prob,u,size))
    assert abs(rows[75][2]-26/27)<1e-12 and all(r[2]==0 for r in rows[:51])
    p.curve([(r[0],r[2]) for r in rows],TEAL)
    p.dots([(.75,26/27)],AMBER,7)
    xx,yy=p.xy(.75,26/27);f.text(xx-12,yy+28,'p=0,75 → S=26/27',13,AMBER,anchor='end')
    f.line(*p.xy(.5,0),*p.xy(.5,1),MUTED,1.5,'5,6')
    f.legend([('Predicción asintótica',TEAL),('Ejemplo p=0,75',AMBER)])
    f.note('Umbral p=1/2; la curva no describe una realización finita (C194.4).')
    return f.bytes(),csv_bytes(['p_retencion','u_no_conexion','S_fraccion_gigante'],rows)


def build195():
    f=Figure(195,'Un enlace adicional detiene una cascada',
        'Cadena A–B–C–D · umbral proporcional 1/2 · semilla A',
        'Conteo exacto de nodos activados por ronda con y sin enlace adicional B–D.')
    p=Plot(f,0,4,0,4.3,'Ronda sincrónica','Nodos activos (de 4)',
           [(i,str(i)) for i in range(5)],[(i,str(i)) for i in range(5)])
    original=[1,2,3,4,4]; extra=[1,1,1,1,1]
    assert original[-1]==4 and extra[-1]==1
    rows=list(zip(range(5),original,extra))
    for arr,col in [(original,TEAL),(extra,AMBER)]:
        p.curve([(i,v) for i,v in enumerate(arr)],col);p.dots([(i,v) for i,v in enumerate(arr)],col)
    f.legend([('Cadena original',TEAL),('Cadena + enlace B–D',AMBER)])
    f.note('La exigencia de B sube de 1 a 2 apoyos cuando su grado pasa de 2 a 3.')
    return f.bytes(),csv_bytes(['ronda','activos_cadena','activos_con_enlace_BD'],rows)


def build196():
    f=Figure(196,'Seis pérdidas por un daño inicial',
        'Capas A y B de cuatro nodos · dependencia mutua y ruta al nodo 1',
        'Cascada determinista causada por eliminar A2; seis fallos secuenciales.')
    p=Plot(f,0,6,0,8,'Paso de la cascada','Nodos fuera de servicio (dos capas)',
           [(i,str(i)) for i in range(7)],[(i,str(i)) for i in range(0,9,2)])
    events=[('A2','daño externo'),('B2','dependencia'),('B3','sin ruta'),('A3','dependencia'),('A4','sin ruta'),('B4','dependencia')]
    rows=[(0,'ninguno','preparación',0,8)]
    rows.extend((i,node,cause,i,8-i) for i,(node,cause) in enumerate(events,1))
    p.curve([(i,i) for i in range(7)],BLUE)
    p.dots([(i,i) for i in range(7)],BLUE,5)
    for i,(node,_) in enumerate(events,1):
        x,y=p.xy(i,i)
        f.text(x+(3 if i<6 else -4),y-13,node,14,BLUE,bold=True,anchor='end' if i==6 else 'middle')
    f.legend([('Pérdidas acumuladas',BLUE)])
    f.note('La capa A aislada conservaría tres nodos tras eliminar A2 (C196.3).')
    assert rows[-1][3:]==(6,2)
    return f.bytes(),csv_bytes(['paso','nodo_perdido','causa','perdidas_acumuladas','nodos_operativos'],rows)


def build197():
    f=Figure(197,'Bloqueo de fase y estabilidad',
        'ω₁=1; ω₂=1,6; κ=0,5 · Δ̇=0,6−sin(Δ)',
        'Campo de fase relativo y equilibrios bloqueados estable e inestable.')
    p=Plot(f,0,math.pi,-.5,.65,'Diferencia de fase Δ (radianes)','Velocidad relativa Δ̇ (rad/tiempo)',
        [(i*math.pi/4,s) for i,s in enumerate(['0','π/4','π/2','3π/4','π'])],
        [(v,f'{v:.1f}') for v in [-.4,-.2,0,.2,.4,.6]])
    rows=[(i*math.pi/200,math.sin(i*math.pi/200),.6-math.sin(i*math.pi/200)) for i in range(201)]
    p.curve([(r[0],r[2]) for r in rows],TEAL)
    p.hline(0)
    stable=math.asin(.6);unstable=math.pi-stable
    assert abs(stable-.6435011087932844)<1e-12 and abs(math.cos(stable)-.8)<1e-12
    for delta,label,color in [(stable,'Estable',BLUE),(unstable,'Inestable',AMBER)]:
        x,y=p.xy(delta,0);f.circle(x,y,6,color);f.text(x,y-15,label,14,color,bold=True,anchor='middle')
    f.legend([('Δ̇=0,6−sin Δ',TEAL),('Bloqueo estable',BLUE),('Inestable',AMBER)])
    f.note('Δ*=arcsin(0,6)=0,6435 rad; recuperación lineal 0,8 por tiempo.')
    return f.bytes(),csv_bytes(['delta_radianes','sin_delta','velocidad_relativa'],rows)


def build198():
    f=Figure(198,'Kuramoto: coherencia no es fracción bloqueada',
        'Distribución lorentziana Δ=1 · umbral de acoplamiento Kc=2',
        'Comparación entre la amplitud de orden r y la fracción de osciladores bloqueados.')
    p=Plot(f,0,6,0,1,'Acoplamiento K (unidades del modelo)','r y fracción bloqueada (adimensional)',
           [(i,str(i)) for i in range(7)],[(i/5,f'{i/5:.1f}') for i in range(6)])
    rows=[]
    for j in range(121):
        k=j/20
        r=0 if k<=2 else math.sqrt(1-2/k)
        fraction=0 if k<=2 else 2/math.pi*math.atan(k*r)
        rows.append((k,r,fraction))
    assert abs(rows[60][1]-1/math.sqrt(3))<1e-12 and abs(rows[60][2]-2/3)<1e-12
    p.curve([(r[0],r[1]) for r in rows],BLUE)
    p.curve([(r[0],r[2]) for r in rows],TEAL)
    p.dots([(3,1/math.sqrt(3))],BLUE,7)
    p.dots([(3,2/3)],TEAL,7)
    f.line(*p.xy(2,0),*p.xy(2,1),MUTED,1.5,'5,6')
    f.legend([('Coherencia r',BLUE),('Fracción bloqueada',TEAL)])
    f.note('Para K=3: r=0,57735 y f_b=2/3; teoría continua, no recuento finito.')
    return f.bytes(),csv_bytes(['K','r_coherencia','fraccion_bloqueada'],rows)


def build199():
    f=Figure(199,'La barrera desaparece cerca del pliegue',
        'ẋ=μ−x² · V(x)=x³/3−μx para μ>0',
        'Barrera local entre el equilibrio estable y el inestable, función de μ.')
    p=Plot(f,0,1,0,1.4,'Parámetro μ (unidades del modelo)','Barrera local ΔV (unidades de potencial)',
           [(i/5,f'{i/5:.1f}') for i in range(6)],[(i*.25,f'{i*.25:.2f}') for i in range(6)])
    rows=[]
    for i in range(101):
        mu=i/100
        x=math.sqrt(mu)
        barrier=4/3*mu**1.5
        rows.append((mu,x,-x,barrier,2*x))
    assert abs(rows[25][3]-1/6)<1e-12 and abs(rows[100][3]-4/3)<1e-12
    p.curve([(r[0],r[3]) for r in rows],BLUE)
    f.legend([('ΔV=(4/3)μ^(3/2)',BLUE)])
    f.note('Barrera de potencial local; el potencial cúbico no confina globalmente.')
    return f.bytes(),csv_bytes(['mu','x_estable','x_inestable','barrera_deltaV','tasa_local'],rows)


def result200(rho):
    # Physical branch selected by repeated iteration starting from q0=rho.
    q=(1-math.sqrt(max(0,(1-9*rho)/(1-rho))))/4 if rho<=1/9 else 1
    s=rho+(1-rho)*(6*q*q-8*q**3+3*q**4)
    return q,s


def build200():
    f=Figure(200,'Cascada por dos apoyos: salto teórico',
        'Red 4-regular · siembra independiente ρ · umbral absoluto de 2 vecinos',
        'Fracción final activa S y transición discontinua en ρc=1/9 del modelo de árbol.')
    p=Plot(f,0,.16,0,1,'Probabilidad inicial ρ','Fracción final activa S',
           [(v,f'{v:.2f}') for v in [0,.04,.08,.12,.16]],[(i/5,f'{i/5:.1f}') for i in range(6)])
    lo=[(i/1000,*result200(i/1000)) for i in range(112) if i/1000<=1/9]
    lo.append((1/9,*result200(1/9)))
    hi=[(1/9,1.,1.),(.16,1.,1.)]
    p.curve([(r,s) for r,q,s in lo],TEAL)
    p.curve([(r,s) for r,q,s in hi],BLUE)
    p.hline(1,color=GRID)
    a=p.xy(1/9,11/32);b=p.xy(1/9,1)
    f.line(*a,*b,AMBER,2,'4,5')
    f.circle(*a,6,TEAL);f.circle(*b,6,BLUE,hollow=True)
    f.text(b[0]+10,b[1]+26,'ρc = 1/9',13,AMBER,bold=True)
    requested=[(.020,.022597),(.050,.068782),(.080,.139416),(.100,.218750),(.105,.250610),(.110,.303917),(.112,1.),(.120,1.),(.150,1.)]
    for rho,expected in requested:
        q,s=result200(rho)
        assert abs(s-expected)<.00000051,(rho,s,expected)
        f.circle(*p.xy(rho,s),4,AMBER)
    rows=[(r,q,s) for r,q,s in lo]+[(r,q,s) for r,q,s in hi[1:]]
    rows += [(rho,*result200(rho)) for rho,_ in requested]
    rows.sort(key=lambda row:row[0])
    assert abs(result200(.1)[1]-7/32)<1e-12 and abs(result200(1/9)[1]-11/32)<1e-12
    f.legend([('Rama de semillas',TEAL),('Activación completa',BLUE),('Tabla original',AMBER)])
    f.note('Datos teóricos del C200; salto del límite izquierdo 11/32 a S=1.')
    return f.bytes(),csv_bytes(['rho_semillas','q_mensaje','S_fraccion_final'],rows)


BUILDERS = {n:globals()[f'build{n}'] for n in range(191,201)}
ARCHIVE='complexa-c191-c200-reproducible.zip'


def generated_assets():
    assets={}
    for n, fn in BUILDERS.items():
        svg,csv=fn()
        assert svg.startswith(b'<svg') and len(svg)>2000 and csv.count(b'\n')>=4
        assets[f'figures/complexa-c{n}-fig01.svg']=svg
        assets[f'repro/c191-200/c{n}.csv']=csv
    return assets


def bundle(assets,root):
    files={f'c{n}.csv':assets[f'repro/c191-200/c{n}.csv'] for n in BUILDERS}
    files.update({f'complexa-c{n}-fig01.svg':assets[f'figures/complexa-c{n}-fig01.svg'] for n in BUILDERS})
    files['make_figures.py']=(root/'repro/c191-200/make_figures.py').read_bytes()
    files['README.md']=(root/'repro/c191-200/README.md').read_bytes()
    b=io.BytesIO()
    with zipfile.ZipFile(b,'w',compression=zipfile.ZIP_DEFLATED,compresslevel=9) as z:
        for name,content in sorted(files.items()):
            meta=zipfile.ZipInfo(name,date_time=(2020,1,1,0,0,0))
            meta.compress_type=zipfile.ZIP_DEFLATED
            meta.external_attr=0o644 << 16
            z.writestr(meta,content,compress_type=zipfile.ZIP_DEFLATED,compresslevel=9)
    return b.getvalue()


def hash_info(name,b):
    return {'path':name,'bytes':len(b),'sha256':hashlib.sha256(b).hexdigest()}


def process(root,check):
    assets=generated_assets()
    archive=bundle(assets,root)
    b64=base64.b64encode(archive)+b'\n'
    assets['zip-data/'+ARCHIVE+'.b64']=b64
    manifest_path=root/'manifest.json'
    manifest=json.loads(manifest_path.read_text(encoding='utf-8'))
    addpaths={*assets,'repro/c191-200/README.md','repro/c191-200/make_figures.py'}
    hashes=[hash_info(name, data) for name,data in assets.items()]
    hashes.extend(hash_info(name,(root/name).read_bytes()) for name in ['repro/c191-200/README.md','repro/c191-200/make_figures.py'])
    manifest['files']=[e for e in manifest['files'] if e['path'] not in addpaths]+sorted(hashes,key=lambda x:x['path'])
    manifest['archives']['C191-C200']={
        'input_path':'zip-data/'+ARCHIVE+'.b64',
        'public_path':'/descargas/'+ARCHIVE,
        'sha256':hashlib.sha256(archive).hexdigest()}
    manifest['publication']='COMPLEXA_reproducibility_2026-10-26_C161-C200'
    # Manifest is intentionally updated only in generate mode; it is validated in check mode.
    if check:
        for name,data in assets.items():
            assert (root/name).read_bytes()==data, f'Not reproducible: {name}'
        index={e['path']:e for e in manifest['files']}
        for entry in hashes:
            assert index[entry['path']]==entry,entry['path']
        assert json.loads(manifest_path.read_text(encoding='utf-8'))==manifest,'Manifest differs'
        print('PASS: 10 SVG, 10 CSV, ZIP, SHA256 and manifest regenerated byte-for-byte')
    else:
        for name,data in assets.items():
            path=root/name;path.parent.mkdir(parents=True,exist_ok=True);path.write_bytes(data)
        manifest_path.write_text(json.dumps(manifest,ensure_ascii=False,indent=2)+'\n',encoding='utf-8')
        print('GENERATED:',len(assets),'assets; ZIP bytes:',len(archive))


def main():
    parser=argparse.ArgumentParser()
    group=parser.add_mutually_exclusive_group(required=True)
    group.add_argument('--generate',action='store_true')
    group.add_argument('--check',action='store_true')
    parser.add_argument('--root',type=Path,default=ROOT)
    args=parser.parse_args()
    process(args.root.resolve(),check=args.check)


if __name__=='__main__':main()