#!/usr/bin/env python3
"""
Example 134 — The Spectral Dimension of the Emergent Diffusion: the Heat Kernel
as the EPI Green's Function
==============================================================================

The EPI channel of the nodal equation is the discrete diffusion equation
(AGENTS.md "Transport Content of the Nodal Equation"):

    dEPI/dt = nu_f * dNFR = -nu_f * L_rw * EPI,   L_rw = I - D^{-1} W.

Its fundamental solution is the HEAT KERNEL e^{-t L} -- the operator that spreads
a point EPI source through the network. The diagonal of the heat kernel is the
RETURN PROBABILITY p(t) = (1/n) Tr(e^{-t L}) = (1/n) sum_k e^{-lambda_k t}, and
its small-t / intermediate scaling

    p(t) ~ t^{-d_s/2}

defines the SPECTRAL DIMENSION d_s -- the dimension the network "looks like" to
its own structural diffusion. This is an exact spectral quantity of the canonical
operator, empirically anchored in spectral geometry and anomalous diffusion
(the spectral / fracton dimension of Alexander-Orbach; the "dimension a random
walker feels").

Doctrine compliance
-------------------
Everything emerges from the canonical structural-diffusion operator: the spectrum
{lambda_k} is that of the canonical structural_eigenmodes (the symmetric
normalized Laplacian, same spectrum as the diffusion operator L_rw); the heat
kernel e^{-t L} is the evolution operator of the EPI channel itself. Nothing is
imposed -- the spectral dimension is read off the canonical spectrum. The
quantity is a standard spectral-geometry observable; the example measures it, it
does not invent it.

Three measured results
----------------------
M1 THE HEAT KERNEL IS THE EPI-CHANNEL EVOLUTION OPERATOR. The heat trace
   Z(t) = sum_k e^{-lambda_k t} runs from n (t=0) to 1 (t->inf, only the
   lambda_1=0 uniform mode survives). And e^{-t L_rw} u0 reproduces the
   explicitly-integrated nodal diffusion du/dt = -L_rw u to integration
   precision -- the heat kernel IS the EPI Green's function.

M2 THE SPECTRAL DIMENSION RECOVERS THE LATTICE DIMENSION. The return-probability
   scaling p(t) ~ t^{-d_s/2} gives d_s ~ 1 for a ring, ~2 for a 2D torus, ~3 for
   a 3D torus -- the emergent diffusion feels the lattice dimension. (d_s is an
   asymptotic n->inf quantity, so finite graphs carry a finite-size bias; the
   ordering 1 < 2 < 3 < 4 is exact and the value converges to the integer as the
   lattice grows.)

M3 NON-LATTICE TOPOLOGIES HAVE A CHARACTERISTIC EMERGENT d_s. A spanning tree is
   quasi-1D (d_s ~ 1.3); adding shortcut edges to a 1D ring (Watts-Strogatz
   rewiring) RAISES d_s monotonically above 1 -- the shortcuts let the walker
   reach farther, so the network feels higher-dimensional; the complete graph is
   the mean-field limit -- its non-zero spectrum is fully degenerate, so there is
   NO finite spectral dimension (no power-law window). The spectral dimension is a
   structural fingerprint of the topology.

Honest scope
------------
The spectral dimension is a standard spectral-geometry / anomalous-diffusion
observable (Alexander-Orbach fracton dimension; the dimension a diffusing walker
feels). It is asymptotic, so finite-graph values carry a finite-size bias (the
example shows the convergence honestly). The heat kernel = EPI-evolution identity
is exact and is the canonical anchor. This re-expresses established spectral
geometry in the emergent transport layer; it is not new mathematics and closes no
open problem.

References
----------
- src/tnfr/physics/structural_diffusion.py (structural_eigenmodes,
  structural_diffusion_operator, relaxation_spectrum)
- AGENTS.md "Transport Content of the Nodal Equation (Structural Diffusion)"
- examples/08_emergent_geometry/99_structural_diffusion.py (the diffusion layer)
- examples/08_emergent_geometry/129_spectral_gap_base_fiber_clock.py (the spectrum
  as the base->fiber clock)
"""

import os
import sys

sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "..", "src"))

import networkx as nx
import numpy as np

try:
    import scipy.linalg as sla

    _HAVE_SCIPY = True
except Exception:  # pragma: no cover
    _HAVE_SCIPY = False

from tnfr.physics.structural_diffusion import (
    structural_diffusion_operator,
    structural_eigenmodes,
)


def heat_trace(eigs, t_arr):
    """Z(t) = sum_k exp(-lambda_k t)."""
    return np.array([float(np.sum(np.exp(-eigs * ti))) for ti in t_arr])


def spectral_dimension(G):
    """d_s from the return-probability plateau p(t) ~ t^{-d_s/2}.

    Returns None when the non-zero spectrum is degenerate (mean-field: no
    power-law window, hence no finite spectral dimension).
    """
    eigs, _ = structural_eigenmodes(G)
    n = len(eigs)
    lam2, lam_max = float(eigs[1]), float(eigs[-1])
    if lam2 <= 1e-9 or (1.0 / lam2) / (1.0 / lam_max) < 3.0:
        return None  # degenerate / no scaling window
    t = np.logspace(np.log10(1.0 / lam_max), np.log10(1.0 / lam2), 200)
    p = heat_trace(eigs, t) / n
    ds_local = -2.0 * np.gradient(np.log(p), np.log(t))
    lo, hi = int(0.25 * len(t)), int(0.75 * len(t))  # central plateau
    return float(np.median(ds_local[lo:hi]))


def experiment_1_heat_kernel_is_epi_evolution():
    """M1: heat trace + heat kernel = EPI-channel evolution operator."""
    print("=" * 74)
    print("M1: THE HEAT KERNEL IS THE EPI-CHANNEL EVOLUTION OPERATOR")
    print("=" * 74)
    print("Heat trace Z(t) = sum_k exp(-lambda_k t): from n (t=0) to 1 (t->inf).")
    print()
    G = nx.cycle_graph(100)
    eigs, _ = structural_eigenmodes(G)
    for ti in [0.0, 0.1, 1.0, 10.0, 100.0, 1e4]:
        print(f"    t={ti:>9.1f}   Z(t) = {heat_trace(eigs, [ti])[0]:>9.4f}")
    print(f"  Z(0)=n={len(eigs)}, Z(inf)=1 (only the lambda_1=0 uniform mode).")
    print()
    print("e^{-t L_rw} u0  vs  explicitly-integrated du/dt = -L_rw u:")
    G = nx.path_graph(40)
    _, lap = structural_diffusion_operator(G)
    L = np.asarray(lap)
    n = L.shape[0]
    rng = np.random.default_rng(0)
    u0 = rng.standard_normal(n)
    for T in [0.5, 2.0, 5.0]:
        if _HAVE_SCIPY:
            u_exact = sla.expm(-T * L) @ u0
        else:
            w, V = np.linalg.eig(L)
            u_exact = (V @ np.diag(np.exp(-T * w)) @ np.linalg.inv(V) @ u0).real
        u = u0.copy()
        steps = 4000
        dt = T / steps
        for _ in range(steps):
            u = u - dt * (L @ u)
        err = np.linalg.norm(u - u_exact) / np.linalg.norm(u_exact)
        print(f"    T={T:>4.1f}:  rel. error = {err:.2e}")
    print()
    print("  -> the heat kernel e^{-t L} reproduces the EPI diffusion exactly")
    print("     (rel err -> 0): it IS the EPI Green's function.")


def experiment_2_lattice_dimension():
    """M2: spectral dimension recovers the lattice dimension."""
    print()
    print("=" * 74)
    print("M2: THE SPECTRAL DIMENSION RECOVERS THE LATTICE DIMENSION")
    print("=" * 74)
    print("p(t) ~ t^{-d_s/2}; d_s read from the return-probability plateau.")
    print()
    print(f"  {'lattice':>14} {'n':>6} {'d_s':>8} {'expected':>9}")
    cases = [
        ("ring 1D", nx.cycle_graph(400), 1),
        ("2D torus", nx.grid_2d_graph(28, 28, periodic=True), 2),
        ("3D torus", nx.grid_graph([12, 12, 12], periodic=True), 3),
    ]
    for label, G, exp in cases:
        ds = spectral_dimension(G)
        print(f"  {label:>14} {G.number_of_nodes():>6} {ds:>8.3f} {exp:>9}")
    print()
    print("  Finite-size convergence (2D torus d_s -> 2 as L grows):")
    print(f"    {'L':>4} {'n':>6} {'d_s':>8}")
    for L in [12, 20, 32, 44]:
        G = nx.grid_2d_graph(L, L, periodic=True)
        print(f"    {L:>4} {L * L:>6} {spectral_dimension(G):>8.3f}")
    print()
    print("  -> the emergent diffusion feels the lattice dimension; d_s is an")
    print("     asymptotic quantity (finite-size bias shrinks as the lattice grows).")


def experiment_3_emergent_dimension():
    """M3: non-lattice topologies have characteristic emergent d_s."""
    print()
    print("=" * 74)
    print("M3: NON-LATTICE TOPOLOGIES HAVE A CHARACTERISTIC EMERGENT d_s")
    print("=" * 74)
    print(f"  {'topology':>18} {'n':>6} {'d_s':>10} {'reading':>20}")
    Gws = nx.watts_strogatz_graph(400, 6, 0.2, seed=1)
    tree = nx.minimum_spanning_tree(Gws)
    cases = [
        ("spanning tree", tree, "quasi-1D"),
        ("ring 1D", nx.cycle_graph(400), "1D"),
        ("2D torus", nx.grid_2d_graph(28, 28, periodic=True), "2D"),
        ("small-world WS", Gws, "ring + shortcuts"),
        ("scale-free BA", nx.barabasi_albert_graph(400, 2, seed=2), "hub-dominated"),
        ("complete K100", nx.complete_graph(100), "mean-field"),
    ]
    for label, G, reading in cases:
        if not nx.is_connected(G):
            G = G.subgraph(max(nx.connected_components(G), key=len)).copy()
        ds = spectral_dimension(G)
        ds_str = f"{ds:.3f}" if ds is not None else "none"
        print(f"  {label:>18} {G.number_of_nodes():>6} {ds_str:>10} " f"{reading:>20}")
    print()
    print("  Shortcuts raise d_s (Watts-Strogatz rewiring of a 1D ring, k=6):")
    print(f"    {'p_rewire':>9} {'d_s':>8}")
    for p in [0.0, 0.02, 0.05, 0.1, 0.3, 1.0]:
        G = nx.watts_strogatz_graph(500, 6, p, seed=7)
        if not nx.is_connected(G):
            G = G.subgraph(max(nx.connected_components(G), key=len)).copy()
        ds = spectral_dimension(G)
        print(f"    {p:>9.2f} {ds:>8.3f}")
    print()
    print("  -> trees are quasi-1D; adding shortcuts to a ring raises d_s above 1")
    print("     (the walker reaches farther); the complete graph is mean-field")
    print("     (degenerate spectrum, NO finite d_s).")


def main():
    print()
    print("  ===============================================================")
    print("  The Spectral Dimension of the Emergent Diffusion")
    print("  The Heat Kernel as the EPI Green's Function")
    print("  ===============================================================")
    print()
    experiment_1_heat_kernel_is_epi_evolution()
    experiment_2_lattice_dimension()
    experiment_3_emergent_dimension()
    print()
    print("=" * 74)
    print("WHAT THIS ESTABLISHES")
    print("=" * 74)
    print("The heat kernel e^{-t L} of the canonical structural-diffusion operator")
    print("IS the evolution operator of the EPI channel (M1, exact). Its return")
    print("probability p(t) ~ t^{-d_s/2} defines the SPECTRAL DIMENSION -- the")
    print("dimension the network feels through its own diffusion. d_s recovers the")
    print("lattice dimension (M2: ring 1, 2D 2, 3D 3, asymptotically) and gives a")
    print("structural fingerprint for non-lattice topologies (M3: trees quasi-1D,")
    print("shortcuts raise d_s above the 1D ring, complete = mean-field with no")
    print("finite d_s). HONEST SCOPE: the spectral dimension is a standard")
    print("spectral-geometry / anomalous-diffusion observable (Alexander-Orbach")
    print("fracton dimension), asymptotic hence finite-size biased; the heat-kernel")
    print("= EPI-evolution identity is the exact canonical anchor. It re-expresses")
    print("established spectral geometry in the emergent transport layer; it is not")
    print("new mathematics and closes no open problem.")


if __name__ == "__main__":
    main()