1"""Calibrated compact-support anomaly score and numerical verification."""
  2import json
  3from pathlib import Path
  4import numpy as np
  5from scipy.special import gammaln
  6from scipy.integrate import quad
  7from sklearn.metrics import roc_auc_score
  8
  9
 10def parameters(d, R2):
 11    if R2 <= d + 2:
 12        raise ValueError("R2 must exceed d+2")
 13    gamma = 2.0 / (R2 - d - 2.0)
 14    a = 1.0 / gamma
 15    logC = (gammaln(d / 2 + 1 + a) - d / 2 * np.log(np.pi)
 16            - d / 2 * np.log(R2) - gammaln(1 + a))
 17    return gamma, a, logC
 18
 19
 20def compact_score(X, mu, cov, R2, barrier=False):
 21    """Negative log density; outside support is +inf (or finite training barrier)."""
 22    X = np.asarray(X)
 23    L = np.linalg.cholesky(cov)
 24    v = np.linalg.solve(L, (X - mu).T).T
 25    r2 = np.sum(v * v, axis=1)
 26    gamma, a, logC = parameters(X.shape[1], R2)
 27    u = 1.0 - r2 / R2
 28    if barrier:
 29        # Stable surrogate specified in the idea for optimization.
 30        return -logC - a * np.log(np.maximum(u, 1e-6)) + np.logaddexp(0., r2 - R2), r2
 31    score = np.full(len(X), np.inf)
 32    inside = r2 < R2
 33    score[inside] = -logC - a * np.log1p(-r2[inside] / R2)
 34    return score, r2
 35
 36
 37def gaussian_score(X, mu, cov):
 38    L = np.linalg.cholesky(cov)
 39    v = np.linalg.solve(L, (X - mu).T).T
 40    return np.sum(v * v, axis=1)
 41
 42
 43def sample_compact(rng, n, d, R2):
 44    """Exact sampler: r2/R2 ~ Beta(d/2, 1+1/gamma), direction uniform."""
 45    _, a, _ = parameters(d, R2)
 46    y = rng.beta(d / 2, a + 1, size=n)
 47    direction = rng.normal(size=(n, d))
 48    direction /= np.linalg.norm(direction, axis=1, keepdims=True)
 49    return direction * np.sqrt(R2 * y)[:, None]
 50
 51
 52def verify(d=5, seed=1013):
 53    rng = np.random.default_rng(seed)
 54    radii = [d + 2.5, d + 6, d + 20, d + 80]
 55    rows = []
 56    # Predictions: direct normalization is 1 and calibrated E[r^2] is d.
 57    for R2 in radii:
 58        gamma, a, logC = parameters(d, R2)
 59        z = sample_compact(rng, 250000, d, R2)
 60        r2 = np.sum(z*z, axis=1)
 61        surface = 2 * np.pi ** (d / 2) / np.exp(gammaln(d / 2))
 62        integral = quad(lambda t: surface * t ** (d - 1) * np.exp(logC)
 63                        * (1 - t * t / R2) ** a, 0, np.sqrt(R2),
 64                        epsabs=1e-10)[0]
 65        rows.append({"R2": R2, "gamma": gamma, "a": a,
 66                     "mean_r2": float(r2.mean()), "target_d": d,
 67                     "max_r2": float(r2.max()),
 68                     "normalization": float(integral),
 69                     "predicted_boundary_slope": -a})
 70    # Prediction: score near boundary has slope -1/gamma against log(epsilon).
 71    boundary_sweep = []
 72    eps = np.logspace(-2, -10, 9)
 73    for R2 in radii:
 74        _, a, logC = parameters(d, R2)
 75        x = np.sqrt(R2 * (1 - eps))[:, None] * np.array([[1.] + [0.] * (d - 1)])
 76        scores, _ = compact_score(x, np.zeros(d), np.eye(d), R2)
 77        observed = np.polyfit(np.log(eps), scores + logC, 1)[0]
 78        boundary_sweep.append({"R2": R2, "predicted_slope": -a,
 79                               "observed_slope": float(observed),
 80                               "abs_error": float(abs(observed + a))})
 81    R2 = d + 10
 82    outside = compact_score(np.array([[np.sqrt(R2) * 1.001] + [0.]*(d-1)]),
 83                            np.zeros(d), np.eye(d), R2)[0][0]
 84    rows_boundary = {"sweep": boundary_sweep,
 85                     "outside_is_infinite": bool(np.isinf(outside))}
 86    # Affine geometry prediction.
 87    A = np.array([[1.5, .2, 0, 0, 0], [.1, .8, .1, 0, 0],
 88                  [0, .1, 1.2, .1, 0], [0, 0, .1, .9, .2],
 89                  [0, 0, 0, .1, 1.1]])
 90    z = sample_compact(rng, 2000, d, d + 10)
 91    cov = A @ A.T
 92    _, r_orig = compact_score(z, np.zeros(d), np.eye(d), d + 10)
 93    _, r_aff = compact_score(z @ A.T, np.zeros(d), cov, d + 10)
 94    affine_err = float(np.max(np.abs(r_orig-r_aff)))
 95    # Secondary equal-feature OOD comparison.
 96    nominal = sample_compact(rng, 4000, d, d + 10)
 97    anomalies = rng.normal(size=(4000, d)) + np.array([2.2, 0, 0, 0, 0])
 98    X = np.vstack([nominal, anomalies])
 99    y = np.r_[np.zeros(len(nominal)), np.ones(len(anomalies))]
100    gs = gaussian_score(X, np.zeros(d), np.eye(d))
101    cs, _ = compact_score(X, np.zeros(d), np.eye(d), d + 10)
102    cs_rank = np.where(np.isinf(cs), np.max(np.where(np.isfinite(cs), cs, 0)) + 1e6, cs)
103    auc_g = roc_auc_score(y, gs); auc_c = roc_auc_score(y, cs_rank)
104    qg = np.quantile(gs[:len(nominal)], .95); qc = np.quantile(cs_rank[:len(nominal)], .95)
105    fpr_g = np.mean(gs[len(nominal):] <= qg); fpr_c = np.mean(cs_rank[len(nominal):] <= qc)
106    return {"calibration_rows": rows, "boundary": rows_boundary,
107            "affine_max_abs_radius_error": affine_err,
108            "ood": {"gaussian_auc": float(auc_g), "compact_auc": float(auc_c),
109                    "gaussian_anomaly_accept_rate_at_95pct_nominal": float(fpr_g),
110                    "compact_anomaly_accept_rate_at_95pct_nominal": float(fpr_c)}}
111
112if __name__ == "__main__":
113    out = verify()
114    Path("results.json").write_text(json.dumps(out, indent=2))
115    print(json.dumps(out, indent=2))