#!/usr/bin/env python3
"""exactas: la solucion exacta de los circuitos del corpus que la tienen.

    siru/corpus/exactas                         escribe NOMBRE.exacta de todos
    siru/corpus/exactas --escribir ARCHIVO ...  escribe NOMBRE.exacta de esos
    siru/corpus/exactas --comparar [--tabla] [--plazo S] [ARCHIVO ...]
                                           corre siru y lo compara
    siru/corpus/exactas --lista                 los circuitos que tiene
    siru/corpus/exactas --buscar                los .d y .ds del corpus que son
                                           circuitos lineales simples

Pedido del autor (04/10/2026): saber no solo que un circuito corre, sino si
da bien. Para los circuitos lineales chicos la solucion se escribe a mano y
no depende de siru ni de SCEPTRE: es la referencia independiente. Las
ecuaciones de cada uno, de donde salen y como se resuelven estan en
siru/corpus/EXACTAS.md, y tambien en el encabezado de cada NOMBRE.exacta.

Los tiempos comunes. La idea del corpus es que siru (el .patron), SCEPTRE
(siru/corpus/sceptre) y SPICE dejen sus valores en los mismos tiempos para
compararlos; la exacta se suma a esos. Los tiempos son los del .patron; si
no hay, los del -Sceptre; si tampoco, el 10, 50 y 90 % del intervalo (como
siru/corpus/sceptre). El .patron lo escribio siru: no dice si siru da bien, dice
si cambio; contra la exacta, si.

NOMBRE.exacta (al lado del .d o del .ds): la exacta en esos tiempos, una
tabla de texto con el tiempo primero y una linea '# TIME ...' con los
nombres, como el .patron y el NOMBRE.sceptre, legible por tabular2/interpola;
el encabezado lleva las ecuaciones. Los nombres son los de SCEPTRE en
mayusculas (VC3, IL23, IR5, QC3); siru y el .patron los escriben en
minusculas, y la comparacion no distingue.

--comparar corre siru sobre una copia del directorio del circuito (como
siru/corpus/convergencia, para no tocar el corpus) y:
  - evalua la formula exacta en todos los tiempos de la salida de siru y
    da, por variable, la mayor diferencia absoluta y la relativa
    |siru - exacta| / max|exacta| (la salida estandar y, con --tabla,
    siru/corpus/EXACTAS-comparacion.md);
  - escribe NOMBRE-Exacta, como el NOMBRE-Sceptre de siru/corpus/sceptre: en los
    tiempos comunes, por variable, la exacta y, por cada fuente que haya
    --siru corrido ahora, el .patron, SCEPTRE del -Sceptre--, su valor y la
    diferencia con la exacta; al final la mayor diferencia relativa de cada
    una. Una columna de SPICE se agregara igual cuando haya un circuito con
    exacta y .cir.

De donde sale cada exacta. Unos pocos circuitos tienen la formula escrita
a mano (el tanque, el pendulo de c01 con funciones elipticas). El resto los
resuelve red_lineal, que arma la solucion de cualquier circuito lineal de
R, L, C, E y J --constantes, controladas, senoidales o por tabla del
tiempo-- sin escribir sus ecuaciones (ver su comentario: los cortes de
bobinas y los lazos de capacitores, el arranque desparejo, las fuentes como
seudoestados). leer_lineal lee esos circuitos de su .d o .ds cuando son
simples (R, L, C, E y J de valor numerico, sin modelos, llaves ni
ecuaciones propias); AUTOMATICOS son los del corpus que se resuelven asi,
y --buscar lista los candidatos.

Las soluciones de los sistemas lineales x' = A x + b son
x(t) = xs + expm(A t) (x0 - xs), con xs = -A^-1 b el estado de reposo y
expm la exponencial de matriz de scipy: exacta a la precision de la maquina,
sin paso de integracion. Requiere numpy y scipy (solo para probar; siru no
los necesita).

Copyright (C) Diego Saravia <dsa@unsa.edu.ar>. GPL version 3 o posterior.
"""
import argparse
import importlib.machinery
import importlib.util
import math
import os
import re
import sys

import numpy as np
from scipy.linalg import expm, null_space

AQUI = os.path.dirname(os.path.abspath(__file__))
CORPUS = os.path.join(AQUI, 'corpus')
TABLA = os.path.join(AQUI, 'EXACTAS-comparacion.md')


def lineal(A, b, x0):
    """t -> x(t) de x' = A x + b, x(0) = x0, por la exponencial de matriz"""
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)
    x0 = np.array(x0, dtype=float)
    xs = -np.linalg.solve(A, b)
    return lambda t: xs + expm(A * t) @ (x0 - xs)


def _equilibrar(M):
    """filas y columnas de M a norma 1: (S, D) con S M D equilibrada"""
    r = np.linalg.norm(M, axis=1); r[r == 0] = 1
    S = 1 / r
    c = np.linalg.norm(M * S[:, None], axis=0); c[c == 0] = 1
    return S, 1 / c


def red_lineal(elementos, iniciales, reposo=False, t0=0.0):
    """la solucion de un circuito lineal de R, L, C, E y J, sin escribir sus
    ecuaciones a mano. elementos: [(nombre, a, b, valor)], el tipo por la
    letra; convencion de SCEPTRE: la corriente va de a a b por el elemento,
    V = Ua - Ub, una E sube el potencial de a a b (V = -E) y una J es la
    corriente. El valor de una E o una J puede ser un numero, una fuente
    controlada ('V', rama, k) o ('I', rama, k) --k por la tension o la
    corriente de otra rama--, una senoidal ('sen', A, w), A sin(w t), o una
    tabla del tiempo ('tabla', [(t, v), ...]), lineal por tramos, con un t
    repetido como salto y el valor del extremo fuera de la tabla (como
    tabular2).
    Estados: las corrientes de L y las tensiones de C.

    Con los estados fijos el circuito es algebraico, M x = r0 + R s. Los
    cortes de solo bobinas y fuentes de corriente, y los lazos de solo
    capacitores y fuentes de tension, son el espacio nulo izquierdo de M:
    restricciones G s = g de los estados. La derivada sale de M x = r0 + R s,
    VL = L dIL/dt, IC = C dVC/dt y G ds = 0, por cuadrados minimos con
    residuo nulo. Un arranque que no cumple G s = g se proyecta conservando
    el flujo de los cortes y la carga de los lazos (la metrica de L y C: el
    impulso de tension, o de corriente, que lo equilibra; lo mismo que hace
    siru). Una senoidal se genera con dos seudoestados, p = A sin(w t) y
    q = A cos(w t), con p' = w q y q' = -w p; una tabla con uno, p, cuya
    derivada es la pendiente del tramo. Asi todo es x' = A x + b en las
    coordenadas libres, con b constante en cada tramo de las tablas, y se
    resuelve por la exponencial de matriz tramo a tramo, sin paso de
    integracion; en un salto de la tabla p cambia en el acto. sol(t) recibe
    el tiempo absoluto; la corrida empieza en t0. Con reposo (RUN INITIAL CONDITIONS) se arranca del
    estado de reposo con las fuentes en su valor de t0."""
    nodos = sorted({n for _, a, b, _ in elementos for n in (a, b)} - {'GND'})
    idx = {n: k for k, n in enumerate(nodos)}
    ne, nn = len(elementos), len(nodos)
    nombres = [e[0] for e in elementos]
    est = [e for e in elementos if e[0][0] in 'LC']
    ns = len(est); pos = {e[0]: j for j, e in enumerate(est)}
    senos = [(n, v[1], v[2]) for n, a, b, v in elementos if isinstance(v, tuple) and v[0] == 'sen']
    tablas = [(n, [(float(x), float(y)) for x, y in v[1]]) for n, a, b, v in elementos
              if isinstance(v, tuple) and v[0] == 'tabla']
    nsen = 2 * len(senos)
    npse = nsen + len(tablas)                  # p y q de cada senoidal, p de cada tabla
    nx = nn + ne

    def ua(a, b):
        u = np.zeros(nx)
        if a != 'GND': u[idx[a]] += 1
        if b != 'GND': u[idx[b]] -= 1
        return u

    def de_rama(r):                             # la fila que da V o I de otra rama
        k = nombres.index(r)
        return k, ua(elementos[k][1], elementos[k][2])
    # sistema estatico: M x = r0 + R [s; pseudo]
    M = np.zeros((nx, nx)); r0 = np.zeros(nx); R = np.zeros((nx, ns + npse))
    for k, (n, a, b, v) in enumerate(elementos):
        if a != 'GND': M[idx[a], nn + k] += 1
        if b != 'GND': M[idx[b], nn + k] -= 1
    js, jt = 0, 0
    for k, (n, a, b, v) in enumerate(elementos):
        f = nn + k; t = n[0]
        if isinstance(v, tuple) and t not in 'EJ':
            raise ValueError('rama %s: solo una E o una J puede ser controlada o senoidal' % n)
        if t == 'R': M[f] = ua(a, b); M[f, nn + k] -= v
        elif t == 'L': M[f, nn + k] = 1; R[f, pos[n]] = 1
        elif t == 'C': M[f] = ua(a, b); R[f, pos[n]] = 1
        elif t == 'E':                          # Ua - Ub = -E
            M[f] = ua(a, b)
            if not isinstance(v, tuple): r0[f] = -v
            elif v[0] == 'sen': R[f, ns + 2 * js] = -1; js += 1
            elif v[0] == 'tabla': R[f, ns + nsen + jt] = -1; jt += 1
            else:
                kr, u = de_rama(v[1])
                if v[0] == 'V': M[f] += v[2] * u
                else: M[f, nn + kr] += v[2]
        elif t == 'J':                          # I = J
            M[f, nn + k] = 1
            if not isinstance(v, tuple): r0[f] = v
            elif v[0] == 'sen': R[f, ns + 2 * js] = 1; js += 1
            elif v[0] == 'tabla': R[f, ns + nsen + jt] = 1; jt += 1
            else:
                kr, u = de_rama(v[1])
                if v[0] == 'V': M[f] -= v[2] * u
                else: M[f, nn + kr] -= v[2]
        else:
            raise ValueError('rama %s: tipo %s' % (n, t))
    # restricciones de los estados: w^T M = 0; con M equilibrada (S M D) el
    # espacio nulo izquierdo es S w', y el corte numerico no confunde una
    # escala chica con una dependencia (a04: C = 5e6 y 5e-4 juntos)
    Sm, Dm = _equilibrar(M)
    W = null_space((M * Sm[:, None] * Dm[None, :]).T, rcond=1e-10)
    W = Sm[:, None] * W if W.size else W
    G = W.T @ R if W.size else np.zeros((0, ns + npse)); g = -(W.T @ r0) if W.size else np.zeros(0)
    keep = np.linalg.norm(G, axis=1) > 1e-12 if len(G) else []
    if len(G) and np.any(~keep) and np.any(np.abs(g[~keep]) > 1e-9):
        raise ValueError('fuentes incompatibles')
    G, g = (G[keep], g[keep]) if len(G) else (G, g)
    if len(G) and np.abs(G[:, ns:]).max(initial=0) > 1e-12:
        raise ValueError('una senoidal o una tabla en un lazo de capacitores o un corte de bobinas')
    G = G[:, :ns]
    peso = np.array([e[3] for e in est])
    # derivadas: [x ; ds] con M x = r0 + R s, VL = L dsL, IC = C dsC, G ds = 0
    Z = np.zeros((nx + ns + len(G), nx + ns))
    Z[:nx, :nx] = M
    for j, (n, a, b, v) in enumerate(est):
        k = nombres.index(n)
        if n[0] == 'L': Z[nx + j, :nx] = ua(a, b) / v; Z[nx + j, nx + j] = -1   # dI/dt = VL/L
        else: Z[nx + j, nn + k] = 1 / v; Z[nx + j, nx + j] = -1                 # dV/dt = IC/C
    Z[nx + ns:, nx:] = G
    Sz, Dz = _equilibrar(Z)
    Ze = Z * Sz[:, None] * Dz[None, :]
    if np.linalg.matrix_rank(Ze, tol=1e-10) < nx + ns: raise ValueError('el sistema no se determina')

    def resolver(sp):
        rhs = np.concatenate([r0 + R @ sp, np.zeros(ns + len(G))])
        ye, *_ = np.linalg.lstsq(Ze, Sz * rhs, rcond=1e-12)
        y = Dz * ye
        res = np.linalg.norm(Z @ y - rhs)
        assert res < 1e-8 * (1 + np.linalg.norm(rhs)), 'inconsistente (residuo %g de %g)' % (res, np.linalg.norm(rhs))
        return y[:nx], y[nx:]
    # el arranque: el estado inicial proyectado sobre G s = g con la metrica L, C
    s0 = np.array([iniciales.get(e[0], 0.0) for e in est])
    if len(G):
        Wi = np.diag(1 / peso)
        s0 = s0 + Wi @ G.T @ np.linalg.solve(G @ Wi @ G.T, g - G @ s0)
    def valor_tabla(pts, t):
        """la tabla en t, lineal entre puntos; en un salto (x repetido) el
        valor de despues; fuera de la tabla, el del extremo"""
        if t < pts[0][0]: return pts[0][1]
        i = max(k for k, (x, _) in enumerate(pts) if x <= t)
        if i == len(pts) - 1: return pts[-1][1]
        (x0, y0), (x1, y1) = pts[i], pts[i + 1]
        return y0 + (y1 - y0) * (t - x0) / (x1 - x0)

    def pendiente(pts, a, b):
        """la pendiente de la tabla en el tramo (a, b), sin quiebres adentro"""
        c = (a + b) / 2
        for (x0, y0), (x1, y1) in zip(pts, pts[1:]):
            if x0 <= c < x1 and x1 > x0:
                return (y1 - y0) / (x1 - x0)
        return 0.0
    pse0 = np.concatenate([[a_ * math.sin(w * t0), a_ * math.cos(w * t0)] for _, a_, w in senos]
                          + [[valor_tabla(pts, t0)] for _, pts in tablas]) if npse else np.zeros(0)
    N = null_space(G) if len(G) else np.eye(ns)
    m = N.shape[1]
    rot = np.zeros((npse, npse))
    for j, (_, a_, w) in enumerate(senos):
        rot[2 * j, 2 * j + 1] = w; rot[2 * j + 1, 2 * j] = -w

    def campo(y):                               # y = [z (m); pseudo (npse)]
        z, pse = y[:m], y[m:]
        _, d = resolver(np.concatenate([s0 + N @ z, pse]))
        return np.concatenate([N.T @ d, rot @ pse])
    nt = m + npse
    b = campo(np.zeros(nt)); A = np.column_stack([campo(np.eye(nt)[j]) - b for j in range(nt)]) if nt else np.zeros((0, 0))
    y0 = np.concatenate([np.zeros(m), pse0])
    if reposo and m:                            # dz = 0 con las fuentes en t0
        z = np.linalg.lstsq(A[:m, :m], -(b[:m] + A[:m, m:] @ pse0), rcond=None)[0]
        y0[:m] = z
    # los quiebres de las tablas, despues de t0; en cada uno se guarda el estado
    quiebres = sorted({x for _, pts in tablas for x, _ in pts if x > t0})

    def avanzar(y, a, b_):
        """y de a a b_, con las pendientes de las tablas en ese tramo"""
        bb = b.copy()
        for j, (_, pts) in enumerate(tablas):
            bb[m + nsen + j] = pendiente(pts, a, b_)
        Aa = np.zeros((nt + 1, nt + 1)); Aa[:nt, :nt] = A; Aa[:nt, nt] = bb
        return (expm(Aa * (b_ - a)) @ np.append(y, 1.0))[:nt]
    guardados = [(t0, y0)]                     # (tiempo, estado justo despues)
    for q in quiebres:
        tq, yq = guardados[-1]
        y = avanzar(yq, tq, q)
        for j, (_, pts) in enumerate(tablas):  # en un salto, el valor de despues
            y[m + nsen + j] = valor_tabla(pts, q)
        guardados.append((q, y))

    def sol(t):
        k = max(i for i, (tq, _) in enumerate(guardados) if tq <= t) if t >= t0 else 0
        tq, yq = guardados[k]
        y = avanzar(yq, tq, t) if t > tq else yq
        sp = np.concatenate([s0 + N @ y[:m], y[m:]])
        x, _ = resolver(sp)
        o = {}
        for k, (n, a, bb, v) in enumerate(elementos):
            u = ua(a, bb); o['V' + n] = u[:nn] @ x[:nn]; o['I' + n] = x[nn + k]
        for e, val in zip(est, sp[:ns]):
            if e[0][0] == 'C': o['Q' + e[0]] = e[3] * val
        return o
    # el estado del arranque (proyectado, o el de reposo) para quien quiera
    # arrancar igual: siru/corpus/exactas --spice se lo da a ngspice como IC
    sp0 = s0 + N @ y0[:m] if m else s0
    sol.estado0 = {e[0]: float(v) for e, v in zip(est, sp0)}
    sol.elementos = elementos
    return sol, A[:m, :m] if m else np.zeros((0, 0)), len(G)


# ---------------------------------------------------------------- circuitos
# Cada circuito: las ecuaciones (texto, va al encabezado del .exacta), el
# intervalo, los quiebres y sol(t) -> {VARIABLE: valor}. Ver EXACTAS.md.

def tanque(c, r, e, j, corte, v0, nc, nr):
    """J (fuente de corriente) carga el capacitor c, que descarga por r hacia
    la fuente e; j vale J hasta t = corte y 0 despues"""
    tau = r * c

    def v(t):
        if t <= corte:
            vinf = e + r * j
            return vinf + (v0 - vinf) * math.exp(-t / tau)
        vc = v(corte)
        return e + (vc - e) * math.exp(-(t - corte) / tau)

    def sol(t):
        x = v(t)
        return {'V' + nc: x, 'I' + nr: (x - e) / r, 'Q' + nc: c * x}
    ecs = ['%s * dV/dt = J(t) - (V - %g)/%s,  V(0) = %g,  %s = %g, %s = %g'
           % (nc, e, nr, v0, nc, c, nr, r),
           'J(t) = %g si t <= %g, 0 despues;  tau = %s*%s = %g s' % (j, corte, nr, nc, tau),
           'V = %g + (%g - %g) exp(-t/tau) hasta el corte; despues V = %g + (V(%g) - %g) exp(-(t-%g)/tau)'
           % (e + r * j, v0, e + r * j, e, corte, e, corte),
           'I%s = (V - %g)/%g' % (nr, e, r)]
    # la misma solucion por red_lineal, que ademas da las ramas (para
    # --spice); se comprueba que coincida con la formula
    el = [('J1', 'GND', '1', ('tabla', [(0, j), (corte, j), (corte, 0), (21600, 0)])),
          (nr, '1', '2', r), (nc, '1', 'GND', c), ('E3', 'GND', '2', e)]
    red, _, _ = red_lineal(el, {nc: v0})
    for x in (1000.0, corte + 500.0):
        assert abs(red(x)['V' + nc] - sol(x)['V' + nc]) < 1e-9 * abs(sol(x)['V' + nc]), 'tanque: red_lineal no da la formula'
    return dict(ecuaciones=ecs, desde=0.0, hasta=21600.0, quiebres=[corte], sol=red)


def seis_elementos():
    """E1 = 1 entre GND y 1; R12 de 1 a 2; L2 de 2 a GND; L23 de 2 a 3; C3 y
    R3 de 3 a GND; todo 1. Estados IL2, IL23, VC3"""
    f = lineal([[-1, -1, 0], [-1, -1, -1], [0, 1, -1]], [1, 1, 0], [0, 0, 0])

    def sol(t):
        il2, il23, vc3 = f(t)
        u2 = 1 - il2 - il23
        return {'IL2': il2, 'IL23': il23, 'VC3': vc3, 'QC3': vc3, 'VL2': u2, 'IR12': il2 + il23,
                'IR3': vc3, 'IC3': il23 - vc3, 'IE1': il2 + il23, 'VR12': il2 + il23,
                'VR3': vc3, 'VL23': u2 - vc3}
    ecs = ['nodo 2: (1 - U2)/R12 = IL2 + IL23  ->  U2 = 1 - IL2 - IL23',
           'L2 dIL2/dt = U2;  L23 dIL23/dt = U2 - VC3;  C3 dVC3/dt = IL23 - VC3/R3',
           'todo 1, IL2(0) = IL23(0) = VC3(0) = 0',
           "x' = A x + b, x = (IL2, IL23, VC3), A = [[-1,-1,0],[-1,-1,-1],[0,1,-1]], b = (1,1,0)",
           'x(t) = xs + expm(A t)(x0 - xs), xs = -A^-1 b = (1, 0, 0)',
           'IR12 = IE1 = IL2 + IL23;  IC3 = IL23 - VC3;  IR3 = VC3;  VL2 = U2;  QC3 = C3*VC3 = VC3']
    el = [('E1', 'GND', '1', 1), ('L2', '2', 'GND', 1), ('C3', '3', 'GND', 1), ('R3', '3', 'GND', 1),
          ('R12', '1', '2', 1), ('L23', '2', '3', 1)]
    red, _, _ = red_lineal(el, {})
    for x in (0.7, 4.0):
        assert all(abs(red(x)[k] - sol(x)[k]) < 1e-12 for k in ('IL2', 'IL23', 'VC3', 'IC3')), 'seis: red_lineal no da la formula'
    return dict(ecuaciones=ecs, desde=0.0, hasta=10.0, quiebres=[], sol=red)


def eqdif():
    # Z' = 0: Z queda en su valor inicial (0) y se resuelve el par X, Y
    g = lineal([[-6, 5], [5, -7]], [10, 0], [0, 0])

    def sol(t):
        x, y = g(t)
        return {'X': x, 'Y': y, 'Z': 0.0}
    ecs = ["X' = -6 X + 5 Y + 10;  Y' = 5 X - 7 Y + 0.2 Z;  Z' = 0;  X(0) = Y(0) = Z(0) = 0",
           'Z = 0 siempre; (X, Y) = xs + expm(A t)(x0 - xs), A = [[-6,5],[5,-7]], b = (10, 0)',
           'xs = (70/17, 50/17); autovalores -1.475 y -11.525',
           'el .ds integra con Euler a paso 0.1: su error contra esta es del orden de 0.1']
    return dict(ecuaciones=ecs, desde=0.0, hasta=1.0, quiebres=[], sol=sol)


def carga_ax():
    def sol(t):
        return {'VE1': -1.0, 'IR1': 1 / 0.9, 'CARGA': t / 0.9}
    ecs = ['E1 = 1 entre GND y 1, R1 = 0.9 de 1 a GND: IR1 = 1/0.9, VE1 = -1',
           "CARGA' = IE1 = IR1, CARGA(0) = 0  ->  CARGA = t/0.9"]
    return dict(ecuaciones=ecs, desde=0.0, hasta=1.0, quiebres=[], sol=sol)


def divisor_seno():
    def sol(t):
        s = math.sin(t)
        return {'U1': s, 'U2': s / 2, 'IR12': s / 2, 'IR4': s / 2, 'VR12': s / 2,
                'VR4': s / 2, 'IESIN1': s / 2, 'VESIN1': -s, 'DVESIN1DT': -math.cos(t)}
    ecs = ['ESIN1 = sin(T) entre GND y 1; R12 = 1 de 1 a 2; R4 = 1 de 2 a GND',
           'U1 = sin(t); IR12 = IR4 = IESIN1 = sin(t)/2; U2 = VR12 = VR4 = sin(t)/2',
           'VESIN1 = -sin(t), DVESIN1DT = -cos(t)']
    return dict(ecuaciones=ecs, desde=0.0, hasta=10.0, quiebres=[], sol=sol)


def rc_salto(como, saltos, hasta):
    """J sobre C = 1 y R = 1 en paralelo, v' = J - v, v(0) = 0; J vale 1 y 0
    alternados, empezando en 1, y cambia en cada salto (los ejemplos
    inventados/salto-*.ds)"""
    cortes = [0.0] + list(saltos) + [hasta]

    def sol(t):
        v = 0.0
        for k in range(len(cortes) - 1):
            a, b = cortes[k], cortes[k + 1]
            j = 1.0 if k % 2 == 0 else 0.0
            fin = min(t, b)
            v = j + (v - j) * math.exp(-(fin - a))
            if t <= b:
                return {'VC1': v}
        return {'VC1': v}
    ecs = ["C dV/dt = J(t) - V/R, C = R = 1, V(0) = 0; J = %s" % como,
           'J vale 1 y 0 alternados desde 1, y cambia en t = %s' % ', '.join('%g' % x for x in saltos),
           'entre saltos J es constante: V = J + (V(t0) - J) exp(-(t - t0))']
    return dict(ecuaciones=ecs, desde=0.0, hasta=float(hasta), quiebres=list(saltos), sol=sol)


NUMERO = r'[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?'


def leer_lineal(ruta):
    """(elementos, iniciales, desde, hasta, reposo) de un .d o .ds que sea un
    circuito lineal simple: solo R, L, C, E y J (con o sin el tipo de siru,
    RB-R1) de valor numerico constante, sin modelos ni llaves. Si no lo es,
    devuelve la razon (una cadena)."""
    import re
    secc, elementos, valores, ini = None, [], {}, {}
    desde, hasta, reposo, fin = 0.0, None, False, None
    for linea in open(ruta, errors='replace'):
        l = re.split(r'[;#]', linea, 1)[0].rstrip()
        if not l.strip():
            continue
        u = l.strip().upper()
        m = re.match(r'^(ELEMENTS|OUTPUTS|FUNCTIONS|INITIAL[ _]CONDITIONS|DEFINED[ _]PARAMETERS|RUN[ _]CONTROLS|'
                     r'CONTROL[ _]EQUATIONS|EQUATIONS|VARIABLES|CONSTANTS|CIRCUIT[ _]DESCRIPTION|MODEL[ _]DESCRIPTION|'
                     r'ELEMENT[ _]DESCRIPTION|END|EXECUTE|RERUN[ _]DESCRIPTION|DATA)\b', u)
        if m:
            secc = m.group(1).replace(' ', '_')
            if secc in ('MODEL_DESCRIPTION', 'RERUN_DESCRIPTION', 'EQUATIONS', 'ELEMENT_DESCRIPTION'):
                if secc != 'EQUATIONS':
                    return 'trae %s' % secc
            resto = u[m.end():].strip()
            if not resto:
                continue
            u = resto
        if secc == 'EQUATIONS':
            return 'trae ecuaciones propias (EQUATIONS)'
        if secc == 'ELEMENTS':
            if u.startswith('DERIVATIVE'):
                return 'DERIVATIVE: fuente variable'
            m = re.match(r'^(?:([A-Z]+)-)?([A-Z][A-Z0-9_]*)\s*,\s*([A-Z0-9_]+)\s*-\s*([A-Z0-9_]+)\s*=\s*(\S+)$', u)
            if not m:
                return 'rama que no se entiende: %s' % u
            tipo, nombre, a, b, v = m.groups()
            letra = (tipo or nombre)[0]
            if letra not in 'RLCEJ' or (tipo and tipo not in ('RB', 'LB', 'CB', 'EB', 'JB')):
                return 'rama %s de tipo %s' % (nombre, tipo or letra)
            if not re.match('^%s$' % NUMERO, v):
                return 'rama %s con valor %s' % (nombre, v)
            if a == b:
                return 'rama %s con los dos extremos en %s' % (nombre, a)
            elementos.append((letra + nombre[1:] if tipo else nombre, a, b, float(v)))
        elif secc in ('INITIAL_CONDITIONS', 'DEFINED_PARAMETERS'):
            for parte in u.split(','):
                m = re.match(r'^\s*([A-Z][A-Z0-9_]*)\s*=\s*(%s)\s*$' % NUMERO, parte)
                if m:
                    valores[m.group(1)] = float(m.group(2))
        elif secc == 'RUN_CONTROLS':
            m = re.match(r'^STOP\s*TIME\s*=\s*(%s)' % NUMERO, u)
            if m:
                hasta = float(m.group(1))
            m = re.match(r'^START\s*TIME\s*=\s*(%s)' % NUMERO, u)
            if m:
                desde = float(m.group(1))
            if re.match(r'^RK\s*=\s*\[\s*-1\s*\]', u):
                return 'RK=[-1]: un solo instante, sin integrar'
            if re.match(r'^RUN\s+INITIAL\s+CONDITIONS', u):
                reposo = True
            if re.match(r'^(RUN\s+(SENSITIVITY|MONTE|WORST|OPTIMI))', u):
                return u.split()[1].lower()
        elif secc == 'CONTROL_EQUATIONS':
            m = re.match(r'^T\s*<=?\s*(%s|[A-Z][A-Z0-9_]*)' % NUMERO, u)
            if m:
                fin = m.group(1)
    if not elementos:
        return 'sin ramas'
    if hasta is None and fin is not None:
        hasta = float(fin) if re.match('^%s$' % NUMERO, fin) else valores.get(fin)
    if hasta is None:
        return 'sin tiempo final'
    nombres = {e[0] for e in elementos}
    for k, v in valores.items():
        if k[0] == 'V' and k[1:] in nombres and k[1] == 'C':
            ini[k[1:]] = v
        elif k[0] == 'I' and k[1:] in nombres and k[1] == 'L':
            ini[k[1:]] = v
    return elementos, ini, desde, hasta, reposo


def de_archivo(rel):
    """un circuito de CIRCUITOS leido de su archivo con leer_lineal"""
    r = leer_lineal(os.path.join(CORPUS, rel))
    if isinstance(r, str):
        raise ValueError('%s: %s' % (rel, r))
    elementos, ini, desde, hasta, reposo = r
    sol0, A, ncort = red_lineal(elementos, ini, reposo, t0=desde)
    ecs = ['circuito lineal leido de %s y resuelto por red_lineal (siru/corpus/exactas, EXACTAS.md):' % os.path.basename(rel),
           '  ramas: ' + ', '.join('%s %s-%s = %g' % e for e in elementos),
           '  iniciales: ' + (', '.join('%s = %g' % kv for kv in sorted(ini.items())) or 'todo 0')
           + ('; arranca del reposo (RUN INITIAL CONDITIONS)' if reposo else ''),
           "  x' = A x + b en %d estados libres%s; autovalores de A: %s" % (
               A.shape[0], ' (%d restricciones de cortes o lazos)' % ncort if ncort else '',
               ', '.join('%.4g' % complex(z).real if abs(complex(z).imag) < 1e-12 else '%.4g%+.4gi' % (z.real, z.imag)
                         for z in np.linalg.eigvals(A)) if A.size else '-')]
    return dict(ecuaciones=ecs, desde=desde, hasta=hasta, quiebres=[], sol=sol0)


# Los circuitos lineales simples del corpus que se resuelven leyendo su
# archivo (siru/corpus/exactas --buscar los lista). Se agregan a mano, despues de
# mirar que la exacta tenga sentido.
AUTOMATICOS = [
    'exprtk_varios_exprtkasades2015/12elemb.d',
    'GLM/seis-elementos-bdf2-directo.d',
    'GLM/seis-elementos-bdf2.d',
    'GLM/seis-elementos-bdf2.ds',
    'GLM/seis-elementos-glm1.d',
    'GLM/seis-elementos-glm1.ds',
    'GLM/seis-elementos-rk0.ds',
    'IMEX/seis-elementos-ars222.d',
    'IMEX/seis-elementos.d',
    'ROS/seis-elementos.d',
    'Simul_corrida/cap-C-salto.d',
    'Simul_corrida/cap-CB-salto.d',
    'Simul_corrida/cap-Q.d',
    'Simul_corrida/cap-QP.d',
    'Simul_corrida/cap-igual.d',
    'Simul_corrida/chopped_circuito.d',
    'Simul_corrida/control-R.d',
    'Simul_corrida/control-S1.d',
    'Simul_corrida/dos-bobinas-salto.d',
    'Simul_corrida/dos-cap-salto.d',
    'Simul_corrida_siru/circuitoB_1.ds',
    'Simul_varios_Siru_perl_examples/circuitoBGF.ds',
    'Simul_varios_Siru_perl_examples/circuitoB_2.ds',
    'Simul_varios_Siru_perl_examples/circuitoCNB.ds',
    'Simul_varios_Siru_perl_examples/circuitoCSLO2.ds',
    'Simul_varios_Siru_perl_examples/circuitoCSLO3.ds',
    'Simul_varios_Siru_perl_examples/circuitoCSLO4.ds',
    'Simul_varios_Siru_perl_examples/circuitoCSLOCLO.ds',
    'Simul_varios_Siru_perl_examples/circuitoCSLO_1.ds',
    'Simul_varios_Siru_perl_examples/circuitoJ.ds',
    'Simul_varios_Siru_perl_examples_otros/11elem_2.d',
    'Simul_varios_Siru_perl_examples_otros/12elem_1.d',
    'Simul_varios_Siru_perl_examples_otros/6elementos_2.d',
    'Simul_varios_Siru_perl_examples_otros/circuito-S1.d',
    'Simul_varios_Siru_perl_examples_otros/sueltos.ds',
    'chopped/6elementos-condiciones.d',
    'chopped/circuito_3.d',
    'exprtk_varios_exprtk2/circuito_4.d',
    'exprtk_varios_exprtk3_tests/circuito_5.d',
    'exprtk_varios_exprtkasades2015/11elem_3.d',
    'exprtk_varios_exprtkasades2015/12elem_2.d',
    'exprtk_varios_exprtkasades2015/ab.d',
    'exprtk_varios_exprtkasades2015/abc.d',
    'exprtk_varios_exprtkasades2015/abc.ds',
    'exprtk_varios_exprtkasades2015/bob.d',
    'exprtk_varios_exprtkasades2015/circpar.d',
    'exprtk_varios_exprtkasades2015/circparmini.d',
    'exprtk_varios_exprtkasades2015/circparminib.d',
    'exprtk_varios_exprtkasades2015/circparminic.d',
    'exprtk_varios_exprtkasades2015/circuito_6.d',
    'exprtk_varios_exprtkasades2015/cond.d',
    'exprtk_varios_exprtkasades2015/xelem.d',
    'exprtk_varios_exprtkasades2015/yelem.d',
    'inventados/caja-rc.ds',
    'inventados/caja-rl.d',
    'manual_para_manual_2013_ejemplo_tests/simple.d',
    'manual_para_manual_2013_simusol_ejemplos_conmodelos/pru_pared5_1.d',
    'packages_examples_models_projects_pru_pared5_tests/pru_pared5_2.d',
    'sceptresimusol_share_sceptresimusol_examples/minimum.d',
    'slex_viejos_PMP_Siru_examples/circuitoB_4.ds',
]

def pendulo(th0, hasta):
    """c01: el pendulo PX'' = -sin(PX) que arranca en reposo en PX = th0. La
    solucion exacta con las funciones elipticas de Jacobi (scipy.special):
    PX = 2 arcsin(k sn(K - t; k)), PY = PX' = -2 k cn(K - t; k),
    k = sin(th0/2), K la integral eliptica completa (parametro m = k^2)"""
    from scipy.special import ellipj, ellipk
    k = math.sin(th0 / 2)
    mm = k * k
    K = ellipk(mm)

    def sol(t):
        sn, cn, dn, _ = ellipj(K - t, mm)
        return {'PX': 2 * math.asin(k * sn), 'PY': -2 * k * cn}
    ecs = ["PX' = PY,  PY' = -sin(PX),  PX(0) = %g,  PY(0) = 0" % th0,
           'PX = 2 arcsin(k sn(K - t; k)),  PY = -2 k cn(K - t; k),  k = sin(PX(0)/2) = %.10g' % k,
           'K = %.10g, la integral eliptica completa de primera especie (m = k^2); periodo 4K = %.10g' % (K, 4 * K),
           'solo la corrida maestra: los RERUN (otros PX(0) y otros integradores) no']
    return dict(ecuaciones=ecs, desde=0.0, hasta=hasta, quiebres=[], sol=sol)


def a04():
    """a04: el amplificador de pequena senal de dos etapas, con el modelo SS1
    (EA = .0005 VR2, R2 = 2000, R1 = .3, JB = 50 IR1) expandido a mano"""
    def ss1(t, b, e, c):
        x = 'X' + t
        return [('EA' + t, e, x, ('V', 'R2' + t, 0.0005)), ('R2' + t, c, e, 2000),
                ('R1' + t, b, x, 0.3), ('JB' + t, c, e, ('I', 'R1' + t, 50.0))]
    el = [('E1', 'GND', '2', ('sen', 0.001, 0.000628)), ('C1', '2', '3', 5e6), ('C2', '4', '5', 5e6),
          ('RB1', '3', 'GND', 100), ('RL1', 'GND', '4', 1), ('RB2', '5', 'GND', 100), ('RL2', 'GND', '6', 1)] \
        + ss1('T1', '3', 'GND', '4') + ss1('T2', '5', 'GND', '6')
    sol, A, _ = red_lineal(el, {})
    ecs = ['E1 = .001 sin(.000628 t) de GND a 2; C1 2-3 y C2 4-5 = 5e6; RB1 3-GND y RB2 5-GND = 100;',
           'RL1 GND-4 y RL2 GND-6 = 1; T1 (3-GND-4) y T2 (5-GND-6) con el modelo SS1 (B-E-C):',
           '  EA E-X = .0005 VR2, R2 C-E = 2000, R1 B-X = .3, JB C-E = 50 IR1',
           'resuelto por red_lineal (fuentes controladas y senoidal); VC1(0) = VC2(0) = 0',
           'autovalores: ' + ', '.join('%.4g' % z.real for z in np.linalg.eigvals(A))]
    return dict(ecuaciones=ecs, desde=0.0, hasta=30000.0, quiebres=[], sol=sol)


def pared():
    """slex/examples/data/Pared.d: una pared de dos capas (CLAI, CLAD) entre el
    aire de afuera, EA = TABLE 1(TIME), y la casa a 20, con la radiacion
    JRAD = Q2P(9, T2(TIME)) = 9 TABLE 2(TIME) sobre la cara de afuera"""
    t1 = [(0, 0), (10800, -1), (18000, -2), (25200, -3), (28800, -3), (32400, 0), (36000, 3.5), (39600, 7),
          (43200, 9.5), (46800, 12), (50400, 14), (54000, 15), (57600, 13), (61200, 11), (64800, 9.5),
          (68400, 8), (72000, 6.5), (75600, 5), (79200, 2.5), (82800, 1), (86400, 0)]
    t2 = [(0, 0), (28800, 0), (29160, 0), (32400, 100), (36000, 360), (39600, 580), (43200, 700), (46800, 700),
          (54000, 580), (57600, 360), (61200, 100), (61560, 0), (86400, 0)]
    el = [('EA', 'GND', 'AMB', ('tabla', t1)), ('JRAD', 'GND', 'LAI', ('tabla', [(x, 9 * y) for x, y in t2])),
          ('ECAS', 'GND', 'CAS', 20), ('RAI', 'AMB', 'LAI', 0.0185185), ('RDC', 'LAD', 'CAS', 0.0138889),
          ('RID', 'LAI', 'LAD', 0.0133333), ('CLAD', 'LAD', 'GND', 2.88e6), ('CLAI', 'LAI', 'GND', 2.88e6)]
    sol, A, _ = red_lineal(el, {'CLAI': 20, 'CLAD': 20})
    ecs = ['EA = TABLE 1(TIME) de GND a AMB; JRAD = 9 TABLE 2(TIME) de GND a LAI; ECAS = 20 de GND a CAS;',
           'RAI AMB-LAI = .0185185, RID LAI-LAD = .0133333, RDC LAD-CAS = .0138889; CLAI y CLAD = 2.88e6 a GND',
           'VCLAI(0) = VCLAD(0) = 20; resuelto por red_lineal con las dos tablas como fuentes por tramos',
           'autovalores: ' + ', '.join('%.4g' % z.real for z in np.linalg.eigvals(A))]
    return dict(ecuaciones=ecs, desde=0.0, hasta=86400.0, quiebres=[], sol=sol)


TANQUE = tanque(41860, 0.1, 5, 500, 10800, 18, 'CA', 'R5')
SEIS = seis_elementos()
CIRCUITOS = {
    'sceptresimusol_share_sceptresimusol_examples/tanque_3.d': TANQUE,
    'simusol-2017.0.2-Nilsa_bin/tanque_4.d': TANQUE,
    'simusol-gui_examples_gui_projects_tanquebgui_tests/tanquebgui.d':
        tanque(41860, 0.1, 5, 500, 14400, 18, 'C2', 'R2'),
    'Simul_perl_data/6elementos_1.d': SEIS,
    'exprtk_varios_exprtkasades2015/circuito-mano.d': SEIS,
    'exprtk_varios_exprtkasades2015/circuito-manov.d': SEIS,
    'Simul_varios_Siru_perl_examples_otros/circuito_2.ds': SEIS,
    'Simul_varios_Siru_perl_examples/eqdif.ds': eqdif(),
    'Simul_varios_Siru_perl_examples/circuitoAX.ds': carga_ax(),
    'Simul_varios_Siru_perl_examples_otros/circuito-SE2.ds': divisor_seno(),
    'Nsceptresimusol_samples/c01.d': pendulo(0.5, 20.0),
    'Nsceptresimusol_samples/a04_1.d': a04(),
    'sceptresimusol_share_sceptresimusol_examples/a04_2.d': a04(),
    'slex_examples_data/Pared.d': pared(),
    'inventados/salto-tabla.ds': rc_salto('TABLE 1(T), tabla adentro del archivo', [1], 3),
    'inventados/salto-tabla-archivo.ds': rc_salto("Table('salto-tabla-archivo.txt',...), archivo aparte", [1], 3),
    'inventados/salto-if.ds': rc_salto('if(T<1,1,0)', [1], 3),
    'inventados/salto-sgn.ds': rc_salto('0.5*(1-sgn(T-1))', [1], 3),
    'inventados/salto-mod.ds': rc_salto('if(mod(T,2)<1,1,0), onda cuadrada', [1, 2, 3, 4], 5),
    'inventados/salto-funcion.ds': rc_salto('FJ(T), FJ(x) = if(x<1,1,0)', [1], 3),
    'inventados/salto-afin.ds': rc_salto('if(4*T-2 < 2,1,0), el tiempo transformado', [1], 3),
    'inventados/salto-glm.ds': rc_salto('TABLE 1(T), integrado con BDF2 (GLM=[20])', [1], 3),
    'inventados/salto-funcion-tabla.ds': rc_salto("FX(T), FX(x) = Table('salto-tabla-archivo.txt',...,x,1)", [1], 3),
}


class _Perezoso(dict):
    """CIRCUITOS con los AUTOMATICOS leidos recien cuando se piden"""
    def __missing__(self, k):
        if k in AUTOMATICOS:
            self[k] = de_archivo(k)
            return self[k]
        raise KeyError(k)

    def __contains__(self, k):
        return dict.__contains__(self, k) or k in AUTOMATICOS

    def __iter__(self):
        return iter(list(dict.keys(self)) + [k for k in AUTOMATICOS if not dict.__contains__(self, k)])


CIRCUITOS = _Perezoso(CIRCUITOS)


def clave(archivo):
    a = os.path.abspath(archivo)
    if a.startswith(CORPUS + os.sep):
        return os.path.relpath(a, CORPUS)
    return archivo


def junto(rel, sufijo):
    """el archivo de rel con otro sufijo: NOMBRE.exacta, NOMBRE.patron, NOMBRE-Sceptre"""
    return os.path.join(CORPUS, os.path.splitext(rel)[0] + sufijo)


def leer_tabla(ruta):
    """(nombres, filas) de una tabla de texto con el tiempo primero: el
    encabezado es la ultima linea '#' antes de los datos con tantas palabras
    como columnas (la regla de tabular2)"""
    comentarios, filas = [], []
    for linea in open(ruta):
        if linea.startswith('#'):
            # un -Sceptre con RERUN trae un bloque por corrida: vale el
            # primero, el de la corrida maestra (c01: siete corridas)
            if filas and 'RERUN' in linea:
                break
            if not filas:
                comentarios.append(linea[1:].split())
            continue
        c = linea.split()
        if not c:
            continue
        try:
            filas.append([float(x) for x in c])
        except ValueError:
            continue
    n = len(filas[0]) if filas else 0
    nombres = next((c for c in reversed(comentarios) if len(c) == n), None)
    return nombres, filas


def tiempos_de(rel):
    """los tiempos en que se comparan todos: los del .patron (lo que dejo
    siru), si no los del -Sceptre (los que eligio siru/corpus/sceptre), si no el
    10, 50 y 90 % del intervalo, como siru/corpus/sceptre. Devuelve (tiempos, de donde)"""
    for suf, que in (('.patron', 'del .patron'), ('-Sceptre', 'del -Sceptre')):
        r = junto(rel, suf)
        if os.path.exists(r):
            _, filas = leer_tabla(r)
            if filas:
                return [f[0] for f in filas], que
    c = CIRCUITOS[rel]
    return [c['desde'] + (c['hasta'] - c['desde']) * q for q in (0.1, 0.5, 0.9)], 'al 10, 50 y 90 % del intervalo'


def escribir(rel):
    c = CIRCUITOS[rel]
    ts, que = tiempos_de(rel)
    nombres = list(c['sol'](ts[0]).keys())
    ruta = junto(rel, '.exacta')
    with open(ruta, 'w') as f:
        f.write('# %s: la solucion exacta en los tiempos %s, escrita por siru/corpus/exactas\n'
                '# (ver siru/corpus/EXACTAS.md). Ecuaciones:\n' % (os.path.basename(rel), que))
        for e in c['ecuaciones']:
            f.write('#   %s\n' % e)
        f.write('# TIME %s\n' % ' '.join(nombres))
        for t in ts:
            v = c['sol'](t)
            f.write('%.10g %s\n' % (t, ' '.join('%.10g' % v[n] for n in nombres)))
    return ruta


def columnas(rel, sufijo, quitar=''):
    """{VARIABLE: [valores en los tiempos]} de NOMBRE.patron o de la columna
    de SCEPTRE (o de siru, con quitar='_siru') del -Sceptre"""
    r = junto(rel, sufijo)
    if not os.path.exists(r):
        return {}
    nombres, filas = leer_tabla(r)
    if not nombres:
        return {}
    res = {}
    for k, n in enumerate(nombres[1:], 1):
        if sufijo == '-Sceptre':
            if quitar and not n.endswith(quitar):
                continue
            if not quitar and (n.endswith('_siru') or n.endswith('_dif')):
                continue
            n = n[:-len(quitar)] if quitar else n
        res[n.upper()] = [f[k] for f in filas]
    return res


def spice(rel, plazo=600):
    """NOMBRE.cir (el circuito para ngspice, desde las ramas de red_lineal)
    y NOMBRE.spice (lo que da ngspice en los tiempos comunes). Devuelve un
    mensaje de error o None."""
    import subprocess
    import tempfile
    c = CIRCUITOS[rel]
    sol = c['sol']
    el = getattr(sol, 'elementos', None)
    if not el:
        return 'sin ramas (la exacta no es de red_lineal)'
    if c['desde'] != 0:
        return 'START TIME distinto de 0'
    # n delante: un nodo 1T8 en una expresion de ngspice es 1e12*8 (T = tera)
    nodo = lambda n: '0' if n == 'GND' else 'n' + n
    nombres = {}
    for n, a, b, v in el:                      # el nombre en SPICE: la letra del tipo
        t = n[0]
        pre = {'E': 'V', 'J': 'I'}.get(t, '')
        if isinstance(v, tuple) and v[0] in ('V', 'I'):
            pre = {('E', 'V'): 'E', ('E', 'I'): 'H', ('J', 'V'): 'G', ('J', 'I'): 'F'}[(t, v[0])]
        nombres[n] = pre + n
    # una rama cuya corriente controla a otra lleva en serie una fuente de 0 V
    sensadas = {v[1] for _, _, _, v in el if isinstance(v, tuple) and v[0] == 'I'}
    lineas = ['* %s: generado por siru/corpus/exactas --spice desde las ramas de la exacta' % os.path.basename(rel),
              '* (el arranque es el de la exacta: proyectado si era desparejo)']
    e0 = sol.estado0
    for n, a, b, v in el:
        sn = nombres[n]
        if n in sensadas:
            lineas.append('VS%s %s nS%s 0' % (n, nodo(a), n))
            a = 'S' + n
        t = n[0]
        if t == 'R':
            lineas.append('%s %s %s %.15g' % (sn, nodo(a), nodo(b), v))
        elif t in 'LC':
            lineas.append('%s %s %s %.15g IC=%.15g' % (sn, nodo(a), nodo(b), v, e0.get(n, 0.0)))
        elif not isinstance(v, tuple):
            if t == 'E': lineas.append('%s %s %s DC %.15g' % (sn, nodo(b), nodo(a), v))
            else: lineas.append('%s %s %s DC %.15g' % (sn, nodo(a), nodo(b), v))
        elif v[0] == 'sen':
            f = v[2] / (2 * math.pi)
            p, m = (nodo(b), nodo(a)) if t == 'E' else (nodo(a), nodo(b))
            lineas.append('%s %s %s SIN(0 %.15g %.15g)' % (sn, p, m, v[1], f))
        elif v[0] == 'tabla':
            pts, prev = [], None
            for x, y in v[1]:
                if prev is not None and x <= prev:
                    x = prev + 1e-9 * max(1.0, abs(prev))   # PWL pide tiempos crecientes
                pts.append('%.15g %.15g' % (x, y)); prev = x
            p, m = (nodo(b), nodo(a)) if t == 'E' else (nodo(a), nodo(b))
            lineas.append('%s %s %s PWL(%s)' % (sn, p, m, ' '.join(pts)))
        else:
            r = [e for e in el if e[0] == v[1]][0]
            if v[0] == 'V':                    # k por la tension de la rama r
                p, m = (nodo(b), nodo(a)) if t == 'E' else (nodo(a), nodo(b))
                lineas.append('%s %s %s %s %s %.15g' % (sn, p, m, nodo(r[1]), nodo(r[2]), v[2]))
            else:                              # k por la corriente de la rama r
                p, m = (nodo(b), nodo(a)) if t == 'E' else (nodo(a), nodo(b))
                lineas.append('%s %s %s VS%s %.15g' % (sn, p, m, v[1], v[2]))
    vectores = []
    for n, a, b, v in el:
        if n[0] == 'C':
            vectores.append(('V' + n, 'v(%s)' % nodo(a) if b == 'GND' else
                             ('-v(%s)' % nodo(b) if a == 'GND' else 'v(%s,%s)' % (nodo(a), nodo(b)))))
        elif n[0] == 'L':
            vectores.append(('I' + n, 'i(%s)' % nombres[n]))
    if not vectores:
        return 'sin capacitores ni bobinas'
    hasta = c['hasta']
    lineas += ['.options reltol=1e-7 abstol=1e-14 vntol=1e-12 chgtol=1e-16 trtol=1',
               '.tran %.6g %.15g 0 %.6g uic' % (hasta / 2000, hasta, hasta / 2000),
               '.control', 'run', 'set wr_singlescale', 'set wr_vecnames']
    # cada vector con su let: wrdata no acepta expresiones como -v(1)
    for k, (_, x) in enumerate(vectores):
        lineas.append('let s%d = %s' % (k, x))
    lineas += ['wrdata salida.txt ' + ' '.join('s%d' % k for k in range(len(vectores))), '.endc', '.end']
    cir = os.path.join(CORPUS, os.path.splitext(rel)[0] + '.cir')
    texto = '\n'.join(lineas) + '\n'
    with tempfile.TemporaryDirectory() as d:
        open(os.path.join(d, 'c.cir'), 'w').write(texto)
        try:
            r = subprocess.run(['ngspice', '-b', 'c.cir'], cwd=d, capture_output=True, text=True, timeout=plazo)
        except subprocess.TimeoutExpired:
            return 'ngspice pasa el plazo'
        out = os.path.join(d, 'salida.txt')
        if not os.path.exists(out):
            err = [l for l in (r.stdout + r.stderr).splitlines() if 'rror' in l]
            return 'ngspice: ' + (err[0] if err else 'sin salida')
        filas = []
        for l in open(out):
            w = l.split()
            try: filas.append([float(x) for x in w])
            except ValueError: continue
    if not filas or any(len(f) != len(vectores) + 1 for f in filas):
        return 'ngspice: la salida no trae %d columnas' % (len(vectores) + 1)
    open(cir, 'w').write(texto)
    # ngspice con uic no escribe la fila de t = 0 (minimum.d empieza en
    # 0.054): se agrega con el arranque, que es el IC que se le dio
    if filas[0][0] > c['desde']:
        cero = []
        for n, _ in vectores:
            cero.append(e0.get(n[1:], 0.0))
        filas.insert(0, [c['desde']] + cero)
    t = [f[0] for f in filas]
    ts, que = tiempos_de(rel)
    conv = None
    with open(junto(rel, '.spice'), 'w') as f:
        f.write('# %s: lo que da ngspice (%s, el .cir de al lado) en los tiempos %s;\n'
                '# escrito por siru/corpus/exactas --spice (ver siru/corpus/EXACTAS.md)\n'
                % (os.path.basename(rel), 'ngspice-42', que))
        f.write('# TIME %s\n' % ' '.join(n for n, _ in vectores))
        for x in ts:
            vals = []
            for k in range(len(vectores)):
                y = [fi[k + 1] for fi in filas]
                vals.append(_interp(t, y, x))
            f.write('%.10g %s\n' % (x, ' '.join('%.12g' % v for v in vals)))
    return None


def _interp(t, y, x):
    import bisect
    if x <= t[0]: return y[0]
    if x >= t[-1]: return y[-1]
    i = bisect.bisect_right(t, x)
    return y[i - 1] + (y[i] - y[i - 1]) * (x - t[i - 1]) / (t[i] - t[i - 1])


def comparar(rel, plazo):
    """corre siru sobre una copia, lo compara con la formula en todos los
    tiempos de siru y escribe NOMBRE-Exacta con la exacta, siru, el .patron
    y SCEPTRE en los mismos tiempos. Devuelve (error, {variable: (abs, rel, t)})"""
    loader = importlib.machinery.SourceFileLoader('convergencia', os.path.join(AQUI, 'convergencia'))
    spec = importlib.util.spec_from_loader('convergencia', loader)
    conv = importlib.util.module_from_spec(spec)
    loader.exec_module(conv)
    t, vals, error = conv.correr(os.path.join(CORPUS, rel), 1, plazo)
    if error:
        return error, {}
    sol = CIRCUITOS[rel]['sol']
    exactas = [sol(x) for x in t]
    por_mayus = {k.upper(): k for k in vals}
    res = {}
    # la escala de cada variable es su maximo, pero no menos de una milesima
    # de la mayor del circuito (sin contar las cargas): una variable que es 0
    # (VL1 de 12elemb) daria una diferencia relativa de 1 por un redondeo de
    # 1e-15
    # sin las cargas: Q = C V lleva la escala de C (en el tanque, 2e6) y
    # achicaba la diferencia relativa de las corrientes
    general = max((abs(e[n]) for e in exactas for n in e if not n.startswith('Q')), default=0.0)
    for nombre in exactas[0]:
        col = por_mayus.get(nombre)
        if col is None:
            continue
        escala = max(max(abs(e[nombre]) for e in exactas), 1e-3 * general)
        peor, cuando = 0.0, t[0]
        for i, x in enumerate(t):
            d = abs(vals[col][i] - exactas[i][nombre])
            if d > peor:
                peor, cuando = d, x
        res[nombre] = (peor, peor / escala if escala > 0 else (0.0 if peor == 0 else math.inf), cuando)

    # NOMBRE-Exacta: todas las fuentes en los mismos tiempos
    ts, que = tiempos_de(rel)
    fuentes = [('siru', {n: [conv.interpolar(t, vals[por_mayus[n]], x) for x in ts]
                         for n in exactas[0] if n in por_mayus})]
    patron = columnas(rel, '.patron')
    if patron:
        fuentes.append(('patron', patron))
    sceptre = columnas(rel, '-Sceptre')
    if sceptre:
        fuentes.append(('sceptre', sceptre))
    spice_ = columnas(rel, '.spice')
    if spice_:
        fuentes.append(('spice', spice_))
    exac = {n: [sol(x)[n] for x in ts] for n in exactas[0]}
    cab, filas, peores = ['TIME'], [[x] for x in ts], []
    for n in exac:
        hay = [(f, v[n]) for f, v in fuentes if n in v and len(v[n]) == len(ts)]
        if not hay:
            continue
        cab.append(n)
        for i in range(len(ts)):
            filas[i].append(exac[n][i])
        # el mismo piso que la comparacion en toda la corrida: una milesima
        # de la mayor variable (sin las cargas); si no, una variable que
        # vale 0 daba diferencias relativas de 1e12 (IL5 de minimum en t = 0)
        escala = max(max(abs(x) for x in exac[n]), 1e-3 * general) or 1.0
        for f, v in hay:
            cab += [n + '_' + f, n + '_' + f + '_dif']
            for i in range(len(ts)):
                filas[i] += [v[i], v[i] - exac[n][i]]
            peores.append('%s_%s %.2g' % (n, f, max(abs(v[i] - exac[n][i]) for i in range(len(ts))) / escala))
    with open(junto(rel, '-Exacta'), 'w') as f:
        f.write('# %s: la solucion exacta contra siru (corrida ahora), el .patron (lo que\n'
                '# dejo siru antes), SCEPTRE (del -Sceptre) y ngspice (del .spice), los que haya,\n# en los tiempos %s;\n'
                '# escrito por siru/corpus/exactas. Por variable: la exacta y, por fuente, su valor\n'
                '# y la diferencia fuente - exacta.\n' % (os.path.basename(rel), que))
        f.write('# %s\n' % '  '.join('%22s' % c for c in cab))
        for fila in filas:
            f.write('  %s\n' % '  '.join('%22.10g' % x for x in fila))
        f.write('# mayor diferencia relativa (|dif| / max|exacta|): %s\n' % ' '.join(peores))
    return None, res


def main():
    p = argparse.ArgumentParser(description='soluciones exactas de circuitos del corpus')
    p.add_argument('archivos', nargs='*')
    p.add_argument('--escribir', action='store_true')
    p.add_argument('--comparar', action='store_true')
    p.add_argument('--tabla', action='store_true', help='con --comparar, escribe siru/corpus/EXACTAS-comparacion.md')
    p.add_argument('--plazo', type=float, default=900)
    p.add_argument('--lista', action='store_true')
    p.add_argument('--buscar', action='store_true', help='los .d y .ds del corpus que leer_lineal entiende')
    p.add_argument('--spice', action='store_true', help='NOMBRE.cir y NOMBRE.spice: el circuito corrido con ngspice')
    op = p.parse_args()

    if op.buscar:
        for r, _, fs in sorted(os.walk(CORPUS)):
            for f in sorted(fs):
                if f.endswith(('.d', '.ds')):
                    rel = os.path.relpath(os.path.join(r, f), CORPUS)
                    x = leer_lineal(os.path.join(r, f))
                    if not isinstance(x, str):
                        print('%s%s' % (rel, '' if rel not in CIRCUITOS else '  (ya esta)'))
        return 0
    rels = [clave(a) for a in op.archivos] or sorted(CIRCUITOS)
    for r in rels:
        if r not in CIRCUITOS:
            sys.exit('%s: no tiene solucion exacta en siru/corpus/exactas (ver --lista)' % r)
    if op.lista:
        for r in sorted(CIRCUITOS):
            print(r)
        return 0
    if op.spice:
        malos = 0
        for r in rels:
            e = spice(r, op.plazo)
            print('%s: %s' % (r, e or 'ngspice OK'))
            malos += bool(e) and not str(e).startswith('sin ')
        return min(malos, 125)
    if not op.comparar:
        for r in rels:
            print(escribir(r))
        return 0

    filas, malos = [], 0
    for r in rels:
        error, res = comparar(r, op.plazo)
        if error:
            malos += 1
            print('%s: no corre (%s)' % (r, error))
            filas.append((r, 'no corre: ' + error, '-', '-', '-'))
            continue
        if not res:
            malos += 1
            print('%s: siru no escribe ninguna de las variables de la exacta' % r)
            filas.append((r, 'sin variables en comun', '-', '-', '-'))
            continue
        peor = max(res, key=lambda v: res[v][1])
        a, rl, cuando = res[peor]
        detalle = ', '.join('%s %.2g' % (v, res[v][1]) for v in res)
        print('%s: peor %s, %.3g absoluta, %.2g relativa, en t=%g (%s)' % (r, peor, a, rl, cuando, detalle))
        filas.append((r, peor, '%.3g' % a, '%.2g' % rl, '%g' % cuando))
    if op.tabla:
        with open(TABLA, 'w') as t:
            t.write('# siru contra la solucion exacta\n\n'
                    'La escribe `siru/corpus/exactas --comparar --tabla`. Por circuito, la variable\n'
                    'que mas se aparta de la exacta, en los tiempos de la salida de siru:\n'
                    'diferencia absoluta y relativa (al maximo de la exacta). Las ecuaciones\n'
                    'y por que cada diferencia es la que es, en `EXACTAS.md`.\n\n'
                    '| archivo | peor variable | absoluta | relativa | en t |\n'
                    '|---|---|---|---|---|\n')
            for fila in filas:
                t.write('| `%s` | %s | %s | %s | %s |\n' % fila)
    return min(malos, 125)


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