From a586a2330f3b7b29972513964390600ba4eb35b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Wed, 5 Nov 2025 15:48:36 +0000 Subject: [PATCH 1/5] initial version --- NAMESPACE | 2 + R/PM_bestdose.R | 359 ++++++++++++++++++++++++++++++ R/PM_result.R | 66 ++++++ R/extendr-wrappers.R | 4 + man/PM_model.Rd | 86 +++---- man/PM_result.Rd | 105 +++++++++ man/interp.Rd | 2 +- src/rust/src/bestdose_executor.rs | 214 ++++++++++++++++++ src/rust/src/lib.rs | 70 ++++++ 9 files changed, 864 insertions(+), 44 deletions(-) create mode 100644 R/PM_bestdose.R create mode 100644 src/rust/src/bestdose_executor.rs diff --git a/NAMESPACE b/NAMESPACE index 07e9b8d43..99e1132dd 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -47,6 +47,7 @@ export(NPreport) export(NPrun) export(PMFortranConfig) export(PM_batch) +export(PM_bestdose) export(PM_build) export(PM_compare) export(PM_cov) @@ -94,6 +95,7 @@ export(add_renal) export(add_shapes) export(add_smooth) export(additive) +export(bestdose) export(build_model) export(build_plot) export(cli_ask) diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R new file mode 100644 index 000000000..5040fccba --- /dev/null +++ b/R/PM_bestdose.R @@ -0,0 +1,359 @@ +#' @title +#' Object to contain BestDose optimization results +#' +#' @description +#' `r lifecycle::badge("experimental")` +#' +#' This object is created after a successful BestDose optimization run. +#' BestDose finds optimal dosing regimens to achieve target drug concentrations +#' or AUC values using Bayesian optimization. +#' +#' @export +PM_bestdose <- R6::R6Class( + "PM_bestdose", + public = list( + #' @field result A list containing: + #' * **doses** Vector of optimal dose amounts (mg) + #' * **objf** Objective function value (lower is better) + #' * **status** Optimization status ("converged" or "max_iterations") + #' * **predictions** Data frame with concentration-time predictions + #' * **auc_predictions** Data frame with AUC predictions (if applicable) + #' * **method** Optimization method used ("posterior" or "uniform") + result = NULL, + + #' @description + #' Create a new PM_bestdose object by running BestDose optimization + #' + #' @param prior A PM_result or PM_final object, or path to theta.csv file + #' @param model A PM_model object or path to compiled model + #' @param past_data Optional: PM_data object or path to CSV with patient history + #' @param target PM_data object or path to CSV with target doses/observations + #' @param dose_range Named list with min and max allowable doses, e.g., + #' list(min = 0, max = 1000). Default: list(min = 0, max = 1000) + #' @param bias_weight Numeric [0,1]. 0 = fully personalized, 1 = population-based. + #' Default = 0.5 (balanced) + #' @param target_type One of "concentration", "auc_from_zero", or "auc_from_last_dose". + #' Default = "concentration" + #' @param time_offset Optional: time offset for past/future concatenation + #' @param settings Optional: list of settings (algorithm, error models, etc.) + #' If NULL, will use defaults from prior object + #' + #' @details + #' + #' ## Target Data Format + #' + #' The target data should be a Pmetrics-formatted CSV with: + #' - Dose events: Set dose amount to 0 for doses to be optimized, + #' or non-zero for fixed doses + #' - Observation events: Set OUT values to the target concentrations or AUCs + #' + #' ## Target Types + #' + #' - **concentration**: Optimize to match target concentrations at observation times + #' - **auc_from_zero**: Optimize to match cumulative AUC from time 0 + #' - **auc_from_last_dose**: Optimize to match interval AUC from last dose + #' + #' ## Bias Weight + #' + #' Controls the balance between patient-specific and population-based optimization: + #' - 0.0: Fully personalized (minimize variance) + #' - 0.5: Balanced approach (recommended default) + #' - 1.0: Population-based (minimize bias from population) + #' + #' @examples + #' \dontrun{ + #' # Load NPAG result + #' result <- PM_load(1) + #' + #' # Create target data + #' target <- PM_data$new("target.csv") + #' + #' # Run BestDose optimization + #' bd <- PM_bestdose$new( + #' prior = result, + #' model = result$model, + #' target = target, + #' dose_range = list(min = 50, max = 500), + #' bias_weight = 0.5 + #' ) + #' + #' # View results + #' print(bd) + #' bd$get_doses() + #' bd$get_predictions() + #' } + initialize = function(prior, + model, + past_data = NULL, + target, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + settings = NULL) { + + # Validate inputs + if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) { + cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose") + } + + if (bias_weight < 0 || bias_weight > 1) { + cli::cli_abort("bias_weight must be between 0 and 1") + } + + if (is.null(dose_range$min) || is.null(dose_range$max)) { + cli::cli_abort("dose_range must have both 'min' and 'max' elements") + } + + if (dose_range$min >= dose_range$max) { + cli::cli_abort("dose_range$min must be less than dose_range$max") + } + + # Parse prior + prior_path <- private$.parse_prior(prior) + + # Parse model + model_info <- private$.parse_model(model) + + # Parse data files + past_data_path <- if (!is.null(past_data)) { + private$.parse_data(past_data) + } else { + NULL + } + + target_data_path <- private$.parse_data(target) + + # Parse settings + if (is.null(settings)) { + settings <- private$.default_bestdose_settings(prior) + } + + # Call Rust function + cli::cli_alert_info("Running BestDose optimization...") + + result <- tryCatch({ + bestdose( + model_path = model_info$path, + prior_path = prior_path, + past_data_path = past_data_path, + target_data_path = target_data_path, + time_offset = time_offset, + dose_min = dose_range$min, + dose_max = dose_range$max, + bias_weight = bias_weight, + target_type = target_type, + params = settings, + kind = model_info$kind + ) + }, error = function(e) { + cli::cli_abort(c( + "x" = "BestDose optimization failed", + "i" = conditionMessage(e) + )) + }) + + self$result <- result + + cli::cli_alert_success("Optimization complete!") + cli::cli_alert_info("Method: {result$method}") + cli::cli_alert_info("Status: {result$status}") + cli::cli_alert_info("Optimal doses: {paste(round(result$doses, 2), collapse = ', ')} mg") + + invisible(self) + }, + + #' @description + #' Print summary of BestDose results + print = function() { + cat("BestDose Optimization Results\n") + cat("==============================\n\n") + cat(sprintf("Optimal doses: %s mg\n", paste(round(self$result$doses, 2), collapse = ", "))) + cat(sprintf("Objective function: %.6f\n", self$result$objf)) + cat(sprintf("Status: %s\n", self$result$status)) + cat(sprintf("Method: %s\n", self$result$method)) + cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions))) + if (!is.null(self$result$auc_predictions)) { + cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions))) + } + invisible(self) + }, + + #' @description + #' Get optimal dose values + #' @return Numeric vector of optimal doses + get_doses = function() { + self$result$doses + }, + + #' @description + #' Get concentration-time predictions + #' @return Data frame with predictions + get_predictions = function() { + self$result$predictions + }, + + #' @description + #' Get AUC predictions (if available) + #' @return Data frame with AUC predictions or NULL + get_auc_predictions = function() { + self$result$auc_predictions + }, + + #' @description + #' Get objective function value + #' @return Numeric objective function value + get_objf = function() { + self$result$objf + }, + + #' @description + #' Get optimization status + #' @return Character string with status + get_status = function() { + self$result$status + }, + + #' @description + #' Get optimization method used + #' @return Character string: "posterior" or "uniform" + get_method = function() { + self$result$method + }, + + #' @description + #' Save results to RDS file + #' @param filename Path to save file. Default: "bestdose_result.rds" + save = function(filename = "bestdose_result.rds") { + saveRDS(self, filename) + cli::cli_alert_success("Results saved to {filename}") + invisible(self) + } + ), + + private = list( + # Helper function to parse prior + .parse_prior = function(prior) { + if (inherits(prior, "PM_result")) { + # Extract theta.csv from PM_result + theta_path <- file.path(prior$rundir, "outputs", "theta.csv") + if (!file.exists(theta_path)) { + cli::cli_abort("theta.csv not found in PM_result outputs") + } + return(theta_path) + } else if (inherits(prior, "PM_final")) { + # Need to write out PM_final to CSV + temp_path <- tempfile(fileext = ".csv") + private$.write_prior_csv(prior, temp_path) + return(temp_path) + } else if (is.character(prior)) { + if (!file.exists(prior)) { + cli::cli_abort("Prior file not found: {prior}") + } + return(prior) + } else { + cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv") + } + }, + + # Helper function to parse model + .parse_model = function(model) { + if (inherits(model, "PM_model")) { + # Check if model is compiled + compiled_path <- model$model_list$compiled_location + if (is.null(compiled_path) || !file.exists(compiled_path)) { + cli::cli_abort("Model must be compiled first. Use model$compile()") + } + + # Determine model type from model_list + kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) { + "analytical" + } else { + "ode" + } + + return(list( + path = compiled_path, + kind = kind + )) + } else if (is.character(model)) { + # Assume it's a path to compiled model + if (!file.exists(model)) { + cli::cli_abort("Model file not found: {model}") + } + # Try to infer type from extension or filename + kind <- if (grepl("analytical", model, ignore.case = TRUE)) { + "analytical" + } else { + "ode" + } + return(list(path = model, kind = kind)) + } else { + cli::cli_abort("model must be PM_model or path to compiled model") + } + }, + + # Helper function to parse data + .parse_data = function(data) { + if (inherits(data, "PM_data")) { + # Write to temp CSV + temp_path <- tempfile(fileext = ".csv") + write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE) + return(temp_path) + } else if (is.character(data)) { + if (!file.exists(data)) { + cli::cli_abort("Data file not found: {data}") + } + return(data) + } else { + cli::cli_abort("data must be PM_data or path to CSV file") + } + }, + + # Helper to create default settings + .default_bestdose_settings = function(prior) { + # Extract settings from prior if it's a PM_result + if (inherits(prior, "PM_result")) { + # Use the settings from the result + settings <- prior$settings + return(settings) + } else { + # Use sensible defaults + settings <- list( + algorithm = "NPAG", + ranges = list(), # Will be filled from prior + error_models = list( + list( + initial = 0.0, + type = "additive", + coeff = c(0.0, 0.2, 0.0, 0.0) + ) + ), + max_cycles = 500, + points = 2028, + seed = 22, + prior = "prior.csv", + idelta = 0.25, # 15 minutes for AUC calculations + tad = 0.0 + ) + } + return(settings) + }, + + # Helper to write PM_final to CSV + .write_prior_csv = function(prior, path) { + # Convert PM_final to CSV with prob column + df <- as.data.frame(prior$popPoints) + df$prob <- prior$popProb + write.csv(df, path, row.names = FALSE, quote = FALSE) + } + ) +) + +#' @export +PM_bestdose$load <- function(filename = "bestdose_result.rds") { + if (!file.exists(filename)) { + cli::cli_abort("File not found: {filename}") + } + readRDS(filename) +} diff --git a/R/PM_result.R b/R/PM_result.R index e0e324523..93474777a 100755 --- a/R/PM_result.R +++ b/R/PM_result.R @@ -196,6 +196,72 @@ PM_result <- R6::R6Class( return(sim) }, + #' @description + #' Run BestDose optimization using this result as the prior + #' @details + #' BestDose finds optimal dosing regimens to achieve target drug concentrations + #' or AUC values. The algorithm uses Bayesian posterior estimation combined with + #' dual optimization to balance patient-specific adaptation and population-level + #' robustness. By default, uses the `$final`, `$model`, and `$data` objects from + #' this result. Most commonly, you will supply a different `target` data object. + #' @param target PM_data object or path to CSV with target doses/observations. + #' Required. This defines the dosing template and target values. Set dose amounts + #' to 0 for doses to be optimized. + #' @param past_data Optional: PM_data object or path to CSV with patient history. + #' If NULL (default), uses the full dataset from this result. + #' @param dose_range Named list with min and max allowable doses. + #' Default: list(min = 0, max = 1000) + #' @param bias_weight Numeric [0,1] controlling personalization level. + #' 0 = fully personalized, 1 = population-based. Default: 0.5 (balanced) + #' @param target_type One of "concentration" (default), "auc_from_zero", or + #' "auc_from_last_dose" + #' @param time_offset Optional: time offset for past/future concatenation + #' @param ... Additional arguments passed to [PM_bestdose] + #' @return A [PM_bestdose] object containing optimal doses and predictions + #' @examples + #' \dontrun{ + #' # Load NPAG result + #' result <- PM_load(1) + #' + #' # Create target data + #' target <- PM_data$new("target.csv") + #' + #' # Run BestDose optimization + #' bd <- result$bestdose( + #' target = target, + #' dose_range = list(min = 50, max = 500), + #' bias_weight = 0.5 + #' ) + #' + #' # View results + #' print(bd) + #' bd$get_doses() + #' } + bestdose = function(target, + past_data = NULL, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + ...) { + # Use this result's data as past_data if not specified + if (is.null(past_data)) { + past_data <- self$data + } + + PM_bestdose$new( + prior = self, + model = self$model, + past_data = past_data, + target = target, + dose_range = dose_range, + bias_weight = bias_weight, + target_type = target_type, + time_offset = time_offset, + ... + ) + }, + #' @description #' Save the current PM_result object to an .Rdata file. #' @details diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R index 3694efead..c65e13970 100755 --- a/R/extendr-wrappers.R +++ b/R/extendr-wrappers.R @@ -40,5 +40,9 @@ model_parameters <- function(model_path, kind) .Call(wrap__model_parameters, mod #'@export setup_logs <- function() .Call(wrap__setup_logs) +#' Run BestDose optimization to find optimal doses +#'@export +bestdose <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) + # nolint end diff --git a/man/PM_model.Rd b/man/PM_model.Rd index 4eb85ea03..685ee7ba1 100755 --- a/man/PM_model.Rd +++ b/man/PM_model.Rd @@ -32,41 +32,41 @@ a \code{donttest} block to avoid automatic compilation. mod_list <- list( pri = c( - CL = ab(10, 200), - V0 = ab(0, 100), - ka = ab(0, 3), - k23 = ab(0, 5), - k32 = ab(0, 5), - lag1 = ab(0, 2) - ), - cov = c( - wt = interp() - ), - sec = function() { - V = V0 * (wt/70) - ke = CL/V # define here to make eqn simpler - }, - eqn = function() { - dx[1] = -ka * x[1] - dx[2] = rateiv[1] + ka * x[1] - (ke + k23) * x[2] + k32 * x[3] - dx[3] = k23 * x[2] - k32 * x[3] - dx[4] = x[1] / V - }, - lag = function() { - tlag[1] = lag1 - }, - out = function() { - y[1] = x[1]/V - y[2] = x[4] # AUC, not fitted to any data, not required - }, - err = c( - proportional(2, c(0.1, 0.15, 0, 0)) # only applies to y[1] - ) - ) + CL = ab(10, 200), + V0 = ab(0, 100), + ka = ab(0, 3), + k23 = ab(0, 5), + k32 = ab(0, 5), + lag1 = ab(0, 2) + ), + cov = c( + wt = interp() + ), + sec = function() { + V <- V0 * (wt / 70) + ke <- CL / V # define here to make eqn simpler + }, + eqn = function() { + dx[1] <- -ka * x[1] + dx[2] <- rateiv[1] + ka * x[1] - (ke + k23) * x[2] + k32 * x[3] + dx[3] <- k23 * x[2] - k32 * x[3] + dx[4] <- x[1] / V + }, + lag = function() { + tlag[1] <- lag1 + }, + out = function() { + y[1] <- x[1] / V + y[2] <- x[4] # AUC, not fitted to any data, not required + }, + err = c( + proportional(2, c(0.1, 0.15, 0, 0)) # only applies to y[1] + ) +) \donttest{ - mod <- PM_model$new(mod_list) - } +mod <- PM_model$new(mod_list) +} } \author{ @@ -115,10 +115,10 @@ and \emph{error models}. These portions of the model have specific and defined creator functions and no additional R code is permissible. They take this form: -\if{html}{\out{
}}\preformatted{block_name = c( - var1 = creator(), +\if{html}{\out{
}}\preformatted{block_name = c( + var1 = creator(), var2 = creator() -) +) }\if{html}{\out{
}} Note the comma separating the creator functions, "\verb{c(}" to open the vector and "\verb{)}" to close the vector. @@ -128,13 +128,13 @@ equations}, \emph{model equations} (e.g. ODEs), \emph{lag time}, \emph{bioavaila and \emph{outputs}. These parts of the model are defined as R functions without arguments, but whose body contains any permissible R code. -\if{html}{\out{
}}\preformatted{block_name = function() \{ +\if{html}{\out{
}}\preformatted{block_name = function() \{ - # any valid R code + # any valid R code # can use primary or secondary parameters and covariates # lines are not separated by commas -\} +\} }\if{html}{\out{
}} Note the absence of arguments between the "\verb{()}", the opening curly brace "\verb{\{}" to start @@ -225,9 +225,9 @@ Note that \code{wt = interp()} is equivalent to \code{wt = interp("lm")}, since are not estimated for these equations but they are available to every other block in the model. For example: -\if{html}{\out{
}}\preformatted{sec = function() \{ - V = V0 * (wt/70) -\} +\if{html}{\out{
}}\preformatted{sec = function() \{ + V = V0 * (wt/70) +\} }\if{html}{\out{
}} Note that the function @@ -361,7 +361,7 @@ but will be specific to the \code{out} block. For example, -\if{html}{\out{
}}\preformatted{out = function() \{ +\if{html}{\out{
}}\preformatted{out = function() \{ V = V0 * wt # only needed if not included in sec block y[1] = x[1]/V #Vp and Vm must be defined in pri or sec blocks diff --git a/man/PM_result.Rd b/man/PM_result.Rd index c6a5e7e56..2648cc753 100755 --- a/man/PM_result.Rd +++ b/man/PM_result.Rd @@ -14,6 +14,31 @@ After a run completes, results are stored on your hard drive. They are loaded back into R with \link{PM_load} to create the \link{PM_result} object, which contains both the results and functions to analyze or plot the result. } +\examples{ + +## ------------------------------------------------ +## Method `PM_result$bestdose` +## ------------------------------------------------ + +\dontrun{ +# Load NPAG result +result <- PM_load(1) + +# Create target data +target <- PM_data$new("target.csv") + +# Run BestDose optimization +bd <- result$bestdose( + target = target, + dose_range = list(min = 50, max = 500), + bias_weight = 0.5 +) + +# View results +print(bd) +bd$get_doses() +} +} \author{ Michael Neely, Julian Otalvaro } @@ -67,6 +92,7 @@ new optimal sampling results.} \item \href{#method-PM_result-nca}{\code{PM_result$nca()}} \item \href{#method-PM_result-report}{\code{PM_result$report()}} \item \href{#method-PM_result-sim}{\code{PM_result$sim()}} +\item \href{#method-PM_result-bestdose}{\code{PM_result$bestdose()}} \item \href{#method-PM_result-save}{\code{PM_result$save()}} \item \href{#method-PM_result-validate}{\code{PM_result$validate()}} \item \href{#method-PM_result-step}{\code{PM_result$step()}} @@ -247,6 +273,85 @@ arguments, e.g. \verb{$sim(include = 1:2, predInt = 1, limits = NA)}.} } \if{html}{\out{
}} } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PM_result-bestdose}{}}} +\subsection{Method \code{bestdose()}}{ +Run BestDose optimization using this result as the prior +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{PM_result$bestdose( + target, + past_data = NULL, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + ... +)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{target}}{PM_data object or path to CSV with target doses/observations. +Required. This defines the dosing template and target values. Set dose amounts +to 0 for doses to be optimized.} + +\item{\code{past_data}}{Optional: PM_data object or path to CSV with patient history. +If NULL (default), uses the full dataset from this result.} + +\item{\code{dose_range}}{Named list with min and max allowable doses. +Default: list(min = 0, max = 1000)} + +\item{\code{bias_weight}}{Numeric \link{0,1} controlling personalization level. +0 = fully personalized, 1 = population-based. Default: 0.5 (balanced)} + +\item{\code{target_type}}{One of "concentration" (default), "auc_from_zero", or +"auc_from_last_dose"} + +\item{\code{time_offset}}{Optional: time offset for past/future concatenation} + +\item{\code{...}}{Additional arguments passed to \link{PM_bestdose}} +} +\if{html}{\out{
}} +} +\subsection{Details}{ +BestDose finds optimal dosing regimens to achieve target drug concentrations +or AUC values. The algorithm uses Bayesian posterior estimation combined with +dual optimization to balance patient-specific adaptation and population-level +robustness. By default, uses the \verb{$final}, \verb{$model}, and \verb{$data} objects from +this result. Most commonly, you will supply a different \code{target} data object. +} + +\subsection{Returns}{ +A \link{PM_bestdose} object containing optimal doses and predictions +} +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +# Load NPAG result +result <- PM_load(1) + +# Create target data +target <- PM_data$new("target.csv") + +# Run BestDose optimization +bd <- result$bestdose( + target = target, + dose_range = list(min = 50, max = 500), + bias_weight = 0.5 +) + +# View results +print(bd) +bd$get_doses() +} +} +\if{html}{\out{
}} + +} + } \if{html}{\out{
}} \if{html}{\out{}} diff --git a/man/interp.Rd b/man/interp.Rd index 0d7eb549e..6a75496a5 100755 --- a/man/interp.Rd +++ b/man/interp.Rd @@ -24,7 +24,7 @@ interpolation between values or not. } \examples{ \dontrun{ -cov = c( +cov <- c( wt = interp(), # same as interp("lm") or interp("linear") visit = interp("none") ) diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs new file mode 100644 index 000000000..e8821f0c7 --- /dev/null +++ b/src/rust/src/bestdose_executor.rs @@ -0,0 +1,214 @@ +use crate::settings::settings; +use extendr_api::prelude::*; +use pmcore::bestdose::{BestDoseProblem, BestDoseResult, DoseRange, Target}; +use pmcore::prelude::{data, ODE}; +use pmcore::routines::initialization::parse_prior; +use std::path::PathBuf; + +/// Helper to parse target type from string +pub(crate) fn parse_target_type(target_str: &str) -> std::result::Result { + match target_str.to_lowercase().as_str() { + "concentration" => Ok(Target::Concentration), + "auc_from_zero" | "auc" => Ok(Target::AUCFromZero), + "auc_from_last_dose" | "auc_interval" => Ok(Target::AUCFromLastDose), + _ => Err(format!( + "Invalid target type: {}. Must be 'concentration', 'auc_from_zero', or 'auc_from_last_dose'", + target_str + )), + } +} + +/// R-compatible prediction row for BestDose output +#[derive(Debug, IntoDataFrameRow)] +pub struct BestDosePredictionRow { + id: String, + time: f64, + observed: f64, + pop_mean: f64, + pop_median: f64, + post_mean: f64, + post_median: f64, + outeq: usize, +} + +impl BestDosePredictionRow { + pub fn from_np_prediction( + pred: &pmcore::routines::output::predictions::NPPredictionRow, + id: &str, + ) -> Self { + Self { + id: id.to_string(), + time: pred.time(), + observed: pred.obs().unwrap_or(0.0), + pop_mean: pred.pop_mean(), + pop_median: pred.pop_median(), + post_mean: pred.post_mean(), + post_median: pred.post_median(), + outeq: pred.outeq(), + } + } +} + +/// R-compatible AUC prediction row +#[derive(Debug, IntoDataFrameRow)] +pub struct BestDoseAucRow { + time: f64, + auc: f64, +} + +/// Convert BestDoseResult to R-compatible list structure +pub(crate) fn convert_bestdose_result_to_r( + result: BestDoseResult, +) -> std::result::Result { + // Extract doses + let doses: Vec = result.doses(); + + // Objective function + let objf = result.objf(); + + // Status + let status_str = format!("{:?}", result.status()); + + // Predictions as data frame + let pred_rows: Vec = result + .predictions() + .predictions() + .iter() + .map(|p| BestDosePredictionRow::from_np_prediction(p, "subject_1")) + .collect(); + let pred_df = pred_rows + .into_dataframe() + .map_err(|e| format!("Failed to create predictions dataframe: {:?}", e))?; + + // AUC predictions (if available) + let auc_val = if let Some(auc_preds) = result.auc_predictions() { + let auc_rows: Vec = auc_preds + .iter() + .map(|(time, auc)| BestDoseAucRow { + time: *time, + auc: *auc, + }) + .collect(); + let auc_df = auc_rows + .into_dataframe() + .map_err(|e| format!("Failed to create AUC dataframe: {:?}", e))?; + Robj::from(auc_df) + } else { + Robj::from(()) // NULL for no AUC + }; + + // Optimization method + let method_str = format!("{}", result.optimization_method()); + + // Build the list using list! macro + let output = list!( + doses = doses, + objf = objf, + status = status_str, + predictions = pred_df, + auc_predictions = auc_val, + method = method_str + ); + + Ok(output.into()) +} + +/// Execute bestdose optimization for ODE models +pub(crate) fn bestdose_ode( + model_path: PathBuf, + prior_path: PathBuf, + past_data_path: Option, + target_data_path: PathBuf, + time_offset: Option, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, +) -> std::result::Result { + // 1. Load the model + let (_lib, (eq, meta)) = + unsafe { pmcore::prelude::pharmsol::exa::load::load::(model_path) }; + + // 2. Parse settings from R list + let settings = settings(params, meta.get_params(), "/tmp/bestdose") + .map_err(|e| format!("Failed to parse settings: {}", e))?; + + // 3. Parse the prior (theta + weights) + let (population_theta, prior_weights) = parse_prior( + &prior_path.to_str().unwrap().to_string(), + &settings, + ) + .map_err(|e| format!("Failed to parse prior: {}", e))?; + + let population_weights = prior_weights.ok_or_else(|| { + "Prior file must contain a 'prob' column with weights".to_string() + })?; + + // 4. Load past data (if provided) + let past_data = if let Some(path) = past_data_path { + let data = data::read_pmetrics(path.to_str().unwrap()) + .map_err(|e| format!("Failed to read past data: {}", e))?; + let subjects = data.subjects(); + if subjects.is_empty() { + return Err("Past data file contains no subjects".to_string()); + } + Some(subjects[0].clone()) + } else { + None + }; + + // 5. Load target data + let target_data = { + let data = data::read_pmetrics(target_data_path.to_str().unwrap()) + .map_err(|e| format!("Failed to read target data: {}", e))?; + let subjects = data.subjects(); + if subjects.is_empty() { + return Err("Target data file contains no subjects".to_string()); + } + subjects[0].clone() + }; + + // 6. Parse target type + let target_enum = parse_target_type(target_type)?; + + // 7. Create dose range + let doserange = DoseRange::new(dose_min, dose_max); + + // 8. Create and solve bestdose problem + let problem = BestDoseProblem::new( + &population_theta, + &population_weights, + past_data, + target_data, + time_offset, + eq, + doserange, + bias_weight, + settings, + target_enum, + ) + .map_err(|e| format!("Failed to create BestDose problem: {}", e))?; + + let result = problem + .optimize() + .map_err(|e| format!("Optimization failed: {}", e))?; + + Ok(result) +} + +/// Execute bestdose optimization for analytical models (placeholder - not yet supported) +pub(crate) fn bestdose_analytical( + _model_path: PathBuf, + _prior_path: PathBuf, + _past_data_path: Option, + _target_data_path: PathBuf, + _time_offset: Option, + _dose_min: f64, + _dose_max: f64, + _bias_weight: f64, + _target_type: &str, + _params: List, +) -> std::result::Result { + Err("BestDose for analytical models is not yet supported".to_string()) +} diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs index 3ee6cd5a4..09af6a5b6 100755 --- a/src/rust/src/lib.rs +++ b/src/rust/src/lib.rs @@ -1,4 +1,5 @@ // mod build; +mod bestdose_executor; mod executor; mod logs; mod settings; @@ -8,6 +9,7 @@ use anyhow::Result; use extendr_api::prelude::*; use pmcore::prelude::{data::read_pmetrics, pharmsol::exa::build, Analytical, ODE}; use simulation::SimulationRow; +use std::path::PathBuf; use std::process::Command; use tracing_subscriber::layer::SubscriberExt; @@ -247,6 +249,73 @@ fn setup_logs() -> anyhow::Result<()> { Ok(()) } +/// Run BestDose optimization to find optimal doses +///@export +#[extendr] +fn bestdose( + model_path: &str, + prior_path: &str, + past_data_path: Nullable, + target_data_path: &str, + time_offset: Nullable, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + kind: &str, +) -> Robj { + RFormatLayer::reset_global_timer(); + let _ = setup_logs(); + + println!("Starting BestDose optimization..."); + + let past_path = match past_data_path.into_option() { + Some(p) => Some(PathBuf::from(p)), + None => None, + }; + + let time_offset_opt = time_offset.into_option(); + + let result = match kind { + "ode" => bestdose_executor::bestdose_ode( + model_path.into(), + prior_path.into(), + past_path, + target_data_path.into(), + time_offset_opt, + dose_min, + dose_max, + bias_weight, + target_type, + params.clone(), + ), + "analytical" => bestdose_executor::bestdose_analytical( + model_path.into(), + prior_path.into(), + past_path, + target_data_path.into(), + time_offset_opt, + dose_min, + dose_max, + bias_weight, + target_type, + params.clone(), + ), + _ => { + return Robj::from(format!("{} is not a supported model type", kind)); + } + }; + + match result { + Ok(bd_result) => match bestdose_executor::convert_bestdose_result_to_r(bd_result) { + Ok(r) => r, + Err(e) => Robj::from(format!("Failed to convert result: {}", e)), + }, + Err(e) => Robj::from(format!("BestDose failed: {}", e)), + } +} + // Macro to generate exports. // This ensures exported functions are registered with R. // See corresponding C code in `entrypoint.c`. @@ -260,6 +329,7 @@ extendr_module! { fn fit; fn model_parameters; fn setup_logs; + fn bestdose; } // To generate the exported function in R, run the following command: From 7304c03bb312c22fe6a9349bef9d3a9ed4806cc1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Wed, 5 Nov 2025 16:54:41 +0000 Subject: [PATCH 2/5] working example --- R/PM_bestdose.R | 702 +++++++++++++++--------------- R/PM_result.R | 108 ++--- src/rust/src/bestdose_executor.rs | 15 +- 3 files changed, 420 insertions(+), 405 deletions(-) diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R index 5040fccba..b18740676 100644 --- a/R/PM_bestdose.R +++ b/R/PM_bestdose.R @@ -1,359 +1,375 @@ #' @title #' Object to contain BestDose optimization results -#' +#' #' @description #' `r lifecycle::badge("experimental")` -#' +#' #' This object is created after a successful BestDose optimization run. #' BestDose finds optimal dosing regimens to achieve target drug concentrations #' or AUC values using Bayesian optimization. -#' +#' #' @export PM_bestdose <- R6::R6Class( - "PM_bestdose", - public = list( - #' @field result A list containing: - #' * **doses** Vector of optimal dose amounts (mg) - #' * **objf** Objective function value (lower is better) - #' * **status** Optimization status ("converged" or "max_iterations") - #' * **predictions** Data frame with concentration-time predictions - #' * **auc_predictions** Data frame with AUC predictions (if applicable) - #' * **method** Optimization method used ("posterior" or "uniform") - result = NULL, - - #' @description - #' Create a new PM_bestdose object by running BestDose optimization - #' - #' @param prior A PM_result or PM_final object, or path to theta.csv file - #' @param model A PM_model object or path to compiled model - #' @param past_data Optional: PM_data object or path to CSV with patient history - #' @param target PM_data object or path to CSV with target doses/observations - #' @param dose_range Named list with min and max allowable doses, e.g., - #' list(min = 0, max = 1000). Default: list(min = 0, max = 1000) - #' @param bias_weight Numeric [0,1]. 0 = fully personalized, 1 = population-based. - #' Default = 0.5 (balanced) - #' @param target_type One of "concentration", "auc_from_zero", or "auc_from_last_dose". - #' Default = "concentration" - #' @param time_offset Optional: time offset for past/future concatenation - #' @param settings Optional: list of settings (algorithm, error models, etc.) - #' If NULL, will use defaults from prior object - #' - #' @details - #' - #' ## Target Data Format - #' - #' The target data should be a Pmetrics-formatted CSV with: - #' - Dose events: Set dose amount to 0 for doses to be optimized, - #' or non-zero for fixed doses - #' - Observation events: Set OUT values to the target concentrations or AUCs - #' - #' ## Target Types - #' - #' - **concentration**: Optimize to match target concentrations at observation times - #' - **auc_from_zero**: Optimize to match cumulative AUC from time 0 - #' - **auc_from_last_dose**: Optimize to match interval AUC from last dose - #' - #' ## Bias Weight - #' - #' Controls the balance between patient-specific and population-based optimization: - #' - 0.0: Fully personalized (minimize variance) - #' - 0.5: Balanced approach (recommended default) - #' - 1.0: Population-based (minimize bias from population) - #' - #' @examples - #' \dontrun{ - #' # Load NPAG result - #' result <- PM_load(1) - #' - #' # Create target data - #' target <- PM_data$new("target.csv") - #' - #' # Run BestDose optimization - #' bd <- PM_bestdose$new( - #' prior = result, - #' model = result$model, - #' target = target, - #' dose_range = list(min = 50, max = 500), - #' bias_weight = 0.5 - #' ) - #' - #' # View results - #' print(bd) - #' bd$get_doses() - #' bd$get_predictions() - #' } - initialize = function(prior, - model, - past_data = NULL, - target, - dose_range = list(min = 0, max = 1000), - bias_weight = 0.5, - target_type = "concentration", - time_offset = NULL, - settings = NULL) { - - # Validate inputs - if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) { - cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose") - } - - if (bias_weight < 0 || bias_weight > 1) { - cli::cli_abort("bias_weight must be between 0 and 1") - } - - if (is.null(dose_range$min) || is.null(dose_range$max)) { - cli::cli_abort("dose_range must have both 'min' and 'max' elements") - } - - if (dose_range$min >= dose_range$max) { - cli::cli_abort("dose_range$min must be less than dose_range$max") - } - - # Parse prior - prior_path <- private$.parse_prior(prior) - - # Parse model - model_info <- private$.parse_model(model) - - # Parse data files - past_data_path <- if (!is.null(past_data)) { - private$.parse_data(past_data) - } else { - NULL - } - - target_data_path <- private$.parse_data(target) - - # Parse settings - if (is.null(settings)) { - settings <- private$.default_bestdose_settings(prior) - } - - # Call Rust function - cli::cli_alert_info("Running BestDose optimization...") - - result <- tryCatch({ - bestdose( - model_path = model_info$path, - prior_path = prior_path, - past_data_path = past_data_path, - target_data_path = target_data_path, - time_offset = time_offset, - dose_min = dose_range$min, - dose_max = dose_range$max, - bias_weight = bias_weight, - target_type = target_type, - params = settings, - kind = model_info$kind - ) - }, error = function(e) { - cli::cli_abort(c( - "x" = "BestDose optimization failed", - "i" = conditionMessage(e) - )) - }) - - self$result <- result - - cli::cli_alert_success("Optimization complete!") - cli::cli_alert_info("Method: {result$method}") - cli::cli_alert_info("Status: {result$status}") - cli::cli_alert_info("Optimal doses: {paste(round(result$doses, 2), collapse = ', ')} mg") - - invisible(self) - }, - - #' @description - #' Print summary of BestDose results - print = function() { - cat("BestDose Optimization Results\n") - cat("==============================\n\n") - cat(sprintf("Optimal doses: %s mg\n", paste(round(self$result$doses, 2), collapse = ", "))) - cat(sprintf("Objective function: %.6f\n", self$result$objf)) - cat(sprintf("Status: %s\n", self$result$status)) - cat(sprintf("Method: %s\n", self$result$method)) - cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions))) - if (!is.null(self$result$auc_predictions)) { - cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions))) - } - invisible(self) - }, - - #' @description - #' Get optimal dose values - #' @return Numeric vector of optimal doses - get_doses = function() { - self$result$doses - }, - - #' @description - #' Get concentration-time predictions - #' @return Data frame with predictions - get_predictions = function() { - self$result$predictions - }, - - #' @description - #' Get AUC predictions (if available) - #' @return Data frame with AUC predictions or NULL - get_auc_predictions = function() { - self$result$auc_predictions - }, - - #' @description - #' Get objective function value - #' @return Numeric objective function value - get_objf = function() { - self$result$objf - }, - - #' @description - #' Get optimization status - #' @return Character string with status - get_status = function() { - self$result$status - }, - - #' @description - #' Get optimization method used - #' @return Character string: "posterior" or "uniform" - get_method = function() { - self$result$method - }, - - #' @description - #' Save results to RDS file - #' @param filename Path to save file. Default: "bestdose_result.rds" - save = function(filename = "bestdose_result.rds") { - saveRDS(self, filename) - cli::cli_alert_success("Results saved to {filename}") - invisible(self) - } - ), - - private = list( - # Helper function to parse prior - .parse_prior = function(prior) { - if (inherits(prior, "PM_result")) { - # Extract theta.csv from PM_result - theta_path <- file.path(prior$rundir, "outputs", "theta.csv") - if (!file.exists(theta_path)) { - cli::cli_abort("theta.csv not found in PM_result outputs") - } - return(theta_path) - } else if (inherits(prior, "PM_final")) { - # Need to write out PM_final to CSV - temp_path <- tempfile(fileext = ".csv") - private$.write_prior_csv(prior, temp_path) - return(temp_path) - } else if (is.character(prior)) { - if (!file.exists(prior)) { - cli::cli_abort("Prior file not found: {prior}") - } - return(prior) - } else { - cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv") - } - }, - - # Helper function to parse model - .parse_model = function(model) { - if (inherits(model, "PM_model")) { - # Check if model is compiled - compiled_path <- model$model_list$compiled_location - if (is.null(compiled_path) || !file.exists(compiled_path)) { - cli::cli_abort("Model must be compiled first. Use model$compile()") - } - - # Determine model type from model_list - kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) { - "analytical" - } else { - "ode" - } - - return(list( - path = compiled_path, - kind = kind - )) - } else if (is.character(model)) { - # Assume it's a path to compiled model - if (!file.exists(model)) { - cli::cli_abort("Model file not found: {model}") - } - # Try to infer type from extension or filename - kind <- if (grepl("analytical", model, ignore.case = TRUE)) { - "analytical" - } else { - "ode" + "PM_bestdose", + public = list( + #' @field result A list containing: + #' * **doses** Vector of optimal dose amounts (mg) + #' * **objf** Objective function value (lower is better) + #' * **status** Optimization status ("converged" or "max_iterations") + #' * **predictions** Data frame with concentration-time predictions + #' * **auc_predictions** Data frame with AUC predictions (if applicable) + #' * **method** Optimization method used ("posterior" or "uniform") + result = NULL, + + #' @description + #' Create a new PM_bestdose object by running BestDose optimization + #' + #' @param prior A PM_result or PM_final object, or path to theta.csv file + #' @param model A PM_model object or path to compiled model + #' @param past_data Optional: PM_data object or path to CSV with patient history + #' @param target PM_data object or path to CSV with target doses/observations + #' @param dose_range Named list with min and max allowable doses, e.g., + #' list(min = 0, max = 1000). Default: list(min = 0, max = 1000) + #' @param bias_weight Numeric [0,1]. 0 = fully personalized, 1 = population-based. + #' Default = 0.5 (balanced) + #' @param target_type One of "concentration", "auc_from_zero", or "auc_from_last_dose". + #' Default = "concentration" + #' @param time_offset Optional: time offset for past/future concatenation + #' @param settings Optional: list of settings (algorithm, error models, etc.) + #' If NULL, will use defaults from prior object + #' + #' @details + #' + #' ## Target Data Format + #' + #' The target data should be a Pmetrics-formatted CSV with: + #' - Dose events: Set dose amount to 0 for doses to be optimized, + #' or non-zero for fixed doses + #' - Observation events: Set OUT values to the target concentrations or AUCs + #' + #' ## Target Types + #' + #' - **concentration**: Optimize to match target concentrations at observation times + #' - **auc_from_zero**: Optimize to match cumulative AUC from time 0 + #' - **auc_from_last_dose**: Optimize to match interval AUC from last dose + #' + #' ## Bias Weight + #' + #' Controls the balance between patient-specific and population-based optimization: + #' - 0.0: Fully personalized (minimize variance) + #' - 0.5: Balanced approach (recommended default) + #' - 1.0: Population-based (minimize bias from population) + #' + #' @examples + #' \dontrun{ + #' # Load NPAG result + #' result <- PM_load(1) + #' + #' # Create target data + #' target <- PM_data$new("target.csv") + #' + #' # Run BestDose optimization + #' bd <- PM_bestdose$new( + #' prior = result, + #' model = result$model, + #' target = target, + #' dose_range = list(min = 50, max = 500), + #' bias_weight = 0.5 + #' ) + #' + #' # View results + #' print(bd) + #' bd$get_doses() + #' bd$get_predictions() + #' } + initialize = function(prior, + model, + past_data = NULL, + target, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + settings = NULL) { + # Validate inputs + if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) { + cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose") + } + + if (bias_weight < 0 || bias_weight > 1) { + cli::cli_abort("bias_weight must be between 0 and 1") + } + + if (is.null(dose_range$min) || is.null(dose_range$max)) { + cli::cli_abort("dose_range must have both 'min' and 'max' elements") + } + + if (dose_range$min >= dose_range$max) { + cli::cli_abort("dose_range$min must be less than dose_range$max") + } + + # Parse prior + prior_path <- private$.parse_prior(prior) + + # Parse model + model_info <- private$.parse_model(model) + + # Parse data files + past_data_path <- if (!is.null(past_data)) { + private$.parse_data(past_data) + } else { + NULL + } + + target_data_path <- private$.parse_data(target) + + # Parse settings + if (is.null(settings)) { + settings <- private$.default_bestdose_settings(prior, model) + } + + # Call Rust function + cli::cli_alert_info("Running BestDose optimization...") + + result <- tryCatch( + { + bestdose( + model_path = model_info$path, + prior_path = prior_path, + past_data_path = past_data_path, + target_data_path = target_data_path, + time_offset = time_offset, + dose_min = dose_range$min, + dose_max = dose_range$max, + bias_weight = bias_weight, + target_type = target_type, + params = settings, + kind = model_info$kind + ) + }, + error = function(e) { + cli::cli_abort(c( + "x" = "BestDose optimization failed", + "i" = conditionMessage(e) + )) + } + ) + + self$result <- result + + cli::cli_alert_success("Optimization complete!") + cat("Result names:", names(result), "\n") + cat("Result structure:\n") + print(str(result)) + if ("method" %in% names(result)) { + cli::cli_alert_info("Method: {result[[\"method\"]]}") + } + if ("status" %in% names(result)) { + cli::cli_alert_info("Status: {result[[\"status\"]]}") + } + if ("doses" %in% names(result)) { + cli::cli_alert_info("Optimal doses: {paste(round(result[[\"doses\"]], 2), collapse = ', ')} mg") + } + + invisible(self) + }, + + #' @description + #' Print summary of BestDose results + print = function() { + cat("BestDose Optimization Results\n") + cat("==============================\n\n") + cat(sprintf("Optimal doses: %s mg\n", paste(round(self$result$doses, 2), collapse = ", "))) + cat(sprintf("Objective function: %.6f\n", self$result$objf)) + cat(sprintf("Status: %s\n", self$result$status)) + cat(sprintf("Method: %s\n", self$result$method)) + cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions))) + if (!is.null(self$result$auc_predictions)) { + cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions))) + } + invisible(self) + }, + + #' @description + #' Get optimal dose values + #' @return Numeric vector of optimal doses + get_doses = function() { + self$result$doses + }, + + #' @description + #' Get concentration-time predictions + #' @return Data frame with predictions + get_predictions = function() { + self$result$predictions + }, + + #' @description + #' Get AUC predictions (if available) + #' @return Data frame with AUC predictions or NULL + get_auc_predictions = function() { + self$result$auc_predictions + }, + + #' @description + #' Get objective function value + #' @return Numeric objective function value + get_objf = function() { + self$result$objf + }, + + #' @description + #' Get optimization status + #' @return Character string with status + get_status = function() { + self$result$status + }, + + #' @description + #' Get optimization method used + #' @return Character string: "posterior" or "uniform" + get_method = function() { + self$result$method + }, + + #' @description + #' Save results to RDS file + #' @param filename Path to save file. Default: "bestdose_result.rds" + save = function(filename = "bestdose_result.rds") { + saveRDS(self, filename) + cli::cli_alert_success("Results saved to {filename}") + invisible(self) } - return(list(path = model, kind = kind)) - } else { - cli::cli_abort("model must be PM_model or path to compiled model") - } - }, - - # Helper function to parse data - .parse_data = function(data) { - if (inherits(data, "PM_data")) { - # Write to temp CSV - temp_path <- tempfile(fileext = ".csv") - write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE) - return(temp_path) - } else if (is.character(data)) { - if (!file.exists(data)) { - cli::cli_abort("Data file not found: {data}") + ), + private = list( + # Helper function to parse prior + .parse_prior = function(prior) { + if (inherits(prior, "PM_result")) { + # Extract theta.csv from PM_result + theta_path <- file.path(prior$rundir, "outputs", "theta.csv") + if (!file.exists(theta_path)) { + cli::cli_abort("theta.csv not found in PM_result outputs") + } + return(theta_path) + } else if (inherits(prior, "PM_final")) { + # Need to write out PM_final to CSV + temp_path <- tempfile(fileext = ".csv") + private$.write_prior_csv(prior, temp_path) + return(temp_path) + } else if (is.character(prior)) { + if (!file.exists(prior)) { + cli::cli_abort("Prior file not found: {prior}") + } + return(prior) + } else { + cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv") + } + }, + + # Helper function to parse model + .parse_model = function(model) { + if (inherits(model, "PM_model")) { + # Check if model is compiled + compiled_path <- model$binary_path + if (is.null(compiled_path) || !file.exists(compiled_path)) { + cli::cli_abort("Model must be compiled first. Use model$compile()") + } + + # Determine model type from model_list + kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) { + "analytical" + } else { + "ode" + } + + return(list( + path = compiled_path, + kind = kind + )) + } else if (is.character(model)) { + # Assume it's a path to compiled model + if (!file.exists(model)) { + cli::cli_abort("Model file not found: {model}") + } + # Try to infer type from extension or filename + kind <- if (grepl("analytical", model, ignore.case = TRUE)) { + "analytical" + } else { + "ode" + } + return(list(path = model, kind = kind)) + } else { + cli::cli_abort("model must be PM_model or path to compiled model") + } + }, + + # Helper function to parse data + .parse_data = function(data) { + if (inherits(data, "PM_data")) { + # Write to temp CSV + temp_path <- tempfile(fileext = ".csv") + write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE) + return(temp_path) + } else if (is.character(data)) { + if (!file.exists(data)) { + cli::cli_abort("Data file not found: {data}") + } + return(data) + } else { + cli::cli_abort("data must be PM_data or path to CSV file") + } + }, + + # Helper to create default settings + .default_bestdose_settings = function(prior, model) { + # Extract settings from prior if it's a PM_result + if (inherits(prior, "PM_result")) { + # Use the settings from the result + settings <- prior$settings + return(settings) + } else { + # Get parameter ranges from model (same format as in fit() method) + param_ranges <- lapply(model$model_list$pri, function(x) { + c(x$min, x$max) + }) + names(param_ranges) <- tolower(names(param_ranges)) + + # Use sensible defaults + settings <- list( + algorithm = "NPAG", + ranges = param_ranges, + error_models = list( + list( + initial = 0.0, + type = "additive", + coeff = c(0.0, 0.2, 0.0, 0.0) + ) + ), + max_cycles = 500, + points = 2028, + seed = 22, + prior = "prior.csv", + idelta = 0.25, # 15 minutes for AUC calculations + tad = 0.0 + ) + } + return(settings) + }, + + # Helper to write PM_final to CSV + .write_prior_csv = function(prior, path) { + # Convert PM_final to CSV with prob column + df <- as.data.frame(prior$popPoints) + df$prob <- prior$popProb + write.csv(df, path, row.names = FALSE, quote = FALSE) } - return(data) - } else { - cli::cli_abort("data must be PM_data or path to CSV file") - } - }, - - # Helper to create default settings - .default_bestdose_settings = function(prior) { - # Extract settings from prior if it's a PM_result - if (inherits(prior, "PM_result")) { - # Use the settings from the result - settings <- prior$settings - return(settings) - } else { - # Use sensible defaults - settings <- list( - algorithm = "NPAG", - ranges = list(), # Will be filled from prior - error_models = list( - list( - initial = 0.0, - type = "additive", - coeff = c(0.0, 0.2, 0.0, 0.0) - ) - ), - max_cycles = 500, - points = 2028, - seed = 22, - prior = "prior.csv", - idelta = 0.25, # 15 minutes for AUC calculations - tad = 0.0 - ) - } - return(settings) - }, - - # Helper to write PM_final to CSV - .write_prior_csv = function(prior, path) { - # Convert PM_final to CSV with prob column - df <- as.data.frame(prior$popPoints) - df$prob <- prior$popProb - write.csv(df, path, row.names = FALSE, quote = FALSE) - } - ) + ) ) #' @export PM_bestdose$load <- function(filename = "bestdose_result.rds") { - if (!file.exists(filename)) { - cli::cli_abort("File not found: {filename}") - } - readRDS(filename) + if (!file.exists(filename)) { + cli::cli_abort("File not found: {filename}") + } + readRDS(filename) } diff --git a/R/PM_result.R b/R/PM_result.R index 93474777a..a2c4aaa64 100755 --- a/R/PM_result.R +++ b/R/PM_result.R @@ -54,7 +54,7 @@ PM_result <- R6::R6Class( #' Use the `$save` method on the augmented `PM_result` object to save it with the #' new optimal sampling results. opt_samp = NULL, - + #' @description #' Create new object populated with data from previous run #' @details @@ -73,9 +73,9 @@ PM_result <- R6::R6Class( if (!inherits(out[[x]], "R6")) { # older save cli::cli_abort(c("x" = "The object was saved in an older format. Please re-run the analysis.")) } else { - if(x == "model"){ - args <- list(x = out[[x]], compile = FALSE) - } else { + if (x == "model") { + args <- list(x = out[[x]], compile = FALSE) + } else { args <- list(out[[x]], path = path, quiet = TRUE) } self[[x]] <- do.call(get(paste0("PM_", x))$new, args = args) # was saved in R6 format, but remake to update if needed @@ -83,24 +83,24 @@ PM_result <- R6::R6Class( } } ) - + # these are diagnostics, not R6 self$errfile <- out$errfile self$success <- out$success - + # add the pop/post data to data if (is.null(self$data$pop) | is.null(self$data$post)) { self$data <- PM_data$new(self$data$data, quiet = TRUE) self$data$pop <- self$pop$data self$data$post <- self$post$data } - + return(self) }, #' @description #' Fit the model to the data #' #' @details - #' This method is used to fit the model in the [PM_result] object to data. + #' This method is used to fit the model in the [PM_result] object to data. #' It calls the `$fit` method of the model stored in the `model` field. #' @param data Optional data to fit. If not provided, the data stored in the #' `data` field of the [PM_result] object will be used. This can be useful to @@ -109,15 +109,15 @@ PM_result <- R6::R6Class( #' @param ... Additional arguments passed to the model's `$fit` method. #' @return Returns an invisible [PM_result]. #' @export - #' - fit = function(data, ...){ + #' + fit = function(data, ...) { if (missing(data)) { data <- self$data } res <- self$model$fit(data = data, ...) return(invisible(res)) }, - + #' @description #' Plot generic function based on type #' @param type Type of plot based on class of object @@ -129,7 +129,7 @@ PM_result <- R6::R6Class( self[[type]]$plot(...) } }, - + #' @description #' Summary generic function based on type #' @param type Type of summary based on class of object @@ -141,7 +141,7 @@ PM_result <- R6::R6Class( self[[type]]$summary(...) } }, - + #' @description #' AUC generic function based on type #' @param type Type of AUC based on class of object @@ -152,7 +152,7 @@ PM_result <- R6::R6Class( } self[[type]]$auc(...) }, - + #' @description #' Perform non-compartmental analysis #' @details @@ -180,22 +180,22 @@ PM_result <- R6::R6Class( if (!"poppar" %in% names(dots)) { dots$poppar <- self$final } - + if (!"data" %in% names(dots)) { dots$data <- self$data } - + if (!"model" %in% names(dots)) { dots$model <- self$model } - + # store copy of the final object bk_final <- self$final$clone() sim <- do.call(PM_sim$new, dots) self$final <- bk_final return(sim) }, - + #' @description #' Run BestDose optimization using this result as the prior #' @details @@ -205,15 +205,15 @@ PM_result <- R6::R6Class( #' robustness. By default, uses the `$final`, `$model`, and `$data` objects from #' this result. Most commonly, you will supply a different `target` data object. #' @param target PM_data object or path to CSV with target doses/observations. - #' Required. This defines the dosing template and target values. Set dose amounts + #' Required. This defines the dosing template and target values. Set dose amounts #' to 0 for doses to be optimized. #' @param past_data Optional: PM_data object or path to CSV with patient history. #' If NULL (default), uses the full dataset from this result. - #' @param dose_range Named list with min and max allowable doses. + #' @param dose_range Named list with min and max allowable doses. #' Default: list(min = 0, max = 1000) #' @param bias_weight Numeric [0,1] controlling personalization level. #' 0 = fully personalized, 1 = population-based. Default: 0.5 (balanced) - #' @param target_type One of "concentration" (default), "auc_from_zero", or + #' @param target_type One of "concentration" (default), "auc_from_zero", or #' "auc_from_last_dose" #' @param time_offset Optional: time offset for past/future concatenation #' @param ... Additional arguments passed to [PM_bestdose] @@ -222,17 +222,17 @@ PM_result <- R6::R6Class( #' \dontrun{ #' # Load NPAG result #' result <- PM_load(1) - #' + #' #' # Create target data #' target <- PM_data$new("target.csv") - #' + #' #' # Run BestDose optimization #' bd <- result$bestdose( #' target = target, #' dose_range = list(min = 50, max = 500), #' bias_weight = 0.5 #' ) - #' + #' #' # View results #' print(bd) #' bd$get_doses() @@ -248,7 +248,7 @@ PM_result <- R6::R6Class( if (is.null(past_data)) { past_data <- self$data } - + PM_bestdose$new( prior = self, model = self$model, @@ -261,7 +261,7 @@ PM_result <- R6::R6Class( ... ) }, - + #' @description #' Save the current PM_result object to an .Rdata file. #' @details @@ -282,7 +282,7 @@ PM_result <- R6::R6Class( #' your current working directory, specify `run = 1` to save the result to the "outputs" #' subfolder of the "1" folder. #' @param file Custom file name. Default is "PMout.Rdata". If `run` is not specified, `file` - #' should be the full path and filename. + #' should be the full path and filename. save = function(run, file = "PMout.Rdata") { if (missing(run)) { cli::cli_inform(c( @@ -309,7 +309,7 @@ PM_result <- R6::R6Class( ) save(PMout, file = paste0(outputfolder, "/", file)) }, - + #' @description #' Validate the result by internal simulation methods. #' @param ... Arguments passed to [PM_valid]. @@ -321,9 +321,8 @@ PM_result <- R6::R6Class( " " = "For example, if your results are in {.code my_run}, use {.code my_run$save(1)} to save back to the outputs folder of run 1." )) return(invisible(self)) - }, - + #' @description #' Conduct stepwise linear regression of Bayesian posterior parameter values #' and covariates. @@ -331,7 +330,7 @@ PM_result <- R6::R6Class( step = function(...) { PM_step(self$cov$data, ...) }, - + #' @description #' Calculate optimal sampling times. #' @@ -346,7 +345,7 @@ PM_result <- R6::R6Class( }) return(invisible(self)) }, - + #' @description #' `r lifecycle::badge("deprecated")` #' @@ -410,7 +409,7 @@ PM_result$load <- function(...) { #' run4 <- PM_load(file = "Pmetrics/MyRuns/4/outputs/PMout.Rdata") # loads from Pmetrics/MyRuns/4/outputs/PMout.Rdata #' run5 <- PM_load() # loads from ./PMout.Rdata #' } -#' +#' #' @author Michael Neely and Julian Otalvaro #' @seealso [PM_final], #' [PM_cycle], [PM_op], [PM_cov], @@ -427,20 +426,22 @@ PM_load <- function(run, path = ".", file = "PMout.Rdata") { names(aux_list) <- names(Out)[i] result <- append(result, aux_list) } - + return(result) } - + found <- "" # initialize - + if (!missing(run)) { filepath <- file.path(path, run, "outputs", file) } else { filepath <- file.path(path, file) - } - - if (file.exists(filepath)) { found <- filepath } - + } + + if (file.exists(filepath)) { + found <- filepath + } + if (found != "") { result <- output2List(Out = get(load(found))) rebuild <- PM_result$new(result, path = dirname(found), quiet = TRUE) @@ -464,31 +465,32 @@ update <- function(res, found) { # start conversion n_cyc <- nrow(dat$mean) n_out <- max(res$op$outeq) - dat$gamlam <- dat$gamlam %>% select(starts_with("add")|starts_with("prop")) + dat$gamlam <- dat$gamlam %>% select(starts_with("add") | starts_with("prop")) if (ncol(gamlam) == 1 & n_out > 1) { gamlam <- cbind(gamlam, replicate((n_out - 1), gamlam[, 1])) } gamlam <- gamlam %>% - pivot_longer( - cols = everything(), - values_to = "value", names_to = c("type", "outeq"), - names_sep = "\\." - ) %>% - mutate(cycle = rep(1:n_cyc, each = n_out)) %>% - select(cycle, value, outeq, type) %>% arrange(cycle, outeq) + pivot_longer( + cols = everything(), + values_to = "value", names_to = c("type", "outeq"), + names_sep = "\\." + ) %>% + mutate(cycle = rep(1:n_cyc, each = n_out)) %>% + select(cycle, value, outeq, type) %>% + arrange(cycle, outeq) if (is.matrix(dat$mean)) { # old fortran format, but not rust format dat$mean <- tibble::tibble(cycle = 1:n_cyc) %>% - dplyr::bind_cols(tidyr::as_tibble(dat$mean)) + dplyr::bind_cols(tidyr::as_tibble(dat$mean)) dat$median <- tibble::tibble(cycle = 1:n_cyc) %>% - dplyr::bind_cols(tidyr::as_tibble(dat$median)) + dplyr::bind_cols(tidyr::as_tibble(dat$median)) dat$sd <- tibble::tibble(cycle = 1:n_cyc) %>% - dplyr::bind_cols(tidyr::as_tibble(dat$sd)) + dplyr::bind_cols(tidyr::as_tibble(dat$sd)) } msg <- c(msg, "cycle") res$cycle <- dat } } - + ####### DONE PROCESSING, INFORM ######### if (!is.null(msg)) { cat( @@ -516,6 +518,6 @@ update <- function(res, found) { cat("Results saved\n") } } - + return(res) } diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs index e8821f0c7..8caee71ed 100644 --- a/src/rust/src/bestdose_executor.rs +++ b/src/rust/src/bestdose_executor.rs @@ -135,15 +135,12 @@ pub(crate) fn bestdose_ode( .map_err(|e| format!("Failed to parse settings: {}", e))?; // 3. Parse the prior (theta + weights) - let (population_theta, prior_weights) = parse_prior( - &prior_path.to_str().unwrap().to_string(), - &settings, - ) - .map_err(|e| format!("Failed to parse prior: {}", e))?; - - let population_weights = prior_weights.ok_or_else(|| { - "Prior file must contain a 'prob' column with weights".to_string() - })?; + let (population_theta, prior_weights) = + parse_prior(&prior_path.to_str().unwrap().to_string(), &settings) + .map_err(|e| format!("Failed to parse prior: {}", e))?; + + let population_weights = prior_weights + .ok_or_else(|| "Prior file must contain a 'prob' column with weights".to_string())?; // 4. Load past data (if provided) let past_data = if let Some(path) = past_data_path { From 15793e42f59b015cd86099b6e98d7ecc6a830343 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Wed, 5 Nov 2025 19:25:53 +0000 Subject: [PATCH 3/5] working example --- R/PM_bestdose.R | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R index b18740676..e455b9cd2 100644 --- a/R/PM_bestdose.R +++ b/R/PM_bestdose.R @@ -158,19 +158,6 @@ PM_bestdose <- R6::R6Class( self$result <- result cli::cli_alert_success("Optimization complete!") - cat("Result names:", names(result), "\n") - cat("Result structure:\n") - print(str(result)) - if ("method" %in% names(result)) { - cli::cli_alert_info("Method: {result[[\"method\"]]}") - } - if ("status" %in% names(result)) { - cli::cli_alert_info("Status: {result[[\"status\"]]}") - } - if ("doses" %in% names(result)) { - cli::cli_alert_info("Optimal doses: {paste(round(result[[\"doses\"]], 2), collapse = ', ')} mg") - } - invisible(self) }, @@ -179,10 +166,11 @@ PM_bestdose <- R6::R6Class( print = function() { cat("BestDose Optimization Results\n") cat("==============================\n\n") - cat(sprintf("Optimal doses: %s mg\n", paste(round(self$result$doses, 2), collapse = ", "))) - cat(sprintf("Objective function: %.6f\n", self$result$objf)) - cat(sprintf("Status: %s\n", self$result$status)) - cat(sprintf("Method: %s\n", self$result$method)) + cat(sprintf("Optimal doses: [%.2f, %.2f] mg\n", self$get_doses()[1], self$get_doses()[2])) + cat(sprintf("Objective function: %.10f\n", self$get_objf())) + cat(sprintf("ln(Objective): %.4f\n", log(self$get_objf()))) + cat(sprintf("Method: %s\n", self$get_method())) + cat(sprintf("Status: %s\n", self$get_status())) cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions))) if (!is.null(self$result$auc_predictions)) { cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions))) From bab4eb2cd477e1800a056037436884abf9314305 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Wed, 5 Nov 2025 19:28:31 +0000 Subject: [PATCH 4/5] missing files --- .gitignore | 3 - Examples/Rscript/examples.R | 778 +++++++++++++++++++ Examples/Runs/.gitkeep | 0 Examples/src/bad.csv | 1 + Examples/src/ex.csv | 260 +++++++ Examples/src/ex_full.csv | 261 +++++++ Examples/src/ptaex1.csv | 10 + Examples/src/simTemp.csv | 10 + inst/Examples/Rscript/bestdose_simple_test.R | 37 + inst/Examples/Runs/bestdose_result.rds | Bin 0 -> 34108 bytes inst/Examples/src/bestdose_past.csv | 9 + inst/Examples/src/bestdose_prior.csv | 47 ++ inst/Examples/src/bestdose_target.csv | 9 + 13 files changed, 1422 insertions(+), 3 deletions(-) create mode 100644 Examples/Rscript/examples.R create mode 100644 Examples/Runs/.gitkeep create mode 100644 Examples/src/bad.csv create mode 100644 Examples/src/ex.csv create mode 100644 Examples/src/ex_full.csv create mode 100644 Examples/src/ptaex1.csv create mode 100644 Examples/src/simTemp.csv create mode 100644 inst/Examples/Rscript/bestdose_simple_test.R create mode 100644 inst/Examples/Runs/bestdose_result.rds create mode 100644 inst/Examples/src/bestdose_past.csv create mode 100644 inst/Examples/src/bestdose_prior.csv create mode 100644 inst/Examples/src/bestdose_target.csv diff --git a/.gitignore b/.gitignore index 6cd3c721e..d6cce0004 100755 --- a/.gitignore +++ b/.gitignore @@ -64,7 +64,6 @@ tests/testthat/Runs/ .vscode/ .Rd2pdf65996/ other/ -Examples Experimental .httr-oauth inst/rust/template/* @@ -73,8 +72,6 @@ inst/options/PMoptions.json inst/doc target/ -*_test.R -*.csv *.txt 1 2 diff --git a/Examples/Rscript/examples.R b/Examples/Rscript/examples.R new file mode 100644 index 000000000..15ccf6921 --- /dev/null +++ b/Examples/Rscript/examples.R @@ -0,0 +1,778 @@ +# INTRODUCTION ------------------------------------------------------------ + +# Lines that start with "#" are comments and ignored by R. Follow the +# directions in them. Execute each non-comment line in this script by +# putting your cursor on it and sending it to the R console. +# You can do this in several ways: +# Windows +# R-studio +# 1) The Run button at the top +# 2) Ctrl-Enter +# R GUI - when the script window is active +# 1) The Run line or selection button at the top +# 2) Ctrl-R +# Mac +# R-studio +# 1) The Run button at the top +# 2) Command-Enter +# R GUI +# 1) Command-Enter + +# This script also serves to introduce several R programming functions +# and techniques. For any function, you can get help by typing ?function_name +# in the R console (the lower left window pane in RStudio). + +# Load Pmetrics into memory. You must include this line at the beginning +# of every script. + +library(Pmetrics) + +# EXERCISE 1 - NPAG RUN ------------------------------------------------ + +# EXAMPLE NPAG RUN - tlag, ka, kel, vol + +# It is useful to annotate your runs as above, so that you can remember +# what you did later! + + +# Tell R where your working directory is going to be. +# Windows users: Make sure that you separate directories with a +# forward slash "/" or double backslashes "\\". Unfortunately, Windows is the only OS that uses +# backslashes "\", so R conforms to Unix/Linux style. + +wd <- "C:/Users/siel3/code/LAPKB/Pmetrics_rust/Examples" + +wd <- glue::glue("{getwd()}/inst/Examples/Runs") + +# change to the working directory to the Examples folder +setwd(wd) + +# DATA OBJECT + +# Pmetrics always needs data and a model to run +# create our first data object + +# create a new data object by reading a file +# set the limit of quantification (loq) to 1: see ?PM_data for help +exData <- PM_data$new(data = "../src/ex.csv", loq = 1) + +# you can look at this file directly by opening it in +# a spreadsheet program like Excel, or a text editor + +# exData is an R6 object, which means that contains both data and methods to +# process that data, for example: +exData$data # contains your original datafile +exData$standard_data # contains the standardized and validated data, +exData$summary() # prints the summary of the data to the terminal, or + +# another way to do that is using the more common S3 framework in R: +summary(exData) + +# To look at the contents of an object: +names(exData) + +# other examples of things that can be done with this object are +exData # view the original data in the viewer +exData$print(standard = TRUE) # view the standardized data in the viewer +exData$print(viewer = FALSE) # view original data in console +exData$plot() #plot the raw data; more on that later + +# MODEL OBJECT +# You can specify a model by reading a file or directly as an object. We'll do both. +# The following code creates the same model as in /src/model.txt file. +# See PMmanual() for details on creating models in R compared to text files. +# The advantage of creating them in R is that one does not need to copy model +# files into folders to provide necessary inputs. + +mod1 <- PM_model$new( + pri = list( + Ka = ab(0.1, 0.9), + Ke = ab(0.001, 0.1), + V = ab(30, 120), + lag1 = ab(0, 4) + ), + cov = list( + wt = interp(), + africa = interp("none"), + age = interp(), + gender = interp("none"), + height = interp() + ), + eqn = function(){ + two_comp_bolus + }, + lag = function(){ + lag[1] = lag1 + }, + out = function(){ + Y[1] = X[2]/V + }, + err = list( + proportional(5, c(0.02, 0.05, -0.0002, 0)) + ) +) + + +# look at it +mod1 + +#plot it +mod1$plot() + + +# in the working directory we have another file "model.txt" that contains the old +# representation of the same model we previously presented, let's take a look at it. +system("cat ../src/model.txt") + +# PM_model$new() also accepts the path to a model file +# create the same model using this file +mod1b <- PM_model$new("../src/model.txt") +mod1b + +# PM_model provides a method to update the different elements of a model, for example: +mod1b$update( + pri = list( + ka = ab(0.001, 5) +)) + +# It is case sensitive, so ka is different from Ka. To remove a parameter, set it to NULL. + + +mod1b$arg_list$pri + +# to copy a model use the $clone() method. +mod1b <- mod1$clone() + +# simply using mod1b <- mod1 will cause mod1b to be changed if mod1 is changed, +# as R6 objects use reference semantics. For more details you can refer to +# https://adv-r.hadley.nz/r6.html, Section 14.4. + +#lastly, use the app! PMmanual() and the article on models for details on this. +build_model() #start from scratch +build_model(exData) #start with data to match covariates +build_model(mod1) #start with a model and update it + + + +# To keep everything tidy, we are working in a folder specific to store the runs + + +run1 <- mod1$fit(data = exData, run = 1, overwrite = TRUE) # execute the fit and return the results to run1 + + +# +# After the run is complete the results are returned to the assigned object, here 'run1' + +# you need get the extracted information back into R. +# They will be sequentially numbered as /1, /2, /3,... in your working directory. + +# One benefit of having this fit object is that it is possible to run multiple +# fittings without needing to move data files around + +getwd() +list.files() + +# Result Object - the result is already returned to run1 above, but if you need to load it later, +# you can use PM_load() +run1 <- PM_load(1) + +# Create a PM_result object by reading a run folder. The "1" in the parentheses tells Pmetrics to +# look in the /1 folder. + +# Plot the raw data using R6 with various options. Type ?plot.PM_data in the R console for help. +run1$data$plot() +run1$data$plot(overlay = FALSE, xlim = c(119, 145)) +run1$data$plot(xlim = c(119, 146), group = "gender", group_names = c("Male", "Female"), + marker = list(color = c("blue","red"), symbol = c("circle","triangle-up"))) + +run1$data$plot(xlim = c(119, 146), group = "gender", group_names = c("Male", "Female"), marker = list(color = "Set2")) + + +# The following are the older S3 method with plot(...) for the first two examples +# You can use R6 or S3 for any Pmetrics object +# We will focus on R6 as the more modern way. +plot(run1$data) + +# here's a summary of the original data file; ?summary.PMmatrix for help +run1$data$summary() + +# Plot some observed vs. predicted data. Type ?plot.PMop in the R console for help. +# Click on points to highlight all points from the same subject. +run1$op$plot() +run1$op$plot(pred.type = "pop") +run1$op$plot(line = list(lm = list(ci = 0, color = "red"), loess = FALSE)) + +# the original op object data can be accessed via +run1$op$data + +# get a summary with bias and imprecision of the population predictions; +# ?summary.PMop for help +run1$op$summary(pred.type = "pop") + +# the S3 way +summary(run1$op, pred.type = "pop") + +# look at the summary for the posterior predictions (default pred.type) based +# on means of parameter values +run1$op$summary(icen = "mean") + +# The OP plot can be disaggregated into a Tidy compatible format from $data (see https://www.tidyverse.org/) +# This allow pre processing in ways more flexible than the default plot method. +library(tidyverse) +run1$op$data %>% plot() +run1$op$data %>% + filter(pred > 5) %>% + filter(pred < 10) %>% + summary() + +# see a header with the first 10 rows of the op object +head(run1$op$data, 10) + + +# Plot final population joint density information. Type ?plot.PMfinal in the R console for help. +run1$final$plot() + +# add a kernel density curve +run1$final$plot(line = list(density = list(color = "red"))) +run1$final$data %>% plot() + + +# A bivariate plot. Plotting formulae in R are of the form 'y~x' +run1$final$plot(ke ~ v, + marker = list(color = "red", symbol = "diamond")) + + + +# or the S3 way +plot(run1$final) + +# The original final object can be accessed using +run1$final$data +names(run1$final$data) + +# see the population points +run1$final$popPoints + +# or +run1$final$data$popPoints + +# see the population mean parameter values +run1$final$popMean + +# see a summary with confidence intervals around the medians +# and the Median Absolute Weighted Difference (MAWD); +# ?summary.PMfinal for help +run1$final$summary() + +# summarize the cycle information; ?summary.PM_cycle for help +run1$cycle$summary() +run1$cycle$data %>% summary() + + +# Plot cycle information +# Type ?plot.PM_cycle in the R console for help. +run1$cycle$plot() + +# names of the cycle object; ?makeCycle for help +names(run1$cycle$data) + +# gamma/lamda value on last 6 cycles +tail(run1$cycle$data$gamlam) + +# Plot covariate information. Type ?plot.PMcov in the R console for help. +# Recall that plotting formulae in R are of the form 'y~x' +run1$cov$plot(v ~ wt) + +run1$cov$data %>% + filter(age > 25) %>% + plot(v ~ wt) + +# Plot +run1$cov$plot(ke ~ age, line = list(loess = FALSE, lm = TRUE), + marker = list(symbol = 3)) + +# Another plot with mean Bayesian posterior parameter and covariate values... +# Remember the 'icen' argument? +run1$cov$plot(v ~ wt, icen = "mean") + +# When time is the x variable, the y variable is aggregated by subject. +# In R plot formulae, calculations on the fly can be included using the I() function +run1$cov$plot(I(v * wt) ~ time) + +# The previous cov object can be seen via: +run1$cov + +# but to access individual elements, use: +run1$cov$data[, 1:3] # for example +names(run1$cov) + +# summarize with mean covariates; ?summary.PMcov for help +run1$cov$summary(icen = "mean") + + +# Look at all possible covariate-parameter relationships by multiple linear regression with forward +# and backward elimination - type ?PM_step in the R console for help. +run1$step() +# or on the cov object directly +run1$cov$step() +# icen works here too.... +run1$step(icen = "mean") +# forward elimination only +run1$step(direction = "forward") + + +# EXERCISE 2 - NPAG WITH COVARIATES --------------------------------------- + +# Again, without copying files, let's create another run object, this time using +# a model that include covariates + +# First clone mod1 +mod2 <- mod1$clone() + +# Then update it +mod2 <- mod2$update( + pri = list( + V0 = ab(30, 120), + V = NULL + ), + sec = function(x) { + V = V0*(WT/55) + }, + err = list( + proportional(2.39, c(0.02, 0.05, -0.0002, 0), fixed = TRUE) + ) +) +# we can also make a model object by loading a file +mod2b <- PM_model$new("../src/model2.txt") + + +run2 <- mod2$fit(data = exData, run = 2, overwrite = TRUE) + +# for future loading +run2 <- PM_load(2) + + + + +# EXERCISE 3 - COMPARING MODELS ------------------------------------------- + + +# Let's compare model 1 and model 2. You can compare any number of models. +# Type ?PM_compare for help. +PM_compare(run1, run2) + + + +# EXERCISE 4 - MODEL VALIDATION ------------------------------------------- + +# MODEL VALIDATION EXAMPLES +# Example of Pmetrics visual predictive check and prediction-corrected visual predictive check +# for model validation - be sure to have executed the NPAG run above +# Type ?makeValid in the R console for help. +# Choose wt as the covariate to bin. Accept all default bin sizes. +run2$validate(limits = c(0, 3)) + +# To see what it contains, use: +run2$valid + +# Default visual predictive check; ?plot.PM_valid for help +run2$valid$plot() + +# or old S3 +plot(run2$valid) + + +# Generate a prediction-corrected visual predictive check; type ?plot.PMvalid in the R console for help. +run2$valid$plot(type = "pcvpc") + +# Create an npde plot +run2$valid$plot(type = "npde") + +# Here is another way to generate a visual predicive check... +npc_2 <- run2$valid$simdata$plot(obs = run2$op, log = FALSE, binSize = 0.5) + +# The jagged appearance of the plot when binSize=0 is because different subjects have +# different doses, covariates, and observation times, which are all combined in one simulation. +# Collapsing simulation times within 1 hour bins (binSize=1) smooths +# the plot, but can change the P-values in the numerical predictive check below. + +npc_2 +# ...and here is a numerical predictive check +# P-values are binomial test of proportion of observations less than +# the respective quantile + + +# EXERCISE 5 - SIMULATOR RUN ---------------------------------------------- + +setwd(wd) +dir.create("../Sim") +setwd("../Sim") + +# The following will simulate 100 sets of parameters/concentrations using the +# first subject in the data file as a template. +# Limits are put on the simulated parameter ranges to be the same as in the model. +# The population parameter values from the NPAG run in exercise 2 are used for the Monte Carlo Simulation. +simdata <- run2$sim(include = 1, limits = NA, nsim = 100) + +# Below is the alternate way to simulate, which is particularly useful if you define +# your own population parameters. See ?SIMrun for details on this as well as +# the article on simulation linked by PMmanual(). +poppar <- list( + wt = 1, + mean = c(0.6, 0.05, 77.8, 1.2), + cov = diag(c(0.07, 0.0004, 830, 0.45)) +) + +simOther <- PM_sim$new(poppar = poppar, data = exData, model = mod1, + include = 1, limits = NA, nsim = 100) + + +# simulate from a model with new data +sim_new <- run2$sim( + data = "../src/ptaex1.csv", + include = 2, limits = NA, + predInt = c(120, 144, 0.5) +) + + + +# Plot them; ?plot.PM_sim for help +simdata$plot() +simOther$plot() +sim_new$plot(log = FALSE) + +# Simulate using multiple subjects as templates +simdata <- run2$sim(include = 1:4, limits = NA, nsim = 100) + +# Plot the third simulation; use include/exclude to specify the ID numbers, +# which are the same as the ID numbers in the template data file +simdata$plot(include = 2) + +# or in S3 +plot(simdata$data, include = 3) + +# Parse and combine multiple files and plot them. Note that combining simulations from templates +# with different simulated observation times can lead to unpredictable plots +simdata2 <- run2$sim(include = 1:4, limits = NA, nsim = 100) +simdata2$plot() + +# simulate with covariates +# in this case we use the covariate-parameter correlations from run 2, which +# are found in the cov.2 object; we re-define the mean weight to be 50 with +# SD of 20, and limits of 10 to 70 kg. We fix africa, gender and height covariates, +# but allow age (the last covariate) to be simulated, using the mean, sd, and +# limits in the original population, since we didn't specify them. +# See ?SIMrun for more help on this and the Pmetrics manual. + +covariate <- list( + cov = run2$cov, + mean = list(wt = 50), + sd = list(wt = 20), + limits = list(wt = c(10, 70)), + fix = c("africa", "gender", "height") +) + +# now simulate with this covariate list object +simdata3 <- run2$sim(include = 1:4, limits = NA, nsim = 100, covariate = covariate) + +# compare difference in simulations without covariates simulated... +# PM_sim's plot function defaults to the first simulation if there +# are multiple simulations and "at" is not specified. +simdata$plot() + +# ...and with covariates simulated +simdata3$plot() + +# Here are the simulated parameters and covariates for the first subject's +# template; note that both wt and age are simulated, using proper covariances +# with simulated PK parameters +simdata3$data$parValues + +# We can summarize simulations too. See ?summary.PM_sim for help. +simdata3$summary(field = "obs", include = 3) + +# look in the working directory and find the "c_simdata.csv" and "c_simmodel.txt" files +# which were made when you simulated with covariates. Compare to original +# "simdata.csv" and "simmoddel.txt" files to note that simulated covariates become +# Primary block variables, and are removed from the template data file. + +# EXERCISE 6 - SAVING PMETRICS OBJECTS ------------------------------------ + +setwd(wd) + +# The following objects have methods to save them to or load them from files: +# PM_fit +# PM_result +# PM_sim +# PM_pta + +# Example - save the PM_result (run2) to the "2" folder +run2$save(file = "2/outputs/run2.rds") # rds is the recommended file extension +list.files("2/outputs") +copy_run2 <- PM_load(file = "2/outputs/run2.rds") +copy_run2 + +# If you want to save multiple objects into one single file, R provides the +# following functionality + +save(exData, mod1, run1, simdata, file = "2/test_drug.Rdata") +list.files("2") +load("2/test_drug.Rdata") + +# or +save.image("2/workspace.Rdata") # This will save all variables in your environment +list.files("2") +load("2/workspace.Rdata") + + +# EXERCISE 7 - CONTINUING RUNS OR EXTERNAL VALIDATIONS -------------------- + +# Example of a run with a non-uniform density +# This is a good way to continue a previous run, +# in this case it continues where run 1 left off + +# note that we can supply a run number to model, data, and prior arguments. The numbers do not +# have to be the same. This will copy the appropriate files from the specified run to be used +# in the current run. By specifying a prior, we are starting with the non-uniform density from the +# end of the specified fun. +run3 <- mod2$fit(data = exData, prior = 2) +run3 <- PM_load(3) + +# We could also generate Bayesian posterior parameter estimates for a new population this way: +# run3 <- mod2$fit(data=PM_data("newPop.csv"), prior = 2, cycles = 0) +# This won't run because we don't have a newPop.csv file, +# but shows you how it could be done. + + + +# EXERCISE 8 - PROBABILITY OF TARGET ATTAINMENT --------------------------- + +# Note: these can be computationally intense and take some time. + +# Examples of probability of target attainment analysis +# Be sure to have executed the NPAG run above and used PM_load(2) in EXERCISE 2 +# Type ?PM_sim, ?PM_pta, or ?plot.PM_pta into the R console for help. + + +# simulate with the template data file that contains different doses +# Look at PM_sim for help on arguments to this function, including predInt, +# seed, limits, nsim. + +simlist1 <- run1$sim( + limits = c(0, 3), data = "../src/ptaex1.csv", + predInt = c(120, 144, 0.5), seed = rep(-17, 4) +) + +# now simulate with covariates; make sure that you defined the covariate +# object first in Exercise 5 above and have loaded the results of Exercise 2 +# with PM_load(2) +simlist2 <- run2$sim( + limits = 5, data = "../src/ptaex1.csv", + predInt = c(120, 144, 0.5), seed = rep(-17, 4), + covariate = covariate +) + +# make the first PMpta object to calculate the time above each target for at +# least 60% of the dosing +# interval from 120 to 144 hours. Include labels for the simulations. +# ?makePTA for help +# define simulation labels first + +simlabels <- c("600 mg daily", "1200 mg daily", "300 mg bid", "600 mg bid") + +pta1_2 <- PM_pta$new( + simdata = simlist1, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), target_type = "time", + success = 0.6, start = 120, end = 144 +) + +pta1b_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), target_type = "time", + success = 0.6, start = 120, end = 144 +) + +# summarize the results +pta1_2$summary() + +# in the summary()$pta, reg_num is the simulation template ID number; +# target in this case is the MIC; prop_success is the proportion of the simulated +# profiles for each dose/MIC that are above the success threshold (0.6); pdi.mean and pdi.sd +# are the mean and standard deviation of the pharmacodynamic index (PDI), in this case proportion of the interval > MIC. +# In the $pdi, target and simnum are the same, but now the median and confidence +# intervals (default 95%) PDI are shown. +# ?summary.PMpta for help + +# Plot the first without covariates. We didn't include simulation +# labels in the makePTA command, so generics are used here, but we move it to +# the bottom left; ?legend for help on arguments to supply to the +# legend list argument to plot.PMpta. +pta1_2$plot(ylab = "Proportion with %T>MIC of at least 60%", grid = TRUE, legend = list(x = "bottomleft")) + +pta1b_2$summary() + +# Plot the second with covariates simulated. Note the regimen labels are included, but we move +# the legend to the bottom left. +pta1b_2$plot( + ylab = "Proportion with %T>MIC of at least 60%", grid = TRUE, + legend = list(x = "bottomleft") +) + +# Now we'll define success as free auc:mic > 100 with a free drug fraction of 50% +pta2_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), + free_fraction = 0.7, + target_type = "auc", success = 100, start = 120, end = 144 +) +summary(pta2_2) +pta2_2$plot( + ylab = "Proportion with AUC/MIC of at least 100", grid = TRUE, + legend = list(x = "bottomleft") +) + +# success is Cmax/MIC >=10 +pta3_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), + target_type = "peak", success = 10, start = 120, end = 144 +) +pta3_2$summary() +pta3_2$plot(ylab = "Proportion with peak/MIC of at least 10", grid = TRUE) + +# success = Cmin:MIC > 1 +pta4_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), + target_type = "min", success = 1, start = 120, end = 144 +) +pta4_2$summary() +pta4_2$plot(ylab = "Proportion with Cmin/MIC of at least 1", grid = TRUE, legend = list(x = "bottomleft")) + +# now plot the PDI (pharmacodynamic index) of each regimen, rather than the proportion +# of successful profiles. A PDI plot is always available for PMpta objects. +pta4_2$plot(at = 1, type = "pdi", ylab = "Cmin:MIC", grid = TRUE) + +# Each regimen has the 90% confidence interval PDI around the median curve, +# in the corresponding, semi-transparent color. Make the CI much narrower... +pta4_2$plot(at = 1, type = "pdi", ci = 0.1) + +# ...or gone altogether, put back the grid, redefine the colors, and make lines narrower +pta4_2$plot( + + at = 1, type = "pdi", ci = 0, grid = TRUE, + line = list( + color = c("blue", "purple", "black", "brown"), + width = 1 + ) +) + +# now let's repeat the analysis but simulate the distribution of MICs +# using susceptibility of Staphylococcus aureus to vancomycin contained +# in the mic1 dataset within Pmetrics (?mic1) + +# see the source with ?mic1 +pta4b_2 <- PM_pta$new( + simdata = simlist2, + simlabels = c("600 mg daily", "1200 mg daily", "300 mg bid", "600 mg bid"), + target = makePTAtarget(mic1), target_type = "min", success = 1, start = 120, end = 144 +) + +pta4b_2$summary() +# plot it +pta4b_2$plot( + grid = TRUE, ylab = "Proportion with Cmin/MIC of at least 1", + marker = list(color = "red"), line = list(color = "black") +) +pta4b_2$plot(type = "pdi", grid = TRUE, ylab = "Proportion with Cmin/MIC of at least 1") + +# note that the plot changes since target MICs are no longer discrete +# since most of the MICs are very low, the regimens all look very similar + +# success = concentration at time 3 hours:MIC > 2 +pta5_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), target_type = 123, success = 2, start = 120, end = 144 +) +pta5_2$summary() +pta5_2$plot(ylab = "Proportion with C3/MIC of at least 1", grid = TRUE, legend = list(x = .3, y = 0.1)) + + +# success is trough >10 +pta6_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = 10, target_type = 144, success = 1, start = 120, end = 144 +) +plot(pta6_2) +pta6_2$summary() + +# EXERCISE 10 - OPTIMAL SAMPLE TIMES -------------------------------------- + +setwd(wd) +dir.create("../MMopt") +setwd("../MMopt") + +# calculate MM-optimal sample times for Run 2, and the 1200 mg once daily dose in the PTA +# By specifying the predInt to start and stop at 120 and 144 hours, with an interval of 1 hour, +# we are sampling at steady state. Including "subject 2", means only the 1200 mg once daily dose +# will serve as a simulation template. + +run2$opt( + data = "../src/ptaex1.csv", + nsamp = 2, predInt = c(120, 140, 1), + limits = NA, + include = 2 +) +# see the optimal sample times and the Bayes Risk of misclassification, +# which is only useful to compare optimal sampling regimens, i.e. the +# absolute value is less helpful, but is the statistic minimized by the +# selected optimal sample times for a given model + +mmopt_2 <- PM_opt$new( + poppar = run2$final, + model = run2$model, + data = "../src/ptaex1.csv", + nsamp = 2, predInt = c(120, 140, 1), + include = 2 +) + +# plot it, with the red lines indicating the optimal sample times. +# see ?plot.MMopt for help +mmopt_2$plot() +plot(mmopt_2) +plot(mmopt_2, line = list(color = "slategrey"), times = list(color = "orange")) + +# EXERCISE 11 - ASSAY ERROR ----------------------------------------------- +# see ?makeErrorPoly for more help +# This will let you choose the best set of C0, C1, C2, C3 for your modeling, +# based on assay validation data which includes the "obs", which are the +# nominal concentrations of the standards, and "sd", which is the standard +# deviation of replicate measurements of each of the standards, i.e. the +# inter-day and/or intra-day standard deviation + +obs <- c(0, 25, 50, 100, 250, 500, 1000, 2000, 5000) +sd <- c(0.5, 6.4, 8.6, 12, 8.6, 37.2, 60.1, 165.7, 483) + +# See plots.pdf, page 50 +makeErrorPoly(obs = obs, sd = sd) + +# choose the one with the best R-squared that will never result in a +# negative value for the SD + + + +# Ancillary functions ----------------------------------------------------- + +# Be sure to check out the help files for the following functions: +# +# makeAUC() - calculate AUC from a variety of inputs +# makeNCA() - non-compartmental analysis +# NM2PM() - convert NONMEM data files to Pmetrics data files +# qgrowth() - CDC growth charts +# zBMI() - CDC Pediatric BMI z-scores and percentiles +# ss.PK() - sample size for Phase 1 PK studies diff --git a/Examples/Runs/.gitkeep b/Examples/Runs/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/Examples/src/bad.csv b/Examples/src/bad.csv new file mode 100644 index 000000000..0d644a159 --- /dev/null +++ b/Examples/src/bad.csv @@ -0,0 +1 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, #ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT 1,1,0,1,600,.,.,1,.,..,.,.,.,.,46.7,1,21,1,160 1,.,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,0,120,.,.,.,.,.,10.44,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,0,121,.,.,.,.,.,12.89,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,122,.,.,.,.,.,14.98,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,125.99,.,.,.,.,.,16.69,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,129,.,.,.,.,.,20.15,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,132,.,.,.,.,.,14.97,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,143.98,.,.,.,.,.,12.57,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 2,1,0,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,24,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,48,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,72,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,96,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,0,120,.,.,.,.,.,3.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,1,120,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,0,120.98,.,.,.,.,.,5.84,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,121.98,.,.,.,.,.,6.54,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,126,.,.,.,.,.,6.14,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,129.02,.,.,.,.,.,6.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,132.02,.,.,.,.,.,4.44,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,144,.,.,.,.,.,3.76,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 3,1,0,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,0,120.08,.,.,.,.,.,4.06,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,121.07,.,.,.,.,.,3.24,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,122.08,.,.,.,.,.,3.09,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,126.08,.,.,.,.,.,7.98,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,129.05,.,.,.,.,.,7.23,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,132.1,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,144.08,.,.,.,.,.,3.82,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 4,1,0,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,24,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,48,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,72,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,96,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,0,120,.,.,.,.,.,2.1,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,1,120,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,0,121,.,.,.,.,.,3.05,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,122.02,.,.,.,.,.,5.21,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,126,.,.,.,.,.,5.09,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,129.03,.,.,.,.,.,4.24,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,132,.,.,.,.,.,3.69,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,144.02,.,.,.,.,.,1.96,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 5,1,0,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,24,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,48,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,72,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,96,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,1,120,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,0,121,.,.,.,.,.,2.64,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,122,.,.,.,.,.,4.8,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,126,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,129.02,.,.,.,.,.,4.13,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,132,.,.,.,.,.,2.81,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,144,.,.,.,.,.,2.21,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 6,1,0,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,24,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,48,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,72,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,96,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,0,120,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,1,120,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,0,121,.,.,.,.,.,6.89,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,121.98,.,.,.,.,.,6.64,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,126,.,.,.,.,.,13.72,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,129,.,.,.,.,.,12.69,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,131.98,.,.,.,.,.,10.58,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,144.98,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 7,1,0,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,24,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,48,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,72,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,96,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,0,120,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,1,120,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,0,121.03,.,.,.,.,.,4.46,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,122.03,.,.,.,.,.,4.54,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,126.02,.,.,.,.,.,12.19,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,129.08,.,.,.,.,.,12.1,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,132.03,.,.,.,.,.,8.61,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,144.03,.,.,.,.,.,6.37,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 8,1,0,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,24,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,48,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,72,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,96,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,0,120,.,.,.,.,.,6.19,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,1,120,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,0,121.03,.,.,.,.,.,6.33,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,122,.,.,.,.,.,6.24,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,125.98,.,.,.,.,.,13.03,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,128.98,.,.,.,.,.,11.86,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,132,.,.,.,.,.,11.45,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,143.98,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 9,1,0,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,24,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,48,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,72,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,96,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,0,120,.,.,.,.,.,2.85,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,1,120,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,0,120.97,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,122,.,.,.,.,.,6.65,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,125.98,.,.,.,.,.,6.81,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,128.98,.,.,.,.,.,6.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,132,.,.,.,.,.,7.48,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,143.98,.,.,.,.,.,4.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 10,1,0,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,24,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,48,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,72,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,96,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,1,120,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,0,121,.,.,.,.,.,4.36,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,122.02,.,.,.,.,.,7.79,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,126,.,.,.,.,.,11.02,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,129,.,.,.,.,.,8.86,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,131.97,.,.,.,.,.,6.09,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,144,.,.,.,.,.,4.15,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 11,1,0,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,24,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,48,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,72,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,96,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,0,120,.,.,.,.,.,2.09,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,1,120,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,0,121.03,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,122,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,125.98,.,.,.,.,.,7.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,129,.,.,.,.,.,6.31,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,132,.,.,.,.,.,5.82,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,144.13,.,.,.,.,.,2.63,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 12,1,0,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,24,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,48,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,72,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,96,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,0,120,.,.,.,.,.,7.09,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,1,120,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,0,121.03,.,.,.,.,.,6.18,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,122.13,.,.,.,.,.,8.66,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,126,.,.,.,.,.,11.16,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,129,.,.,.,.,.,9.51,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,132,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,144,.,.,.,.,.,7.89,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 13,1,0,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,24,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,48,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,72,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,96,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,0,120,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,1,120,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,0,121,.,.,.,.,.,3.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,122,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,126,.,.,.,.,.,10.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,129.02,.,.,.,.,.,12.84,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,132,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,144,.,.,.,.,.,8.06,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 14,1,0,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,24,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,48,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,72,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,96,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,0,120,.,.,.,.,.,3.63,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,1,120,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,0,121,.,.,.,.,.,4.49,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,122,.,.,.,.,.,5.5,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,126,.,.,.,.,.,7.28,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,129,.,.,.,.,.,5.27,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,132,.,.,.,.,.,4.89,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,144,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 15,1,0,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,24,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,48,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,72,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,96,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,0,120,.,.,.,.,.,5.53,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,1,120,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,0,121,.,.,.,.,.,4.81,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,122,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,126,.,.,.,.,.,9.96,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,129,.,.,.,.,.,8.55,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,132.05,.,.,.,.,.,7.54,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,144.05,.,.,.,.,.,5.74,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 16,1,0,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,24,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,48,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,72,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,96,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,0,120,.,.,.,.,.,5.48,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,1,120,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,0,121,.,.,.,.,.,6.59,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,122,.,.,.,.,.,8.91,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,126,.,.,.,.,.,10.57,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,129,.,.,.,.,.,9.52,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,132,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,143.97,.,.,.,.,.,4.96,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 17,1,0,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,24,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,48,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,72,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,96,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,0,120,.,.,.,.,.,2.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,1,120,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,0,121.02,.,.,.,.,.,1.86,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,122.02,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,126,.,.,.,.,.,9.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,129,.,.,.,.,.,6.96,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,132,.,.,.,.,.,5.64,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,144.08,.,.,.,.,.,3.59,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 18,1,0,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,24,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,48,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,72,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,96,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,0,120,.,.,.,.,.,7.95,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,1,120,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,0,120.98,.,.,.,.,.,7.47,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,121.98,.,.,.,.,.,8.67,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,126,.,.,.,.,.,13.83,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,129.17,.,.,.,.,.,14.01,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,132.17,.,.,.,.,.,8.97,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,143.97,.,.,.,.,.,8.4,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 19,1,0,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,24,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,48,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,72,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,96,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,0,120,.,.,.,.,.,5.42,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,1,120,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,0,121,.,.,.,.,.,7.08,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,122,.,.,.,.,.,7.27,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,125.98,.,.,.,.,.,20.07,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,128.98,.,.,.,.,.,18.24,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,132,.,.,.,.,.,15.36,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,144,.,.,.,.,.,10.92,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 20,1,0,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,24,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,48,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,72,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,96,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,0,120,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,1,120,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,0,120.77,.,.,.,.,.,4.5,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,121.75,.,.,.,.,.,3.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,125.67,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,128.67,.,.,.,.,.,11.56,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,143.67,.,.,.,.,.,6.45,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 \ No newline at end of file diff --git a/Examples/src/ex.csv b/Examples/src/ex.csv new file mode 100644 index 000000000..392a45499 --- /dev/null +++ b/Examples/src/ex.csv @@ -0,0 +1,260 @@ +#ID,TIME,DOSE,OUT,WT,AFRICA,AGE,GENDER,HEIGHT +1,0,600,.,46.7,1,21,1,160 +1,24,600,.,46.7,1,21,1,160 +1,48,600,.,46.7,1,21,1,160 +1,72,600,.,46.7,1,21,1,160 +1,96,600,.,46.7,1,21,1,160 +1,120,.,10.44,46.7,1,21,1,160 +1,120,600,.,46.7,1,21,1,160 +1,121,.,12.89,46.7,1,21,1,160 +1,122,.,14.98,46.7,1,21,1,160 +1,125.99,.,16.69,46.7,1,21,1,160 +1,129,.,20.15,46.7,1,21,1,160 +1,132,.,14.97,46.7,1,21,1,160 +1,143.98,.,12.57,46.7,1,21,1,160 +2,0,600,.,66.5,1,30,1,174 +2,24,600,.,66.5,1,30,1,174 +2,48,600,.,66.5,1,30,1,174 +2,72,600,.,66.5,1,30,1,174 +2,96,600,.,66.5,1,30,1,174 +2,120,.,3.56,66.5,1,30,1,174 +2,120,600,.,66.5,1,30,1,174 +2,120.98,.,5.84,66.5,1,30,1,174 +2,121.98,.,6.54,66.5,1,30,1,174 +2,126,.,6.14,66.5,1,30,1,174 +2,129.02,.,6.56,66.5,1,30,1,174 +2,132.02,.,4.44,66.5,1,30,1,174 +2,144,.,3.76,66.5,1,30,1,174 +3,0,600,.,46.7,1,24,0,164 +3,24,600,.,46.7,1,24,0,164 +3,48,600,.,46.7,1,24,0,164 +3,72,600,.,46.7,1,24,0,164 +3,96,600,.,46.7,1,24,0,164 +3,120,600,.,46.7,1,24,0,164 +3,120.08,.,4.06,46.7,1,24,0,164 +3,121.07,.,3.24,46.7,1,24,0,164 +3,122.08,.,3.09,46.7,1,24,0,164 +3,126.08,.,7.98,46.7,1,24,0,164 +3,129.05,.,7.23,46.7,1,24,0,164 +3,132.1,.,4.71,46.7,1,24,0,164 +3,144.08,.,3.82,46.7,1,24,0,164 +4,0,600,.,50.8,1,25,1,165 +4,24,600,.,50.8,1,25,1,165 +4,48,600,.,50.8,1,25,1,165 +4,72,600,.,50.8,1,25,1,165 +4,96,600,.,50.8,1,25,1,165 +4,120,.,2.1,50.8,1,25,1,165 +4,120,600,.,50.8,1,25,1,165 +4,121,.,3.05,50.8,1,25,1,165 +4,122.02,.,5.21,50.8,1,25,1,165 +4,126,.,5.09,50.8,1,25,1,165 +4,129.03,.,4.24,50.8,1,25,1,165 +4,132,.,3.69,50.8,1,25,1,165 +4,144.02,.,1.96,50.8,1,25,1,165 +5,0,600,.,65.8,1,22,1,181 +5,24,600,.,65.8,1,22,1,181 +5,48,600,.,65.8,1,22,1,181 +5,72,600,.,65.8,1,22,1,181 +5,96,600,.,65.8,1,22,1,181 +5,120,.,2.93,65.8,1,22,1,181 +5,120,600,.,65.8,1,22,1,181 +5,121,.,2.64,65.8,1,22,1,181 +5,122,.,4.8,65.8,1,22,1,181 +5,126,.,3.7,65.8,1,22,1,181 +5,129.02,.,4.13,65.8,1,22,1,181 +5,132,.,2.81,65.8,1,22,1,181 +5,144,.,2.21,65.8,1,22,1,181 +6,0,600,.,65,1,23,1,177 +6,24,600,.,65,1,23,1,177 +6,48,600,.,65,1,23,1,177 +6,72,600,.,65,1,23,1,177 +6,96,600,.,65,1,23,1,177 +6,120,.,6.92,65,1,23,1,177 +6,120,600,.,65,1,23,1,177 +6,121,.,6.89,65,1,23,1,177 +6,121.98,.,6.64,65,1,23,1,177 +6,126,.,13.72,65,1,23,1,177 +6,129,.,12.69,65,1,23,1,177 +6,131.98,.,10.58,65,1,23,1,177 +6,144.98,.,6.62,65,1,23,1,177 +7,0,600,.,51.7,1,27,0,161 +7,24,600,.,51.7,1,27,0,161 +7,48,600,.,51.7,1,27,0,161 +7,72,600,.,51.7,1,27,0,161 +7,96,600,.,51.7,1,27,0,161 +7,120,.,5.41,51.7,1,27,0,161 +7,120,600,.,51.7,1,27,0,161 +7,121.03,.,4.46,51.7,1,27,0,161 +7,122.03,.,4.54,51.7,1,27,0,161 +7,126.02,.,12.19,51.7,1,27,0,161 +7,129.08,.,12.1,51.7,1,27,0,161 +7,132.03,.,8.61,51.7,1,27,0,161 +7,144.03,.,6.37,51.7,1,27,0,161 +8,0,600,.,51.2,1,22,1,163 +8,24,600,.,51.2,1,22,1,163 +8,48,600,.,51.2,1,22,1,163 +8,72,600,.,51.2,1,22,1,163 +8,96,600,.,51.2,1,22,1,163 +8,120,.,6.19,51.2,1,22,1,163 +8,120,600,.,51.2,1,22,1,163 +8,121.03,.,6.33,51.2,1,22,1,163 +8,122,.,6.24,51.2,1,22,1,163 +8,125.98,.,13.03,51.2,1,22,1,163 +8,128.98,.,11.86,51.2,1,22,1,163 +8,132,.,11.45,51.2,1,22,1,163 +8,143.98,.,7.83,51.2,1,22,1,163 +9,0,600,.,55,1,23,1,174 +9,24,600,.,55,1,23,1,174 +9,48,600,.,55,1,23,1,174 +9,72,600,.,55,1,23,1,174 +9,96,600,.,55,1,23,1,174 +9,120,.,2.85,55,1,23,1,174 +9,120,600,.,55,1,23,1,174 +9,120.97,.,3.7,55,1,23,1,174 +9,122,.,6.65,55,1,23,1,174 +9,125.98,.,6.81,55,1,23,1,174 +9,128.98,.,6.51,55,1,23,1,174 +9,132,.,7.48,55,1,23,1,174 +9,143.98,.,4.51,55,1,23,1,174 +10,0,600,.,52.1,1,32,1,163 +10,24,600,.,52.1,1,32,1,163 +10,48,600,.,52.1,1,32,1,163 +10,72,600,.,52.1,1,32,1,163 +10,96,600,.,52.1,1,32,1,163 +10,120,.,2.93,52.1,1,32,1,163 +10,120,600,.,52.1,1,32,1,163 +10,121,.,4.36,52.1,1,32,1,163 +10,122.02,.,7.79,52.1,1,32,1,163 +10,126,.,11.02,52.1,1,32,1,163 +10,129,.,8.86,52.1,1,32,1,163 +10,131.97,.,6.09,52.1,1,32,1,163 +10,144,.,4.15,52.1,1,32,1,163 +11,0,600,.,56.5,1,34,1,165 +11,24,600,.,56.5,1,34,1,165 +11,48,600,.,56.5,1,34,1,165 +11,72,600,.,56.5,1,34,1,165 +11,96,600,.,56.5,1,34,1,165 +11,120,.,2.09,56.5,1,34,1,165 +11,120,600,.,56.5,1,34,1,165 +11,121.03,.,2.68,56.5,1,34,1,165 +11,122,.,4.71,56.5,1,34,1,165 +11,125.98,.,7.71,56.5,1,34,1,165 +11,129,.,6.31,56.5,1,34,1,165 +11,132,.,5.82,56.5,1,34,1,165 +11,144.13,.,2.63,56.5,1,34,1,165 +12,0,600,.,47.9,1,54,0,160 +12,24,600,.,47.9,1,54,0,160 +12,48,600,.,47.9,1,54,0,160 +12,72,600,.,47.9,1,54,0,160 +12,96,600,.,47.9,1,54,0,160 +12,120,.,7.09,47.9,1,54,0,160 +12,120,600,.,47.9,1,54,0,160 +12,121.03,.,6.18,47.9,1,54,0,160 +12,122.13,.,8.66,47.9,1,54,0,160 +12,126,.,11.16,47.9,1,54,0,160 +12,129,.,9.51,47.9,1,54,0,160 +12,132,.,8.14,47.9,1,54,0,160 +12,144,.,7.89,47.9,1,54,0,160 +13,0,600,.,60.5,1,24,1,180 +13,24,600,.,60.5,1,24,1,180 +13,48,600,.,60.5,1,24,1,180 +13,72,600,.,60.5,1,24,1,180 +13,96,600,.,60.5,1,24,1,180 +13,120,.,6.62,60.5,1,24,1,180 +13,120,600,.,60.5,1,24,1,180 +13,121,.,3.18,60.5,1,24,1,180 +13,122,.,5.41,60.5,1,24,1,180 +13,126,.,10.18,60.5,1,24,1,180 +13,129.02,.,12.84,60.5,1,24,1,180 +13,132,.,12.35,60.5,1,24,1,180 +13,144,.,8.06,60.5,1,24,1,180 +14,0,600,.,59.2,1,26,1,174 +14,24,600,.,59.2,1,26,1,174 +14,48,600,.,59.2,1,26,1,174 +14,72,600,.,59.2,1,26,1,174 +14,96,600,.,59.2,1,26,1,174 +14,120,.,3.63,59.2,1,26,1,174 +14,120,600,.,59.2,1,26,1,174 +14,121,.,4.49,59.2,1,26,1,174 +14,122,.,5.5,59.2,1,26,1,174 +14,126,.,7.28,59.2,1,26,1,174 +14,129,.,5.27,59.2,1,26,1,174 +14,132,.,4.89,59.2,1,26,1,174 +14,144,.,2.68,59.2,1,26,1,174 +15,0,450,.,43,1,19,0,150 +15,24,450,.,43,1,19,0,150 +15,48,450,.,43,1,19,0,150 +15,72,450,.,43,1,19,0,150 +15,96,450,.,43,1,19,0,150 +15,120,.,5.53,43,1,19,0,150 +15,120,450,.,43,1,19,0,150 +15,121,.,4.81,43,1,19,0,150 +15,122,.,8.14,43,1,19,0,150 +15,126,.,9.96,43,1,19,0,150 +15,129,.,8.55,43,1,19,0,150 +15,132.05,.,7.54,43,1,19,0,150 +15,144.05,.,5.74,43,1,19,0,150 +16,0,600,.,64.4,1,25,1,173 +16,24,600,.,64.4,1,25,1,173 +16,48,600,.,64.4,1,25,1,173 +16,72,600,.,64.4,1,25,1,173 +16,96,600,.,64.4,1,25,1,173 +16,120,.,5.48,64.4,1,25,1,173 +16,120,600,.,64.4,1,25,1,173 +16,121,.,6.59,64.4,1,25,1,173 +16,122,.,8.91,64.4,1,25,1,173 +16,126,.,10.57,64.4,1,25,1,173 +16,129,.,9.52,64.4,1,25,1,173 +16,132,.,7.83,64.4,1,25,1,173 +16,143.97,.,4.96,64.4,1,25,1,173 +17,0,600,.,54.8,1,23,1,170 +17,24,600,.,54.8,1,23,1,170 +17,48,600,.,54.8,1,23,1,170 +17,72,600,.,54.8,1,23,1,170 +17,96,600,.,54.8,1,23,1,170 +17,120,.,2.11,54.8,1,23,1,170 +17,120,600,.,54.8,1,23,1,170 +17,121.02,.,1.86,54.8,1,23,1,170 +17,122.02,.,6.92,54.8,1,23,1,170 +17,126,.,9.11,54.8,1,23,1,170 +17,129,.,6.96,54.8,1,23,1,170 +17,132,.,5.64,54.8,1,23,1,170 +17,144.08,.,3.59,54.8,1,23,1,170 +18,0,450,.,44.3,1,20,0,164 +18,24,450,.,44.3,1,20,0,164 +18,48,450,.,44.3,1,20,0,164 +18,72,450,.,44.3,1,20,0,164 +18,96,450,.,44.3,1,20,0,164 +18,120,.,7.95,44.3,1,20,0,164 +18,120,450,.,44.3,1,20,0,164 +18,120.98,.,7.47,44.3,1,20,0,164 +18,121.98,.,8.67,44.3,1,20,0,164 +18,126,.,13.83,44.3,1,20,0,164 +18,129.17,.,14.01,44.3,1,20,0,164 +18,132.17,.,8.97,44.3,1,20,0,164 +18,143.97,.,8.4,44.3,1,20,0,164 +19,0,600,.,50,1,36,1,168 +19,24,600,.,50,1,36,1,168 +19,48,600,.,50,1,36,1,168 +19,72,600,.,50,1,36,1,168 +19,96,600,.,50,1,36,1,168 +19,120,.,5.42,50,1,36,1,168 +19,120,600,.,50,1,36,1,168 +19,121,.,7.08,50,1,36,1,168 +19,122,.,7.27,50,1,36,1,168 +19,125.98,.,20.07,50,1,36,1,168 +19,128.98,.,18.24,50,1,36,1,168 +19,132,.,15.36,50,1,36,1,168 +19,144,.,10.92,50,1,36,1,168 +20,0,600,.,59,1,31,1,170 +20,24,600,.,59,1,31,1,170 +20,48,600,.,59,1,31,1,170 +20,72,600,.,59,1,31,1,170 +20,96,600,.,59,1,31,1,170 +20,120,.,4.71,59,1,31,1,170 +20,120,600,.,59,1,31,1,170 +20,120.77,.,4.5,59,1,31,1,170 +20,121.75,.,3.35,59,1,31,1,170 +20,125.67,.,12.35,59,1,31,1,170 +20,128.67,.,11.56,59,1,31,1,170 +20,143.67,.,6.45,59,1,31,1,170 \ No newline at end of file diff --git a/Examples/src/ex_full.csv b/Examples/src/ex_full.csv new file mode 100644 index 000000000..0f3b8679c --- /dev/null +++ b/Examples/src/ex_full.csv @@ -0,0 +1,261 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, +#ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT +1,1,0,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,120,.,.,.,.,.,10.44,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,121,.,.,.,.,.,12.89,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,122,.,.,.,.,.,14.98,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,125.99,.,.,.,.,.,16.69,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,129,.,.,.,.,.,20.15,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,132,.,.,.,.,.,14.97,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,143.98,.,.,.,.,.,12.57,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +2,1,0,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,24,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,48,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,72,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,96,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,0,120,.,.,.,.,.,3.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,1,120,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,0,120.98,.,.,.,.,.,5.84,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,121.98,.,.,.,.,.,6.54,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,126,.,.,.,.,.,6.14,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,129.02,.,.,.,.,.,6.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,132.02,.,.,.,.,.,4.44,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,144,.,.,.,.,.,3.76,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +3,1,0,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,0,120.08,.,.,.,.,.,4.06,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,121.07,.,.,.,.,.,3.24,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,122.08,.,.,.,.,.,3.09,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,126.08,.,.,.,.,.,7.98,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,129.05,.,.,.,.,.,7.23,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,132.1,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,144.08,.,.,.,.,.,3.82,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +4,1,0,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,24,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,48,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,72,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,96,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,0,120,.,.,.,.,.,2.1,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,1,120,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,0,121,.,.,.,.,.,3.05,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,122.02,.,.,.,.,.,5.21,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,126,.,.,.,.,.,5.09,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,129.03,.,.,.,.,.,4.24,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,132,.,.,.,.,.,3.69,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,144.02,.,.,.,.,.,1.96,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +5,1,0,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,24,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,48,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,72,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,96,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,1,120,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,0,121,.,.,.,.,.,2.64,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,122,.,.,.,.,.,4.8,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,126,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,129.02,.,.,.,.,.,4.13,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,132,.,.,.,.,.,2.81,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,144,.,.,.,.,.,2.21,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +6,1,0,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,24,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,48,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,72,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,96,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,0,120,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,1,120,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,0,121,.,.,.,.,.,6.89,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,121.98,.,.,.,.,.,6.64,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,126,.,.,.,.,.,13.72,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,129,.,.,.,.,.,12.69,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,131.98,.,.,.,.,.,10.58,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,144.98,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +7,1,0,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,24,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,48,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,72,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,96,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,0,120,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,1,120,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,0,121.03,.,.,.,.,.,4.46,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,122.03,.,.,.,.,.,4.54,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,126.02,.,.,.,.,.,12.19,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,129.08,.,.,.,.,.,12.1,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,132.03,.,.,.,.,.,8.61,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,144.03,.,.,.,.,.,6.37,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +8,1,0,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,24,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,48,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,72,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,96,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,0,120,.,.,.,.,.,6.19,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,1,120,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,0,121.03,.,.,.,.,.,6.33,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,122,.,.,.,.,.,6.24,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,125.98,.,.,.,.,.,13.03,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,128.98,.,.,.,.,.,11.86,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,132,.,.,.,.,.,11.45,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,143.98,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +9,1,0,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,24,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,48,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,72,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,96,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,0,120,.,.,.,.,.,2.85,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,1,120,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,0,120.97,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,122,.,.,.,.,.,6.65,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,125.98,.,.,.,.,.,6.81,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,128.98,.,.,.,.,.,6.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,132,.,.,.,.,.,7.48,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,143.98,.,.,.,.,.,4.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +10,1,0,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,24,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,48,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,72,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,96,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,1,120,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,0,121,.,.,.,.,.,4.36,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,122.02,.,.,.,.,.,7.79,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,126,.,.,.,.,.,11.02,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,129,.,.,.,.,.,8.86,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,131.97,.,.,.,.,.,6.09,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,144,.,.,.,.,.,4.15,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +11,1,0,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,24,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,48,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,72,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,96,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,0,120,.,.,.,.,.,2.09,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,1,120,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,0,121.03,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,122,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,125.98,.,.,.,.,.,7.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,129,.,.,.,.,.,6.31,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,132,.,.,.,.,.,5.82,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,144.13,.,.,.,.,.,2.63,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +12,1,0,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,24,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,48,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,72,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,96,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,0,120,.,.,.,.,.,7.09,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,1,120,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,0,121.03,.,.,.,.,.,6.18,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,122.13,.,.,.,.,.,8.66,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,126,.,.,.,.,.,11.16,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,129,.,.,.,.,.,9.51,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,132,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,144,.,.,.,.,.,7.89,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +13,1,0,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,24,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,48,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,72,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,96,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,0,120,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,1,120,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,0,121,.,.,.,.,.,3.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,122,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,126,.,.,.,.,.,10.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,129.02,.,.,.,.,.,12.84,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,132,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,144,.,.,.,.,.,8.06,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +14,1,0,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,24,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,48,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,72,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,96,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,0,120,.,.,.,.,.,3.63,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,1,120,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,0,121,.,.,.,.,.,4.49,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,122,.,.,.,.,.,5.5,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,126,.,.,.,.,.,7.28,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,129,.,.,.,.,.,5.27,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,132,.,.,.,.,.,4.89,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,144,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +15,1,0,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,24,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,48,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,72,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,96,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,0,120,.,.,.,.,.,5.53,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,1,120,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,0,121,.,.,.,.,.,4.81,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,122,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,126,.,.,.,.,.,9.96,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,129,.,.,.,.,.,8.55,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,132.05,.,.,.,.,.,7.54,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,144.05,.,.,.,.,.,5.74,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +16,1,0,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,24,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,48,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,72,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,96,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,0,120,.,.,.,.,.,5.48,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,1,120,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,0,121,.,.,.,.,.,6.59,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,122,.,.,.,.,.,8.91,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,126,.,.,.,.,.,10.57,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,129,.,.,.,.,.,9.52,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,132,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,143.97,.,.,.,.,.,4.96,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +17,1,0,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,24,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,48,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,72,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,96,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,0,120,.,.,.,.,.,2.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,1,120,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,0,121.02,.,.,.,.,.,1.86,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,122.02,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,126,.,.,.,.,.,9.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,129,.,.,.,.,.,6.96,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,132,.,.,.,.,.,5.64,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,144.08,.,.,.,.,.,3.59,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +18,1,0,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,24,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,48,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,72,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,96,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,0,120,.,.,.,.,.,7.95,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,1,120,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,0,120.98,.,.,.,.,.,7.47,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,121.98,.,.,.,.,.,8.67,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,126,.,.,.,.,.,13.83,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,129.17,.,.,.,.,.,14.01,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,132.17,.,.,.,.,.,8.97,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,143.97,.,.,.,.,.,8.4,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +19,1,0,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,24,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,48,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,72,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,96,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,0,120,.,.,.,.,.,5.42,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,1,120,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,0,121,.,.,.,.,.,7.08,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,122,.,.,.,.,.,7.27,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,125.98,.,.,.,.,.,20.07,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,128.98,.,.,.,.,.,18.24,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,132,.,.,.,.,.,15.36,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,144,.,.,.,.,.,10.92,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +20,1,0,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,24,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,48,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,72,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,96,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,0,120,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,1,120,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,0,120.77,.,.,.,.,.,4.5,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,121.75,.,.,.,.,.,3.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,125.67,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,128.67,.,.,.,.,.,11.56,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,143.67,.,.,.,.,.,6.45,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 \ No newline at end of file diff --git a/Examples/src/ptaex1.csv b/Examples/src/ptaex1.csv new file mode 100644 index 000000000..93a207c9a --- /dev/null +++ b/Examples/src/ptaex1.csv @@ -0,0 +1,10 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, +#ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT +1,1,0,0,600,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 +2,1,0,0,1200,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +2,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 +3,1,0,0,300,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +3,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 +4,1,0,0,600,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +4,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 \ No newline at end of file diff --git a/Examples/src/simTemp.csv b/Examples/src/simTemp.csv new file mode 100644 index 000000000..a4d379b8a --- /dev/null +++ b/Examples/src/simTemp.csv @@ -0,0 +1,10 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, +#ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT +1,1,0,0,500,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. +2,1,0,0,1000,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +2,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. +3,1,0,0,250,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +3,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. +4,1,0,0,500,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +4,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. \ No newline at end of file diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R new file mode 100644 index 000000000..032d36c51 --- /dev/null +++ b/inst/Examples/Rscript/bestdose_simple_test.R @@ -0,0 +1,37 @@ +library(Pmetrics) + +setwd("inst/Examples/Runs") + +mod_onecomp <- PM_model$new( + pri = list( + ke = ab(0.001, 3.0), + v = ab(25.0, 250.0) + ), + eqn = function() { + dx[1] <- -ke * X[1] + B[1] + }, + out = function() { + Y[1] <- X[1] / v + }, + err = list( + additive(1, c(0, 0.20, 0, 0)) + ) +) + + +past_file <- "../src/bestdose_past.csv" +target_file <- "../src/bestdose_target.csv" +prior_file <- "../src/bestdose_prior.csv" + + +bd <- PM_bestdose$new( + prior = prior_file, + model = mod_onecomp, + past_data = past_file, + target = target_file, + dose_range = list(min = 0, max = 300), + bias_weight = 0.0, + target_type = "concentration" +) + +bd$print() diff --git a/inst/Examples/Runs/bestdose_result.rds b/inst/Examples/Runs/bestdose_result.rds new file mode 100644 index 0000000000000000000000000000000000000000..ead529e5b16b3f8b3bdc1122b36ffc7cca3d006b GIT binary patch literal 34108 zcmXV1c_7pO|930*5kgEPSFD2++Q?B+8H$|iR+J@*G8;;+*X#Lu9*@W4`2NwS&SgcRpudUh531_>}kr$>)`jhnK$}TptOe6z?mmDcq}+w8b$wMOx2M ze`SyRO;*>Vslcv;^UmTKI9HhTBBdo=s`1^TzkjI%r2&;Zdc`Z#uN~qWpE1mQwgh{a z|Adi~frWbykNg->{OGkm`ynj#?S+3Yl|qjP)-HUDwCZ#8lGCrf>wZKpz58{{^X^lh z4AkfI+IwzvI&9mXedNH=Gh^7x7Y}xD<8>BEEL$04`_)G_cC+Tr%Tr}L2ti||0cN^K zUsT>MIgp*K?_uK3=x>-)Ib>g7s#06IfLvNonfs??@~E-a`HvTMauV~mWr?>socYmz z?fS{qpQjv-d`uSO)DR!8zA0Y$r*~h*(IXR^=NHZx*IzeODSPnXC*upt;zu~v@SOgI zcj>X^pAVnAP_|!Rb1^ml?qkL|5~{;>+~jeeZ5_NX>~)|{->+X`>bmRrp#J5X$eE)ymt+r$=ehRm?D<^pSq|57@o_V;=>DnflE2|3(b#VF_=;K> zr!70m4Q=xy>!w+%+?jg&cag6N#(k8QW2Mrs4;>jV&D&$5bz-fx+Y@v14Pm(3KL5nP z+X{Ziw*^%9z|F|_x6G|}E4g^L289TWoh)d(K|+>xY@quA{myKlw;-d|MiS z$A|}gL+l-O@W~z;w0t~Ms@cBwvUl(9@P``Dd}^$pz1`_)9VK$xefG-T`zfE~H3!pu zOl40!&D#@HIU*7A;=z|0hS?mwJ47~Ur7q^j*u@W)A(;Yb~zb6xVd(~-M1pR%T^Vog3gizKLMX}?`dxh>z^nC#5sx1lf}`HwQj)9xO1 zX`-W&=Z@QFPux#;sKXc@oaw{8eUV|-88U+)iF9*;h4eLv8B$MgmD2YHj$9ojWUH;gom5YhSYo?;`7mokHA_HWaN zc|@FqCFSQWn)msxRV{2Cd2u!J9|@{(>&Im`Dfzc|A4nv8b0UTn|69|B`g+B9xVPhb zJ>|A{cgsoz&aBblP}Oy>1>ZB@+8 zR7?Ljz0dpZIdY|VVP@va^~%pJ8?SpFyd0w@wiH}aEw|P_aZ4+ug_HKCxL)tx7W$z+ zeXZ*z@Zd_Tqgj_9#44+P-)9sO!?GwH*&nR`Eg|8-5sB=QA+3D=YtiNN=;!%o{j%38 zJY6Efb)4_n5>*RKKDgC%8$7tsM{;=B+9pA3`_XmTf?n8_xlLoCX6LhQK_4f2p2+r} zYr0VL?CPW_u*e(xHL7!SCb9DD z+1{S?SBftJYww)usTg1)2DORj%kRx$Yo(8yRxSuXex_rdf7?egA8OyxFG(GI5WYMs@Y{(W7Ipx8$Rce`H?F|FP93{R!&K z3%uB+ee#N9lX8vLo7>QwLY8FXwVL7KYX@?UhN(;1lzljB=~I|?nZvugKlPu!__a;y zMK?DWUdz9A_U09xOO+pF?$m9wR!z|lzwMKw@o2S~bLGR-bDabKZoIuE`)40kJ2%z0 zphEA$jn9!Q5q>w!Upx%!d#f0?O1<5!G^n0-=S#-6*DAHWk01V6IUdXkcyulBKa^-7 zN9Mz%$b!7;P)E+Riu$#o6d!w@rPj5dK?th@Hlq)eZ)UE23Ok}im9m6owGSNHvRJY9 zhxW;%@8ZoLhpUcyT9=#V-oeL|V7%1>isBLTpAh#L$IOBgI_}?EShhG({v!D5jpjOB z&?ZVq61-MaZ#(Bm^<723;bD`5b1!~$o;J&tRP))pg4uFp#+Eu-`_;R8*5bq(b^pCq z)8ao)>|C{ChQCvQ)B1YSCr^>u>H64M`upHy{Bh>iMi1h@K3}>rpWbY@L-mdx=V@0 zzjl< zRMjHQ-I>htG%HC;)=s#5;G#HRJVmzto7rl*p=o+e!WGS5$Df-Iieq+WefPYk^YuvV zyP=>`t25I0+2RwL#kQ}8_YyGT{P@gY`}0TMs=uT<{BUTv86uvOdNHp8`Qms>_P77G zAUc1P*QXaGxjib+8c6W&I9Pf6Yk!~L%`ah6XtVu4R9shN)o+%5Q{_Ns&7_Z4_c!f@ zmCV@QkO-K1G!&A*@XAi~zx?f=2Z|fWw!-Y2D4EX0+j^)|Za+#aLb+)Pm3|tSLyz}< zVi{jBei7#16|vlvv2pT<-_M6h{-=nd@Wsy}(9jRf6n&R-5z)Q*Uf)h#XAc~F563?L zQMA|DdPQ`|F?>z_PYRnrZIL+ttz7G-sjTFumlB^VD=o#MD-}-MH5m(MHhYfjnvE+p zn^BtJ=0iRDrf;?5pzM+Jl`N+js}sA_!!F0m)hm+9pKm0)QfQd{5@Hw34EDy|y71uW z^$$`NW4=clD)(3{qQaiJ^6e~+oQSHQf?aRzS%FuN4)(9EnAVJD1)0hIEO%{+tba3H zWO_~JC*G^&>iBZJqKiZDgJbsAP31bF(KE*V2Se`Ny7Ql=4afWQqR(GAW1;cL zlD^S{T`GpRw@zMwdufheyRUkwKD2Mpe)(ECvQMd_%400-Nbng9dBi+fJmIn9(YNPC zqfd-Pjt^pw96z*cclI~`wr3Kw;jJXs-3l(Ps?`(!Vy?_4b-Pc*Jh`fCwbkWTcM+A} zHNGYHsp^-^oqkjL=i1+e-+Ytcdfj>~h_+>ML1xyfVkk(;HRp;(w(n8dPOD+B7ukJ9 zyK9z8O>g;7tu`|KZob<}E+6~UQ9>bIY2TF@7QA#ZV?Ou#H!ri+NVzGW95=J5BQNXE zG#yk4h7BE?>Si2(1s{`3-Er5q?(x7bwes&ij*K7s%y%njon+kT9lC5$8Gw;poBq7* z<^20&<5z4{mb4QqZn(A`cvOYUwiF+X8KJxi%y*KI%=@N{MK@>}b0yXPd}@1idRv*z zkqwEETFI)fSS`w|n$e8xb1lq;K78g*+iGNC_NR{i%paFiJSO=rcPvvv_eO>sZ*9QN`^?kb zH<}Q28S37b7XN+nYda%6+`x1k>w5h_#LYPWhqC$?%-`=Asx3N_>Ka}6*7x>$_K%^Y z=wSLo(w-kL%p~sy(%-uF+iOPhi=WCN-lo5s>THb-^9ZZZY$|$rmG$QI2H=m5L zC;F_1u4kt1Jus9yJH1en{7W%5*U~!g;D74Ei%C{*)UR0ovDdbi$}v)p!^@nG2y91; zj#2$WX}^&N&bvCp*O40fUM05n?mvo;T#^6Te~9&H1cYHj*Dm1=XU`mKP!U1`u;b)yV|qLR5#9@KRz;D3fq}JQQmOx>4?2*-@($@ zbIC>#iMib}f!kJUs1uW5pNsoL znU{-NI>duYEen6ki&s-8`V96oQ}(*W7CvQ+*&6-&H&8;>u2CKS@!DI{`0RZN7LRzR zbU%NeJG(!%Z{%O)`n9}Szq3oBPSb7(w?l{WLeJl83K^uUt-TwDU1(h_XvTi;Qa=@u zmMBVT4rco=7k^3LJ@a}rp7ee1UBljyye49E&Ghq!dY>Qe-W+sIHi%9>dw-#lvo(ux z_`c<@T_oJX$F>6po>_7Je7RV=e1maDJX}J!zIW0oBui8LcY${qX00jkx!>1ab=S9S zD_hwt9r54IKd^t7hFLZ7VYReBs%u9URIhuEqV-Pnl;*tek9(p=7STzgJ&+A+uZ?c( zyEeDeV5e(;<0yB=q+mfk+H)c56-j61XMn7k%z$~ts&=jMK;NH<=4i%65_$N8&l=w_Q<#hTEampJWuj^_#q%4Bt{t zFf})7dyd?;5|}+z`6!qkZuN5dl>I#BSpF5s3oSpt8=1V|DXx_!jGV=ykrU5vS!LLKd75ajab<|v`220{8Qm89e`QwdBxHQM&t6hOx)pm9S^hnF z-QdCF(zTqzXjH^Z!`NbRfxArLE{WZ@87D7kN}z+y_VvZ)&YrE6+i=>m6rJ=@O`rL) zsBSobbxK@*Urp=Y8z$NY#vE_q<%bhl1uUz;7vru6T-f_}pZsE0u6EsK&rjW%Nz=X@ z)l-$!E7&A|n~{?Sck{lDm=%01VeB_K82(|Hq-D21;@Yq18&O%eR@?rmzg0?p6t!nM z>{Nf`mdF^r3qAjBMTtoH#z}GETO$Tn4l2`bG(Nkz^+clt$B&?M(6W%gSWw86vx;91 zdIY`PW`-rze~C(QjS`J8pZ*qbr%jJ#o)pDbL$Tl*QY5;=m+W^KBHq>2kIrVRZevYg%aw%199Q_ga&T z(V6?+rbYW>Va!X^>OVahC*^k1K>jBCz%HDq& zDr#IFP!D?c27=b6+}L$L|Z;t zk8S+A=h_v&ZT;J<^9Duug<^-!RzHW5E;LU_PZde;WO%) zd%18<**0&&#P7<8MYUFs&L?Z?6z&9h?9~BsT-^_!-qN{;505Yz=wx}nnVQqyz9;3U zvt?@zSA2Wl!EKy8`hrtixAUY+llGmw9ZyZAx_<-|xHVm_ z_z;$eDUw!iX`EJ%`1AG7*MnJm2`<(5!s$BSc2^#V3YsX%Pe&uIG+Zuk!_Qu?dKE&t zh1gptIj&VWJlM)R^de~bYUMELo_<0`lYJ9sv8$`dlD+!kTF0Lc4Vh0G%FSXe&JaAF z4rThPC;X@1|4{sk-R|W|;r$cVd$}8F(K^|C6=*OMK^}QK?S3LY9xi-% z0LgPZTUFGEqr$bcw!O7^ySsGCHxQ3G!0W=?BWM3;B3J7;cEPqA4U(5eOHziYCtVha z^g+#6$OX4QA;);PtPZsj*3X2V6*tFmgu|Xy??9IZEmYdBe|&UZWh`Heo!P5)W4*wbl0h?zm4Sgr3Wm!sgTU8~5!OMeAYT z(es5l6^)sjDOJj>pZDQvIiWE6zIb9x4O|lUnLgLmZp&Z3PgEkcgZtH_(5Bdmi0i_c zPl;ER?xc=x8T)loFX_+znImtsPc%@a`jT)aW9OYsbqPBoU(i*@ltv4rI~7Te1F!YZ z;^ih?{Z_3V5NaQ7SA61$>7R>ZnW4LODdVX}*;~BH-;X&(f5%tTJ?;}k_|u<=A|#D} z?v!wXMlG(gtGL02#ma1R+hUrJRQ2I|q-U6y^li&3aUHq=*zE1tV*10E4G~gk1DlD< zqn4YHCnCQdYQQM*&BxR!qPQ~p&-lP0dG4=9iUQ8b7M}$+NKj{tXY$CREFard=;Ua@ z=1zqi&c?xv+v<_4OJ-fIZx~@E^i_AnS8o)b*GN&t-Oj~H-Y4kU=&dc;H&Rq`7}!y` zLICbc``6+I^bRg|Y}$HnCjJFo`hK8|3^%_Kz6E!#8A#bFlx!P+ zk}peQOkGWsW`fP5KlE={-%Q(KGodt^wuyX()glb}8*!D!GR5YQI@mJ$!$1*&0?Co} z%$6f7O*6KcxY}S4k2dCZ$KZTemm1;9FcbVK*8Z_Uk&V`Pv^L|n_QSN``xtPFYIN#m z@?n;W`;jSv88(6bFd+v2mTnS{RyU-T#iR9YCbo?dOt?W~UXm*q1zalqVWtR#9ivK7 zsx(<@Qsk47Uqu^ehTPdkf_59E-$>YQGjVMcy_5UAk)Ym%kEg$E+!r9ieSaSf=idEh z&(69})8RfyL?X~8Tqj_53I63V%gtmX6v_8K(k1fiYy;Vi^}P`zN3WD3J4SwGH_*($ ziR?2IznW;^oOzyWG~k8eI$I8%(+i9wVBWuT~foP(H~03K)Lj&_`nYm{HDew zWgNNz2~)w5=?}L;Ff~4Qj3VFMY7mdq;l{@B3*wQ6+%L8qoF~h^5hj%Zg*C!db3$Lx zZJyzUS+g-UIF>GbpLj1+*9R$GKNzrsDGshhDX`+>QeRWTdNX4&B0Z}uUiZniEAUM- zNHMH(?i_L#H~uxHUuJ_QMpF6bPL2sA$+{FbfRE&@Y|1|(i*-0Z<<^B3E}Y^#qpKZb zn=#ulnM+X(aEVjAC24%=40;>mw?QAYo$exBfkTRrPN3QYi4^?tS-j(jaN#0H*%`lo zCYT4Qv7W}g*9+htlH>Z=rjiLFcq^UdIv86$pP(7}9h%O?OXO}WseuRyiDB|2SSMn9 zNHh_$H;Ss}N1{1qqb+2a;L1Jzep@qG>lmcYt(9GYWELy6Q1k0*V|Tt+qfirB6BDoo zxM&ZHWkBx`w?I_kus(Z(Mk}y+xin$JsS6+3E^n}ztd7>GGnr@!!>NhARot6jRvXo< zqkFUq%+hb4?%{SEUdiWv2wu)~W}C@o5H(nrz)CbmGo_t*i=Tc3l$^qNYsFY6uMk{>d--=030Ea}vI@vs5;(mXv|J&B{9K7RTcQx{Gz<&lm%|ME&6u-s zcX1co2Bf&h?fhvs2Qi`~xkyu3sCbkENdYy1k)%Hx)hU-}*PYF5Q;wMvW6{nV-0v$FvJrO^e|N0lz#@?)tIg{-#fva^KrbUDuv#-F_+nN9W~Nst2#S6{pYxhJ4SA138lw1n zD>PjCzVlGDaK3ADohpjAe0X$TY?b@sj!*5Dx-*ljF-3U2eCpPf3M^a#_c?bB`3DJq zwX!$eQsh%~fOTC)w~5R~5Ok@*aS@X%Ka(fekUz3BVXE2ZZ6 zk6SaTBD9@_ye!_jl`~O|9b~(@^gUxFsEf))%X2lpq@o~S!+N*LETLX9CO7=O9eRICHz*5!n*Y^=ItBbqn~>5(qaU@VQ`6C<2?EB051ALH^x!m8|+>5 zTXA=xn2BX!XM!5r&o&jJxNW?)LwBs-8J*s2Cd5JSp1HtL@}hriGqoE={y5-G*B{js zA8g>?%OL48r(%5gaF2}_cIahDioO`9bryOr!ISHPG%_&yxnHUF5>Ot$sevfb!(Tt) zjD29Cv+=MW(tfpaJ>6BUnZIrbc2^2TOaK4wxJSY6EJ!5LI{hl-_$qnHeS_3Skvs)q z;lLLl1pQVnk(<_xtP3m@22ZT|pDzTCL)+8$HL~x1fs69sg^=%9f+J#E=n{UAv5e2) z8#3?ORiKZ}%BLSaR|<{=Qe;)fr8MvfhCIhf*%%E7|8~Xl9%%rQrl-c$;HKAj@9kjJ z>1GwxWUP4Z*H&8Wl0pv$?#IWY?jD{D7Pdg})n60YNMWvDms(D>45W+RwV1gDq0Xql zpNpp(jd>`dIjWu>FLc3+uh0&68)ajFmlOygy%Yv2Zy4uWrCSO)b5#4>t-j|H*XCm| zTkfx!W?$|~Si4k$nlcf}kCrYFrETA|=zH{n+t{vjqId@QT<}hc2+nNA1W6r5Pp)Eb zI7B*Ck-hqLxNc+W=Wx6Hj=%GsTo${9?_6{Vo@(rn*?589N0$~i{W1Od0xZCoDPg;c zF`=i%h@Z4p>-MR+@vnCk7iPzFaOBXKI z{>^w>nXPNSqO10&IKi@fJKn7Zg0SIo@YwSV7B*?U>Q-D0d19lupxg}v-2k*7Wz~^O zWMjl|pLMdAR73X8G_g&?;;kJ=#0nQH)mbfZcM+7Sf_kSh&-6;sBB)aicyFZ$T5_X( zX$z#4c7Q=NT*i#uX?#xu-g~GUv5R0_xK6mUFd)qhx4lm0Ww3JkHXG44r%iK|<0TEX z%;_&6yB}ny_WZ&rGAi=B*PUMS*74w1v1Iez2r+iYMQ-)U$%i-auRdnE6Wi1_&9Sx9 zkcTqd1zBI3tNZ2hM_67DoT+M?x`c2J1>UeENji<{Bi_cs_RJ6^qQ|9(s%#DLsYxu_ zq>vJLlYQ0S^ME09HAaIbx_#pXq(e`!MYG?_@*NTSi|gTI1Z~D4vU{W-D{vEK6Q}a< zqp^r2kfG-BS|H+X@K3ev*^o)6iWVWB7SJVt(INE)ifVZSnet}kS~^_h6ofDVbp|>3 z53`@CD&9QI^bPcce_xzVic#7h|JQc>q`k*M!x)3rOBbeO$`tCJlm7LJcaR@LuV_-Ao z5ecXN6uZTpUWXL>s6v zBqmG|8q(p)Dq{H4EcTXsN+Px z`!DG5CmZ_?V^w=a$y-tLNKMvdyIrti@08z3Se4x0NK3Fhq$?}ag&m-6NYfov6AfZ( z#9&mLnd!y4y}sgwQ^8U{s7zB;z6#ii4G4huKv!i@Ip3!YN`%ssR<@ph$}PN3pK0w8!(-0iy=qgNJYF9BY$0H zOH2`D?%WCzJfLl4&$*UDPFPN@nJ5mC`*~@q0lsBOXv06_wJvQ_A>)=Ka3Hubcm55Q zD4ju+XR(37aHDEcjl9g)dKo2K*iyEu=?$=#oM&$F}V zAHbnaS8ETBH)p;ca1|dw-|0o}WVA4!19Y`*93w$5!oOUB;hYNPP_2YODI)LWonHLs z6^gHJC*RF&86R?%yJ?Q5z*v$$jlMIE<^0sa=PlrAhBv?ZTtU8orvw-=t!!Vyu8@rH zRtR8|66;IxH5fNOf5D#sPZ*cJ+cC_>LTu;Tc+Sc zRB{wMU6V%~W^bocrIPUr42CBl{Nt=h&Hqm=!}`F6iv)}WfV1mJ8$?}vx0WtVsF1PR zC&(;vqro6~*ES<7H{&%~2^|rZ`4K4~TX#Ta&S-OBga}czL6ySXy0-Ir@{mXmV>L%1 zx&jZU%f{6pe>3X^d1?d%9GTGF3cxYQRF#{?`Q!I85WASSZLbsLr6wV87DKr^>P~aFDvIQDt<@BNh(*CXtdrGGJbf@j)VB8&3sE zQJ(cvEHEY-{2JR9C{o00HRSaG^nn-R?y{>wo>{SXMyj&0wyU(X!{C8wvTTYZzW_xR zj|qYTC7=ejVc2JgrfZ`&1Q~B=F~a_b9uEo7BR*kb8ZzX$O@@ix;yH48KR_gjhr;f3w}w+Ux5`V3a3XBh@n9ixk-5c9DqEDjC-ExKU{Q=IDG3 zmC0E4HpYl#tQnPvfme(`(+s+0$Vpb0ptp(l{@=c!mPG+=5(U7jGkaY$wnitvoUR&!WGzQkD7kq8uEgwFrC>VX{4pfcDKDw zv)>L}`y<#4Ud^8O3?!ohb#BT=rVx1>ibI{sD0xf4s6{(1puJFYir{n>&~BrFQs<#L ziT~e?S5n<()L?-cK!5fM$UV+HU6{lvXeA!(hWB?jUP<^7#sBT*pZ`JoAUZxZ7HQ0z z4@swKF^UN(uBaN%$$h{GRMg$iz~Z3D?6Qj_RkEy{T?PDQm9{EaxtY$inOIA6|x^I>OR?E zydK@eHHinZQ&j!h^z!<+mnc%8??holNGE})QujIP#141np>&Ew&swx8Jq(<6fF#$( zc9CYO3@s+|tle%`(=-{2m=w&b6{mC%;CxK(iswDj#(HW4b#uwXNZ#^^j1XnEv8^do zELie>d;a?v{QI7z=teReg?vFtrkOE)Y>Qz8*zg4qHNb05Ds&nozeB3gOqmrR^dlaC zj8H4fZeojt-RAvK=96rBGf@nlE>XN^1LUM70k_2gLL3;&lSgg*(>4!u2}vVfH#l6R z828#2IXW%DDEic-iEsrTzXQNc{1#TP0Dd5k&HmIyih}%SeZ4~_tyn-ym-T%e^ER-U_^Urnl zV|Y2>%yNGNpgg_#k0=K%!QVC7%lFFVa{!+uOL2v5Kcf#1QH`RT8s76TF?>6auW9Dx z6uT}+wH(?GFivGG*_Qq((GrP;vCKfh&=y4)_iHjXTICISM<@3HL(1E)RB(}dnJ`WJ00Z5ZxCLW4oS?A-Ev8K6SDM&$%M8+chKkge}BPumM)=eh@ z+~j>YIPJ00JyL(5erwf95yVy`E zeIEM&r5mzADm?APw0lAfs zumgcKw9)%I?U<3X6U0v{0YD}Uh#yx`5I+;kO3s-{0=uD|XR8GHLxOwFZWw(78*p+~ zN`SZ-zf=IjAqtT^6xP)A%b%IdDzkUmnj(v*od1$VNv=0YmWHn?&K&cb^cA+?KSI^C z1N%=jV3hCk9KSu4i0()4_Ij+qt_F-k2spd81m|7@LA-9j)s_d4C`)TEYe}Tb)sEx0 zsc_?}J=bs3M(cJ2pnN>#2Q1R26=EiTdHb097`J}XI|Bh1>>({s6fN1t)g8~VpX~*$ zm%?d}50Ia;Wg>Qe?f?XLCYsL#vw?*l<&GdH)}plY)e^D7xu2nBqxv?mHUK3Lf7$sF zkRY5gkhGXe>2Qf2E+FSgqZ%Kf^P|&3;5#C#0Oh(FjUjA5pteAc?ap#Q$fFj? z>D`h7G<2eYqSnKB(mL&HG;GE^5o3j*sRA2gBynY31bW*;H5 z(MnP9CD~Q72x@|JNCNT~-dMV~mLyWdkAA>PisMUeN7M?RhU_1({z&lBjzR3B&FTad zZ$j%I*Uf)XD6o7xOVHmW_~7tt zur5KMp#>#^6!8MoMKbN5jnXs|pi9%N2H$n^7BDe4ryC}_F$Y1S zxEG_lK0=cB5G9l(6Lty+N(%@njPzBY!tfF#34o*3-a&!+ye$?8XnWlTD3SQX+zzHQ zu=@{Nj>L*L^*7RhNpZw7ynS|Zr{e?t&Is^Q&$q?KoZS%T0@TwNwMF@Pu=(6)wtqy!(O9&Y*c&Kym)9&ZGM9~bE zp*|tA50o2dUo2rKNPj>>l4pGLoAm=iNf5s*mSEbyX>D2vEm}7~Fj>W;23+07%fesL z=PajAG(HYqE2en<15Sy#9C^CBG4WI@0MOW&ap?r_BWe`eGw;2iIN^ z!)s&w66RvTE}sI%)U&Mr$W`OLOGI*=I%*=!f^Hh4K?sfkbdqr_nBg8EkJz=iexn&i zDDbfI00!28m0va#sOjK?1*e)VNPB3t1rGLV)&tosMsSi}kST65a95d~286 zoM;oU_X4+?Q-2qmiK-dOCGyEz0JYjm$W&k1gGe$-6%>%=B@x``QJar=aX^*c^9p;Q z!M|(WGuDe~e~%;owW0wQ-9NyrHUPS^9hz`jM zOstbbO-#&Y~@{m;eoVNv<4& zf`%KR!;|Dh*2FSm3~yv}L65#J4VjFoM~FgpV|Y`W3%8n6DDt96aj+mg(0&hMA3z|T zXEzH10*jMj{fc`w6UK?&#*(o8t){%bc{EFGiknF};_^QPr;M6_Jhp>{x{8vuQJhxM z48Ce8|K!7y+b6+8aDm-G6f~_zRy{hEXAOpYE~a9f|(mX|nxfeoQyHJ^4O#Z)T6D;}E0l z!#0+Ntva&9I0i_~^P)lIy)R7KiyUzpx#bx3sngL30!ho13(FtDKV74L6*Z)`mI?cTGt;h5%aVDvC?b9U@FJ$k(hu^2?9zxGY2yF`Nif zZ@I_@y)(;G!D3_lo9d}gr=|t!aEGNN0H9TPKc`Z^{piUu;nRe*^Osy?HwF;PKo;6G z|77om-{mAM2qRh}z)u)q^3faK&06I4{#F-<5s5;+|2A-3nOo;1De`-!H!uft2LQs6 zvb_#;fZ9FhAuZoQgx(O|a9)8=Y=kVvHOcKW1D`B(2?mApfq#788!hewNP;RvMxm(C z5GmN`r5Vt~uS2&o)&#Y0R)8hh+cOD$x|d!Gq{7wf>cyAZ1ua#a0o?@b1^4UjJ;V?# zu+x+#XjG>*yoYthv0hhyxc0uyQnZNwOK{SjqZHL4*)=vdIJ0dXB!-vDV9^ok1A+4N zC>cg97m-9JZ>AJZM3{z38v5%oL6pKIV+q1-Hmj9;5y;xZll#Pr6V>nuM*fCON1#xM zIA>}CZ-}TpV571AhNiFt$FonXL6c(ce<22+9$FPr~&!*xmFxD`i8%;YgtDR2x z^J*SizYpB+&z=k@d~|6aOP=ERGO@Nf$!hhvM zg?41|l_ZlW)MT-C&Ek6KXr|O>z+vrXJ#^7yuLhAwLF)`J1>O}H$CEUS4w!r?_+GU@ zbGkD+50pPi@S6C*X7*HYhPEnum%Kk@<1p(BC=#@rz8a_6voxPUGGaPv*Q|1!1c~qM zc6>ITq`x{5F&bnOKApMrjio#BsN+5VLMQ8DM}gEhmnjMcqreZnabSdJ2_+Ct2Aoe0 zav!AZoR(7e&EAVZfOCM{*#R8@>W=Rh2${bUmJVO-T7@VAc1e+C!^Qbj6#0^N%^F^m zdsBXx0j)PzT3>iROlnbt_@na86lHceSSweNrLJAWxrbCl8#15AX%IfR&OiD@PR6Du zltd4>W*(8*@PWvJ*q~&7gsu*It;wfbM$K$YFA1Up_JDg1|6jH*SyjFlcPI7>vJkzA z>n>jr-R@!-TAP`OJJwu@G%T=s=9TP6yvMVZ*JvFJS;NbsBWh zb4=IXnN9U(yiLTuPbk3$2?iF78Yn{1Qw2$QXMLRbEv^y4sBQza=ZKOJRh@OG8`*vz zp2GUNbkisJ8C>bF@o#9B3$l&w=)CwCJLzUn3Ug;VWxZ<^eE}FBvv(2+W2wG_N#0Ba zb{)7+W!T^8#bh$#r#!^Zam~=e+;I6}*4u*>Sk36-f0%pbct$0#V#f|Dkl1n`K?{Qd zk*#O>SjnKzU(rz6iv$Kl4uuacg^!wkNzoxua(j%JDd58Ne?S~2E3>S+yl_ug1?~d{ zd%F+-#PPwjMdcS3W;~*i-_X1iw+u)FJ|_$ZM!mKe+k|&idUaM#4GjWNN)AJWbeJ6K zS@mbL`O6#`-aR=qN{>x-z^==;!&Py~bPrt;M{+7t2>H|AAId~aVe55)2-X!?qf&BLlq%QkP9#ICwI_UZxMg4ZI#ml+X_)K=# zg!8Ng)q}XdN4``3c^HrS`&BSxNM8iO5)sGpEGWLlsmLK|k*R_)yu11Ksa_Zp%81xVEIQAO_Vv=$>~V zjo4F>#C*YMM#s7p7$VFOBx*x#E z00KT46+cgCrAD~l&1{3Obgj~c1@*7?Z^@vU@W&9o)30htk!1QD8f19Uap?JLO5t>b z4GDK(H%@F70M%V%>bH0e47QOoL6Lpj?mAK%k2LLgl?$?v3`l^h{2QhWLcY1Nz*39L zLiFWB;hd#u!)O;U^OU87K)BAP=%N?EK}ric2*GUfCeVp|iT{N}zh1eN4qxtiiPV0x za$1m3_ky;6@E5+9%to*sp^BN{Nb0n=gowiUn$gQ2AqBb%z}@opyvoR#uMkz14<-`R zY&q9dhIX=lu1=-7Z-nMXyFtDUOA0n>?a*~+nhG1NX%&*oUlK@&5wSGW$gK-_N-lmY z$kx!SHaNZOJG5sMzNJoN0TNi1DcQDzF~8jdS~VlnyebI>=m$zH<`M?z6$y=ujbqFlGyr$5H^&iM8JPf zXYT2Iz+7ez6J+atLxDdCCQ}4WLx^^t%AJ;Ihvex?*F=#TTVbL$VQVaF09_SnY!}DP z`FqnTOjIi|+1)u6R01@puR0$f<11-iZbkxQuf2;=kp>zlAVuQiH1PlqonQ7^gLP1Q z4*Zm0-HEFqPks6sKz!{9%F-PF-FYArJR^Ft;@wH{vCQv4tN|;$EP`~su{|Sg;C0;U zE(TegVL$cZUoFtS%u=Qk=S;}918Gw~ZNeAfGUa4u5zUa{7?xj;2OY+eE2;tO0)qB_ zVF{oDg-){ZGZCb}!XNah8=5@SWSLF)MDRgSylz0RS3q9e63yZ$XK&+1$@x5b zYNMqSmeZ76X4vA07`CV4vtkKjuS$rD)K-9{*Vw zah}0HUYvSz1+chr`Ist3iaCpii-GX?*3oduVA}ux;E%EdQ}@bn5{Y9w#g??epbvB%bEJ6Qfl z=-t*g0-=7uRS0M#;C=QPBf5X1Yo81k`^B+~?*-=A>b}@}0glBM^n#>VPIl3V*g@op z*}^uj$eB*Eo**w7(oMwA;XOiS^Fi+Ew5Reb1Xry_0mWiXRWCkoz<}VR5yZK2b_ele zQK2f7w=3I3K@U)Jg|gbU2%}iAU~91dsu!bj5=!{eF|>(S)}*1Mv%M#G_487G#H~eb z&d*c_)1Q}K(~KD=GizioDI(j1xIuizECrgx8M4=q*xxunyUgz2s5#YMX+G`?E*_s% z1+Ss$Zgjlic8E;jpMcq!Pqu7~DAc*1j@i7jqyu@XBM0tdQdrcb2=^Knfj4oc@g5r) zKry9>gFq;u(N503)Fn)V#G;jf5Ono`Y)8udD!)#sjm6Nr7LiLLq+<%yrKPvwX>GnB zDXm8dMikbC$j3iwwMW3{TG};JQAqlGu~2R`=${IbO(hhqBN$MyyN4fRTmBwevjz^v zZ)&>KIpqAn`MY_EUr?Q;sJFGN`(C(7u0&D`N%1*IBO?ZAj^dA0xaK3ODaIPe{>p1Q2gWS#eKq(bO@M#_xFC%~< z%!8m7&QjqdqT7rauRX`c>$7VYD;8Ew5W5(QP6TV6nv>CxzozW%^v3*q#d#OZOP^)@4*5y#q(9x?;z+f_% ziOTimm&MTZA`uy2R!;H|0CTderZ_cw_GgcQbASz6+d@~hYj_t5w)1l&Hl~kZ!Pvas z0^FGWV%ZVVVQ)!V_oxwQ?i5KGF+fXMAboRhk}=R0alf1EC%efvCTE3(<>@eq{lt~c=_O+p()a!J`{hYBkK|LUQ> ztFUnZ(BMdZ-3;@7nL{lN>jJ)l>kd@|%?aFXHVn`4FmAZ*b#xxZTSm>CaG{_L>y1qf& z56FE1Jhk}ff5#9Bc3=_$C=#^H@m}0`^AZ+2hGvkMlC1MzqzB;U63Ly^0Aq4IY-s}> z3_M2xTGhJ}FYfJa>gBLhGg>n`ka=j4-VIuGR|P%#yLYR1*PC91R5g!Ddih^`UIxF& z$B?bm$=^n~VkQ=s{3S0P{_^f#!uL8b_hG%jYc3xiDFAZMF25z&@FhTKau+ek0TxC( z5`>_K5#3OcEC5@CEOMOMOJh?>6=W2Ila&jl(Xf7q(;Eh7#w!>Ae-Qc1zysiM0bBrDWBegUy?=^h{^ zWg3A!Q`d}?Ufc(MFHLsR-MlniuAy_jRhOR9DgZYCxK7#P~1@Nyn zR_nNAI(^aPJqS=xVsszfGypmBBdl;0NeW|@=#Ufy^)PGhNrT2E>%e^4+6VkVJvNCY z(JwjGx?WVM#87rI;y|C0S|o5+q`wBZ>%X^IyX?WBeUfCie!*b?$JE|yZIhe=)mTg@ z-(7PK3E&ZDFl@8+skcJ;{e7s}NGO`s%@IWfENF~@Fio)rYMN+vC*0= zAQ*cE2*!ZGiJ~X&gs5Wvwq5?d>N9>qu@rP;!4q%@Y~6VgP?q_RVBw-0qW1;nhd6D3X%$OyyTVwp@tM^c@s-^k!CPn&zcM1dWDw0qmLL?P#GPB?^`7L-L4}8jkmK7C5nt5H2;PR zmN@_&=nV#^a{%EO1MScazlNjGKj2V?MP#G_EcFXeTIs_8BKBQ5+3Z$LyyIP+zHrh4 zEUSL6s-w^xa3~{eE2(f18ypPa2uB=1*tQlqhY*l=BKL^|NVi4ho$fnHUusC}i`9X$ z=U*T?n3DL#*C#08WFRt6#qS8Y)T8GSpzKBH@)bt#2w*(+AkP7cZWh#h;RCRD6hLiB zo7wvB4$89LNUOmj9~`)m91#J9Thntu^C@n;av8vtb1k)an3amT2Fkme3w%gLJYqGh z4zqsZ|L{@<7_vXbWfu@AuSwr{cUAW#jy4zq;^d1 zUW$w66?fk%E2Ey>QXGLu%|a0~L+2wC7_Cgqo3g8=xZ}V{oAXPi?Oh3S3S*#QMBjf2 zNEm_qfC@iKitv{)gt&n*G~Rb99C9BVoS6Zsgg`#-XyY37U$Bh;5LzE*{0Q4bYy5_Nn3h6XqT=%)Z0L}r1( z0RsQCn+jM?9L63}mz=D|NzWk$u-zxiTfX@lMVvT52|WbigK3S^7#hAj0hpefe#!sY zP^B%0Dm|)mI(xRf<+GjW^ycUhc9wsU9qR;u09uLv`(gO528@k=qIIE&mjdf43q)2if~)1(`7<6G*H*1%Rpg^x2VSo ztfG2RHlgr!3Mssy*@`oi1gPV98YD*8#IiFm6`X@i6SNExgLhHKhO+>4r#=z=uV#5Z z!I!52=;XzK1W@vS5am^*a8p1vNezM?zwuU|t9jA*lhN5n%O!}3PYzJ@j%6%)L0I!! zbuR;YTP0RYNQ45s9zuYv4HT_~!}sq%^3CCcfQ-yiNIgQ`Km*W30hRH&C-U}i@1+uyD(UPNuu zPJg1DDJme210AU+Uy;z=Jr;4pyZqgqw0EZy9S?-EJ=~g_0Fb7^T_~@?q}`F=2EfnT zz3f>-A_I*#vqV9y_bxYr*h03>aG(+`J%Bc!>j7F5SoOe$CjKC+cPOCOkOz#+i1-G# zeNYvqt?HWtO#hZArvNn+LsG`m*-btIDe0s?lk~fN>b|w5|8)xW9*`Ae577TpK_L`a z8&OMC$@IZjb+TXxecyz3hS2S`I*#P()Av)i>BLT4r*?q=o+?C;!`j-G4#4d7dfqa+ z#?#2Tmh9|SjHthIL{tUA>NQCIBC>f_(hKY(kprAae#~`X8EG-jHT~P|sUVOGKhvF5 zW{oX2mL4c|Gh2rsqTqyHjrch=%T ztS3r$K`Z>Re+h;;aA`lQDF8)yoTe?m8dsrl!YrW$q!~d5u8WZp%+?p6nX&1EIP#HS zdhC+}^!w_%5nQTQKYOgaw`Va7IiW;XdaYddCDc_fWNGM0)UG0Wo_lX!A^Pw505+gG zb>k0oT}FTs31HRUqn+-x_5fuhuu1CCRSG$p@06={co^cZy*?3CiCNuxTGUE(^9;~D z{<#uu)4o@>!B~35%mCd0n2;qbn33x6mQX;e$_^d_5G41Y-FFN;{XJu#MTL>N;bj1L zf|y~6r$*?a=3UH#8)MUD2qSte|46nSV}MywzY6Xt0ScH{jUE}^%7!$xiQ@JQapL)jur$U41HfD$nEfwr=#t8$EDuq<*uo&klN;H zI?)KUMo**TrN+BFP>VVc6Fwmi``7mnY(s^H|VEopuK z$SakB*-v)X>vO_;n(`*M8&3?xE?mqvYN=OxP|3_VzF&2(X1hQ=_w6x64O?-Mt7zo-4|%z4>V(bAfyBSB@(NqHSHw)c=v9E-TSK=LbMf8Vxt(?WVGZ{fS`W0c5-dwYG z`N70MeH**u-&=Ryhio&4RctD7-J-^Vvyf=CFWP&9vg zo8HqffND(Fdu>|sZjkYiTg>S&t7v8n;>*j6meJmpe$O}m{2A;L@)+1ujv@1XK9$eM zRFvmH8?kve$MF&wC*A6|m&0e>*ctxNuVnw$#Y9~R(x2QaJm*5b_)_Xmg=cp>?ll>GFaHWwX&dRl?R341>6Q{qTr)%{2 zwyqw$4%vT!a@uSCC=Btl$*XkfUFK^?xV^ns_Led#R`SPZ14+=JB)Vl-o zQZULJGcE)476V=ir%>_Kj84W5^)fGX(?T-+(j!4;8+;#m4x`s7W{Qk`&4pwXWzf;m z9{)aj@C|EyA!Z6a*;C*W<2;s#!Zbp&hHD5^PCS@579!$m2Cpt^+h zY7)2oBMwilQWt$TFijOuyp?#S@iV)GN3^O!hykhNUs@0Zpkse%VZC@*px3&|Q{7~!1XXA^6ij3+6x=P<~Pj;-uOwwr; za_We3OeaU2ML-m-uTx@3n2OEW`FXx8+wk$*X606nF6tCAQ_^oflAc%m|9YdU>YDXl zevUe!vbJtv=vY)%6LWu+QTZ(`Xf!XwoINGxOm7tC&Tb=J@=s1hg*{=Hn(`srp4@N4 zfdfXjLuvnGDW|>$K9SF7#5t?s0fjme&zZY8!_H~KyJO(R?tC-!3i>U*KLR`zND>5W zQYWvtn5!f|+V!uLOpb21$+25RPs4{|03+Dd(XVA<2b2lL=6OE=R~{Q$=JWr&QPr}l zJqQHZlT=6(Xms+RBBCDkBL8EON?@Qj!a#3J)#e&p=@8hj5GvbRA{(Jx4KE14l6fNTX*PFTr-%D3=99B(wR`v2v)6^KhNtRxMyXYcp@FM% zC5gCiPDUwee4MZZ?K(c6dTXZ@3pp%i!fPWc{xeG58X~Wg9{*ItMM)c_+r7HOc>6_H zxuJ_)^|AU^Si!cOinga4jhS=E#1ll-?~aatdZUMdI_1{E6{kP47C(#AvUahY{Z9P0 zlJ#tY0~GK6M)YUc1rZb8f?lFX$uH7JH{Z{s~CyJ+pM8t89cMhr3Bo-05V}>pBVb}KBQg|XWX&k-PwP(vX z%Dizels}mPr_ONg>-56fOaeD6p1O#yiJOSbQMa1jf3Tc?uQA)KaaB~k};YfT^R*;k9q+UG=1_Qk}+OnTbs z7jwxwj8g_g>xe9q75_g5sipsyL24xY_Rb*n6|1N7O^S|FYR%7ez&?i&OSyY5My+~p zM0ANW``OnjCu6%BcUipO|1|yhuO9^u9m+!~Vsx*caMKk~fzgdwgSN}Cx&O-_b>Pk) zwau(F(d2z5;r-xWql1#NQ>wJ4TyyCAd)}^hn}xqKdf+m9-=FrLf3FXQ=UeCgIy?7mPpw3a{pXIF z(CG7TM#YRCN5$3$S~9%%@PI``)ZyjReh$CotX}BrOy#}IQ~F;&>-{{Y7i-s`t)5NU z&$==wD(%ZgH1Q)Ohzay_0E1MX8nTVX1?NrXCWgIO#U&*q z>hN)h^Co}j{u6e?UEcrvQK`9yQsyVzUcExHNE5|FgVo$-Uf23Bbl5im_Na%Yd{88x z6<##$60Papd?rVo{x2f;{?*pnSI*cY$a=z#dgD&ZYN#MyoY{91WG+ z4981JTrQrNE>nEknxwpMX5z;ns%-C0LgEjw9mU|%j#6D&KiNBb)VEdM@$QA5clM}> z(?}gB2WAd|^yA1RtlpP;O?n)?6iO=9DF+*2ZdbqGc+N(!33s5|ANbB3*R!QhwdC`6 za7_$d@!xr)3W)Y|9#~jKzbYl0A@r=M@ne~_J6LG7$lK7@9biziJRs`#F`~}g)6=rt z6zH#I*n7CW^?~c9-k%uS%@pgJ(Xr()ADB}QUsIca-{)86-Hh6UxkZhBU6rMV_Ir)^ z)uJuD{~Y28cD}sr-b_pQ>Uw93dh|V~Yl($E-oEN%a<<^=v&PAE$*Hktg(AETv}vG& z8Rh3LotD0T5Ai2)T=X(gbw2aC8dXw_Idkbud-kDx$EVh@HXJ1*iXYm1Cvt@M_B=3> zD~U046?6uciF)WXAGE<@Pl*l4prF}au( z_!wU#FW>4#>7g?T;TySWvV-#&I-VH9bAja(FXykXRgDg&`6w z*n#3WfB}Ah{pv=fjBp9Pj|Hugj`#c1%gc0~aqENF*XrDl z0ICf<^OqC!W^~Ldtu%3zzd9ZNOx@%9?Hqz_egN2uKEvrxQeFit_BxLnA{7^@^=tkM zJsN6fuPn!3QAKIx$L}k&Hn%Da19%E)bjV`Km}{t*149hy##Q07&I1BoXetfbR=dsW ze%#rGGY=g5m{C-4Wb?5{UDS(YW&P#twtEgbHtU4;Kj*XN`9xGIWVZr!c%seZWmc%E zh|NX6IXOAeyD6S83HP}B_=xzJvO&kc^4Cv+1N}^XJ;Js6{&m?)koeb{mCPJ@un|1w z6<{)@Xr_L66|RInsc!7Ik2G=z}zxzi+Vi3M<42a< zqX{OLubcAQ^>J+5$@!8Imlx8gM1$?y#{8EVCWVOtm4KJ(qNLr3HDp^gPIs#!iN;K@ zO=k__@9 z19$%I+9ph(EQM!QFf@|LUV(bc7N|d!b&#VKNzJLGm_N;1J^i1tDx%$E^K(TfNS!F! zuabY1z;eHn_qp)~(UAImJ-u=K>E7e>OMy6HPVM4v;-1H%C3fki2Clt1YS{|ZoK5>z zagC3-oSZYBq3X^i-Zx!>jD`6{P_Q>ubyt|uQjuZG`StKu=M-GirX0OMk>V7&C7x=r zRij&^H)s8rh9ZnO?~+)_Ep=kvw}pQi+N#aA8lHWu(d_raUHN2?Sa@z&3@RW|J1(_~ z7L(E~_|e{OeS$6Va#^|gR$uD2!hRp;tI5J+P&4OAuHRDH8ewArKJlUw5p~sV6CZaw zVV_#98&6}{$yxhn{j}fp4SY8NAHr#)Og82BUktLEyXvdfCnYS8KcXDfD7PIb-6_4Q z$;QUYO{UMMgd~v=ydI-MQt?1J|Gc&_o=VD*e!V(j8UXAxC|P~x#h{ax@<%>F1lbp7 zz+w1E(WoXZrz=jcB&1qo;Xfl*Xncs@uEy}9{h1=*#Hs{;tD1nK*%PXLCtUo|4MdzF z$B~Ph(Lf7Pa{bi^x1Ewl3uP|YFdhL zYlS3d`gh{CV0%OUH}kTl)o*=@zbW6crlj97Vc!XToa7JJe#)^!deQ$mu%gqe*l~ZQDTaTEAgqU0qABg8cFoQFhFGRjNt* z%^q~k_>#=Cu2UaAv-eZ@&LVh?lH^@OtdSsL^dFS){E=)3^}JrxPL#oUgr;YRX*X-y zhg9ZqZst7$pMB9Zg_bM5@RtELku=s_PcJ7;+bq^V+?*b-OXEJF?PEF8Z!$j3WR}@& zj&f5`m$YzOv53h?hgHzGOi83`Gt~SVt7u=@r_D~|3ItqP6oy~5YSB|z$y3A0`pv|> zHRc?(;gsWzJ+3;REiUY3qN~Jf<2bxs1^E`I3u5wS{oD%@^l#SbXr9vYz<(insW`|5 z_Unw9R-gEG;4Foj9%DOBM343fYYFCsQXz}4ei04nYh2lvi-Kw!(r0gYyl#~#F@b8e z3HYF8E8ho#p$Xx~l^IJ>w29@@V187N5ZJIU=%Tf zXAEKGrBgr1g{x$2^N7_ZS%kcud|6pxoi%Y`WJxjgmItQz`hYfWg2-6uOs+|sst<1szru2+?D349=GG?xozo`ChhpuXuJlyHbQ%d+ zOs-rdQDez&477-FOED>^i=bKWPl(B(qql~jw>(jMq$^dk0#*KJv%(lR&)&uvHUaTv z+e1NNasHjBn;;jJzE!3?D;Te=*6v4rH6Xeqlc@;i6Rj76g(YpQK+fNy)5P_%oHs%IYAj)aLH2pSR}rLheW0XPw-O_&w9C z-Pcbgdv*L&{vEY-W?_Za;n6<;7Zu!e0&?3s8g;P3g3g@hUAMltVcH_V%8YnO;9w7LwNx-&r%o@{q zI0uWu?>+Hwre5+FxbuxtE1gp<-M4#Gu$c=n<5TbB8kG%t9<8VWVb*-Qor3n&H>l`BR>Hu?$a- zr59UEB(;!pk!+=iNz2>Jq-&&Ehvc!8NEM~FN3jUwzr6w-cC&2#`Zxa`F~OOc3e!Zn zVE(c)$#7OUp+vq6Ckfi3LK;AjH8%>Cj5ZlB)}L4T{l*eM84uE%`R&e^P2RxQC-wY0 zN+p`|R`;tN*Oq5BO3Qs3y#Sul@5-f(5vAsEST`;g0~}iCDW~$>LK{(MH9y{8BqjRQ z2wKK`;UCoy<#L`oVAkwtPVwhg=N98@rf2jYm^_m=zhuID137Z7o1S^32FK~u@ah-d zgj4k;xZLw7JV@w1Ue}`dFgAoi;gLnizAchr@fL;vK{Mo&Y&{Y8&IBRem6N&kE#i8*b(sa!RcDPey_4M6OB zyRC{+xE%b(@=34E{4)-IkA`8K%zE-4hu?_56H1>4Uy^B~R@tj!(tT$T?}?YPk6w#>F^@p_#Op25Jh2kuFikjzsYfI- zef+W*x6HGwP{h;R`h}+|=fpsGEryKH^_?&@-P+xJt$$v*Pt|MV@gIeEJy?m19!xF% z-}W9L%{FmxlwzDGt=s4ZB8+}z9I+43o`JZ&4gsYMe5198%F(ETCxdaN4&B{C9Rt5etk>2KzGZQ2V^|&TROHjEF zl(&%6yk75o2~J}nc;=8DnzH}fp}z@n@EtbAh*o~mubY@RVN~%~NTqnl4d#}j6PqD! zyA7ND{3C42A8A`IJR;r2E+N6&P^75J>LmLTmQRtBBrp4XEDhhf_S#iH_GV**ODe;lnSUqry! zsyK1$1IESwbuM)qt4!@W+T^u-9siS_E9{!MVdDYSv6z3ORVT6}ykT?pMdyG#@~6&F z<;O5q+N-y@aPiY;IQ;nn5Z$icdiGStArcNJ@x=9~@CDm>eUK95DJ;G`5CT!RVQ!(ZM3$;UzCLlv0V+4Vn8O*jHzGm>$Xsad$V0N zH-}@n@`4?I4C^cHo%8tI>xfJ&-V^XtNexZ0zbNjV+nHC{#a8k9-%Eo`K>(yw3>riY zv8ohV0ZT*L)Z+)RG+tr%{Q5EAG53qve*t$mfKuzzVYW|T!8>TrqX7xm2Tv%1sR-Fx zO4?nq>B`)Q?Bi?&DnBfD{$YK=O9Y-Hq^8_wT6Ev`(->qk=S-dD%@Iwb+@`bI95*kW zP{wZ{BE@KBjsE~#od^%pmxlwEQ=TMh-snC@H5EDLHPY!+?szI6a z?fIPe*Dce{A7|zy8qn_(R%h7o??Nd5HPdV)IQf6Y_$E|27kYi)VDEmDn3k zs=f1jIoso#&8!N?gK*dwZke_!;LuO5YA+x$Z|I>3)EuvJc9DA8}XcNP(QdHI1H~(cq_wFIGD}wA8nla>M z9{1QQlHr%x6hP2TyEnN`e2Jk^60}BFC%1&4?ioD_SKVPEvB`@7f(}QS~1uDh=E)lI21{~Kyh4U z9xnWrxNXaxd85BV?{J5_nZ$wJu9Tc4)=Ii7TRj#CMx!=bi!XJ#UiMmll+zpXQjUB4 zvx?}Yu-}K5;o26pgE|2lC;$BN!&_rp4N513SGD0Yt5Y*~uUpnmuNf&@%c%K))|2x8 zN3-83-X`e~@YW0p1_|xHisSU0bt}(aceRF(&%2WFA731+c$H4bbo=8(GTutR%Vz8! zx2@9AB`=@d2hP~+Yb}c7FJ*%t;~tyuuPu+6qIwO+QY;Zn6yjzcR*rzk(zbH0jN3^C zIj+WwT_*xPNS#npPQg~OyOCxv@}l#)QL^65fSzfSGQf#Ts!rR#Z=N|S3AgILoxI0c zf|%>P3W#nK&1CQ4epkM~a|dsq*L zW3A?5-JP8g7B+@o$e{>k559d01c^De-kUj>G$c*^u@L78^Nx%r>5>Z#5c@GMX&^O` zNtZOkAo4cy%iBD#8Ow=P9@C>YjIcGK#47tIrgGvG;ojXYK0eOHYxOtKb#`JQYn)|* z%iFEyh1;pHaYM?Bo-OecbeWGnI#uUNwzICeW&_hRw=2;_U=|E^{Eon+MLn*VV|qJ~ z?9S}%`mHr6wSgySNqR^~s>j>g2!cH$W|`27K&AX0MUVE$!4>%?eQ#D>in|XU>)fw@ z0K;$~MwVUdzAYP}nzdmVf3**ilA=#sl5o~Y4#F26_ML5|GT!Pkm%|>`F4oIbxcYnX zD(tPjfKL_U7(80)UZ|Kq3#|@}u16$T_Kw>I>R7~EVNsy&5;%E;rVO||XXInt+u_kS z3}z)7?v=O|4lm|76FspxLd5ce7e|$XcDMb^TCzFzLq$%LVJhb#9RJJUKHm| zfAsNfe6?yXOnF?K*ZDDS*(N$(u5->cMpW&!=VaVjH!SNT;exAy5_|%L0FsZbZ&0=rUEgI2jwxl*~Lv!)bg#@7mcI-l!>xLgpo5gCL^5_1e>=uXj>7NlnGzo#y)QuRWSfKkF&l zGrCQGj$mx**^NNJQcvlwX0m7{|7agqKsg#BPG{-C+Y4w=ODeJW=e@I~x}RrSO4nuA zx`;dyiL!_@&{awMQpMu)R-u~CjIHPa4tUy}znq?DCIp&NmjN54*& zwzpXjPSs_ZmR$j8ETcdX+ppKRz6`hF0)RY__hz*1d+V1x$PYH;O~D_l&n)d^EIcs2 z$B{&`JoNWCY0bPu`6;QRh`o*CFuN=>zWAm15NK{WfO zJISU-(*YR_AB-AH>N4 zf1Q0kLD34N)-5j_>q{BPni|>U2WJuA7j*mVie06fljPSI0(L?o$&;TEa5Z38GVFJ= zx40n!7YBAJTV$bT*_KXb_e*I8ne-zv*^M}Ml;OX%r7MTOwYySykDst0HSYyX*^}{K z7w(hpB7S=&OA#UuocVUX5^Lp)2?(WuR-}NNtpvPE5qVO-2s0)-jX;k$s_KY5Ryf`J z6+cXW4u<_oG8~dX5-% zs)=L1qdBC67LRsgMyX4Kcr5&BD;W|qXIRPV@OX@0&f3%@M&#r?H!PH~T}ipSA7~6r z+~1v&p`chaW4orr_>FGT`t;`}#c070%ZMmaar++E^v{CZx~)V4+coV!DX1tLwS29o zmxHzZT(_1)CS(!6vCf0NZjkw-Up+pj-}uj-?{g$NG0DrI&lO3!Ha6lvtA|9_Y4R<4 z>?{fK)*4N#e_1Ov(>;lZD zP>JU;E5!ljn0bT0%>8Y1IJJL&Z+-^{yzcQm$@f!cOsRYN$3ZR0{tuk0dV&|j|1syP zjOI15I^}?9)mP-W38jMCmQ(UpbE*h76v_#OKBp*wQzbwJsh+sLM6Gr3y%tT96s(YB ztB)a^zDN-4+2XmHI`Z_iNNU;^VwLnKrnDS@Jk8kp%aaD8Xn8YfAu;KMm{Z3(rt6ua zFeNjPYf)X@*Ug@hWOFmN=urxtu**s>BY`Woy&$GG>spJHA7M<`kOrY~T(v208IO+l1d6Tloh$n+3V?Y3ZI>3h0N70iy zhWlSLpR7a-|KGGzo-Im9mPOC#p^lG%NX89Zv@t51#zbG$47{Q1WZ(Li`=m)g6#;mU zrNQm|Vrb1gYqCGR9#>d#!v_acfa%yy%5ZdtvbEaaXSmctkJd#SZD-!j2K{*M>6q=- zX!DhlZMT$YH^vQH6Y{wl+mEHxxRmy6l8aKDSs?=}JHYP>6>=mzMYdc8FN8MjuP+a{ zq>H|l$!#ccODA;G+lomrQQn&`Q2ypECLG8fEcxNB_4ttwHuYZkzGraYJTufq zM&T9@S*r@djcCb1?XRnsjwqbcqRd99*f~_ek8oOR-}@a1LOkt1w}xKspI@wL zj~y>FcK7nKgLrDQKuOhJ9@2i7-EYon6*;rK-FnIs4x{HJoGI(gmFSqjlAzIb$4u73 zRrn?=-F|FN%M(}0>5Kl88Bl%QlhRe*Zf}H(+w~vCp*Kb}D(pk2BdHJ2xh6wELMb?G05}j{7(?sDuhNRzrF$kE9 zP!FR*u9iK_F56O*)4Q>;hWBesNFTLOT(thP2lu@}W`Z(yU7^zPkx_Lx1m2ZDYW^8! z|2SSCW;vma-itFw^zg7O2VIAOVQl@GU;O2q=32+!63>U6SDn&cSZda$s;Xd@hJRu2 zxBSbujI(u9-q^CBnpeE0^LgOW)ySioz^sg-522%J!0r70gr?J2x^<;9&yX<*H?BHL z{&L-2m?6*Cm4Qxu0PVqN$V&TG8OJkFHP(Ch5LlNjmpNOq19Z3%OV_>1QCiZJdC9l+ zb$4Z~CTMfX-a1>Gy3axka$mH5UxM@>^hxh+y>!o$_S9mNpNk2zN@&5?*_X(w`dJoR?f@%+Hm>hS&^8ITe9V{kOj{yCC^Gt>XA+ig|JCLP_;fHiD7} zG$a-uMjPv9@Kq^DIOg1}W`|AXLn;i=oZl%eT_= z#dm_EGX6`}iIxLDy*Ot57}RfX-chFQJwQ}H%lj>B_OdnS#Li_8!8pYAnj0em~Rwm%73tg?1ef)6r30e@Sq)T*6=byM> z_H)tbI^f@Q9n&;6=vLpwk`g>M&GD@i!{PAX_2-1}@Up0mADHG6rJaslk9Lz9K4h1Q z(`q?n%>PBAqF0WLudivrGZwS(O}cV34hxrzmYcD>3p`&el#f17N}RgOMS=OntxdVU zCU6huMju|U(r(tlKc~FMYl`jNw7LQ5{nKUySPh&H46v8zheZh>Z+>*wzXK>LghYBg zi$0dW?O6t+ugVC}Udi%d!0ZntHJWmdM z%m|eq_u8a-%!n@BJ6qGJW!}44?cdW0D_m8d_L&L|Uuq)NN!W@z_-%kv{1%(SkSu^1 zJH8J`Q?cs+-P0n%suhXriqLjYFpAQ%;Wt*Ic!Z~3r26D%NKU&?FIL3=_SO*Zt8{=2A`8BEo8{n(hrKxiJ$MtXUqf$E!l5r9nw==9;z)F2r z?EF~yNGa#R5GK_O9iqIDdQ2nNO_+gWv~Mrm%&%}G=FKxZ_E?l&%73lET#QtuJAKY@ z6K88B%saF z#0oq*%hzdqy&RZ*_x`DHGWgkvPYbw})2V3qvBEGx{2ArlfnhSXDE$qj7V37sxn%K0 zq9cC%8SdRD3PzB@rB8f*fFWj!T25pNQT5kN6TmENQDR>p51X1xi_e>w17_+IRI?Op zDWh&eFp>$5LyAb8Oic$4)rIAVKHxqgzx$}V@B|F8j(_Uim<>%gXzv$!2?NM`kwc5Z zMccm)n6upZ3cmnEtqE+W!uNA0NZP6E-wJ(i307gL7wBpJ;>3JJGU8$QRDqb1GkMeK z?dP(zZ))GHk}Ekce6(6=j$AJn6TBW7le#2cM&YB)@otOYMVPQZB`ZHKP<|makMM{VU&8z1aQP$rW^H zwfIKuN@$GDlVb>Buj_63#+7ULv|}wd)A`?C<@7$Q>$tdW###%*##(!c0?WbZf5o?r zEa2e6>%E>oC5X~oB%TV7WE{{33OK4RwM7qd;$2NwY*HqrzMegv`%K*_y+we?N4V zjSSRzV*!MXr{qSD#T|)52d2;mtv^X`^zBCWuR_B(gMx4{zr$xL(b{xPtKqP${ltr) zjb$j0fL_LG1{F@9_@wZzHvDO0^O-o-5#eP3A7VuO^-A9~Cvn-z&q6jaV$w|89VGcLtwT{lxv{b7~Q z_)o(t^mdEB(OF#;tZ0q&&ZYonKiJdQOVA%)b|7`V?=E%VMUEw?Fl$Lcx59MF8pw#? zJhfqake#DiNN5(f8~!fR!v9@B#Ko89D0nVDT;%Ofxz4BPL_a6*hlfOUb#_lg?~ymi>c_TR}jM_xDehx0i8bfyP82M8Hk?OFO7_7dvCCk%QB9H_Ix#r zlkpb`Tr72WkdHscKAep8^|*w7Vm{uIl}+!2w({@3y1hC3#Ef!ft(zjR-apzRHQ-5c z-%pTBw0O056m$BrH;&{J9r5NcWPnj5+4YKu_niJp;Ge)HUYoy`ylL+rw6>{yRcm4* zzk$>5Y6@p=a4Rdq*+2U0)*$KD`5a;Ulu+fzhf`D3f6vb1V80XMxNGN0tHp+ej&%G# z-KMU@C}XTM_5Ym{GWo877wGixWUn_qe)s$EP2QMo2$ySsyqVRQ>#Oo;&WzKs5&2O} zZ!mqEUHLO!e-*?LRGU65$1L)OJUSNfPbl_E=e*3o!}oa^_p{MFo$vh}@5eiR-=x#* ze2jf%mVUnkPf^l<)R3-6++4}EiFJHzS;flyA5+rdVCc!gmxCiThmRMcPTre)b!Ji@ zAn6ybW8^k>H)49X7ng`5<9)99hHHk{GLj6&maA1eERqcAxjj=V=BE8GEEybb1Z>p1 zJTGr+o98Z-ZKg!A9Kd;Tud|}hCP`|K>z&Yw#nET)xyz!{$y44@H<2i-{0T<|Ot`jf z9HC2_Q!fw3ilaJCd~{%~{%R>Yhz^sdvDcvauPNSE%O3~q4)8ej_w(5F=e6oS+MWv7 zY1*acs7UKkmuw0YHJtg{l`qi{{Zirb>D8|I;2&g($>8*2WzYMCH(Bu)xMhIrQ@*yw H2lxILHP*tE literal 0 HcmV?d00001 diff --git a/inst/Examples/src/bestdose_past.csv b/inst/Examples/src/bestdose_past.csv new file mode 100644 index 000000000..f31bbc5e9 --- /dev/null +++ b/inst/Examples/src/bestdose_past.csv @@ -0,0 +1,9 @@ +"id","evid","time","dur","dose","addl","ii","input","out","outeq","c0","c1","c2","c3" +1,1,0,0,150,0,0,1,.,.,.,.,.,. +1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,. +1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,. +1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,. +1,1,12,0,75,0,0,1,.,.,.,.,.,. +1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,. +1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,. +1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,. diff --git a/inst/Examples/src/bestdose_prior.csv b/inst/Examples/src/bestdose_prior.csv new file mode 100644 index 000000000..0559dea12 --- /dev/null +++ b/inst/Examples/src/bestdose_prior.csv @@ -0,0 +1,47 @@ +ke,v,prob +0.08736658442020416,104.1576635837555,0.06411826509758521 +0.1305751524925232,97.43967413902283,0.03921568628597478 +0.3540655417442322,86.90321147441864,0.03930537754808413 +0.2908282660484314,113.02810430526733,0.05369391213529181 +0.10566180405616761,140.3564077615738,0.01960784315457728 +0.9795492498397828,221.82963848114014,0.019607843139348134 +0.3214192095041275,68.54407012462616,0.019607843027666164 +0.04363290762901306,85.5460512638092,0.01960721617788871 +0.3493796042442322,75.69715678691864,0.019607073213094923 +0.09591740503311158,70.4919171333313,0.019606975069368423 +0.3318302191734314,92.54958868026733,0.019607488410309255 +0.3091094713449478,92.16600060462952,0.01969419594801783 +0.3319534166574478,118.26951622962952,0.019631823530513178 +0.06250557525157929,107.17031538486481,0.019670607599452994 +0.06885235497951508,76.57066643238068,0.01970716827239974 +0.104416601395607,91.60142242908478,0.019644091781941257 +0.11730292952060699,78.54966461658478,0.019597452684180432 +0.06250557525157929,95.70058882236481,0.019558918237699178 +0.013994081234931946,178.82485330104828,0.019739518261060587 +0.3406021231412888,99.99475717544556,0.019495129816772868 +0.07190197534561157,67.0641827583313,0.019509009848450342 +0.02050767450332642,140.61931788921356,0.019494247900270656 +0.283145404958725,125.85132718086243,0.021689589949674865 +0.28171865940093993,140.98521947860718,0.020065175979272123 +0.2902425238609314,94.48318243026733,0.02049255770567427 +0.09111055810451509,93.09410393238068,0.019720100726323925 +0.2925854926109314,128.45290899276733,0.019127402830562693 +0.2925854926109314,105.60134649276733,0.02235654024528764 +0.3019573676109314,127.61794805526733,0.01954540376390492 +0.044804392004013066,99.6524965763092,0.021082697423336363 +0.28490263152122497,92.62867093086243,0.018708872021706954 +0.08112040762901307,98.0704653263092,0.021075404370655606 +0.036757444095611574,109.3835186958313,0.018114822039674315 +0.09826037378311157,111.2731671333313,0.021727337750962878 +0.09298869409561158,113.6022686958313,0.01713443002609149 +0.3132096666574478,115.89646935462952,0.00474084568678399 +0.07658791284561157,113.3385968208313,0.002945117861059519 +0.08678084223270416,87.5463354587555,0.02933991222454541 +0.08736658442020416,128.0199682712555,0.019952065962974042 +0.05615650746822357,130.38859486579895,0.000015108929159013186 +0.056742249655723574,131.13566517829895,0.012605671648896763 +0.08736658442020416,87.5463354587555,0.007152585492468555 +0.07658791284561157,113.3825421333313,0.03210569535127533 +0.3132096666574478,115.94041466712952,0.03478763091025865 +0.05615650746822357,130.43254017829895,0.04616814785864839 +0.08736658442020416,127.9760229587555,0.00001919610085464287 diff --git a/inst/Examples/src/bestdose_target.csv b/inst/Examples/src/bestdose_target.csv new file mode 100644 index 000000000..29d8788a6 --- /dev/null +++ b/inst/Examples/src/bestdose_target.csv @@ -0,0 +1,9 @@ +"id","evid","time","dur","dose","addl","ii","input","out","outeq","c0","c1","c2","c3" +1,1,0,0,0,0,0,1,.,.,.,.,.,. +1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,. +1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,. +1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,. +1,1,12,0,0,0,0,1,.,.,.,.,.,. +1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,. +1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,. +1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,. From 2236e6cad358066ae1e6a08abf91b374fc465b9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Tue, 18 Nov 2025 14:24:20 +0000 Subject: [PATCH 5/5] bestdose is subdivided in two steps, posterior calculation and optimization --- Cargo.lock | 10 +- NAMESPACE | 1 + R/PM_bestdose.R | 508 ++++++++++--------- R/extendr-wrappers.R | 4 + inst/.gitignore | 2 + inst/Examples/Rscript/bestdose_simple_test.R | 19 +- man/makeNCA.Rd | 2 +- man/plot.PM_cov.Rd | 2 +- man/plot.PM_cycle.Rd | 2 +- man/plot.PM_final.Rd | 2 +- man/plot.PM_op.Rd | 2 +- man/plot.PM_pta.Rd | 12 +- man/plot.PM_sim.Rd | 2 +- src/rust/Cargo.toml | 4 +- src/rust/src/bestdose_executor.rs | 345 ++++++++++--- src/rust/src/executor.rs | 2 +- src/rust/src/lib.rs | 40 ++ 17 files changed, 626 insertions(+), 333 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 0a6d1402e..a4064bec1 100755 --- a/Cargo.lock +++ b/Cargo.lock @@ -1286,9 +1286,9 @@ dependencies = [ [[package]] name = "pharmsol" -version = "0.20.0" +version = "0.21.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "98b8e2ab3a0e91cd4b20c28544cb3676e8df31aa490cf5680ec0531259b5fa4e" +checksum = "2fc25564d039d0cd5701013aa3785a339b14cf0b51409d7b817320bc360dc944" dependencies = [ "argmin", "argmin-math", @@ -1320,6 +1320,7 @@ version = "0.1.0" dependencies = [ "anyhow", "extendr-api", + "libloading", "pmcore", "rayon", "tracing", @@ -1328,13 +1329,12 @@ dependencies = [ [[package]] name = "pmcore" -version = "0.21.0" +version = "0.22.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "74646008c207f79cdfd54152cefae76ebcba895fd6005b23de42dfb838a2ca95" +checksum = "3866100507aa3bcba475381af3102d84b5e503bce82cb82f56cf0fa46cc1e408" dependencies = [ "anyhow", "argmin", - "argmin-math", "csv", "faer", "faer-ext", diff --git a/NAMESPACE b/NAMESPACE index 99e1132dd..63feb0401 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -48,6 +48,7 @@ export(NPrun) export(PMFortranConfig) export(PM_batch) export(PM_bestdose) +export(PM_bestdose_problem) export(PM_build) export(PM_compare) export(PM_cov) diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R index e455b9cd2..684c04e02 100644 --- a/R/PM_bestdose.R +++ b/R/PM_bestdose.R @@ -1,3 +1,103 @@ +bestdose_parse_prior <- function(prior) { + if (inherits(prior, "PM_result")) { + theta_path <- file.path(prior$rundir, "outputs", "theta.csv") + if (!file.exists(theta_path)) { + cli::cli_abort("theta.csv not found in PM_result outputs") + } + theta_path + } else if (inherits(prior, "PM_final")) { + temp_path <- tempfile(fileext = ".csv") + bestdose_write_prior_csv(prior, temp_path) + temp_path + } else if (is.character(prior)) { + if (!file.exists(prior)) { + cli::cli_abort("Prior file not found: {prior}") + } + prior + } else { + cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv") + } +} + +bestdose_write_prior_csv <- function(prior, path) { + df <- as.data.frame(prior$popPoints) + df$prob <- prior$popProb + write.csv(df, path, row.names = FALSE, quote = FALSE) +} + +bestdose_parse_model <- function(model) { + if (inherits(model, "PM_model")) { + compiled_path <- model$binary_path + if (is.null(compiled_path) || !file.exists(compiled_path)) { + cli::cli_abort("Model must be compiled first. Use model$compile()") + } + + kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) { + "analytical" + } else { + "ode" + } + + list(path = compiled_path, kind = kind, model = model) + } else if (is.character(model)) { + if (!file.exists(model)) { + cli::cli_abort("Model file not found: {model}") + } + kind <- if (grepl("analytical", model, ignore.case = TRUE)) { + "analytical" + } else { + "ode" + } + list(path = model, kind = kind, model = NULL) + } else { + cli::cli_abort("model must be PM_model or path to compiled model") + } +} + +bestdose_parse_data <- function(data) { + if (inherits(data, "PM_data")) { + temp_path <- tempfile(fileext = ".csv") + write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE) + temp_path + } else if (is.character(data)) { + if (!file.exists(data)) { + cli::cli_abort("Data file not found: {data}") + } + data + } else { + cli::cli_abort("data must be PM_data or path to CSV file") + } +} + +bestdose_default_settings <- function(prior, model) { + if (inherits(prior, "PM_result")) { + return(prior$settings) + } + + param_ranges <- lapply(model$model_list$pri, function(x) { + c(x$min, x$max) + }) + names(param_ranges) <- tolower(names(param_ranges)) + + list( + algorithm = "NPAG", + ranges = param_ranges, + error_models = list( + list( + initial = 0.0, + type = "additive", + coeff = c(0.0, 0.2, 0.0, 0.0) + ) + ), + max_cycles = 500, + points = 2028, + seed = 22, + prior = "prior.csv", + idelta = 0.25, + tad = 0.0 + ) +} + #' @title #' Object to contain BestDose optimization results #' @@ -12,152 +112,44 @@ PM_bestdose <- R6::R6Class( "PM_bestdose", public = list( - #' @field result A list containing: - #' * **doses** Vector of optimal dose amounts (mg) - #' * **objf** Objective function value (lower is better) - #' * **status** Optimization status ("converged" or "max_iterations") - #' * **predictions** Data frame with concentration-time predictions - #' * **auc_predictions** Data frame with AUC predictions (if applicable) - #' * **method** Optimization method used ("posterior" or "uniform") result = NULL, - - #' @description - #' Create a new PM_bestdose object by running BestDose optimization - #' - #' @param prior A PM_result or PM_final object, or path to theta.csv file - #' @param model A PM_model object or path to compiled model - #' @param past_data Optional: PM_data object or path to CSV with patient history - #' @param target PM_data object or path to CSV with target doses/observations - #' @param dose_range Named list with min and max allowable doses, e.g., - #' list(min = 0, max = 1000). Default: list(min = 0, max = 1000) - #' @param bias_weight Numeric [0,1]. 0 = fully personalized, 1 = population-based. - #' Default = 0.5 (balanced) - #' @param target_type One of "concentration", "auc_from_zero", or "auc_from_last_dose". - #' Default = "concentration" - #' @param time_offset Optional: time offset for past/future concatenation - #' @param settings Optional: list of settings (algorithm, error models, etc.) - #' If NULL, will use defaults from prior object - #' - #' @details - #' - #' ## Target Data Format - #' - #' The target data should be a Pmetrics-formatted CSV with: - #' - Dose events: Set dose amount to 0 for doses to be optimized, - #' or non-zero for fixed doses - #' - Observation events: Set OUT values to the target concentrations or AUCs - #' - #' ## Target Types - #' - #' - **concentration**: Optimize to match target concentrations at observation times - #' - **auc_from_zero**: Optimize to match cumulative AUC from time 0 - #' - **auc_from_last_dose**: Optimize to match interval AUC from last dose - #' - #' ## Bias Weight - #' - #' Controls the balance between patient-specific and population-based optimization: - #' - 0.0: Fully personalized (minimize variance) - #' - 0.5: Balanced approach (recommended default) - #' - 1.0: Population-based (minimize bias from population) - #' - #' @examples - #' \dontrun{ - #' # Load NPAG result - #' result <- PM_load(1) - #' - #' # Create target data - #' target <- PM_data$new("target.csv") - #' - #' # Run BestDose optimization - #' bd <- PM_bestdose$new( - #' prior = result, - #' model = result$model, - #' target = target, - #' dose_range = list(min = 50, max = 500), - #' bias_weight = 0.5 - #' ) - #' - #' # View results - #' print(bd) - #' bd$get_doses() - #' bd$get_predictions() - #' } - initialize = function(prior, - model, + problem = NULL, + bias_weight = NULL, + initialize = function(prior = NULL, + model = NULL, past_data = NULL, - target, + target = NULL, dose_range = list(min = 0, max = 1000), bias_weight = 0.5, target_type = "concentration", time_offset = NULL, - settings = NULL) { - # Validate inputs - if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) { - cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose") + settings = NULL, + result = NULL, + problem_obj = NULL, + bias_override = NULL) { + if (!is.null(result)) { + private$.set_result(result, problem_obj, bias_override) + return(invisible(self)) } - if (bias_weight < 0 || bias_weight > 1) { - cli::cli_abort("bias_weight must be between 0 and 1") - } - - if (is.null(dose_range$min) || is.null(dose_range$max)) { - cli::cli_abort("dose_range must have both 'min' and 'max' elements") - } - - if (dose_range$min >= dose_range$max) { - cli::cli_abort("dose_range$min must be less than dose_range$max") + if (is.null(target)) { + cli::cli_abort("target must be supplied when computing a new BestDose result") } - # Parse prior - prior_path <- private$.parse_prior(prior) - - # Parse model - model_info <- private$.parse_model(model) - - # Parse data files - past_data_path <- if (!is.null(past_data)) { - private$.parse_data(past_data) - } else { - NULL - } - - target_data_path <- private$.parse_data(target) - - # Parse settings - if (is.null(settings)) { - settings <- private$.default_bestdose_settings(prior, model) - } - - # Call Rust function - cli::cli_alert_info("Running BestDose optimization...") - - result <- tryCatch( - { - bestdose( - model_path = model_info$path, - prior_path = prior_path, - past_data_path = past_data_path, - target_data_path = target_data_path, - time_offset = time_offset, - dose_min = dose_range$min, - dose_max = dose_range$max, - bias_weight = bias_weight, - target_type = target_type, - params = settings, - kind = model_info$kind - ) - }, - error = function(e) { - cli::cli_abort(c( - "x" = "BestDose optimization failed", - "i" = conditionMessage(e) - )) - } + problem <- PM_bestdose_problem$new( + prior = prior, + model = model, + past_data = past_data, + target = target, + dose_range = dose_range, + bias_weight = bias_weight, + target_type = target_type, + time_offset = time_offset, + settings = settings ) - self$result <- result - - cli::cli_alert_success("Optimization complete!") + raw <- problem$optimize_raw(bias_weight = bias_weight) + private$.set_result(raw$result, problem, raw$bias_weight) invisible(self) }, @@ -171,6 +163,9 @@ PM_bestdose <- R6::R6Class( cat(sprintf("ln(Objective): %.4f\n", log(self$get_objf()))) cat(sprintf("Method: %s\n", self$get_method())) cat(sprintf("Status: %s\n", self$get_status())) + if (!is.null(self$bias_weight)) { + cat(sprintf("Bias weight (lambda): %.2f\n", self$bias_weight)) + } cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions))) if (!is.null(self$result$auc_predictions)) { cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions))) @@ -230,126 +225,149 @@ PM_bestdose <- R6::R6Class( } ), private = list( - # Helper function to parse prior - .parse_prior = function(prior) { - if (inherits(prior, "PM_result")) { - # Extract theta.csv from PM_result - theta_path <- file.path(prior$rundir, "outputs", "theta.csv") - if (!file.exists(theta_path)) { - cli::cli_abort("theta.csv not found in PM_result outputs") - } - return(theta_path) - } else if (inherits(prior, "PM_final")) { - # Need to write out PM_final to CSV - temp_path <- tempfile(fileext = ".csv") - private$.write_prior_csv(prior, temp_path) - return(temp_path) - } else if (is.character(prior)) { - if (!file.exists(prior)) { - cli::cli_abort("Prior file not found: {prior}") - } - return(prior) - } else { - cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv") + .set_result = function(result, problem, bias_weight) { + self$result <- result + self$problem <- problem + self$bias_weight <- bias_weight + } + ) +) + +#' @title +#' Prepare a reusable BestDose optimization problem +#' +#' @description +#' `r lifecycle::badge("experimental")` +#' +#' Use `PM_bestdose_problem` to mirror the Rust workflow: compute posterior +#' support points once, inspect them in R, and solve for multiple bias weights +#' without repeating the expensive initialization step. +#' +#' @export +PM_bestdose_problem <- R6::R6Class( + "PM_bestdose_problem", + public = list( + handle = NULL, + theta = NULL, + theta_dim = NULL, + param_names = NULL, + posterior_weights = NULL, + population_weights = NULL, + bias_weight = NULL, + target_type = NULL, + dose_range = NULL, + model_info = NULL, + settings = NULL, + initialize = function(prior, + model, + past_data = NULL, + target, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + settings = NULL) { + if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) { + cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose") } - }, - # Helper function to parse model - .parse_model = function(model) { - if (inherits(model, "PM_model")) { - # Check if model is compiled - compiled_path <- model$binary_path - if (is.null(compiled_path) || !file.exists(compiled_path)) { - cli::cli_abort("Model must be compiled first. Use model$compile()") - } - - # Determine model type from model_list - kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) { - "analytical" - } else { - "ode" - } - - return(list( - path = compiled_path, - kind = kind - )) - } else if (is.character(model)) { - # Assume it's a path to compiled model - if (!file.exists(model)) { - cli::cli_abort("Model file not found: {model}") - } - # Try to infer type from extension or filename - kind <- if (grepl("analytical", model, ignore.case = TRUE)) { - "analytical" - } else { - "ode" - } - return(list(path = model, kind = kind)) - } else { - cli::cli_abort("model must be PM_model or path to compiled model") + if (bias_weight < 0 || bias_weight > 1) { + cli::cli_abort("bias_weight must be between 0 and 1") } - }, - # Helper function to parse data - .parse_data = function(data) { - if (inherits(data, "PM_data")) { - # Write to temp CSV - temp_path <- tempfile(fileext = ".csv") - write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE) - return(temp_path) - } else if (is.character(data)) { - if (!file.exists(data)) { - cli::cli_abort("Data file not found: {data}") - } - return(data) - } else { - cli::cli_abort("data must be PM_data or path to CSV file") + if (is.null(dose_range$min) || is.null(dose_range$max)) { + cli::cli_abort("dose_range must have both 'min' and 'max' elements") } - }, - # Helper to create default settings - .default_bestdose_settings = function(prior, model) { - # Extract settings from prior if it's a PM_result - if (inherits(prior, "PM_result")) { - # Use the settings from the result - settings <- prior$settings - return(settings) - } else { - # Get parameter ranges from model (same format as in fit() method) - param_ranges <- lapply(model$model_list$pri, function(x) { - c(x$min, x$max) - }) - names(param_ranges) <- tolower(names(param_ranges)) - - # Use sensible defaults - settings <- list( - algorithm = "NPAG", - ranges = param_ranges, - error_models = list( - list( - initial = 0.0, - type = "additive", - coeff = c(0.0, 0.2, 0.0, 0.0) - ) - ), - max_cycles = 500, - points = 2028, - seed = 22, - prior = "prior.csv", - idelta = 0.25, # 15 minutes for AUC calculations - tad = 0.0 - ) + if (dose_range$min >= dose_range$max) { + cli::cli_abort("dose_range$min must be less than dose_range$max") + } + + prior_path <- bestdose_parse_prior(prior) + model_info <- bestdose_parse_model(model) + past_data_path <- if (!is.null(past_data)) bestdose_parse_data(past_data) else NULL + target_data_path <- bestdose_parse_data(target) + + if (is.null(settings)) { + model_for_settings <- if (!is.null(model_info$model)) model_info$model else model + settings <- bestdose_default_settings(prior, model_for_settings) } - return(settings) + + prep <- bestdose_prepare( + model_path = model_info$path, + prior_path = prior_path, + past_data_path = past_data_path, + target_data_path = target_data_path, + time_offset = time_offset, + dose_min = dose_range$min, + dose_max = dose_range$max, + bias_weight = bias_weight, + target_type = target_type, + params = settings, + kind = model_info$kind + ) + + if (is.character(prep)) { + cli::cli_abort(prep) + } + + dim <- as.integer(prep$theta_dim) + theta_matrix <- matrix(prep$theta_values, nrow = dim[1], ncol = dim[2]) + colnames(theta_matrix) <- prep$param_names + + self$handle <- prep$handle + self$theta <- theta_matrix + self$theta_dim <- dim + self$param_names <- prep$param_names + self$posterior_weights <- prep$posterior_weights + self$population_weights <- prep$population_weights + self$bias_weight <- prep$bias_weight + self$target_type <- prep$target_type + self$dose_range <- dose_range + self$model_info <- model_info + self$settings <- settings + + cli::cli_alert_success("BestDose problem prepared with %d support points", dim[1]) }, + finalize = function() { + self$handle <- NULL + }, + + #' @description + #' Run optimization and return raw list (doses, objf, predictions) + optimize_raw = function(bias_weight = NULL) { + private$.run_optimize(bias_weight) + }, + + #' @description + #' Run optimization and return a `PM_bestdose` result object + optimize = function(bias_weight = NULL) { + raw <- self$optimize_raw(bias_weight) + PM_bestdose$new( + result = raw$result, + problem_obj = self, + bias_override = raw$bias_weight + ) + } + ), + private = list( + .run_optimize = function(bias_weight) { + if (is.null(self$handle)) { + cli::cli_abort("BestDose problem handle has been released") + } + + bw <- if (is.null(bias_weight)) self$bias_weight else bias_weight + + if (bw < 0 || bw > 1) { + cli::cli_abort("bias_weight must be between 0 and 1") + } + + res <- bestdose_optimize(self$handle, bw) + if (is.character(res)) { + cli::cli_abort(res) + } - # Helper to write PM_final to CSV - .write_prior_csv = function(prior, path) { - # Convert PM_final to CSV with prob column - df <- as.data.frame(prior$popPoints) - df$prob <- prior$popProb - write.csv(df, path, row.names = FALSE, quote = FALSE) + list(result = res, bias_weight = bw) } ) ) diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R index c65e13970..e199aec1d 100755 --- a/R/extendr-wrappers.R +++ b/R/extendr-wrappers.R @@ -44,5 +44,9 @@ setup_logs <- function() .Call(wrap__setup_logs) #'@export bestdose <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) +bestdose_prepare <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) + +bestdose_optimize <- function(handle, bias_weight) .Call(wrap__bestdose_optimize, handle, bias_weight) + # nolint end diff --git a/inst/.gitignore b/inst/.gitignore index 28f483682..4472f674a 100755 --- a/inst/.gitignore +++ b/inst/.gitignore @@ -1,3 +1,5 @@ .DS_Store /.quarto/ + +template/ diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R index 032d36c51..7cab124fa 100644 --- a/inst/Examples/Rscript/bestdose_simple_test.R +++ b/inst/Examples/Rscript/bestdose_simple_test.R @@ -24,14 +24,27 @@ target_file <- "../src/bestdose_target.csv" prior_file <- "../src/bestdose_prior.csv" -bd <- PM_bestdose$new( +# Prepare the problem once (posterior + handle to optimized model) +problem <- PM_bestdose_problem$new( prior = prior_file, model = mod_onecomp, past_data = past_file, target = target_file, dose_range = list(min = 0, max = 300), bias_weight = 0.0, - target_type = "concentration" + target_type = "concentration" # "concentration", "auc_from_zero", "auc_from_last_dose" ) -bd$print() +cat("\nPosterior support points:\n") +print(head(problem$theta)) + +# Reuse the same problem for different bias weights +bias_weights <- seq(0, 1, by = 0.25) +results <- lapply(bias_weights, function(lambda) { + problem$optimize(bias_weight = lambda) +}) + +for (i in seq_along(results)) { + cat("\n=== Bias weight:", bias_weights[i], "===\n") + results[[i]]$print() +} diff --git a/man/makeNCA.Rd b/man/makeNCA.Rd index 4dc0ffcba..737122886 100755 --- a/man/makeNCA.Rd +++ b/man/makeNCA.Rd @@ -25,7 +25,7 @@ makeNCA( \itemize{ \item It can be the run number, e.g. 3, that has been previously loaded with \link{PM_load}. Either the mdata file from the run (NPAG or IT2B) can be used (default) or -the post object can be used (NPAG only) by specifying \code{postPred = T} below. If \link{x} is a run number that corresponds +the post object can be used (NPAG only) by specifying \code{postPred = T} below. If \link[ggplot2:aes_position]{ggplot2::x} is a run number that corresponds to both an NPAG and IT2B run which have been previously loaded into memory with \link{PM_load}, the NPAG run will be used. \item It can be the run number of a run that has \emph{not} been previously loaded with \link{PM_load}. In this case, the current working directory should be the Runs folder as \link{makeNCA} will call \link{PM_load}. diff --git a/man/plot.PM_cov.Rd b/man/plot.PM_cov.Rd index ea54d695a..fedc8dc5b 100755 --- a/man/plot.PM_cov.Rd +++ b/man/plot.PM_cov.Rd @@ -171,7 +171,7 @@ NPex$cov$plot(V ~ wt, } } \seealso{ -\link{PM_cov}, \link{PM_result}, \link{schema} +\link{PM_cov}, \link{PM_result}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cycle}()}, diff --git a/man/plot.PM_cycle.Rd b/man/plot.PM_cycle.Rd index 67acb7e6b..f91068602 100755 --- a/man/plot.PM_cycle.Rd +++ b/man/plot.PM_cycle.Rd @@ -114,7 +114,7 @@ NPex$cycle$plot( } } \seealso{ -\link{PM_result}, \link{schema} +\link{PM_result}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/man/plot.PM_final.Rd b/man/plot.PM_final.Rd index 9e2374ae1..5b9129fab 100755 --- a/man/plot.PM_final.Rd +++ b/man/plot.PM_final.Rd @@ -161,7 +161,7 @@ ITex$final$plot(Ke ~ V) } } \seealso{ -\link{PM_final}, \link{schema} +\link{PM_final}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/man/plot.PM_op.Rd b/man/plot.PM_op.Rd index d3696144b..ed683317a 100755 --- a/man/plot.PM_op.Rd +++ b/man/plot.PM_op.Rd @@ -187,7 +187,7 @@ NPex$op$plot(stats = list(x = 0.5, y = 0.2, font = list(size = 7, color = "blue" } } \seealso{ -\link{PM_result}, \link{PM_op}, \link{schema} +\link{PM_result}, \link{PM_op}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/man/plot.PM_pta.Rd b/man/plot.PM_pta.Rd index 833e52de3..06168a61f 100755 --- a/man/plot.PM_pta.Rd +++ b/man/plot.PM_pta.Rd @@ -61,7 +61,7 @@ the \code{line} argument should be a list of the options for group based plottin where each group corresponds to a simulated regimen. The possible elements of the \code{line} list should be exactly named: \itemize{ -\item color Maps to the \link{plot_ly} \code{colors} argument to override default colors +\item color Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{colors} argument to override default colors applied to the lines for each regimen. This can be a named palette, which can be obtained with \code{RColorBrewer::display.brewer.all()} or a vector of hexadecimal color names. One way to ensure reliable color palettes is to use the @@ -72,12 +72,12 @@ of JavaScript on the ColorBrewer website. The default is "Set1". Palettes with fewer colors than regimens will be recycled. A color can also be a character vector of color names, recycled as needed. For example, a print-friendly choice is \code{line = list(color = "black")}. -\item width Maps to the \link{plot_ly} \code{width} argument to override default widths +\item width Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{width} argument to override default widths applied to the lines for each regimen. All lines will have the same width. The default value is 2. -\item dash Maps to the \link{plot_ly} \code{linetypes} argument to override default styles +\item dash Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{linetypes} argument to override default styles applied to the lines for each regimen. If numeric, will map to \code{lty} \link{par} values. -It can also be a character vector of dash names as listed in \link{plot_ly}. +It can also be a character vector of dash names as listed in \link[plotly:plot_ly]{plotly::plot_ly}. Example: \code{line = list(color = "Blues", width = 1, dash = 2)}, whicb will result in dotted lines (dash = 2) all with width 1 but in different shades of blue. }} @@ -97,12 +97,12 @@ where each group corresponds to a simulated regimen. The possible elements of th marker color does not need to also be specified. Even if line plotting is suppressed with \code{line = F}, the default color value of "Set1" will be applied to markers, unless specified, e.g. \code{marker = list(color = "Blues")}. -\item symbol Maps to the \link{plot_ly} \code{symbols} argument to override default symbols +\item symbol Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{symbols} argument to override default symbols applied to the markers for each regimen. If only one value is supplied for this, it will be recycled for each regimen, i.e. all will have the same symbol. See \code{plotly::schema()}, traces > scatter > attributes > marker > symbol > values for options. -\item size Maps to the \link{plot_ly} \code{size} argument to override default size +\item size Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{size} argument to override default size applied to the markers for each regimen. All markers will have the same size. The default value is 12. }} diff --git a/man/plot.PM_sim.Rd b/man/plot.PM_sim.Rd index 2cf39f45c..7278bdb05 100755 --- a/man/plot.PM_sim.Rd +++ b/man/plot.PM_sim.Rd @@ -191,7 +191,7 @@ simEx$plot(log = FALSE, line = list(color = "orange")) } } \seealso{ -\link{PM_sim}, \link{plot_ly}, \link{schema} +\link{PM_sim}, \link[plotly:plot_ly]{plotly::plot_ly}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml index bca68233b..d8fe6025a 100755 --- a/src/rust/Cargo.toml +++ b/src/rust/Cargo.toml @@ -9,8 +9,8 @@ name = 'pm_rs' [dependencies] extendr-api = '*' -pmcore = {version ="=0.21.0", features = ["exa"]} -# pmcore = { path = "/Users/jotalvaro/code/LAPKB/PMcore", features = ["exa"] } +pmcore = { version = "=0.22.1", features = ["exa"] } +libloading = "0.8" rayon = "1.10.0" anyhow = "1.0.97" diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs index 8caee71ed..2d83731f6 100644 --- a/src/rust/src/bestdose_executor.rs +++ b/src/rust/src/bestdose_executor.rs @@ -1,4 +1,4 @@ -use crate::settings::settings; +use crate::{logs::RFormatLayer, settings::settings}; use extendr_api::prelude::*; use pmcore::bestdose::{BestDoseProblem, BestDoseResult, DoseRange, Target}; use pmcore::prelude::{data, ODE}; @@ -112,8 +112,102 @@ pub(crate) fn convert_bestdose_result_to_r( Ok(output.into()) } +/// Opaque handle that keeps the dynamic model library alive while reusing the +/// prepared `BestDoseProblem` for multiple optimization runs. +pub struct BestDoseProblemHandle { + problem: BestDoseProblem, + #[allow(dead_code)] + library: libloading::Library, +} + +impl BestDoseProblemHandle { + #[allow(clippy::too_many_arguments)] + pub fn new( + model_path: PathBuf, + prior_path: PathBuf, + past_data_path: Option, + target_data_path: PathBuf, + time_offset: Option, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + ) -> std::result::Result { + let (library, (eq, meta)) = + unsafe { pmcore::prelude::pharmsol::exa::load::load::(model_path) }; + + let settings = settings(params, meta.get_params(), "/tmp/bestdose") + .map_err(|e| format!("Failed to parse settings: {}", e))?; + + let (population_theta, prior_weights) = + parse_prior(&prior_path.to_str().unwrap().to_string(), &settings) + .map_err(|e| format!("Failed to parse prior: {}", e))?; + + let population_weights = prior_weights + .ok_or_else(|| "Prior file must contain a 'prob' column with weights".to_string())?; + + let past_data = if let Some(path) = past_data_path { + let data = data::read_pmetrics(path.to_str().unwrap()) + .map_err(|e| format!("Failed to read past data: {}", e))?; + let subjects = data.subjects(); + if subjects.is_empty() { + return Err("Past data file contains no subjects".to_string()); + } + Some(subjects[0].clone()) + } else { + None + }; + + let target_data = { + let data = data::read_pmetrics(target_data_path.to_str().unwrap()) + .map_err(|e| format!("Failed to read target data: {}", e))?; + let subjects = data.subjects(); + if subjects.is_empty() { + return Err("Target data file contains no subjects".to_string()); + } + subjects[0].clone() + }; + + let target_enum = parse_target_type(target_type)?; + let doserange = DoseRange::new(dose_min, dose_max); + + let problem = BestDoseProblem::new( + &population_theta, + &population_weights, + past_data, + target_data, + time_offset, + eq, + doserange, + bias_weight, + settings, + target_enum, + ) + .map_err(|e| format!("Failed to create BestDose problem: {}", e))?; + + Ok(Self { problem, library }) + } + + pub fn optimize( + &self, + bias_weight: Option, + ) -> std::result::Result { + let configured_problem = match bias_weight { + Some(weight) => self.problem.clone().with_bias_weight(weight), + None => self.problem.clone(), + }; + + configured_problem + .optimize() + .map_err(|e| format!("Optimization failed: {}", e)) + } + + pub fn problem(&self) -> &BestDoseProblem { + &self.problem + } +} -/// Execute bestdose optimization for ODE models pub(crate) fn bestdose_ode( model_path: PathBuf, prior_path: PathBuf, @@ -126,72 +220,20 @@ pub(crate) fn bestdose_ode( target_type: &str, params: List, ) -> std::result::Result { - // 1. Load the model - let (_lib, (eq, meta)) = - unsafe { pmcore::prelude::pharmsol::exa::load::load::(model_path) }; - - // 2. Parse settings from R list - let settings = settings(params, meta.get_params(), "/tmp/bestdose") - .map_err(|e| format!("Failed to parse settings: {}", e))?; - - // 3. Parse the prior (theta + weights) - let (population_theta, prior_weights) = - parse_prior(&prior_path.to_str().unwrap().to_string(), &settings) - .map_err(|e| format!("Failed to parse prior: {}", e))?; - - let population_weights = prior_weights - .ok_or_else(|| "Prior file must contain a 'prob' column with weights".to_string())?; - - // 4. Load past data (if provided) - let past_data = if let Some(path) = past_data_path { - let data = data::read_pmetrics(path.to_str().unwrap()) - .map_err(|e| format!("Failed to read past data: {}", e))?; - let subjects = data.subjects(); - if subjects.is_empty() { - return Err("Past data file contains no subjects".to_string()); - } - Some(subjects[0].clone()) - } else { - None - }; - - // 5. Load target data - let target_data = { - let data = data::read_pmetrics(target_data_path.to_str().unwrap()) - .map_err(|e| format!("Failed to read target data: {}", e))?; - let subjects = data.subjects(); - if subjects.is_empty() { - return Err("Target data file contains no subjects".to_string()); - } - subjects[0].clone() - }; - - // 6. Parse target type - let target_enum = parse_target_type(target_type)?; - - // 7. Create dose range - let doserange = DoseRange::new(dose_min, dose_max); - - // 8. Create and solve bestdose problem - let problem = BestDoseProblem::new( - &population_theta, - &population_weights, - past_data, - target_data, + let handle = BestDoseProblemHandle::new( + model_path, + prior_path, + past_data_path, + target_data_path, time_offset, - eq, - doserange, + dose_min, + dose_max, bias_weight, - settings, - target_enum, - ) - .map_err(|e| format!("Failed to create BestDose problem: {}", e))?; - - let result = problem - .optimize() - .map_err(|e| format!("Optimization failed: {}", e))?; + target_type, + params, + )?; - Ok(result) + handle.optimize(None) } /// Execute bestdose optimization for analytical models (placeholder - not yet supported) @@ -209,3 +251,176 @@ pub(crate) fn bestdose_analytical( ) -> std::result::Result { Err("BestDose for analytical models is not yet supported".to_string()) } + +pub(crate) struct PosteriorSummary { + theta_values: Vec, + theta_dim: (i32, i32), + param_names: Vec, + posterior_weights: Vec, + population_weights: Vec, + bias_weight: f64, + target_type: Target, +} + +fn summarize_problem(problem: &BestDoseProblem) -> PosteriorSummary { + let theta = problem.posterior_theta(); + let matrix = theta.matrix(); + let nrows = matrix.nrows() as i32; + let ncols = matrix.ncols() as i32; + let mut theta_values = vec![0.0; (nrows * ncols) as usize]; + + for col in 0..ncols as usize { + for row in 0..nrows as usize { + theta_values[row + col * nrows as usize] = *matrix.get(row, col); + } + } + + PosteriorSummary { + theta_values, + theta_dim: (nrows, ncols), + param_names: theta.param_names(), + posterior_weights: problem.posterior_weights().to_vec(), + population_weights: problem.population_weights().to_vec(), + bias_weight: problem.bias_weight(), + target_type: problem.target_type(), + } +} + +fn vec_to_doubles(values: Vec, label: &str) -> std::result::Result { + Doubles::try_from(values) + .map_err(|e| format!("Failed to convert {} to doubles: {:?}", label, e)) +} + +fn dims_to_integers(dim: (i32, i32)) -> std::result::Result { + Integers::try_from(vec![dim.0, dim.1]) + .map_err(|e| format!("Failed to convert dims to integers: {:?}", e)) +} + +fn names_to_strings(names: &[String]) -> Strings { + Strings::from_values(names.iter().map(|s| s.as_str())) +} + +pub(crate) fn prepare_bestdose_problem( + model_path: PathBuf, + prior_path: PathBuf, + past_data_path: Option, + target_data_path: PathBuf, + time_offset: Option, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, +) -> std::result::Result<(BestDoseProblemHandle, PosteriorSummary), String> { + let handle = BestDoseProblemHandle::new( + model_path, + prior_path, + past_data_path, + target_data_path, + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params, + )?; + + let summary = summarize_problem(handle.problem()); + Ok((handle, summary)) +} + +pub(crate) fn bestdose_prepare_internal( + model_path: &str, + prior_path: &str, + past_data_path: Nullable, + target_data_path: &str, + time_offset: Nullable, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + kind: &str, +) -> Robj { + RFormatLayer::reset_global_timer(); + let _ = crate::setup_logs(); + + let past_path = past_data_path.into_option().map(PathBuf::from); + let time_offset = time_offset.into_option(); + + let preparation = match kind { + "ode" => prepare_bestdose_problem( + PathBuf::from(model_path), + PathBuf::from(prior_path), + past_path, + PathBuf::from(target_data_path), + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params.clone(), + ), + "analytical" => Err("BestDose for analytical models is not yet supported".to_string()), + other => Err(format!("{} is not a supported model type", other)), + }; + + match preparation { + Ok((handle, summary)) => { + let theta_values = match vec_to_doubles(summary.theta_values, "theta_values") { + Ok(values) => values, + Err(e) => return Robj::from(e), + }; + let theta_dim = match dims_to_integers(summary.theta_dim) { + Ok(dim) => dim, + Err(e) => return Robj::from(e), + }; + let posterior_weights = + match vec_to_doubles(summary.posterior_weights, "posterior_weights") { + Ok(values) => values, + Err(e) => return Robj::from(e), + }; + let population_weights = + match vec_to_doubles(summary.population_weights, "population_weights") { + Ok(values) => values, + Err(e) => return Robj::from(e), + }; + let param_names = names_to_strings(&summary.param_names); + let handle_ptr = ExternalPtr::new(handle); + + let output = list!( + handle = handle_ptr, + theta_values = theta_values, + theta_dim = theta_dim, + param_names = param_names, + posterior_weights = posterior_weights, + population_weights = population_weights, + bias_weight = summary.bias_weight, + target_type = format!("{:?}", summary.target_type), + nspp = summary.theta_dim.0, + n_parameters = summary.theta_dim.1 + ); + + output.into() + } + Err(e) => Robj::from(format!("BestDose prepare failed: {}", e)), + } +} + +pub(crate) fn bestdose_optimize_internal( + handle: ExternalPtr, + bias_weight: Nullable, +) -> Robj { + let weight = bias_weight.into_option(); + + match handle.try_addr() { + Ok(inner) => match inner.optimize(weight) { + Ok(result) => match convert_bestdose_result_to_r(result) { + Ok(robj) => robj, + Err(e) => Robj::from(format!("Failed to convert result: {}", e)), + }, + Err(e) => Robj::from(format!("BestDose optimization failed: {}", e)), + }, + Err(e) => Robj::from(format!("Invalid BestDose handle: {}", e)), + } +} diff --git a/src/rust/src/executor.rs b/src/rust/src/executor.rs index 2be93ee94..8cf5461f0 100755 --- a/src/rust/src/executor.rs +++ b/src/rust/src/executor.rs @@ -42,7 +42,7 @@ pub(crate) fn fit( let data = data::read_pmetrics(data.to_str().unwrap()).expect("Failed to read data"); //dbg!(&data); let mut algorithm = dispatch_algorithm(settings, eq, data)?; - let result = algorithm.fit()?; + let mut result = algorithm.fit()?; result.write_outputs()?; Ok(()) } diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs index 09af6a5b6..ffa43b6fc 100755 --- a/src/rust/src/lib.rs +++ b/src/rust/src/lib.rs @@ -13,6 +13,7 @@ use std::path::PathBuf; use std::process::Command; use tracing_subscriber::layer::SubscriberExt; +use crate::bestdose_executor::BestDoseProblemHandle; use crate::logs::RFormatLayer; fn validate_paths(data_path: &str, model_path: &str) { @@ -316,6 +317,43 @@ fn bestdose( } } +#[extendr] +fn bestdose_prepare( + model_path: &str, + prior_path: &str, + past_data_path: Nullable, + target_data_path: &str, + time_offset: Nullable, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + kind: &str, +) -> Robj { + bestdose_executor::bestdose_prepare_internal( + model_path, + prior_path, + past_data_path, + target_data_path, + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params, + kind, + ) +} + +#[extendr] +fn bestdose_optimize( + handle: ExternalPtr, + bias_weight: Nullable, +) -> Robj { + bestdose_executor::bestdose_optimize_internal(handle, bias_weight) +} + // Macro to generate exports. // This ensures exported functions are registered with R. // See corresponding C code in `entrypoint.c`. @@ -330,6 +368,8 @@ extendr_module! { fn model_parameters; fn setup_logs; fn bestdose; + fn bestdose_prepare; + fn bestdose_optimize; } // To generate the exported function in R, run the following command: