#!/usr/bin/env python3
"""sceptre: corre SCEPTRE sobre un .d, corre siru sobre el mismo .d, y compara.

    siru/corpus/sceptre [-v] [-k] [--umbral 0.01] [--solo-sceptre]
                   [--calcula EJECUTABLE] [--bibliotecas DIR] CIRCUITO.d ...

Hasta el 30/09/2026 ningun .d del corpus estaba comparado contra SCEPTRE
porque no habia salidas suyas. Esto las produce y compara. Para cada .d deja
al lado dos archivos, los dos tablas de texto con el tiempo primero y una
linea # de nombres, legibles por tabular2/interpola:

  NOMBRE.sceptre   lo que calcula SCEPTRE: la primera corrida de su .sim
  NOMBRE-Sceptre   SCEPTRE contra siru --sobre NOMBRE.ds si lo hay, porque el
                   .d es la entrada de SCEPTRE y el .ds la de siru; si no
                   sobre el .d-- en los tiempos del NOMBRE.patron (si
                   lo hay; si no al 10, 50 y 90 % del intervalo), interpolados
                   con tabular2: por variable de OUTPUTS, SCEPTRE, siru y la
                   diferencia siru - SCEPTRE; al final, la mayor diferencia
                   relativa de cada variable (|dif| / max|SCEPTRE|)

y avisa por la salida de errores de cada variable cuya mayor diferencia
relativa pasa el umbral (--umbral, 0.01 por omision). Estado de salida: la
cantidad de .d que difieren o que no se pudieron comparar. regenera lo corre
con --sceptre.

Y lleva la tabla siru/corpus/SCEPTRE.md: una fila por .d del corpus con si SCEPTRE
lo corre (y si no, su primer error), si siru da igual o distinto, la peor
diferencia, y una columna 'por que' que se escribe a mano y se conserva. Cada
corrida actualiza las filas de los .d que corrio; --tabla la rehace entera con
lo que hay al lado de cada .d (NOMBRE.sceptre, NOMBRE-Sceptre), sin correr
nada.

El SCEPTRE es sceptresimusol (siru/sceptresimusol, el de Simusol: Novender
2017 con la salida .sim). Se compila la primera vez en siru/corpus/sceptre-build,
exportado del arbol para no llenar la fuente de .o. No se usa Nsceptresimusol,
cuya fase 1 escribe nombres de tabla corruptos en el Fortran que genera; ni
novender/SC1-21354, que anda pero no escribe .sim. Como bin/sceptre.sh de
sceptresimusol: el .d en mayusculas (esta version solo entiende tarjetas en
mayusculas), fase 1 (sceptre1: el Fortran del circuito), compilar con
sceptre2.o y correr. Con -k lo de SCEPTRE (listados, Fortran, ejecutable)
queda en el temporal; con -v se muestra lo que dice. --solo-sceptre deja el
NOMBRE.sceptre y no compara.

Fortran agregado: lo que haya en sceptre-fortran/ junto al .d (.f, .f90,
.for) se compila y se enlaza con el circuito. Ahi va lbits.f, que el
analisis AC de la fase 2 llama y libsceptre.a no trae, y las funciones
Fortran externas de los .d. Ver README.md, "Contra SCEPTRE".

Copyright (C) Diego Saravia <dsa@unsa.edu.ar>. GPL version 3 o posterior.
"""
import glob
import os
import re
import shutil
import subprocess
import sys
import tempfile

AQUI = os.path.dirname(os.path.abspath(__file__))
SIRU = os.path.join(AQUI, '..', 'bin', 'siru')
FUENTE = os.path.join(AQUI, '..', '..', 'sceptresimusol')
BUILD = os.path.join(AQUI, 'sceptre-build')
S = os.path.join(BUILD, 'src', 'packages', 'sceptresimusol')
TABLA = os.path.join(AQUI, 'SCEPTRE.md')
CORPUS = os.path.join(AQUI, 'corpus')
ERROR = {}   # .d -> por que no corrio SCEPTRE o siru, para la tabla
FC = ['gfortran', '-fno-automatic', '-finit-local-zero', '-O0', '-std=legacy',
      '-fallow-argument-mismatch']


def compilar_sceptre():
    """sceptresimusol en siru/corpus/sceptre-build, una sola vez"""
    if os.access(os.path.join(S, 'phase1', 'sceptre1'), os.X_OK) \
            and os.path.exists(os.path.join(S, 'phase2', 'libsceptre.a')):
        return True
    if not os.path.isdir(FUENTE):
        print('no esta siru/sceptresimusol', file=sys.stderr); return False
    if not shutil.which('gfortran'):
        print('falta gfortran', file=sys.stderr); return False
    print('compilando sceptresimusol en %s (una sola vez)' % BUILD, file=sys.stderr)
    shutil.rmtree(BUILD, ignore_errors=True)
    if subprocess.run(['svn', 'export', '-q', FUENTE, BUILD]).returncode != 0:
        shutil.copytree(FUENTE, BUILD)
    with open(os.path.join(BUILD, 'make.log'), 'w') as log:
        subprocess.run(['make', 'FFLAGS=' + ' '.join(FC[1:])], cwd=S, stdout=log, stderr=log)
    if not os.access(os.path.join(S, 'phase1', 'sceptre1'), os.X_OK):
        print('no se pudo compilar sceptresimusol: ver %s/make.log' % BUILD, file=sys.stderr)
        return False
    return True


def correr_sceptre(d, verboso, guardar):
    """NOMBRE.sceptre al lado del .d; devuelve su ruta o None"""
    dir_, nom = os.path.split(os.path.abspath(d))
    nom = nom[:-2]
    T = tempfile.mkdtemp(prefix='sceptre-')
    with open(d, encoding='latin-1') as f, open(os.path.join(T, nom + '.d'), 'w', encoding='latin-1') as g:
        g.write(f.read().upper())           # 1. mayusculas
    pasos = [
        ([os.path.join(S, 'phase1', 'sceptre1'), '-o', '.', '-m',
          os.path.join(S, 'extras', 'modellib.dat'), nom], nom + '.ls0'),
    ]
    ok = True
    for orden, log in pasos:
        with open(os.path.join(T, log), 'w') as l:
            ok = subprocess.run(orden, cwd=T, stdout=l, stderr=subprocess.STDOUT,
                                stdin=subprocess.DEVNULL).returncode == 0
    f_for = os.path.join(T, nom + '.for')
    if ok and os.path.getsize(f_for) > 0 if os.path.exists(f_for) else False:
        with open(f_for, encoding='latin-1') as f, open(os.path.join(T, nom + '.f'), 'w', encoding='latin-1') as g:
            g.write(f.read().replace('$', '_'))
        # Fortran agregado: lo que haya en sceptre-fortran/, junto al .d, se
        # compila y se enlaza con el circuito (pedido del autor: LBITS, que
        # el analisis AC de la fase 2 llama y que solo esta en la fase 1, y
        # para probar funciones Fortran externas). No toca el SCEPTRE.
        agregados = sorted(glob.glob(os.path.join(dir_, 'sceptre-fortran', '*.f'))
                           + glob.glob(os.path.join(dir_, 'sceptre-fortran', '*.f90'))
                           + glob.glob(os.path.join(dir_, 'sceptre-fortran', '*.for')))
        # Tablas de tabular2 (ETABLE(N, X, Y, V, M)): si el circuito la llama,
        # se compila tabular2/fortran/etable.f90 y se enlaza con libtabular2
        # (la del arbol, src/.libs, si esta). ETABLE lee etable.txt del
        # directorio de la corrida: se copian el etable.txt que esta junto al
        # .d y los archivos que nombra.
        tabular = []
        with open(os.path.join(T, nom + '.f'), encoding='latin-1') as f:
            # formato fijo: los blancos no cuentan (SCEPTRE escribe 'ETABL  E(')
            usa_etable = 'ETABLE(' in re.sub(r'\s', '', f.read()).upper()
        if usa_etable:
            T2 = os.path.join(AQUI, '..', '..', 'tabular2')
            libs = os.path.join(T2, 'src', '.libs')
            tabular = [os.path.join(T2, 'fortran', 'etable.f90')]
            if os.path.isdir(libs):
                tabular += ['-L' + libs, '-Wl,-rpath,' + os.path.abspath(libs)]
            tabular += ['-ltabular2']
            et = os.path.join(dir_, 'etable.txt')
            if os.path.exists(et):
                shutil.copy(et, T)
                for l in open(et, errors='replace'):
                    c = l.split()
                    if len(c) >= 2 and not l.lstrip().startswith('#'):
                        a = os.path.join(dir_, c[1])
                        if os.path.exists(a):
                            shutil.copy(a, T)
        with open(os.path.join(T, nom + '.ls4'), 'w') as l:
            ok = subprocess.run(FC + ['-o', nom + '_run', nom + '.f'] + agregados + [os.path.join(S, 'phase2', 'sceptre2.o'),
                                      os.path.join(S, 'phase2', 'libsceptre.a'),
                                      os.path.join(S, 'auxpro', 'libsceptreaux.a')] + tabular,
                                cwd=T, stdout=l, stderr=subprocess.STDOUT).returncode == 0
        if ok:
            with open(os.path.join(T, nom + '.ls5'), 'w') as l:
                try:
                    ok = subprocess.run(['./' + nom + '_run'], cwd=T, stdout=l, stderr=subprocess.STDOUT,
                                        stdin=subprocess.DEVNULL, timeout=600).returncode == 0
                except subprocess.TimeoutExpired:
                    ok = False
    else:
        ok = False
    if verboso:
        for log in ('.ls0', '.ls4', '.ls5'):
            p = os.path.join(T, nom + log)
            if os.path.exists(p):
                sys.stderr.write(open(p, errors='replace').read())
    sim = os.path.join(T, nom + '.sim')
    # RUN INITIAL CONDITIONS ONLY: SCEPTRE calcula solo el punto inicial, lo
    # deja en el listado (.ls2, '+I/C TRANSIENT VALUES AT TIME EQUALS ZERO:',
    # una linea 'NOMBRE = 1.23D-01' por variable) y no escribe .sim. Corrio
    # bien: el punto es una tabla de una fila con TIME = 0.
    ls2 = os.path.join(T, nom + '.ls2')
    if ok and (not os.path.exists(sim) or os.path.getsize(sim) == 0) and os.path.exists(ls2):
        texto = open(ls2, errors='replace').read()
        if 'INITIAL CONDITIONS ONLY RUN COMPLETED' in texto:
            # La primera corrida del listado (hasta el segundo 'I/C TRANSIENT
            # VALUES', que es el primer rerun), con sus dos formas: los
            # estados calculados, 'NOMBRE = 1.2D-01' bajo 'I/C TRANSIENT
            # VALUES AT TIME EQUALS ZERO' (si el circuito tiene estados), y
            # los OUTPUTS, '0NOMBRE   1.2D-01' bajo 'RESULTS OF INITIAL
            # CONDITION COMPUTATIONS'. Asi se compara con lo que saca siru,
            # que son los OUTPUTS.
            # Con RERUN hay un bloque por corrida, cada uno con su marca; van
            # al .sceptre como los bloques del .simul de siru: una linea
            # '#  RERUN k/N' antes de cada uno (0 es la maestra).
            marca = 'I/C TRANSIENT VALUES AT TIME EQUALS ZERO:'
            partes = texto.split(marca)
            bloques = []
            for parte in (partes[1:] if len(partes) > 1 else [texto]):
                pares, vistos = [], set()
                for l in parte.split('\n'):
                    m = re.match(r'^0?\s*([A-Za-z]\w*)\s*=\s*([-+0-9.]+[DE][-+]?\d+)\s*$', l) \
                        or re.match(r'^0([A-Za-z]\w*)\s+([-+0-9.]+[DE][-+]?\d+)\s*$', l)
                    if m and m.group(1) not in vistos:
                        vistos.add(m.group(1))
                        pares.append((m.group(1), float(m.group(2).replace('D', 'E'))))
                if pares:
                    bloques.append(pares)
            if bloques:
                ref = os.path.join(dir_, nom + '.sceptre')
                with open(ref, 'w') as s:
                    s.write('# %s.d: lo que calcula SCEPTRE (sceptresimusol); escrito por siru/corpus/sceptre\n' % nom)
                    s.write('# RUN INITIAL CONDITIONS ONLY: solo el punto inicial, del listado\n')
                    for k, pares in enumerate(bloques):
                        if len(bloques) > 1:
                            s.write('#  RERUN %d/%d\n' % (k, len(bloques) - 1))
                        s.write('# TIME ' + ' '.join(n for n, _ in pares) + '\n')
                        s.write('0 ' + ' '.join('%.10g' % v for _, v in pares) + '\n')
                sensibilidad_sceptre(T, nom, dir_)
                montecarlo_sceptre(T, nom, dir_)
                peorcaso_sceptre(T, nom, dir_)
                if not guardar:
                    shutil.rmtree(T, ignore_errors=True)
                return ref
    if not ok or not os.path.exists(sim) or os.path.getsize(sim) == 0:
        print('%s: SCEPTRE no dio salida%s' % (d, ' (queda en %s)' % T if guardar else ''), file=sys.stderr)
        ls1 = os.path.join(T, nom + '.ls1')
        ERROR[os.path.abspath(d)] = 'sin salida'
        if os.path.exists(ls1):
            vistos = set()
            lineas = open(ls1, errors='replace').read().splitlines()
            for i, l in enumerate(lineas):
                if 'ERROR MESSAGE' in l and i + 1 < len(lineas):
                    m = lineas[i + 1].strip('0 ')
                    if m not in vistos and len(vistos) < 5:
                        if not vistos: ERROR[os.path.abspath(d)] = m.strip().lower()
                        vistos.add(m); print('   ' + m, file=sys.stderr)
        if not guardar:
            shutil.rmtree(T, ignore_errors=True)
        return None
    # 4. el .sim como tabla, un bloque por corrida. En el .sim las filas de
    # cada corrida vienen primero y despues su cola: '#RUN/RERUN NUMBER n',
    # '#nam ...' (los nombres), '#ini', '#fin', '#min', '#max'.
    bloques, actual = [], []
    with open(sim, errors='replace') as f:
        for l in f:
            if l.startswith('#RUN/RERUN NUMBER'):
                bloques.append({'filas': actual, 'nombres': None}); actual = []
            elif l.startswith('#nam') and bloques and bloques[-1]['nombres'] is None:
                bloques[-1]['nombres'] = l[4:].split()
            elif l.startswith('#'):
                continue
            elif l.strip():
                actual.append(l.rstrip('\n'))
    if actual:
        bloques.append({'filas': actual, 'nombres': None})
    bloques = [b for b in bloques if b['filas'] and b['nombres']]
    nombres = bloques[0]['nombres'] if bloques else None
    filas = bloques[0]['filas'] if bloques else None
    if not nombres or not filas:
        print('%s: el .sim de SCEPTRE no trae nombres o filas' % d, file=sys.stderr)
        if not guardar:
            shutil.rmtree(T, ignore_errors=True)
        return None
    ref = os.path.join(dir_, nom + '.sceptre')
    with open(ref, 'w') as s:
        s.write('# %s.d: lo que calcula SCEPTRE (sceptresimusol); escrito por siru/corpus/sceptre\n' % nom)
        for k, b in enumerate(bloques):
            if len(bloques) > 1:
                s.write('#  RERUN %d/%d\n' % (k, len(bloques) - 1))
            s.write('# ' + ' '.join(b['nombres']) + '\n')
            s.write('\n'.join(b['filas']) + '\n')
    sensibilidad_sceptre(T, nom, dir_)
    montecarlo_sceptre(T, nom, dir_)
    peorcaso_sceptre(T, nom, dir_)
    if guardar:
        print('%s: lo de SCEPTRE queda en %s' % (d, T), file=sys.stderr)
    else:
        shutil.rmtree(T, ignore_errors=True)
    return ref


def sensibilidad_sceptre(T, nom, dir_):
    """la tabla SENSITIVITY CALCULATION(S) del listado (.ls2), si la hay, a
    NOMBRE.sceptre-dc, con las columnas de NOMBRE.dc de siru. Ojo: en a08.d
    esa tabla no coincide con las propias corridas de SCEPTRE (ver
    siru/PLAN_DC_RERUN.md, paso 3); la de siru si."""
    ls2 = os.path.join(T, nom + '.ls2')
    if not os.path.exists(ls2):
        return
    texto = open(ls2, errors='replace').read()
    if 'SENSITIVITY CALCULATION' not in texto:
        return
    bloque = texto.split('SENSITIVITY CALCULATION', 1)[1].split('\n')
    filas, conj, dep = [], 0, None
    for l in bloque[1:]:
        if re.match(r'^1', l):      # salto de pagina: fin de la tabla
            break
        m = re.match(r'^\s*(\d+)\.\s+(\w+)\s+(\w+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)', l)
        if m:
            conj, dep = int(m.group(1)), m.group(2)
            filas.append((conj, dep, m.group(3), m.group(4), m.group(5))); continue
        m = re.match(r'^\s+(\w+)\s+(\w+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)\s*$', l)
        if m and dep is not None:
            dep = m.group(1)
            filas.append((conj, dep, m.group(2), m.group(3), m.group(4))); continue
        m = re.match(r'^\s+(\w+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)\s*$', l)
        if m and dep is not None:
            filas.append((conj, dep, m.group(1), m.group(2), m.group(3)))
    if not filas:
        return
    with open(os.path.join(dir_, nom + '.sceptre-dc'), 'w') as s:
        s.write('# %s.d: la tabla de sensibilidad de SCEPTRE (SENSITIVITY CALCULATION(S) del listado);\n'
                '# escrito por siru/corpus/sceptre. No coincide con sus propias corridas (a08): ver\n'
                '# siru/PLAN_DC_RERUN.md, paso 3.\n' % nom)
        s.write('#%7s %-12s %-12s %22s %22s\n' % ('CONJ', 'DEPENDIENTE', 'INDEPENDIENTE', 'DERIVADA', 'NORMALIZADA'))
        for f in filas:
            s.write('%8d %-12s %-12s %22s %22s\n' % f)


def montecarlo_sceptre(T, nom, dir_):
    """las OBSERVED STATISTICS del Monte Carlo del listado (.ls2), si las hay, a
    NOMBRE.sceptre-mc, con las columnas del NOMBRE.mc de siru"""
    ls2 = os.path.join(T, nom + '.ls2')
    if not os.path.exists(ls2):
        return
    texto = open(ls2, errors='replace').read()
    if 'OBSERVED STATISTICS' not in texto:
        return
    filas, conj = [], 0
    semilla = re.findall(r'FINAL RANDOM NUMBER\s*=\s*(\d+)', texto)
    for bloque in texto.split('OBSERVED STATISTICS')[1:]:
        conj += 1
        tipo = None
        for l in bloque.split('\n'):
            if re.match(r'^\s*INDEPENDENT', l): tipo = 'IND'; continue
            if re.match(r'^\s*DEPENDENT', l): tipo = 'DEP'; continue
            if re.match(r'^1', l): break
            m = re.match(r'^\s+(\w+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)\s*$', l)
            if m and tipo:
                filas.append((conj, tipo) + m.groups())
    if not filas:
        return
    with open(os.path.join(dir_, nom + '.sceptre-mc'), 'w') as s:
        s.write('# %s.d: las OBSERVED STATISTICS del Monte Carlo de SCEPTRE (del listado); escrito por siru/corpus/sceptre.\n' % nom)
        if semilla: s.write('# semilla final %s\n' % semilla[-1])
        s.write('#%5s %-5s %-12s %20s %20s %20s %20s %20s\n' % ('CONJ', 'TIPO', 'VARIABLE', 'NOMINAL', 'MINIMO', 'MAXIMO', 'MEDIA', 'SIGMA'))
        for f in filas:
            s.write('%6d %-5s %-12s %20s %20s %20s %20s %20s\n' % f)


def peorcaso_sceptre(T, nom, dir_):
    """las WORST CASE COMPUTATION del listado (.ls2), si las hay, a
    NOMBRE.sceptre-wc: por dependiente el valor nominal, HIGH y LOW de SCEPTRE y,
    por independiente, su valor y el gradiente de SCEPTRE en cada punto. No se
    compara con NOMBRE.wc de siru: el gradiente (adjunto) de SCEPTRE esta mal
    (PLAN_DC_RERUN, paso 3), asi que sus puntos extremos son otros; lo de siru
    se verifica corriendo SCEPTRE en los puntos que siru elige (a09wch, a09wcl)."""
    ls2 = os.path.join(T, nom + '.ls2')
    if not os.path.exists(ls2):
        return
    texto = open(ls2, errors='replace').read()
    if 'WORST CASE COMPUTATION' not in texto:
        return
    filas, dep, lado = [], None, None
    for l in texto.split('\n'):
        m = re.match(r'^\+WORST CASE COMPUTATION -\s+(\w+) VALUE', l)
        if m:
            lado = m.group(1); continue
        m = re.match(r'^0OBJECTIVE FUNCTION\s+(\w+)\s*=\s*([-+0-9.]+E[-+]?\d+)', l)
        if m and lado:
            if lado == 'NOMINAL':
                dep = {'nombre': m.group(1), 'NOMINAL': m.group(2), 'ind': {}}
                filas.append(dep)
            if dep:
                dep[lado] = m.group(2)
            continue
        m = re.match(r'^\s+(\w+)\s+([-+0-9.]+E[-+]?\d+)\s+([-+0-9.]+E[-+]?\d+)', l)
        if m and dep and lado:
            dep['ind'].setdefault(m.group(1), {})[lado] = (m.group(2), m.group(3))
    if not filas:
        return
    with open(os.path.join(dir_, nom + '.sceptre-wc'), 'w') as s:
        s.write('# %s.d: las WORST CASE COMPUTATION de SCEPTRE (del listado); escrito por siru/corpus/sceptre.\n' % nom)
        s.write('# Ojo: el gradiente de SCEPTRE es su adjunto, que esta mal (PLAN_DC_RERUN paso 3); sus puntos\n')
        s.write('# HIGH y LOW no son los de siru (NOMBRE.wc) y no se comparan.\n')
        s.write('#%-12s %18s %18s %18s %-12s %18s %18s %18s %18s %18s %18s\n' % ('DEPENDIENTE', 'NOMINAL', 'HIGH', 'LOW',
                'INDEP', 'VALOR', 'GRADIENTE', 'VALOR_HIGH', 'GRAD_HIGH', 'VALOR_LOW', 'GRAD_LOW'))
        for d in filas:
            for n, v in d['ind'].items():
                g = lambda lado, j: v.get(lado, ('-', '-'))[j]
                s.write(' %-12s %18s %18s %18s %-12s %18s %18s %18s %18s %18s %18s\n' % (
                    d['nombre'], d['NOMINAL'], d.get('HIGH', '-'), d.get('LOW', '-'), n,
                    g('NOMINAL', 0), g('NOMINAL', 1), g('HIGH', 0), g('HIGH', 1), g('LOW', 0), g('LOW', 1)))


def leer_tabla(ruta):
    """nombres (la ultima linea # antes de los datos) y filas de numeros"""
    nombres, filas = None, []
    with open(ruta, encoding='utf-8', errors='replace') as f:
        for linea in f:
            s = linea.strip()
            if not s:
                continue
            if s.startswith('#'):
                if not filas:
                    nombres = s.lstrip('#').split()
                continue
            try:
                filas.append([float(x) for x in s.split()])
            except ValueError:
                continue
    if filas and nombres and len(nombres) != len(filas[0]):
        nombres = None
    return nombres, filas


def leer_tabla_simul(ruta):
    """un .simul como tabla: los nombres de sceptrevars y las filas"""
    nombres = nombres_simul(ruta)
    _, filas = leer_tabla(ruta)
    return nombres, filas


def nombres_simul(ruta):
    """la linea de sceptrevars de un .simul: como el .d llama a cada columna"""
    with open(ruta, encoding='utf-8', errors='replace') as f:
        lineas = f.readlines()
    for i, l in enumerate(lineas):
        m = re.match(r'#\s*para Simusol:\s*(.*)', l)
        if m:
            palabras = m.group(1).split()
            for cual in ('sceptrevars', 'variables'):
                if cual in palabras:
                    return lineas[i + 1 + palabras.index(cual)].lstrip('#').split()
    return None


def interpolar(tabla, tiempos, columnas):
    """con tabular2/interpola, lineal: {tiempo: {columna: valor}}"""
    ordenes = ''.join('interpolar %.17g %s\n' % (t, ' '.join(columnas)) for t in tiempos)
    r = subprocess.run(['interpola', tabla], input=ordenes, capture_output=True, text=True)
    valores, t = {}, None
    for l in r.stdout.splitlines():
        m = re.match(r'(\S+)\s*=\s*([-+0-9.eE]+)', l)
        if m:
            t = float(m.group(2)); valores[t] = {}; continue
        m = re.match(r'\s+(\S+)\s+(\S+)\s*$', l)
        if m and t is not None:
            try:
                valores[t][m.group(1)] = float(m.group(2))
            except ValueError:
                pass
    return valores


def bloques_de(ruta, simul):
    """los bloques (corridas) de una tabla: [(titulo, lineas)]. En un .simul de
    siru cada bloque empieza en su cabecera ('#   Archivo ...'); en un .sceptre,
    en su linea '#  RERUN k/N'. Sin reruns hay un solo bloque."""
    lineas = open(ruta, encoding='utf-8', errors='replace').read().split('\n')
    con_datos = lambda ls: any(x.strip() and not x.startswith('#') for x in ls)
    bloques, actual, titulo = [], [], ''
    for l in lineas:
        if (simul and l.startswith('#   Archivo ')) or (not simul and re.match(r'#\s+RERUN \d+/\d+', l)):
            if con_datos(actual):
                bloques.append((titulo, actual))
            actual, titulo = [], ''
        m = re.match(r'#\s+(RERUN \d+/\d+.*)', l)
        if m:
            titulo = m.group(1).strip()
        actual.append(l)
    if con_datos(actual):
        bloques.append((titulo, actual))
    return bloques


def comparar(d, ref, calcula, bibliotecas, umbral):
    """NOMBRE-Sceptre; devuelve las variables que difieren, o None si no se pudo.
    Con RERUN, bloque por bloque: la maestra con la maestra, el rerun k con el
    rerun k; la mayor diferencia de cada variable es sobre todos los bloques."""
    dir_, nom = os.path.split(os.path.abspath(d))
    nom = nom[:-2]
    secciones, peor, escala_global = [], {}, {}
    comunes_todos, faltan_todos = [], set()
    with tempfile.TemporaryDirectory() as T:
        salida = os.path.join(T, nom + '.simul')
        orden = ['perl', SIRU, '--salida', salida]
        if calcula:
            orden += ['--calcula', calcula]
        if bibliotecas:
            orden += ['--bibliotecas', bibliotecas]
        # El .d es la entrada de SCEPTRE y el .ds la de siru (decision del
        # autor, 01/10/2026): si hay NOMBRE.ds al lado, siru corre ese.
        ds = os.path.join(dir_, nom + '.ds')
        fuente_siru = ds if os.path.exists(ds) else d
        # CALCULA_EXTERNAS=1: los archivos del corpus con funciones en Fortran
        # (EXTERNAL) son de confianza; sin el permiso Calcula no las carga
        r = subprocess.run(orden + [fuente_siru], capture_output=True, text=True,
                           env=dict(os.environ, CALCULA_EXTERNAS='1'))
        if r.returncode != 0 or not os.path.exists(salida):
            msg = (r.stdout + r.stderr).strip().splitlines()
            print('%s: siru no corrio: %s' % (d, msg[-1] if msg else ''), file=sys.stderr)
            ERROR[os.path.abspath(d)] = 'siru no corrio: ' + (msg[-1] if msg else '')
            return None
        # Lo que siru escribe al lado de su salida sobre el punto DC (NOMBRE.mc,
        # las estadisticas de Monte Carlo; NOMBRE.wc, el peor caso) queda al
        # lado del .d, para compararlo con lo que siru/corpus/sceptre saco del
        # listado de SCEPTRE.
        for ext in ('.mc', '.wc'):
            if os.path.exists(os.path.join(T, nom + ext)):
                shutil.copy(os.path.join(T, nom + ext), os.path.join(dir_, nom + ext))
        b_ref = bloques_de(ref, False)
        b_sim = bloques_de(salida, True)
        if not b_ref:
            print('%s: %s sin tabla' % (d, ref), file=sys.stderr); return None
        if not b_sim:
            print('%s: el .simul no trae filas' % d, file=sys.stderr); return None
        # Los bloques se emparejan por su numero de RERUN (un rerun que no
        # corrio en uno de los dos no desplaza a los demas); sin numero, por
        # posicion.
        def numero(titulo, k):
            m = re.match(r'RERUN (\d+)/', titulo or '')
            return int(m.group(1)) if m else k
        por_num_sim = {numero(t, k): (t, ls) for k, (t, ls) in enumerate(b_sim)}
        pares = [(numero(t, k), (t, ls), por_num_sim.get(numero(t, k))) for k, (t, ls) in enumerate(b_ref)]
        faltan_sim = [k for k, _, b in pares if b is None]
        if faltan_sim or len(b_ref) != len(b_sim):
            print('%s: SCEPTRE trae %d corridas y siru %d%s' % (d, len(b_ref), len(b_sim),
                  '; siru no tiene la ' + ', '.join(str(k) for k in faltan_sim) if faltan_sim else ''), file=sys.stderr)
        patron = os.path.join(dir_, nom + '.patron')
        tiempos_patron = None
        if os.path.exists(patron):
            _, filas_p = leer_tabla(patron)
            tiempos_patron = [f[0] for f in filas_p] or None
        for k, (titulo_ref, lineas_ref), bs in pares:
            if bs is None:
                continue
            titulo_sim, lineas_sim = bs
            ref_k = os.path.join(T, 'ref%d.sceptre' % k)
            sim_k = os.path.join(T, 'sim%d.simul' % k)
            open(ref_k, 'w').write('\n'.join(lineas_ref) + '\n')
            open(sim_k, 'w').write('\n'.join(lineas_sim) + '\n')
            nombres_ref, filas_ref = leer_tabla(ref_k)
            if not nombres_ref or not filas_ref:
                print('%s: el bloque %d de %s no trae tabla' % (d, k, ref), file=sys.stderr); continue
            nombres_sim = nombres_simul(sim_k)
            if not nombres_sim:
                print('%s: el bloque %d del .simul no trae nombres' % (d, k), file=sys.stderr); continue
            por_mayus = {n.upper(): n for n in nombres_sim[1:]}
            comunes = [(n, por_mayus[n.upper()]) for n in nombres_ref[1:] if n.upper() in por_mayus]
            faltan = [n for n in nombres_ref[1:] if n.upper() not in por_mayus]
            faltan_todos.update(faltan)
            if not comunes:
                print('%s: ninguna variable en comun (SCEPTRE %s, siru %s)'
                      % (d, nombres_ref[1:], nombres_sim[1:]), file=sys.stderr)
                if k == 0:
                    return None
                continue
            if k == 0:
                comunes_todos = comunes
            t0, t1 = filas_ref[0][0], filas_ref[-1][0]
            tiempos = tiempos_patron or [t0 + (t1 - t0) * p for p in (0.1, 0.5, 0.9)]
            # un rerun puede tener otro STOP TIME: solo los tiempos que estan
            # en los dos bloques
            tiempos = [t for t in tiempos if t0 - 1e-9 * abs(t1 - t0) <= t <= t1 + 1e-9 * abs(t1 - t0)] or [t0 + (t1 - t0) * p for p in (0.1, 0.5, 0.9)]
            # Con una sola fila (RUN INITIAL CONDITIONS ONLY, un punto) no hay
            # nada que interpolar: se compara la fila misma.
            def valores(tabla, columnas):
                nombres_t, filas_t = leer_tabla_simul(tabla) if tabla.endswith('.simul') else leer_tabla(tabla)
                if nombres_t and len(filas_t) == 1:
                    fila = dict(zip(nombres_t, filas_t[0]))
                    return {fila[nombres_t[0]]: {c: fila[c] for c in columnas if c in fila}}
                return interpolar(tabla, tiempos, columnas)
            if len(filas_ref) == 1:
                tiempos = [filas_ref[0][0]]
            v_ref = valores(ref_k, [n for n, _ in comunes])
            v_sim = valores(sim_k, [m for _, m in comunes])
            for n, _ in comunes:
                jn = nombres_ref.index(n)
                escala_global[n] = max(escala_global.get(n, 0.0), max(abs(f[jn]) for f in filas_ref))
                peor.setdefault(n, 0.0)
            # las filas de este bloque: (titulo, k, fila); y la mayor
            # diferencia de cada variable, sobre todos los bloques
            for t in tiempos:
                tr = min(v_ref, key=lambda x: abs(x - t)) if v_ref else None
                ts = min(v_sim, key=lambda x: abs(x - t)) if v_sim else None
                fila = ['%22.15g' % t]
                for n, m in comunes:
                    a = v_ref.get(tr, {}).get(n) if tr is not None else None
                    b = v_sim.get(ts, {}).get(m) if ts is not None else None
                    if a is None or b is None:
                        fila += ['%22s' % '-'] * 3
                        continue
                    dif = b - a
                    fila += ['%22.10g' % a, '%22.10g' % b, '%22.3g' % dif]
                    peor[n] = max(peor[n], abs(dif) / (escala_global[n] or 1.0))
                secciones.append((titulo_ref or titulo_sim, k, fila))
    comunes = comunes_todos
    if not comunes:
        return None
    with open(os.path.join(dir_, nom + '-Sceptre'), 'w', encoding='utf-8') as s:
        s.write('# %s.d en SCEPTRE (sceptresimusol) contra %s en siru, interpolados con tabular2\n'
                '# en los tiempos %s; escrito por siru/corpus/sceptre.\n'
                '# Por variable: lo de SCEPTRE, lo de siru y la diferencia siru - SCEPTRE.\n'
                % (nom, os.path.basename(fuente_siru),
                   'del patron' if tiempos_patron else 'al 10, 50 y 90 % del intervalo'))
        if faltan_todos:
            s.write('# Sin comparar, siru no las saca: %s\n' % ' '.join(sorted(faltan_todos)))
        cab = ['TIME']
        for n, _ in comunes:
            cab += [n, n + '_siru', n + '_dif']
        ultimo = None
        corridas = len(set(x[1] for x in secciones))
        for titulo, k, fila in secciones:
            if k != ultimo:
                if corridas > 1:
                    s.write('#  %s\n' % (titulo or 'RERUN %d' % k))
                s.write('# ' + ' '.join('%22s' % c for c in cab).lstrip() + '\n')
                ultimo = k
            s.write(' '.join(fila) + '\n')
        s.write('# mayor diferencia relativa (|dif| / max|SCEPTRE|):'
                + ''.join(' %s %.2g' % (n, peor[n]) for n, _ in comunes) + '\n')
    difieren = [n for n, _ in comunes if peor[n] > umbral]
    mc = comparar_mc(d, dir_, nom, umbral)
    un_punto = all(len([x for x in secciones if x[1] == k]) == 1 for k in set(x[1] for x in secciones))
    if difieren:
        print('%s: difiere de SCEPTRE en %s%s' % (d, ' '.join('%s (%.2g)' % (n, peor[n]) for n in difieren),
              '; ojo, tras RUN MONTE CARLO los OUTPUTS de SCEPTRE son los de la ultima muestra, los de siru los nominales'
              if mc is not None and un_punto else ''), file=sys.stderr)
    else:
        print('%s: igual a SCEPTRE (umbral %g%s)' % (d, umbral, ', %d corridas' % corridas if corridas > 1 else ''))
    return difieren


def comparar_mc(d, dir_, nom, umbral):
    """NOMBRE.mc (siru) contra NOMBRE.sceptre-mc (SCEPTRE), fila por fila: por
    variable el nominal, el minimo, el maximo, la media y el sigma. Devuelve la
    lista de lo que difiere, o None si falta alguna de las dos tablas."""
    a, b = os.path.join(dir_, nom + '.mc'), os.path.join(dir_, nom + '.sceptre-mc')
    if not (os.path.exists(a) and os.path.exists(b)):
        return None
    def leer(ruta):
        filas = {}
        for l in open(ruta, encoding='utf-8'):
            c = l.split()
            if not c or c[0].startswith('#') or len(c) < 8:
                continue
            filas[(c[0], c[1], c[2].upper())] = [float(x) for x in c[3:8]]
        return filas
    fa, fb = leer(a), leer(b)
    columnas = ['NOMINAL', 'MINIMO', 'MAXIMO', 'MEDIA', 'SIGMA']
    difieren = []
    for k in fb:
        if k not in fa:
            difieren.append('%s (siru no la da)' % k[2]); continue
        for j, col in enumerate(columnas):
            escala = max(abs(fb[k][0]), abs(fb[k][j]), 1e-300)
            r = abs(fa[k][j] - fb[k][j]) / escala
            if r > umbral:
                difieren.append('%s.%s (%.2g)' % (k[2], col, r))
    if difieren:
        print('%s: Monte Carlo difiere de SCEPTRE en %s' % (d, ' '.join(difieren)), file=sys.stderr)
    else:
        print('%s: Monte Carlo igual a SCEPTRE (%d filas, umbral %g)' % (d, len(fb), umbral))
    return difieren


def tabla_leer():
    """las filas de siru/corpus/SCEPTRE.md: {archivo: [sceptre, siru, peor, por que]}"""
    filas = {}
    if not os.path.exists(TABLA):
        return filas
    for l in io_abrir(TABLA):
        c = [x.strip() for x in l.strip().strip('|').split('|')]
        if len(c) == 5 and c[0].endswith('.d`'):
            filas[c[0].strip('`')] = c[1:]
    return filas


def io_abrir(ruta):
    with open(ruta, encoding='utf-8') as f:
        return f.readlines()


def tabla_fila(d, umbral, anterior):
    """la fila de un .d con lo que hay a su lado; anterior conserva el error
    de SCEPTRE (que no queda en ningun archivo) y el por que"""
    base = d[:-2]
    sc, si, peor = (anterior or ['sin correr', '-', '-', ''])[:3]
    por_que = (anterior or [''] * 4)[3]
    err = ERROR.get(os.path.abspath(d))
    if os.path.exists(base + '.sceptre'):
        sc = 'corre'
    elif err and not err.startswith('siru'):
        sc = 'no corre: ' + err
    if os.path.exists(base + '-Sceptre') and os.path.exists(base + '.sceptre'):
        ultima = io_abrir(base + '-Sceptre')[-1]
        pares = re.findall(r'(\S+) ([0-9.eE+-]+)', ultima.split(':', 1)[-1])
        if pares:
            n, v = max(pares, key=lambda p: float(p[1]))
            peor = '%s %s' % (n, v)
            si = 'igual' if float(v) <= umbral else 'distinto'
    elif sc != 'corre':
        si, peor = '-', '-'
    if err and err.startswith('siru'):
        si, peor = err.replace('|', '/'), '-'
    return [sc.replace('|', '/'), si, peor, por_que]


def tabla_escribir(filas):
    corren = sum(1 for f in filas.values() if f[0] == 'corre')
    iguales = sum(1 for f in filas.values() if f[1] == 'igual')
    distintos = sum(1 for f in filas.values() if f[1] == 'distinto')
    with open(TABLA, 'w', encoding='utf-8') as s:
        s.write('# Los .d del corpus contra SCEPTRE\n\n'
                'La escribe `siru/corpus/sceptre` (cada corrida actualiza sus filas; `--tabla` la\n'
                'rehace con lo que hay al lado de cada .d). La columna **por que** se escribe\n'
                'a mano y se conserva. Lo que se entendio de los casos, en `README.md`,\n'
                '\u00abContra SCEPTRE\u00bb. SCEPTRE corre el `.d`; siru, el `.ds` si lo hay.\n\n'
                '%d .d: SCEPTRE corre %d; siru igual en %d, distinto en %d.\n\n'
                % (len(filas), corren, iguales, distintos))
        s.write('| archivo | SCEPTRE | siru | peor diferencia | por que |\n|---|---|---|---|---|\n')
        for k in sorted(filas):
            s.write('| `%s` | %s |\n' % (k, ' | '.join(filas[k])))


def tabla_actualizar(archivos, umbral, solo_estos=False):
    filas = tabla_leer()
    if solo_estos:   # --tabla: fuera las filas de los .d que ya no estan en el corpus
        hay = {os.path.relpath(os.path.abspath(d), CORPUS) for d in archivos}
        filas = {k: v for k, v in filas.items() if k in hay}
    for d in archivos:
        k = os.path.relpath(os.path.abspath(d), CORPUS)
        filas[k] = tabla_fila(d, umbral, filas.get(k))
    tabla_escribir(filas)


def main(args):
    calcula = bibliotecas = None
    verboso = guardar = solo = False
    umbral = 0.01
    archivos = []
    i = 0
    while i < len(args):
        a = args[i]
        if a == '--calcula':
            calcula = args[i + 1]; i += 2
        elif a == '--bibliotecas':
            bibliotecas = args[i + 1]; i += 2
        elif a == '--umbral':
            umbral = float(args[i + 1]); i += 2
        elif a == '-v':
            verboso = True; i += 1
        elif a == '-k':
            guardar = True; i += 1
        elif a == '--solo-sceptre':
            solo = True; i += 1
        elif a == '--tabla':
            todos = sorted(os.path.join(r, f) for r, _, fs in os.walk(CORPUS) for f in fs if f.endswith('.d'))
            tabla_actualizar(todos, umbral, solo_estos=True)
            return 0
        elif a in ('-h', '--help'):
            print(__doc__); return 0
        else:
            archivos.append(a); i += 1
    if not archivos:
        print(__doc__, file=sys.stderr); return 2
    if not compilar_sceptre():
        return 77
    mal = 0
    for d in archivos:
        if not d.endswith('.d') or not os.path.exists(d):
            print('%s: no es un .d que exista' % d, file=sys.stderr); mal += 1; continue
        ref = correr_sceptre(d, verboso, guardar)
        if ref is None:
            mal += 1; continue
        if solo:
            continue
        r = comparar(d, ref, calcula, bibliotecas, umbral)
        if r is None or r:
            mal += 1
    # solo los del corpus van a la tabla
    del_corpus = [d for d in archivos if os.path.abspath(d).startswith(CORPUS + os.sep)]
    if del_corpus:
        tabla_actualizar(del_corpus, umbral)
    return min(mal, 125)


if __name__ == '__main__':
    sys.exit(main(sys.argv[1:]))
