1"""Numerical MVP for eigenmode-targeted hidden-state actuator selection."""
  2import json
  3import numpy as np
  4
  5SEED = 2944
  6
  7def make_family(rng, lam):
  8    Q, _ = np.linalg.qr(rng.normal(size=(8, 8)))
  9    scales = np.array([1.0, .8, 1.25, .9, 1.1, .75, 1.3, .95])
 10    V = Q @ np.diag(scales)
 11    eigs = np.array([lam, .62, .55, .48, .40, .32, .25, .15])
 12    return V @ np.diag(eigs) @ np.linalg.inv(V), V
 13
 14def mode_data(A):
 15    vals, V = np.linalg.eig(A)
 16    j = int(np.argmax(np.abs(vals)))
 17    lam = float(vals[j].real)
 18    v = V[:, j].real
 19    v /= np.linalg.norm(v)
 20    w = np.real(np.linalg.inv(V)[j, :])
 21    w /= w @ v
 22    return lam, v, w
 23
 24def resolvent_gain(A, C, S, z=1.0):
 25    B = np.eye(A.shape[0])[:, S]
 26    R = C @ np.linalg.solve(np.eye(A.shape[0]) - z * A, B)
 27    return float(np.linalg.svd(R, compute_uv=False)[0])
 28
 29def finite_gain(A, C, S, horizon=300):
 30    B = np.eye(A.shape[0])[:, S]
 31    P, total = np.eye(A.shape[0]), np.zeros((C.shape[0], len(S)))
 32    for _ in range(horizon):
 33        total += C @ P @ B
 34        P = A @ P
 35    return float(np.linalg.svd(total, compute_uv=False)[0])
 36
 37def main():
 38    rng = np.random.default_rng(SEED)
 39    # A readout with nonzero exposure to the dominant mode.
 40    C = np.array([[1., -.7, .4, .2, -.3, .5, .1, -.2]])
 41
 42    A, _ = make_family(rng, .97)
 43    lam, v, w = mode_data(A)
 44    eig_resid = np.linalg.norm(A @ v - lam * v)
 45    left_resid = np.linalg.norm(w @ A - lam * w)
 46    exposure = float(np.linalg.norm(C @ v))
 47
 48    # Fixed eigenbasis: only lambda changes, isolating the predicted 1/(1-lambda)
 49    # growth of the DC resolvent.
 50    base_rng = np.random.default_rng(123)
 51    _, V = make_family(base_rng, .5)
 52    eigvals = np.array([.0, .62, .55, .48, .40, .32, .25, .15])
 53    deltas = np.array([.20, .12, .08, .05, .03, .02, .01])
 54    gains = []
 55    for d in deltas:
 56        D = eigvals.copy(); D[0] = 1.0 - d
 57        Ad = V @ np.diag(D) @ np.linalg.inv(V)
 58        gains.append(resolvent_gain(Ad, C, np.arange(8), z=1.0))
 59    slope = float(np.polyfit(np.log(1 / deltas), np.log(gains), 1)[0])
 60
 61    # Compare the literal prompt score (which is inverse footprint) with the
 62    # modal input-leverage score |w_i| and random coordinate subsets.
 63    rows = []
 64    for seed in [11, 22, 33]:
 65        rr = np.random.default_rng(seed)
 66        At, _ = make_family(rr, .985)
 67        lt, vt, wt = mode_data(At)
 68        k = 2
 69        literal = np.argsort(1.0 / (1e-8 + vt**2))[-k:][::-1]
 70        right_rank = np.argsort(vt**2)[-k:][::-1]
 71        left_rank = np.argsort(wt**2)[-k:][::-1]
 72        sets = [rr.choice(8, k, replace=False) for _ in range(500)]
 73        rg = np.array([finite_gain(At, C, s) for s in sets])
 74        lg = finite_gain(At, C, literal)
 75        wg = finite_gain(At, C, right_rank)
 76        mg = finite_gain(At, C, left_rank)
 77        rows.append({
 78            "seed": seed, "lambda": lt,
 79            "literal_inverse_footprint_set": literal.tolist(),
 80            "right_eigenvector_top_set": right_rank.tolist(),
 81            "left_eigenvector_modal_leverage_set": left_rank.tolist(),
 82            "literal_gain": lg, "right_rank_gain": wg, "modal_left_rank_gain": mg,
 83            "random_mean_gain": float(rg.mean()),
 84            "modal_ratio_to_random_mean": float(mg / rg.mean()),
 85            "literal_ratio_to_random_mean": float(lg / rg.mean()),
 86        })
 87    result = {
 88        "seed": SEED,
 89        "eigenpair_residual": float(eig_resid),
 90        "left_eigenpair_residual": float(left_resid),
 91        "biorthogonal_residual": float(abs(w @ v - 1)),
 92        "mode_exposure": exposure,
 93        "resolvent_deltas": deltas.tolist(), "resolvent_gains": gains,
 94        "resolvent_log_slope": slope, "selection_rows": rows,
 95        "mean_modal_ratio": float(np.mean([r["modal_ratio_to_random_mean"] for r in rows])),
 96        "mean_literal_ratio": float(np.mean([r["literal_ratio_to_random_mean"] for r in rows]))
 97    }
 98    print(json.dumps(result, indent=2))
 99    with open("results.json", "w") as f: json.dump(result, f, indent=2)
100
101if __name__ == "__main__": main()