From fe7cfa476a219c9b6951303404f0c2ae923b8a51 Mon Sep 17 00:00:00 2001 From: Markus <66058642+mhovd@users.noreply.github.com> Date: Fri, 24 Apr 2026 14:14:16 +0200 Subject: [PATCH 1/2] refactor: Rename support points to parameters Simulator is agnostic to NP or P approach, API should reflect that --- examples/compare_solvers.rs | 4 +- examples/gendata.rs | 10 +-- src/data/structs.rs | 14 ++-- src/lib.rs | 2 +- src/optimize/mod.rs | 2 +- src/optimize/{spp.rs => parameters.rs} | 18 ++--- src/simulator/cache.rs | 4 +- src/simulator/equation/analytical/mod.rs | 42 ++++++------ src/simulator/equation/mod.rs | 50 +++++++------- src/simulator/equation/ode/mod.rs | 83 ++++++++++++------------ src/simulator/equation/sde/mod.rs | 50 +++++++------- src/simulator/likelihood/matrix.rs | 30 ++++----- src/simulator/likelihood/mod.rs | 4 +- src/simulator/likelihood/subject.rs | 2 +- tests/test_solvers.rs | 4 +- 15 files changed, 162 insertions(+), 157 deletions(-) rename src/optimize/{spp.rs => parameters.rs} (78%) diff --git a/examples/compare_solvers.rs b/examples/compare_solvers.rs index 4d534c29..69c57424 100644 --- a/examples/compare_solvers.rs +++ b/examples/compare_solvers.rs @@ -51,7 +51,7 @@ fn main() { .missing_observation(24.0, 0) .build(); - let spp = vec![0.1, 0.05, 0.03, 50.0]; // ke, kcp, kpc, V + let parameters = vec![0.1, 0.05, 0.03, 50.0]; // ke, kcp, kpc, V // Run each solver and collect predictions let bdf = two_cpt(OdeSolver::Bdf); @@ -69,7 +69,7 @@ fn main() { // ── Run all solvers and collect results ─────────────────────── let mut rows: Vec<(&str, u128, Vec)> = Vec::new(); for (name, ode) in &results { - let (preds, us) = timed(|| ode.estimate_predictions(&subject, &spp).unwrap()); + let (preds, us) = timed(|| ode.estimate_predictions(&subject, ¶meters).unwrap()); let preds: Vec = preds.flat_predictions().to_vec(); rows.push((name, us, preds)); } diff --git a/examples/gendata.rs b/examples/gendata.rs index 067639da..e7a6a346 100644 --- a/examples/gendata.rs +++ b/examples/gendata.rs @@ -45,8 +45,8 @@ fn main() { // let v_dist = rand_distr::Normal::new(50.0, 10.0).unwrap(); // let ske_dist = rand_distr::Normal::new(0.1, 0.01).unwrap(); - let mut support_points = vec![]; - let mut file = File::create("spp.csv").unwrap(); + let mut parameters = vec![]; + let mut file = File::create("parameters.csv").unwrap(); for _ in 0..100 { let ke = ke_dist.sample(&mut rand::rng()); // let ke = 1.2; @@ -55,14 +55,14 @@ fn main() { // let ske = ske_dist.sample(&mut rand::thread_rng()); // let v = v_dist.sample(&mut rand::thread_rng()); let v = 50.0; - support_points.push(vec![ke]); + parameters.push(vec![ke]); println!("{ke}, {ske}, {v}"); writeln!(file, "{}, {}, {}", ke, ske, v).unwrap(); } let mut data = vec![]; - for (i, spp) in support_points.iter().enumerate() { - let trajectories = sde.estimate_predictions(&subject, spp).unwrap(); + for (i, parameter) in parameters.iter().enumerate() { + let trajectories = sde.estimate_predictions(&subject, parameter).unwrap(); let trajectory = trajectories.row(0); // dbg!(&trajectory); let mut sb = Subject::builder(format!("id{}", i)).bolus(0.0, 20.0, 0); diff --git a/src/data/structs.rs b/src/data/structs.rs index 2acb676a..396cc04c 100644 --- a/src/data/structs.rs +++ b/src/data/structs.rs @@ -592,12 +592,12 @@ impl Occasion { } fn add_lagtime(&mut self, reorder: Option<(&Fa, &Lag, &[f64], &Covariates)>) { - if let Some((_, fn_lag, spp, covariates)) = reorder { - let spp = nalgebra::DVector::from_vec(spp.to_vec()); + if let Some((_, fn_lag, parameters, covariates)) = reorder { + let parameters = nalgebra::DVector::from_vec(parameters.to_vec()); for event in self.events.iter_mut() { let time = event.time(); if let Event::Bolus(bolus) = event { - let lagtime = fn_lag(&spp.clone().into(), time, covariates); + let lagtime = fn_lag(¶meters.clone().into(), time, covariates); if let Some(l) = lagtime.get(&bolus.input()) { *bolus.mut_time() += l; } @@ -609,12 +609,12 @@ impl Occasion { fn add_bioavailability(&mut self, reorder: Option<(&Fa, &Lag, &[f64], &Covariates)>) { // If lagtime is empty, return early - if let Some((fn_fa, _, spp, covariates)) = reorder { - let spp = nalgebra::DVector::from_vec(spp.to_vec()); + if let Some((fn_fa, _, parameters, covariates)) = reorder { + let parameters = nalgebra::DVector::from_vec(parameters.to_vec()); for event in self.events.iter_mut() { let time = event.time(); if let Event::Bolus(bolus) = event { - let fa = fn_fa(&spp.clone().into(), time, covariates); + let fa = fn_fa(¶meters.clone().into(), time, covariates); if let Some(f) = fa.get(&bolus.input()) { bolus.set_amount(bolus.amount() * f); } @@ -656,7 +656,7 @@ impl Occasion { /// /// # Arguments /// - /// * `reorder` - Optional tuple containing references to (Fa, Lag, support point, covariates) for adjustments + /// * `reorder` - Optional tuple containing references to (Fa, Lag, parameter, covariates) for adjustments /// * `ignore` - If true, filter out events marked as ignore /// * `mappings` - Optional reference to an [equation::Mapper] for input remapping /// diff --git a/src/lib.rs b/src/lib.rs index 9dd8e195..cf1c1b0a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -14,7 +14,7 @@ pub use crate::data::Interpolation::*; pub use crate::data::*; pub use crate::equation::*; pub use crate::optimize::effect::get_e2; -pub use crate::optimize::spp::SppOptimizer; +pub use crate::optimize::parameters::parametersOptimizer; pub use crate::simulator::equation::{ self, ode::{ExplicitRkTableau, OdeSolver, SdirkTableau}, diff --git a/src/optimize/mod.rs b/src/optimize/mod.rs index a37d90e7..14f2cd50 100644 --- a/src/optimize/mod.rs +++ b/src/optimize/mod.rs @@ -1,2 +1,2 @@ pub mod effect; -pub mod spp; +pub mod parameters; diff --git a/src/optimize/spp.rs b/src/optimize/parameters.rs similarity index 78% rename from src/optimize/spp.rs rename to src/optimize/parameters.rs index ddebb297..796cd5ef 100644 --- a/src/optimize/spp.rs +++ b/src/optimize/parameters.rs @@ -7,28 +7,28 @@ use ndarray::{Array1, Axis}; use crate::{prelude::simulator::log_likelihood_matrix, AssayErrorModels, Data, Equation}; -pub struct SppOptimizer<'a, E: Equation> { +pub struct parametersOptimizer<'a, E: Equation> { equation: &'a E, data: &'a Data, sig: &'a AssayErrorModels, pyl: &'a Array1, } -impl CostFunction for SppOptimizer<'_, E> { +impl CostFunction for parametersOptimizer<'_, E> { type Param = Vec; type Output = f64; - fn cost(&self, spp: &Self::Param) -> Result { - let theta = Array1::from(spp.clone()).insert_axis(Axis(0)); + fn cost(&self, parameters: &Self::Param) -> Result { + let theta = Array1::from(parameters.clone()).insert_axis(Axis(0)); let log_psi = log_likelihood_matrix(self.equation, self.data, &theta, self.sig, false)?; let psi = log_psi.mapv(f64::exp); if psi.ncols() > 1 { - tracing::error!("Psi in SppOptimizer has more than one column"); + tracing::error!("Psi in parametersOptimizer has more than one column"); } if psi.nrows() != self.pyl.len() { tracing::error!( - "Psi in SppOptimizer has {} rows, but spp has {}", + "Psi in parametersOptimizer has {} rows, but parameters has {}", psi.nrows(), self.pyl.len() ); @@ -42,7 +42,7 @@ impl CostFunction for SppOptimizer<'_, E> { } } -impl<'a, E: Equation> SppOptimizer<'a, E> { +impl<'a, E: Equation> parametersOptimizer<'a, E> { pub fn new( equation: &'a E, data: &'a Data, @@ -56,8 +56,8 @@ impl<'a, E: Equation> SppOptimizer<'a, E> { pyl, } } - pub fn optimize_point(self, spp: Array1) -> Result, Error> { - let simplex = create_initial_simplex(&spp.to_vec()); + pub fn optimize_point(self, parameters: Array1) -> Result, Error> { + let simplex = create_initial_simplex(¶meters.to_vec()); let solver: NelderMead, f64> = NelderMead::new(simplex).with_sd_tolerance(1e-2)?; let res = Executor::new(self, solver) .configure(|state| state.max_iters(5)) diff --git a/src/simulator/cache.rs b/src/simulator/cache.rs index a9d18311..eaa2f332 100644 --- a/src/simulator/cache.rs +++ b/src/simulator/cache.rs @@ -28,10 +28,10 @@ use crate::simulator::likelihood::SubjectPredictions; /// Default maximum number of entries per cache. pub const DEFAULT_CACHE_SIZE: u64 = 100_000; -/// Cache key: (subject_hash, support_point_hash) +/// Cache key: (subject_hash, parameters_hash) pub(crate) type PredictionKey = (u64, u64); -/// Cache key for SDE: (subject_hash, support_point_hash, error_model_hash) +/// Cache key for SDE: (subject_hash, parameters_hash, error_model_hash) pub(crate) type SdeKey = (u64, u64, u64); /// Thread-safe LRU cache for subject predictions. diff --git a/src/simulator/equation/analytical/mod.rs b/src/simulator/equation/analytical/mod.rs index 72376202..dd71b0ce 100644 --- a/src/simulator/equation/analytical/mod.rs +++ b/src/simulator/equation/analytical/mod.rs @@ -13,7 +13,7 @@ pub use three_compartment_models::*; pub use two_compartment_cl_models::*; pub use two_compartment_models::*; -use super::spphash; +use super::parametershash; use crate::data::error_model::AssayErrorModels; use crate::simulator::cache::{PredictionCache, DEFAULT_CACHE_SIZE}; @@ -151,13 +151,13 @@ impl EquationPriv for Analytical { // } // #[inline(always)] - // fn get_lag(&self, spp: &[f64]) -> Option> { - // Some((self.lag)(&V::from_vec(spp.to_owned()))) + // fn get_lag(&self, parameters: &[f64]) -> Option> { + // Some((self.lag)(&V::from_vec(parameters.to_owned()))) // } // #[inline(always)] - // fn get_fa(&self, spp: &[f64]) -> Option> { - // Some((self.fa)(&V::from_vec(spp.to_owned()))) + // fn get_fa(&self, parameters: &[f64]) -> Option> { + // Some((self.fa)(&V::from_vec(parameters.to_owned()))) // } #[inline(always)] @@ -188,7 +188,7 @@ impl EquationPriv for Analytical { fn solve( &self, x: &mut Self::S, - support_point: &[f64], + parameters: &[f64], covariates: &Covariates, infusions: &[Infusion], ti: f64, @@ -217,7 +217,7 @@ impl EquationPriv for Analytical { // 2) March over each sub-interval let mut current_t = ts[0]; - let mut sp = V::from_vec(support_point.to_vec(), NalgebraContext); + let mut sp = V::from_vec(parameters.to_vec(), NalgebraContext); let mut rateiv = V::zeros(self.get_ndrugs(), NalgebraContext); for &next_t in &ts[1..] { @@ -253,7 +253,7 @@ impl EquationPriv for Analytical { #[inline(always)] fn process_observation( &self, - support_point: &[f64], + parameters: &[f64], observation: &Observation, error_models: Option<&AssayErrorModels>, _time: f64, @@ -266,7 +266,7 @@ impl EquationPriv for Analytical { let out = &self.out; (out)( x, - &V::from_vec(support_point.to_vec(), NalgebraContext), + &V::from_vec(parameters.to_vec(), NalgebraContext), observation.time(), covariates, &mut y, @@ -280,12 +280,12 @@ impl EquationPriv for Analytical { Ok(()) } #[inline(always)] - fn initial_state(&self, spp: &[f64], covariates: &Covariates, occasion_index: usize) -> V { + fn initial_state(&self, parameters: &[f64], covariates: &Covariates, occasion_index: usize) -> V { let init = &self.init; let mut x = V::zeros(self.get_nstates(), NalgebraContext); if occasion_index == 0 { (init)( - &V::from_vec(spp.to_vec(), NalgebraContext), + &V::from_vec(parameters.to_vec(), NalgebraContext), 0.0, covariates, &mut x, @@ -554,19 +554,19 @@ impl Equation for Analytical { fn estimate_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - _estimate_likelihood(self, subject, support_point, error_models) + _estimate_likelihood(self, subject, parameters, error_models) } fn estimate_log_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - let ypred = _subject_predictions(self, subject, support_point)?; + let ypred = _subject_predictions(self, subject, parameters)?; ypred.log_likelihood(error_models) } @@ -579,28 +579,28 @@ impl Equation for Analytical { fn _subject_predictions( analytical: &Analytical, subject: &Subject, - support_point: &[f64], + parameters: &[f64], ) -> Result { if let Some(cache) = &analytical.cache { - let key = (subject.hash(), spphash(support_point)); + let key = (subject.hash(), parametershash(parameters)); if let Some(cached) = cache.get(&key) { return Ok(cached); } - let result = analytical.simulate_subject(subject, support_point, None)?.0; + let result = analytical.simulate_subject(subject, parameters, None)?.0; cache.insert(key, result.clone()); Ok(result) } else { - Ok(analytical.simulate_subject(subject, support_point, None)?.0) + Ok(analytical.simulate_subject(subject, parameters, None)?.0) } } fn _estimate_likelihood( ode: &Analytical, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - let ypred = _subject_predictions(ode, subject, support_point)?; + let ypred = _subject_predictions(ode, subject, parameters)?; Ok(ypred.log_likelihood(error_models)?.exp()) } diff --git a/src/simulator/equation/mod.rs b/src/simulator/equation/mod.rs index e462537c..c6c0955c 100644 --- a/src/simulator/equation/mod.rs +++ b/src/simulator/equation/mod.rs @@ -131,7 +131,7 @@ pub(crate) trait EquationPriv: EquationTypes { fn solve( &self, state: &mut Self::S, - support_point: &[f64], + parameters: &[f64], covariates: &Covariates, infusions: &[Infusion], start_time: f64, @@ -148,7 +148,7 @@ pub(crate) trait EquationPriv: EquationTypes { #[allow(clippy::too_many_arguments)] fn process_observation( &self, - support_point: &[f64], + parameters: &[f64], observation: &Observation, error_models: Option<&AssayErrorModels>, time: f64, @@ -160,7 +160,7 @@ pub(crate) trait EquationPriv: EquationTypes { fn initial_state( &self, - support_point: &[f64], + parameters: &[f64], covariates: &Covariates, occasion_index: usize, ) -> Self::S; @@ -168,7 +168,7 @@ pub(crate) trait EquationPriv: EquationTypes { #[allow(clippy::too_many_arguments)] fn simulate_event( &self, - support_point: &[f64], + parameters: &[f64], event: &Event, next_event: Option<&Event>, error_models: Option<&AssayErrorModels>, @@ -193,7 +193,7 @@ pub(crate) trait EquationPriv: EquationTypes { } Event::Observation(observation) => { self.process_observation( - support_point, + parameters, observation, error_models, event.time(), @@ -208,7 +208,7 @@ pub(crate) trait EquationPriv: EquationTypes { if let Some(next_event) = next_event { self.solve( x, - support_point, + parameters, covariates, infusions, event.time(), @@ -232,7 +232,7 @@ pub(crate) trait EquationPriv: EquationTypes { /// is provided for backward compatibility. #[allow(private_bounds)] pub trait Equation: EquationPriv + 'static + Clone + Sync { - /// Estimate the likelihood of the subject given the support point and error model. + /// Estimate the likelihood of the subject given the parameter and error model. /// /// **Deprecated**: Use [`estimate_log_likelihood`](Self::estimate_log_likelihood) instead /// for better numerical stability, especially with many observations or extreme parameter values. @@ -242,7 +242,7 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { /// /// # Parameters /// - `subject`: The subject data - /// - `support_point`: The parameter values + /// - `parameters`: The parameter values /// - `error_model`: The error model /// /// # Returns @@ -254,11 +254,11 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { fn estimate_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result; - /// Estimate the log-likelihood of the subject given the support point and error model. + /// Estimate the log-likelihood of the subject given the parameter and error model. /// /// This function calculates the log of how likely the observed data is given the model /// parameters and error model. It is numerically more stable than `estimate_likelihood` @@ -269,7 +269,7 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { /// /// # Parameters /// - `subject`: The subject data - /// - `support_point`: The parameter values + /// - `parameters`: The parameter values /// - `error_models`: The error model /// /// # Returns @@ -277,7 +277,7 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { fn estimate_log_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result; @@ -287,16 +287,16 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { /// /// # Parameters /// - `subject`: The subject data - /// - `support_point`: The parameter values + /// - `parameters`: The parameter values /// /// # Returns /// Predicted concentrations fn estimate_predictions( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], ) -> Result { - Ok(self.simulate_subject(subject, support_point, None)?.0) + Ok(self.simulate_subject(subject, parameters, None)?.0) } /// Get the number of output equations in the model. @@ -313,7 +313,7 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { /// /// # Parameters /// - `subject`: The subject data - /// - `support_point`: The parameter values + /// - `parameters`: The parameter values /// - `error_model`: The error model (optional) /// /// # Returns @@ -321,7 +321,7 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { fn simulate_subject( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: Option<&AssayErrorModels>, ) -> Result<(Self::P, Option), PharmsolError> { let mut output = Self::P::new(self.nparticles()); @@ -329,15 +329,13 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync { for occasion in subject.occasions() { let covariates = occasion.covariates(); - let mut x = self.initial_state(support_point, covariates, occasion.index()); + let mut x = self.initial_state(parameters, covariates, occasion.index()); let mut infusions = Vec::new(); - let events = occasion.process_events( - Some((self.fa(), self.lag(), support_point, covariates)), - true, - ); + let events = occasion + .process_events(Some((self.fa(), self.lag(), parameters, covariates)), true); for (index, event) in events.iter().enumerate() { self.simulate_event( - support_point, + parameters, event, events.get(index + 1), error_models, @@ -372,12 +370,12 @@ impl EqnKind { } } -/// Hash support points to a u64 for cache key generation. +/// Hash parameters to a u64 for cache key generation. #[inline(always)] -fn spphash(spp: &[f64]) -> u64 { +fn parametershash(parameters: &[f64]) -> u64 { use std::hash::{Hash, Hasher}; let mut hasher = ahash::AHasher::default(); - for &value in spp { + for &value in parameters { // Normalize -0.0 to 0.0 for consistent hashing let bits = if value == 0.0 { 0u64 } else { value.to_bits() }; bits.hash(&mut hasher); diff --git a/src/simulator/equation/ode/mod.rs b/src/simulator/equation/ode/mod.rs index ae69facc..c1343d16 100644 --- a/src/simulator/equation/ode/mod.rs +++ b/src/simulator/equation/ode/mod.rs @@ -8,7 +8,7 @@ use crate::{ Event, Observation, PharmsolError, Subject, }; -use super::spphash; +use super::parametershash; use crate::simulator::cache::{PredictionCache, DEFAULT_CACHE_SIZE}; use crate::simulator::equation::Predictions; use closure::PMProblem; @@ -161,10 +161,10 @@ impl State for V { fn _estimate_likelihood( ode: &ODE, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - let ypred = _subject_predictions(ode, subject, support_point)?; + let ypred = _subject_predictions(ode, subject, parameters)?; Ok(ypred.log_likelihood(error_models)?.exp()) } @@ -172,19 +172,19 @@ fn _estimate_likelihood( fn _subject_predictions( ode: &ODE, subject: &Subject, - support_point: &[f64], + parameters: &[f64], ) -> Result { if let Some(cache) = &ode.cache { - let key = (subject.hash(), spphash(support_point)); + let key = (subject.hash(), parametershash(parameters)); if let Some(cached) = cache.get(&key) { return Ok(cached); } - let result = ode.simulate_subject(subject, support_point, None)?.0; + let result = ode.simulate_subject(subject, parameters, None)?.0; cache.insert(key, result.clone()); Ok(result) } else { - Ok(ode.simulate_subject(subject, support_point, None)?.0) + Ok(ode.simulate_subject(subject, parameters, None)?.0) } } @@ -195,15 +195,15 @@ impl EquationTypes for ODE { impl EquationPriv for ODE { //#[inline(always)] - // fn get_lag(&self, spp: &[f64]) -> Option> { - // let spp = DVector::from_vec(spp.to_vec()); - // Some((self.lag)(&spp)) + // fn get_lag(&self, parameters: &[f64]) -> Option> { + // let parameters = DVector::from_vec(parameters.to_vec()); + // Some((self.lag)(¶meters)) // } // #[inline(always)] - // fn get_fa(&self, spp: &[f64]) -> Option> { - // let spp = DVector::from_vec(spp.to_vec()); - // Some((self.fa)(&spp)) + // fn get_fa(&self, parameters: &[f64]) -> Option> { + // let parameters = DVector::from_vec(parameters.to_vec()); + // Some((self.fa)(¶meters)) // } #[inline(always)] fn lag(&self) -> &Lag { @@ -232,7 +232,7 @@ impl EquationPriv for ODE { fn solve( &self, _state: &mut Self::S, - _support_point: &[f64], + _parameters: &[f64], _covariates: &Covariates, _infusions: &[Infusion], _start_time: f64, @@ -243,7 +243,7 @@ impl EquationPriv for ODE { #[inline(always)] fn process_observation( &self, - _support_point: &[f64], + _parameters: &[f64], _observation: &Observation, _error_models: Option<&AssayErrorModels>, _time: f64, @@ -256,12 +256,17 @@ impl EquationPriv for ODE { } #[inline(always)] - fn initial_state(&self, spp: &[f64], covariates: &Covariates, occasion_index: usize) -> V { + fn initial_state( + &self, + parameters: &[f64], + covariates: &Covariates, + occasion_index: usize, + ) -> V { let init = &self.init; let mut x = V::zeros(self.get_nstates(), NalgebraContext); if occasion_index == 0 { - let spp = DVector::from_vec(spp.to_vec()); - (init)(&spp.into(), 0.0, covariates, &mut x); + let parameters = DVector::from_vec(parameters.to_vec()); + (init)(¶meters.into(), 0.0, covariates, &mut x); } x } @@ -274,7 +279,7 @@ impl ODE { &self, solver: &mut S, events: &[Event], - spp_v: &V, + parameters_v: &V, covariates: &Covariates, error_models: Option<&AssayErrorModels>, bolus_v: &mut V, @@ -305,7 +310,7 @@ impl ODE { (self.diffeq)( solver.state().y, - spp_v, + parameters_v, event.time(), state_without_bolus, zero_bolus, @@ -315,7 +320,7 @@ impl ODE { (self.diffeq)( solver.state().y, - spp_v, + parameters_v, event.time(), state_with_bolus, bolus_v, @@ -333,7 +338,7 @@ impl ODE { y_out.fill(0.0); (self.out)( solver.state().y, - spp_v, + parameters_v, observation.time(), covariates, y_out, @@ -394,19 +399,19 @@ impl Equation for ODE { fn estimate_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - _estimate_likelihood(self, subject, support_point, error_models) + _estimate_likelihood(self, subject, parameters, error_models) } fn estimate_log_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - let ypred = _subject_predictions(self, subject, support_point)?; + let ypred = _subject_predictions(self, subject, parameters)?; ypred.log_likelihood(error_models) } @@ -417,7 +422,7 @@ impl Equation for ODE { fn simulate_subject( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: Option<&AssayErrorModels>, ) -> Result<(Self::P, Option), PharmsolError> { let mut output = Self::P::new(self.nparticles()); @@ -436,7 +441,7 @@ impl Equation for ODE { let zero_bolus = V::zeros(ndrugs, NalgebraContext); let zero_rateiv = V::zeros(ndrugs, NalgebraContext); let mut bolus_v = V::zeros(ndrugs, NalgebraContext); - let spp_v: V = DVector::from_vec(support_point.to_vec()).into(); + let parameters_v: V = DVector::from_vec(parameters.to_vec()).into(); // Pre-allocate output vector for observations let mut y_out = V::zeros(self.get_nouteqs(), NalgebraContext); @@ -445,26 +450,24 @@ impl Equation for ODE { for occasion in subject.occasions() { let covariates = occasion.covariates(); let infusions = occasion.infusions_ref(); - let events = occasion.process_events( - Some((self.fa(), self.lag(), support_point, covariates)), - true, - ); + let events = occasion + .process_events(Some((self.fa(), self.lag(), parameters, covariates)), true); let problem = OdeBuilder::::new() .atol(vec![self.atol]) .rtol(self.rtol) .t0(occasion.initial_time()) .h0(1e-3) - .p(support_point.to_vec()) + .p(parameters.to_vec()) .build_from_eqn(PMProblem::with_params_v( self.diffeq, nstates, ndrugs, - support_point.to_vec(), - spp_v.clone(), + parameters.to_vec(), + parameters_v.clone(), covariates, infusions.as_slice(), - self.initial_state(support_point, covariates, occasion.index()) + self.initial_state(parameters, covariates, occasion.index()) .into(), )?)?; @@ -475,7 +478,7 @@ impl Equation for ODE { self, &mut solver, &events, - &spp_v, + ¶meters_v, covariates, error_models, &mut bolus_v, @@ -494,7 +497,7 @@ impl Equation for ODE { self, &mut solver, &events, - &spp_v, + ¶meters_v, covariates, error_models, &mut bolus_v, @@ -513,7 +516,7 @@ impl Equation for ODE { self, &mut solver, &events, - &spp_v, + ¶meters_v, covariates, error_models, &mut bolus_v, @@ -532,7 +535,7 @@ impl Equation for ODE { self, &mut solver, &events, - &spp_v, + ¶meters_v, covariates, error_models, &mut bolus_v, diff --git a/src/simulator/equation/sde/mod.rs b/src/simulator/equation/sde/mod.rs index 704165e5..e9bcb87f 100644 --- a/src/simulator/equation/sde/mod.rs +++ b/src/simulator/equation/sde/mod.rs @@ -14,7 +14,7 @@ use crate::{ Subject, }; -use super::spphash; +use super::parametershash; use crate::simulator::cache::{SdeLikelihoodCache, DEFAULT_CACHE_SIZE}; use diffsol::VectorCommon; @@ -33,7 +33,7 @@ use super::{Equation, EquationPriv, EquationTypes, Predictions, State}; /// * `drift` - Function defining the deterministic component of the SDE /// * `difussion` - Function defining the stochastic component of the SDE /// * `x` - Current state vector -/// * `support_point` - Parameter vector for the model +/// * `parameters` - Parameter vector for the model /// * `cov` - Covariates that may influence the system dynamics /// * `infusions` - Infusion events to be applied during simulation /// * `ti` - Starting time @@ -47,7 +47,7 @@ pub(crate) fn simulate_sde_event( drift: &Drift, difussion: &Diffusion, x: V, - support_point: &[f64], + parameters: &[f64], cov: &Covariates, infusions: &[Infusion], ndrugs: usize, @@ -61,7 +61,7 @@ pub(crate) fn simulate_sde_event( let mut sde = em::EM::new( *drift, *difussion, - DVector::from_column_slice(support_point), + DVector::from_column_slice(parameters), x.inner().clone(), cov.clone(), infusions.to_vec(), @@ -252,13 +252,13 @@ impl EquationPriv for SDE { // } // #[inline(always)] - // fn get_lag(&self, spp: &[f64]) -> Option> { - // Some((self.lag)(&V::from_vec(spp.to_owned()))) + // fn get_lag(&self, parameters: &[f64]) -> Option> { + // Some((self.lag)(&V::from_vec(parameters.to_owned()))) // } // #[inline(always)] - // fn get_fa(&self, spp: &[f64]) -> Option> { - // Some((self.fa)(&V::from_vec(spp.to_owned()))) + // fn get_fa(&self, parameters: &[f64]) -> Option> { + // Some((self.fa)(&V::from_vec(parameters.to_owned()))) // } #[inline(always)] @@ -289,7 +289,7 @@ impl EquationPriv for SDE { fn solve( &self, state: &mut Self::S, - support_point: &[f64], + parameters: &[f64], covariates: &Covariates, infusions: &[Infusion], ti: f64, @@ -301,7 +301,7 @@ impl EquationPriv for SDE { &self.drift, &self.diffusion, particle.clone().into(), - support_point, + parameters, covariates, infusions, ndrugs, @@ -323,7 +323,7 @@ impl EquationPriv for SDE { #[inline(always)] fn process_observation( &self, - support_point: &[f64], + parameters: &[f64], observation: &crate::Observation, error_models: Option<&AssayErrorModels>, _time: f64, @@ -338,7 +338,7 @@ impl EquationPriv for SDE { let mut y = V::zeros(self.get_nouteqs(), NalgebraContext); (self.out)( &x[i].clone().into(), - &V::from_vec(support_point.to_vec(), NalgebraContext), + &V::from_vec(parameters.to_vec(), NalgebraContext), observation.time(), covariates, &mut y, @@ -373,7 +373,7 @@ impl EquationPriv for SDE { #[inline(always)] fn initial_state( &self, - support_point: &[f64], + parameters: &[f64], covariates: &Covariates, occasion_index: usize, ) -> Self::S { @@ -382,7 +382,7 @@ impl EquationPriv for SDE { let mut state: V = DVector::zeros(self.get_nstates()).into(); if occasion_index == 0 { (self.init)( - &V::from_vec(support_point.to_vec(), NalgebraContext), + &V::from_vec(parameters.to_vec(), NalgebraContext), 0.0, covariates, &mut state, @@ -400,7 +400,7 @@ impl Equation for SDE { /// # Arguments /// /// * `subject` - Subject data containing observations - /// * `support_point` - Parameter vector for the model + /// * `parameters` - Parameter vector for the model /// * `error_model` - Error model to use for likelihood calculations /// /// # Returns @@ -409,21 +409,21 @@ impl Equation for SDE { fn estimate_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { - _estimate_likelihood(self, subject, support_point, error_models) + _estimate_likelihood(self, subject, parameters, error_models) } fn estimate_log_likelihood( &self, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { // For SDE, the particle filter computes likelihood in regular space. // We compute it directly and then take the log. - let lik = _estimate_likelihood(self, subject, support_point, error_models)?; + let lik = _estimate_likelihood(self, subject, parameters, error_models)?; if lik > 0.0 { Ok(lik.ln()) @@ -441,21 +441,25 @@ impl Equation for SDE { fn _estimate_likelihood( sde: &SDE, subject: &Subject, - support_point: &[f64], + parameters: &[f64], error_models: &AssayErrorModels, ) -> Result { if let Some(cache) = &sde.cache { - let key = (subject.hash(), spphash(support_point), error_models.hash()); + let key = ( + subject.hash(), + parametershash(parameters), + error_models.hash(), + ); if let Some(cached) = cache.get(&key) { return Ok(cached); } - let ypred = sde.simulate_subject(subject, support_point, Some(error_models))?; + let ypred = sde.simulate_subject(subject, parameters, Some(error_models))?; let result = ypred.1.unwrap(); cache.insert(key, result); Ok(result) } else { - let ypred = sde.simulate_subject(subject, support_point, Some(error_models))?; + let ypred = sde.simulate_subject(subject, parameters, Some(error_models))?; Ok(ypred.1.unwrap()) } } diff --git a/src/simulator/likelihood/matrix.rs b/src/simulator/likelihood/matrix.rs index fdb3e5c1..aa25f134 100644 --- a/src/simulator/likelihood/matrix.rs +++ b/src/simulator/likelihood/matrix.rs @@ -1,7 +1,7 @@ //! Population-level log-likelihood matrix computation. //! //! This module provides functions for computing log-likelihood matrices -//! across populations of subjects and parameter support points. +//! across populations of subjects and parameter parameters. use ndarray::{Array2, Axis, ShapeBuilder}; use rayon::prelude::*; @@ -11,7 +11,7 @@ use crate::{Data, Equation, PharmsolError}; use super::progress::ProgressTracker; -/// Calculate the log-likelihood matrix for all subjects and support points. +/// Calculate the log-likelihood matrix for all subjects and parameters. /// /// This function computes log-likelihoods directly in log-space, which is numerically /// more stable than computing likelihoods and then taking logarithms. This is especially @@ -21,12 +21,12 @@ use super::progress::ProgressTracker; /// # Parameters /// - `equation`: The equation to use for simulation /// - `subjects`: The subject data -/// - `support_points`: The support points to evaluate (rows = support points, cols = parameters) +/// - `parameters`: The parameters to evaluate (rows = parameters, cols = parameters) /// - `error_models`: The error models to use (observation-based sigma) /// - `progress`: Whether to display a progress bar during computation` /// /// # Returns -/// A 2D array of log-likelihoods with shape (n_subjects, n_support_points) +/// A 2D array of log-likelihoods with shape (n_subjects, n_parameters) /// /// # Example /// ```ignore @@ -35,7 +35,7 @@ use super::progress::ProgressTracker; /// let log_liks = log_likelihood_matrix( /// &equation, /// &data, -/// &support_points, +/// ¶meters, /// &error_models, /// false /// )?; @@ -43,20 +43,20 @@ use super::progress::ProgressTracker; pub fn log_likelihood_matrix( equation: &impl Equation, subjects: &Data, - support_points: &Array2, + parameters: &Array2, error_models: &AssayErrorModels, progress: bool, ) -> Result, PharmsolError> { - let mut log_psi: Array2 = Array2::default((subjects.len(), support_points.nrows()).f()); + let mut log_psi: Array2 = Array2::default((subjects.len(), parameters.nrows()).f()); let subjects_vec = subjects.subjects(); let progress_tracker = if progress { - let total = subjects_vec.len() * support_points.nrows(); + let total = subjects_vec.len() * parameters.nrows(); println!( - "Computing log-likelihood matrix: {} subjects × {} support points...", + "Computing log-likelihood matrix: {} subjects × {} parameters...", subjects_vec.len(), - support_points.nrows() + parameters.nrows() ); Some(ProgressTracker::new(total)) } else { @@ -75,7 +75,7 @@ pub fn log_likelihood_matrix( let subject = subjects_vec.get(i).unwrap(); match equation.estimate_log_likelihood( subject, - &support_points.row(j).to_vec(), + ¶meters.row(j).to_vec(), error_models, ) { Ok(log_likelihood) => { @@ -110,11 +110,11 @@ pub fn log_likelihood_matrix( pub fn log_psi( equation: &impl Equation, subjects: &Data, - support_points: &Array2, + parameters: &Array2, error_models: &AssayErrorModels, progress: bool, ) -> Result, PharmsolError> { - log_likelihood_matrix(equation, subjects, support_points, error_models, progress) + log_likelihood_matrix(equation, subjects, parameters, error_models, progress) } /// Calculate the likelihood matrix (deprecated). @@ -132,12 +132,12 @@ pub fn log_psi( pub fn psi( equation: &impl Equation, subjects: &Data, - support_points: &Array2, + parameters: &Array2, error_models: &AssayErrorModels, progress: bool, ) -> Result, PharmsolError> { let log_psi_matrix = - log_likelihood_matrix(equation, subjects, support_points, error_models, progress)?; + log_likelihood_matrix(equation, subjects, parameters, error_models, progress)?; // Exponentiate to get likelihoods (may underflow to 0 for extreme values) Ok(log_psi_matrix.mapv(f64::exp)) diff --git a/src/simulator/likelihood/mod.rs b/src/simulator/likelihood/mod.rs index c703dee7..09df9e28 100644 --- a/src/simulator/likelihood/mod.rs +++ b/src/simulator/likelihood/mod.rs @@ -18,7 +18,7 @@ //! ## For Non-Parametric Algorithms //! //! Use [`log_likelihood_matrix`] to compute a matrix of log-likelihoods across -//! all subjects and support points: +//! all subjects and parameters: //! //! ```ignore //! use pharmsol::prelude::simulator::{log_likelihood_matrix, LikelihoodMatrixOptions}; @@ -26,7 +26,7 @@ //! let log_liks = log_likelihood_matrix( //! &equation, //! &data, -//! &support_points, +//! ¶meters, //! &error_models, //! LikelihoodMatrixOptions::new().with_progress(), //! )?; diff --git a/src/simulator/likelihood/subject.rs b/src/simulator/likelihood/subject.rs index 77d8963b..e1ec07b7 100644 --- a/src/simulator/likelihood/subject.rs +++ b/src/simulator/likelihood/subject.rs @@ -138,7 +138,7 @@ impl From> for SubjectPredictions { /// Container for predictions across a population of subjects. /// /// This struct holds predictions for multiple subjects organized in a 2D array -/// where rows represent subjects and columns represent support points (or +/// where rows represent subjects and columns represent parameters (or /// other groupings). pub struct PopulationPredictions { /// 2D array of subject predictions diff --git a/tests/test_solvers.rs b/tests/test_solvers.rs index 2ce08f3a..091b53e9 100644 --- a/tests/test_solvers.rs +++ b/tests/test_solvers.rs @@ -40,9 +40,9 @@ fn one_cpt(solver: OdeSolver) -> equation::ODE { fn preds(solver: OdeSolver) -> Vec { let sub = subject(); - let spp = vec![0.1, 50.0]; + let parameters = vec![0.1, 50.0]; one_cpt(solver) - .estimate_predictions(&sub, &spp) + .estimate_predictions(&sub, ¶meters) .unwrap() .flat_predictions() .to_vec() From 2e712879e691b9c2b9bd53200783e91a3504b17f Mon Sep 17 00:00:00 2001 From: Markus <66058642+mhovd@users.noreply.github.com> Date: Fri, 24 Apr 2026 14:16:58 +0200 Subject: [PATCH 2/2] Format --- src/simulator/equation/analytical/mod.rs | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/simulator/equation/analytical/mod.rs b/src/simulator/equation/analytical/mod.rs index dd71b0ce..3b485092 100644 --- a/src/simulator/equation/analytical/mod.rs +++ b/src/simulator/equation/analytical/mod.rs @@ -280,7 +280,12 @@ impl EquationPriv for Analytical { Ok(()) } #[inline(always)] - fn initial_state(&self, parameters: &[f64], covariates: &Covariates, occasion_index: usize) -> V { + fn initial_state( + &self, + parameters: &[f64], + covariates: &Covariates, + occasion_index: usize, + ) -> V { let init = &self.init; let mut x = V::zeros(self.get_nstates(), NalgebraContext); if occasion_index == 0 {