Merge pull request #53 from CCPBioSim/7-nan-error-handling

Add error handling for NaN and invalid values in calculations:

Author: harryswift01

CodeEntropy/EntropyFunctions.py

@@ -40,6 +40,14 @@ def frequency_calculation(lambdas, temp):
     pi = np.pi
     # get kT in Joules from given temperature
     kT = UAC.get_KT2J(temp)
+
+    lambdas = np.array(lambdas)  # Ensure input is a NumPy array
+
+    # Check for negatives and raise an error if any are found
+    if np.any(lambdas < 0):
+        raise ValueError(f"Negative eigenvalues encountered: {lambdas[lambdas < 0]}")
+
+    # Compute frequencies safely
     frequencies = 1 / (2 * pi) * np.sqrt(lambdas / kT)
 
     return frequencies

CodeEntropy/GeometricFunctions.py

@@ -181,6 +181,26 @@ def get_sphCoord_axes(arg_r):
     x2y2 = arg_r[0] ** 2 + arg_r[1] ** 2
     r2 = x2y2 + arg_r[2] ** 2
 
+    # Check for division by zero
+    if r2 == 0.0:
+        raise ValueError("r2 is zero, cannot compute spherical coordinates.")
+
+    if x2y2 == 0.0:
+        raise ValueError("x2y2 is zero, cannot compute sin_phi and cos_phi.")
+
+    # Check for non-negative values inside the square root
+    if x2y2 / r2 < 0:
+        raise ValueError(
+            f"Negative value encountered for sin_theta calculation: {x2y2 / r2}. "
+            f"Cannot take square root."
+        )
+
+    if x2y2 < 0:
+        raise ValueError(
+            f"Negative value encountered for sin_phi and cos_phi calculation: {x2y2}. "
+            f"Cannot take square root."
+        )
+
     if x2y2 != 0.0:
         sin_theta = np.sqrt(x2y2 / r2)
         cos_theta = arg_r[2] / np.sqrt(r2)
@@ -259,6 +279,13 @@ def get_weighted_forces(
     # divide the sum of forces by the mass of the bead to get the weighted forces
     mass = bead.total_mass()
 
+    # Check that mass is positive to avoid division by 0 or negative values inside sqrt
+    if mass <= 0:
+        raise ValueError(
+            f"Invalid mass value: {mass}. Mass must be positive to compute the square "
+            f"root."
+        )
+
     weighted_force = forces_trans / np.sqrt(mass)
 
     return weighted_force
@@ -309,14 +336,18 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5)
     moment_of_inertia = bead.moment_of_inertia()
 
     for dimension in range(3):
-        # cannot divide by zero
-        if np.isclose(moment_of_inertia[dimension, dimension], 0):
-            weighted_torque[dimension] = torques[dimension]
-        else:
-            weighted_torque[dimension] = torques[dimension] / np.sqrt(
-                moment_of_inertia[dimension, dimension]
+        # Check if the moment of inertia is valid for square root calculation
+        inertia_value = moment_of_inertia[dimension, dimension]
+
+        if np.isclose(inertia_value, 0):
+            raise ValueError(
+                f"Invalid moment of inertia value: {inertia_value}. "
+                f"Cannot compute weighted torque."
             )
 
+        # compute weighted torque if moment of inertia is valid
+        weighted_torque[dimension] = torques[dimension] / np.sqrt(inertia_value)
+
     return weighted_torque
 
 

CodeEntropy/poseidon/extractData/generalFunctions.py

@@ -18,6 +18,16 @@ def distance(x0, x1, dimensions):
     x1 = np.array(x1)
     delta = np.abs(x1 - x0)
     delta = np.where(delta > 0.5 * dimensions, delta - dimensions, delta)
+
+    if np.any(delta < 0):
+        negative_indices = np.where(delta < 0)
+        negative_values = delta[negative_indices]
+        raise ValueError(
+            f"Negative values encountered in 'delta' at indices {negative_indices} "
+            f"with values {negative_values}. "
+            f"Cannot take square root of negative values."
+        )
+
     dist = np.sqrt((delta**2).sum(axis=-1))
     return dist
 
@@ -173,8 +183,38 @@ def angle(a, b, c, dimensions):
     ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba)
     bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc)
     ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac)
+
+    if np.any(ba < 0):
+        negative_indices = np.where(ba < 0)
+        negative_values = ba[negative_indices]
+        raise ValueError(
+            f"Negative values encountered in 'ba' at indices {negative_indices} "
+            f"with values {negative_values}. "
+            f"Cannot take square root of negative values."
+        )
+
     dist_ba = np.sqrt((ba**2).sum(axis=-1))
+
+    if np.any(bc < 0):
+        negative_indices = np.where(bc < 0)
+        negative_values = bc[negative_indices]
+        raise ValueError(
+            f"Negative values encountered in 'bc' at indices {negative_indices} "
+            f"with values {negative_values}. "
+            f"Cannot take square root of negative values."
+        )
+
     dist_bc = np.sqrt((bc**2).sum(axis=-1))
+
+    if np.any(ac < 0):
+        negative_indices = np.where(ac < 0)
+        negative_values = ac[negative_indices]
+        raise ValueError(
+            f"Negative values encountered in 'ac' at indices {negative_indices} "
+            f"with values {negative_values}. "
+            f"Cannot take square root of negative values."
+        )
+
     dist_ac = np.sqrt((ac**2).sum(axis=-1))
 
     cosine_angle = (dist_ac**2 - dist_bc**2 - dist_ba**2) / (-2 * dist_bc * dist_ba)