pmcore/bestdose/
predictions.rs

1//! AUC / dense-grid helpers used by the cost function.
2//!
3//! For the AUC targets, concentrations are simulated on a dense time grid and
4//! integrated with the trapezoidal rule:
5//!
6//! ```text
7//! AUC(t) = Σᵢ (C[i] + C[i-1]) / 2 × (t[i] - t[i-1])
8//! ```
9
10use pharmsol::prelude::*;
11
12/// Find the time of the last dose (bolus or infusion) before a given observation
13/// time. Returns `0.0` if no dose exists before `obs_time`.
14pub fn find_last_dose_time_before(subject: &Subject, obs_time: f64) -> f64 {
15    let mut last_dose_time = 0.0;
16
17    for occasion in subject.occasions() {
18        for event in occasion.events() {
19            let event_time = match event {
20                Event::Bolus(b) => Some(b.time()),
21                Event::Infusion(i) => Some(i.time()),
22                Event::Observation(_) => None,
23            };
24
25            if let Some(t) = event_time {
26                if t < obs_time && t > last_dose_time {
27                    last_dose_time = t;
28                }
29            }
30        }
31    }
32
33    last_dose_time
34}
35
36/// Generate a dense time grid for AUC calculations.
37///
38/// The grid contains the observation times plus regular points spaced `idelta`
39/// apart (in the model's time units), sorted and deduplicated. A non-positive
40/// `idelta` disables the regular grid, leaving only the observation times.
41pub fn calculate_dense_times(
42    start_time: f64,
43    end_time: f64,
44    obs_times: &[f64],
45    idelta: f64,
46) -> Vec<f64> {
47    let mut times = Vec::new();
48
49    times.extend_from_slice(obs_times);
50
51    if idelta > 0.0 {
52        let mut t = start_time;
53        while t <= end_time {
54            times.push(t);
55            t += idelta;
56        }
57    }
58
59    if !times.contains(&end_time) {
60        times.push(end_time);
61    }
62
63    times.sort_by(|a, b| a.partial_cmp(b).unwrap());
64
65    let tolerance = 1e-10;
66    let mut unique_times = Vec::new();
67    let mut last_time = f64::NEG_INFINITY;
68
69    for &t in &times {
70        if (t - last_time).abs() > tolerance {
71            unique_times.push(t);
72            last_time = t;
73        }
74    }
75
76    unique_times
77}
78
79/// Calculate cumulative AUC at target times using the trapezoidal rule.
80///
81/// Integrates from the first dense time point and extracts the cumulative AUC at
82/// each of `target_times`.
83pub fn calculate_auc_at_times(
84    dense_times: &[f64],
85    dense_predictions: &[f64],
86    target_times: &[f64],
87) -> Vec<f64> {
88    assert_eq!(dense_times.len(), dense_predictions.len());
89
90    let mut target_aucs = Vec::with_capacity(target_times.len());
91    let mut auc = 0.0;
92    let mut target_idx = 0;
93    let tolerance = 1e-10;
94
95    for i in 1..dense_times.len() {
96        let dt = dense_times[i] - dense_times[i - 1];
97        let avg_conc = (dense_predictions[i] + dense_predictions[i - 1]) / 2.0;
98        auc += avg_conc * dt;
99
100        if target_idx < target_times.len()
101            && (dense_times[i] - target_times[target_idx]).abs() < tolerance
102        {
103            target_aucs.push(auc);
104            target_idx += 1;
105        }
106    }
107
108    target_aucs
109}
110
111/// Calculate interval AUC for each observation independently.
112///
113/// For each observation at time `t`, integrates from the last dose before `t` to
114/// `t` (e.g. dosing-interval AUCτ at steady state).
115pub fn calculate_interval_auc_per_observation(
116    subject: &Subject,
117    dense_times: &[f64],
118    dense_predictions: &[f64],
119    obs_times: &[f64],
120) -> Vec<f64> {
121    assert_eq!(dense_times.len(), dense_predictions.len());
122
123    let mut interval_aucs = Vec::with_capacity(obs_times.len());
124    let tolerance = 1e-10;
125
126    for &obs_time in obs_times {
127        let last_dose_time = find_last_dose_time_before(subject, obs_time);
128
129        let start_idx = dense_times
130            .iter()
131            .position(|&t| (t - last_dose_time).abs() < tolerance || t > last_dose_time)
132            .unwrap_or(0);
133
134        let end_idx = dense_times
135            .iter()
136            .position(|&t| (t - obs_time).abs() < tolerance || t > obs_time)
137            .unwrap_or(dense_times.len() - 1);
138
139        let mut auc = 0.0;
140        for i in (start_idx + 1)..=end_idx.min(dense_times.len() - 1) {
141            let dt = dense_times[i] - dense_times[i - 1];
142            let avg_conc = (dense_predictions[i] + dense_predictions[i - 1]) / 2.0;
143            auc += avg_conc * dt;
144        }
145
146        interval_aucs.push(auc);
147    }
148
149    interval_aucs
150}
pmcore/bestdose/predictions.rs

pmcore/bestdose/
predictions.rs
