From 6622fd5775bac0c91f4e8cda07b5160f17ba1ae1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eloy=20Pe=C3=B1a-Asensio?= Date: Wed, 22 Jul 2026 21:20:04 +0200 Subject: [PATCH] Fix silent likelihood/reporting corruption in the Dynesty log-likelihood pipeline Four bugs in DynestyMetSim.py let a bad model draw leak a finite, misleading value into the nested-sampling likelihood surface (or crash the results summary) instead of being rejected/surfaced correctly: 1. np.interp() in logLikelihoodDynesty (simulated_time, simulated_lc_intensity, simulated_lag) clamped out-of-range height queries to the simulation's edge value instead of returning NaN. A simulation that terminates before covering the full observed height range (e.g. the meteoroid fully ablates or hits h_kill/v_kill early) therefore slipped past the NaN-count guard that was supposed to reject it, and got scored against a duplicated final value instead of being rejected. Fixed by passing left=np.nan, right=np.nan, matching the pattern already used correctly elsewhere in the file (_worker_simulate_and_interp). 2. The frame-integration branch of the luminosity comparison (the common case, since camera frame time is normally longer than the simulation step) never checked that integrateLuminosity()'s output covered the same points as the plain-interpolation branch did. Added the missing NaN-count guard to that branch. 3. runSimulationDynesty caught ZeroDivisionError from the simulation and silently retried with a brand new, unrelated nominal Constants() object, discarding the sampled parameters entirely. dynesty would then score that live point using a completely different model without knowing it. Fixed by letting the exception propagate and having logLikelihoodDynesty reject the point with -inf, the same way it already handles a timeout. 4. summaryResultsTable's marginal-mode histogram (Mode_{Ndim} column) ran on the raw, NaN-containing samples/weights instead of the already-masked x_valid/w_valid used for every other statistic in the same loop. np.min/ np.max of a NaN-containing array is NaN, and np.histogram raises ValueError on a NaN range, so this crashed the entire results/plotting stage whenever any parameter had a NaN sample. The same histogram-plus-mode logic was duplicated (correctly, with masking) in the posterior distribution plot code a few hundred lines down, which is how the divergence was caught. Extracted a shared _weightedHistogramMode() helper used by both call sites, so there's now a single place responsible for this computation. Test plan: - Added wmpl/Dynesty/Tests/test_DynestyMetSim.py: 6 regression tests covering all four bugs. runSimulationDynesty/runSimulation/constructConstants are mocked out with hand-built SimpleNamespace stand-ins, so the suite runs in milliseconds with no real MetSim/dynesty evaluation. - Verified each test fails against the pre-fix code (via git stash) and passes with the fixes applied. - Run with: python -m wmpl.Dynesty.Tests.test_DynestyMetSim --- wmpl/Dynesty/DynestyMetSim.py | 92 +++++++------ wmpl/Dynesty/Tests/test_DynestyMetSim.py | 168 +++++++++++++++++++++++ 2 files changed, 221 insertions(+), 39 deletions(-) create mode 100644 wmpl/Dynesty/Tests/test_DynestyMetSim.py diff --git a/wmpl/Dynesty/DynestyMetSim.py b/wmpl/Dynesty/DynestyMetSim.py index 5f6ffbd8..6d28604f 100644 --- a/wmpl/Dynesty/DynestyMetSim.py +++ b/wmpl/Dynesty/DynestyMetSim.py @@ -937,6 +937,30 @@ def _quantiles_from_samples(arr_2d, qs): out[name] = np.nanquantile(arr_2d, q, axis=0, method='linear') return out +def _weightedHistogramMode(x_valid, w_valid, smooth): + """ Weighted, optionally dynesty-style-smoothed histogram of already NaN-filtered samples. + + Arguments: + x_valid: [ndarray] 1D sample values, with any NaNs already masked out by the caller. + w_valid: [ndarray] Importance weights matching x_valid (need not be renormalized). + smooth: [float or int] Number of bins (int) or smoothing fraction used to derive + nbins = round(10/smooth) and apply a Gaussian KDE-style smoothing (float). + + Return: + centers: [ndarray] Bin centers. + hist: [ndarray] (Smoothed) weighted histogram counts, same length as centers. + + """ + lo, hi = np.min(x_valid), np.max(x_valid) + if isinstance(smooth, int): + hist, edges = np.histogram(x_valid, bins=smooth, weights=w_valid, range=(lo, hi)) + else: + nbins = int(round(10./smooth)) + hist, edges = np.histogram(x_valid, bins=nbins, weights=w_valid, range=(lo, hi)) + hist = norm_kde(hist, 10.0) + centers = 0.5*(edges[1:] + edges[:-1]) + return centers, hist + def _maybe_integrate_luminosity(sim, obs_data): """ If (1/fps_lum) > sim.const.dt, integrate luminosity over fps window and @@ -2088,14 +2112,7 @@ def summaryResultsTable(results, mode_value = mode_raw[i] # mode via corner logic - lo, hi = np.min(x), np.max(x) - if isinstance(smooth, int): - hist, edges = np.histogram(x, bins=smooth, weights=w, range=(lo,hi)) - else: - nbins = int(round(10./smooth)) - hist, edges = np.histogram(x, bins=nbins, weights=w, range=(lo,hi)) - hist = norm_kde(hist, 10.0) - centers = 0.5*(edges[1:] + edges[:-1]) + centers, hist = _weightedHistogramMode(x_valid, w_valid, smooth) mode_Ndim = centers[np.argmax(hist)] # now apply your log & unit transforms *after* computing stats @@ -2552,15 +2569,7 @@ def log_transf(v): continue # Compute histogram - lo, hi = np.min(x_valid), np.max(x_valid) - if isinstance(smooth, int): - hist, edges = np.histogram(x_valid, bins=smooth, weights=w_valid, range=(lo, hi)) - else: - nbins = int(round(10./smooth)) - hist, edges = np.histogram(x_valid, bins=nbins, weights=w_valid, range=(lo, hi)) - hist = norm_kde(hist, 10.0) # dynesty-style smoothing - - centers = 0.5*(edges[1:] + edges[:-1]) + centers, hist = _weightedHistogramMode(x_valid, w_valid, smooth) # Fill under the curve ax.fill_between(centers, hist, color='blue', alpha=0.6) @@ -5858,20 +5867,14 @@ def runSimulationDynesty(parameter_guess, real_event, var_names, fix_var): """ - # build the const to run the + # build the const to run the const_nominal = constructConstants(parameter_guess, real_event, var_names, fix_var) - try: - # Run the simulation - frag_main, results_list, wake_results = runSimulation(const_nominal, compute_wake=False) - simulation_MetSim_object = SimulationResults(const_nominal, frag_main, results_list, wake_results) - except ZeroDivisionError as e: - print(f"Error during simulation: {e}") - # run again with the nominal values to avoid the error - const_nominal = Constants() - # Run the simulation - frag_main, results_list, wake_results = runSimulation(const_nominal, compute_wake=False) - simulation_MetSim_object = SimulationResults(const_nominal, frag_main, results_list, wake_results) + # Run the simulation. A ZeroDivisionError means this parameter draw is unphysical; + # let it propagate so the caller (logLikelihoodDynesty) can reject the point with + # -np.inf instead of silently swapping in an unrelated nominal simulation. + frag_main, results_list, wake_results = runSimulation(const_nominal, compute_wake=False) + simulation_MetSim_object = SimulationResults(const_nominal, frag_main, results_list, wake_results) return simulation_MetSim_object @@ -6112,7 +6115,10 @@ def logLikelihoodDynesty(guess_var, obs_metsim_obj, flags_dict, fix_var, timeout # check if the OS is not Linux if os.name != 'posix': # If not Linux, run the simulation without timeout - simulation_results = runSimulationDynesty(guess_var, obs_metsim_obj, var_names, fix_var) + try: + simulation_results = runSimulationDynesty(guess_var, obs_metsim_obj, var_names, fix_var) + except ZeroDivisionError: + return -np.inf # unphysical parameter draw, reject the point else: # Set timeout handler signal.signal(signal.SIGALRM, timeout_handler) @@ -6123,14 +6129,17 @@ def logLikelihoodDynesty(guess_var, obs_metsim_obj, flags_dict, fix_var, timeout except TimeoutException: print('timeout') return -np.inf # immediately return -np.inf if times out + except ZeroDivisionError: + return -np.inf # unphysical parameter draw, reject the point finally: signal.alarm(0) # Cancel alarm ### LUM CALC ### - simulated_time = np.interp(obs_metsim_obj.height_lum, - np.flip(simulation_results.leading_frag_height_arr), - np.flip(simulation_results.time_arr)) + simulated_time = np.interp(obs_metsim_obj.height_lum, + np.flip(simulation_results.leading_frag_height_arr), + np.flip(simulation_results.time_arr), + left=np.nan, right=np.nan) # check if the length of the lag_sim is the same as the length of the obs_metsim_obj.lag if np.sum(~np.isnan(simulated_time)) != np.sum(~np.isnan(obs_metsim_obj.time_lum)): return -np.inf @@ -6140,11 +6149,15 @@ def logLikelihoodDynesty(guess_var, obs_metsim_obj, flags_dict, fix_var, timeout # find the integral of the luminosity in time in between FPS but not valid for CAMO narrowfield cameras as there is no smearing becuse it follows the meteor if (1/obs_metsim_obj.fps_lum > simulation_results.const.dt): # and (not any('1T' in station for station in obs_metsim_obj.stations_lum) or not any('2T' in station for station in obs_metsim_obj.stations_lum)): # FPS is lower than the simulation time step need to integrate the luminosity simulated_lc_intensity, _ = integrateLuminosity(all_simulated_time,obs_metsim_obj.time_lum,simulation_results.luminosity_arr,simulation_results.const.dt,obs_metsim_obj.fps_lum,obs_metsim_obj.P_0m) + # check if the length of the simulated_lc_intensity is the same as the length of the obs_metsim_obj.luminosity + if np.sum(~np.isnan(simulated_lc_intensity)) != np.sum(~np.isnan(obs_metsim_obj.luminosity)): + return -np.inf else: # too high frame rate, just interpolate the luminosity - simulated_lc_intensity = np.interp(obs_metsim_obj.height_lum, - np.flip(simulation_results.leading_frag_height_arr), - np.flip(simulation_results.luminosity_arr)) + simulated_lc_intensity = np.interp(obs_metsim_obj.height_lum, + np.flip(simulation_results.leading_frag_height_arr), + np.flip(simulation_results.luminosity_arr), + left=np.nan, right=np.nan) # check if the length of the simulated_lc_intensity is the same as the length of the obs_metsim_obj.luminosity if np.sum(~np.isnan(simulated_lc_intensity)) != np.sum(~np.isnan(obs_metsim_obj.luminosity)): return -np.inf @@ -6153,9 +6166,10 @@ def logLikelihoodDynesty(guess_var, obs_metsim_obj, flags_dict, fix_var, timeout lag_sim = simulation_results.leading_frag_length_arr - (obs_metsim_obj.v_init*simulation_results.time_arr) - simulated_lag = np.interp(obs_metsim_obj.height_lag, - np.flip(simulation_results.leading_frag_height_arr), - np.flip(lag_sim)) + simulated_lag = np.interp(obs_metsim_obj.height_lag, + np.flip(simulation_results.leading_frag_height_arr), + np.flip(lag_sim), + left=np.nan, right=np.nan) lag_sim = simulated_lag - simulated_lag[0] diff --git a/wmpl/Dynesty/Tests/test_DynestyMetSim.py b/wmpl/Dynesty/Tests/test_DynestyMetSim.py new file mode 100644 index 00000000..c5ce5d6b --- /dev/null +++ b/wmpl/Dynesty/Tests/test_DynestyMetSim.py @@ -0,0 +1,168 @@ +""" Regression tests for the three logLikelihoodDynesty/runSimulationDynesty bugs found in +DynestyMetSim.py, all of which let a bad model draw leak a finite log-likelihood into the +nested-sampling likelihood surface instead of being rejected with -inf: + + - out-of-range np.interp queries used to clamp to the simulation's edge value instead of + returning NaN, so a simulation that terminates before covering the observed height range + slipped past the NaN-count guard that was supposed to reject it + - the frame-integration branch of the luminosity comparison (the common case, since camera + frame time is normally longer than the simulation step) never had that NaN-count guard at + all - only the plain-interpolation branch did + - a ZeroDivisionError during the simulation used to be swallowed and silently replaced with an + unrelated nominal Constants() simulation, so dynesty would score a completely different model + at that live point without any indication it happened + - the marginal-mode histogram in the results summary table (and, near-identically, in the + posterior distribution plot) ran on the raw, NaN-containing samples/weights instead of the + already NaN-masked x_valid/w_valid used for every other statistic, so np.min/np.max silently + became NaN whenever a parameter had any NaN sample - which np.histogram then turns into a hard + ValueError, crashing the whole results/plotting stage. Both call sites now share + _weightedHistogramMode(), which only ever operates on pre-masked input. + +Every simulation result here is a hand-built stand-in (SimpleNamespace) with runSimulationDynesty +patched out, so these run in milliseconds instead of paying for a real MetSim/dynesty evaluation. + +Run under pytest, or directly: + + python -m wmpl.Dynesty.Tests.test_DynestyMetSim +""" + +import warnings +from types import SimpleNamespace +from unittest.mock import patch + +import numpy as np + +import wmpl.Dynesty.DynestyMetSim as DynestyMetSim + + +def _makeObs(**overrides): + """ Minimal ObservationData-like stand-in with just the attributes logLikelihoodDynesty reads. """ + obs = SimpleNamespace( + height_lum=np.array([98000.0, 90000.0, 82000.0, 75000.0, 70000.0]), + time_lum=np.array([0.1, 1.0, 2.0, 3.0, 4.0]), + luminosity=np.array([100.0, 90.0, 80.0, 70.0, 60.0]), + fps_lum=30.0, + P_0m=840.0, + height_lag=np.array([98000.0, 90000.0, 82000.0]), + lag=np.array([0.0, 1.0, 2.0]), + v_init=16000.0, + noise_lum=1.0, + noise_lag=1.0, + ) + for key, value in overrides.items(): + setattr(obs, key, value) + return obs + + +def _makeSim(height_max=100000.0, height_min=80000.0, n=5, dt=1.0): + """ Minimal SimulationResults-like stand-in: a straight descent from height_max to height_min, + sampled every dt seconds (dt=1.0 by default since most tests never reach the integration branch + that actually depends on realistic time spacing). """ + height = np.linspace(height_max, height_min, n) # decreasing, like a real leading_frag_height_arr + time = np.arange(n)*dt + return SimpleNamespace( + leading_frag_height_arr=height, + time_arr=time, + leading_frag_length_arr=time*1000.0, + luminosity_arr=np.linspace(50.0, 150.0, n), + const=SimpleNamespace(dt=dt), + ) + + +def test_truncated_simulation_is_rejected(): + """ Bug 1: obs heights below the simulation's lowest reached height (75000/70000, while the sim + stops at 80000) must make the NaN-count guard fire, since np.interp is no longer allowed to clamp + them to the simulation's edge time/value. """ + sim = _makeSim(height_min=80000.0) + obs = _makeObs() + + with patch.object(DynestyMetSim, "runSimulationDynesty", return_value=sim): + result = DynestyMetSim.logLikelihoodDynesty([1.0], obs, {"v_init": []}, {}) + + assert result == -np.inf, "truncated simulation was not rejected (np.interp clamping regression)" + + +def test_integration_branch_checks_coverage(): + """ Bug 2: force the frame-integration branch (1/fps_lum > dt) with a realistic, dt-spaced time + grid, and make the last observed frame time (100.0) fall far outside every simulated time sample. + integrateLuminosity() returns NaN there while obs.luminosity has a real value at that index, so + the branch must now reject the draw instead of silently nansum-ing over the mismatch. """ + sim = _makeSim(height_max=100000.0, height_min=80000.0, n=500, dt=0.01) # 1/30 s frame >> 0.01 s dt + obs = _makeObs( + height_lum=np.array([98000.0, 95000.0, 92000.0, 90000.0]), + time_lum=np.array([0.5, 1.5, 2.5, 100.0]), # last point far beyond any simulated time + luminosity=np.array([100.0, 90.0, 80.0, 70.0]), + ) + + with patch.object(DynestyMetSim, "runSimulationDynesty", return_value=sim), \ + warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) # expected empty-slice mean at i=3 + result = DynestyMetSim.logLikelihoodDynesty([1.0], obs, {"v_init": []}, {}) + + assert result == -np.inf, "integration branch let a coverage mismatch through" + + +def test_zero_division_rejected_not_swallowed(): + """ Bug 3: a ZeroDivisionError raised while running the simulation must reject the point (-inf), + not silently fall back to an unrelated nominal Constants() simulation. """ + obs = _makeObs() + + with patch.object(DynestyMetSim, "runSimulationDynesty", side_effect=ZeroDivisionError("boom")): + result = DynestyMetSim.logLikelihoodDynesty([1.0], obs, {"v_init": []}, {}) + + assert result == -np.inf, "ZeroDivisionError was not converted to a -inf rejection" + + +def test_run_simulation_dynesty_propagates_zero_division(): + """ Bug 3, lower level: runSimulationDynesty itself must not swallow ZeroDivisionError into a + fallback nominal simulation - it must propagate so the caller can decide how to reject the draw. """ + dummy_const = SimpleNamespace() + with patch.object(DynestyMetSim, "constructConstants", return_value=dummy_const), \ + patch.object(DynestyMetSim, "runSimulation", side_effect=ZeroDivisionError("boom")): + try: + DynestyMetSim.runSimulationDynesty([1.0], object(), ["v_init"], {}) + except ZeroDivisionError: + pass + else: + raise AssertionError("runSimulationDynesty swallowed ZeroDivisionError into a fallback") + + +def test_weighted_histogram_mode_recovers_true_peak(): + """ Bug 4: _weightedHistogramMode (shared by summaryResultsTable's Mode_{Ndim} column and the + posterior distribution plot) must reproduce a sane weighted mode from pre-masked samples. """ + rng = np.random.default_rng(0) + x_valid = rng.normal(loc=5.0, scale=1.0, size=5000) + w_valid = np.full(x_valid.size, 1.0/x_valid.size) + + centers, hist = DynestyMetSim._weightedHistogramMode(x_valid, w_valid, smooth=0.02) + mode = centers[np.argmax(hist)] + + assert np.isfinite(mode), "mode is not finite on clean, pre-masked input" + assert abs(mode - 5.0) < 0.5, "recovered mode is far from the true peak at 5.0" + + +def test_unmasked_nan_input_raises_documenting_why_callers_must_mask(): + """ Bug 4 root cause, pinned down: feeding _weightedHistogramMode raw (unmasked) samples containing + a NaN reproduces the original crash (np.min/np.max of a NaN-containing array is NaN, and + np.histogram rejects a NaN range outright). This is why both call sites now mask into + x_valid/w_valid before calling the helper - if this assertion ever stops raising, np.histogram's + behavior changed and the masking discipline in the two callers should be re-checked. """ + x = np.array([1.0, 2.0, np.nan, 3.0]) + w = np.full(4, 0.25) + + try: + DynestyMetSim._weightedHistogramMode(x, w, smooth=0.02) + except ValueError: + pass + else: + raise AssertionError("expected ValueError from an unmasked NaN range - masking may be needed elsewhere too") + + +if __name__ == "__main__": + test_truncated_simulation_is_rejected() + test_integration_branch_checks_coverage() + test_zero_division_rejected_not_swallowed() + test_run_simulation_dynesty_propagates_zero_division() + test_weighted_histogram_mode_recovers_true_peak() + test_unmasked_nan_input_raises_documenting_why_callers_must_mask() + print("All DynestyMetSim log-likelihood regression checks passed.")