1import json
  2import numpy as np
  3
  4
  5def steering(points, omega, kappa):
  6    points = np.asarray(points)
  7    return np.exp(1j * kappa * (omega @ points.T)) / np.sqrt(len(omega))
  8
  9
 10def orthonormal_columns(A):
 11    Q, _ = np.linalg.qr(A)
 12    return Q[:, :A.shape[1]]
 13
 14
 15def adversarial_rotation(U, eps, rng):
 16    """Rotate every signal direction toward an orthogonal complement by angle asin(eps)."""
 17    M, s = U.shape
 18    R = rng.standard_normal((M, s)) + 1j * rng.standard_normal((M, s))
 19    R = R - U @ (U.conj().T @ R)
 20    W = orthonormal_columns(R)
 21    theta = np.arcsin(eps)
 22    V = np.cos(theta) * U + np.sin(theta) * W
 23    V = orthonormal_columns(V)
 24    return V
 25
 26
 27def scores(U, phi):
 28    corr = U.conj().T @ phi
 29    power = np.sum(np.abs(corr) ** 2, axis=0)
 30    residual = 1.0 - power
 31    # A simple cosine-style control: retain only the strongest learned direction.
 32    cosine = 1.0 - np.max(np.abs(corr) ** 2, axis=0)
 33    return residual.real, cosine.real
 34
 35
 36def nearest_local_minimum(values, grid, target, radius=0.18):
 37    d = np.linalg.norm(grid - target[None, :], axis=1)
 38    eligible = d <= radius
 39    idx = np.where(eligible)[0][np.argmin(values[eligible])]
 40    return grid[idx], float(values[idx])
 41
 42
 43def main():
 44    rng = np.random.default_rng(3137)
 45    M, s, kappa = 64, 3, 18.0
 46    omega = rng.normal(size=(M, 2))
 47    omega /= np.linalg.norm(omega, axis=1, keepdims=True)
 48    targets = np.array([[-0.42, -0.22], [0.08, 0.31], [0.43, -0.34]])
 49    U = orthonormal_columns(steering(targets, omega, kappa))
 50    P = U @ U.conj().T
 51
 52    # Cheap direct verification of |q_tilde-q| <= ||P_tilde-P||.
 53    check_rows = []
 54    for eps in [0.0, 0.02, 0.05, 0.10, 0.20, 0.35]:
 55        V = adversarial_rotation(U, eps, rng)
 56        Pt = V @ V.conj().T
 57        projector_distance = np.linalg.norm(Pt - P, 2)
 58        test_points = rng.uniform(-0.8, 0.8, size=(300, 2))
 59        ph = steering(test_points, omega, kappa)
 60        q, _ = scores(U, ph)
 61        qt, _ = scores(V, ph)
 62        max_delta = float(np.max(np.abs(qt - q)))
 63        check_rows.append({"requested_eps": eps, "projector_distance": projector_distance,
 64                           "max_score_change": max_delta,
 65                           "bound_slack": projector_distance - max_delta})
 66
 67    # Localization on a common dense candidate grid, using all three known wells.
 68    axis = np.linspace(-0.8, 0.8, 81)
 69    xx, yy = np.meshgrid(axis, axis)
 70    grid = np.stack([xx.ravel(), yy.ravel()], axis=1)
 71    grid_phi = steering(grid, omega, kappa)
 72    eps = 0.20
 73    V = adversarial_rotation(U, eps, rng)
 74    q_idea, c_idea = scores(V, grid_phi)
 75    # Baseline uses the same estimated subspace but a max-cosine similarity rather than
 76    # accumulated projection residual, which is a common single-component control.
 77    idea_err, base_err = [], []
 78    idea_vals, base_vals = [], []
 79    for t in targets:
 80        pi, vi = nearest_local_minimum(q_idea, grid, t)
 81        pb, vb = nearest_local_minimum(c_idea, grid, t)
 82        idea_err.append(float(np.linalg.norm(pi - t)))
 83        base_err.append(float(np.linalg.norm(pb - t)))
 84        idea_vals.append(vi)
 85        base_vals.append(vb)
 86
 87    # Repeat rotations to assess graceful degradation, keeping the same Fourier setup.
 88    degradation = []
 89    for eps2 in [0.0, 0.05, 0.10, 0.20, 0.35, 0.50]:
 90        V2 = adversarial_rotation(U, eps2, rng)
 91        q2, c2 = scores(V2, grid_phi)
 92        ei, eb = [], []
 93        for t in targets:
 94            pi, _ = nearest_local_minimum(q2, grid, t)
 95            pb, _ = nearest_local_minimum(c2, grid, t)
 96            ei.append(np.linalg.norm(pi - t)); eb.append(np.linalg.norm(pb - t))
 97        degradation.append({"eps": eps2, "projector_distance": float(np.linalg.norm(V2@V2.conj().T-P,2)),
 98                            "idea_mean_localization_error": float(np.mean(ei)),
 99                            "cosine_mean_localization_error": float(np.mean(eb))})
100
101    out = {"config": {"M": M, "rank": s, "kappa": kappa, "targets": targets.tolist()},
102           "math_check": check_rows,
103           "localization_eps_0p2": {"idea_mean_error": float(np.mean(idea_err)),
104                                     "cosine_mean_error": float(np.mean(base_err)),
105                                     "idea_errors": idea_err, "cosine_errors": base_err,
106                                     "idea_scores": idea_vals, "cosine_scores": base_vals},
107           "degradation": degradation}
108    print(json.dumps(out, indent=2))
109
110
111if __name__ == "__main__":
112    main()