Skip to content

🧱 Mason: Consolidating Reaction-Diffusion Logic (Fixes SRP/DRY) #362

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

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

🧱 Mason: Consolidating Reaction-Diffusion Logic (Fixes SRP/DRY) #362

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
4 changes: 4 additions & 0 deletions .jules/mason.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,3 +26,7 @@
## 2026-06-25 - Decoupling LoraHub Algorithms
**Violation:** Open/Closed Principle (OCP). The `LoraEnsemble` hardcoded `combine` (weighted sum) and `evaluate_objective` (L1 Regularization), preventing extension for advanced merging techniques like SLERP or TIES.
**Remedy:** Strategy Pattern. Extracted `CombinationStrategy` and `ObjectiveStrategy` traits. Implemented `LinearCombinationStrategy` and `L1RegularizationStrategy`. Refactored `LoraEnsemble` to use these strategies via dependency injection.

## 2026-07-22 - Consolidating Reaction-Diffusion Logic
**Violation:** Single Responsibility Principle (SRP) and Don't Repeat Yourself (DRY). The `TuringSystem` duplicated complex stencil iteration logic in both its inherent `step` method (Euler integration) and `derivative_in_place` implementation (OdeSystem). This coupled the spatial discretization tightly with the time integration method.
**Remedy:** Extracted Method with Closure. Extracted `apply_reaction_diffusion_stencil` which handles the spatial discretization and delegates the update rule to a closure. This centralizes the logic and decouples the "what" (physics) from the "how" (state update).
181 changes: 76 additions & 105 deletions math_explorer/src/biology/morphogenesis.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -235,26 +235,28 @@ impl<K: ReactionKinetics> TuringSystem<K> {
self.state.v_mut()
}

/// Updates the grid using a finite-difference Laplacian and reaction kinetics.
pub fn step(&mut self, dt: f64) {
let n = self.state.len();
/// Applies the reaction-diffusion stencil to the grid, delegating the update rule to a closure.
fn apply_reaction_diffusion_stencil<F>(
state: &TuringState,
kinetics: &K,
d_u: f64,
d_v: f64,
dx: f64,
mut writer: F,
) where
F: FnMut(usize, f64, f64, f64, f64), // i, u_curr, v_curr, du_dt, dv_dt
{
let n = state.len();
if n == 0 {
return;
}

// Ensure buffers are the right size
if self.next_state.len() != n {
self.next_state = TuringState::new(n);
}

let dx_sq = self.dx * self.dx;
let dx_sq = dx * dx;
let inv_dx_sq = 1.0 / dx_sq;

// Optimization: Lift boundary checks out of the loop and use slices
let u = &self.state.u;
let v = &self.state.v;
let next_u = &mut self.next_state.u;
let next_v = &mut self.next_state.v;
let u = &state.u;
let v = &state.v;

// 1. Handle i = 0
{
Expand All @@ -274,12 +276,12 @@ impl<K: ReactionKinetics> TuringSystem<K> {
let lap_u = (u_next - 2.0 * u_curr + u_prev) * inv_dx_sq;
let lap_v = (v_next - 2.0 * v_curr + v_prev) * inv_dx_sq;

let (reaction_u, reaction_v) = self.kinetics.reaction(u_curr, v_curr);
let (reaction_u, reaction_v) = kinetics.reaction(u_curr, v_curr);

unsafe {
*next_u.get_unchecked_mut(i) = u_curr + dt * (self.d_u * lap_u + reaction_u);
*next_v.get_unchecked_mut(i) = v_curr + dt * (self.d_v * lap_v + reaction_v);
}
let du_dt = d_u * lap_u + reaction_u;
let dv_dt = d_v * lap_v + reaction_v;

writer(i, u_curr, v_curr, du_dt, dv_dt);
}

// 2. Handle i = 1..n-1 (Hot Path)
Expand All @@ -298,10 +300,12 @@ impl<K: ReactionKinetics> TuringSystem<K> {
let lap_u = (u_next - 2.0 * u_curr + u_prev) * inv_dx_sq;
let lap_v = (v_next - 2.0 * v_curr + v_prev) * inv_dx_sq;

let (reaction_u, reaction_v) = self.kinetics.reaction(u_curr, v_curr);
let (reaction_u, reaction_v) = kinetics.reaction(u_curr, v_curr);

*next_u.get_unchecked_mut(i) = u_curr + dt * (self.d_u * lap_u + reaction_u);
*next_v.get_unchecked_mut(i) = v_curr + dt * (self.d_v * lap_v + reaction_v);
let du_dt = d_u * lap_u + reaction_u;
let dv_dt = d_v * lap_v + reaction_v;

writer(i, u_curr, v_curr, du_dt, dv_dt);

// Shift window
u_prev = u_curr;
Expand All @@ -326,12 +330,46 @@ impl<K: ReactionKinetics> TuringSystem<K> {
let lap_u = (u_next - 2.0 * u_curr + u_prev) * inv_dx_sq;
let lap_v = (v_next - 2.0 * v_curr + v_prev) * inv_dx_sq;

let (reaction_u, reaction_v) = self.kinetics.reaction(u_curr, v_curr);
let (reaction_u, reaction_v) = kinetics.reaction(u_curr, v_curr);

let du_dt = d_u * lap_u + reaction_u;
let dv_dt = d_v * lap_v + reaction_v;

*next_u.get_unchecked_mut(i) = u_curr + dt * (self.d_u * lap_u + reaction_u);
*next_v.get_unchecked_mut(i) = v_curr + dt * (self.d_v * lap_v + reaction_v);
writer(i, u_curr, v_curr, du_dt, dv_dt);
}
}
}

/// Updates the grid using a finite-difference Laplacian and reaction kinetics.
pub fn step(&mut self, dt: f64) {
let n = self.state.len();
if n == 0 {
return;
}

// Ensure buffers are the right size
if self.next_state.len() != n {
self.next_state = TuringState::new(n);
}

let next_u = &mut self.next_state.u;
let next_v = &mut self.next_state.v;

// Split borrow: pass immutable references to parameters and state
Self::apply_reaction_diffusion_stencil(
&self.state,
&self.kinetics,
self.d_u,
self.d_v,
self.dx,
|i, u_curr, v_curr, du_dt, dv_dt| {
// Forward Euler Step
unsafe {
*next_u.get_unchecked_mut(i) = u_curr + dt * du_dt;
*next_v.get_unchecked_mut(i) = v_curr + dt * dv_dt;
}
},
);

// Swap buffers (states)
std::mem::swap(&mut self.state, &mut self.next_state);
Expand All @@ -356,91 +394,24 @@ impl<K: ReactionKinetics> OdeSystem<TuringState> for TuringSystem<K> {
*out = TuringState::new(n);
}

let dx_sq = self.dx * self.dx;
let inv_dx_sq = 1.0 / dx_sq;

// Optimization: Lift boundary checks out of the loop and use slices
let u = &state.u;
let v = &state.v;
let out_u = &mut out.u;
let out_v = &mut out.v;

// 1. Handle i = 0
{
let i = 0;
// Safety: n > 0 checked above
let u_curr = unsafe { *u.get_unchecked(i) };
let v_curr = unsafe { *v.get_unchecked(i) };

let u_prev = u_curr;
let v_prev = v_curr;
let (u_next, v_next) = if n > 1 {
unsafe { (*u.get_unchecked(1), *v.get_unchecked(1)) }
} else {
(u_curr, v_curr)
};

let lap_u = (u_next - 2.0 * u_curr + u_prev) * inv_dx_sq;
let lap_v = (v_next - 2.0 * v_curr + v_prev) * inv_dx_sq;

let (reaction_u, reaction_v) = self.kinetics.reaction(u_curr, v_curr);

unsafe {
*out_u.get_unchecked_mut(i) = self.d_u * lap_u + reaction_u;
*out_v.get_unchecked_mut(i) = self.d_v * lap_v + reaction_v;
}
}

// 2. Handle i = 1..n-1 (Hot Path)
if n > 2 {
// Optimization: Sliding Window / Register Rotation
unsafe {
let mut u_prev = *u.get_unchecked(0);
let mut u_curr = *u.get_unchecked(1);
let mut v_prev = *v.get_unchecked(0);
let mut v_curr = *v.get_unchecked(1);

for i in 1..n - 1 {
let u_next = *u.get_unchecked(i + 1);
let v_next = *v.get_unchecked(i + 1);

let lap_u = (u_next - 2.0 * u_curr + u_prev) * inv_dx_sq;
let lap_v = (v_next - 2.0 * v_curr + v_prev) * inv_dx_sq;

let (reaction_u, reaction_v) = self.kinetics.reaction(u_curr, v_curr);

*out_u.get_unchecked_mut(i) = self.d_u * lap_u + reaction_u;
*out_v.get_unchecked_mut(i) = self.d_v * lap_v + reaction_v;

// Shift window
u_prev = u_curr;
u_curr = u_next;
v_prev = v_curr;
v_curr = v_next;
// Use the shared stencil logic
Self::apply_reaction_diffusion_stencil(
state,
&self.kinetics,
self.d_u,
self.d_v,
self.dx,
|i, _u_curr, _v_curr, du_dt, dv_dt| {
// Just store the derivative
unsafe {
*out_u.get_unchecked_mut(i) = du_dt;
*out_v.get_unchecked_mut(i) = dv_dt;
}
}
}

// 3. Handle i = n-1
if n > 1 {
let i = n - 1;
unsafe {
let u_curr = *u.get_unchecked(i);
let v_curr = *v.get_unchecked(i);
let u_prev = *u.get_unchecked(i - 1);
let v_prev = *v.get_unchecked(i - 1);
let u_next = u_curr;
let v_next = v_curr;

let lap_u = (u_next - 2.0 * u_curr + u_prev) * inv_dx_sq;
let lap_v = (v_next - 2.0 * v_curr + v_prev) * inv_dx_sq;

let (reaction_u, reaction_v) = self.kinetics.reaction(u_curr, v_curr);

*out_u.get_unchecked_mut(i) = self.d_u * lap_u + reaction_u;
*out_v.get_unchecked_mut(i) = self.d_v * lap_v + reaction_v;
}
}
},
);
}
}

Expand Down