import numpy as np


def cyclic_triangle_incidence(n):
    A = np.zeros((n, n), dtype=float)
    for j in range(n):
        A[j % n, j] = 1.0
        A[(j + 1) % n, j] = 1.0
        A[(j + 2) % n, j] = 1.0
    return A


def local_operator(A, self_loop=0.0):
    up = A.sum(axis=1)
    down = A.sum(axis=0)
    P = (A / up[:, None]) @ np.diag(1.0 / down) @ A.T
    if self_loop:
        P = (1.0 - self_loop) * P + self_loop * np.eye(P.shape[0])
    return P


def teleport(P, alpha):
    n = P.shape[0]
    return (1.0 - alpha) * P + alpha * np.ones((n, n)) / n


def nonprincipal_radius(P):
    vals = np.linalg.eigvals(P)
    k = np.argmin(np.abs(vals - 1.0))
    return float(np.max(np.abs(np.delete(vals, k))))


def mfpt_to_target(P, target):
    n = P.shape[0]
    keep = [i for i in range(n) if i != target]
    Q = P[np.ix_(keep, keep)]
    t = np.linalg.solve(np.eye(n - 1) - Q, np.ones(n - 1))
    out = np.zeros(n)
    out[keep] = t
    return out


def ridge_accuracy(X, y, train, reg=1e-2):
    Xtr = X[train]
    Y = np.eye(2)[y[train]]
    W = np.linalg.solve(Xtr.T @ Xtr + reg * np.eye(X.shape[1]), Xtr.T @ Y)
    pred = np.argmax(X @ W, axis=1)
    return float(np.mean(pred == y)), pred


def run_experiment(seed=7):
    rng = np.random.default_rng(seed)
    n = 60
    A = cyclic_triangle_incidence(n)
    P0 = local_operator(A, self_loop=0.1)
    alphas = np.array([0.0, 0.01, 0.03, 0.1, 0.3, 0.6])

    # Prediction 1: uniform teleportation contracts every nonprincipal eigenvalue
    # by exactly 1-alpha on the regular complex.
    rho0 = nonprincipal_radius(P0)
    spectral = []
    for a in alphas:
        rho = nonprincipal_radius(teleport(P0, a))
        spectral.append((float(a), rho, rho / rho0, 1.0 - a))

    # Prediction 2: first-passage times decrease as global jumps become more likely.
    target = n // 2
    mfpt = []
    for a in alphas:
        t = mfpt_to_target(teleport(P0, a), target)
        mfpt.append((float(a), float(np.mean(t))))

    # Small node classification test: two noisy community signals on the cycle.
    y = np.zeros(n, dtype=int)
    y[n // 2:] = 1
    X = np.zeros((n, 8))
    X[:, 0] = 2 * y - 1
    X[:, 1] = np.cos(2 * np.pi * np.arange(n) / n)
    X[:, 2:] = rng.normal(0, 0.8, size=(n, 6))
    train = np.r_[np.arange(0, 8), np.arange(n // 2, n // 2 + 8)]
    # Repeated propagation isolates the oversmoothing/mixing tradeoff.
    cls = []
    for a in alphas:
        P = teleport(P0, a)
        H = X.copy()
        for _ in range(8):
            H = P @ H
        acc, _ = ridge_accuracy(H, y, train)
        cls.append((float(a), acc, float(np.std(H[:, 0]))))

    return {"n": n, "rho0": rho0, "spectral": spectral, "mfpt": mfpt,
            "classification": cls}


if __name__ == '__main__':
    import json
    print(json.dumps(run_experiment(), indent=2))