Refactor: Epidemiology safety and Physics determinism

.jules/refactor.md

@@ -13,3 +13,9 @@ Updated `li_ma_significance` to return `Result`.
 Updated `solve_linear_system` and `solve_normal_equation` to return `Result`.
 ## 2026-01-20 - [Physics/Chaos, Applied/Isosurface] **Debt:** Stringly Typed Errors **Risk:** Panic/Unpredictable Behavior
 Replaced `Result<T, String>` with typed error enums `ChaosError` and `IsosurfaceError`. This improves error handling safety and allows consumers to match on specific error cases instead of parsing strings.
+
+## 2026-06-21 - [Epidemiology/Stochastic] **Debt:** Panic in Library Code **Risk:** Panic risk
+`SIRModel::react` panics on invalid reaction index. This causes the entire simulation to crash instead of returning an error.
+
+## 2026-06-21 - [Physics/StatMech] **Debt:** Hardcoded RNG **Risk:** Architectural limitation
+`SpinLattice` and its methods hardcode `rand::thread_rng()`, making deterministic testing and reseeding impossible without mocking.

math_explorer/src/epidemiology/error.rs

@@ -14,6 +14,8 @@ pub enum EpidemiologyError {
         v_rows: usize,
         v_cols: usize,
     },
+    /// Invalid reaction index for the system.
+    InvalidReactionIndex(usize),
 }
 
 impl fmt::Display for EpidemiologyError {
@@ -36,6 +38,7 @@ impl fmt::Display for EpidemiologyError {
                 "Matrix dimensions mismatch: F=({}, {}), V=({}, {})",
                 f_rows, f_cols, v_rows, v_cols
             ),
+            Self::InvalidReactionIndex(idx) => write!(f, "Invalid reaction index: {}", idx),
         }
     }
 }

math_explorer/src/epidemiology/stochastic.rs

@@ -1,4 +1,5 @@
 use crate::epidemiology::compartmental::{SIRModel, SIRState};
+use crate::epidemiology::error::EpidemiologyError;
 use rand::Rng;
 
 /// A trait for systems that can be simulated stochastically.
@@ -10,7 +11,7 @@ pub trait StochasticSystem<State> {
     fn propensities(&self, state: &State) -> Vec<f64>;
 
     /// Updates the state according to the reaction that occurred.
-    fn react(&self, state: &mut State, reaction_index: usize);
+    fn react(&self, state: &mut State, reaction_index: usize) -> Result<(), EpidemiologyError>;
 }
 
 /// A solver for stochastic simulation using the Gillespie Algorithm (SSA).
@@ -34,7 +35,7 @@ pub trait StochasticSystem<State> {
 ///
 /// // Run for 10 time units
 /// while time < 10.0 {
-///     let dt = solver.step(&model, &mut state);
+///     let dt = solver.step(&model, &mut state).unwrap();
 ///     if dt.is_infinite() { break; }
 ///     time += dt;
 /// }
@@ -55,16 +56,20 @@ impl<R: Rng> GillespieSolver<R> {
     /// Performs one step of the Gillespie algorithm.
     ///
     /// Returns the time elapsed for this step.
-    /// Returns `f64::INFINITY` if no reactions can occur (total propensity is 0).
-    pub fn step<S, State>(&mut self, system: &S, state: &mut State) -> f64
+    /// Returns `Ok(f64::INFINITY)` if no reactions can occur (total propensity is 0).
+    pub fn step<S, State>(
+        &mut self,
+        system: &S,
+        state: &mut State,
+    ) -> Result<f64, EpidemiologyError>
     where
         S: StochasticSystem<State>,
     {
         let rates = system.propensities(state);
         let total_rate: f64 = rates.iter().sum();
 
         if total_rate <= 0.0 {
-            return f64::INFINITY;
+            return Ok(f64::INFINITY);
         }
 
         // 1. Determine time step tau
@@ -92,9 +97,9 @@ impl<R: Rng> GillespieSolver<R> {
         }
 
         // 3. Update state
-        system.react(state, reaction_index);
+        system.react(state, reaction_index)?;
 
-        tau
+        Ok(tau)
     }
 }
 
@@ -111,7 +116,7 @@ impl StochasticSystem<SIRState> for SIRModel {
         vec![infection_rate, recovery_rate]
     }
 
-    fn react(&self, state: &mut SIRState, reaction_index: usize) {
+    fn react(&self, state: &mut SIRState, reaction_index: usize) -> Result<(), EpidemiologyError> {
         match reaction_index {
             0 => {
                 // Infection: S decreases by 1, I increases by 1
@@ -123,8 +128,9 @@ impl StochasticSystem<SIRState> for SIRModel {
                 state.i -= 1.0;
                 state.r += 1.0;
             }
-            _ => panic!("Invalid reaction index for SIR Model"),
+            _ => return Err(EpidemiologyError::InvalidReactionIndex(reaction_index)),
         }
+        Ok(())
     }
 }
 
@@ -195,7 +201,7 @@ mod tests {
         let initial_i = state.i;
 
         // Take one step
-        let dt = solver.step(&model, &mut state);
+        let dt = solver.step(&model, &mut state).unwrap();
 
         assert!(dt > 0.0, "Time step should be positive");
         assert!(dt.is_finite());

math_explorer/src/physics/stat_mech/ising.rs

@@ -21,7 +21,7 @@ use rand::Rng;
 /// Simulating a ferromagnetic phase transition (Low Temperature):
 ///
 /// ```
-/// use math_explorer::physics::stat_mech::ising::SpinLattice;
+/// use math_explorer::physics::stat_mech::ising::{SpinLattice, MetropolisSolver, IsingSolver};
 /// use math_explorer::physics::stat_mech::KB;
 ///
 /// // 1. Setup system parameters
@@ -39,9 +39,8 @@ use rand::Rng;
 ///
 /// // 4. Run Metropolis Simulation to reach equilibrium
 /// // Note: This is a probabilistic process.
-/// for _ in 0..100_000 {
-///     lattice.metropolis_step(temp, j_coupling, h_field);
-/// }
+/// let mut solver = MetropolisSolver::new(rand::thread_rng());
+/// solver.evolve(&mut lattice, 200_000, temp, j_coupling, h_field);
 ///
 /// // 5. Check Magnetization
 /// let m = lattice.magnetization();
@@ -62,6 +61,11 @@ impl SpinLattice {
     /// Creates a new random spin lattice.
     pub fn new(width: usize, height: usize) -> Self {
         let mut rng = rand::thread_rng();
+        Self::new_with_rng(width, height, &mut rng)
+    }
+
+    /// Creates a new random spin lattice using a provided RNG.
+    pub fn new_with_rng<R: Rng + ?Sized>(width: usize, height: usize, rng: &mut R) -> Self {
         let count = width * height;
         let spins = (0..count)
             .map(|_| if rng.gen_bool(0.5) { 1 } else { -1 })
@@ -123,55 +127,98 @@ impl SpinLattice {
     /// 1. Select a random site.
     /// 2. Calculate energy cost to flip Delta E.
     /// 3. Accept flip if Delta E < 0 OR with probability e^(-beta * Delta E).
+    #[deprecated(since = "0.2.0", note = "Use MetropolisSolver::step instead")]
     pub fn metropolis_step(&mut self, temperature: f64, j_coupling: f64, h_field: f64) {
-        let mut rng = rand::thread_rng();
-        let x = rng.gen_range(0..self.width);
-        let y = rng.gen_range(0..self.height);
-
-        let s = self.get(x, y);
-
-        // Neighbors for energy difference (all 4 neighbors needed for Delta E)
-        let left_x = (x + self.width - 1) % self.width;
-        let right_x = (x + 1) % self.width;
-        let up_y = (y + self.height - 1) % self.height;
-        let down_y = (y + 1) % self.height;
-
-        let sum_neighbors = self.get(left_x, y) as f64
-            + self.get(right_x, y) as f64
-            + self.get(x, up_y) as f64
-            + self.get(x, down_y) as f64;
-
-        // Delta E = E_new - E_old
-        // E_old_local = -J * s * sum_neighbors - h * s
-        // E_new_local = -J * (-s) * sum_neighbors - h * (-s)
-        // Delta E = 2 * J * s * sum_neighbors + 2 * h * s
+        let rng = rand::thread_rng();
+        let mut solver = MetropolisSolver::new(rng);
+        solver.step(self, temperature, j_coupling, h_field);
+    }
+
+    /// Performs multiple Metropolis steps efficiently using a lookup table and batching.
+    ///
+    /// This method is significantly faster than calling `metropolis_step` in a loop.
+    #[deprecated(since = "0.2.0", note = "Use MetropolisSolver::evolve instead")]
+    pub fn evolve(&mut self, steps: usize, temperature: f64, j_coupling: f64, h_field: f64) {
+        let rng = rand::thread_rng();
+        let mut solver = MetropolisSolver::new(rng);
+        solver.evolve(self, steps, temperature, j_coupling, h_field);
+    }
+}
+
+/// A trait for strategies that evolve the Ising Model.
+pub trait IsingSolver {
+    /// Performs one step of the simulation.
+    fn step(&mut self, lattice: &mut SpinLattice, temperature: f64, j_coupling: f64, h_field: f64);
+
+    /// Evolves the simulation for a number of steps.
+    fn evolve(
+        &mut self,
+        lattice: &mut SpinLattice,
+        steps: usize,
+        temperature: f64,
+        j_coupling: f64,
+        h_field: f64,
+    );
+}
+
+/// A Metropolis-Hastings solver for the Ising Model.
+///
+/// Wraps an RNG to allow for deterministic simulations.
+pub struct MetropolisSolver<R> {
+    rng: R,
+}
+
+impl<R: Rng> MetropolisSolver<R> {
+    /// Creates a new solver with the given RNG.
+    pub fn new(rng: R) -> Self {
+        Self { rng }
+    }
+}
+
+impl<R: Rng> IsingSolver for MetropolisSolver<R> {
+    fn step(&mut self, lattice: &mut SpinLattice, temperature: f64, j_coupling: f64, h_field: f64) {
+        let x = self.rng.gen_range(0..lattice.width);
+        let y = self.rng.gen_range(0..lattice.height);
+
+        let s = lattice.get(x, y);
+
+        // Neighbors
+        let left_x = (x + lattice.width - 1) % lattice.width;
+        let right_x = (x + 1) % lattice.width;
+        let up_y = (y + lattice.height - 1) % lattice.height;
+        let down_y = (y + 1) % lattice.height;
+
+        let sum_neighbors = lattice.get(left_x, y) as f64
+            + lattice.get(right_x, y) as f64
+            + lattice.get(x, up_y) as f64
+            + lattice.get(x, down_y) as f64;
+
         let delta_e = 2.0 * s as f64 * (j_coupling * sum_neighbors + h_field);
 
         let should_flip = if delta_e < 0.0 {
             true
         } else {
             let beta = 1.0 / (KB * temperature);
             let prob = (-beta * delta_e).exp();
-            rng.r#gen::<f64>() < prob
+            self.rng.r#gen::<f64>() < prob
         };
 
         if should_flip {
-            self.set(x, y, -s);
+            lattice.set(x, y, -s);
         }
     }
 
-    /// Performs multiple Metropolis steps efficiently using a lookup table and batching.
-    ///
-    /// This method is significantly faster than calling `metropolis_step` in a loop because:
-    /// 1. It precomputes Boltzmann factors (`exp(-beta * dE)`) into a lookup table.
-    /// 2. It reuses the random number generator, avoiding TLS overhead.
-    pub fn evolve(&mut self, steps: usize, temperature: f64, j_coupling: f64, h_field: f64) {
-        let mut rng = rand::thread_rng();
+    fn evolve(
+        &mut self,
+        lattice: &mut SpinLattice,
+        steps: usize,
+        temperature: f64,
+        j_coupling: f64,
+        h_field: f64,
+    ) {
         let beta = 1.0 / (KB * temperature);
 
-        // Precompute probabilities: lut[s_idx][sum_idx]
-        // s_idx: 0 for s=-1, 1 for s=1
-        // sum_idx: 0..5 for sum in {-4, -2, 0, 2, 4} mapped via (sum + 4) / 2
+        // Precompute probabilities
         let mut lut = [[0.0; 5]; 2];
 
         for (s_idx, lut_row) in lut.iter_mut().enumerate() {
@@ -181,49 +228,40 @@ impl SpinLattice {
                 let delta_e = 2.0 * s_val * (j_coupling * sum_val + h_field);
 
                 if delta_e <= 0.0 {
-                    *entry = 1.1; // Always accept (value > 1.0 ensures check passes)
+                    *entry = 1.1; // Always accept
                 } else {
                     *entry = (-beta * delta_e).exp();
                 }
             }
         }
 
-        let width = self.width;
-        let height = self.height;
+        let width = lattice.width;
+        let height = lattice.height;
 
         for _ in 0..steps {
-            let x = rng.gen_range(0..width);
-            let y = rng.gen_range(0..height);
+            let x = self.rng.gen_range(0..width);
+            let y = self.rng.gen_range(0..height);
 
-            // Manual inline of get() to help optimizer
             let idx = y * width + x;
-            let s = self.spins[idx];
+            let s = lattice.spins[idx];
 
-            // Neighbors
-            // Use wrapping arithmetic or simple checks to avoid modulo if possible,
-            // but for random access, modulo is robust.
-            // Optimizing modulo:
             let left_x = if x == 0 { width - 1 } else { x - 1 };
             let right_x = if x == width - 1 { 0 } else { x + 1 };
             let up_y = if y == 0 { height - 1 } else { y - 1 };
             let down_y = if y == height - 1 { 0 } else { y + 1 };
 
-            let neighbor_sum = self.spins[y * width + left_x] as i32
-                + self.spins[y * width + right_x] as i32
-                + self.spins[up_y * width + x] as i32
-                + self.spins[down_y * width + x] as i32;
+            let neighbor_sum = lattice.spins[y * width + left_x] as i32
+                + lattice.spins[y * width + right_x] as i32
+                + lattice.spins[up_y * width + x] as i32
+                + lattice.spins[down_y * width + x] as i32;
 
-            // Map s (-1 or 1) to 0 or 1
             let s_idx = if s == -1 { 0 } else { 1 };
-            // Map neighbor_sum (-4..4) to 0..4
             let sum_idx = ((neighbor_sum + 4) / 2) as usize;
 
             let prob = lut[s_idx][sum_idx];
 
-            // If prob > 1.0, it was delta_e < 0, so flip.
-            // Else compare with random.
-            if prob > 1.0 || rng.r#gen::<f64>() < prob {
-                self.spins[idx] = -s;
+            if prob > 1.0 || self.rng.r#gen::<f64>() < prob {
+                lattice.spins[idx] = -s;
             }
         }
     }
@@ -234,6 +272,9 @@ mod tests {
     use super::*;
     use crate::physics::stat_mech::KB;
 
+    use rand::SeedableRng;
+    use rand::rngs::StdRng;
+
     #[test]
     fn test_ising_disordered() {
         // High temperature limit. T >> J/k_B.
@@ -251,9 +292,8 @@ mod tests {
 
         // Run many steps to equilibrate
         let steps = 10_000;
-        for _ in 0..steps {
-            lattice.metropolis_step(t_high, j_val, h_val);
-        }
+        let mut solver = MetropolisSolver::new(rand::thread_rng());
+        solver.evolve(&mut lattice, steps, t_high, j_val, h_val);
 
         let m = lattice.magnetization();
         let max_m = (width * height) as f64;
@@ -268,4 +308,33 @@ mod tests {
             avg_m_per_spin
         );
     }
+
+    #[test]
+    fn test_ising_deterministic() {
+        let mut rng = StdRng::seed_from_u64(42);
+        let width = 10;
+        let height = 10;
+        let mut lattice = SpinLattice::new_with_rng(width, height, &mut rng);
+
+        // Capture initial state
+        let initial_spins = lattice.spins.clone();
+
+        let mut solver = MetropolisSolver::new(rng);
+        // Run a few steps
+        solver.evolve(&mut lattice, 100, 1.0, 1.0, 0.0);
+
+        // Ensure state changed
+        assert_ne!(lattice.spins, initial_spins);
+
+        // Re-run with same seed
+        let mut rng2 = StdRng::seed_from_u64(42);
+        let mut lattice2 = SpinLattice::new_with_rng(width, height, &mut rng2);
+        let mut solver2 = MetropolisSolver::new(rng2);
+        solver2.evolve(&mut lattice2, 100, 1.0, 1.0, 0.0);
+
+        assert_eq!(
+            lattice.spins, lattice2.spins,
+            "Simulation should be deterministic"
+        );
+    }
 }