Safety Refactor: Climate, Game Theory, and Radar Gating

math_explorer/src/applied/algorithms/sorting/linear.rs

@@ -15,7 +15,7 @@ impl Sorter<u64> for RadixSorter {
         }
 
         // Safe unwrap: sorted_data is not empty (checked above)
-        let max_val = *sorted_data.iter().max().expect("sorted_data is not empty");
+        let max_val = *sorted_data.iter().max().unwrap_or(&0);
         // Comparisons to find max? technically yes, but usually considered O(n) overhead not part of sort loop logic strictly.
         // We'll count them for rigorousness.
         stats.comparisons += sorted_data.len() as u64 - 1;

math_explorer/src/applied/game_theory/error.rs

@@ -7,6 +7,8 @@ pub enum GameTheoryError {
     NonSquarePayoffMatrix { rows: usize, cols: usize },
     /// Invalid parameter value (e.g. frequency outside [0, 1]).
     InvalidParameter { name: String, value: f64 },
+    /// Population has gone extinct (sum close to zero), preventing normalization.
+    PopulationExtinction,
 }
 
 impl fmt::Display for GameTheoryError {
@@ -20,6 +22,9 @@ impl fmt::Display for GameTheoryError {
             Self::InvalidParameter { name, value } => {
                 write!(f, "Invalid parameter {}: {}", name, value)
             }
+            Self::PopulationExtinction => {
+                write!(f, "Population extinction: sum of proportions is zero.")
+            }
         }
     }
 }

math_explorer/src/applied/game_theory/evolutionary.rs

@@ -52,7 +52,7 @@ impl ReplicatorDynamics {
         initial_population: DVector<f64>,
         time_horizon: f64,
         dt: f64,
-    ) -> Vec<(f64, DVector<f64>)> {
+    ) -> Result<Vec<(f64, DVector<f64>)>, GameTheoryError> {
         self.simulate_with_strategy(initial_population, time_horizon, dt, &RungeKutta4)
     }
 
@@ -72,7 +72,7 @@ impl ReplicatorDynamics {
         time_horizon: f64,
         dt: f64,
         solver: &S,
-    ) -> Vec<(f64, DVector<f64>)>
+    ) -> Result<Vec<(f64, DVector<f64>)>, GameTheoryError>
     where
         S: Solver<DVector<f64>>,
     {
@@ -89,14 +89,17 @@ impl ReplicatorDynamics {
 
             // Normalize to prevent numerical drift from simplex
             let sum = current_x.sum();
+            if sum.abs() < 1e-12 {
+                return Err(GameTheoryError::PopulationExtinction);
+            }
             if (sum - 1.0).abs() > 1e-9 {
                 current_x /= sum;
             }
 
             trajectory.push((current_t, current_x.clone()));
         }
 
-        trajectory
+        Ok(trajectory)
     }
 }
 
@@ -142,7 +145,7 @@ mod tests {
 
         // Start off-center. Should cycle (closed orbits for zero-sum RPS).
         let init = DVector::from_vec(vec![0.4, 0.3, 0.3]);
-        let trajectory = system.simulate(init, 10.0, 0.01);
+        let trajectory = system.simulate(init, 10.0, 0.01).unwrap();
 
         // Check that we didn't leave the simplex significantly
         let last_state = &trajectory.last().unwrap().1;
@@ -163,7 +166,9 @@ mod tests {
         let init = DVector::from_vec(vec![0.4, 0.3, 0.3]);
 
         // Inject Euler solver
-        let trajectory = system.simulate_with_strategy(init, 5.0, 0.001, &Euler);
+        let trajectory = system
+            .simulate_with_strategy(init, 5.0, 0.001, &Euler)
+            .unwrap();
 
         let last_state = &trajectory.last().unwrap().1;
         assert!((last_state.sum() - 1.0).abs() < 1e-6);

math_explorer/src/climate/loss.rs

@@ -53,8 +53,8 @@ pub fn earth_movers_distance(z1: &DMatrix<f32>, z2: &DMatrix<f32>) -> f32 {
 
         // The EMD for 1D distributions is the L1 norm of the difference
         // between the sorted samples.
-        col1.sort_by(|a, b| a.partial_cmp(b).unwrap());
-        col2.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        col1.sort_by(|a, b| a.total_cmp(b));
+        col2.sort_by(|a, b| a.total_cmp(b));
 
         let emd_i: f32 = col1
             .iter()
@@ -126,4 +126,17 @@ mod tests {
         let emd = earth_movers_distance(&z1, &z2);
         assert!((emd - 4.0).abs() < 1e-6);
     }
+
+    #[test]
+    fn test_earth_movers_distance_nan() {
+        let z1 = DMatrix::from_row_slice(2, 1, &[f32::NAN, 1.0]);
+        let z2 = DMatrix::from_row_slice(2, 1, &[2.0, 1.0]);
+
+        // This should not panic.
+        // total_cmp sorts NaN to the end (or beginning, depending on sign bit, but consistently).
+        // z1 sorted: [1.0, NaN] (likely) or [NaN, 1.0]
+        // z2 sorted: [1.0, 2.0]
+        // It computes distance.
+        let _emd = earth_movers_distance(&z1, &z2);
+    }
 }

math_explorer/src/physics/medical/radar_gating/error.rs

@@ -9,6 +9,8 @@ pub enum RadarError {
     InsufficientSnapshots { required: usize, actual: usize },
     /// Signal subspace dimension error (e.g. >= samples).
     InvalidSignalSubspace { samples: usize, subspace: usize },
+    /// Eigenvalue decomposition failed or produced NaNs.
+    NumericalInstability,
 }
 
 impl fmt::Display for RadarError {
@@ -29,6 +31,9 @@ impl fmt::Display for RadarError {
                 "Signal subspace dimension {} equals or exceeds sample size {}",
                 subspace, samples
             ),
+            Self::NumericalInstability => {
+                write!(f, "Eigenvalue decomposition failed or produced NaNs.")
+            }
         }
     }
 }

math_explorer/src/physics/medical/radar_gating/super_resolution.rs

@@ -121,8 +121,13 @@ impl MusicEstimator {
             .enumerate()
             .map(|(i, &v)| (v, i))
             .collect();
+
+        if pairs.iter().any(|(v, _)| v.is_nan()) {
+            return Err(RadarError::NumericalInstability);
+        }
+
         // Sort descending (largest first)
-        pairs.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal));
+        pairs.sort_by(|a, b| b.0.total_cmp(&a.0));
 
         // 3. Extract Noise Subspace En
         // The last (N - P) eigenvectors correspond to noise.
@@ -201,3 +206,27 @@ impl MusicEstimator {
         Ok(spectrum)
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use nalgebra::Complex;
+
+    #[test]
+    fn test_music_estimator_nan_handling() {
+        let n = 4;
+        let mut estimator = MusicEstimator::new(n, 2, 1);
+
+        // Add valid snapshots
+        let snap1 = vec![Complex::new(1.0, 0.0); n];
+        let snap2 = vec![Complex::new(0.0, 1.0); n];
+        estimator.add_snapshot(&snap1).unwrap();
+        estimator.add_snapshot(&snap2).unwrap();
+
+        // Normally we can't easily force eigen decomposition to produce NaNs without
+        // invalid input or matrix properties.
+        // However, we can assert that with valid input it doesn't return NumericalInstability.
+        let result = estimator.compute_spectrum(0.0, 1.0, 0.1, 1e9, 3e8);
+        assert!(result.is_ok());
+    }
+}

math_explorer/tests/test_game_theory.rs

@@ -34,7 +34,9 @@ fn test_mean_field_integration() {
 fn test_evolutionary_integration() {
     let payoff = DMatrix::from_row_slice(2, 2, &[3.0, 0.0, 5.0, 1.0]); // Prisoner's Dilemma-ish
     let system = ReplicatorDynamics::new(payoff).unwrap();
-    let traj = system.simulate(DVector::from_vec(vec![0.5, 0.5]), 1.0, 0.1);
+    let traj = system
+        .simulate(DVector::from_vec(vec![0.5, 0.5]), 1.0, 0.1)
+        .unwrap();
     assert!(traj.len() > 0);
 }