pmcore/estimation/nonparametric/
predictions.rs

1use std::path::Path;
2
3use anyhow::{bail, Result};
4use pharmsol::{prelude::simulator::Prediction, Censor, Data, Predictions as PredTrait};
5use rayon::prelude::*;
6use serde::{Deserialize, Serialize};
7
8use crate::{
9    estimation::nonparametric::{theta::Theta, weights::Weights},
10    estimation::nonparametric::{weighted_median, Posterior},
11};
12
13#[derive(Debug, Clone, Serialize, Deserialize)]
14pub struct NPPredictionRow {
15    id: String,
16    time: f64,
17    outeq: usize,
18    block: usize,
19    obs: Option<f64>,
20    cens: Censor,
21    pop_mean: f64,
22    pop_median: f64,
23    post_mean: f64,
24    post_median: f64,
25}
26
27impl NPPredictionRow {
28    pub fn id(&self) -> &str {
29        &self.id
30    }
31    pub fn time(&self) -> f64 {
32        self.time
33    }
34    pub fn outeq(&self) -> usize {
35        self.outeq
36    }
37    pub fn block(&self) -> usize {
38        self.block
39    }
40    pub fn obs(&self) -> Option<f64> {
41        self.obs
42    }
43    pub fn pop_mean(&self) -> f64 {
44        self.pop_mean
45    }
46    pub fn pop_median(&self) -> f64 {
47        self.pop_median
48    }
49    pub fn post_mean(&self) -> f64 {
50        self.post_mean
51    }
52    pub fn post_median(&self) -> f64 {
53        self.post_median
54    }
55
56    pub fn censoring(&self) -> Censor {
57        self.cens
58    }
59}
60
61#[derive(Debug, Clone, Serialize, Deserialize)]
62pub struct NPPredictions {
63    predictions: Vec<NPPredictionRow>,
64}
65
66impl IntoIterator for NPPredictions {
67    type Item = NPPredictionRow;
68    type IntoIter = std::vec::IntoIter<NPPredictionRow>;
69
70    fn into_iter(self) -> Self::IntoIter {
71        self.predictions.into_iter()
72    }
73}
74
75impl Default for NPPredictions {
76    fn default() -> Self {
77        NPPredictions::new()
78    }
79}
80
81impl NPPredictions {
82    pub fn new() -> Self {
83        NPPredictions {
84            predictions: Vec::new(),
85        }
86    }
87
88    pub fn add(&mut self, row: NPPredictionRow) {
89        self.predictions.push(row);
90    }
91
92    pub fn predictions(&self) -> &[NPPredictionRow] {
93        &self.predictions
94    }
95
96    /// Write the predictions to a CSV file readable by Pmetrics.
97    ///
98    /// The parent directory is created if it does not already exist. The header
99    /// matches the fields of [`NPPredictionRow`]:
100    /// `id,time,outeq,block,obs,cens,pop_mean,pop_median,post_mean,post_median`.
101    pub fn write(&self, path: &Path) -> Result<()> {
102        tracing::debug!("Writing predictions...");
103
104        super::create_parent_dir(path)?;
105
106        let mut writer = csv::WriterBuilder::new()
107            .has_headers(true)
108            .from_path(path)?;
109
110        for row in &self.predictions {
111            writer.serialize(row)?;
112        }
113
114        writer.flush()?;
115        Ok(())
116    }
117
118    pub fn calculate(
119        equation: &(impl pharmsol::prelude::simulator::Equation + Sync),
120        data: &Data,
121        theta: &Theta,
122        w: &Weights,
123        posterior: &Posterior,
124        idelta: f64,
125        tad: f64,
126    ) -> Result<Self> {
127        let data = data.clone().expand(idelta, tad);
128        let subjects = data.subjects();
129
130        if subjects.len() != posterior.matrix().nrows() {
131            bail!("Number of subjects and number of posterior means do not match");
132        };
133
134        let support_points: Vec<Vec<f64>> = theta
135            .matrix()
136            .row_iter()
137            .map(|spp| spp.iter().cloned().collect())
138            .collect();
139
140        let per_subject: Vec<Vec<NPPredictionRow>> = subjects
141            .par_iter()
142            .enumerate()
143            .map(|(subject_index, subject)| -> Result<Vec<NPPredictionRow>> {
144                let predictions: Vec<Vec<Prediction>> = support_points
145                    .par_iter()
146                    .map(|spp| {
147                        Ok(equation
148                            .simulate_subject_dense(subject, spp, None)?
149                            .0
150                            .get_predictions())
151                    })
152                    .collect::<Result<Vec<_>>>()?;
153
154                let Some(first_spp_preds) = predictions.first() else {
155                    return Ok(Vec::new());
156                };
157                let n_points = first_spp_preds.len();
158
159                let mut pop_mean: Vec<f64> = vec![0.0; n_points];
160                for (i, outer_pred) in predictions.iter().enumerate() {
161                    for (j, pred) in outer_pred.iter().enumerate() {
162                        pop_mean[j] += pred.prediction() * w[i];
163                    }
164                }
165
166                let mut pop_median: Vec<f64> = Vec::with_capacity(n_points);
167                for j in 0..n_points {
168                    let mut values: Vec<f64> = Vec::new();
169                    let mut weights: Vec<f64> = Vec::new();
170
171                    for (i, outer_pred) in predictions.iter().enumerate() {
172                        values.push(outer_pred[j].prediction());
173                        weights.push(w[i]);
174                    }
175
176                    pop_median.push(weighted_median(&values, &weights));
177                }
178
179                let mut posterior_mean: Vec<f64> = vec![0.0; n_points];
180                for (i, outer_pred) in predictions.iter().enumerate() {
181                    for (j, pred) in outer_pred.iter().enumerate() {
182                        posterior_mean[j] +=
183                            pred.prediction() * posterior.matrix()[(subject_index, i)];
184                    }
185                }
186
187                let mut posterior_median: Vec<f64> = Vec::with_capacity(n_points);
188                for j in 0..n_points {
189                    let mut values: Vec<f64> = Vec::new();
190                    let mut weights: Vec<f64> = Vec::new();
191
192                    for (i, outer_pred) in predictions.iter().enumerate() {
193                        values.push(outer_pred[j].prediction());
194                        weights.push(posterior.matrix()[(subject_index, i)]);
195                    }
196
197                    posterior_median.push(weighted_median(&values, &weights));
198                }
199
200                Ok(first_spp_preds
201                    .iter()
202                    .enumerate()
203                    .map(|(j, p)| NPPredictionRow {
204                        id: subject.id().clone(),
205                        time: p.time(),
206                        outeq: p.outeq(),
207                        block: p.occasion(),
208                        obs: p.observation(),
209                        cens: p.censoring(),
210                        pop_mean: pop_mean[j],
211                        pop_median: pop_median[j],
212                        post_mean: posterior_mean[j],
213                        post_median: posterior_median[j],
214                    })
215                    .collect())
216            })
217            .collect::<Result<Vec<_>>>()?;
218
219        let mut container = NPPredictions::new();
220        for rows in per_subject {
221            for row in rows {
222                container.add(row);
223            }
224        }
225
226        Ok(container)
227    }
228}
pmcore/estimation/nonparametric/predictions.rs

pmcore/estimation/nonparametric/
predictions.rs
