Skip to content

⚙️ Core: Strategy Pattern for Ising Model #368

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
fderuiter wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

⚙️ Core: Strategy Pattern for Ising Model #368

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions .jules/systems_core.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,3 +58,8 @@
**Violation:** The `gillespie_step_time` function in `epidemiology::stochastic` had hardcoded `rand::thread_rng()` dependency (Side Effect) and hardcoded 2-reaction logic, violating Dependency Inversion and Open/Closed Principle.
**Refactor:** Applied the **Strategy Pattern**. Defined `StochasticSystem` trait for reaction networks and `GillespieSolver` struct with injected RNG (`R: Rng`).
**Trade-off:** Increased complexity (Generics, Trait definition) compared to a single function, but enabled deterministic testing (seeded RNG) and support for any reaction network (Composability).

## 2027-04-10 - Strategy Pattern for Ising Model
**Violation:** `SpinLattice` in `physics::stat_mech` was coupled to the Metropolis algorithm via hardcoded `metropolis_step` and `evolve` methods, violating OCP and hindering the addition of advanced solvers (e.g., Wolff).
**Refactor:** Applied the **Strategy Pattern**. Extracted `IsingSolver` trait and implemented `MetropolisSolver`. Refactored `SpinLattice` to delegate to these solvers, maintaining backward compatibility via deprecated methods.
**Trade-off:** Minimal boilerplate for trait/struct definition. Enables swapping update rules and deterministic testing (via injected RNG), while preserving high-performance batch updates.
1 change: 1 addition & 0 deletions math_explorer/Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,3 +25,4 @@ path = "tests/algorithmic_information.rs"

[dev-dependencies]
approx = "0.5.1"
rand = "0.8.5"
11 changes: 7 additions & 4 deletions math_explorer/examples/bench_ising_custom.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
use math_explorer::physics::stat_mech::KB;
use math_explorer::physics::stat_mech::ising::SpinLattice;
use math_explorer::physics::stat_mech::ising::{IsingSolver, MetropolisSolver, SpinLattice};
use std::time::Instant;

fn main() {
Expand Down Expand Up @@ -31,14 +31,17 @@ fn main() {
old_speed
);

// --- After: evolve batch ---
// --- After: MetropolisSolver Strategy ---
let mut lattice = SpinLattice::new(width, height);
// Explicitly using the new Strategy Pattern
let mut solver = MetropolisSolver::new(rand::thread_rng());

let start = Instant::now();
lattice.evolve(iterations, temp, j_coupling, h_field);
solver.evolve(&mut lattice, iterations, temp, j_coupling, h_field);
let duration = start.elapsed();
let new_speed = (iterations as f64 / duration.as_secs_f64()) / 1_000_000.0;
println!(
"After (evolve): {:.4}s, {:.2} M/s",
"After (MetropolisSolver): {:.4}s, {:.2} M/s",
duration.as_secs_f64(),
new_speed
);
Expand Down
265 changes: 162 additions & 103 deletions math_explorer/src/physics/stat_mech/ising.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,152 @@
//! - **$T > T_c$**: Paramagnetic phase (Disordered spins).

use super::KB;
use rand::Rng;
use rand::{thread_rng, Rng};

/// Defines the strategy for evolving the Ising Model.
pub trait IsingSolver {
/// Performs a single update step on the lattice.
fn step(&mut self, lattice: &mut SpinLattice, temperature: f64, j_coupling: f64, h_field: f64);

/// Evolves the lattice for a given number of steps.
///
/// Default implementation loops over `step`, but solvers may override for optimization.
fn evolve(
&mut self,
lattice: &mut SpinLattice,
steps: usize,
temperature: f64,
j_coupling: f64,
h_field: f64,
) {
for _ in 0..steps {
self.step(lattice, temperature, j_coupling, h_field);
}
}
}

/// A Metropolis-Hastings solver for the Ising Model.
///
/// This solver implements the standard Metropolis algorithm.
/// It holds an internal RNG to ensure reproducibility if seeded,
/// and reuses it for batch operations to maximize performance.
pub struct MetropolisSolver<R: Rng> {
rng: R,
}

impl<R: Rng> MetropolisSolver<R> {
/// Creates a new MetropolisSolver 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 for energy difference (all 4 neighbors needed for Delta E)
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;

// Delta E = 2 * J * s * sum_neighbors + 2 * h * s
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();
self.rng.r#gen::<f64>() < prob
};

if should_flip {
lattice.set(x, y, -s);
}
}

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
let mut lut = [[0.0; 5]; 2];

for (s_idx, lut_row) in lut.iter_mut().enumerate() {
let s_val = if s_idx == 0 { -1.0 } else { 1.0 };
for (sum_idx, entry) in lut_row.iter_mut().enumerate() {
let sum_val = (sum_idx as f64 * 2.0) - 4.0;
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)
} else {
*entry = (-beta * delta_e).exp();
}
}
}

let width = lattice.width;
let height = lattice.height;

for _ in 0..steps {
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 = lattice.spins[idx];

// Neighbors
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 = 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 || self.rng.r#gen::<f64>() < prob {
lattice.spins[idx] = -s;
}
}
}
}

/// A 2D Spin Lattice for the Ising Model.
///
Expand Down Expand Up @@ -59,9 +204,13 @@ pub struct SpinLattice {
}

impl SpinLattice {
/// Creates a new random spin lattice.
/// Creates a new random spin lattice using the thread-local RNG.
pub fn new(width: usize, height: usize) -> Self {
let mut rng = rand::thread_rng();
Self::new_with_rng(width, height, &mut rand::thread_rng())
}

/// Creates a new random spin lattice using a provided RNG.
pub fn new_with_rng<R: 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 })
Expand Down Expand Up @@ -120,112 +269,22 @@ impl SpinLattice {

/// Performs one Metropolis algorithm step (one attempt).
///
/// 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
/// Use `MetropolisSolver` or the `IsingSolver` trait instead.
#[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 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
};

if should_flip {
self.set(x, y, -s);
}
let mut solver = MetropolisSolver::new(thread_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 because:
/// 1. It precomputes Boltzmann factors (`exp(-beta * dE)`) into a lookup table.
/// 2. It reuses the random number generator, avoiding TLS overhead.
/// # Deprecated
/// Use `MetropolisSolver` or the `IsingSolver` trait instead.
#[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 mut rng = rand::thread_rng();
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
let mut lut = [[0.0; 5]; 2];

for (s_idx, lut_row) in lut.iter_mut().enumerate() {
let s_val = if s_idx == 0 { -1.0 } else { 1.0 };
for (sum_idx, entry) in lut_row.iter_mut().enumerate() {
let sum_val = (sum_idx as f64 * 2.0) - 4.0;
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)
} else {
*entry = (-beta * delta_e).exp();
}
}
}

let width = self.width;
let height = self.height;

for _ in 0..steps {
let x = rng.gen_range(0..width);
let y = rng.gen_range(0..height);

// Manual inline of get() to help optimizer
let idx = y * width + x;
let s = self.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;

// 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;
}
}
let mut solver = MetropolisSolver::new(thread_rng());
solver.evolve(self, steps, temperature, j_coupling, h_field);
}
}

Expand Down
Loading