Fix silent likelihood/reporting corruption in the Dynesty log-likelihood pipeline

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]
 

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