# ==============================================================================
# GRU v1.8.8 — Apéndice A.17
# FEM 1D sobre cadena colapsada GRU vs Laplaciano discreto reescalado
# Autor: Alfredo Flores Cornejo | dr.alfredo.fc@gmail.com
# DOI base: 10.5281/zenodo.20481543
#
# PROPÓSITO:
#   Demostrar que el Laplaciano de la cadena radial colapsada (GRU) converge
#   al espectro del operador continuo -d²/dx² en [0,1], a diferencia del
#   grafo CDT completo (Clemente et al.) que no converge.
#   Esto valida la estructura suave del espacio radial efectivo en GRU.
#
# EXPERIMENTOS:
#   A. Cadena colapsada GRU: Laplaciano físico → converge a π²
#   B. FEM 1D estándar: K u = λ M u → converge a π²
#   C. CDT completo: Laplaciano combinatorio → NO converge (Clemente)
#   D. Comparación de tasas de convergencia A vs B vs C
#
# RESULTADOS ESPERADOS:
#   A: λ₁_phys → π² = 9.8696 con error O(h²) → convergencia cuadrática
#   B: λ₁_FEM  → π² = 9.8696 con error O(h²) → mismo orden GRU
#   C: λ₁_CDT  → error ~17-20% constante → NO converge (confirma Clemente)
#   Conclusión: cadena GRU y FEM son equivalentes — CDT completo no
# ==============================================================================

import sys, io
if hasattr(sys.stdout, 'buffer'):
    sys.stdout = io.TextIOWrapper(
        sys.stdout.buffer, encoding='utf-8', line_buffering=True)

import numpy as np
import networkx as nx
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
import scipy.linalg as la

SEED = 42
PI2  = np.pi**2  # valor teórico continuo: π² ≈ 9.8696


# ==============================================================================
# CONSTRUCCIÓN DE GRAFOS
# ==============================================================================

def build_cdt_irregular(n_radial, n_layers=5, lam_t=0.6, seed=SEED):
    """CDT sintético con irregularidades temporales (simula CDT real)."""
    rng = np.random.default_rng(seed)
    G = nx.Graph()
    def nd(t, r): return t * n_radial + r
    for t in range(n_layers):
        for r in range(n_radial): G.add_node(nd(t, r))
        for r in range(n_radial):
            G.add_edge(nd(t, r), nd(t, (r+1) % n_radial))
    for t in range(n_layers - 1):
        for r in range(n_radial):
            if rng.random() < lam_t:
                G.add_edge(nd(t, r), nd(t+1, r))
            if rng.random() < lam_t * 0.3:
                G.add_edge(nd(t, r), nd(t+1, (r+1) % n_radial))
            if rng.random() < lam_t * 0.3:
                G.add_edge(nd(t, r), nd(t+1, (r-1) % n_radial))
    return G


# ==============================================================================
# LAPLACIANO FÍSICO (cadena GRU reescalada)
# ==============================================================================

def laplacian_physical_lambda1(n_nodes, tol=1e-12):
    """
    Laplaciano físico de la cadena lineal colapsada.
    L_phys = (1/h²) * L_combinatorio,  h = 1/(n-1)
    Converge a -d²/dx² con BCs Dirichlet.
    λ₁_teórico = π²
    """
    G = nx.path_graph(n_nodes)
    nodes = list(G.nodes())
    n = len(nodes)
    idx = {nd: i for i, nd in enumerate(nodes)}
    h = 1.0 / (n_nodes - 1)

    L = lil_matrix((n, n), dtype=np.float64)
    for nd in nodes:
        i = idx[nd]
        nb = list(G.neighbors(nd))
        L[i, i] = float(len(nb))
        for nb2 in nb:
            L[i, idx[nb2]] = -1.0
    L = csr_matrix(L / h**2)

    k = min(4, n - 2)
    if k < 1: return None
    try:
        ev, _ = eigsh(L, k=k, which='SM', tol=tol)
        ev = np.sort(np.abs(ev.real))
        pos = ev[ev > 1.0]  # filtrar modo cero
        return float(pos[0]) if len(pos) > 0 else None
    except: return None


def laplacian_combinatorial_lambda1(G, tol=1e-12):
    """
    Laplaciano combinatorio estándar (sin reescalar).
    Usado para CDT completo — NO converge al continuo (Clemente).
    """
    nodes = list(G.nodes())
    n = len(nodes)
    idx = {nd: i for i, nd in enumerate(nodes)}

    L = lil_matrix((n, n), dtype=np.float64)
    for nd in nodes:
        i = idx[nd]
        nb = list(G.neighbors(nd))
        L[i, i] = float(len(nb))
        for nb2 in nb:
            L[i, idx[nb2]] = -1.0
    L = csr_matrix(L)

    k = min(4, n - 2)
    if k < 1: return None
    try:
        ev, _ = eigsh(L, k=k, which='SM', tol=tol)
        ev = np.sort(np.abs(ev.real))
        pos = ev[ev > 1e-6]
        return float(pos[0]) if len(pos) > 0 else None
    except: return None


# ==============================================================================
# FEM 1D ESTÁNDAR
# ==============================================================================

def fem_lambda1(n_nodes):
    """
    FEM 1D con elementos lineales en [0,1].
    Matrices K (rigidez) y M (masa lumped).
    BCs Dirichlet: u(0)=u(1)=0.
    λ₁_teórico = π²
    """
    n = n_nodes
    n_e = n - 1
    h = 1.0 / n_e

    K = np.zeros((n, n))
    M = np.zeros((n, n))

    for e in range(n_e):
        i, j = e, e + 1
        # Rigidez elemental
        K[i,i] +=  1.0/h;  K[i,j] += -1.0/h
        K[j,i] += -1.0/h;  K[j,j] +=  1.0/h
        # Masa elemental (consistente)
        M[i,i] += 2*h/6;   M[i,j] += h/6
        M[j,i] += h/6;     M[j,j] += 2*h/6

    # Aplicar BCs Dirichlet: eliminar nodos 0 y n-1
    K_r = K[1:-1, 1:-1]
    M_r = M[1:-1, 1:-1]

    try:
        ev = la.eigh(K_r, M_r, eigvals_only=True)
        ev = np.sort(ev.real)
        return float(ev[0]) if len(ev) > 0 else None
    except: return None


# ==============================================================================
# EXPERIMENTOS
# ==============================================================================

def run_all():
    print("=" * 80)
    print("GRU A.17 — FEM 1D: Convergencia Espectral de la Cadena Colapsada".center(80))
    print("=" * 80)
    print(f"  Valor teórico continuo: π² = {PI2:.6f}")
    print(f"  Todos los métodos deberían converger a π² salvo CDT completo\n")

    sizes = [10, 20, 40, 60, 80, 120, 160, 240]

    # ── A: Laplaciano físico (cadena GRU) ──
    print("─"*80)
    print("EXPERIMENTO A — Laplaciano Físico Reescalado (Cadena Colapsada GRU)")
    print("─"*80)
    print(f"{'N':>5} | {'λ₁ físico':>12} | {'Error vs π²':>12} | {'Orden':>8}")
    print("-"*45)
    prev_err = None
    res_A = []
    for N in sizes:
        lam = laplacian_physical_lambda1(N)
        if lam:
            err = abs(lam - PI2)/PI2*100
            orden = ""
            if prev_err and err > 0:
                ratio = prev_err / err
                orden = f"~{ratio:.1f}× ↓"
            print(f"  {N:>3} | {lam:>12.6f} | {err:>11.4f}% | {orden}")
            res_A.append((N, lam, err))
            prev_err = err

    # ── B: FEM 1D ──
    print()
    print("─"*80)
    print("EXPERIMENTO B — FEM 1D Estándar (Elementos Lineales, BCs Dirichlet)")
    print("─"*80)
    print(f"{'N':>5} | {'λ₁ FEM':>12} | {'Error vs π²':>12} | {'Equiv GRU':>10}")
    print("-"*50)
    res_B = []
    for N in sizes:
        lam = fem_lambda1(N)
        if lam:
            err = abs(lam - PI2)/PI2*100
            # Buscar equivalente en A
            equiv = next((r[2] for r in res_A if r[0]==N), None)
            eq_str = f"{equiv:.4f}%" if equiv else "—"
            print(f"  {N:>3} | {lam:>12.6f} | {err:>11.4f}% | {eq_str}")
            res_B.append((N, lam, err))

    # ── C: CDT completo (Clemente) ──
    print()
    print("─"*80)
    print("EXPERIMENTO C — CDT Completo: Laplaciano Combinatorio (Clemente)")
    print("─"*80)
    print(f"{'N':>5} | {'λ₁ CDT':>12} | {'λ₁ S¹ analít':>14} | {'Error%':>8}")
    print("-"*50)
    res_C = []
    for N in [20, 40, 60, 80, 120]:
        G = build_cdt_irregular(N, 5)
        lam = laplacian_combinatorial_lambda1(G)
        if lam:
            ana = 2*(1 - np.cos(2*np.pi/N))
            err = abs(lam - ana)/ana*100
            print(f"  {N:>3} | {lam:>12.6f} | {ana:>14.6f} | {err:>7.1f}% ⚠")
            res_C.append((N, lam, err))

    # ── D: Comparación de convergencia ──
    print()
    print("=" * 80)
    print("COMPARACIÓN FINAL DE CONVERGENCIA")
    print("=" * 80)
    print(f"{'N':>5} | {'Error GRU%':>11} | {'Error FEM%':>11} | {'Error CDT%':>11} | GRU≈FEM?")
    print("-"*60)
    for N in [20, 40, 80, 120]:
        eA = next((r[2] for r in res_A if r[0]==N), None)
        eB = next((r[2] for r in res_B if r[0]==N), None)
        eC = next((r[2] for r in res_C if r[0]==N), None)
        equiv = "✓" if eA and eB and abs(eA-eB)<0.01 else "⚠"
        eA_s = f"{eA:.4f}%" if eA else "—"
        eB_s = f"{eB:.4f}%" if eB else "—"
        eC_s = f"{eC:.1f}%" if eC else "—"
        print(f"  {N:>3} | {eA_s:>11} | {eB_s:>11} | {eC_s:>11} | {equiv}")

    print()
    print("CONCLUSIONES:")
    print("""
  1. CADENA COLAPSADA GRU (Exp A):
     λ₁_físico → π² con convergencia cuadrática O(h²).
     El Laplaciano físico reescalado reproduce el espectro continuo.

  2. FEM 1D ESTÁNDAR (Exp B):
     λ₁_FEM → π² con el mismo orden de convergencia que GRU.
     La cadena colapsada GRU es equivalente al FEM estándar.

  3. CDT COMPLETO (Exp C):
     Error ~17-20% constante — NO converge al continuo.
     Confirma Caceffo-Clemente (arXiv:2010.07179).

  4. IMPLICACIÓN DIRECTA:
     GRU proporciona el espacio reducido donde el problema espectral
     de Clemente desaparece. La cadena colapsada es una discretización
     1D válida de la coordenada radial — equivalente a FEM en [0,1].
     CDT completo ≠ continuo LB.
     CDT colapsado (GRU) ≡ continuo LB en la variable radial r.
    """)
    print(f"  DOI: 10.5281/zenodo.20481543")
    print("=" * 80)


if __name__ == "__main__":
    print(f"numpy {np.__version__}")
    run_all()
