From 19e4655748b16b60c124f7200845f63fe33b8e86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Tue, 8 Oct 2024 00:31:59 +0100 Subject: [PATCH 1/2] support for addl --- src/data/event.rs | 7 ++ src/data/parse_pmetrics.rs | 178 +++++++++++++++++++++++------------ src/tests/data/addl_test.csv | 5 + 3 files changed, 128 insertions(+), 62 deletions(-) create mode 100644 src/tests/data/addl_test.csv diff --git a/src/data/event.rs b/src/data/event.rs index 2c30882..33ddaa6 100644 --- a/src/data/event.rs +++ b/src/data/event.rs @@ -17,6 +17,13 @@ impl Event { Event::Observation(observation) => observation.time, } } + pub(crate) fn inc_time(&mut self, dt: f64) { + match self { + Event::Bolus(bolus) => bolus.time += dt, + Event::Infusion(infusion) => infusion.time += dt, + Event::Observation(observation) => observation.time += dt, + } + } } /// An instantaenous input of drug diff --git a/src/data/parse_pmetrics.rs b/src/data/parse_pmetrics.rs index 239a5e9..dbdfaf2 100644 --- a/src/data/parse_pmetrics.rs +++ b/src/data/parse_pmetrics.rs @@ -97,8 +97,14 @@ pub fn read_pmetrics(path: impl Into) -> Result { // Parse events for row in rows.clone() { - let event: Event = Event::try_from(row.clone())?; - events.push(event); + match row.parse_events() { + Ok(ev) => events.extend(ev), + Err(e) => { + // dbg!(&row); + // dbg!(&e); + return Err(e); + } + } } // Parse covariates @@ -191,7 +197,8 @@ pub fn read_pmetrics(path: impl Into) -> Result { covariates.add_covariate(name, covariate) } // Create the block - let occasion = Occasion::new(events, covariates, block_index); + let mut occasion = Occasion::new(events, covariates, block_index); + occasion.sort(); occasions.push(occasion); } @@ -226,8 +233,8 @@ struct Row { #[serde(deserialize_with = "deserialize_option_isize")] addl: Option, /// Dosing interval - #[serde(deserialize_with = "deserialize_option_isize")] - ii: Option, + #[serde(deserialize_with = "deserialize_option_f64")] + ii: Option, /// Input compartment #[serde(deserialize_with = "deserialize_option_usize")] input: Option, @@ -262,68 +269,83 @@ impl Row { _ => None, } } -} - -impl TryFrom for Event { - type Error = PmetricsError; + fn parse_events(self) -> Result, PmetricsError> { + let mut events: Vec = Vec::new(); - fn try_from(row: Row) -> Result { - match row.evid { - 0 => Ok(Event::Observation(Observation::new( - row.time, - row.out + match self.evid { + 0 => events.push(Event::Observation(Observation::new( + self.time, + self.out .ok_or_else(|| PmetricsError::MissingObservationOut { - id: row.id.clone(), - time: row.time, + id: self.id.clone(), + time: self.time, })?, - row.outeq + self.outeq .ok_or_else(|| PmetricsError::MissingObservationOuteq { - id: row.id.clone(), - time: row.time, + id: self.id.clone(), + time: self.time, })? - 1, - row.get_errorpoly(), - row.out == Some(-99.0), + self.get_errorpoly(), + self.out == Some(-99.0), ))), 1 | 4 => { - if row.dur.unwrap_or(0.0) > 0.0 { - Ok(Event::Infusion(Infusion::new( - row.time, - row.dose.ok_or_else(|| PmetricsError::MissingInfusionDose { - id: row.id.clone(), - time: row.time, - })?, - row.input + let event = if self.dur.unwrap_or(0.0) > 0.0 { + Event::Infusion(Infusion::new( + self.time, + self.dose + .ok_or_else(|| PmetricsError::MissingInfusionDose { + id: self.id.clone(), + time: self.time, + })?, + self.input .ok_or_else(|| PmetricsError::MissingInfusionInput { - id: row.id.clone(), - time: row.time, + id: self.id.clone(), + time: self.time, })? - 1, - row.dur.ok_or_else(|| PmetricsError::MissingInfusionDur { - id: row.id.clone(), - time: row.time, + self.dur.ok_or_else(|| PmetricsError::MissingInfusionDur { + id: self.id.clone(), + time: self.time, })?, - ))) + )) } else { - Ok(Event::Bolus(Bolus::new( - row.time, - row.dose.ok_or_else(|| PmetricsError::MissingBolusDose { - id: row.id.clone(), - time: row.time, + Event::Bolus(Bolus::new( + self.time, + self.dose.ok_or_else(|| PmetricsError::MissingBolusDose { + id: self.id.clone(), + time: self.time, })?, - row.input.ok_or_else(|| PmetricsError::MissingBolusInput { - id: row.id, - time: row.time, + self.input.ok_or_else(|| PmetricsError::MissingBolusInput { + id: self.id, + time: self.time, })? - 1, - ))) + )) + }; + if self.addl.is_some() && self.ii.is_some() { + if self.addl.unwrap_or(0) != 0 && self.ii.unwrap_or(0.0) > 0.0 { + let mut ev = event.clone(); + let interval = &self.ii.unwrap().abs(); + let repetitions = &self.addl.unwrap().abs(); + let direction = &self.addl.unwrap().signum(); + + for _ in 0..*repetitions { + ev.inc_time((*direction as f64) * interval); + events.push(ev.clone()); + } + } } + events.push(event); } - _ => Err(PmetricsError::UnknownEvid { - evid: row.evid, - id: row.id.clone(), - time: row.time, - }), - } + _ => { + return Err(PmetricsError::UnknownEvid { + evid: self.evid, + id: self.id.clone(), + time: self.time, + }); + } + }; + Ok(events) } } @@ -409,18 +431,50 @@ where deserializer.deserialize_map(CovsVisitor) } -// #[cfg(test)] -// mod tests { -// use super::*; +#[cfg(test)] +mod tests { + + use super::*; + + #[test] + fn test_addl() { + let data = read_pmetrics("src/tests/data/addl_test.csv"); + + assert!(data.is_ok(), "Failed to parse data"); + + let data = data.unwrap(); + let subjects = data.get_subjects(); + let first_subject = subjects.get(0).unwrap(); + let second_subject = subjects.get(1).unwrap(); + let s1_occasions = first_subject.occasions(); + let s2_occasions = second_subject.occasions(); + let first_scenario = s1_occasions.get(0).unwrap(); + let second_scenario = s2_occasions.get(0).unwrap(); + + let s1_times = first_scenario + .events() + .iter() + .map(|e| e.get_time()) + .collect::>(); + + // Negative ADDL, observations shifted forward + + assert_eq!( + s1_times, + vec![-120.0, -108.0, -96.0, -84.0, -72.0, -60.0, -48.0, -36.0, -24.0, -12.0, 0.0, 9.0] + ); -// #[test] -// fn test_read_pmetrics() { -// let path = std::path::Path::new("examples/data/bimodal_ke_blocks.csv"); -// let data = read_pmetrics(&path).unwrap(); + let s2_times = second_scenario + .events() + .iter() + .map(|e| e.get_time()) + .collect::>(); -// assert_eq!(data.nsubjects(), 1); -// assert_eq!(data.nobs(), 30); + // Positive ADDL, no shift in observations -// println!("{}", data); -// } -// } + assert_eq!( + s2_times, + vec![0.0, 9.0, 12.0, 24.0, 36.0, 48.0, 60.0, 72.0, 84.0, 96.0, 108.0, 120.0] + ); + } +} diff --git a/src/tests/data/addl_test.csv b/src/tests/data/addl_test.csv new file mode 100644 index 0000000..d059e00 --- /dev/null +++ b/src/tests/data/addl_test.csv @@ -0,0 +1,5 @@ +ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3 +1,1,0,0,600,-10,12,1,.,.,.,.,.,. +1,0,9,.,.,.,.,.,100,100,.,.,.,. +2,1,0,0,600,10,12,1,.,.,.,.,.,. +2,0,9,.,.,.,.,.,100,100,.,.,.,. \ No newline at end of file From d20e9a5b02b679bf195e3b0d213deb45eed874e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Wed, 9 Oct 2024 12:20:48 +0100 Subject: [PATCH 2/2] fixed a bug that caused the program to generate NAs if two covariates where given at the same time --- src/data/parse_pmetrics.rs | 34 +++++++++++++++++---------- src/simulator/equation/ode/closure.rs | 1 + src/simulator/equation/ode/mod.rs | 3 ++- 3 files changed, 25 insertions(+), 13 deletions(-) diff --git a/src/data/parse_pmetrics.rs b/src/data/parse_pmetrics.rs index dbdfaf2..55f52c7 100644 --- a/src/data/parse_pmetrics.rs +++ b/src/data/parse_pmetrics.rs @@ -168,19 +168,29 @@ pub fn read_pmetrics(path: impl Into) -> Result { value: value.unwrap(), }, )); - } else if let Some(next) = next_occurrence { - // Linear interpolation for non-fixed covariates - let (next_time, next_value) = next; + } else if let Some((next_time, next_value)) = next_occurrence { if let Some(current_value) = value { - let slope = (next_value.unwrap() - current_value) / (next_time - time); - covariate.add_segment(CovariateSegment::new( - time, - *next_time, - InterpolationMethod::Linear { - slope, - intercept: current_value - slope * time, - }, - )); + if *next_time == time { + covariate.add_segment(CovariateSegment::new( + time, + *next_time, + InterpolationMethod::CarryForward { + value: current_value, + }, + )); + } else { + let slope = + (next_value.unwrap() - current_value) / (next_time - time); + covariate.add_segment(CovariateSegment::new( + time, + *next_time, + InterpolationMethod::Linear { + slope, + intercept: current_value - slope * time, + }, + )); + } + last_value = Some((next_time, next_value)); } } else if let Some((last_time, last_value)) = last_value { diff --git a/src/simulator/equation/ode/closure.rs b/src/simulator/equation/ode/closure.rs index 92cd7a5..ede1d30 100644 --- a/src/simulator/equation/ode/closure.rs +++ b/src/simulator/equation/ode/closure.rs @@ -5,6 +5,7 @@ use diffsol::{ op::{NonLinearOp, Op, OpStatistics}, vector::Vector, }; + use std::{cell::RefCell, rc::Rc}; use crate::data::{Covariates, Infusion}; diff --git a/src/simulator/equation/ode/mod.rs b/src/simulator/equation/ode/mod.rs index b10c35e..f7b60e3 100644 --- a/src/simulator/equation/ode/mod.rs +++ b/src/simulator/equation/ode/mod.rs @@ -135,9 +135,10 @@ impl EquationPriv for ODE { start_time: f64, end_time: f64, ) { - if start_time == end_time { + if f64::abs(start_time - end_time) < 1e-8 { return; } + // dbg!(start_time, end_time); let problem = build_pm_ode::( self.diffeq, |_p: &V, _t: T| state.clone(),