From 6c89a476e3f307fa8254fb08cc7ac84202b428ad Mon Sep 17 00:00:00 2001 From: Markus <66058642+mhovd@users.noreply.github.com> Date: Fri, 16 Oct 2026 12:02:02 +0200 Subject: [PATCH 1/2] fix: cover negative times in Data::expand grid The expand grid started at 0 and keyed times as u64 microseconds, so negative times produced no grid and collided at 0. Use signed microseconds and start at the first event when it is before 0, then continue on multiples of idelta from 0. Also fix two clippy lints flagged by a newer toolchain. Co-Authored-By: Claude Opus 5.5 --- src/data/parser/pmetrics/tests.rs | 4 +- src/data/structs.rs | 171 ++++++++++++++++++++++++++++-- src/optimize/parameters.rs | 1 - 3 files changed, 163 insertions(+), 13 deletions(-) diff --git a/src/data/parser/pmetrics/tests.rs b/src/data/parser/pmetrics/tests.rs index 619c29a0..0191d06f 100644 --- a/src/data/parser/pmetrics/tests.rs +++ b/src/data/parser/pmetrics/tests.rs @@ -906,7 +906,7 @@ fn programmatic_covariate_name_is_lowercased_for_export() { .build()]); let bytes = data.to_pmetrics_csv_bytes().unwrap(); let mut reader = ::csv::Reader::from_reader(bytes.as_slice()); - assert_eq!(reader.headers().unwrap().iter().last(), Some("wt")); + assert_eq!(reader.headers().unwrap().iter().next_back(), Some("wt")); let parsed = Data::from_pmetrics_csv_bytes(&bytes).unwrap(); let weight = parsed.subjects()[0].occasions()[0] @@ -927,7 +927,7 @@ fn covariate_key_and_name_may_differ_by_ascii_case() { let bytes = Data::new(vec![subject]).to_pmetrics_csv_bytes().unwrap(); let mut reader = ::csv::Reader::from_reader(bytes.as_slice()); - assert_eq!(reader.headers().unwrap().iter().last(), Some("wt")); + assert_eq!(reader.headers().unwrap().iter().next_back(), Some("wt")); } #[test] diff --git a/src/data/structs.rs b/src/data/structs.rs index 40434cb5..5d7d03f0 100644 --- a/src/data/structs.rs +++ b/src/data/structs.rs @@ -149,6 +149,13 @@ impl Data { /// to zero and the dataset is returned unchanged. /// * `tad` - Additional time to add after the last dose (time after dose) /// + /// The grid for each occasion starts at 0. If the occasion's earliest event + /// (dose or observation) is before 0, the grid instead starts at that event + /// and continues from the next multiple of `idelta`, so the simulation + /// still starts at the first real event and every later grid time is a + /// multiple of `idelta` measured from t = 0. The grid ends at the last dose + /// (or infusion end) of the occasion plus `tad`, inclusive. + /// /// # Returns /// /// A new `Data` object with expanded observations @@ -160,7 +167,7 @@ impl Data { // Work in integer microseconds so the grid always makes progress. // A sub-microsecond `idelta` would otherwise round back to the previous // time on every iteration and spin forever. - let step_us = (idelta * 1e6).round() as u64; + let step_us = (idelta * 1e6).round() as i64; if step_us == 0 { return self.clone(); } @@ -194,14 +201,27 @@ impl Data { .unwrap_or(0.0) + tad; + // Earliest event of this occasion. The grid starts here when it is + // before 0, so dosing histories with negative times are covered. + let first_time = occasion + .events + .iter() + .map(|event| match event { + Event::Bolus(bolus) => bolus.time(), + Event::Infusion(infusion) => infusion.time(), + Event::Observation(obs) => obs.time(), + }) + .min_by(|a, b| a.partial_cmp(b).unwrap()) + .unwrap_or(0.0); + // Create a set of existing (time, outeq) pairs for fast lookup - let existing_obs: std::collections::HashSet<(u64, OutputLabel)> = + let existing_obs: std::collections::HashSet<(i64, OutputLabel)> = old_events .iter() .filter_map(|event| match event { Event::Observation(obs) => { // Convert to microseconds for consistent comparison - let time_key = (obs.time() * 1e6).round() as u64; + let time_key = (obs.time() * 1e6).round() as i64; Some((time_key, obs.outeq().clone())) } _ => None, @@ -211,16 +231,12 @@ impl Data { // Generate new observation times, stepping in integer // microseconds to guarantee forward progress. let mut new_events = Vec::new(); - let last_time_us = (last_time * 1e6).round() as u64; - let mut time_key = 0u64; - while time_key <= last_time_us { - let time = time_key as f64 / 1e6; - + let mut add_grid_point = |time_key: i64| { for outeq in &outeq_values { // Only add if this (time, outeq) combination doesn't exist if !existing_obs.contains(&(time_key, outeq.clone())) { let obs = Observation::new( - time, + time_key as f64 / 1e6, None, outeq.clone(), None, @@ -230,7 +246,26 @@ impl Data { new_events.push(Event::Observation(obs)); } } - + }; + + // Signed microseconds, so times before 0 neither saturate nor + // collide. When the first event is before 0, the grid starts at + // that event (so the simulation does not start earlier than the + // data) and continues from the next multiple of `step_us`, to + // stay aligned with the grid after 0. + let last_time_us = (last_time * 1e6).round() as i64; + let first_time_us = (first_time * 1e6).round() as i64; + let mut time_key = if first_time_us < 0 { + let aligned = -(-first_time_us).div_euclid(step_us) * step_us; + if aligned != first_time_us && first_time_us <= last_time_us { + add_grid_point(first_time_us); + } + aligned + } else { + 0 + }; + while time_key <= last_time_us { + add_grid_point(time_key); time_key += step_us; } @@ -1736,4 +1771,120 @@ mod tests { // Occasion 1: last dose = 10, grid at 0,5,10 -> 2 generated + 1 real = 3 assert_eq!(count_obs(&subject.occasions()[1]), 3); } + + fn sorted_obs_times(occasion: &Occasion) -> Vec { + let mut times: Vec = occasion + .events() + .iter() + .filter_map(|e| match e { + Event::Observation(o) => Some(o.time()), + _ => None, + }) + .collect(); + times.sort_by(|a, b| a.partial_cmp(b).unwrap()); + times + } + + #[test] + fn expand_grid_covers_negative_times() { + // Dosing history before t = 0 with a trough sample just before the + // dose at 0; the grid ends at that last dose + tad. + let subject = Subject::builder("s1") + .bolus(-24.0, 100.0, 0) + .observation(-0.5, 5.0, 0) + .bolus(0.0, 100.0, 0) + .build(); + let data = Data::from(subject); + + let expanded = data.expand(0.5, 12.0); + let occasion = &expanded.subjects()[0].occasions()[0]; + let times = sorted_obs_times(occasion); + + assert_eq!(times.first(), Some(&-24.0)); + assert_eq!(times.last(), Some(&12.0)); + // -24..=12 in steps of 0.5 + assert_eq!(times.len(), 73); + assert!(times.windows(2).all(|w| w[0] < w[1]), "no duplicate times"); + + // -0.5 appears once, and it is the original observation + let at_minus_half: Vec<_> = occasion + .events() + .iter() + .filter_map(|e| match e { + Event::Observation(o) if o.time() == -0.5 => Some(o), + _ => None, + }) + .collect(); + assert_eq!(at_minus_half.len(), 1); + assert_eq!(at_minus_half[0].value(), Some(5.0)); + + assert!(times.contains(&0.0)); + } + + #[test] + fn expand_negative_observation_keeps_grid_point_at_zero() { + let subject = Subject::builder("s1") + .bolus(0.0, 100.0, 0) + .observation(-0.5, 5.0, 0) + .build(); + let data = Data::from(subject); + + let expanded = data.expand(1.0, 2.0); + let times = sorted_obs_times(&expanded.subjects()[0].occasions()[0]); + + // Grid starts at the observation at -0.5, then 0 is still added + assert_eq!(times, vec![-0.5, 0.0, 1.0, 2.0]); + } + + #[test] + fn expand_negative_grid_starts_at_first_event_then_aligns() { + let subject = Subject::builder("s1") + .bolus(-23.3, 100.0, 0) + .bolus(0.0, 100.0, 0) + .observation(1.0, 5.0, 0) + .build(); + let data = Data::from(subject); + + let expanded = data.expand(0.5, 1.0); + let times = sorted_obs_times(&expanded.subjects()[0].occasions()[0]); + + // First point at the first dose, then the grid continues from -23.0 + assert_eq!(times[..3], [-23.3, -23.0, -22.5]); + for t in times.iter().skip(1) { + let k = t / 0.5; + assert_eq!(k, k.round(), "{t} is not a multiple of idelta"); + } + assert!(times.contains(&0.0)); + assert_eq!(times.last(), Some(&1.0)); + } + + #[test] + fn expand_negative_observation_starts_grid_without_extra_point() { + // No doses: the grid starts at the first observation and ends at 0 + tad + let subject = Subject::builder("s1") + .observation(-0.5, 5.0, 0) + .observation(12.0, 3.0, 0) + .build(); + let data = Data::from(subject); + + let expanded = data.expand(1.0, 3.0); + let times = sorted_obs_times(&expanded.subjects()[0].occasions()[0]); + + assert_eq!(times, vec![-0.5, 0.0, 1.0, 2.0, 3.0, 12.0]); + } + + #[test] + fn expand_last_time_before_zero() { + // An observation is needed so the dataset has an output equation + let subject = Subject::builder("s1") + .bolus(-10.0, 100.0, 0) + .observation(-9.0, 5.0, 0) + .build(); + let data = Data::from(subject); + + let expanded = data.expand(1.0, 2.0); + let times = sorted_obs_times(&expanded.subjects()[0].occasions()[0]); + + assert_eq!(times, vec![-10.0, -9.0, -8.0]); + } } diff --git a/src/optimize/parameters.rs b/src/optimize/parameters.rs index fcf21a14..7d564549 100644 --- a/src/optimize/parameters.rs +++ b/src/optimize/parameters.rs @@ -74,7 +74,6 @@ impl<'a, E: Equation> ParameterOptimizer<'a, E> { } /// Optimize the parameters to minimize the negative log-likelihood against the data. - 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)?; From d66463a597bed8988d8e8070f60a4940b2681a0c Mon Sep 17 00:00:00 2001 From: Markus <66058642+mhovd@users.noreply.github.com> Date: Fri, 16 Oct 2026 12:10:57 +0200 Subject: [PATCH 2/2] fix: keep exact first-event time for expand's seed point The seed point at a negative first event used the rounded microsecond key, which could move it across the first dose or before the simulation start. Use the exact time and start the aligned grid strictly after it. Co-Authored-By: Claude Opus 5.5 --- src/data/structs.rs | 71 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 59 insertions(+), 12 deletions(-) diff --git a/src/data/structs.rs b/src/data/structs.rs index 5d7d03f0..94c321d9 100644 --- a/src/data/structs.rs +++ b/src/data/structs.rs @@ -231,12 +231,12 @@ impl Data { // Generate new observation times, stepping in integer // microseconds to guarantee forward progress. let mut new_events = Vec::new(); - let mut add_grid_point = |time_key: i64| { + let mut add_grid_point = |time: f64, time_key: i64| { for outeq in &outeq_values { // Only add if this (time, outeq) combination doesn't exist if !existing_obs.contains(&(time_key, outeq.clone())) { let obs = Observation::new( - time_key as f64 / 1e6, + time, None, outeq.clone(), None, @@ -249,23 +249,27 @@ impl Data { }; // Signed microseconds, so times before 0 neither saturate nor - // collide. When the first event is before 0, the grid starts at - // that event (so the simulation does not start earlier than the - // data) and continues from the next multiple of `step_us`, to - // stay aligned with the grid after 0. + // collide. When the first event is before 0, the grid starts with + // a point at that event's exact time, so it is never moved across + // the first dose or before the simulation start. It then continues + // from the first multiple of `step_us` strictly after that event, + // to stay aligned with the grid after 0. let last_time_us = (last_time * 1e6).round() as i64; - let first_time_us = (first_time * 1e6).round() as i64; - let mut time_key = if first_time_us < 0 { - let aligned = -(-first_time_us).div_euclid(step_us) * step_us; - if aligned != first_time_us && first_time_us <= last_time_us { - add_grid_point(first_time_us); + let mut time_key = if first_time < 0.0 { + let first_time_us = (first_time * 1e6).round() as i64; + if first_time <= last_time { + add_grid_point(first_time, first_time_us); + } + let mut aligned = -(-first_time_us).div_euclid(step_us) * step_us; + if aligned as f64 / 1e6 <= first_time { + aligned += step_us; } aligned } else { 0 }; while time_key <= last_time_us { - add_grid_point(time_key); + add_grid_point(time_key as f64 / 1e6, time_key); time_key += step_us; } @@ -1887,4 +1891,47 @@ mod tests { assert_eq!(times, vec![-10.0, -9.0, -8.0]); } + + /// Observation times after expanding a subject with doses at `first_dose` + /// and 0, an observation at 1, `idelta = 0.5` and `tad = 1`. + fn expanded_times_with_first_dose(first_dose: f64) -> Vec { + let subject = Subject::builder("s1") + .bolus(first_dose, 100.0, 0) + .bolus(0.0, 100.0, 0) + .observation(1.0, 5.0, 0) + .build(); + let expanded = Data::from(subject).expand(0.5, 1.0); + sorted_obs_times(&expanded.subjects()[0].occasions()[0]) + } + + #[test] + fn expand_first_point_keeps_exact_time_when_rounding_up() { + // -1.1000004 rounds to -1.1 µs-wise, which would put the first point + // after the first dose + let times = expanded_times_with_first_dose(-1.1000004); + assert_eq!(times[..3], [-1.1000004, -1.0, -0.5]); + } + + #[test] + fn expand_first_point_keeps_exact_time_when_rounding_down() { + // -1.1000006 rounds to -1.100001, which would start the simulation + // before the first dose + let times = expanded_times_with_first_dose(-1.1000006); + assert_eq!(times[..3], [-1.1000006, -1.0, -0.5]); + } + + #[test] + fn expand_first_point_rounding_onto_grid_multiple() { + // Rounds down onto -1.0: the next grid point must still be after the dose + let times = expanded_times_with_first_dose(-0.9999996); + assert_eq!(times[..3], [-0.9999996, -0.5, 0.0]); + + // Rounds up onto -1.0: -1.0 is after the dose and stays on the grid + let times = expanded_times_with_first_dose(-1.0000004); + assert_eq!(times[..3], [-1.0000004, -1.0, -0.5]); + + // Exactly on a grid multiple: one point there, no duplicate + let times = expanded_times_with_first_dose(-1.0); + assert_eq!(times[..3], [-1.0, -0.5, 0.0]); + } }