"""Forward-intersection spectral latent dynamics: numerical verification.

Toy: G=span(e0,...,e3) in R^5.  e0/e1/e2 are supported stable modes;
e3 is an unstable nuisance mode whose image leaks into e4.  We compare
exact nullspace intersections with a finite-tolerance principal-angle filter.
"""
import json
import numpy as np


def orth(A, tol=1e-10):
    if A.size == 0:
        return np.zeros((A.shape[0], 0))
    u, s, _ = np.linalg.svd(A, full_matrices=False)
    scale = s[0] if len(s) else 1.0
    return u[:, :int(np.sum(s > tol * scale))]


def intersection_basis(Q, R, tol=1e-8):
    """Basis for range(Q) intersect range(R) from ker([Q,-R])."""
    p, q = Q.shape[1], R.shape[1]
    _, s, vh = np.linalg.svd(np.concatenate([Q, -R], axis=1), full_matrices=True)
    scale = s[0] if len(s) else 1.0
    rank = int(np.sum(s > tol * scale))
    return orth(Q @ vh[rank:, :].T[:p, :], tol)


def principal_filter(Q, R, tau):
    """Approximate intersection: retain Q directions with sin(angle)<=tau."""
    # Principal angles require orthonormal bases on both sides.
    Rorth = orth(R)
    U, c, _ = np.linalg.svd(Q.T @ Rorth, full_matrices=False)
    keep = c >= np.sqrt(max(0., 1. - tau * tau))
    return Q @ U[:, keep], c


def make_K(alpha):
    K = np.zeros((5, 5))
    K[0, 0], K[1, 1], K[2, 2] = .90, .70, .80
    K[3, 3] = 1.30                 # unsupported unstable mode
    K[4, 3] = alpha                # forward incompatibility/leakage
    K[4, 4] = 1.0
    return K


def run(seed=7):
    np.random.default_rng(seed)       # fixed-seed API; this toy is deterministic
    Q0 = np.eye(5)[:, :4]
    tau = .20
    predicted_boundary = 1.30 * tau / np.sqrt(1 - tau * tau)
    alphas = np.array([0., .10, .20, .25, .263, .27, .35, .60, 1.0])
    rows = []
    for alpha in alphas:
        K = make_K(alpha)
        Q1, cosines = principal_filter(Q0, K @ Q0, tau)
        Af, A0 = Q1.T @ K @ Q1, Q0.T @ K @ Q0
        ef, e0 = np.linalg.eigvals(Af), np.linalg.eigvals(A0)
        supported = [.90, .70, .80]
        support_err = max(min(abs(z - x) for z in ef) for x in supported)
        sin_angle = alpha / np.sqrt(1.30**2 + alpha**2)
        rows.append({
            "alpha": float(alpha),
            "predicted_nuisance_retained": bool(alpha <= predicted_boundary),
            "observed_nuisance_retained": bool(np.any(abs(ef - 1.30) < 1e-6)),
            "predicted_dim": 4 if alpha <= predicted_boundary else 3,
            "observed_dim": int(Q1.shape[1]),
            "sin_principal_angle": float(sin_angle),
            "spectral_radius_baseline": float(max(abs(e0))),
            "spectral_radius_filtered": float(max(abs(ef)) if len(ef) else 0.),
            "supported_eigenvalue_error": float(support_err),
            "cosines": [float(x) for x in cosines],
        })

    exact_dims = {str(a): int(intersection_basis(Q0, make_K(a) @ Q0).shape[1])
                  for a in [0., .1, .5]}

    # Equal-step rollout from a nuisance-contaminated latent state.
    alpha = .60
    K = make_K(alpha)
    Q1, _ = principal_filter(Q0, K @ Q0, tau)
    A0, Af = Q0.T @ K @ Q0, Q1.T @ K @ Q1
    z0 = np.array([1., -1., .5, 1.])
    base, filt = z0.copy(), Q1.T @ (Q0 @ z0)
    for _ in range(30):
        base, filt = A0 @ base, Af @ filt
    rollout = {"steps": 30, "baseline_norm": float(np.linalg.norm(base)),
               "filtered_norm": float(np.linalg.norm(filt)),
               "norm_ratio_filtered_over_baseline": float(np.linalg.norm(filt) / np.linalg.norm(base))}

    observed = {
        "P1_boundary_observed_transition": "see rows; dimension/nuisance switch",
        "P2_max_supported_error": max(r["supported_eigenvalue_error"] for r in rows),
        "P3_exact_dimensions": exact_dims, "rollout": rollout}
    passed = (observed["P2_max_supported_error"] < 1e-10
              and exact_dims == {"0.0": 4, "0.1": 3, "0.5": 3}
              and all(r["observed_dim"] == r["predicted_dim"] for r in rows))
    return {"tau": tau,
            "predicted": {"P1_boundary_alpha": predicted_boundary,
                           "P2_supported_eigenpairs_preserved_error": "0 (floating point)",
                           "P3_exact_intersection_dimension": "4 at alpha=0; 3 for any nonzero alpha"},
            "observed": observed, "rows": rows, "pass": bool(passed)}


if __name__ == "__main__":
    print(json.dumps(run(), indent=2))