1"""Small reachable-set risk head for Gaussian linearized dynamics."""
 2import numpy as np
 3from scipy.special import ndtr
 4
 5
 6def propagate(mu, cov, jacobian, process_cov):
 7    """One Gaussian covariance recursion step."""
 8    mu = np.asarray(mu, dtype=float)
 9    cov = np.asarray(cov, dtype=float)
10    J = np.asarray(jacobian, dtype=float)
11    Q = np.asarray(process_cov, dtype=float)
12    return mu, J @ cov @ J.T + Q
13
14
15def halfspace_probability(mu, cov, c, b, eps=1e-12):
16    """P(c^T X >= b), X ~ N(mu,cov)."""
17    c = np.asarray(c, dtype=float)
18    variance = max(float(c @ cov @ c), eps)
19    margin = (float(b - c @ mu)) / np.sqrt(variance)
20    return float(ndtr(-margin)), margin
21
22
23def reachable_risk(mu0, cov0, dynamics, controls, unsafe_halfspaces, process_cov,
24                   return_trace=False):
25    """Propagate a locally linear Gaussian model and union risk over H steps.
26
27    dynamics(mu, u) returns (next_mean, jacobian). Halfspaces are (c, b).
28    """
29    mu, cov = np.asarray(mu0, float), np.asarray(cov0, float)
30    survival = 1.0
31    trace = []
32    for u in controls:
33        next_mu, J = dynamics(mu, u)
34        mu, cov = propagate(next_mu, cov, J, process_cov)
35        qs, margins = [], []
36        for c, b in unsafe_halfspaces:
37            q, m = halfspace_probability(mu, cov, c, b)
38            qs.append(q); margins.append(m)
39            survival *= (1.0 - q)
40        trace.append({"mu": mu.copy(), "cov": cov.copy(), "q": qs, "margin": margins})
41    risk = 1.0 - survival
42    return (risk, trace) if return_trace else risk
43
44
45def fixed_false_warning_recall(scores, labels, false_positive_rate=0.05):
46    """Recall at an empirical threshold allowing at most target false warnings."""
47    scores, labels = np.asarray(scores), np.asarray(labels).astype(bool)
48    negatives = scores[~labels]
49    threshold = np.inf if len(negatives) == 0 else np.quantile(negatives, 1-false_positive_rate)
50    return float(np.mean(scores[labels] >= threshold)), float(threshold)