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.")