"""Numerical MVP for eigenmode-targeted hidden-state actuator selection."""
import json
import numpy as np

SEED = 2944

def make_family(rng, lam):
    Q, _ = np.linalg.qr(rng.normal(size=(8, 8)))
    scales = np.array([1.0, .8, 1.25, .9, 1.1, .75, 1.3, .95])
    V = Q @ np.diag(scales)
    eigs = np.array([lam, .62, .55, .48, .40, .32, .25, .15])
    return V @ np.diag(eigs) @ np.linalg.inv(V), V

def mode_data(A):
    vals, V = np.linalg.eig(A)
    j = int(np.argmax(np.abs(vals)))
    lam = float(vals[j].real)
    v = V[:, j].real
    v /= np.linalg.norm(v)
    w = np.real(np.linalg.inv(V)[j, :])
    w /= w @ v
    return lam, v, w

def resolvent_gain(A, C, S, z=1.0):
    B = np.eye(A.shape[0])[:, S]
    R = C @ np.linalg.solve(np.eye(A.shape[0]) - z * A, B)
    return float(np.linalg.svd(R, compute_uv=False)[0])

def finite_gain(A, C, S, horizon=300):
    B = np.eye(A.shape[0])[:, S]
    P, total = np.eye(A.shape[0]), np.zeros((C.shape[0], len(S)))
    for _ in range(horizon):
        total += C @ P @ B
        P = A @ P
    return float(np.linalg.svd(total, compute_uv=False)[0])

def main():
    rng = np.random.default_rng(SEED)
    # A readout with nonzero exposure to the dominant mode.
    C = np.array([[1., -.7, .4, .2, -.3, .5, .1, -.2]])

    A, _ = make_family(rng, .97)
    lam, v, w = mode_data(A)
    eig_resid = np.linalg.norm(A @ v - lam * v)
    left_resid = np.linalg.norm(w @ A - lam * w)
    exposure = float(np.linalg.norm(C @ v))

    # Fixed eigenbasis: only lambda changes, isolating the predicted 1/(1-lambda)
    # growth of the DC resolvent.
    base_rng = np.random.default_rng(123)
    _, V = make_family(base_rng, .5)
    eigvals = np.array([.0, .62, .55, .48, .40, .32, .25, .15])
    deltas = np.array([.20, .12, .08, .05, .03, .02, .01])
    gains = []
    for d in deltas:
        D = eigvals.copy(); D[0] = 1.0 - d
        Ad = V @ np.diag(D) @ np.linalg.inv(V)
        gains.append(resolvent_gain(Ad, C, np.arange(8), z=1.0))
    slope = float(np.polyfit(np.log(1 / deltas), np.log(gains), 1)[0])

    # Compare the literal prompt score (which is inverse footprint) with the
    # modal input-leverage score |w_i| and random coordinate subsets.
    rows = []
    for seed in [11, 22, 33]:
        rr = np.random.default_rng(seed)
        At, _ = make_family(rr, .985)
        lt, vt, wt = mode_data(At)
        k = 2
        literal = np.argsort(1.0 / (1e-8 + vt**2))[-k:][::-1]
        right_rank = np.argsort(vt**2)[-k:][::-1]
        left_rank = np.argsort(wt**2)[-k:][::-1]
        sets = [rr.choice(8, k, replace=False) for _ in range(500)]
        rg = np.array([finite_gain(At, C, s) for s in sets])
        lg = finite_gain(At, C, literal)
        wg = finite_gain(At, C, right_rank)
        mg = finite_gain(At, C, left_rank)
        rows.append({
            "seed": seed, "lambda": lt,
            "literal_inverse_footprint_set": literal.tolist(),
            "right_eigenvector_top_set": right_rank.tolist(),
            "left_eigenvector_modal_leverage_set": left_rank.tolist(),
            "literal_gain": lg, "right_rank_gain": wg, "modal_left_rank_gain": mg,
            "random_mean_gain": float(rg.mean()),
            "modal_ratio_to_random_mean": float(mg / rg.mean()),
            "literal_ratio_to_random_mean": float(lg / rg.mean()),
        })
    result = {
        "seed": SEED,
        "eigenpair_residual": float(eig_resid),
        "left_eigenpair_residual": float(left_resid),
        "biorthogonal_residual": float(abs(w @ v - 1)),
        "mode_exposure": exposure,
        "resolvent_deltas": deltas.tolist(), "resolvent_gains": gains,
        "resolvent_log_slope": slope, "selection_rows": rows,
        "mean_modal_ratio": float(np.mean([r["modal_ratio_to_random_mean"] for r in rows])),
        "mean_literal_ratio": float(np.mean([r["literal_ratio_to_random_mean"] for r in rows]))
    }
    print(json.dumps(result, indent=2))
    with open("results.json", "w") as f: json.dump(result, f, indent=2)

if __name__ == "__main__": main()