Detection Theory Curves¶

import numpy as np

from doppler.detection import (
    det_dwell,
    det_pd,
    det_threshold,
    det_threshold_f,
    det_threshold_noncoherent,
)

PFA = 1e-5
PD_TARGET = 0.9
SNR_DB = [0, 3, 6, 10]  # curves for left panel
MAX_DWELL = 64
DWELL_X = np.arange(1, MAX_DWELL + 1)

ETA = det_threshold(PFA)  # threshold is Pfa-only; computed once

# Left panel: Pd vs dwell for each SNR.
snr_amps = [10 ** (db / 20) for db in SNR_DB]
pd_curves = [[det_pd(snr, int(d), ETA) for d in DWELL_X] for snr in snr_amps]

# Right panel: minimum dwell achieving Pd = 0.9 vs SNR from -3 to 15 dB.
snr_db_sweep = np.linspace(-3, 15, 300)
snr_amp_sweep = 10 ** (snr_db_sweep / 20)
min_dwell = [
    det_dwell(float(s), PD_TARGET, PFA, MAX_DWELL) for s in snr_amp_sweep
]

# Mask SNRs where det_dwell returned -1 (not achievable within MAX_DWELL).
valid = np.array(min_dwell)
mask = valid > 0

# Third panel: what an ESTIMATED noise reference costs.
#
# Every curve above prices a statistic normalised by a KNOWN noise power.
# A burst detector has one burst and must estimate the noise from it, so
# the exact H0 law becomes R^2 = n*F(n, n) -- a fatter tail than the
# chi-square gate prices. det_threshold_f is the exact gate; the penalty
# below is what using the chi-square one instead really delivers.
#
# The degrees of freedom are the trap, and the comparator is DERIVED
# rather than picked to match: under H0 the burst's n prompts give
# sum Re^2 ~ s^2*chi2(n), so a caller treating sum Im^2/n as exactly s^2
# believes R^2 ~ chi2(n) and gates at that quantile. Since
# det_threshold_noncoherent(pfa, M) prices chi2(2M), that is M = n/2.
# Pricing it at n gives 4.8x instead of 41x at n=16.
#
# EVEN n only -- n // 2 would price chi2(n-1) otherwise, which is the same
# gate as n-1 and overstates the penalty by about a fifth. doppler cannot
# express an odd-dof known-noise gate; det_threshold_f handles odd n
# deliberately, because a burst's prompt count is whatever it contained.
PFA_BURST = 1e-3
DOFS = [2, 4, 8, 16, 32, 64, 128]
assert all(n % 2 == 0 for n in DOFS)  # see the note above


def realized_pfa(quantile, n):
    """The pfa a given F(n,n) quantile really buys.

    det_threshold_f is monotone DECREASING in pfa, so bisect on the
    geometric mean of the bracket rather than inverting analytically.
    """
    lo, hi = 1e-15, 0.9999
    for _ in range(200):
        mid = float(np.sqrt(lo * hi))
        if det_threshold_f(mid, n) > quantile:
            lo = mid
        else:
            hi = mid
    return float(np.sqrt(lo * hi))


chi_gates = [
    det_threshold_noncoherent(PFA_BURST, max(n // 2, 1)) ** 2 / n for n in DOFS
]
f_gates = [det_threshold_f(PFA_BURST, n) for n in DOFS]
penalty = [realized_pfa(g, n) / PFA_BURST for g, n in zip(chi_gates, DOFS)]

# The TWO prices, and the estimate quality that causes both.
#
# `penalty` is a multiplier on the false-alarm RATE, paid by a detector
# that ignores the estimation. It is NOT "what an estimated reference
# costs" -- quoted as a cost it names no unit. The cost of doing it RIGHT
# is a higher threshold: detection sensitivity, in dB, which is the number
# a link budget wants.
#
# Both are the same sqrt(2/n) seen from two sides. sigma_hat^2 =
# sum Im^2 / n is sigma^2 * chi2(n)/n, whose relative standard deviation
# is sqrt(2/n) -- 35% at n=16. The floor the detector divides by is known
# only to a third, which is why the correct gate must sit so much higher.
sens_db = [10 * np.log10(f / g) for f, g in zip(f_gates, chi_gates)]
est_rel = [np.sqrt(2.0 / n) for n in DOFS]

# Self-checks: the theory functions must be internally consistent.
# Coherent integration only helps: Pd is non-decreasing in dwell.
for snr_db, pds in zip(SNR_DB, pd_curves):
    assert np.all(np.diff(pds) >= -1e-12), (
        f"Pd not monotone in dwell at {snr_db} dB"
    )

# det_dwell() must return the *minimum* dwell: Pd first crosses the
# target at M, i.e. Pd(M) >= target and Pd(M-1) < target.
for snr_db, snr_amp in zip(SNR_DB, snr_amps):
    m = det_dwell(snr_amp, PD_TARGET, PFA, MAX_DWELL)
    assert m > 0, f"Pd={PD_TARGET} unreachable at {snr_db} dB"
    assert det_pd(snr_amp, m, ETA) >= PD_TARGET, "det_dwell undershoots"
    assert m == 1 or det_pd(snr_amp, m - 1, ETA) < PD_TARGET, (
        "det_dwell is not minimal"
    )
    print(f"SNR {snr_db:+3d} dB: minimum dwell M = {m}")

# Integration gain compensates SNR: a stronger signal never needs a
# longer dwell, so the right-panel curve is non-increasing.
assert np.all(np.diff(valid[mask]) <= 0), "min dwell not monotone in SNR"

# The estimated-noise penalty is real, large, and shrinks as the estimate
# hardens. 41x at n=16 is the number detection_core.h quotes and the
# certification re-derives two independent ways (analytically here, and
# by 2e6 Monte-Carlo draws in the models characterization).
assert all(f > g for f, g in zip(f_gates, chi_gates)), (
    "the F gate must always be the stricter one"
)
assert abs(penalty[DOFS.index(16)] - 41.0) < 1.0, (
    f"penalty at n=16 is {penalty[DOFS.index(16)]:.1f}x, expected ~41x"
)
assert all(a > b for a, b in zip(penalty, penalty[1:])), (
    "the penalty must shrink monotonically as the noise estimate hardens"
)
assert penalty[-1] > 5.0, (
    "even at n=128 the chi-square gate is still materially mispriced"
)

# Price two: doing it correctly costs SENSITIVITY, not false alarms.
i16 = DOFS.index(16)
assert abs(sens_db[i16] - 3.27) < 0.05, (
    f"the correct gate costs {sens_db[i16]:.2f} dB at n=16, expected 3.27"
)
assert all(a > b for a, b in zip(sens_db, sens_db[1:])), (
    "more prompts in the reference must cost less sensitivity"
)
# Doubling n tightens sqrt(2/n) by sqrt(2) and roughly halves the dB cost.
assert all(0.4 < b / a < 0.75 for a, b in zip(sens_db, sens_db[1:])), (
    "the dB cost should roughly halve per doubling of n"
)
assert abs(est_rel[i16] - 0.3536) < 1e-3, "sqrt(2/16) is 35.4%"

for n, pen, db, rel in zip(DOFS, penalty, sens_db, est_rel):
    print(
        f"n={n:<4} noise known to {rel * 100:5.1f}%  ->  "
        f"ignore it: {pen:5.1f}x the priced pfa   |   "
        f"price it: {db:5.2f} dB of sensitivity"
    )