Saved notebook contents
Code, text, and outputs are displayed as saved in the repository. This documentation build does not execute or revalidate this study.
This notebook has saved text outputs but no embedded figures. The paper PDFs referenced by its code are not included in this checkout.
Checking jaxglitches against published results
The literature on LISA glitch inference reports numbers that this package should be able to reproduce. Only one of those pipelines is public, so most of the comparisons below are against published values rather than against running code. Four are attempted here, and they are the four that need nothing we do not already have:
| Reference | What is checked | |
|---|---|---|
| T2 | Muratore et al. 2025, Table I | the optimal SNR of seven LPF-drawn glitches |
| T4 | Baghi et al. 2022, Sec. VI C | the LPF\(\to\)LISA SNR distribution |
| T6 | Boumerdassi et al. 2026 | the 50/75/90 percentiles of a two-year population |
| T7 | Sauter et al. 2025, Sec. VI | which parameters degenerate for a sub-sample glitch |
Two of these (T2, T7) come out cleanly and settle an ambiguity in the paper being checked. Two (T4, T6) agree in the median and disagree in the tail, and the notebook says as precisely as it can where the disagreement lives.
Everything is computed from the optimal matched-filter SNR
\[\rho^2 \;=\; 4\int_{f_{\min}}^{f_{\max}}\sum_{c\in\{A,E,T\}}
\frac{|\tilde h_c(f)|^2}{S_c(f)}\,\mathrm{d}f ,\]
which is what all four references use. It is evaluated with jaxglitches'
clean_signal_f for \(\tilde h_c\); the noise PSD is not part of the package and comes
from notebooks/noise.py.
import os
import sys
from pathlib import Path
# Pinned to the CPU: this notebook writes a figure into paper/figures and GPU
# reductions are not bit-reproducible, so the file would change between runs even
# though none of the numbers do. Same reason as unequal_vs_equal_arm.ipynb.
os.environ.setdefault("JAX_PLATFORMS", "cpu")
REPO = next(p for p in [Path.cwd(), *Path.cwd().parents] if (p / "pyproject.toml").exists())
sys.path.insert(0, str(REPO / "notebooks"))
sys.path.insert(0, str(REPO / "paper" / "validation"))
import numpy as np
import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
import jax.random as jr
import matplotlib.pyplot as plt
from jaxglitches import clean_signal_f, fisher_matrix, snr as jg_snr
from jaxglitches.catalog_generator import _load_lpf_catalog, run_catalog, RATE_ORDINARY
import noise as ns
from _style import C, FULL_IN, save, set_source
set_source("notebooks/glitch_only/compare_other_codes.ipynb") # for figures/MANIFEST.json
# Two noise settings are used below. `SCIRD` is the package default (the LISA red-book
# amplitudes used everywhere else in the paper); `LDC` is the pair Muratore et al. quote
# for their Spritz run, which is also the closest thing to what Baghi et al. would have
# got from the LDC software. They differ by 1.25x in acceleration and 1.9x in optical
# path, which is a factor of a few in SNR, so the choice is always stated.
SCIRD = dict(A=ns._A_DEF, P=ns._P_DEF) # 3.0e-15 m/s^2/rHz, 15e-12 m/rHz
LDC = dict(A=2.4e-15, P=7.7e-12) # Muratore et al., Spritz noise table
def snr_unit(tau, *, A, P, fmin=1e-5, fmax=1.0, n=4001):
"""Optimal SNR of a glitch of decay time `tau` and unit velocity kick.
The waveform is linear in Delta_v, so rho(Delta_v, tau) = |Delta_v| * snr_unit(tau)
exactly. Computing the unit response once on a grid of tau and scaling is what makes
a 3000-event catalogue cheap, and it keeps the Delta_v dependence explicit -- which
matters below, because two of the papers checked here appear to mean different
things by their amplitude column.
The integral is a trapezoidal sum on a log-spaced grid rather than the uniform-grid
sum of `jaxglitches.snr`; the two are cross-checked against each other in T2.
"""
f = jnp.asarray(np.geomspace(fmin, fmax, n))
S = jnp.stack(ns.psd_tdi1(f, A=A, P=P), axis=-1) # (n, 3) one-sided [1/Hz]
def one(t):
h = clean_signal_f(jnp.array([0.0, 1.0, t]), f)
return jnp.sqrt(4.0 * jnp.trapezoid(jnp.sum(jnp.abs(h) ** 2 / S, axis=-1), f))
return np.asarray(jax.vmap(one)(jnp.atleast_1d(jnp.asarray(tau, dtype=float))))
def snr_of(dv, tau, **kw):
"""Optimal SNR of glitches with velocity kicks `dv` and decay times `tau`."""
return np.abs(np.asarray(dv)) * snr_unit(tau, **kw)
print(f"jax {jax.__version__} on {jax.default_backend()}")
print(f"LPF catalogue: {sum(len(_load_lpf_catalog(r)[0]) for r in ('ordinary','cold'))} "
f"fitted events shipped with the package")
jax 0.10.1 on cpu
LPF catalogue: 530 fitted events shipped with the package
T2 — the seven glitch SNRs of Muratore et al. (2025)
Their Table I lists seven glitches drawn from the LPF catalogue, injected into a 30-day data set, with an optimal SNR quoted for each in second-generation TDI \(AET\). The SNR is generation-independent as long as the PSD matches the waveform — \(\tilde h^{(2)} = (1-D^4)\tilde h^{(1)}\) and \(S^{(2)} = |1-D^4|^2 S^{(1)}\), so \(|\tilde h|^2/S\) is identical — so this is computed in TDI-1.
Two things have to be pinned down before the comparison means anything.
What the amplitude column is. It is headed \(\ln(\Delta v/[\mathrm{m\,s^{-2}}])\) — an acceleration — while \(\Delta v\) is a velocity, and their Eq. (7) gives \(\Delta v = 2A_0\beta_0\). So the tabulated number is either \(\ln\Delta v\) directly (reading A) or \(\ln A_0\), the shapelet amplitude, in which case \(\Delta v = 2A_0\tau\) (reading B). For glitch #3 the two differ by a factor \(2\tau = 369\).
The noise. Their prior table gives \(N_{\rm TM}\in(7,8)\times10^{-12}\,
\mathrm{m\,s^{-2}\,Hz^{-1/2}}\) and \(N_{\rm OMS}\in(2,3)\times10^{-15}\,
\mathrm{m\,Hz^{-1/2}}\), which are interchanged: those are an optical-path amplitude and
an acceleration amplitude respectively, and their own Spritz table has them the other
way round (\(N_{\rm TM}=2.4\times10^{-15}\), \(N_1=7.77\times10^{-12}\)). We use the Spritz
values. Their Eqs. (17)–(20) are then algebraically identical to noise.psd_tdi1, which
is checked below rather than assumed.
There is a test that needs neither of these decisions. Glitches #3 and #118 have the same tabulated amplitude and different \(\tau\), so their SNR ratio is independent of the amplitude convention's overall scale and of the noise normalisation — it depends only on the spectral shape. That ratio alone separates the two readings.
# Muratore et al. 2025, Table I. (name, t0 [s], ln(amplitude), tau [s], quoted SNR)
MURATORE_TABLE_I = [
("#3", 2.97e5, -25.5, 184.47, 72),
("#118", 7.50e5, -25.5, 4.23, 230),
("#90", 1.25e6, -24.4, 43.21, 479),
("#72", 1.75e6, -24.5, 1.86, 637),
("#179", 2.20e6, -23.9, 4.77, 1172),
("#248", 2.25e6, -23.6, 13.96, 1544),
("#41", 2.50e6, -28.3, 5.57, 21),
]
T_OBS_MUR = 2_592_000.0 # 30 days
# --- first, that their noise model really is ours ----------------------------------
# Muratore Eq. (19): S_A = 8 sin^2 x [2 S_TM (3 + 2 cos x + cos 2x) + S_oms (2 + cos x)]
# ours (noise.psd_tdi1): S_A = 8 sin^2 x [4 (1 + cos x + cos^2 x) S_pm + (2 + cos x) S_op]
# and 2(3 + 2 cos x + cos 2x) = 4(1 + cos x + cos^2 x) identically. Likewise Eq. (20)
# against S_T, using sin^4(x/2) = (1 - cos x)^2 / 4.
_f = jnp.asarray(np.geomspace(1e-5, 1.0, 2001))
_x = 2 * jnp.pi * _f * ns._ARM_m / ns._C_SI
_Spm, _Sop = ns._single_link_psd(_f, A=LDC["A"], P=LDC["P"])
_SA_mur = 8 * jnp.sin(_x)**2 * (2 * _Spm * (3 + 2*jnp.cos(_x) + jnp.cos(2*_x))
+ _Sop * (2 + jnp.cos(_x)))
_ST_mur = (16 * _Sop * (1 - jnp.cos(_x)) * jnp.sin(_x)**2
+ 128 * _Spm * jnp.sin(_x)**2 * jnp.sin(_x/2)**4)
_SA, _SE, _ST = ns.psd_tdi1(_f, **LDC)
print("Muratore Eqs. (17)-(20) vs noise.psd_tdi1:")
print(f" max rel. difference, S_A : {float(jnp.max(jnp.abs(_SA_mur/_SA - 1))):.2e}")
print(f" max rel. difference, S_T : {float(jnp.max(jnp.abs(_ST_mur/_ST - 1))):.2e}")
# --- and that the log-grid integral agrees with jaxglitches.snr on a uniform grid ---
_fu = jnp.asarray(np.fft.rfftfreq(int(T_OBS_MUR / 5.0), 5.0)) # dt = 5 s
_psd = ns.psd_tdi1_array(jnp.where(_fu > 0, _fu, 1.0), t_obs=T_OBS_MUR, **LDC)
_p = jnp.array([1e5, 1e-12, 50.0])
_uni = float(jg_snr(clean_signal_f(_p, _fu), _psd))
_log = float(snr_of(1e-12, 50.0, **LDC, fmax=0.1)[0])
print(f"\njaxglitches.snr on the uniform grid : {_uni:.6f}")
print(f"trapezoid on the log grid : {_log:.6f} "
f"(rel. {abs(_log/_uni - 1):.1e})")
Muratore Eqs. (17)-(20) vs noise.psd_tdi1:
max rel. difference, S_A : 1.37e-12
max rel. difference, S_T : 3.46e-11
jaxglitches.snr on the uniform grid : 18.861791
trapezoid on the log grid : 18.861786 (rel. 2.4e-07)
tau_tab = np.array([g[3] for g in MURATORE_TABLE_I])
lnx = np.array([g[2] for g in MURATORE_TABLE_I])
quoted = np.array([g[4] for g in MURATORE_TABLE_I], dtype=float)
# The two readings of the amplitude column.
dv_A = np.exp(lnx) # the column is ln(Delta v)
dv_B = 2 * np.exp(lnx) * tau_tab # the column is ln(A_0), and Delta v = 2 A_0 tau
rho_A = snr_of(dv_A, tau_tab, **LDC)
rho_B = snr_of(dv_B, tau_tab, **LDC)
print(f"{'glitch':>7} {'tau [s]':>9} {'quoted':>8} {'A: dv=e^x':>11} {'ratio':>7} "
f"{'B: dv=2 e^x tau':>16} {'ratio':>8}")
for (name, _, _, tau, q), a, b in zip(MURATORE_TABLE_I, rho_A, rho_B):
print(f"{name:>7} {tau:9.2f} {q:8.0f} {a:11.4g} {a/q:7.2f} {b:16.4g} {b/q:8.1f}")
ok = np.array([n != "#41" for n, *_ in MURATORE_TABLE_I])
print(f"\nreading A, excluding #41: ratios {rho_A[ok]/quoted[ok]}")
print(f" median {np.median(rho_A[ok]/quoted[ok]):.3f}, "
f"worst {np.max(np.abs(rho_A[ok]/quoted[ok] - 1))*100:.1f}% off")
print(f"reading B: ratios span {rho_B.min()/quoted[np.argmin(rho_B)]:.0f} to "
f"{(rho_B/quoted).max():.0f}")
# The shape-only test: #3 and #118 share the tabulated amplitude, so this ratio is
# independent of both the amplitude convention's scale and the noise normalisation.
i3, i118 = 0, 1
print(f"\nSNR(#118)/SNR(#3), a pure test of the spectral shape:")
print(f" quoted {quoted[i118]/quoted[i3]:.3f}")
print(f" reading A {rho_A[i118]/rho_A[i3]:.3f}")
print(f" reading B {rho_B[i118]/rho_B[i3]:.3f}")
# #41 is inconsistent with their own table: relative to #179 (nearly the same tau) the
# quoted SNRs differ by less than the amplitudes do.
i41, i179 = 6, 4
print(f"\n#179 and #41 have tau = {tau_tab[i179]:.2f} and {tau_tab[i41]:.2f} s, so their"
f" SNRs should\nbe in almost the ratio of their amplitudes, "
f"e^({lnx[i179]} - {lnx[i41]}) = {np.exp(lnx[i179]-lnx[i41]):.1f}:")
print(f" reading A gives {rho_A[i179]/rho_A[i41]:.1f}")
print(f" their table gives {quoted[i179]/quoted[i41]:.1f}")
# f_max: does the analysis band matter? The shortest glitch has f_knee = 0.086 Hz.
print(f"\nsensitivity to the upper edge of the band (reading A):")
for fmax in (0.1, 0.5, 2.0):
r = snr_of(dv_A, tau_tab, **LDC, fmax=fmax)
print(f" f_max = {fmax:4.1f} Hz -> {np.array2string(r, precision=1)}")
glitch tau [s] quoted A: dv=e^x ratio B: dv=2 e^x tau ratio
#3 184.47 72 70.07 0.97 2.585e+04 359.0
#118 4.23 230 241.8 1.05 2046 8.9
#90 43.21 479 508.1 1.06 4.391e+04 91.7
#72 1.86 637 662.8 1.04 2465 3.9
#179 4.77 1172 1195 1.02 1.14e+04 9.7
#248 13.96 1544 1507 0.98 4.208e+04 27.3
#41 5.57 21 14.61 0.70 162.7 7.7
reading A, excluding #41: ratios [0.97315935 1.05150686 1.06081563 1.04044612 1.01948661 0.97620413]
median 1.030, worst 6.1% off
reading B: ratios span 8 to 359
SNR(#118)/SNR(#3), a pure test of the spectral shape:
quoted 3.194
reading A 3.452
reading B 0.079
#179 and #41 have tau = 4.77 and 5.57 s, so their SNRs should
be in almost the ratio of their amplitudes, e^(-23.9 - -28.3) = 81.5:
reading A gives 81.8
their table gives 55.8
sensitivity to the upper edge of the band (reading A):
f_max = 0.1 Hz -> [ 70.1 241.8 508.1 662.8 1194.8 1507.3 14.6]
f_max = 0.5 Hz -> [ 70.1 241.8 508.1 662.8 1194.8 1507.3 14.6]
f_max = 2.0 Hz -> [ 70.1 241.8 508.1 662.8 1194.8 1507.3 14.6]
T4 — the LPF\(\to\)LISA SNR distribution of Baghi et al. (2022)
Their Sec. VI C projects each catalogued LPF transient onto a single LISA test mass and integrates their Eq. (28) over \(f\in[10^{-5},1]\,\)Hz. They report SNRs "ranging between \(10^{-2}\) and \(10^{4}\), with 50% of events having SNRs larger than 10".
The catalogue they produced is the one shipped inside jaxglitches — the
effective_glitch_parameters files that lisaglitch's sampler is built on — so this is
a check of the same events through our waveform and our noise, not of a resampled
population. Their noise came from "the LDC software", so LDC is the fairer of the two
settings; SCIRD is shown alongside to size the sensitivity.
beta_cat, level_cat, run_of = [], [], []
for run in ("ordinary", "cold"):
b, l = _load_lpf_catalog(run)
beta_cat.append(b); level_cat.append(l); run_of += [run] * len(b)
beta_cat = np.concatenate(beta_cat)
level_cat = np.concatenate(level_cat)
run_of = np.array(run_of)
print(f"catalogue: {len(beta_cat)} events "
f"({np.sum(run_of=='ordinary')} ordinary, {np.sum(run_of=='cold')} cold)")
print(f" |Delta v| percentiles (1/50/99): "
f"{np.array2string(np.percentile(np.abs(level_cat), [1, 50, 99]), precision=2)} m/s")
print(f" tau percentiles (1/50/99): "
f"{np.array2string(np.percentile(beta_cat, [1, 50, 99]), precision=2)} s")
rho_cat = {}
for lab, kw in (("LDC", LDC), ("SciRD", SCIRD)):
r = snr_of(level_cat, beta_cat, **kw)
rho_cat[lab] = r
print(f"\n{lab:6s}: median {np.median(r):6.2f} frac(SNR>10) {np.mean(r > 10):.2f}")
print(f" 1-99 percentile range [{np.percentile(r,1):.2g}, {np.percentile(r,99):.2g}]")
print(f" full min-max [{r.min():.2g}, {r.max():.2g}]")
print("\nBaghi et al. quote: 1e-2 to 1e4, with 50% of events above SNR 10")
catalogue: 530 events (377 ordinary, 153 cold)
|Delta v| percentiles (1/50/99): [4.52e-15 3.29e-13 9.22e-12] m/s
tau percentiles (1/50/99): [1.00e-01 3.86e+00 4.30e+04] s
LDC : median 8.62 frac(SNR>10) 0.44
1-99 percentile range [0.0016, 1.7e+02]
full min-max [2.3e-15, 3.6e+05]
SciRD : median 5.95 frac(SNR>10) 0.33
1-99 percentile range [0.0013, 1.2e+02]
full min-max [1.8e-15, 2.7e+05]
Baghi et al. quote: 1e-2 to 1e4, with 50% of events above SNR 10
T6 — the two-year population percentiles of Boumerdassi et al. (2026)
They resample the LPF catalogue at roughly 730 glitches per two-year observation and quote glitch-SNR thresholds of \(\{8, 90, 400\}\) as the 50th, 75th and 90th percentiles of that population. Percentiles do not depend on the event rate, so the rate difference (ours is the LPF interarrival rate, about four times theirs) does not enter the comparison — it is reported below only because it is a visible difference in modelling.
Two things are checked separately: whether our resampler reproduces the catalogue it resamples from, and whether the resulting population reproduces their percentiles. The first is ours to get right; the second is a comparison between two papers.
T_2YR = 2 * 365.25 * 86400.0
cat = run_catalog(T_2YR, key=jr.PRNGKey(0))
dv_res, tau_res = np.asarray(cat["Deltav"]), np.asarray(cat["tau"])
print(f"our LPF rate: {RATE_ORDINARY*86400:.2f}/day -> {cat['n_glitches']} events in "
f"2 yr (they quote ~730, i.e. ~1/day)")
BOUM = {50: 8.0, 75: 90.0, 90: 400.0}
Q = [25, 50, 75, 90, 95]
pct = {}
def report(label, rho, key=None):
q = np.percentile(rho, Q)
if key is not None:
pct[key] = q
line = " ".join(f"p{k}={v:>8.3g}" for k, v in zip(Q, q))
rat = ", ".join(f"p{k}: {BOUM[k]/v:.2f}" for k, v in zip(Q, q) if k in BOUM)
print(f"{label:<36s} {line}\n{'':36s} theirs/ours = {rat}")
# The comparison proper is against the catalogue itself: both papers start from the same
# 530 fitted events, so a resampler (theirs a normalising flow, ours a Gaussian KDE in
# (log tau, log|Delta v|)) is a difference in method we should factor out, not in.
print("\nthe catalogue both papers resample:")
report(" Delta v = level", rho_cat["LDC"], key="raw_lvl")
report(" Delta v = 2 alpha tau", snr_of(2*level_cat*beta_cat, beta_cat, **LDC),
key="raw_2at")
# Then our resampler, against the catalogue it resamples: if it were exact these two
# would coincide.
print("\nour two-year resampled population:")
report(" Delta v = level", snr_of(dv_res, tau_res, **LDC), key="res_lvl")
report(" Delta v = 2 alpha tau", snr_of(2*dv_res*tau_res, tau_res, **LDC),
key="res_2at")
print(f"\nBoumerdassi et al.: "
f"{'':9s}p50={BOUM[50]:>8.3g} p75={BOUM[75]:>8.3g} p90={BOUM[90]:>8.3g}")
for seed in (1, 2):
c = run_catalog(T_2YR, key=jr.PRNGKey(seed))
q = np.percentile(snr_of(np.asarray(c["Deltav"]), np.asarray(c["tau"]), **LDC),
[50, 75, 90])
print(f" (seed {seed}: N={c['n_glitches']}, 50/75/90 = {np.array2string(q, precision=0)})")
# Is the mismatch a constant rescaling (same shape, different noise or band) or a
# different shape? A flat ratio across quantiles means the former.
for key, lab in (("raw_lvl", "level"), ("raw_2at", "2 alpha tau")):
r = np.array([BOUM[k] / v for k, v in zip(Q, pct[key]) if k in BOUM])
print(f"\ntheirs/ours across the three quantiles, catalogue read as {lab}: "
f"{np.array2string(r, precision=2)}")
print(f" spread about the mean: {r.max()/r.min():.1f}x"
f" -> {'same shape, scale off by ' + format(1/r.mean(), '.1f') + 'x' if r.max()/r.min() < 1.5 else 'different shape'}")
our LPF rate: 4.30/day -> 3056 events in 2 yr (they quote ~730, i.e. ~1/day)
the catalogue both papers resample:
Delta v = level p25= 2.89 p50= 8.62 p75= 20.8 p90= 43.2 p95= 78.8
theirs/ours = p50: 0.93, p75: 4.33, p90: 9.25
Delta v = 2 alpha tau p25= 4.86 p50= 30.9 p75= 393 p90=1.68e+03 p95=3.23e+03
theirs/ours = p50: 0.26, p75: 0.23, p90: 0.24
our two-year resampled population:
Delta v = level p25= 1.94 p50= 6.63 p75= 20.9 p90= 55.8 p95= 97.7
theirs/ours = p50: 1.21, p75: 4.30, p90: 7.16
Delta v = 2 alpha tau p25= 1.93 p50= 17.3 p75= 245 p90=1.67e+03 p95=4.38e+03
theirs/ours = p50: 0.46, p75: 0.37, p90: 0.24
Boumerdassi et al.: p50= 8 p75= 90 p90= 400
(seed 1: N=3199, 50/75/90 = [ 7. 19. 47.])
(seed 2: N=3209, 50/75/90 = [ 7. 21. 54.])
theirs/ours across the three quantiles, catalogue read as level: [0.93 4.33 9.25]
spread about the mean: 10.0x -> different shape
theirs/ours across the three quantiles, catalogue read as 2 alpha tau: [0.26 0.23 0.24]
spread about the mean: 1.1x -> same shape, scale off by 4.1x
T7 — which parameters degenerate for a short glitch (Sauter et al. 2025)
They report that a maximum-likelihood subtraction leaves most glitches below SNR 2, and that the failures fall in two families. The second is a claim about the likelihood surface: as \(\beta\) approaches the \(4\,\)Hz sampling interval, "the TDI response will fail to capture a change in the shape of the signal, only affecting the overall scale. However, \(\Delta v\) changes also result in a change to the response's scale, creating a degeneracy of parameter choices." — i.e. they attribute the failure to a \(\Delta v\)–\(\beta\) degeneracy.
That is testable directly, and much more cheaply than by sampling: the degeneracy has to
appear as a near-singular block of the Fisher matrix, which jaxglitches.fisher_matrix
returns by automatic differentiation. The expansion of the waveform below the knee
predicts something different from what they say. With \(2\pi f\tau\ll1\),
\[\frac{\Delta\tilde\nu_g}{\nu_0}(f) \;\simeq\;
\frac{-i\,\Delta v\;e^{-2\pi i f(t_0+2\tau)}}{c\,(2\pi f)} ,\]
so \(\tau\) survives only inside the phase, \(\partial_\tau \tilde h = 2\,\partial_{t_0} \tilde h\) exactly, and the degenerate pair is \((t_0,\tau)\) along \(t_0+2\tau = \mathrm{const}\) — while \(\Delta v\), the one parameter they expect to lose, is the one that stays measurable.
# Their sampling rate, and a band that matches it.
DT_SAUTER = 0.25 # 4 Hz
T_OBS_S = 1.0e5
f_s = jnp.asarray(np.fft.rfftfreq(int(T_OBS_S / DT_SAUTER), DT_SAUTER))
psd_s = ns.psd_tdi1_array(jnp.where(f_s > 0, f_s, 1.0), t_obs=T_OBS_S, **SCIRD)
print(f"grid: {len(f_s)} bins, df = {1/T_OBS_S:.1e} Hz, f_max = {float(f_s[-1]):.2f} Hz")
TAU_SCAN = np.geomspace(1e-2, 1e2, 25)
corr, cond, sig_dv = [], [], []
for tau in TAU_SCAN:
G = np.asarray(fisher_matrix(jnp.array([500.0, 1e-12, tau]), f_s, psd_s))
Cv = np.linalg.inv(G)
s = np.sqrt(np.diag(Cv))
R = Cv / np.outer(s, s)
d = np.sqrt(np.diag(G))
corr.append([abs(R[0, 2]), abs(R[1, 2]), abs(R[0, 1])]) # (t0,tau) (dv,tau) (t0,dv)
cond.append(np.linalg.cond(G / np.outer(d, d)))
sig_dv.append(s[1] / 1e-12)
corr, cond, sig_dv = np.array(corr), np.array(cond), np.array(sig_dv)
print(f"\n{'tau [s]':>9} {'f_knee [Hz]':>12} {'|c(t0,tau)|':>12} {'|c(dv,tau)|':>12} "
f"{'|c(t0,dv)|':>11} {'sigma_dv/dv':>12}")
for i in range(0, len(TAU_SCAN), 3):
t = TAU_SCAN[i]
print(f"{t:9.3g} {1/(2*np.pi*t):12.4g} {corr[i,0]:12.4f} {corr[i,1]:12.4f} "
f"{corr[i,2]:11.4f} {sig_dv[i]:12.3g}")
def corr_at(tau):
G = np.asarray(fisher_matrix(jnp.array([500.0, 1e-12, float(tau)]), f_s, psd_s))
Cv = np.linalg.inv(G); sd = np.sqrt(np.diag(Cv))
return Cv / np.outer(sd, sd), sd
R_dt, s_dt = corr_at(DT_SAUTER)
print(f"\nexactly at tau = {DT_SAUTER} s, their sampling interval:")
print(f" |corr(t0, tau)| = {abs(R_dt[0,2]):.4f} <- the degenerate pair")
print(f" |corr(Delta v, tau)| = {abs(R_dt[1,2]):.4f} <- the pair they name")
print(f" Delta v is still measured to {100*s_dt[1]/1e-12:.1f}% at this SNR")
# where does the (t0, tau) correlation cross 0.99?
# corr(t0, tau) falls monotonically with tau, so reverse both for np.interp
tau_99 = float(np.interp(0.99, corr[::-1, 0], TAU_SCAN[::-1]))
print(f" |corr(t0, tau)| exceeds 0.99 below tau = {tau_99:.2f} s "
f"(f_knee = {1/(2*np.pi*tau_99):.3f} Hz)")
grid: 200001 bins, df = 1.0e-05 Hz, f_max = 2.00 Hz
tau [s] f_knee [Hz] |c(t0,tau)| |c(dv,tau)| |c(t0,dv)| sigma_dv/dv
0.01 15.92 1.0000 0.0450 0.0450 0.0501
0.0316 5.033 0.9999 0.0479 0.0479 0.0501
0.1 1.592 0.9995 0.0663 0.0663 0.0502
0.316 0.5033 0.9982 0.1150 0.1148 0.0504
1 0.1592 0.9946 0.2004 0.1993 0.0511
3.16 0.05033 0.9851 0.3482 0.3430 0.0536
10 0.01592 0.9595 0.5286 0.5072 0.0611
31.6 0.005033 0.8987 0.6694 0.6016 0.0826
100 0.001592 0.7867 0.7314 0.5754 0.149
exactly at tau = 0.25 s, their sampling interval:
|corr(t0, tau)| = 0.9986 <- the degenerate pair
|corr(Delta v, tau)| = 0.1025 <- the pair they name
Delta v is still measured to 5.0% at this SNR
|corr(t0, tau)| exceeds 0.99 below tau = 1.97 s (f_knee = 0.081 Hz)
# The direction of the degeneracy, deep in the regime.
p_deg = jnp.array([500.0, 1e-12, 0.05])
G = np.asarray(fisher_matrix(p_deg, f_s, psd_s))
d = np.sqrt(np.diag(G))
w, V = np.linalg.eigh(G / np.outer(d, d))
v_phys = (V[:, 0] / d)
v_phys = v_phys / np.abs(v_phys[0])
s = np.sqrt(np.diag(np.linalg.inv(G)))
print(f"tau = {float(p_deg[2])} s, normalised-Fisher eigenvalues: "
f"{np.array2string(w, precision=4)}")
print(f"softest direction, physical [dt0 : d(Delta v) : dtau]")
print(f" = {v_phys[0]:+.4f} : {v_phys[1]:+.3e} : {v_phys[2]:+.4f}")
print(f" dt0/dtau = {v_phys[0]/v_phys[2]:+.4f} (the expansion predicts -2)")
print(f" sigma_t0/sigma_tau = {s[0]/s[2]:.4f} (predicts 2)")
print(f" Delta v carries {abs(v_phys[1])/abs(v_phys[0]):.1e} of the degenerate "
f"direction: it is not in it at all")
tau = 0.05 s, normalised-Fisher eigenvalues: [2.1877e-04 1.0000e+00 1.9998e+00]
softest direction, physical [dt0 : d(Delta v) : dtau]
= -1.0000 : +2.384e-17 : +0.5001
dt0/dtau = -1.9996 (the expansion predicts -2)
sigma_t0/sigma_tau = 1.9996 (predicts 2)
Delta v carries 2.4e-17 of the degenerate direction: it is not in it at all
The figure
One four-panel figure, fig_literature.pdf, used in the appendix of the paper.
fig, axes = plt.subplots(2, 2, figsize=(FULL_IN, 5.0))
(ax_a, ax_b), (ax_c, ax_d) = axes
# ---- (a) T2: computed against quoted -----------------------------------------------
lim = [8, 4e4]
ax_a.plot(lim, lim, color=C["grey"], lw=0.8, ls="--", zorder=1)
ax_a.fill_between(lim, [0.9*x for x in lim], [1.1*x for x in lim],
color=C["grey"], alpha=0.15, lw=0, zorder=0)
ax_a.scatter(quoted, rho_A, s=22, color=C["blue"], zorder=3,
label=r"$\Delta v = e^{x}$")
ax_a.scatter(quoted, rho_B, s=22, color=C["orange"], marker="s", zorder=3,
label=r"$\Delta v = 2\,e^{x}\tau$")
ax_a.scatter(quoted[~ok], rho_A[~ok], s=70, facecolor="none",
edgecolor=C["red"], lw=1.0, zorder=4)
ax_a.annotate("#41", (quoted[~ok][0], rho_A[~ok][0]), textcoords="offset points",
xytext=(7, -9), fontsize=6, color=C["red"])
ax_a.set(xscale="log", yscale="log", xlim=lim, ylim=(8, 1e5),
xlabel="SNR quoted by Muratore et al., Table I",
ylabel="SNR computed here")
ax_a.set_title(r"(a) seven LPF glitches, $\pm10\%$ band", fontsize=8)
ax_a.legend(loc="upper left", fontsize=6)
# ---- (b) T4: the LPF -> LISA SNR distribution ---------------------------------------
bins = np.geomspace(1e-4, 1e6, 61)
ax_b.hist(rho_cat["LDC"], bins=bins, color=C["blue"], alpha=0.75, label="LDC noise")
ax_b.hist(rho_cat["SciRD"], bins=bins, histtype="step", color=C["red"], lw=1.0,
label="SciRD noise")
ax_b.axvline(10, color=C["grey"], ls="--", lw=0.8)
ax_b.axvline(np.median(rho_cat["LDC"]), color=C["blue"], ls=":", lw=1.0)
ax_b.annotate("Baghi: half the\nevents above this", xy=(10, ax_b.get_ylim()[1]*0.62),
xytext=(4e2, ax_b.get_ylim()[1]*0.80), fontsize=6, color=C["grey"],
ha="left", va="center",
arrowprops=dict(arrowstyle="->", lw=0.6, color=C["grey"]))
ax_b.set(xscale="log", xlabel="optimal SNR in LISA", ylabel="LPF events per bin")
ax_b.set_title(f"(b) the {len(rho_cat['LDC'])} catalogued events; "
f"{100*np.mean(rho_cat['LDC']>10):.0f}% above 10", fontsize=8)
ax_b.legend(loc="upper left", fontsize=6)
# ---- (c) T6: the two-year population ------------------------------------------------
def cdf(ax, rho, **kw):
r = np.sort(np.asarray(rho))
ax.plot(r, np.linspace(0, 1, len(r)), **kw)
cdf(ax_c, rho_cat["LDC"], color=C["blue"], label=r"catalogue, $\Delta v=$ level")
cdf(ax_c, snr_of(2*level_cat*beta_cat, beta_cat, **LDC), color=C["orange"],
label=r"catalogue, $\Delta v=2\alpha\tau$")
cdf(ax_c, snr_of(dv_res, tau_res, **LDC), color=C["blue"], ls="--", lw=0.9,
label="resampled, 2 yr")
for k, v in BOUM.items():
ax_c.plot(v, k/100, marker="*", ms=9, color=C["red"], ls="none", zorder=5)
ax_c.plot([], [], marker="*", ms=9, color=C["red"], ls="none",
label="Boumerdassi et al.")
ax_c.set(xscale="log", xlim=(1e-2, 1e6), ylim=(0, 1),
xlabel="optimal SNR in LISA", ylabel="cumulative fraction")
ax_c.set_title("(c) the LPF population, two ways of reading it", fontsize=8)
ax_c.legend(loc="upper left", fontsize=6)
# ---- (d) T7: the Fisher correlations ------------------------------------------------
ax_d.semilogx(TAU_SCAN, corr[:, 0], color=C["blue"], label=r"$|r(t_0,\tau)|$")
ax_d.semilogx(TAU_SCAN, corr[:, 1], color=C["orange"], label=r"$|r(\Delta v,\tau)|$")
ax_d.semilogx(TAU_SCAN, corr[:, 2], color=C["green"], ls="--",
label=r"$|r(t_0,\Delta v)|$")
ax_d.axvline(DT_SAUTER, color=C["grey"], ls=":", lw=0.9)
ax_d.text(DT_SAUTER*1.15, 0.45, "4 Hz sampling", fontsize=6, color=C["grey"], rotation=90)
ax_d.set(xlabel=r"glitch decay time $\tau$ [s]",
ylabel="Fisher correlation coefficient", ylim=(0, 1.05))
ax_d.set_title(r"(d) the degeneracy is $(t_0,\tau)$, not $(\Delta v,\tau)$", fontsize=8)
ax_d.legend(loc="center left", fontsize=6)
save(fig, "fig_literature",
note="published values are hard-coded from the references; the only input is the "
"LPF catalogue shipped inside the package")
saved paper/figures/fig_literature.pdf
PosixPath('/home/giorgio/Desktop/jaxglitches/paper/figures/fig_literature.pdf')
Summary
print("T2 Muratore et al. 2025, Table I")
print(f" reading A reproduces 6 of 7 quoted SNRs to better than "
f"{np.max(np.abs(rho_A[ok]/quoted[ok]-1))*100:.0f}%; reading B is off by 10-100x")
print(f" and gets the shape ratio wrong ({rho_B[1]/rho_B[0]:.2f} against {quoted[1]/quoted[0]:.2f}).")
print(f" -> the column is ln(Delta v), a velocity, despite its [m s^-2] label.")
print(f" #41 is inconsistent with their own table by a factor "
f"{quoted[6]/rho_A[6]:.2f}.")
print()
print("T4 Baghi et al. 2022")
print(f" median SNR {np.median(rho_cat['LDC']):.1f} and {100*np.mean(rho_cat['LDC']>10):.0f}% above 10, "
f"against their '50% above 10'.")
print(f" the bulk (1-99 pct) spans [{np.percentile(rho_cat['LDC'],1):.1g}, "
f"{np.percentile(rho_cat['LDC'],99):.1g}]; they quote 1e-2 to 1e4.")
print()
print("T6 Boumerdassi et al. 2026")
print(f" catalogue as 'level' 50/75/90 = {np.array2string(pct['raw_lvl'][1:4], precision=0)}")
print(f" catalogue as '2 alpha tau' 50/75/90 = {np.array2string(pct['raw_2at'][1:4], precision=0)}")
print(f" theirs 50/75/90 = [ 8. 90. 400.]")
_rl = np.array([BOUM[k]/v for k, v in zip(Q, pct['raw_lvl']) if k in BOUM])
_r2 = np.array([BOUM[k]/v for k, v in zip(Q, pct['raw_2at']) if k in BOUM])
print(f" under 'level' the median matches but the ratio runs {_rl.min():.1f}-{_rl.max():.1f};")
print(f" under '2 alpha tau' the ratio is flat at {_r2.mean():.2f} +- "
f"{_r2.std():.2f}, i.e. the same shape")
print(f" with the overall SNR scale {1/_r2.mean():.1f}x lower than ours. Undecidable")
print(f" without their event list.")
print()
print("T7 Sauter et al. 2025")
print(f" at tau = {DT_SAUTER} s: |r(t0,tau)| = {abs(R_dt[0,2]):.4f}, "
f"|r(Delta v,tau)| = {abs(R_dt[1,2]):.4f}.")
print(f" the soft direction is dt0 = {v_phys[0]/v_phys[2]:.3f} dtau, i.e. t0 + 2 tau,")
print(f" with no Delta v component. Their diagnosis names the wrong pair.")
T2 Muratore et al. 2025, Table I
reading A reproduces 6 of 7 quoted SNRs to better than 6%; reading B is off by 10-100x
and gets the shape ratio wrong (0.08 against 3.19).
-> the column is ln(Delta v), a velocity, despite its [m s^-2] label.
#41 is inconsistent with their own table by a factor 1.44.
T4 Baghi et al. 2022
median SNR 8.6 and 44% above 10, against their '50% above 10'.
the bulk (1-99 pct) spans [0.002, 2e+02]; they quote 1e-2 to 1e4.
T6 Boumerdassi et al. 2026
catalogue as 'level' 50/75/90 = [ 9. 21. 43.]
catalogue as '2 alpha tau' 50/75/90 = [ 31. 393. 1679.]
theirs 50/75/90 = [ 8. 90. 400.]
under 'level' the median matches but the ratio runs 0.9-9.3;
under '2 alpha tau' the ratio is flat at 0.24 +- 0.01, i.e. the same shape
with the overall SNR scale 4.1x lower than ours. Undecidable
without their event list.
T7 Sauter et al. 2025
at tau = 0.25 s: |r(t0,tau)| = 0.9986, |r(Delta v,tau)| = 0.1025.
the soft direction is dt0 = -2.000 dtau, i.e. t0 + 2 tau,
with no Delta v component. Their diagnosis names the wrong pair.