From 3a8585cb90b76c3841a563fa12c9de778cae3c10 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Fri, 30 Jan 2026 19:15:40 +0000 Subject: [PATCH 1/2] Refactor Dose Calculation to use Strategy Pattern and Domain Types - Defined `DoseError` enum for typed error handling. - Introduced `TermaInput` and `KernelInput` parameter objects with strict validation. - Extracted `TermaModel` and `DoseKernel` traits to decouple algorithms. - Implemented `ExponentialTerma` and `PointSpreadKernel` strategies. - Deprecated legacy functions `calculate_terma` and `point_kernel`, delegating to new strategies. - Added comprehensive unit tests for validation and strategy logic. This refactor addresses Primitive Obsession and improves Composability/OCP. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- .jules/systems_core.md | 5 + .../src/physics/medical/dose_calculation.rs | 160 ++++++++++++++++-- 2 files changed, 152 insertions(+), 13 deletions(-) diff --git a/.jules/systems_core.md b/.jules/systems_core.md index 183b3ca..f57d2c2 100644 --- a/.jules/systems_core.md +++ b/.jules/systems_core.md @@ -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-02-14 - Strategy Pattern for Dose Calculation +**Violation:** The `dose_calculation.rs` module suffered from primitive obsession (raw `f64` arguments) and lacked extensibility, violating the Open/Closed Principle. +**Refactor:** Applied **Strategy Pattern** and **Parameter Object**. Extracted `TermaModel` and `DoseKernel` traits. Introduced `TermaInput` and `KernelInput` for strict validation. +**Trade-off:** Increased verbosity for simple calculations, but enabled type-safe, extensible dose algorithms and centralized validation logic. diff --git a/math_explorer/src/physics/medical/dose_calculation.rs b/math_explorer/src/physics/medical/dose_calculation.rs index 7857533..7a33ff3 100644 --- a/math_explorer/src/physics/medical/dose_calculation.rs +++ b/math_explorer/src/physics/medical/dose_calculation.rs @@ -1,4 +1,113 @@ //! Dose Calculation Algorithms. +//! +//! # Refactoring Note +//! +//! This module has been refactored to use the Strategy Pattern for Dose Calculation. +//! Legacy functions are preserved but deprecated. + +use std::fmt; + +// --- Domain Errors --- + +#[derive(Debug, Clone, PartialEq)] +pub enum DoseError { + InvalidInput(String), + CalculationError(String), +} + +impl fmt::Display for DoseError { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + match self { + DoseError::InvalidInput(msg) => write!(f, "Invalid Input: {}", msg), + DoseError::CalculationError(msg) => write!(f, "Calculation Error: {}", msg), + } + } +} + +impl std::error::Error for DoseError {} + +// --- Parameter Objects --- + +pub struct TermaInput { + pub incident_fluence: f64, + pub mu: f64, + pub depth: f64, +} + +impl TermaInput { + pub fn new(incident_fluence: f64, mu: f64, depth: f64) -> Result { + if incident_fluence < 0.0 { + return Err(DoseError::InvalidInput("Fluence must be non-negative".into())); + } + if mu < 0.0 { + return Err(DoseError::InvalidInput("Attenuation coefficient (mu) must be non-negative".into())); + } + if depth < 0.0 { + return Err(DoseError::InvalidInput("Depth must be non-negative".into())); + } + Ok(Self { incident_fluence, mu, depth }) + } +} + +pub struct KernelInput { + pub radius: f64, +} + +impl KernelInput { + pub fn new(radius: f64) -> Result { + if radius < 0.0 { + return Err(DoseError::InvalidInput("Radius must be non-negative".into())); + } + // Singularity check + if radius.abs() < 1e-6 { + return Err(DoseError::InvalidInput("Radius cannot be zero (singularity at r=0)".into())); + } + Ok(Self { radius }) + } +} + +// --- Strategies --- + +pub trait TermaModel { + fn calculate(&self, input: &TermaInput) -> Result; +} + +pub trait DoseKernel { + fn evaluate(&self, input: &KernelInput) -> Result; +} + +// --- Implementations --- + +pub struct ExponentialTerma; + +impl TermaModel for ExponentialTerma { + fn calculate(&self, input: &TermaInput) -> Result { + // T = mu * Psi0 * e^(-mu * d) + Ok(input.mu * input.incident_fluence * (-input.mu * input.depth).exp()) + } +} + +pub struct PointSpreadKernel { + pub amplitude: f64, + pub beta: f64, +} + +impl PointSpreadKernel { + pub fn new(amplitude: f64, beta: f64) -> Self { + Self { amplitude, beta } + } +} + +impl DoseKernel for PointSpreadKernel { + fn evaluate(&self, input: &KernelInput) -> Result { + let r = input.radius; + // K = (A / r^2) * e^(-beta * r) + let val = (self.amplitude / (r * r)) * (-self.beta * r).exp(); + Ok(val) + } +} + +// --- Legacy Functions --- /// Calculates the Total Energy Released per Mass (TERMA) for a ray segment. /// @@ -18,12 +127,15 @@ /// # Formula /// /// $T = \mu \Psi_0 e^{-\mu d}$ +#[deprecated(since = "0.2.0", note = "Use TermaModel strategy (ExponentialTerma) instead")] pub fn calculate_terma(incident_fluence: f64, mu: f64, depth: f64) -> f64 { - if incident_fluence < 0.0 || mu < 0.0 || depth < 0.0 { - // Physical quantities should be non-negative, but we return 0.0 or handle gracefully. - return 0.0; - } - mu * incident_fluence * (-mu * depth).exp() + let input = match TermaInput::new(incident_fluence, mu, depth) { + Ok(i) => i, + Err(_) => return 0.0, // Legacy behavior + }; + + let model = ExponentialTerma; + model.calculate(&input).unwrap_or(0.0) } /// Calculates a simplified analytical Point Spread Function (Kernel). @@ -47,16 +159,12 @@ pub fn calculate_terma(incident_fluence: f64, mu: f64, depth: f64) -> f64 { /// /// *Note*: This is a singular kernel at r=0. In practice, finite voxel size integration is used. /// Here we return an error or handle the singularity if r is too close to 0. +#[deprecated(since = "0.2.0", note = "Use DoseKernel strategy (PointSpreadKernel) instead")] pub fn point_kernel(radius: f64, amplitude: f64, beta: f64) -> Result { - if radius.abs() < 1e-6 { - return Err("Radius cannot be zero (singularity at r=0)".to_string()); - } - if radius < 0.0 { - return Err("Radius must be non-negative".to_string()); - } + let input = KernelInput::new(radius).map_err(|e| e.to_string())?; - let val = (amplitude / (radius * radius)) * (-beta * radius).exp(); - Ok(val) + let kernel = PointSpreadKernel::new(amplitude, beta); + kernel.evaluate(&input).map_err(|e| e.to_string()) } #[cfg(test)] @@ -64,6 +172,7 @@ mod tests { use super::*; #[test] + #[allow(deprecated)] fn test_terma_calculation() { // Simple case: no attenuation (mu=0) -> T = 0 assert_eq!(calculate_terma(100.0, 0.0, 10.0), 0.0); @@ -74,6 +183,7 @@ mod tests { } #[test] + #[allow(deprecated)] fn test_point_kernel() { // Error on r=0 assert!(point_kernel(0.0, 1.0, 1.0).is_err()); @@ -86,6 +196,30 @@ mod tests { let k = point_kernel(r, a, b).unwrap(); assert!((k - (-1.0_f64).exp()).abs() < 1e-5); } + + // New tests for refactored code + #[test] + fn test_terma_input_validation() { + assert!(TermaInput::new(-1.0, 0.1, 1.0).is_err()); + assert!(TermaInput::new(100.0, -0.1, 1.0).is_err()); + assert!(TermaInput::new(100.0, 0.1, -1.0).is_err()); + } + + #[test] + fn test_exponential_terma_strategy() { + let input = TermaInput::new(100.0, 0.1, 5.0).unwrap(); + let model = ExponentialTerma; + let result = model.calculate(&input).unwrap(); + // 0.1 * 100 * exp(-0.1 * 5) = 10 * exp(-0.5) + let expected = 10.0 * (-0.5_f64).exp(); + assert!((result - expected).abs() < 1e-6); + } + + #[test] + fn test_kernel_input_validation() { + assert!(KernelInput::new(-1.0).is_err()); + assert!(KernelInput::new(0.0).is_err()); // Singularity + } } /// Calculates the average energy for the Beam Loading Line. From 9aa344cf2746bd106ecbb016c3592ed131524d0f Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Fri, 30 Jan 2026 19:21:05 +0000 Subject: [PATCH 2/2] Fix formatting in dose_calculation.rs Ran `cargo fmt` to address CI failure. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- .../src/physics/medical/dose_calculation.rs | 32 +++++++++++++++---- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/math_explorer/src/physics/medical/dose_calculation.rs b/math_explorer/src/physics/medical/dose_calculation.rs index 7a33ff3..60980d5 100644 --- a/math_explorer/src/physics/medical/dose_calculation.rs +++ b/math_explorer/src/physics/medical/dose_calculation.rs @@ -37,15 +37,23 @@ pub struct TermaInput { impl TermaInput { pub fn new(incident_fluence: f64, mu: f64, depth: f64) -> Result { if incident_fluence < 0.0 { - return Err(DoseError::InvalidInput("Fluence must be non-negative".into())); + return Err(DoseError::InvalidInput( + "Fluence must be non-negative".into(), + )); } if mu < 0.0 { - return Err(DoseError::InvalidInput("Attenuation coefficient (mu) must be non-negative".into())); + return Err(DoseError::InvalidInput( + "Attenuation coefficient (mu) must be non-negative".into(), + )); } if depth < 0.0 { return Err(DoseError::InvalidInput("Depth must be non-negative".into())); } - Ok(Self { incident_fluence, mu, depth }) + Ok(Self { + incident_fluence, + mu, + depth, + }) } } @@ -56,11 +64,15 @@ pub struct KernelInput { impl KernelInput { pub fn new(radius: f64) -> Result { if radius < 0.0 { - return Err(DoseError::InvalidInput("Radius must be non-negative".into())); + return Err(DoseError::InvalidInput( + "Radius must be non-negative".into(), + )); } // Singularity check if radius.abs() < 1e-6 { - return Err(DoseError::InvalidInput("Radius cannot be zero (singularity at r=0)".into())); + return Err(DoseError::InvalidInput( + "Radius cannot be zero (singularity at r=0)".into(), + )); } Ok(Self { radius }) } @@ -127,7 +139,10 @@ impl DoseKernel for PointSpreadKernel { /// # Formula /// /// $T = \mu \Psi_0 e^{-\mu d}$ -#[deprecated(since = "0.2.0", note = "Use TermaModel strategy (ExponentialTerma) instead")] +#[deprecated( + since = "0.2.0", + note = "Use TermaModel strategy (ExponentialTerma) instead" +)] pub fn calculate_terma(incident_fluence: f64, mu: f64, depth: f64) -> f64 { let input = match TermaInput::new(incident_fluence, mu, depth) { Ok(i) => i, @@ -159,7 +174,10 @@ pub fn calculate_terma(incident_fluence: f64, mu: f64, depth: f64) -> f64 { /// /// *Note*: This is a singular kernel at r=0. In practice, finite voxel size integration is used. /// Here we return an error or handle the singularity if r is too close to 0. -#[deprecated(since = "0.2.0", note = "Use DoseKernel strategy (PointSpreadKernel) instead")] +#[deprecated( + since = "0.2.0", + note = "Use DoseKernel strategy (PointSpreadKernel) instead" +)] pub fn point_kernel(radius: f64, amplitude: f64, beta: f64) -> Result { let input = KernelInput::new(radius).map_err(|e| e.to_string())?;