#!/usr/bin/env python3
"""convergencia: corre siru con el paso normal y con la mitad, y compara.

    siru/corpus/convergencia [--tres] [--umbral 0.01] [--plazo SEGUNDOS]
                        [--tabla] CIRCUITO.d|CIRCUITO.ds ...

Pedido del autor (04/10/2026): saber no solo que un circuito corre, sino si
da bien. Esto no dice si el modelo es el correcto (para eso, SCEPTRE en
siru/corpus/sceptre, o una solucion exacta), pero si si el resultado depende del
paso: si al partir el paso a la mitad cambia mas que el umbral, el paso es
grande para ese circuito y la curva no es de fiar.

Para cada circuito:
  - copia su directorio a un temporal, a la misma profundidad que en el
    corpus (los include son relativos: ../../../../slex/...), con slex, siru
    y Simul enlazados; asi el corpus no se toca;
  - corre siru con --paso-factor 1 y 0.5 (y 0.25 con --tres), cada uno en
    su copia, con --salida a un .csv;
  - interpola linealmente la corrida fina en los tiempos de la gruesa y da,
    por variable, la mayor diferencia relativa |gruesa - fina| / max|fina|;
  - con --tres, el orden observado p = log2(|d(1, 1/2)| / |d(1/2, 1/4)|),
    con las diferencias en norma maxima: si el metodo es de orden p, la
    diferencia baja 2^p al partir el paso.

A paso adaptivo siru achica las tolerancias en lugar del paso (ver
Siru::Din, SIRU_PASO_FACTOR). Un circuito sin estados (un punto o un
barrido) no tiene paso que achicar: se informa "sin paso".

Cerca de un salto (una llave que conmuta, una fuente por tramos) la
interpolacion lineal exagera la diferencia: un salto que se corre medio
paso da una diferencia del orden del salto. La columna de la peor variable
dice en que tiempo esta la peor diferencia, para mirarlo.

Salida: una linea por circuito en la salida estandar; --tabla escribe
ademas siru/corpus/CONVERGENCIA.md. Estado de salida: la cantidad de circuitos
que no convergen o no se pudieron correr.

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

AQUI = os.path.dirname(os.path.abspath(__file__))
SIRU = os.path.abspath(os.path.join(AQUI, '..', '..'))
CORPUS = os.path.join(AQUI, 'corpus')
TABLA = os.path.join(AQUI, 'CONVERGENCIA.md')


def copia(fuente, raiz):
    """el directorio de fuente en raiz/<misma ruta respecto de siru>; devuelve
    la ruta del circuito en la copia"""
    fuente = os.path.abspath(fuente)
    rel = os.path.relpath(fuente, SIRU)
    if rel.startswith('..'):
        rel = os.path.join('siru', 'corpus', 'corpus', 'fuera', os.path.basename(fuente))
    destino = os.path.join(raiz, rel)
    shutil.copytree(os.path.dirname(fuente), os.path.dirname(destino),
                    symlinks=True, dirs_exist_ok=True)
    # todo lo de la raiz del repositorio menos siru, que es donde esta el
    # corpus y en la copia tiene solo el directorio del circuito: los include
    # y las rutas relativas de los archivos llegan a slex, siruviejo o
    # tabular2 (con_FOR_1.ds usa ../../../../tabular2/fortran/etable.f90)
    for d in sorted(os.listdir(SIRU)):
        if d == 'siru' or d.startswith('.') or not os.path.isdir(os.path.join(SIRU, d)):
            continue
        if not os.path.lexists(os.path.join(raiz, d)):
            os.symlink(os.path.join(SIRU, d), os.path.join(raiz, d))
    return destino


def correr(fuente, factor, plazo):
    """corre siru con el paso por factor; devuelve (tiempos, {variable: valores}, error)"""
    raiz = tempfile.mkdtemp(prefix='convergencia-')
    try:
        circuito = copia(fuente, raiz)
        csvf = os.path.join(raiz, 'salida.csv')
        siru = os.path.join(SIRU, 'siru', 'bin', 'siru')
        try:
            # SIRU_CALCULA y SIRU_BIBLIOTECAS, como las pruebas de siru/tests:
            # sin ellas se usa el calcula instalado, y entonces un circuito que
            # un calcula recien compilado arregla sigue figurando como que no
            # corre (paso con pru_form y calcula 6.0-36)
            donde = []
            if os.environ.get('SIRU_CALCULA'):
                donde = ['--calcula', os.environ['SIRU_CALCULA']]
                if os.environ.get('SIRU_BIBLIOTECAS'):
                    donde += ['--bibliotecas', os.environ['SIRU_BIBLIOTECAS']]
            # CALCULA_EXTERNAS=1, como siru/corpus/sceptre: los archivos del corpus
            # con funciones en Fortran (fortran/, con_FOR) son propios, y sin
            # el permiso calcula no carga su codigo y siru no arranca
            r = subprocess.run(['perl', siru] + donde + ['--paso-factor', str(factor),
                                '--salida', csvf, os.path.basename(circuito)],
                               cwd=os.path.dirname(circuito), capture_output=True, text=True,
                               timeout=plazo, env=dict(os.environ, CALCULA_EXTERNAS='1'))
        except subprocess.TimeoutExpired:
            return None, None, 'pasa el plazo de %g s' % plazo
        if r.returncode != 0 or not os.path.exists(csvf):
            ult = [l for l in (r.stderr + r.stdout).splitlines() if l.strip()]
            return None, None, 'siru: ' + (ult[-1] if ult else 'estado %d' % r.returncode)
        with open(csvf) as f:
            filas = list(csv.reader(f))
        if len(filas) < 2:
            return None, None, 'sin filas'
        cab = filas[0]
        datos = []
        for fila in filas[1:]:
            try:
                datos.append([float(x) for x in fila])
            except ValueError:
                continue
        # el tiempo es la columna t (o time, tiempo); casi siempre la
        # primera, pero un barrido sin estados (circuito-SE2) escribe todas
        # las variables en orden alfabetico
        it = next((k for k, c in enumerate(cab) if c.lower() in ('t', 'time', 'tiempo')), 0)
        t = [d[it] for d in datos]
        vals = {cab[k]: [d[k] for d in datos] for k in range(len(cab)) if k != it}
        return t, vals, None
    finally:
        shutil.rmtree(raiz, ignore_errors=True)


def interpolar(t, y, x):
    """y(t) lineal en x; t creciente (con tiempos repetidos en un salto,
    toma el ultimo)"""
    lo, hi = 0, len(t) - 1
    if x <= t[0]:
        return y[0]
    if x >= t[-1]:
        return y[-1]
    while hi - lo > 1:
        m = (lo + hi) // 2
        if t[m] <= x:
            lo = m
        else:
            hi = m
    if t[hi] == t[lo]:
        return y[hi]
    return y[lo] + (y[hi] - y[lo]) * (x - t[lo]) / (t[hi] - t[lo])


def diferencia(a, b, tiempos=None):
    """{variable: (relativa, absoluta, tiempo)} de a (gruesa) contra b (fina),
    en los tiempos de a o en los que se den (para el orden, las dos
    diferencias se miden en los mismos tiempos, los de la corrida gruesa)"""
    ta, va = a
    tb, vb = b
    res = {}
    # la escala de cada variable no baja de una milesima de la mayor del
    # circuito (sin las cargas, que llevan la escala de C): una variable que
    # es 0 (IL23 de circuitoBloli, por una simetria) daba una diferencia
    # relativa de 0.83 dividiendo un redondeo por otro, como en siru/corpus/exactas
    general = max((abs(x) for k, xs in vb.items() if not k.lower().startswith('q')
                   for x in xs if math.isfinite(x)), default=0.0)
    for v in va:
        if v not in vb:
            continue
        # los puntos que no son finitos no se comparan: un 0/0 del circuito
        # (PREND = PENCA/PENJ en t = 0, tanque_int) es NaN con cualquier paso
        escala = max(max((abs(x) for x in vb[v] if math.isfinite(x)), default=0.0), 1e-3 * general)
        peor, cuando = 0.0, ta[0]
        lista = tiempos if tiempos is not None else ta
        for i, x in enumerate(lista):
            # En un salto hay dos filas con el mismo tiempo: el estado antes
            # del equilibrio y el de despues (una llave que conmuta). La
            # interpolada de la corrida fina solo puede dar un valor para ese
            # tiempo, el de despues, asi que la fila de antes se saltea: si no,
            # se comparaba el valor previo con el posterior y la diferencia
            # salia del tamanio del salto, 1 en los ocho de llaves_equilibrio,
            # con las dos corridas conmutando en el mismo instante.
            if i + 1 < len(lista) and lista[i + 1] == x:
                continue
            ya = va[v][i] if tiempos is None else interpolar(ta, va[v], x)
            yb = interpolar(tb, vb[v], x)
            if not (math.isfinite(ya) and math.isfinite(yb)):
                continue
            d = abs(ya - yb)
            if d > peor:
                peor, cuando = d, x
        res[v] = (peor / escala if escala > 0 else (0.0 if peor == 0 else math.inf), peor, cuando)
    return res


def main():
    p = argparse.ArgumentParser(description='siru con el paso y con la mitad del paso')
    p.add_argument('archivos', nargs='*')
    p.add_argument('--registros', nargs='+', metavar='REGISTRO',
                   help='no corre nada: arma siru/corpus/CONVERGENCIA.md desde la salida de corridas anteriores '
                        '(varias en paralelo, por ejemplo); si un archivo aparece mas de una vez vale la ultima')
    p.add_argument('--tres', action='store_true', help='tambien un cuarto del paso, y el orden observado')
    p.add_argument('--umbral', type=float, default=0.01)
    p.add_argument('--plazo', type=float, default=900, help='segundos por corrida de siru')
    p.add_argument('--tabla', action='store_true', help='escribir siru/corpus/CONVERGENCIA.md')
    op = p.parse_args()
    if op.registros:
        filas = de_registros(op.registros)
        escribir_tabla(filas, op.umbral)
        print('%s: %d archivos' % (TABLA, len(filas)))
        return 0
    if not op.archivos:
        p.error('faltan los archivos (o --registros)')

    filas, malos = [], 0
    for f in op.archivos:
        nombre = os.path.relpath(os.path.abspath(f), CORPUS) if os.path.abspath(f).startswith(CORPUS + os.sep) else f
        corridas = []
        error = None
        for factor in ((1, 0.5, 0.25) if op.tres else (1, 0.5)):
            t, v, error = correr(f, factor, op.plazo)
            if error:
                break
            corridas.append((t, v))
        if error:
            malos += 1
            print('%s: no corre (%s)' % (nombre, error))
            filas.append((nombre, 'no corre: ' + error, '-', '-', '-'))
            continue
        (t1, v1), (t2, v2) = corridas[0], corridas[1]
        if t1 == t2 and v1 == v2:
            # las mismas filas y los mismos valores: no hay paso que achicar
            # (un punto, un barrido sin estados, o el paso no llega a siru)
            print('%s: sin paso (las dos corridas son identicas, %d filas)' % (nombre, len(t1)))
            filas.append((nombre, 'sin paso', '-', '-', '-'))
            continue
        d12 = diferencia(corridas[0], corridas[1])
        if not d12:
            malos += 1
            print('%s: sin variables en comun' % nombre)
            filas.append((nombre, 'sin variables', '-', '-', '-'))
            continue
        peor = max(d12, key=lambda v: d12[v][0])
        rel, _, cuando = d12[peor]
        orden = '-'
        if op.tres:
            d24 = diferencia(corridas[1], corridas[2], tiempos=t1)
            a = max(x[1] for x in d12.values())
            b = max(x[1] for x in d24.values()) if d24 else 0
            if a > 0 and b > 0:
                orden = '%.1f' % math.log2(a / b)
        ok = rel <= op.umbral
        if not ok:
            malos += 1
        estado = 'converge' if ok else 'no converge'
        print('%s: %s, %s %.2g en t=%g%s' % (nombre, estado, peor, rel, cuando,
                                              '' if orden == '-' else ', orden %s' % orden))
        filas.append((nombre, estado, '%s %.2g' % (peor, rel), '%g' % cuando, orden))

    if op.tabla:
        escribir_tabla(filas, op.umbral)
    return min(malos, 125)


def de_registros(registros):
    """las filas de la tabla desde la salida de corridas anteriores"""
    import re
    filas = {}
    for r in registros:
        for linea in open(r):
            linea = linea.rstrip('\n')
            m = re.match(r'^(\S+): (converge|no converge), (\S+) (\S+) en t=([^,]+)(?:, orden (\S+))?$', linea)
            if m:
                filas[m.group(1)] = (m.group(1), m.group(2), '%s %s' % (m.group(3), m.group(4)),
                                     m.group(5), m.group(6) or '-')
                continue
            m = re.match(r'^(\S+): no corre \((.*)\)$', linea)
            if m:
                filas[m.group(1)] = (m.group(1), 'no corre: ' + m.group(2), '-', '-', '-')
                continue
            m = re.match(r'^(\S+): sin paso', linea)
            if m:
                filas[m.group(1)] = (m.group(1), 'sin paso', '-', '-', '-')
    return [filas[k] for k in sorted(filas)]


def escribir_tabla(filas, umbral):
    cuenta = {}
    for f in filas:
        clave = f[1].split(':')[0]
        cuenta[clave] = cuenta.get(clave, 0) + 1
    with open(TABLA, 'w') as t:
        t.write('# El paso y la mitad del paso\n\n'
                'La escribe `siru/corpus/convergencia` (ver su encabezado): siru con el paso\n'
                'del circuito, con la mitad y con un cuarto, la mayor diferencia relativa\n'
                'por variable (umbral %g) y el orden observado. "converge" no dice que el\n'
                'modelo sea el correcto, solo que el resultado no depende del paso.\n'
                'Lo que significa cada caso que no converge, en `EXACTAS.md`.\n\n'
                '%d archivos: %s.\n\n'
                '| archivo | resultado | peor variable | en t | orden |\n'
                '|---|---|---|---|---|\n'
                % (umbral, len(filas), ', '.join('%d %s' % (n, k) for k, n in
                                                 sorted(cuenta.items(), key=lambda x: -x[1]))))
        for fila in filas:
            t.write('| `%s` | %s | %s | %s | %s |\n' % fila)


if __name__ == '__main__':
    sys.exit(main())
