From 541b7f6704ac2e65c54f078f7f91e769e1004458 Mon Sep 17 00:00:00 2001 From: Jasper Seehofer Date: Fri, 3 Jul 2026 10:22:07 +0200 Subject: [PATCH 01/31] test(regression): pin pre-change values for campaign depth + PV marginalization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Physics-change protocol requirement (assert OLD numerical results before the change so the [PHYSICS] diffs are verifiable): - NEW test_bayesian_statistics_host_z_kernel.py: exact pins of single_host_likelihood (volume_deconv + local_ratio, with/without BH mass, low-z window clamp) via stubbed worker globals — the first tests to numerically execute the production host-z kernel. - test_constants.py: pin HOST_DRAW_Z_MAX == 0.5 and GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT == 0.55 (pre-#20 values) + the host-draw <= population-depth ordering constraint. - physical_relations_test.py: d_L anchors at z = 0.5 / 1.5 incl. the closure-truth case dist(1.5, h=0.67) > dist(1.5, h=0.73) that the current fiducial-h parameter-space cap violates. - test_pp_coverage.py: exact-float pins of the tiny-config harness output for both kernels (the determinism test does not pin across code changes). Prep for issues #20 (depth 1.5) and #16 (sigma_v marginalization). Co-Authored-By: Claude Fable 5 --- .../test_bayesian_statistics_host_z_kernel.py | 162 ++++++++++++++++++ .../physical_relations_test.py | 19 ++ master_thesis_code_test/test_constants.py | 24 +++ .../validation/test_pp_coverage.py | 21 +++ 4 files changed, 226 insertions(+) create mode 100644 master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py diff --git a/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py b/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py new file mode 100644 index 00000000..b23f79a3 --- /dev/null +++ b/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py @@ -0,0 +1,162 @@ +"""Numeric regression pins for the host-z kernel in ``single_host_likelihood``. + +These tests pin the exact numerical output of the production Pipeline-B +per-host likelihood at fixed synthetic inputs, so any change to the host +redshift kernel (window bounds, ``norm(loc=host_z, scale=...)`` width, +volume-deconvolution weight) shows up as a deliberate pin update in the same +diff that changes the physics — the /physics-change regression requirement. + +The worker globals normally installed by ``init_worker`` are set directly on +the module with a stub detection-probability object whose grid is wide enough +that the outside-grid quadrature weights are exactly zero. +""" + +import numpy as np +import pytest + +import master_thesis_code.bayesian_inference.bayesian_statistics as bs + +# Synthetic detection: z ~ 0.1 event at h = 0.73 with 5% distance error. +_DET_D_L = 0.47 # Gpc +_DET_D_L_UNC = 0.0235 # Gpc (5%) +_DET_M = 3.3e5 # M_sun, observer-frame (redshifted) mass +_HOST_PHI = 1.2 +_HOST_THETA = 1.0 + + +class _StubDetectionProbability: + """Deterministic, smooth p_det stub with a grid spanning any test window.""" + + def __init__(self) -> None: + self._dl_centers = np.linspace(0.01, 60.0, 300) + + def detection_probability_without_bh_mass_interpolated_zero_fill( + self, + d_L: np.ndarray, + phi: np.ndarray, + theta: np.ndarray, + h: float, + ) -> np.ndarray: + return np.exp(-np.asarray(d_L, dtype=np.float64) / 5.0) + + def detection_probability_with_bh_mass_interpolated( + self, + d_L: np.ndarray, + M_z: np.ndarray, + phi: np.ndarray, + theta: np.ndarray, + h: float, + ) -> np.ndarray: + d = np.asarray(d_L, dtype=np.float64) + m = np.asarray(M_z, dtype=np.float64) + return np.exp(-d / 5.0) * np.exp(-((np.log10(m) - 5.5) ** 2)) + + def _get_or_build_grid(self, h: float) -> tuple: + class _Interp: + grid = (self._dl_centers,) + + return None, _Interp() + + +def _install_worker_globals() -> None: + """Install the module-level worker state single_host_likelihood reads.""" + cov_3d = np.diag([0.02**2, 0.02**2, 0.05**2]) + cov_4d = np.diag([0.02**2, 0.02**2, 0.05**2, 0.1**2]) + + bs.det_index_to_slot = {0: 0} + bs.det_d_L_arr = np.array([_DET_D_L]) + bs.det_d_L_unc_arr = np.array([_DET_D_L_UNC]) + bs.det_M_arr = np.array([_DET_M]) + bs.det_phi_arr = np.array([_HOST_PHI]) + bs.det_theta_arr = np.array([_HOST_THETA]) + + bs.means_3d = np.array([[_HOST_PHI, _HOST_THETA, 1.0]]) + bs.cov_inv_3d = np.array([np.linalg.inv(cov_3d)]) + bs.log_norm_3d = np.array([-0.5 * (3 * np.log(2 * np.pi) + np.linalg.slogdet(cov_3d)[1])]) + + bs.means_4d = np.array([[_HOST_PHI, _HOST_THETA, 1.0, 1.0]]) + # Diagonal 4D covariance: conditional variance of the M_z fraction is its + # marginal variance and the projection vector vanishes. + bs.sigma2_cond_arr = np.array([0.1**2]) + bs.proj_arr = np.zeros((1, 3)) + bs.cov_inv_4d = np.array([np.linalg.inv(cov_4d)]) + bs.log_norm_4d = np.array([-0.5 * (4 * np.log(2 * np.pi) + np.linalg.slogdet(cov_4d)[1])]) + + bs.detection_probability = _StubDetectionProbability() + + +def _run_case( + host_z: float, + host_z_error: float, + normalization_mode: str, + evaluate_with_bh_mass: bool, +) -> list[float]: + _install_worker_globals() + return bs.single_host_likelihood( + host_phiS=_HOST_PHI, + host_qS=_HOST_THETA, + host_z=host_z, + host_z_error=host_z_error, + host_M=3.0e5, + host_M_error=3.0e4, + detection_index=0, + h=0.73, + evaluate_with_bh_mass=evaluate_with_bh_mass, + normalization_mode=normalization_mode, + base_seed=42, + ) + + +def test_kernel_pin_volume_deconv_without_bh_mass() -> None: + """Spec-z-like host (sigma_z = 0.0015) — volume-deconvolved kernel.""" + num, den, w_num, w_den = _run_case(0.10, 0.0015, "volume_deconv", False) + assert num == pytest.approx(PIN_VD_NUM, rel=1e-9) + assert den == pytest.approx(PIN_VD_DEN, rel=1e-9) + assert w_num == 0.0 + assert w_den == 0.0 + + +def test_kernel_pin_local_ratio_without_bh_mass() -> None: + """Same host — bare photo-z Gaussian kernel (local_ratio mode).""" + num, den, w_num, w_den = _run_case(0.10, 0.0015, "local_ratio", False) + assert num == pytest.approx(PIN_LR_NUM, rel=1e-9) + assert den == pytest.approx(PIN_LR_DEN, rel=1e-9) + assert w_num == 0.0 + assert w_den == 0.0 + + +def test_kernel_pin_volume_deconv_with_bh_mass() -> None: + """With-BH-mass path incl. the seeded MC denominator (base_seed = 42).""" + vals = _run_case(0.10, 0.0015, "volume_deconv", True) + assert vals[0] == pytest.approx(PIN_VD_NUM, rel=1e-9) + assert vals[1] == pytest.approx(PIN_VD_DEN, rel=1e-9) + assert vals[2] == pytest.approx(PIN_VD_BH_NUM, rel=1e-9) + assert vals[3] == pytest.approx(PIN_VD_BH_DEN, rel=1e-9) + assert vals[4] == 0.0 + assert vals[5] == 0.0 + + +def test_kernel_pin_low_z_window_clamp() -> None: + """Low-z host (z_g < 4 sigma_z): lower window bound clamps to z = 1e-6. + + The numerator is exactly zero (the event window near z ~ 0.1 has no + overlap with the host kernel at z = 0.004); the denominator and its + outside-grid quadrature weight both depend directly on the kernel width, + making the weight a sensitive canary for any sigma_z change. + """ + num, den, w_num, w_den = _run_case(0.004, 0.0015, "volume_deconv", False) + assert num == 0.0 + assert den == pytest.approx(PIN_CLAMP_DEN, rel=1e-9) + assert w_num == 0.0 + assert w_den == pytest.approx(PIN_CLAMP_W_DEN, rel=1e-9) + + +# ── Pinned values (captured from the code as of the commit adding this file) ── +PIN_VD_NUM = 1622.0066615957417 +PIN_VD_DEN = 0.9153058114326034 +PIN_LR_NUM = 1607.5871614112384 +PIN_LR_DEN = 0.9152832361769563 +PIN_VD_BH_NUM = 4574.227429970933 +PIN_VD_BH_DEN = 0.913118914446828 +PIN_CLAMP_DEN = 0.9958992939448512 +PIN_CLAMP_W_DEN = 0.24151081256750923 diff --git a/master_thesis_code_test/physical_relations_test.py b/master_thesis_code_test/physical_relations_test.py index f28b14ba..9d4cf73c 100644 --- a/master_thesis_code_test/physical_relations_test.py +++ b/master_thesis_code_test/physical_relations_test.py @@ -213,3 +213,22 @@ def test_prescreen_bound_zero_redshift_limit() -> None: """Analytic limit: dist(0) = 0, so the bound vanishes at z_max = 0 (up to float residue of the hypergeometric evaluation, ~1e-15 Gpc).""" assert luminosity_distance_prescreen_gpc(0.0, h=0.73) == pytest.approx(0.0, abs=1e-12) + + +def test_dist_depth_anchor_old_host_draw() -> None: + """Pin d_L at the pre-#20 host-draw depth (z = 0.5, fiducial cosmology).""" + assert dist(0.5, h=0.73) == pytest.approx(2.7428208171409043, rel=1e-12) + + +def test_dist_depth_anchor_campaign() -> None: + """Pin d_L at the Phase-2 campaign depth (z = 1.5): the population-reach + anchor for the pre-screen and the parameter-space d_L bound.""" + assert dist(1.5, h=0.73) == pytest.approx(10.686188506544777, rel=1e-12) + + +def test_dist_campaign_depth_exceeds_fiducial_cap_at_low_h() -> None: + """At closure-truth h = 0.67 the campaign-depth d_L (11.643 Gpc) exceeds + dist(1.5, h=0.73) = 10.686 Gpc — the parameter-space upper bound must + therefore be computed at the lowest campaign h, not the fiducial h.""" + assert dist(1.5, h=0.67) == pytest.approx(11.643160611608486, rel=1e-12) + assert dist(1.5, h=0.67) > dist(1.5, h=0.73) diff --git a/master_thesis_code_test/test_constants.py b/master_thesis_code_test/test_constants.py index 7d94ecb9..f7b9c8fb 100644 --- a/master_thesis_code_test/test_constants.py +++ b/master_thesis_code_test/test_constants.py @@ -3,7 +3,9 @@ import numpy as np from master_thesis_code.constants import ( + GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT, GPC_TO_MPC, + HOST_DRAW_Z_MAX, KM_TO_M, OMEGA_DE, OMEGA_M, @@ -35,3 +37,25 @@ def test_km_to_m() -> None: def test_radian_to_degree() -> None: """360 degrees = 2π radians.""" assert abs(RADIAN_TO_DEGREE * 2 * np.pi - 360.0) < 1e-10 + + +def test_host_draw_depth_pin() -> None: + """Pin the campaign population depth — any change is a /physics-change. + + Pre-#20 value: 0.5 (pre-dt² horizon justification). The Phase-2 campaign + decision (issue #20, 2026-07-03) deliberately flips this to 1.5 in a + [PHYSICS] commit that updates this pin in the same diff. + """ + assert HOST_DRAW_Z_MAX == 0.5 + + +def test_galaxy_catalog_depth_pin() -> None: + """Pin the documented catalogue depth bound (currently unwired in code).""" + assert GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT == 0.55 + + +def test_host_draw_within_population_model() -> None: + """Ordering constraint: the host draw must not exceed the population model + depth (Model1CrossCheck.max_redshift = 1.5, cosmological_model.py) — the + d_L pre-screen derivation relies on this ordering.""" + assert HOST_DRAW_Z_MAX <= 1.5 diff --git a/master_thesis_code_test/validation/test_pp_coverage.py b/master_thesis_code_test/validation/test_pp_coverage.py index f0addf93..17ec5d75 100644 --- a/master_thesis_code_test/validation/test_pp_coverage.py +++ b/master_thesis_code_test/validation/test_pp_coverage.py @@ -83,3 +83,24 @@ def test_medium_config_calibration_band() -> None: assert 0.5 <= volume["coverage"]["68"] <= 0.85 assert bare["coverage"]["68"] < volume["coverage"]["68"] assert abs(volume["map_bias"]) < abs(bare["map_bias"]) + + +def test_tiny_config_exact_value_pins(tiny_bare: dict, tiny_volume: dict) -> None: + """Exact-float regression pins of the harness output (both kernels). + + Unlike the determinism test (which compares two same-code runs), these + pins freeze the CURRENT numerical behaviour so any host-z kernel change + (e.g. a peculiar-velocity sigma_z term) shows up as a deliberate pin + update in the same diff. + """ + bare = tiny_bare["results"]["0.7200"] + assert bare["map_mean"] == pytest.approx(0.6965000000000001, rel=1e-12) + assert bare["map_bias"] == pytest.approx(-0.023499999999999854, rel=1e-9) + assert bare["coverage"]["68"] == pytest.approx(0.25, rel=1e-12) + assert bare["rail_fraction"] == 0.0 + + volume = tiny_volume["results"]["0.7200"] + assert volume["map_mean"] == pytest.approx(0.7185000000000001, rel=1e-12) + assert volume["map_bias"] == pytest.approx(-0.0014999999999998348, rel=1e-6) + assert volume["coverage"]["68"] == pytest.approx(0.625, rel=1e-12) + assert volume["rail_fraction"] == 0.0 From b52ff8d21ea2051b686fd8c674932ccbb5521398 Mon Sep 17 00:00:00 2001 From: Jasper Seehofer Date: Fri, 3 Jul 2026 10:27:46 +0200 Subject: [PATCH 02/31] [PHYSICS] deepen campaign population to z = 1.5 (fix #20) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit User decision on issue #20 (2026-07-03): 'first go for z=1.5 and then see the results and HPC performance'. The pre-dt2 justification (horizon z ~ 0.18, truncation exact) is retired — post-dt2 the EMRI horizon reaches z ~ 1.5+, so the depth is a deliberate population-model choice matching Model1CrossCheck.max_redshift = 1.5. - constants.py: HOST_DRAW_Z_MAX 0.5 -> 1.5; GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT 0.55 -> 1.55 (documented as UNWIRED: the reduced CSV is full-depth, load-time depth is max_redshift via _get_pruned_galaxy_catalog). - main.py injection_campaign: z_cut is now HOST_DRAW_Z_MAX instead of a hardcoded 0.5 — the P_det injection grid must span the host-draw volume or the selection function is blind above the cut. - cosmological_model.py: parameter-space d_L cap computed at the LOWEST campaign h (H_MIN/100 = 0.60 -> 13.0 Gpc) instead of the fiducial h = 0.73 (10.686 Gpc), which silently dropped z >~ 1.35 events in closure runs at h_true = 0.67 (dist(1.5, 0.67) = 11.643 Gpc) via ParameterOutOfBoundsError; plus an explicit ordering guard HOST_DRAW_Z_MAX <= max_redshift (the d_L pre-screen derivation relies on it). - handler.py: rewrite the stale 'truncation is exact' draw docstrings. - main.py fig23: sky-averaged completeness figure now rendered to the campaign depth. - tests: constants pins flipped 0.5->1.5 / 0.55->1.55 in the same diff (physics-change protocol); NEW campaign-depth machinery validation (Schechter/gammaincc finite, bounded, monotone over z in [0.5, 1.5] on synthetic + committed frozen m_th map; frozen-map f_bar(1.0) ~ 0); pixelated dark-draw W_k test pinned to an explicit z_max = 0.5 window (sampler is depth-agnostic; W_k contrast lives at low z — at 1.5 the near-uniform weights drown the Pearson check in multinomial noise). Refs: population reach anchors dist(1.5, h=0.73) = 10.686 Gpc, dist(1.5, h=0.60) = 13.002 Gpc (physical_relations_test.py pins). Full fast suite: 735 passed. Closes #20. Co-Authored-By: Claude Fable 5 --- master_thesis_code/constants.py | 21 +++++++++++++------ master_thesis_code/cosmological_model.py | 21 +++++++++++++++++-- .../galaxy_catalogue/handler.py | 13 +++++++----- master_thesis_code/main.py | 12 ++++++++--- master_thesis_code_test/test_constants.py | 10 ++++----- .../test_dark_event_injection.py | 11 ++++++++-- .../test_pixel_completeness.py | 19 +++++++++++++++++ 7 files changed, 84 insertions(+), 23 deletions(-) diff --git a/master_thesis_code/constants.py b/master_thesis_code/constants.py index 43853d32..1841311a 100644 --- a/master_thesis_code/constants.py +++ b/master_thesis_code/constants.py @@ -66,12 +66,21 @@ FRACTIONAL_MEASURED_MASS_ERROR: float = 1e-8 # fractional error on measured redshifted mass SKY_LOCALIZATION_ERROR: float = 2 / 180 * np.pi # rad, EMRI sky localization error (2 degrees) GALAXY_CATALOG_REDSHIFT_LOWER_LIMIT: float = 0.00001 # minimum redshift for galaxy catalog -GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT: float = 0.55 # maximum redshift for galaxy catalog -# Redshift-horizon margin for the in-catalog host draw (draw_uniform_hosts). The -# EMRI detection horizon is z ≈ 0.18 and the injection campaign already uses -# z_cut = 0.5, so z < 0.5 is safely beyond the horizon -> the truncation is EXACT -# (p_det = 0 beyond) and only removes never-detectable host candidates. -HOST_DRAW_Z_MAX: float = 0.5 +# Documented catalogue depth bound. NOTE: currently UNWIRED — the reduced +# catalogue CSV is written full-depth (no z cut in parse_to_reduced_catalog) +# and the effective load-time depth is Model1CrossCheck.max_redshift via +# _get_pruned_galaxy_catalog. Kept as documentation of the depth the pipeline +# is validated for; raised alongside HOST_DRAW_Z_MAX (issue #20). +GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT: float = 1.55 # maximum redshift for galaxy catalog +# [PHYSICS] Campaign population depth for the in-catalog host draw (issue #20, +# user decision 2026-07-03: "first go for z=1.5 and then see the results and +# HPC performance"). The pre-dt² justification ("horizon z ≈ 0.18, truncation +# EXACT") is retired: after the dt² fix the EMRI horizon reaches z ~ 1.5+, so +# this is a deliberate population-model choice matching +# Model1CrossCheck.max_redshift = 1.5 (cosmological_model.py), NOT a claim +# that p_det = 0 beyond. The injection-campaign z_cut derives from this +# constant (main.py) so the P_det grid always spans the host-draw volume. +HOST_DRAW_Z_MAX: float = 1.5 LUMINOSITY_DISTANCE_THRESHOLD_GPC: float = 1.55 # Gpc, LISA detection horizon for EMRIs # Multiplicative safety margin on the population-derived d_L pre-screen bound # (physical_relations.luminosity_distance_prescreen_gpc). Placeholder pending diff --git a/master_thesis_code/cosmological_model.py b/master_thesis_code/cosmological_model.py index 3e057199..06902e6c 100644 --- a/master_thesis_code/cosmological_model.py +++ b/master_thesis_code/cosmological_model.py @@ -19,7 +19,7 @@ import numpy.typing as npt from scipy.stats import truncnorm -from master_thesis_code.constants import SNR_THRESHOLD +from master_thesis_code.constants import H_MIN, HOST_DRAW_Z_MAX, SNR_THRESHOLD from master_thesis_code.datamodels.detection import ( Detection as Detection, ) @@ -183,7 +183,24 @@ def _apply_model_assumptions(self) -> None: self.parameter_space.e0.upper_limit = 0.2 self.max_redshift = 1.5 - self.parameter_space.luminosity_distance.upper_limit = dist(redshift=self.max_redshift) + # The d_L pre-screen derivation (physical_relations. + # luminosity_distance_prescreen_gpc) and the host draws rely on the + # population model being at least as deep as the host-draw volume. + if HOST_DRAW_Z_MAX > self.max_redshift: + raise ValueError( + f"HOST_DRAW_Z_MAX = {HOST_DRAW_Z_MAX} exceeds the population-model " + f"depth max_redshift = {self.max_redshift}; the host draw would " + "sample beyond the event population." + ) + # [PHYSICS] Parameter-space d_L cap at the LOWEST campaign h, not the + # fiducial h: d_L(z, h) ~ 1/h, so a cap of dist(1.5, h=0.73) = 10.686 Gpc + # silently drops z ≳ 1.35 events in closure runs at h_true = 0.67 + # (dist(1.5, 0.67) = 11.643 Gpc) via ParameterOutOfBoundsError. + # H_MIN/100 = 0.60 is the lower edge of the inference grid, so the cap + # contains the full population for every campaign truth value. + self.parameter_space.luminosity_distance.upper_limit = dist( + redshift=self.max_redshift, h=H_MIN / 100.0 + ) self.luminostity_detection_threshold = 1.55 # as in Hitchikers Guide def emri_distribution(self, M: float, redshift: float) -> float: diff --git a/master_thesis_code/galaxy_catalogue/handler.py b/master_thesis_code/galaxy_catalogue/handler.py index 744e5d53..27e75816 100644 --- a/master_thesis_code/galaxy_catalogue/handler.py +++ b/master_thesis_code/galaxy_catalogue/handler.py @@ -615,9 +615,11 @@ def draw_uniform_hosts( Sampling is WITH REPLACEMENT: hosts are i.i.d. draws from the uniform distribution over the eligible rows, so the same galaxy may be returned more - than once. The truncation ``z < z_max`` is exact for the inference because the - EMRI detection horizon (z ≈ 0.18) lies far below ``z_max`` = 0.5, so the removed - galaxies have p_det = 0 and never contribute a detectable event. + than once. The truncation ``z < z_max`` is a deliberate population-depth + choice matching ``Model1CrossCheck.max_redshift`` (issue #20): after the + dt² fix the EMRI horizon reaches z ~ 1.5+, so p_det is NOT assumed zero + beyond ``z_max`` — the injection-campaign ``z_cut`` derives from the same + constant so draw, injections, and inference share one depth. Args: number_of_hosts: Number of host galaxies to draw (i.i.d., with replacement). @@ -696,8 +698,9 @@ def draw_rate_weighted_hosts( once. As in :meth:`draw_uniform_hosts`, each returned :class:`HostGalaxy` carries z / sky / M / errors straight from its catalog row — there is NO nearest-neighbour snap and NO overwrite of catalog quantities. The - truncation ``z < z_max`` is exact for the inference because the EMRI - detection horizon (z ≈ 0.18) lies far below ``z_max`` = 0.5. + truncation ``z < z_max`` is a deliberate population-depth choice + matching ``Model1CrossCheck.max_redshift`` (issue #20; see + :meth:`draw_uniform_hosts` for the shared-depth rationale). Args: number_of_hosts: Number of host galaxies to draw (i.i.d., with diff --git a/master_thesis_code/main.py b/master_thesis_code/main.py index aaa53fb5..c703da2a 100644 --- a/master_thesis_code/main.py +++ b/master_thesis_code/main.py @@ -660,7 +660,7 @@ def injection_campaign( rng: Random number generator for reproducibility. use_gpu: Whether to use GPU acceleration. """ - from master_thesis_code.constants import INJECTION_CSV_PATH + from master_thesis_code.constants import HOST_DRAW_Z_MAX, INJECTION_CSV_PATH from master_thesis_code.galaxy_catalogue.handler import ParameterSample from master_thesis_code.memory_management import MemoryManagement from master_thesis_code.parameter_estimation.parameter_estimation import ( @@ -700,7 +700,11 @@ def _alarm_handler(signum: int, frame: object) -> None: iteration = 0 parameter_samples_iter: Iterator[ParameterSample] = iter([]) - z_cut = 0.5 # generous margin above max observed detection z ≈ 0.18 + # [PHYSICS] Injection population depth = host-draw depth (issue #20): the + # P_det grid built from these injections must span the full host-draw + # volume, otherwise the selection function is blind above z_cut (the + # pre-#20 hardcoded 0.5 capped the grid while hosts now reach z = 1.5). + z_cut = HOST_DRAW_Z_MAX skipped_high_z = 0 _EMCEE_BATCH = 1000 # large batch to amortize MCMC overhead (93.5% z-rejected) _LOG_INTERVAL = 100 # log every N successful events @@ -1520,12 +1524,14 @@ def _gen_completeness_fk_skymap() -> tuple[object, object] | None: def _gen_sky_averaged_completeness() -> tuple[object, object] | None: try: + from master_thesis_code.constants import HOST_DRAW_Z_MAX from master_thesis_code.galaxy_catalogue.pixel_completeness import from_cache_or_build from master_thesis_code.plotting.completeness_plots import ( plot_sky_averaged_completeness, ) - return plot_sky_averaged_completeness(from_cache_or_build()) + # Track the campaign depth so the figure shows the full host volume. + return plot_sky_averaged_completeness(from_cache_or_build(), z_max=HOST_DRAW_Z_MAX) except (FileNotFoundError, ValueError): _ROOT_LOGGER.info("fig23 skipped: no m_th map / catalog available") return None diff --git a/master_thesis_code_test/test_constants.py b/master_thesis_code_test/test_constants.py index f7b9c8fb..2228c3e0 100644 --- a/master_thesis_code_test/test_constants.py +++ b/master_thesis_code_test/test_constants.py @@ -42,16 +42,16 @@ def test_radian_to_degree() -> None: def test_host_draw_depth_pin() -> None: """Pin the campaign population depth — any change is a /physics-change. - Pre-#20 value: 0.5 (pre-dt² horizon justification). The Phase-2 campaign - decision (issue #20, 2026-07-03) deliberately flips this to 1.5 in a - [PHYSICS] commit that updates this pin in the same diff. + Phase-2 campaign value: 1.5 (issue #20, user decision 2026-07-03), + replacing the pre-dt² 0.5. Flipped from 0.5 in the same [PHYSICS] diff + that changed the constant, per the physics-change protocol. """ - assert HOST_DRAW_Z_MAX == 0.5 + assert HOST_DRAW_Z_MAX == 1.5 def test_galaxy_catalog_depth_pin() -> None: """Pin the documented catalogue depth bound (currently unwired in code).""" - assert GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT == 0.55 + assert GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT == 1.55 def test_host_draw_within_population_model() -> None: diff --git a/master_thesis_code_test/test_dark_event_injection.py b/master_thesis_code_test/test_dark_event_injection.py index e2eae95b..d9c4ab05 100644 --- a/master_thesis_code_test/test_dark_event_injection.py +++ b/master_thesis_code_test/test_dark_event_injection.py @@ -374,12 +374,19 @@ def test_pixelated_dark_draw_pixel_frequency_tracks_W_k() -> None: pc = _pixel_completeness_gradient() rng = np.random.default_rng(202) n = 40000 - hosts = draw_dark_hosts(n, rng, pc, _M_MIN, _M_MAX, h=H) + # Explicit shallow window: the sampler is depth-agnostic, but the W_k + # contrast between pixels lives at low z where f_k differs. At the + # campaign depth (HOST_DRAW_Z_MAX = 1.5) the volume integral is dominated + # by the f ~ 0 shell, W_k becomes near-uniform, and the occupancy/W_k + # Pearson r drowns in multinomial noise — a power issue, not a sampler + # bug. z_max = 0.5 preserves the discriminating power of this check. + z_max = 0.5 + hosts = draw_dark_hosts(n, rng, pc, _M_MIN, _M_MAX, h=H, z_max=z_max) pix = np.array([pc.ang2pix(float(h.phiS), float(h.qS)) for h in hosts]) observed = np.bincount(pix, minlength=pc.npix).astype(np.float64) # Expected per-pixel weight (the categorical the sampler draws from). - z_grid = np.linspace(1e-6, HOST_DRAW_Z_MAX, 4096) + z_grid = np.linspace(1e-6, z_max, 4096) dVc = np.asarray(comoving_volume_element(z_grid, h=H), dtype=np.float64) p_pop = dVc / (1.0 + z_grid) w_k = pc.pixel_dark_weights(z_grid, p_pop, H) diff --git a/master_thesis_code_test/test_pixel_completeness.py b/master_thesis_code_test/test_pixel_completeness.py index 3831b653..f3f44632 100644 --- a/master_thesis_code_test/test_pixel_completeness.py +++ b/master_thesis_code_test/test_pixel_completeness.py @@ -69,6 +69,25 @@ def test_f_k_decreases_to_zero_at_high_redshift() -> None: assert f[-1] < 0.05 +def test_campaign_depth_machinery_valid_to_z_1p5() -> None: + """Issue #20 validation: the per-pixel Schechter machinery stays finite, + bounded, and monotone over the deepened host volume z ∈ [0.5, 1.5]. + + Out there GLADE+ completeness is far below 0.5 (pure-completion regime): + x_th grows to ~1e2 and gammaincc must underflow smoothly to 0 with no + NaN/overflow, on both a synthetic map and the committed frozen map. + """ + z = np.linspace(0.5, 1.5, 40) + for pc in (_mixed_map(), from_cache_or_build()): + f_bar = np.asarray(pc.f_bar(z, _H)) + assert np.all(np.isfinite(f_bar)) + assert np.all((f_bar >= 0.0) & (f_bar <= 1.0)) + assert np.all(np.diff(f_bar) <= 1e-12) + # Frozen GLADE+ map: the deepened shell is pure completion (f ~ 0). + frozen = from_cache_or_build() + assert float(np.asarray(frozen.f_bar(np.array([1.0]), _H))[0]) < 1e-6 + + def test_f_k_in_unit_interval() -> None: """0 <= f_k <= 1 for every pixel (valid and empty) over a z range.""" pc = _mixed_map() From 8568d9fcaa22c59bb6b49deb8dface3e829da4fc Mon Sep 17 00:00:00 2001 From: Jasper Seehofer Date: Fri, 3 Jul 2026 10:31:33 +0200 Subject: [PATCH 03/31] [PHYSICS] marginalize residual host peculiar velocity into the host-z kernel (issue #16) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit User decision on issue #16 (2026-07-03): 'for the next campaign we marginalize it and test the isolated effect in parallel'. Inference-side only (re-evaluate tier) — no re-simulation, no catalogue rebuild. Old kernel (bayesian_statistics.single_host_likelihood): N(z; z_g, sigma_z_cat) [+/- 4 sigma_z_cat window] New kernel: N(z; z_g, sigma_z_eff), sigma_z_eff^2 = sigma_z_cat^2 + sigma_z_pv^2 sigma_z_pv = (1 + z_g) * SIGMA_V_PEC_KM_S / c, SIGMA_V_PEC_KM_S = 200 References: - (1+z) factor: Davis et al. (2011), arXiv:1012.2912, Eqs. (1)/(A1) (z_obs = z_cos + (1 + z_cos) v_pec/c); Davis & Scrimgeour (2014), arXiv:1405.0105, Eq. (9). - Quadrature with catalogue sigma_z: Mastrogiovanni et al. (2023), arXiv:2305.10488, Sec. IV (icarogw/GLADE+ practice). - sigma_v = 200 km/s: Fishbach et al. (2019), arXiv:1807.05667, Sec. 2.2; Chen et al. (2018), arXiv:1712.06531. LISA-EMRI precedent (Laghi et al. 2021, arXiv:2102.01708) uses 500 km/s with the (1+z) factor — kept as a systematics-budget row. Dimensional analysis: [km/s] / [km/s] = dimensionless redshift error; quadrature of two dimensionless sigmas. Limiting cases: sigma_v -> 0 recovers the old kernel exactly; z -> 0 gives sigma_z_pv = sigma_v/c = 6.67e-4 (the low-z LVK convention); at z = 1.5 the factor 2.5 gives 1.67e-3 ~ the catalogue's parse-time PV floor. Applied ONCE at function entry — all eight downstream consumption points (window bounds, Z_g renorm, prior pdf, per-host D_g, with-BH numerator/ denominator, MC proposal + sampling_pdf) flow through the single norm() object, so no double counting inside the likelihood. The catalogue z_error's parse-time PV-CORRECTION floor (0.0015) is a DISTINCT term (correction uncertainty vs residual dispersion) — documented in constants.py. Ball-tree candidate window / pruning intentionally keep the bare catalogue z_error (second-order candidate-list effect). Also: integration-testing twin mirrors the quadrature; pp_coverage gains an inert sigma_z_pv knob (default 0.0 — committed anchor runs stay bit-identical, proven by the unchanged exact-value pins). Kernel regression pins updated in this diff (protocol): at z_g = 0.1, sigma_z_cat = 0.0015 the kernel widens ~11%, integrals shift 0.2-0.5% (e.g. volume_deconv numerator 1622.007 -> 1629.370); pre-change values in the parent commit. Full fast suite: 735 passed. #16 stays OPEN pending the parallel isolated PV value-correction test. Co-Authored-By: Claude Fable 5 --- .../bayesian_inference/bayesian_statistics.py | 34 ++++++++++++++++--- master_thesis_code/constants.py | 13 +++++++ master_thesis_code/validation/pp_coverage.py | 13 ++++++- .../test_bayesian_statistics_host_z_kernel.py | 23 ++++++++----- 4 files changed, 69 insertions(+), 14 deletions(-) diff --git a/master_thesis_code/bayesian_inference/bayesian_statistics.py b/master_thesis_code/bayesian_inference/bayesian_statistics.py index 1fac3b40..697fc017 100644 --- a/master_thesis_code/bayesian_inference/bayesian_statistics.py +++ b/master_thesis_code/bayesian_inference/bayesian_statistics.py @@ -33,7 +33,9 @@ CRAMER_RAO_BOUNDS_OUTPUT_PATH, INJECTION_DATA_DIR, PREPARED_CRAMER_RAO_BOUNDS_PATH, + SIGMA_V_PEC_KM_S, SNR_THRESHOLD, + SPEED_OF_LIGHT_KM_S, H, ) from master_thesis_code.cosmological_model import LamCDMScenario, Model1CrossCheck @@ -1880,6 +1882,25 @@ def single_host_likelihood( integration_limit_sigma_multiplier = 4.0 + # [PHYSICS] Issue #16 (user decision 2026-07-03): marginalize the residual + # host peculiar-velocity dispersion into the host-z kernel. + # sigma_z_pv = (1 + z_g) * sigma_v / c + # Davis et al. (2011), arXiv:1012.2912, Eqs. (1)/(A1) for the (1+z) factor + # (z_obs = z_cos + (1 + z_cos) v_pec / c); added in quadrature to the + # catalogue redshift error per standard practice (Mastrogiovanni et al. + # 2023, arXiv:2305.10488, Sec. IV; EMRI precedent with the (1+z) factor: + # Laghi et al. 2021, arXiv:2102.01708). The catalogue z_error already + # carries GLADE+'s PV-CORRECTION error (or the 0.0015 parse-time floor); + # SIGMA_V_PEC_KM_S is the residual (uncorrected/nonlinear) dispersion on + # top of it. Applied ONCE here: every downstream consumer (window bounds, + # Z_g renormalization, prior pdf, D_g, MC proposal + sampling_pdf) flows + # through this single sigma and the one norm() object below, so the term + # cannot double-count inside the likelihood. The ball-tree candidate + # window and catalogue pruning (handler.py) intentionally keep the bare + # catalogue z_error — a ±1σ, second-order candidate-list effect. + sigma_z_pv = (1.0 + host_z) * SIGMA_V_PEC_KM_S / SPEED_OF_LIGHT_KM_S + host_z_error_eff = float(np.sqrt(host_z_error**2 + sigma_z_pv**2)) + numerator_integration_upper_redshift_limit = dist_to_redshift( _det_d_L + integration_limit_sigma_multiplier * _det_d_L_unc, h=h ) @@ -1887,18 +1908,18 @@ def single_host_likelihood( _det_d_L - integration_limit_sigma_multiplier * _det_d_L_unc, h=h ) denominator_integration_upper_redshift_limit = ( - host_z + integration_limit_sigma_multiplier * host_z_error + host_z + integration_limit_sigma_multiplier * host_z_error_eff ) # [PHYSICS] clamp to z >= 0: for low-z photo-z hosts (z_g < 4 sigma_z) the window # would extend to unphysical z < 0 where comoving_volume_element still returns # positive values, silently adding prior mass to Z_g / D_g (G2b derivation note, # docs/derivations/G2b_host_z_volume_prior.md). Matches B_num's and D(h)'s z_min. denominator_integration_lower_redshift_limit = max( - host_z - integration_limit_sigma_multiplier * host_z_error, 1e-6 + host_z - integration_limit_sigma_multiplier * host_z_error_eff, 1e-6 ) # construct normal distribution for redshift and mass for host galaxy - galaxy_redshift_normal_distribution = norm(loc=host_z, scale=host_z_error) + galaxy_redshift_normal_distribution = norm(loc=host_z, scale=host_z_error_eff) # [PHYSICS] De-rail fix #1 (commission, 2026-07-01): in-catalogue host-redshift prior. # "global"/"local_ratio" use the BARE photo-z Gaussian N(z; z_g, sigma_z) (unchanged @@ -2210,7 +2231,12 @@ def single_host_likelihood_integration_testing( ABS_ERROR = 1e-20 # construct normal distribution for redshift and mass for host galaxy - galaxy_redshift_normal_distribution = norm(loc=possible_host.z, scale=possible_host.z_error) + # [PHYSICS] Issue #16: mirror the production kernel's residual-PV quadrature + # (see single_host_likelihood) so the integration-testing twin stays a + # faithful cross-check of the production path. + _sigma_z_pv = (1.0 + possible_host.z) * SIGMA_V_PEC_KM_S / SPEED_OF_LIGHT_KM_S + _z_error_eff = float(np.sqrt(possible_host.z_error**2 + _sigma_z_pv**2)) + galaxy_redshift_normal_distribution = norm(loc=possible_host.z, scale=_z_error_eff) # Sky localization weight (phi, theta) is inside the GW likelihood Gaussian. # Verified correct by Phase 14 derivation (Sec. 2.7) -- not a source of error. diff --git a/master_thesis_code/constants.py b/master_thesis_code/constants.py index 1841311a..87a19455 100644 --- a/master_thesis_code/constants.py +++ b/master_thesis_code/constants.py @@ -65,6 +65,19 @@ FRACTIONAL_BLACK_HOLE_MASS_CATALOG_ERROR: float = 0.1 # fractional BH mass catalog uncertainty FRACTIONAL_MEASURED_MASS_ERROR: float = 1e-8 # fractional error on measured redshifted mass SKY_LOCALIZATION_ERROR: float = 2 / 180 * np.pi # rad, EMRI sky localization error (2 degrees) +# [PHYSICS] Residual host peculiar-velocity dispersion, marginalized into the +# host-z kernel at inference time (issue #16 decision 2026-07-03): +# sigma_z_pv = (1 + z_g) * SIGMA_V_PEC_KM_S / c, added in quadrature to the +# catalogue sigma_z in bayesian_statistics.single_host_likelihood. +# (1+z) factor: Davis et al. (2011), arXiv:1012.2912, Eqs. (1)/(A1); quadrature +# convention: Mastrogiovanni et al. (2023), arXiv:2305.10488, Sec. IV. +# 200 km/s follows Fishbach et al. (2019), arXiv:1807.05667, Sec. 2.2 and +# Chen et al. (2018), arXiv:1712.06531; the LISA-EMRI precedent (Laghi et al. +# 2021, arXiv:2102.01708, Sec. 4) uses 500 km/s — kept as a systematics-budget +# row, not the default. Distinct from (residual on top of) the GLADE+ +# PV-CORRECTION error already folded into the catalogue z_error at parse time +# (handler.parse_to_reduced_catalog, 0.0015 floor for rows without it). +SIGMA_V_PEC_KM_S: float = 200.0 GALAXY_CATALOG_REDSHIFT_LOWER_LIMIT: float = 0.00001 # minimum redshift for galaxy catalog # Documented catalogue depth bound. NOTE: currently UNWIRED — the reduced # catalogue CSV is written full-depth (no z cut in parse_to_reduced_catalog) diff --git a/master_thesis_code/validation/pp_coverage.py b/master_thesis_code/validation/pp_coverage.py index f471e9a6..14ae1554 100644 --- a/master_thesis_code/validation/pp_coverage.py +++ b/master_thesis_code/validation/pp_coverage.py @@ -191,6 +191,13 @@ class PPCoverageConfig: n_realizations: int = 120 n_events: int = 250 sigma_z: float = 0.035 + # Flat peculiar-velocity redshift-error term, added in quadrature to + # sigma_z for BOTH the generative truth scatter and the inference kernel + # (the calibrated case). Default 0.0 keeps the committed anchor runs + # bit-identical. Issue #16: the production kernel uses + # (1+z) * SIGMA_V_PEC_KM_S / c; this harness knob is flat because its + # sigma_z is flat too. + sigma_z_pv: float = 0.0 sigma_dl_frac: float = 0.05 injected_truths: list[float] = field(default_factory=lambda: [0.62, 0.72, 0.84]) seed: int = 20260701 @@ -223,7 +230,9 @@ def _run_realization( so only the host-z kernel differs between the two estimator variants. """ - sigma_z = config.sigma_z + # One effective sigma for the truth-scatter draw AND the kernel keeps the + # generative model and the inference consistent (calibrated case). + sigma_z = float(np.hypot(config.sigma_z, config.sigma_z_pv)) z_host = _sample_detected_redshifts(h_true, config.n_events, rng) dL_host = comoving_amplitude_of_z(z_host) / h_true dL_obs = np.clip(dL_host + rng.normal(0.0, config.sigma_dl_frac * dL_host), 1e-3, None) @@ -326,6 +335,7 @@ def main(argv: list[str] | None = None) -> None: parser.add_argument("--n-realizations", type=int, default=120) parser.add_argument("--n-events", type=int, default=250) parser.add_argument("--sigma-z", type=float, default=0.035) + parser.add_argument("--sigma-z-pv", type=float, default=0.0) parser.add_argument("--sigma-dl-frac", type=float, default=0.05) parser.add_argument("--truths", type=float, nargs="+", default=[0.62, 0.72, 0.84]) parser.add_argument("--seed", type=int, default=20260701) @@ -337,6 +347,7 @@ def main(argv: list[str] | None = None) -> None: n_realizations=args.n_realizations, n_events=args.n_events, sigma_z=args.sigma_z, + sigma_z_pv=args.sigma_z_pv, sigma_dl_frac=args.sigma_dl_frac, injected_truths=list(args.truths), seed=args.seed, diff --git a/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py b/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py index b23f79a3..0e9678f1 100644 --- a/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py +++ b/master_thesis_code_test/bayesian_inference/test_bayesian_statistics_host_z_kernel.py @@ -151,12 +151,17 @@ def test_kernel_pin_low_z_window_clamp() -> None: assert w_den == pytest.approx(PIN_CLAMP_W_DEN, rel=1e-9) -# ── Pinned values (captured from the code as of the commit adding this file) ── -PIN_VD_NUM = 1622.0066615957417 -PIN_VD_DEN = 0.9153058114326034 -PIN_LR_NUM = 1607.5871614112384 -PIN_LR_DEN = 0.9152832361769563 -PIN_VD_BH_NUM = 4574.227429970933 -PIN_VD_BH_DEN = 0.913118914446828 -PIN_CLAMP_DEN = 0.9958992939448512 -PIN_CLAMP_W_DEN = 0.24151081256750923 +# ── Pinned values ───────────────────────────────────────────────────────────── +# Updated in the [PHYSICS] issue-#16 commit: the host-z kernel now uses +# sigma_z_eff = sqrt(sigma_z_cat^2 + ((1+z_g) SIGMA_V_PEC_KM_S / c)^2), which +# at z_g = 0.1, sigma_z_cat = 0.0015 widens the kernel by ~11% (sigma_z_pv = +# 7.34e-4) and shifts these integrals by 0.2-0.5%. Pre-change values are in +# the parent commit of that diff (physics-change protocol). +PIN_VD_NUM = 1629.3700900543863 +PIN_VD_DEN = 0.9152972692189939 +PIN_LR_NUM = 1611.8260385718838 +PIN_LR_DEN = 0.9152831657144014 +PIN_VD_BH_NUM = 4594.733503494528 +PIN_VD_BH_DEN = 0.9130872910698521 +PIN_CLAMP_DEN = 0.995760331092859 +PIN_CLAMP_W_DEN = 0.2283619655845074 From cbf4ddbe688c13903e6aa78196db8a6f724137f3 Mon Sep 17 00:00:00 2001 From: Jasper Seehofer Date: Fri, 3 Jul 2026 10:33:42 +0200 Subject: [PATCH 04/31] docs: CHANGELOG entries for #20/#16 physics changes + catalogue provenance in DATA_INVENTORY MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - CHANGELOG [Unreleased]: campaign depth z=1.5 (incl. z_cut wiring, low-h d_L cap fix, ordering guard) and residual-PV marginalization (sigma_v = 200 km/s, (1+z) factor, references). - DATA_INVENTORY: new 'Galaxy Catalogue (reduced GLADE+)' section — the z_cmb 8-col full-depth CSV was previously untracked here; records schema, frame provenance (18e9608), rebuild recipe + append-mode gotcha, retired backups, and the frozen m_th map C1 coupling (unchanged by the depth constants since no CSV rebuild occurred). Co-Authored-By: Claude Fable 5 --- CHANGELOG.md | 20 ++++++++++++++++++++ DATA_INVENTORY.md | 17 +++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e954c47f..e9a36f20 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,26 @@ Format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ## [Unreleased] ### Changed +- **[PHYSICS] Campaign population deepened to z = 1.5 (issue #20, decision 2026-07-03):** + `HOST_DRAW_Z_MAX` 0.5 → 1.5 (pre-dt² "horizon z ≈ 0.18, truncation exact" justification + retired); `injection_campaign` `z_cut` now derives from `HOST_DRAW_Z_MAX` (was hardcoded + 0.5 — would have left the P_det grid blind above z = 0.5); parameter-space d_L cap computed + at the lowest campaign h (0.60 → 13.0 Gpc) instead of fiducial h = 0.73 (10.686 Gpc), which + silently dropped z ≳ 1.35 events in h_true = 0.67 closure sims; explicit + `HOST_DRAW_Z_MAX ≤ max_redshift` ordering guard. `GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT` + 0.55 → 1.55 (documented as unwired — the reduced CSV is full-depth, load-time depth is + `Model1CrossCheck.max_redshift`). Completeness machinery validated over z ∈ [0.5, 1.5] + (finite/bounded/monotone; frozen-map f̄(1.0) ≈ 0 — pure-completion regime). Pre-dt² + injection pools remain RETIRED; regeneration at depth 1.5 is part of the Phase-2 campaign. +- **[PHYSICS] Residual host peculiar-velocity dispersion marginalized into the host-z kernel + (issue #16, decision 2026-07-03):** `single_host_likelihood` now uses + σ_z,eff² = σ_z,cat² + ((1+z_g)·σ_v/c)² with `SIGMA_V_PEC_KM_S = 200` (Davis et al. 2011 + Eqs. 1/A1 for the (1+z) factor; Mastrogiovanni et al. 2023 §IV quadrature convention; + Laghi et al. 2021's 500 km/s kept as a systematics-budget row). Inference-side only — + no re-simulation; distinct from the GLADE+ PV-correction error already folded into the + catalogue z_error at parse time. `pp_coverage` gains an inert `--sigma-z-pv` knob + (default 0.0; committed anchor runs bit-identical). Issue #16 stays open for the isolated + PV value-correction impact test. - **Inference is now deterministic (G4):** the with-BH-mass MC denominator draws from a per-host stream derived from `(base_seed, detection_index, host_z, host_M)`; `--seed` reaches the inference layer via `evaluate(..., base_seed=...)` (default 0). diff --git a/DATA_INVENTORY.md b/DATA_INVENTORY.md index 5066fd0f..974994be 100644 --- a/DATA_INVENTORY.md +++ b/DATA_INVENTORY.md @@ -53,6 +53,23 @@ corresponding tier before reporting results. --- +## Galaxy Catalogue (reduced GLADE+) + +The single on-disk input the whole pipeline shares — previously untracked here. + +| Property | Value | +|----------|-------| +| **File** | `master_thesis_code/galaxy_catalogue/reduced_galaxy_catalogue.csv` (headerless, **1.68 GB**, 22 641 048 rows) | +| **Schema** | **8 columns** (order = `_reduced_catalog_column_names()`): RA_deg, Dec_deg, B_mag, **z_cmb**, z_error (PV-correction error folded in quadrature, 0.0015 floor), stellar_mass, stellar_mass_err, z_flag (1=photo-z, 3=spec-z; trailing) | +| **Frame** | **z_cmb** (CMB frame) since `18e9608` (2026-07-02 rebuild; 99.9% rows shifted, median \|Δz\| 6e-4 — `.planning/gate/GATE_SIGNOFF.md:27`) | +| **Depth** | **Full-depth** (no z cut in the writer; max z ≈ 7.03). Effective load-time depth = `Model1CrossCheck.max_redshift` = 1.5 via `_get_pruned_galaxy_catalog`. `GALAXY_CATALOG_REDSHIFT_UPPER_LIMIT` is documentation-only | +| **Rebuild** | `results/commission_20260701/scratch/rebuild_catalog.py` from repo root; **move the old CSV aside first** (writer appends, `mode="a"`); ~77 s full GLADE+ pass on the dev box | +| **Source** | `master_thesis_code/galaxy_catalogue/GLADE+.txt` (6.4 GB, dev box ONLY — cluster cannot rebuild; staging is rsync of the reduced CSV per `/cluster` skill) | +| **Superseded** | `.zhelio_20260702` (z_helio 8-col, 2026-07-01), `.stale6col_mar28` (6-col) — backups next to the live file, RETIRED | +| **Coupled artifact** | `m_th_map_nside32.npy` (frozen per-pixel m_th, C1: byte-identical on injection + inference sides). Built from the full flag-{1,3} catalogue → **unchanged by the 2026-07-03 depth constants** (no CSV rebuild occurred); MUST be regenerated atomically on both sides if the CSV content ever changes | + +--- + ## Dataset Registry ### phase45-seed200-20260501 *(current canonical, post-Tier-3 fix)* From 3273fa5985bb3e762c96af98edcb8d8bf51f4f81 Mon Sep 17 00:00:00 2001 From: Jasper Seehofer Date: Fri, 3 Jul 2026 11:11:12 +0200 Subject: [PATCH 05/31] [PHYSICS] campaign readiness guards: stale-pool gates, M_z edge clamp, injection provenance MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Readiness-sweep fixes (adversarially verified, 2026-07-03) closing the window between catalogue deepening (#20) and new-injection arrival: Stale-pool gates (A2-STALE-POOL-GATE, campaign-blocking): - SimulationDetectionProbability gains expected_z_max/allow_shallow_pool: rejects pools with max z < 0.9 x expected depth and pools MIXING provenance (different z_cut stamps, or stamped + legacy files — the partial-rsync signature). Default None keeps synthetic/test pools ungated. Production constructors (evaluate + combine paths) pass expected_z_max=HOST_DRAW_Z_MAX; the combine path previously had NO check at all. - BayesianStatistics.evaluate: validate_coverage is now a HARD gate (raise below 95% coverage; --allow_low_pdet_coverage escape) instead of a warning buried in 1 of 41 task logs. - Injection writer stamps z_cut + code_rev into every row: h_inj alone cannot discriminate pre-/post-dt2 or shallow/deep pools (0.73 in every era, identical filenames). - Local 80-file pre-dt2 z_cut=0.5 pool archived to simulations/injections_RETIRED_predt2_zcut0p5_20260703/ (filename collision with the regenerated pool); DATA_INVENTORY retirement note. test_sky_selection data-gated tests now self-skip until the depth-1.5 pool lands. [PHYSICS] M_z out-of-grid policy (A2-EXTRAP): p_det with-BH-mass queries outside the injected M_z support were silently LINEARLY extrapolated (scipy fill_value=None semantics) while four comments claimed 'nearest'. Old: linear extrapolation, clipped to [0,1]. New: clamp M_z to the grid edge — true nearest, the documented Phase 44 boundary intent. Affects only out-of-support queries (rare once the pool covers the campaign range); regression test pins edge equality. Injection-loop consistency (A1/A3): - _TIMEOUT_S 30 -> 90 s, aligned with the simulation loop: timed-out events are DROPPED from the pool, so a timeout-rate correlation with (d_L, M_z) at depth 1.5 would bias the p_det grid; counter logged for smoke-test binning. - Symmetric M_z <= M.upper_limit truncation (the CRB path structurally excludes that corner via the Fisher bounds guard; the pool must match). - Pre-dt2 comments retired (93.5% z-rejection, 24/69500 horizon claim). Ops/instrumentation: - --prescreen_audit CLI flag: bypass the quick-SNR skip while logging PRESCREEN_AUDIT (quick, full, params) lines — enables the smoke-run false-negative measurement for PRE_SCREEN_SNR_FACTOR + issue #19. - --combine fast path now writes run_metadata_combine.json (TC-10: previously no git_commit/args recorded for the combine stage). - ParameterSpace d_L default 7 -> 13.1 Gpc (dist(1.5, h=0.60) = 13.0015; the old default silently rejected z >~ 0.9 events in bare constructions); dead luminostity_detection_threshold removed; completeness-plot default z_max tracks HOST_DRAW_Z_MAX. - Integration fixture pool deepened to z <= 1.5 + provenance-stamped (mirrors the production writer; exercises the gates end-to-end). mypy clean (144 files); full fast suite 734 passed / 15 data-gated skips. Co-Authored-By: Claude Fable 5 --- DATA_INVENTORY.md | 14 ++ master_thesis_code/arguments.py | 31 ++++ .../bayesian_inference/bayesian_statistics.py | 22 ++- .../posterior_combination.py | 5 + .../simulation_detection_probability.py | 77 +++++++-- master_thesis_code/cosmological_model.py | 1 - .../datamodels/parameter_space.py | 7 +- master_thesis_code/main.py | 116 ++++++++++--- .../plotting/completeness_plots.py | 4 +- .../test_simulation_detection_probability.py | 162 ++++++++++++++++++ .../integration/test_evaluation_pipeline.py | 7 +- 11 files changed, 409 insertions(+), 37 deletions(-) diff --git a/DATA_INVENTORY.md b/DATA_INVENTORY.md index 974994be..0ce33358 100644 --- a/DATA_INVENTORY.md +++ b/DATA_INVENTORY.md @@ -53,6 +53,20 @@ corresponding tier before reporting results. --- +## ⚠️ Injection pools RETIRED — 2026-07-03 (depth-1.5 campaign prep) + +All pre-dt² / z_cut = 0.5 injection pools are **RETIRED** by the issue-#20 depth change: +local `simulations/injections/` (80 files) moved to +`simulations/injections_RETIRED_predt2_zcut0p5_20260703/`; cluster pools +(`seed43000_Mz`, `seed700`) marked retired in `cluster/datasets.yaml`. The campaign +regenerates a single-h (h_ref = 0.73) pool at z_cut = 1.5 with the **same filenames** — +never mix the eras. Guards: injection rows now carry `z_cut` + `code_rev` provenance +columns; `SimulationDetectionProbability` rejects shallow/mixed pools +(`expected_z_max`, readiness sweep A2-STALE-POOL-GATE); `--evaluate` hard-fails below +95% P_det grid coverage (`--allow_low_pdet_coverage` to override deliberately). + +--- + ## Galaxy Catalogue (reduced GLADE+) The single on-disk input the whole pipeline shares — previously untracked here. diff --git a/master_thesis_code/arguments.py b/master_thesis_code/arguments.py index 87f0a0f9..d7376575 100644 --- a/master_thesis_code/arguments.py +++ b/master_thesis_code/arguments.py @@ -147,6 +147,16 @@ def catalog_only(self) -> bool: """Skip completion integral: set f_i=1, L_comp=0 (catalog-only diagnostic).""" return bool(self._parsed_arguments.catalog_only) + @property + def allow_low_pdet_coverage(self) -> bool: + """Escape hatch for the hard P_det grid-coverage / shallow-pool gate.""" + return bool(self._parsed_arguments.allow_low_pdet_coverage) + + @property + def prescreen_audit(self) -> bool: + """Bypass the quick-SNR early skip while logging (quick, full) SNR pairs.""" + return bool(self._parsed_arguments.prescreen_audit) + @property def combine(self) -> bool: """Indicates whether to combine per-event posteriors into joint H0 posterior.""" @@ -281,6 +291,27 @@ def _parse_arguments(arguments: list[str]) -> argparse.Namespace: default=False, help="Skip completion integral in evaluation: set f_i=1, L_comp=0 (catalog-only diagnostic).", ) + parser.add_argument( + "--allow_low_pdet_coverage", + action="store_true", + default=False, + help=( + "Proceed despite <95%% P_det grid coverage or a shallow injection pool " + "(stale-pool gate). Only for deliberate re-evaluations of archived " + "shallow baselines." + ), + ) + parser.add_argument( + "--prescreen_audit", + action="store_true", + default=False, + help=( + "Audit mode for the quick-SNR pre-screen: compute the full SNR even " + "when the quick gate would skip, and log PRESCREEN_AUDIT lines with " + "(quick_snr, full_snr, params). Smoke-test use only (issue #19 / " + "PRE_SCREEN_SNR_FACTOR re-validation)." + ), + ) parser.add_argument( "--pdet_dl_bins", type=int, diff --git a/master_thesis_code/bayesian_inference/bayesian_statistics.py b/master_thesis_code/bayesian_inference/bayesian_statistics.py index 697fc017..1590236f 100644 --- a/master_thesis_code/bayesian_inference/bayesian_statistics.py +++ b/master_thesis_code/bayesian_inference/bayesian_statistics.py @@ -31,6 +31,7 @@ ) from master_thesis_code.constants import ( CRAMER_RAO_BOUNDS_OUTPUT_PATH, + HOST_DRAW_Z_MAX, INJECTION_DATA_DIR, PREPARED_CRAMER_RAO_BOUNDS_PATH, SIGMA_V_PEC_KM_S, @@ -844,6 +845,7 @@ def evaluate( fisher_cond_threshold: float = 1e16, normalization_mode: str = "volume_deconv", base_seed: int = 0, + allow_low_pdet_coverage: bool = False, ) -> None: self.catalog_only = catalog_only # G4: deterministic seed for the with-BH-mass MC denominator (threaded to @@ -912,6 +914,11 @@ def evaluate( dl_bins=pdet_dl_bins, mass_bins=pdet_mass_bins, estimator=pdet_estimator, # type: ignore[arg-type] + # Stale-pool gate (issue #20): the pool must span the host-draw + # volume; a z_cut = 0.5-era pool at depth 1.5 yields p_det = 0 + # for essentially all events — silent garbage posteriors. + expected_z_max=HOST_DRAW_Z_MAX, + allow_shallow_pool=allow_low_pdet_coverage, ) _LOGGER.debug("Detection probability functions created.") @@ -920,8 +927,19 @@ def evaluate( detection_probability._get_or_build_grid(h_value) _LOGGER.debug("P_det grid pre-warmed for h=%.4f.", h_value) - # Validate P_det grid coverage for observed events - detection_probability.validate_coverage(h_value, self.cramer_rao_bounds) + # Validate P_det grid coverage for observed events — HARD gate + # (readiness sweep A2-STALE-POOL-GATE, 2026-07-03): a warning buried + # in one of 38 per-task logs does not stop a campaign from burning + # its cpu-h budget on p_det = 0 posteriors. + coverage_fraction = detection_probability.validate_coverage(h_value, self.cramer_rao_bounds) + if coverage_fraction < 0.95 and not allow_low_pdet_coverage: + msg = ( + f"P_det grid covers only {coverage_fraction:.1%} of events' " + "4-sigma d_L windows (< 95%). The injection pool is likely stale " + "or too shallow for this event set. Regenerate the pool, or pass " + "--allow_low_pdet_coverage to proceed deliberately." + ) + raise RuntimeError(msg) # Gray et al. (2020), arXiv:1908.06050, Eq. 9 + Gray-Messenger-Veitch 2022, # arXiv:2111.04629 (Change 5): per-HEALPix-pixel completeness f_k(z,Omega,h), diff --git a/master_thesis_code/bayesian_inference/posterior_combination.py b/master_thesis_code/bayesian_inference/posterior_combination.py index c5a11501..0d1708dc 100644 --- a/master_thesis_code/bayesian_inference/posterior_combination.py +++ b/master_thesis_code/bayesian_inference/posterior_combination.py @@ -570,6 +570,7 @@ def combine_posteriors( SimulationDetectionProbability, ) from master_thesis_code.constants import ( # noqa: PLC0415 + HOST_DRAW_Z_MAX, INJECTION_DATA_DIR, OMEGA_DE, OMEGA_M, @@ -579,6 +580,10 @@ def combine_posteriors( detection_probability = SimulationDetectionProbability( injection_data_dir=INJECTION_DATA_DIR, snr_threshold=SNR_THRESHOLD, + # Same stale-pool depth gate as BayesianStatistics.evaluate — the + # combine path recomputes D(h) from the pool and had NO check at + # all (readiness sweep A2-STALE-POOL-GATE, 2026-07-03). + expected_z_max=HOST_DRAW_Z_MAX, ) for h in h_values: detection_probability._get_or_build_grid(h) diff --git a/master_thesis_code/bayesian_inference/simulation_detection_probability.py b/master_thesis_code/bayesian_inference/simulation_detection_probability.py index 3a904b35..f8a8276c 100644 --- a/master_thesis_code/bayesian_inference/simulation_detection_probability.py +++ b/master_thesis_code/bayesian_inference/simulation_detection_probability.py @@ -169,6 +169,8 @@ def __init__( estimator: _Estimator = _DEFAULT_ESTIMATOR, n_sky_bands: int = _DEFAULT_N_SKY_BANDS, _force_unit_weights: bool = False, + expected_z_max: float | None = None, + allow_shallow_pool: bool = False, ) -> None: self._dl_bins = dl_bins self._mass_bins = mass_bins @@ -258,7 +260,7 @@ def __init__( # Validate required columns. "qS" (ecliptic colatitude) is now required: # it carries the response-anisotropy (ecliptic-latitude) axis of p_det # (Change 1). It is already written by every injection campaign - # (main.py:602), so this is a no-op for the canonical CSVs. + # (main.py:injection_campaign), so this is a no-op for the canonical CSVs. required_cols = {"z", "M", "SNR", "h_inj", "luminosity_distance", "qS"} missing = required_cols - set(self._pooled_df.columns) if missing: @@ -276,6 +278,54 @@ def __init__( "luminosity_distance" ].values.astype(np.float64) + # ── Depth / provenance gates (issue #20 stale-pool hazard, 2026-07-03) ── + # Deepening the host draw (HOST_DRAW_Z_MAX -> 1.5) outdates every + # z_cut = 0.5-era pool: deep hosts reach d_L ~ 13 Gpc while a shallow + # pool's survival grid tops out below ~1 Gpc, so p_det = 0 for + # essentially all events — silently valid-looking garbage posteriors. + # The regenerated campaign writes the SAME filenames, so a partial + # rsync / leftover task file mixes eras undetectably by name alone. + # Production constructors pass expected_z_max=HOST_DRAW_Z_MAX; + # tests and synthetic pools leave it None (no depth gate). + if "z_cut" in self._pooled_df.columns: + n_missing = int(self._pooled_df["z_cut"].isna().sum()) + z_cuts = sorted(float(z) for z in self._pooled_df["z_cut"].dropna().unique()) + if n_missing > 0 or len(z_cuts) > 1: + msg = ( + f"Injection pool mixes provenance: z_cut values {z_cuts} plus " + f"{n_missing} rows lacking the column (legacy files). Leftover " + "task files from a retired pool or a partial rsync poison the " + f"survival grid — purge/archive '{injection_data_dir}' and use " + "one consistently-generated pool." + ) + raise ValueError(msg) + else: + logger.warning( + "Injection pool has no provenance columns (z_cut/code_rev) — " + "pre-2026-07-03 writer. Depth is still gated below if " + "expected_z_max is set." + ) + if "code_rev" in self._pooled_df.columns and self._pooled_df["code_rev"].nunique() > 1: + logger.warning( + "Injection pool spans %d code revisions (%s) — legitimate for " + "straggler resubmits after a non-physics fix, but verify none of " + "them changed SNR semantics.", + self._pooled_df["code_rev"].nunique(), + ", ".join(str(c)[:8] for c in self._pooled_df["code_rev"].unique()), + ) + if expected_z_max is not None: + pool_z_max = float(np.max(self._z_arr)) if len(self._z_arr) else 0.0 + if pool_z_max < 0.9 * float(expected_z_max) and not allow_shallow_pool: + msg = ( + f"Injection pool is SHALLOW: max injected z = {pool_z_max:.3f} " + f"< 0.9 x expected_z_max = {0.9 * float(expected_z_max):.3f}. " + "The survival grid cannot cover the host-draw volume " + f"(HOST_DRAW_Z_MAX-era depth mismatch). Regenerate the pool at " + f"the campaign depth, or pass allow_shallow_pool=True (e.g. for " + "a deliberate re-evaluation of an archived shallow baseline)." + ) + raise ValueError(msg) + # h-invariant detection horizon for each injection. # p_det = survival function of the detection horizon, P(d_hor >= d_L), # with d_hor = SNR·d_L/threshold. @@ -761,7 +811,8 @@ def _build_grid_1d( # Exact survival at each center (monotone by construction). p_det_1d = self._survival_at(dl_centers) - # fill_value=None → nearest extrapolation; the public accessors below + # fill_value=None → LINEAR extrapolation outside the grid (scipy + # semantics); harmless here because the public accessors below # override out-of-grid behavior with the EXACT searchsorted survival. return RegularGridInterpolator( (dl_centers,), @@ -876,10 +927,10 @@ def _build_grid_2d( p_det_grid = np.clip(p_det_grid, 0.0, 1.0) - # fill_value=None → nearest-neighbor extrapolation outside grid. Used - # because high-SNR events' integration bounds can exceed the grid; - # the public 2D accessor clamps d_L below first center / above last - # center to keep monotonicity and boundedness. + # fill_value=None → LINEAR extrapolation outside the grid (scipy + # semantics, NOT nearest). Tolerated only because the public 2D + # accessor clamps BOTH axes before querying: d_L below first center / + # above last center, and M_z to the grid edges (true nearest). return RegularGridInterpolator( (dl_centers, M_centers), p_det_grid, @@ -995,7 +1046,9 @@ def detection_probability_with_bh_mass_interpolated( * d_L below the first center → clamp to the first center (survival ≈ 1 there, since the grid starts near d_L = 0). * d_L above the last center → 0 (no injection's horizon reaches there). - * M_z outside the grid range → nearest (``fill_value=None``). + * M_z outside the grid range → clamped to the nearest grid edge. + (``fill_value=None`` alone would LINEARLY extrapolate — made-up but + plausible-looking values; the explicit clip enforces true nearest.) The result is monotone non-increasing in d_L and bounded in [0, 1]. @@ -1030,10 +1083,14 @@ def detection_probability_with_bh_mass_interpolated( dl_max = float(dl_centers[-1]) # d_L below first center → clamp to first center (survival ≈ 1 there); - # d_L above last center → 0. M_z outside range → nearest via - # fill_value=None on the interpolator. + # d_L above last center → 0. M_z outside range → clamp to the grid + # edge (true nearest): RegularGridInterpolator with fill_value=None + # would silently LINEAR-extrapolate outside the M_z axis, inventing + # p_det values beyond the injected mass support (readiness sweep + # A2-EXTRAP, 2026-07-03). dl_query = np.clip(dl_arr, dl_min, dl_max) - result = np.clip(interp_2d(np.column_stack([dl_query, M_arr])), 0.0, 1.0) + M_query = np.clip(M_arr, float(M_centers[0]), float(M_centers[-1])) # noqa: N806 + result = np.clip(interp_2d(np.column_stack([dl_query, M_query])), 0.0, 1.0) # Above the last center the survival is exactly 0. result = np.where(dl_arr > dl_max, 0.0, result) diff --git a/master_thesis_code/cosmological_model.py b/master_thesis_code/cosmological_model.py index 06902e6c..261d61b5 100644 --- a/master_thesis_code/cosmological_model.py +++ b/master_thesis_code/cosmological_model.py @@ -201,7 +201,6 @@ def _apply_model_assumptions(self) -> None: self.parameter_space.luminosity_distance.upper_limit = dist( redshift=self.max_redshift, h=H_MIN / 100.0 ) - self.luminostity_detection_threshold = 1.55 # as in Hitchikers Guide def emri_distribution(self, M: float, redshift: float) -> float: return self.dN_dz_of_mass(M, redshift) * self.R_emri(M) diff --git a/master_thesis_code/datamodels/parameter_space.py b/master_thesis_code/datamodels/parameter_space.py index 7f237ea0..723b5fd4 100644 --- a/master_thesis_code/datamodels/parameter_space.py +++ b/master_thesis_code/datamodels/parameter_space.py @@ -111,7 +111,12 @@ class ParameterSpace: symbol="luminosity_distance", unit="Gpc", lower_limit=0.0, - upper_limit=7, + # dist(HOST_DRAW_Z_MAX=1.5, h=H_MIN/100=0.60) = 13.0015 Gpc — the + # campaign population reach at the lowest grid h. Model1CrossCheck + # recomputes this exactly; the literal here protects bare + # ParameterSpace() constructions from a sub-horizon cap (the old + # 7 Gpc default silently rejected z >~ 0.9 events). + upper_limit=13.1, derivative_epsilon=1e-4, # ~3e-4 × 1 Gpc ≈ 3e-4; use 1e-4 Gpc (= 0.1 Mpc) ) ) # luminosity distance diff --git a/master_thesis_code/main.py b/master_thesis_code/main.py index c703da2a..230fdaef 100644 --- a/master_thesis_code/main.py +++ b/master_thesis_code/main.py @@ -52,6 +52,17 @@ def main() -> None: if arguments.combine: from master_thesis_code.bayesian_inference.posterior_combination import combine_posteriors + # Provenance for the fast-path combine stage (readiness sweep TC-10): + # the early return below skips _write_run_metadata, so the combine + # stage previously recorded no git_commit/args at all. Distinct + # filename avoids colliding with simulation task metadata. + _write_run_metadata( + arguments.working_directory, + arguments.seed, + arguments, + filename="run_metadata_combine.json", + ) + for variant_dir in ["posteriors", "posteriors_with_bh_mass"]: posteriors_dir = os.path.join(arguments.working_directory, variant_dir) if os.path.isdir(posteriors_dir): @@ -106,6 +117,7 @@ def main() -> None: arguments.h_value, rng=rng, use_gpu=arguments.use_gpu, + prescreen_audit=arguments.prescreen_audit, ) if arguments.evaluate: @@ -122,6 +134,7 @@ def main() -> None: normalization_mode=arguments.normalization_mode, # G4: --seed now reaches the inference layer (deterministic MC denominator). base_seed=seed, + allow_low_pdet_coverage=arguments.allow_low_pdet_coverage, ) if arguments.snr_analysis: @@ -236,7 +249,9 @@ def _get_git_commit() -> str: return "unknown" -def _write_run_metadata(working_directory: str, seed: int, arguments: Arguments) -> None: +def _write_run_metadata( + working_directory: str, seed: int, arguments: Arguments, filename: str | None = None +) -> None: metadata = { "git_commit": _get_git_commit(), "timestamp": datetime.datetime.now().isoformat(), @@ -265,11 +280,12 @@ def _write_run_metadata(working_directory: str, seed: int, arguments: Arguments) if slurm_info: metadata["slurm"] = slurm_info - index = arguments.simulation_index - if index > 0 or "SLURM_ARRAY_TASK_ID" in os.environ: - filename = f"run_metadata_{index}.json" - else: - filename = "run_metadata.json" + if filename is None: + index = arguments.simulation_index + if index > 0 or "SLURM_ARRAY_TASK_ID" in os.environ: + filename = f"run_metadata_{index}.json" + else: + filename = "run_metadata.json" metadata_path = os.path.join(working_directory, filename) with open(metadata_path, "w") as f: json.dump(metadata, f, indent=2) @@ -332,6 +348,7 @@ def data_simulation( rng: np.random.Generator | None = None, *, use_gpu: bool = False, + prescreen_audit: bool = False, ) -> None: # conditional imports because they require GPU from master_thesis_code.memory_management import MemoryManagement @@ -499,8 +516,13 @@ def _sigterm_handler(signum: int, frame: object) -> None: # SNR scales as √T for stationary sources; EMRIs chirp so the # 1-yr / 5-yr ratio can be even lower. Factor 0.3 is a conservative - # compromise between the √T bound (0.447) and chirp margin. - if quick_snr < cosmological_model.snr_threshold * PRE_SCREEN_SNR_FACTOR: + # compromise between the √T bound (0.447) and chirp margin — + # calibrated pre-dt² at z <= 0.5; the smoke run re-measures the + # false-negative rate at depth 1.5 via --prescreen_audit. + _quick_gate_failed = ( + quick_snr < cosmological_model.snr_threshold * PRE_SCREEN_SNR_FACTOR + ) + if _quick_gate_failed and not prescreen_audit: signal.alarm(0) _ROOT_LOGGER.info( f"Quick SNR threshold check failed: {np.round(quick_snr, 3)} < {cosmological_model.snr_threshold * PRE_SCREEN_SNR_FACTOR}." @@ -511,6 +533,23 @@ def _sigterm_handler(signum: int, frame: object) -> None: snr = parameter_estimation.compute_signal_to_noise_ratio() signal.alarm(0) warnings.resetwarnings() + if prescreen_audit: + # Greppable audit record for the smoke run (issue #19 + + # PRE_SCREEN_SNR_FACTOR re-validation): quick/full SNR pairs + # with the parameters that drive waveform cost. A false + # negative is full_snr >= threshold while the quick gate + # would have skipped. + _ROOT_LOGGER.info( + "PRESCREEN_AUDIT quick_snr=%.4f full_snr=%.4f gate_would_skip=%s " + "d_L=%.4f M=%.6e e0=%.4f p0=%.4f", + quick_snr, + float(snr), + _quick_gate_failed, + parameter_estimation.parameter_space.luminosity_distance.value, + parameter_estimation.parameter_space.M.value, + parameter_estimation.parameter_space.e0.value, + parameter_estimation.parameter_space.p0.value, + ) except Warning as e: if "Mass ratio" in str(e): _ROOT_LOGGER.warning( @@ -631,7 +670,7 @@ def _sigterm_handler(signum: int, frame: object) -> None: _INJECTION_COLUMNS = ["z", "M", "phiS", "qS", "SNR", "h_inj", "luminosity_distance"] -def _flush_injection_results(results: list[dict[str, float]], csv_path: str) -> None: +def _flush_injection_results(results: list[dict[str, float | str]], csv_path: str) -> None: """Write injection results to CSV (overwrites previous flush).""" import pandas as pd @@ -695,7 +734,7 @@ def _alarm_handler(signum: int, frame: object) -> None: f"index={simulation_index}, output={csv_path}" ) - results: list[dict[str, float]] = [] + results: list[dict[str, float | str]] = [] counter = 0 iteration = 0 parameter_samples_iter: Iterator[ParameterSample] = iter([]) @@ -706,11 +745,24 @@ def _alarm_handler(signum: int, frame: object) -> None: # pre-#20 hardcoded 0.5 capped the grid while hosts now reach z = 1.5). z_cut = HOST_DRAW_Z_MAX skipped_high_z = 0 - _EMCEE_BATCH = 1000 # large batch to amortize MCMC overhead (93.5% z-rejected) + skipped_high_mz = 0 + timeout_count = 0 + # Provenance stamped into every injection row (stale-pool gate, + # readiness sweep A2, 2026-07-03): h_inj alone cannot discriminate + # pre-/post-dt² or shallow/deep pools (0.73 in every era). + code_rev = _get_git_commit() + _EMCEE_BATCH = 1000 # large batch to amortize MCMC overhead (z-rejection ~0 now + # that z_cut = HOST_DRAW_Z_MAX matches the sampler depth; the pre-#20 93.5% + # figure was measured at z_cut = 0.5) _LOG_INTERVAL = 100 # log every N successful events _GPU_FREE_INTERVAL = 50 # free GPU memory every N waveform computations _FLUSH_INTERVAL = 2000 # flush to disk every N events - _TIMEOUT_S = 30 # SNR-only is fast; 30s is generous + # Aligned with the main simulation loop's 90 s alarm (readiness sweep A1, + # 2026-07-03): the injection SNR uses the FULL 5-yr generator and depth 1.5 + # lifts M_z into corners never timing-profiled at the old 30 s budget; + # timed-out events are DROPPED from the pool, so a timeout-rate correlation + # with (d_L, M_z) would bias the p_det grid. Smoke test bins the counter. + _TIMEOUT_S = 90 while counter < simulation_steps: # Sample events from population model @@ -721,10 +773,10 @@ def _alarm_handler(signum: int, frame: object) -> None: parameter_samples_iter = iter(samples_list) sample = next(parameter_samples_iter) - # Importance sampling: skip events beyond the detection horizon. - # All 24/69500 detections in the initial campaign were at z < 0.18. - # Events at z > z_cut have P_det ≈ 0 and waste GPU time on waveforms - # that will never produce detectable SNR. + # Population-depth consistency cut: z_cut = HOST_DRAW_Z_MAX so the + # P_det grid spans exactly the host-draw volume. NOT a p_det = 0 + # claim — post-dt² the horizon reaches z ~ 1.5+ (issue #20 retired + # the pre-dt² "24/69500 detections at z < 0.18" justification). if sample.redshift > z_cut: skipped_high_z += 1 continue @@ -751,6 +803,17 @@ def _alarm_handler(signum: int, frame: object) -> None: # re-lifts). Maggiore (2008) GW Vol. 1 §4.1.4; Babak et al. (2017) arXiv:1703.09722. # Eq. (4.7) in Maggiore (2008) GW Vol. 1 §4.1.4: M_z = M_source·(1+z) redshifted_M = redshifted_mass(sample.M, sample.redshift) # M_z = M·(1+z) + + # Symmetric M_z truncation (readiness sweep A3, 2026-07-03): the CRB + # path structurally excludes M_z > M.upper_limit (Fisher stencil raises + # ParameterOutOfBoundsError at value + 2ε), so the p_det pool must + # exclude the same corner or the selection function includes a + # population the event set cannot contain. Rate-suppressed (mass + # function dies above ~4e6 M_sun source-frame), expected count ~0 — + # the counter makes that measurable. + if redshifted_M > parameter_estimation.parameter_space.M.upper_limit: + skipped_high_mz += 1 + continue parameter_estimation.parameter_space.M.value = redshifted_M # Set luminosity distance with candidate h value (injection pipeline does not use @@ -806,11 +869,14 @@ def _alarm_handler(signum: int, frame: object) -> None: ) continue except TimeoutError: - # G9 gate: the injection alarm budget is _TIMEOUT_S = 30 s, NOT 90 s - # (the old message was wrong); params logged for timeout binning. + # G9 gate: params logged for timeout binning (smoke test checks + # for (d_L, M_z) correlation before full-campaign sizing). + timeout_count += 1 _ROOT_LOGGER.warning( - "Injection waveform/SNR computation timed out (>%ss). Skipping event... params=%s", + "Injection waveform/SNR computation timed out (>%ss, %d total). " + "Skipping event... params=%s", _TIMEOUT_S, + timeout_count, parameter_estimation.parameter_space._parameters_to_dict(), ) continue @@ -828,6 +894,10 @@ def _alarm_handler(signum: int, frame: object) -> None: "SNR": float(snr), "h_inj": h_value, "luminosity_distance": luminosity_distance, + # Provenance (stale-pool gate): z_cut discriminates pool depth + # eras, code_rev ties rows to the generating commit. + "z_cut": z_cut, + "code_rev": code_rev, } ) counter += 1 @@ -839,7 +909,11 @@ def _alarm_handler(signum: int, frame: object) -> None: # Final write _flush_injection_results(results, csv_path) - _ROOT_LOGGER.info(f"Injection campaign complete: {len(results)} events stored to {csv_path}") + _ROOT_LOGGER.info( + f"Injection campaign complete: {len(results)} events stored to {csv_path} " + f"(skipped: {skipped_high_z} high-z, {skipped_high_mz} M_z > bound, " + f"{timeout_count} timeouts @ {_TIMEOUT_S}s)" + ) def evaluate( @@ -855,6 +929,7 @@ def evaluate( fisher_cond_threshold: float = 1e16, normalization_mode: str = "volume_deconv", base_seed: int | None = None, + allow_low_pdet_coverage: bool = False, ) -> None: from master_thesis_code.bayesian_inference.bayesian_statistics import BayesianStatistics @@ -871,6 +946,7 @@ def evaluate( fisher_cond_threshold=fisher_cond_threshold, normalization_mode=normalization_mode, base_seed=base_seed if base_seed is not None else 0, + allow_low_pdet_coverage=allow_low_pdet_coverage, ) diff --git a/master_thesis_code/plotting/completeness_plots.py b/master_thesis_code/plotting/completeness_plots.py index 4b8bf249..ff188b40 100644 --- a/master_thesis_code/plotting/completeness_plots.py +++ b/master_thesis_code/plotting/completeness_plots.py @@ -18,7 +18,7 @@ from matplotlib.axes import Axes from matplotlib.figure import Figure -from master_thesis_code.constants import H +from master_thesis_code.constants import HOST_DRAW_Z_MAX, H from master_thesis_code.galaxy_catalogue.pixel_completeness import PixelCompleteness from master_thesis_code.plotting._colors import REFERENCE, SEQUENTIAL_CMAP, VARIANT_NO_MASS from master_thesis_code.plotting._helpers import _fig_from_ax, get_figure @@ -121,7 +121,7 @@ def plot_sky_averaged_completeness( comp: PixelCompleteness, *, h: float = H, - z_max: float = 0.5, + z_max: float = HOST_DRAW_Z_MAX, n_z: int = 120, ax: Axes | None = None, ) -> tuple[Figure, Axes]: diff --git a/master_thesis_code_test/bayesian_inference/test_simulation_detection_probability.py b/master_thesis_code_test/bayesian_inference/test_simulation_detection_probability.py index 01bed8d4..ca98baef 100644 --- a/master_thesis_code_test/bayesian_inference/test_simulation_detection_probability.py +++ b/master_thesis_code_test/bayesian_inference/test_simulation_detection_probability.py @@ -1713,3 +1713,165 @@ def test_survival_matches_detected_fraction(self, tmp_path: object) -> None: # The accessor IS the exact survival, so they coincide to machine # precision; assert within a loose tolerance per spec. assert p == pytest.approx(direct, abs=1e-6) + + +class TestStalePoolGates: + """Depth/provenance gates (readiness sweep A2-STALE-POOL-GATE, 2026-07-03).""" + + def test_shallow_pool_raises_with_expected_z_max(self, tmp_path: object) -> None: + """A z <= 1.0 pool must be rejected when the host draw expects z_max = 1.5.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + _create_synthetic_injection_csv(d, h_value=0.73, seed=42) # z in [0.01, 1.0] + with pytest.raises(ValueError, match="SHALLOW"): + SimulationDetectionProbability(d, snr_threshold=20.0, expected_z_max=1.5) + + def test_shallow_pool_escape_hatch(self, tmp_path: object) -> None: + """allow_shallow_pool=True permits deliberate shallow-baseline re-evals.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + _create_synthetic_injection_csv(d, h_value=0.73, seed=42) + sdp = SimulationDetectionProbability( + d, snr_threshold=20.0, expected_z_max=1.5, allow_shallow_pool=True + ) + assert sdp is not None + + def test_no_expected_z_max_no_gate(self, tmp_path: object) -> None: + """Default expected_z_max=None leaves synthetic/test pools ungated.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + _create_synthetic_injection_csv(d, h_value=0.73, seed=42) + assert SimulationDetectionProbability(d, snr_threshold=20.0) is not None + + def test_deep_pool_passes_gate(self, tmp_path: object) -> None: + """A pool spanning the host-draw depth passes the gate.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + z = np.linspace(0.02, 1.48, 100) + _create_controlled_injection_csv(d, 0.73, z, np.full(100, 3e5), np.linspace(50.0, 5.0, 100)) + sdp = SimulationDetectionProbability(d, snr_threshold=20.0, expected_z_max=1.5) + assert sdp is not None + + def test_mixed_z_cut_provenance_raises(self, tmp_path: object) -> None: + """Two files with different z_cut stamps = mixed eras -> hard error.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + z = np.linspace(0.02, 1.48, 50) + for z_cut, suffix in ((0.5, "task_001"), (1.5, "task_002")): + n = len(z) + df = pd.DataFrame( + { + "z": z, + "M": np.full(n, 3e5), + "phiS": np.zeros(n), + "qS": np.zeros(n), + "SNR": np.linspace(50.0, 5.0, n), + "h_inj": 0.73, + "luminosity_distance": dist_vectorized(z, h=0.73), + "z_cut": z_cut, + "code_rev": "deadbeef", + } + ) + df.to_csv(f"{d}/injection_h_0p73_{suffix}.csv", index=False) + with pytest.raises(ValueError, match="mixes provenance"): + SimulationDetectionProbability(d, snr_threshold=20.0) + + def test_partial_provenance_raises(self, tmp_path: object) -> None: + """One stamped + one legacy (unstamped) file = partial-rsync signature.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + z = np.linspace(0.02, 1.48, 50) + n = len(z) + base = { + "z": z, + "M": np.full(n, 3e5), + "phiS": np.zeros(n), + "qS": np.zeros(n), + "SNR": np.linspace(50.0, 5.0, n), + "h_inj": 0.73, + "luminosity_distance": dist_vectorized(z, h=0.73), + } + pd.DataFrame(base).to_csv(f"{d}/injection_h_0p73_task_001.csv", index=False) + pd.DataFrame({**base, "z_cut": 1.5, "code_rev": "deadbeef"}).to_csv( + f"{d}/injection_h_0p73_task_002.csv", index=False + ) + with pytest.raises(ValueError, match="mixes provenance"): + SimulationDetectionProbability(d, snr_threshold=20.0) + + def test_uniform_provenance_passes(self, tmp_path: object) -> None: + """Consistently stamped pool constructs cleanly and passes the depth gate.""" + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + z = np.linspace(0.02, 1.48, 50) + n = len(z) + for suffix in ("task_001", "task_002"): + df = pd.DataFrame( + { + "z": z, + "M": np.full(n, 3e5), + "phiS": np.zeros(n), + "qS": np.zeros(n), + "SNR": np.linspace(50.0, 5.0, n), + "h_inj": 0.73, + "luminosity_distance": dist_vectorized(z, h=0.73), + "z_cut": 1.5, + "code_rev": "deadbeef", + } + ) + df.to_csv(f"{d}/injection_h_0p73_{suffix}.csv", index=False) + sdp = SimulationDetectionProbability(d, snr_threshold=20.0, expected_z_max=1.5) + assert sdp is not None + + +class TestMzEdgeClamp: + """M_z outside the grid must clamp to the edge, not linear-extrapolate.""" + + def test_out_of_range_mz_equals_edge_value(self, tmp_path: object) -> None: + from master_thesis_code.bayesian_inference.simulation_detection_probability import ( + SimulationDetectionProbability, + ) + + d = str(tmp_path) + z = np.linspace(0.02, 1.0, 200) + M = np.linspace(1e5, 5e5, 200) # noqa: N806 + _create_controlled_injection_csv(d, 0.73, z, M, np.linspace(80.0, 5.0, 200)) + sdp = SimulationDetectionProbability(d, snr_threshold=20.0) + interp_2d, _ = sdp._get_or_build_grid(0.73) + m_lo = float(interp_2d.grid[1][0]) + m_hi = float(interp_2d.grid[1][-1]) + d_q = 0.3 # well inside the d_L grid + p_at_hi_edge = sdp.detection_probability_with_bh_mass_interpolated( + d_q, m_hi, 0.0, 0.0, h=0.73 + ) + p_beyond_hi = sdp.detection_probability_with_bh_mass_interpolated( + d_q, 10.0 * m_hi, 0.0, 0.0, h=0.73 + ) + p_at_lo_edge = sdp.detection_probability_with_bh_mass_interpolated( + d_q, m_lo, 0.0, 0.0, h=0.73 + ) + p_below_lo = sdp.detection_probability_with_bh_mass_interpolated( + d_q, 0.1 * m_lo, 0.0, 0.0, h=0.73 + ) + assert p_beyond_hi == pytest.approx(p_at_hi_edge, rel=1e-12) + assert p_below_lo == pytest.approx(p_at_lo_edge, rel=1e-12) diff --git a/master_thesis_code_test/integration/test_evaluation_pipeline.py b/master_thesis_code_test/integration/test_evaluation_pipeline.py index 410660cc..24971ea0 100644 --- a/master_thesis_code_test/integration/test_evaluation_pipeline.py +++ b/master_thesis_code_test/integration/test_evaluation_pipeline.py @@ -39,7 +39,10 @@ def _create_synthetic_injection_csvs(directory: str) -> None: for h_value in [0.70, 0.73, 0.80]: n = 200 h_label = f"{h_value:.2f}".replace(".", "p") - z = rng.uniform(0.01, 1.0, size=n) + # Span the campaign host-draw depth (HOST_DRAW_Z_MAX = 1.5) so the + # production stale-pool depth gate passes, and stamp provenance + # columns exactly like the production injection writer does. + z = rng.uniform(0.01, 1.5, size=n) df = pd.DataFrame( { "z": z, @@ -53,6 +56,8 @@ def _create_synthetic_injection_csvs(directory: str) -> None: ), "h_inj": h_value, "luminosity_distance": z * 4.0, + "z_cut": 1.5, + "code_rev": "fixture", } ) df.to_csv(f"{directory}/injection_h_{h_label}_task_001.csv", index=False) From 97b1571b187ef3171f172a8ce4bdd0c3c43d4ede Mon Sep 17 00:00:00 2001 From: Jasper Seehofer Date: Fri, 3 Jul 2026 11:11:43 +0200 Subject: [PATCH 06/31] =?UTF-8?q?ops(cluster):=20campaign-shape=20job=20te?= =?UTF-8?q?mplates=20=E2=80=94=20walltimes,=20private=20CWDs,=20safe=20res?= =?UTF-8?q?ubmit,=20provenance=20probes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Readiness-sweep fixes TC-01..TC-15 (adversarially verified where campaign-blocking; measured anchors 2026-07-03, jobs 5732036/5735964/5735965): - evaluate.sbatch: 128 cpu/15 min -> 16 cpu/6 h placeholder (anchor: 56-76 min per h-value @ 3355 events / 16 cpus; the old pair guaranteed 38/38 TIMEOUTs + a dead chained combine); partition cpu,cpu_il; 41-value hybrid h-grid (0.655/0.665/0.675 added so the 0.005-dense window covers the 0.67 closure truth); per-task --seed (EVAL_SEED + task id) + --simulation_index metadata; per-h skip-if-output idempotency guard. - combine.sbatch: 45 -> 90 min; cluster-side figures OFF by default (RENDER_FIGURES=1 to opt in) — the figures phase is what killed job 5735965; posteriors-only anchor ~20 min. - ALL array templates (simulate/evaluate/combine/inject): private per-run CWD (/cwd with simulations + master_thesis_code symlinks) replaces the shared /simulations symlink dance — concurrent multi-seed pipelines could silently cross-contaminate CRB/ posterior writes (TC-03); merge.sbatch already --workdir-safe. - submit_pipeline.sh: --array derived from the evaluate.sbatch H_VALUES line (single grid source, TC-14); --h_true flag (threads H_VALUE, suffixes RUN_DIR) and REQUIRED --injection_pool flag (stages the pool into RUN_DIR/simulations/injections; loud failure on empty dir); submit-time posterior archiving (replaces the racy task-0 in-job archive, TC-06). - resubmit_failed.sh: TIMEOUT removed from the default state filter (gpu_h100_short tasks TIMEOUT by design — the old filter would delete ALL outputs of a healthy run); h_value recovered from run_metadata before cleanup with conflict-abort (closure-truth contamination fix); abort if the merged CRB already exists (append-duplication guard). - preflight.sh + cluster.env: O(1) catalogue fingerprint (row-1 z field: 0.001733 = z_cmb, 0.00099... = STALE z_helio — both are 8-col, so the column count cannot discriminate) + shallow-injection-pool depth probe. - datasets.yaml: seed43000_Mz/seed700 RETIRED (pre-dt2, z_cut=0.5); depth15_campaign placeholder pending regeneration. - submit_resimulate_phase50.sh: stale INJECTION_SOURCE default removed (must be set explicitly). - LAUNCHING_JOBS.md / README.md: private-CWD mental model, 2026-07-03 walltime anchors, new flags + resubmit signature. bash -n clean on all 10 shell files; 41-label convention verified against bayesian_statistics posterior filenames. Co-Authored-By: Claude Fable 5 --- cluster/LAUNCHING_JOBS.md | 56 ++++++--- cluster/README.md | 24 +++- cluster/cluster.env | 7 ++ cluster/combine.sbatch | 37 ++++-- cluster/datasets.yaml | 16 ++- cluster/evaluate.sbatch | 115 ++++++++--------- cluster/inject.sbatch | 14 ++- cluster/merge.sbatch | 2 +- cluster/preflight.sh | 32 +++++ cluster/resubmit_failed.sh | 155 +++++++++++++++++++---- cluster/simulate.sbatch | 14 ++- cluster/submit_pipeline.sh | 177 +++++++++++++++++++++++---- cluster/submit_resimulate_phase50.sh | 25 ++-- 13 files changed, 510 insertions(+), 164 deletions(-) diff --git a/cluster/LAUNCHING_JOBS.md b/cluster/LAUNCHING_JOBS.md index cf10368b..3cbe8627 100644 --- a/cluster/LAUNCHING_JOBS.md +++ b/cluster/LAUNCHING_JOBS.md @@ -16,15 +16,17 @@ There are **two distinct directories**, and confusing them causes most failures: | **RUN_DIR** | this run's **output** (logs, CSVs, posteriors) | `$WORKSPACE/run_YYYYMMDD_seedS/` | The package uses **relative paths from the current working directory**: -- it reads the catalog from `./master_thesis_code/galaxy_catalogue/`, +- it reads the catalog from `./master_thesis_code/galaxy_catalogue/` (handler.py:24), - it reads/writes `./simulations/…`. -So every job does the same dance: **`cd $PROJECT_ROOT`** then -**`ln -sfn $RUN_DIR/simulations $PROJECT_ROOT/simulations`** — run from the code, -but redirect `./simulations` to this run's output. Output therefore lands in -`$RUN_DIR/simulations/`. (Because there is one shared symlink, avoid running two -jobs from the same PROJECT_ROOT with different RUN_DIRs *interactively* at once; -SLURM tasks each re-point it at start, and all point to the same RUN_DIR per job.) +So every batch job runs from a **private per-run CWD** (TC-03): it `cd`s into +`$RUN_DIR/cwd/`, which holds two symlinks — +`simulations → $RUN_DIR/simulations` and +`master_thesis_code → $PROJECT_ROOT/master_thesis_code`. Code and catalog come +from the one repo; output lands in `$RUN_DIR/simulations/`. Because each run +owns its CWD, **concurrent runs with different RUN_DIRs are safe** — there is +no shared `$PROJECT_ROOT/simulations` symlink to fight over anymore. +(`merge.sbatch` needs no CWD tricks — it uses absolute `--workdir` paths.) **Env threading:** submit wrappers pass everything the sbatch needs via `sbatch --export=ALL,RUN_DIR=…,BASE_SEED=…`. The sbatch validates them and falls @@ -38,11 +40,11 @@ job `source cluster/modules.sh` (which also exports `$WORKSPACE`, `$PROJECT_ROOT ## 2. Partition cheat-sheet -| Partition | Use | Limits | +| Partition | Use | Limits / anchors (2026-07-03) | |---|---|---| | `gpu_h100_short` | production GPU sim (tasks are time-capped, backfills fast) | 30-min wall, 1 GPU/task | -| `gpu_a100_short` | GPU smoke tests | 5-min wall | -| `cpu_il` | inference / merge / combine | up to 128 cpus, ~15 min/h-value | +| `gpu_a100_short` | injection campaigns (`inject.sbatch`) + GPU smoke tests | `inject.sbatch` requests a 30-min wall; the smoke test uses 5 min | +| `cpu,cpu_il` | inference / merge / combine | evaluate: **56–76 min per h-value @ 3355 events / 16 cpus** (jobs 5732036, volume_deconv; 6h pre-smoke budget, re-size after smoke); combine: **~20 min posteriors** + 90-min budget (job 5735965 anchor); figures rendered locally (`RENDER_FIGURES=0`) | | `dev_gpu_h100` / `dev_*` | quick queue for testing | short wall, fast start | Seed convention everywhere: **per-task seed = `BASE_SEED + SLURM_ARRAY_TASK_ID`** @@ -57,12 +59,19 @@ Run all of these **from `~/MasterThesisCode` after `source cluster/modules.sh`.* ### 3a. Simulation → CRB (+ auto merge → evaluate → combine) One command chains simulate (GPU) → merge (CPU) → evaluate (CPU) → combine: ```bash -bash cluster/submit_pipeline.sh --tasks 100 --steps 50 --seed 42 +bash cluster/submit_pipeline.sh --tasks 100 --steps 50 --seed 42 \ + --injection_pool "$WORKSPACE/injection__seed/simulations/injections" # creates $WORKSPACE/run_YYYYMMDD_seed42/ ; prints all job IDs + a sacct line ``` - `--tasks` = GPU array size, `--steps` = EMRI iterations/task, `--seed` = base seed. +- `--injection_pool` (required unless `--no_injections`) links the pool's + `injection_h_*.csv` into `RUN_DIR/simulations/injections/` at submit time, so + evaluate's p_det grid uses exactly the intended pool (see `cluster/datasets.yaml`). +- `--h_true V` sets the injected truth for closure runs (default 0.73); a + non-default truth is embedded in the run-dir name (`run_YYYYMMDD_seedS_h0p67`). - Dependency chain: simulate → merge (`afterany`, tolerates task timeouts) → - evaluate (`afterok`, 38-point h-grid 0.60–0.86) → combine (`afterok`). + evaluate (`afterok`, h-grid parsed from `evaluate.sbatch` — currently 41 + points 0.60–0.86) → combine (`afterok`). - **Test small first:** `--tasks 2 --steps 10`. ### 3b. Injection campaign → P_det pool @@ -79,7 +88,8 @@ bash cluster/submit_injection.sh --tasks_per_h 80 --steps 900 --seed 12345 Usually part of 3a, but to (re-)evaluate an existing run: ```bash RUN=$WORKSPACE/run_20260516_seed400_phase50 -sbatch --parsable --array=0-37 \ +# --array must match the H_VALUES count in evaluate.sbatch (currently 41 → 0-40) +sbatch --parsable --array=0-40 \ --output="$RUN/logs/evaluate_%A_%a.out" --error="$RUN/logs/evaluate_%A_%a.err" \ --export=ALL,RUN_DIR="$RUN" cluster/evaluate.sbatch # then combine: @@ -166,13 +176,25 @@ inline wrapper that makes `RUN_DIR` and passes `--export`. ## 6. Re-run safety & idempotency - **Skip-if-output** per unit (evaluate/combine already do this): guard on the - target file so resubmits don't redo finished work. -- **Archive-then-write**: `evaluate.sbatch` task 0 archives existing - `posteriors*/` to `simulations/archive/eval_/` before a fresh sweep — so a + target file so resubmits don't redo finished work. `evaluate.sbatch` exits 0 + per-task if its `h_