diff --git a/.gitignore b/.gitignore index 98e91c40..1ee158ab 100644 --- a/.gitignore +++ b/.gitignore @@ -27,4 +27,5 @@ settings.json /benches/*.json log.txt op.csv -*results.txt \ No newline at end of file +*results.txt +covs.csv diff --git a/Cargo.toml b/Cargo.toml index e520f20b..ceb8fdf7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -27,15 +27,16 @@ tokio = { version = "1.32.0", features = ["sync", "rt"] } argmin = { version = "0.10.0", features = [] } argmin-math = { version = "0.4.0", features = ["ndarray_v0_15-nolinalg"] } tracing = "0.1.40" -tracing-subscriber = { version = "0.3.17", features = ["env-filter", "fmt", "time"] } +tracing-subscriber = { version = "0.3.17", features = [ + "env-filter", + "fmt", + "time", +] } chrono = "0.4" -config = {version = "0.14", features=["preserve_order"]} +config = { version = "0.14", features = ["preserve_order"] } faer = "0.19.3" -faer-ext = {version = "0.2.0", features = ["nalgebra", "ndarray"]} -pharmsol = "0.7.0" -# pharmsol = {path = "../pharmsol"} -# pharmsol = {git ="https://github.com/LAPKB/pharmsol", branch="net"} -# pharmsol = {git ="https://github.com/LAPKB/pharmsol", branch="main"} +faer-ext = { version = "0.2.0", features = ["nalgebra", "ndarray"] } +pharmsol = "0.7.1" toml = "0.8.14" #REMOVE rand = "0.8.5" anyhow = "1.0.86" diff --git a/src/routines/output.rs b/src/routines/output.rs index 8eb75d0a..353f6494 100644 --- a/src/routines/output.rs +++ b/src/routines/output.rs @@ -5,6 +5,7 @@ use ndarray::parallel::prelude::*; use ndarray::{Array, Array1, Array2, Axis}; use pharmsol::prelude::data::*; use pharmsol::prelude::simulator::Equation; +use serde::Serialize; // use pharmsol::Cache; use settings::Settings; use std::fs::{create_dir_all, File, OpenOptions}; @@ -82,6 +83,7 @@ impl NPResult { .context("Failed to write observed-predicted file")?; self.write_pred(idelta, tad) .context("Failed to write predictions")?; + self.write_covs().context("Failed to write covariates")?; if self.w.len() > 0 { //TODO: find a better way to indicate that the run failed self.write_posterior() @@ -95,6 +97,19 @@ impl NPResult { pub fn write_obspred(&self) -> Result<()> { tracing::debug!("Writing observations and predictions..."); + #[derive(Debug, Clone, Serialize)] + struct Row { + id: String, + time: f64, + outeq: usize, + block: usize, + obs: f64, + pop_mean: f64, + pop_median: f64, + post_mean: f64, + post_median: f64, + } + let theta: Array2 = self.theta.clone(); let w: Array1 = self.w.clone(); let psi: Array2 = self.psi.clone(); @@ -119,69 +134,56 @@ impl NPResult { .has_headers(true) .from_writer(&outputfile.file); - // Create the headers - writer.write_record([ - "id", - "time", - "outeq", - "obs", - "popMean", - "popMedian", - "postMean", - "postMedian", - ])?; - for (i, subject) in subjects.iter().enumerate() { - // Population predictions - let pop_mean_pred = self - .equation - .simulate_subject(subject, &pop_mean.to_vec(), None) - .0 - .get_predictions() - .clone(); - let pop_median_pred = self - .equation - .simulate_subject(subject, &pop_median.to_vec(), None) - .0 - .get_predictions() - .clone(); - - // Posterior predictions - let post_mean_spp: Vec = post_mean.row(i).to_vec(); - let post_mean_pred = self - .equation - .simulate_subject(subject, &post_mean_spp, None) - .0 - .get_predictions() - .clone(); - let post_median_spp: Vec = post_median.row(i).to_vec(); - let post_median_pred = self - .equation - .simulate_subject(subject, &post_median_spp, None) - .0 - .get_predictions() - .clone(); - - // Write the records - pop_mean_pred - .iter() - .zip(pop_median_pred.iter()) - .zip(post_mean_pred.iter()) - .zip(post_median_pred.iter()) - .for_each(|(((pop_mean, pop_median), post_mean), post_median)| { - writer - .write_record([ - subject.id(), - &format!("{:.3}", pop_mean.time()), - &format!("{}", pop_mean.outeq()), - &format!("{:.4}", pop_mean.observation()), - &format!("{:.4}", pop_mean.prediction()), - &format!("{:.4}", pop_median.prediction()), - &format!("{:.4}", post_mean.prediction()), - &format!("{:.4}", post_median.prediction()), - ]) - .expect("Failed to write record"); - }); + for occasion in subject.occasions() { + let id = subject.id(); + let occ = occasion.index(); + + let subject = Subject::from_occasions(id.clone(), vec![occasion.clone()]); + + // Population predictions + let pop_mean_pred = self + .equation + .simulate_subject(&subject, &pop_mean.to_vec(), None) + .0 + .get_predictions() + .clone(); + let pop_median_pred = self + .equation + .simulate_subject(&subject, &pop_median.to_vec(), None) + .0 + .get_predictions() + .clone(); + + // Posterior predictions + let post_mean_spp: Vec = post_mean.row(i).to_vec(); + let post_mean_pred = self + .equation + .simulate_subject(&subject, &post_mean_spp, None) + .0 + .get_predictions() + .clone(); + let post_median_spp: Vec = post_median.row(i).to_vec(); + let post_median_pred = self + .equation + .simulate_subject(&subject, &post_median_spp, None) + .0 + .get_predictions() + .clone(); + + let row = Row { + id: id.clone(), + time: pop_mean_pred[0].time(), + outeq: pop_mean_pred[0].outeq(), + block: occ, + obs: pop_mean_pred[0].observation(), + pop_mean: pop_mean_pred[0].prediction(), + pop_median: pop_median_pred[0].prediction(), + post_mean: post_mean_pred[0].prediction(), + post_median: post_median_pred[0].prediction(), + }; + writer.serialize(row)?; + } } writer.flush()?; tracing::info!( @@ -190,6 +192,7 @@ impl NPResult { ); Ok(()) } + /// Writes theta, which contains the population support points and their associated probabilities /// Each row is one support point, the last column being probability pub fn write_theta(&self) -> Result<()> { @@ -315,7 +318,7 @@ impl NPResult { .context("Failed to calculate posterior mean and median")?; let (pop_mean, pop_median) = population_mean_median(&theta, &w) - .context("Failed to calculate posterior mean and median")?; + .context("Failed to calculate population mean and median")?; let subjects = data.get_subjects(); if subjects.len() != post_mean.nrows() { @@ -327,66 +330,76 @@ impl NPResult { .has_headers(true) .from_writer(&outputfile.file); - // Create the headers - writer.write_record([ - "id", - "time", - "outeq", - "popMean", - "popMedian", - "postMean", - "postMedian", - ])?; + #[derive(Debug, Clone, Serialize)] + struct Row { + id: String, + time: f64, + outeq: usize, + block: usize, + pop_mean: f64, + pop_median: f64, + post_mean: f64, + post_median: f64, + } for (i, subject) in subjects.iter().enumerate() { - // Population predictions - let pop_mean_pred = self - .equation - .simulate_subject(subject, &pop_mean.to_vec(), None) - .0 - .get_predictions() - .clone(); - let pop_median_pred = self - .equation - .simulate_subject(subject, &pop_median.to_vec(), None) - .0 - .get_predictions() - .clone(); - - // Posterior predictions - let post_mean_spp: Vec = post_mean.row(i).to_vec(); - let post_mean_pred = self - .equation - .simulate_subject(subject, &post_mean_spp, None) - .0 - .get_predictions() - .clone(); - let post_median_spp: Vec = post_median.row(i).to_vec(); - let post_median_pred = self - .equation - .simulate_subject(subject, &post_median_spp, None) - .0 - .get_predictions() - .clone(); - - pop_mean_pred - .iter() - .zip(pop_median_pred.iter()) - .zip(post_mean_pred.iter()) - .zip(post_median_pred.iter()) - .for_each(|(((pop_mean, pop_median), post_mean), post_median)| { - writer - .write_record([ - subject.id(), - &format!("{:.3}", pop_mean.time()), - &format!("{}", pop_mean.outeq()), - &format!("{:.4}", pop_mean.prediction()), - &format!("{:.4}", pop_median.prediction()), - &format!("{:.4}", post_mean.prediction()), - &format!("{:.4}", post_median.prediction()), - ]) - .expect("Failed to write record"); - }); + for occasion in subject.occasions() { + let id = subject.id(); + let block = occasion.index(); + + // Create a new subject with only the current occasion + let subject = Subject::from_occasions(id.clone(), vec![occasion.clone()]); + + // Population predictions + let pop_mean_pred = self + .equation + .simulate_subject(&subject, &pop_mean.to_vec(), None) + .0 + .get_predictions() + .clone(); + let pop_median_pred = self + .equation + .simulate_subject(&subject, &pop_median.to_vec(), None) + .0 + .get_predictions() + .clone(); + + // Posterior predictions + let post_mean_spp: Vec = post_mean.row(i).to_vec(); + let post_mean_pred = self + .equation + .simulate_subject(&subject, &post_mean_spp, None) + .0 + .get_predictions() + .clone(); + let post_median_spp: Vec = post_median.row(i).to_vec(); + let post_median_pred = self + .equation + .simulate_subject(&subject, &post_median_spp, None) + .0 + .get_predictions() + .clone(); + + // Write predictions for each time point + for (((pop_mean, pop_median), post_mean), post_median) in pop_mean_pred + .iter() + .zip(pop_median_pred.iter()) + .zip(post_mean_pred.iter()) + .zip(post_median_pred.iter()) + { + let row = Row { + id: id.clone(), + time: pop_mean.time(), + outeq: pop_mean.outeq(), + block, + pop_mean: pop_mean.prediction(), + pop_median: pop_median.prediction(), + post_mean: post_mean.prediction(), + post_median: post_median.prediction(), + }; + writer.serialize(row)?; + } + } } writer.flush()?; tracing::info!( @@ -395,6 +408,79 @@ impl NPResult { ); Ok(()) } + + /// Writes the covariates + pub fn write_covs(&self) -> Result<()> { + tracing::debug!("Writing covariates..."); + let outputfile = OutputFile::new(&self.settings.output.path, "covs.csv")?; + let mut writer = WriterBuilder::new() + .has_headers(true) + .from_writer(&outputfile.file); + + // Collect all unique covariate names + let mut covariate_names = std::collections::HashSet::new(); + for subject in self.data.get_subjects() { + for occasion in subject.occasions() { + if let Some(cov) = occasion.get_covariates() { + let covmap = cov.covariates(); + for (cov_name, _) in &covmap { + covariate_names.insert(cov_name.clone()); + } + } + } + } + let mut covariate_names: Vec = covariate_names.into_iter().collect(); + covariate_names.sort(); // Ensure consistent order + + // Write the header row: id, time, block, covariate names + let mut headers = vec!["id", "time", "block"]; + headers.extend(covariate_names.iter().map(|s| s.as_str())); + writer.write_record(&headers)?; + + // Write the data rows + for subject in self.data.get_subjects() { + for occasion in subject.occasions() { + if let Some(cov) = occasion.get_covariates() { + let covmap = cov.covariates(); + + for event in occasion.get_events(&None, &None, false) { + let time = match event { + Event::Bolus(bolus) => bolus.time(), + Event::Infusion(infusion) => infusion.time(), + Event::Observation(observation) => observation.time(), + }; + + let mut row: Vec = Vec::new(); + row.push(subject.id().clone()); + row.push(time.to_string()); + row.push(occasion.index().to_string()); + + // Add covariate values to the row + for cov_name in &covariate_names { + if let Some(cov) = covmap.get(cov_name) { + if let Some(value) = cov.interpolate(time) { + row.push(value.to_string()); + } else { + row.push(String::new()); + } + } else { + row.push(String::new()); + } + } + + writer.write_record(&row)?; + } + } + } + } + + writer.flush()?; + tracing::info!( + "Covariates written to {:?}", + &outputfile.get_relative_path() + ); + Ok(()) + } } /// An [NPCycle] object contains the summary of a cycle