From 9b147cc03fb3627ae6402c14f46826a9f73327b7 Mon Sep 17 00:00:00 2001
From: JeffreyCovington <22208551+JeffreyCovington@users.noreply.github.com>
Date: Wed, 26 Feb 2025 18:06:40 -0700
Subject: [PATCH] Fixed how the particle filter handles anomolous
 floating-point calculations. (#238)

---
 epymorph/parameter_fitting/utils/resampler.py | 25 +++++++++++++------
 1 file changed, 17 insertions(+), 8 deletions(-)

diff --git a/epymorph/parameter_fitting/utils/resampler.py b/epymorph/parameter_fitting/utils/resampler.py
index 94842365..ae2081d4 100644
--- a/epymorph/parameter_fitting/utils/resampler.py
+++ b/epymorph/parameter_fitting/utils/resampler.py
@@ -62,7 +62,7 @@ def compute_weights(
         for i in range(self.N):
             expected = expected_observations[i]
             try:
-                log_weights[i] = np.sum(
+                log_weights[i] = np.nansum(
                     np.log(self.likelihood_fn.compute(current_obs_data, expected))
                 )
             except (ValueError, FloatingPointError):
@@ -170,18 +170,27 @@ def compute_weights(
         """
         n_nodes = expected_observations[0].shape[0]
         log_weights = np.zeros(shape=(n_nodes, self.N))
+        shifted_log_weights = np.zeros(shape=(n_nodes, self.N))
 
         for i_node in range(n_nodes):
-            for i_particle in range(self.N):
-                expected = expected_observations[i_particle][i_node]
+            if not np.isnan(current_obs_data[i_node]):
+                for i_particle in range(self.N):
+                    expected = expected_observations[i_particle][i_node]
+                    try:
+                        log_weights[i_node, i_particle] = np.log(
+                            self.likelihood_fn.compute(
+                                current_obs_data[i_node], expected
+                            )
+                        )
+                    except (ValueError, FloatingPointError):
+                        log_weights[i_node, i_particle] = -np.inf
                 try:
-                    log_weights[i_node, i_particle] = np.log(
-                        self.likelihood_fn.compute(current_obs_data[i_node], expected)
+                    shifted_log_weights[i_node, :] = log_weights[i_node, :] - np.max(
+                        log_weights[i_node, :], keepdims=True
                     )
-                except (ValueError, FloatingPointError):
-                    log_weights[i_node, i_particle] = -np.inf
+                except FloatingPointError:
+                    shifted_log_weights[i_node, :] = 0
 
-        shifted_log_weights = log_weights - np.max(log_weights, axis=1, keepdims=True)
         underflow_lower_bound = -(10**2)
         clipped_log_weights = np.clip(
             shifted_log_weights, a_min=underflow_lower_bound, a_max=None