Hardcoded constants audit
Written after tracking down why the 300-run PE campaign’s phenomd_gpry leg
was silently producing wrong results twice over: PhaseDistanceMarginalLikelihood
defaulted dist_bounds to (1000, 8000) Mpc — the class’s own generic
default — regardless of the run’s actual configured distance prior. For a
campaign whose injections sat as close as 200 Mpc, the distance-marginalization
quadrature grid silently never covered the true distance, producing a wrong
(not just noisier — actively misleading) marginal likelihood over the
intrinsic parameters. It was caught by luck, not by design (see
campaign/plan.md amendment 13’s dist_bounds fix and amendment 16’s
GPryEngine noise-floor/classifier-threshold fix for the sibling bug it
surfaced next to).
The failure shape is: a numeric constant that functions as a physical or algorithmic boundary, tolerance, or resolution is baked into library code instead of being a required input (function argument, constructor argument, or config-file value with no silent library-side fallback). When the baked-in default is wrong for the case at hand, the result is not a crash — it’s a plausible-looking wrong answer. That is strictly worse than a loud failure, because nothing downstream (a corner plot, a PP-plot, a credible interval) looks obviously broken.
This document is the full audit that followed: every file in jaxpe/ was
checked for constants of this shape. It is organized by category, most severe
first within each. Entries marked no override have no argument or config
path at all — editing source is the only way to change them. Entries marked
kwarg exists have an override, but nothing forces a caller (in particular
the CLI, which is supposed to be the single place physical choices are made)
to actually supply it from configuration, so the dangerous default is still
live for any caller that doesn’t know to ask for it.
Status legend: 🔴 not fixed · 🟡 patched at one call site, landmine still live in the library · 🟢 fixed at the source.
Worse than the bug we fixed
IMRPhenomT / IMRPhenomTHM are not real waveform models 🟢 both fixed
jaxpe/gw/cbc_models/phenomt.py— fixed.IMRPhenomT(the dominant (2,2) mode) has been faithfully reimplemented against LALSuite’sLALSimIMRPhenomTHM_internals.c/_fits.c(Estelles et al. 2020, arXiv:2004.08302), replacing the old_compute_phenom_coefficientsplaceholder (fixed dict, ignored itseta/chi1/chi2arguments entirely) with the real closed-form construction: 3.5PN TaylorT3 inspiral + 6 NR-fitted higher-order corrections (6×6jnp.linalg.solve), an arcsinh-parametrized merger frequency ansatz (3×3 solve), and Damour&Nagar 2014’s ringdown ansatz, with phase as the exact analytic integral of each region (no numerical integration). Verified with realtests/test_phenomt.pyassertions against LALSuite’s ownIMRPhenomT(not shape-only): mismatch 1.4e-6–4.6e-6 across a parameter grid spanning mass ratio 1:1–5:1, spins to 0.8, and varying inclination/reference phase — the same precision level as the already-correctIMRPhenomD(see below). Also fixed along the way: the reference-phase convention (LALSuite anchorsphaseatf_ref, not at merger — omitting this left frequency/amplitude correct but phase off by a mass/spin-dependent constant) and thephi_ref/inclination mode-combination convention (phi = π/2 - phiReffed to the spin-weightedYlm, not a separate per-modeexp(i·m·phiRef)rotation — the two aren’t interchangeable away from face-on/edge-on).jaxpe/gw/cbc_models/phenomthm.py— fixed (Phase 2). The 4 higher modes ((2,1),(3,3),(4,4),(5,5)) are now built from the real ansatz, reusing Phase 1’s (2,2) coefficient-solving machinery generalized to arbitrary(l,m): mode-specific PN amplitude coefficients (_pn_amp_coefficients, transcribed fromIMRPhenomTSetHMAmplitudeCoefficients’s per-mode blocks), the same merger/ringdown frequency+phase ansätze as the (2,2) mode with per-mode QNM frequency/damping substitution (_QNM_FRING/_QNM_FDAMP/_QNM_FDAMP_N2), and the fixed per-mode phase offsets (π/2, −π/2, π, π/2 for (2,1)/(3,3)/(4,4)/(5,5)). Verified with realtests/test_phenomthm.pyassertions against LALSuite’sIMRPhenomTHM: polarization-level mismatch 1.4e-10–2.5e-6 across the same parameter grid asIMRPhenomT, plus a direct per-mode comparison againstSimIMRPhenomTHM_Modesand an explicit check that odd-mmodes vanish exactly at the equal-mass/equal-spin symmetry point (a real physical symmetry LALSuite itself enforces, not a numerical coincidence). One latent Phase 1 bug was caught and fixed while building this:_solve_22_amplitude_coefficientswas evaluating the merger-amplitude boundary frequency at the wrong point (the eta-dependenttcut22via the inspiral-only ansatz, instead of the fixedt=-150Mvia the full piecewise_omega22dispatch LALSuite actually uses) — fixed with zero regression on the existing (2,2)-only test suite. Also caught and fixed during Phase 2 itself: the odd-m-vanishes mask was initially applied to all higher modes including the even(4,4), zeroing it out incorrectly at the symmetry point.
tests/test_lalsuite_comparison.py’s xfail(strict=True) marker has been
removed now that IMRPhenomTHM passes for real. New tests/test_phenomd.py
gives IMRPhenomD (previously untested against LAL at all, despite being a
careful, correct transcription) the same real-comparison treatment, and a
new standalone jaxpe/gw/match.py (PSD-weighted mismatch/match,
jax.jit/jax.vmap-composable) backs all of these instead of each test
reinventing the comparison. jaxpe/cli.py’s --waveform phenomthm no
longer refuses to run — it dispatches to the now-validated IMRPhenomTHM
under --domain td, same as --waveform phenomt.
The dist_bounds bug, unfixed at the source 🟢 fixed
jaxpe/gw/likelihood/fd_marginal.py:107—PhaseDistanceMarginalLikelihood.__init__’sdist_boundsis now a required keyword-only argument (no default). Omitting it is an immediateTypeErrorat construction, not a silently wrong answer.jaxpe/gw/likelihood/modes.py:291-292,419-420— the same fix applied to the mode-based marginalizer that backsMarginalizedIntrinsicLikelihood(the ESIGMA+GPry path):dist_min/dist_maxare now required keyword-only arguments (no default) on bothlog_marginal_likelihoodandlog_marginal_likelihood_full.jaxpe/gw/likelihood/marginalized_intrinsic.py—MarginalizedIntrinsicLikelihood.__init__now eagerly validates that itssettingsdict containsdist_min/dist_max, raisingValueErrorat construction rather than deferring the (now-mandatory)TypeErrorto the first, possibly expensive,__call__.jaxpe/cli.py— the default that used to live inside the library classes now lives here instead, as a named, documented constant (DISTANCE_MARGINAL_BOUNDS_FALLBACK_MPC = (100.0, 8000.0)) behind a small helper (_dist_bounds_for_marginalization). Bothrun_pebranches (marginalized_phase_distance,marginalized_intrinsic) call it: it returns the run’s actual resolvedluminosity_distanceprior box whenever that box is non-degenerate, and only substitutes the(100, 8000)Mpc fallback (with a loudwarnings.warn) for the edge case of a degenerate box (e.g. a"fixed"distance prior, which is a nonsensical config for a distance-marginalized likelihood in the first place) — it never silently overrides a real, deliberately-narrow user-configured prior.- All existing callers (
tests/test_fd_marginal.py,tests/test_marginalized.py,tests/test_surrogate.py,examples/07_td_higher_mode_route_comparison.py,examples/08_fd_dominant_mode_route_comparison.py,bin/run_direct_ns.py) already passed these explicitly, so the tightened signatures needed no call-site changes there;tests/test_fd_marginal.py,tests/test_marginalized.py, and the non-GPry/ESIGMA-dependent parts oftests/test_cli.pywere re-run green after the change.
Physics (waveforms / likelihood / data / detectors / priors)
| File:line | Constant | Controls | Exposure | Risk | Status |
|---|---|---|---|---|---|
phenomt.py:33,35,38-39 |
clip to q∈[0.01,1], η∈[1e-4,.25], χ∈[-.99,.99] | silently clips out-of-domain params to the boundary value instead of erroring | no warning on clip | plausible-but-wrong waveform at the prior edges | 🔴 |
esigma.py:109 |
self-force spin-table grid ±0.995, range not exposed (only point count is a kwarg) |
valid spin domain for the horizon-flux term | jnp.interp silently clamps beyond it |
wrong (frozen) horizon flux near-extremal spin, no error | 🔴 |
esigma.py:146,192 |
1.5× safety margin on the Newtonian time-to-merger estimate, no override |
ODE grid duration | none | can truncate high-eccentricity inspirals before merger, silently | 🔴 |
phenomd.py:79 |
fM_CUT = 0.2, module-level, zero override |
dimensionless FD validity cutoff | none | amplitude forced to exactly 0 above cutoff (a truncation, at least louder than extrapolation) | 🔴 |
nrsur7dq4.py |
stated NR-training domain (q≤4, spins≤0.8) in the docstring only | surrogate validity | never checked/clipped/warned in __call__ |
silent extrapolation outside the trained region | 🔴 |
jaxpe/gw/priors.py:27-32,69-75 |
bbh_priors/ebbh_priors: chirp_mass=(20,40), mass_ratio=(0.25,1.0), luminosity_distance=(100,2000), time_width=0.1, eccentricity=(0,0.4) |
default prior support used for inference | kwargs, but nothing forces config-sourcing | true parameter outside the default box → sampler explores a support that excludes the truth; the closest analog in the codebase to the dist_bounds bug |
🔴 |
jaxpe/gw/detectors.py:121-128 |
GMST linear-fit fallback (gps_2000, gmst_2000, valid “around 2020” per docstring) used when lal isn’t importable |
Greenwich Mean Sidereal Time → antenna pattern / sky localization in every likelihood call | silent except ImportError, no warning which branch ran |
silently degrades sky-localization accuracy; the stated validity window is now years stale and unchecked | 🔴 |
jaxpe/gw/psd.py:144 |
noverlap = nperseg // 2 inside welch_psd, zero override |
Welch PSD estimator bias/variance tradeoff | none | biases the noise curve every likelihood evaluation is normalized against | 🔴 |
jaxpe/gw/psd.py:38 vs :90 |
f_low=10.0 (aligo_zdhp_psd) vs f_low=5.0 (lalsim_psd), independent defaults |
PSD low-frequency floor | separate, no shared source of truth | can leak asymmetric near-inf PSD bins across detectors if callers don’t align them |
🔴 |
jaxpe/gw/data.py |
f_min=20.0 / sampling_rate=2048.0 / post_trigger=2.0 / tukey_alpha=0.1 kwargs, repeated across make_injection/make_injections/likelihood_from_strain; duration=8.0 kwarg default now has a config-level adaptive alternative (data.adaptive_duration + jaxpe.gw.estimate_segment_duration, opt-in via examples/configs/campaign_10to80_aligned_adaptive.json, default off) but the kwarg itself is still just a kwarg |
analysis segment/grid and low-frequency cutoff | kwargs, not config-forced (partially addressed for duration) |
forgetting to pass f_min for e.g. a third-generation-detector or low-mass analysis silently drops real signal power below 20 Hz |
🟡 |
likelihood/fd_marginal.py:191,206 |
delta = 0.3, 1e-3 * max\|a\|, inside _check_dominant_mode |
the one safety guard that catches a higher-mode model being misused with this closed-form (dominant-mode-only) likelihood | no override | if poorly matched to a given waveform, the guard can pass while the closed-form is actually invalid | 🔴 |
likelihood/modes.py:288,416 |
n_phi: int = 512 (docstring: covers “network SNR ~ 20”) |
phase-marginalization integral resolution | kwarg, but no assertion/check tying it to the run’s actual SNR | a louder real event silently under-resolves the phase marginal — biased, not obviously wrong | 🔴 |
likelihood/modes.py:290,418 |
tc_half_samples: int = 205 (“±0.1s @ 2048 Hz”) |
time-marginalization window | constant baked to an assumed sampling rate, not derived from it | wrong window width at any other sampling rate | 🔴 |
relative_binning_fd.py:124-125, relative_binning_td.py:59,143,260,431 |
chi=1.0, epsilon=0.5, phase_per_bin=0.5 |
heterodyne/relative-binning resolution | kwargs, adequacy never checked against the actual signal | biased likelihood for a fiducial far from truth or fast phase evolution, with no signal anything’s wrong | 🔴 |
relative_binning_td.py:478 |
1e6 * in_band_max substituted for inf |
PSD masking before the Toeplitz/Gohberg-Semencul solve | none | can leak non-negligible correlated-noise weight for PSDs with wider dynamic range than this factor assumes | 🔴 |
likelihood/base.py:114 |
tukey_alpha: float = 0.1, dataclass default shared by every NetworkLikelihood/ModesNetworkLikelihood |
window taper width → spectral leakage | overridable, but duration/f_min-independent | same “kwarg exists but still ships a plausible default” shape as dist_bounds |
🔴 |
Lower severity / largely self-correcting, noted for completeness:
conditioning.py:36 post_trigger=2.0; conditioning.py:67 0.9 × Nyquist
factor in resolve_f_max; modes.py importance-sampling internals (15.0
e-fold cut, width clip (0.01, 0.25), defense=0.2, qmc_seed=7 — covered by
the escalating-rounds/ESS-target self-healing machinery);
marginalized_intrinsic.py:90 effective_sample_size_floor: float = 0.0
(quality floor defaults off unless a caller opts in); waveform.py
(f_start=20.0, amp_hm=0.6) is an explicitly pedagogical toy model, not
used in production. diagnostics/plots.py:45 quantiles=[0.16,0.5,0.84],
bins=40, smooth=0.9 in corner_plot — display-only, fully overridable,
worth noting only because the 1σ quantile choice doubles as a de facto
credible-interval default read off the plot.
Not flagged (checked, found benign): true physical/mathematical constants
(MTSUN_SI, MPC_SI, C_SI, MRSUN_SI, c, G, π); trivial
indexing/loop literals; PN-series and NR-fit coefficients, which are the
physics itself, not tunable boundaries; core/priors.py’s Sine/Cosine
defaults (0..π, -π/2..π/2) are the full physical range, not a truncation;
core/problem.py/core/transforms.py have no meaningful numeric literals;
diagnostics/stats.py never applies a convergence threshold in code (only in
docstring commentary); flows/trainer.py hyperparameters (n_epochs=8,
batch_size=1024, lr=1e-3) are ordinary, fully-exposed ML defaults with no
physics-correctness failure mode; esigma.py’s inspiral_end_radius=4.0,
ode_eps=1e-8, taper_on/off_seconds, f_lower=20.0 and
pyseobnr_model.py/teobresums_model.py’s d_ref_mpc=1000.0, f_low=20.0,
taper_seconds are all constructor kwargs, not silent in-body defaults.
Sampler / other algorithms (MCMC, GP surrogate, normalizing flow)
| File:line | Constant | Controls | Exposure | Risk | Status |
|---|---|---|---|---|---|
cli.py --gpry-noise-level / --gpry-svm-threshold / --gpry-trust-region-threshold |
default=None |
the amendment-16 GPry high-SNR fix (regressor noise floor, infinities-classifier thresholds) | opt-in only — when unset, the flags are simply omitted rather than computed from the run’s SNR, so GPryEngine falls back to GPry’s own noise_level=1e-2 / ~3σ thresholds |
any future invocation (a new campaign variant, a script, a forgotten flag) silently reproduces the exact GPAcquisitionError crash amendment 16 spent a day diagnosing |
🟡 (fix exists, not a forcing function) |
kernels/adaptation.py:37 |
gamma: float = 1.0 (Robbins-Monro step-size adaptation rate) |
how aggressively the HMC/MALA step size adapts | kwarg exists, but the only caller (global_local.py:414) never threads it from config/CLI |
can oscillate rather than converge under noisy short-window acceptance-rate estimates, with no diagnostic surfaced | 🔴 |
kernels/adaptation.py:38-39,68,78 |
lo=1e-8, hi=1e3 step clip; 1e-8 ensemble-scale floor; 1e-6 ensemble-cov jitter |
HMC/MALA adaptation numerical stability | same pattern — unreachable from config/CLI | usually fails loud (degenerate Cholesky) rather than silently, but still hardcoded in every real run | 🔴 |
sampler/global_local.py:137, flows/interface.py:61 |
interval: float = 5.0 (normalizing-flow rational-quadratic-spline domain) |
the whitened-space range the flow can nonlinearly reshape; effectively linear/identity outside [-interval, interval] |
overridable kwarg, never tuned per-problem or from the actual buffered-sample spread | GW posteriors are routinely multimodal/heavy-tailed (sky ring, distance-inclination degeneracy); a mode or tail beyond ±5σ silently gets a bad global proposal, no error | 🔴 |
sampler/postprocessing.py:92 |
max_tau = 100 fallback when autocorrelation-time estimation fails |
burn-in/thinning basis | warns, but continues using this arbitrary number | plausible but unjustified effective-sample-size accounting | 🔴 |
sampler/global_local.py:502 |
batch_size = 100_000 inside to_physical, zero override (contrast postprocessing.py:122, which has the same constant as an overridable kwarg) |
GPU-memory chunking of the vmapped prior transform | none | OOM risk on smaller GPUs; not a correctness bug (chunking is exact), just inflexible | 🔴 |
surrogate/jax_acquisition.py:166 |
_bucket_size(n, bucket=64), no override at its call site |
JIT-cache padding granularity for the JAX NORA predictive | none | performance/recompilation only — padding is proven exact, not a correctness risk | 🟢 (benign) |
Clean: jaxpe/kernels/hmc.py, mala.py, mmala.py, grw.py, uld.py,
base.py — step_size, scale, n_leapfrog, friction, metric_fn, cov
are all properly required or config-threaded constructor arguments; no
orphaned algorithmic literals found. surrogate/gpry_backend.py and
multifidelity.py pass options/**kwargs through to gpry.Runner
verbatim with no hardcoded GP defaults of their own — the landmine above
lives entirely in the CLI wiring around the surrogate library, not in the
surrogate wrapper itself.
Hardware limits
| File:line | Constant | Controls | Exposure | Status |
|---|---|---|---|---|
bin/condor/generate_campaign_dag.py:47 |
PE_REQUEST_CPUS = 8, module constant |
condor request_cpus + the taskset core mask for all 300 jobs |
no CLI flag | 🔴 |
bin/condor/generate_campaign_dag.py (per-variant dict) |
request_memory = "16GB" (HMC) / "4GB" (GPry variants) |
condor memory request | hardcoded per variant, no override | 🔴 |
bin/condor/generate_campaign_dag.py:169 |
request_cpus = 2 for the POSTPROCESS job |
condor CPU request for the report-generation node | hardcoded | 🔴 |
sampler/global_local.py:502 |
batch_size = 100_000 (listed above too) |
GPU chunking | no override | 🔴 |
--maxjobs, --request-disk, --n-injections, --config are already
proper CLI flags on generate_campaign_dag.py — only the CPU/memory sizing
constants above are not.
How to use this document
Each 🔴 row is a candidate for one of:
- Make the argument required (no default) so a missing value is a loud
TypeError/ConfigErrorinstead of a silent wrong answer — the right fix for anything that gates correctness (distance/mass/spin bounds, PSD floors, resolution parameters tied to SNR). - Source it from
cfg/CLI with no library-side fallback — the right fix for anything that’s legitimately run-dependent but currently only reachable by editing source (hardware sizing,gamma/adaptation constants). - Add a runtime domain check (raise or warn) rather than silently
clipping/extrapolating — the right fix for calibration-range constants
(
phenomt.pyclipping,nrsur7dq4.py’s unchecked NR domain,esigma.py’s spin-table clamp).
Update the status column (🔴 → 🟡 → 🟢) as each is addressed, and note the
call site if a fix only patches one caller rather than the library source
(as dist_bounds currently is at cli.py:363).