diff --git a/.jules/refactor.md b/.jules/refactor.md index 0d5c13e..1273d8a 100644 --- a/.jules/refactor.md +++ b/.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` 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. diff --git a/math_explorer/src/epidemiology/error.rs b/math_explorer/src/epidemiology/error.rs index 7ccdede..05f43ee 100644 --- a/math_explorer/src/epidemiology/error.rs +++ b/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), } } } diff --git a/math_explorer/src/epidemiology/stochastic.rs b/math_explorer/src/epidemiology/stochastic.rs index e9a0bd0..26b8552 100644 --- a/math_explorer/src/epidemiology/stochastic.rs +++ b/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 { fn propensities(&self, state: &State) -> Vec; /// 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 { /// /// // 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,8 +56,12 @@ impl GillespieSolver { /// 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(&mut self, system: &S, state: &mut State) -> f64 + /// Returns `Ok(f64::INFINITY)` if no reactions can occur (total propensity is 0). + pub fn step( + &mut self, + system: &S, + state: &mut State, + ) -> Result where S: StochasticSystem, { @@ -64,7 +69,7 @@ impl GillespieSolver { 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 GillespieSolver { } // 3. Update state - system.react(state, reaction_index); + system.react(state, reaction_index)?; - tau + Ok(tau) } } @@ -111,7 +116,7 @@ impl StochasticSystem 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 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()); diff --git a/math_explorer/src/physics/stat_mech/ising.rs b/math_explorer/src/physics/stat_mech/ising.rs index d58f7d0..e032fd4 100644 --- a/math_explorer/src/physics/stat_mech/ising.rs +++ b/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(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,28 +127,72 @@ 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 { + rng: R, +} + +impl MetropolisSolver { + /// Creates a new solver with the given RNG. + pub fn new(rng: R) -> Self { + Self { rng } + } +} + +impl IsingSolver for MetropolisSolver { + 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 { @@ -152,26 +200,25 @@ impl SpinLattice { } else { let beta = 1.0 / (KB * temperature); let prob = (-beta * delta_e).exp(); - rng.r#gen::() < prob + self.rng.r#gen::() < 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::() < prob { - self.spins[idx] = -s; + if prob > 1.0 || self.rng.r#gen::() < 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" + ); + } }