From e932c6a83e400cbefe40a6d474ec71ea3fd08be3 Mon Sep 17 00:00:00 2001 From: Michael Neely Date: Tue, 25 Nov 2025 15:36:49 -0800 Subject: [PATCH 1/3] wip: adding simulation from theta --- R/PM_model.R | 17 +++----- R/PM_sim.R | 119 ++++++++++++++++++++++++++++++++++++++------------- 2 files changed, 95 insertions(+), 41 deletions(-) diff --git a/R/PM_model.R b/R/PM_model.R index 5be502db5..2ed018d99 100644 --- a/R/PM_model.R +++ b/R/PM_model.R @@ -401,7 +401,7 @@ PM_model <- R6::R6Class( ) if (!is.null(x)) { - model_sections <- c("pri", "cov", "sec", "eqn", "lag", "fa", "ini", "out", "err") + #model_sections <- c("pri", "cov", "sec", "eqn", "lag", "fa", "ini", "out", "err") if (is.character(x) && length(x) == 1) { # x is a filename if (!file.exists(x)) { cli::cli_abort(c( @@ -411,11 +411,9 @@ PM_model <- R6::R6Class( } self$arg_list <- private$R6fromFile(x) # read file and populate fields } else if (is.list(x)) { # x is a list in R - purrr::walk(model_sections, \(s) { - if (s %in% names(x)) { - self$arg_list[[s]] <- x[[s]] - } - }) + self$arg_list <- utils::modifyList(self$arg_list, x$arg_list) + self$arg_list$x <- NULL + } else if (inherits(x, "PM_model")) { # x is a PM_model object if (!"arg_list" %in% names(x)) { cli::cli_abort(c( @@ -424,11 +422,7 @@ PM_model <- R6::R6Class( )) } - purrr::walk(model_sections, \(s) { - if (s %in% names(x$arg_list)) { - self$arg_list[[s]] <- x$arg_list[[s]] - } - }) + self$arg_list <- utils::modifyList(self$arg_list, x$arg_list) self$arg_list$x <- NULL } else { cli::cli_abort(c( @@ -740,6 +734,7 @@ PM_model <- R6::R6Class( self$compile() } } else { # default is to compile + browser() self$compile() } }, diff --git a/R/PM_sim.R b/R/PM_sim.R index 68f99ba8b..4c96d3081 100755 --- a/R/PM_sim.R +++ b/R/PM_sim.R @@ -556,6 +556,7 @@ PM_sim <- R6::R6Class( dots <- list(...) old_noise <- c("obsNoise", "obsTimeNoise", "doseNoise", "doseTimeNoise") %in% names(dots) msg <- NULL + simWithTheta <- FALSE # initialize if (any(old_noise)) { cli::cli_abort(c( @@ -665,18 +666,36 @@ PM_sim <- R6::R6Class( # not returning because going on to simulate below - ### This is for loading a saved simulation from file + } else if (is.data.frame (poppar) | tibble::is_tibble (poppar)){ + + theta_ok <- check_valid_theta(poppar) # function at end of file + if (length(theta_ok)>0){ + cli::cli_abort(c("x" = "Your {.arg poppar} object is not valid.", "i" = theta_ok)) + } else { + case <- 8 + simWithTheta <- TRUE + sim <- poppar + class(poppar) <- c(class(poppar), "PM_sim_theta") + final <- list(popPoints = poppar) + } - # CASE 8 - last option, poppar is filename } else { + # CASE 9 - last option, poppar is filename if (file.exists(poppar)) { if (grepl("rds$", poppar, perl = TRUE)) { # poppar is rds filename sim <- readRDS(poppar) - } else { # poppar is a simout.txt name - cli::cli_abort(c( - "x" = "{poppar} is not a compatible file.", - "i" = "Please remake your simulation." - )) + } else { # poppar is another filename + # check if the file has a valid theta format + sim <- readr::read_csv(poppar, show_col_types = FALSE) + theta_ok <- check_valid_theta(poppar) # function at end of file + if (length(theta_ok)>0){ + cli::cli_abort(c("x" = "Your {.arg poppar} object is not valid.", "i" = theta_ok)) + } else { + case <- 9 + simWithTheta <- TRUE + class(sim) <- c(class(sim), "PM_sim_theta") + final <- list(popPoints = sim) + } } # now populate if (inherits(sim, c("PMsim", "PM_sim_data"))) { @@ -685,10 +704,12 @@ PM_sim <- R6::R6Class( private$populate(sim, type = "simlist") } else if (inherits(sim, "PM_sim")) { private$populate(sim, type = "R6sim") + } else if (inherits(sim, "PM_sim_theta")) { + private$populate(sim, type = "theta") } else { cli::cli_abort(c("x" = "{poppar} is not a Pmetrics simulation.")) } - return(self) # end, we have loaded a prior sim + if (!simWithTheta) {return(self)} # end, we have loaded a prior sim } else { cli::cli_abort(c("x" = "{poppar} does not exist in the current working directory.")) } @@ -697,7 +718,7 @@ PM_sim <- R6::R6Class( # If we reach this point, we are creating a new simulation # check model and data - if(case %in% c(2, 3, 7)) { # need model and data if not from PM_result + if(case %in% c(2, 3, 7, 8, 9)) { # need model and data if not from PM_result if (missing(model)) { model <- "model.txt" } # try the default model <- PM_model$new(model) # will abort if can't make PM_model @@ -788,7 +809,7 @@ PM_sim <- R6::R6Class( include = include, exclude = exclude, nsim = nsim, predInt = predInt, covariate = covariate, usePost = usePost, - seed = seed, ode = ode, + seed = seed, ode = ode, simWithTheta = simWithTheta, noise = noise, makecsv = makecsv, outname = outname, clean = clean, quiet = quiet, @@ -898,6 +919,7 @@ PM_sim <- R6::R6Class( include, exclude, nsim, predInt, covariate, usePost, seed, ode, + simWithTheta, noise, makecsv, outname, clean, quiet, nocheck, overwrite, msg) { @@ -906,7 +928,7 @@ PM_sim <- R6::R6Class( ###### POPPAR - npar <- nrow(poppar$popCov) + npar <- ncol(poppar$popPoints) - 1 # number of parameters, minus prob column ###### MODEL @@ -952,7 +974,7 @@ PM_sim <- R6::R6Class( # POSTERIORS ---------------------------------------------- # check if simulating with the posteriors and if so, get all subject IDs - if (usePost) { + if (usePost & !simWithTheta) { if (length(poppar$postMean) == 0) { cli::cli_abort(c( "x" = "Posterior parameters not found.", @@ -982,7 +1004,11 @@ PM_sim <- R6::R6Class( # PARAMETER LIMITS -------------------------------------------------------- - + if (simWithTheta & !all(is.null(limits))){ + cli::cli_warn(c("i" = "Parameter limits are ignored when simulating with fixed theta values.")) + limits <- NULL + } + if (all(is.null(limits))) { # limits are omitted altogether parLimits <- tibble::tibble(par = 1:npar , min = rep(-Inf, npar), max = rep(Inf, npar)) } else if (!any(is.na(limits)) & is.vector(limits)) { # no limit is NA and specified as vector of length 1 or 2 @@ -1014,13 +1040,13 @@ PM_sim <- R6::R6Class( cli::cli_abort(c("x" = "A manual {.arg limits} argument must be a data frame or list with elements {.val par}, {.val min}, and {.val max}.")) } } - + if (!is.null(parLimits) && nrow(parLimits) != npar) { cli::cli_abort(c("x" = "The number of rows in {.arg limits} must match the number of parameters in the model.")) } - - - + + + # COVARIATES ---------------------------------------------------- @@ -1029,6 +1055,11 @@ PM_sim <- R6::R6Class( simWithCov <- FALSE # default is no covariates + if (simWithTheta & !all(is.null(covariate))){ + cli::cli_warn(c("i" = "Simulation with covariates ignored when simulating with fixed theta values.", " " = "Add covariates to theta and to model as primary parameters to simulate.")) + limits <- NULL + } + if(!is.null(covariate)) { simWithCov <- TRUE @@ -1294,6 +1325,11 @@ PM_sim <- R6::R6Class( template <- PM_data$new(template, quiet = TRUE) mod <- PM_model$new(arg_list, quiet = TRUE) + if (simWithTheta & length(postToUse) > 0){ + cli::cli_warn(c("i" = "Simulation with posteriors not possible when simulating with fixed theta values.")) + postToUse <- character(0) + } + if (length(postToUse) > 0) { # simulating from posteriors, each posterior matched to a subject # need to set theta as the posterior mean or median for each subject @@ -1354,18 +1390,22 @@ PM_sim <- R6::R6Class( } else { # postToUse is false - # set theta as nsim rows drawn from prior - thisPrior <- private$getSimPrior( - i = 1, - poppar = poppar, - split = split, - postToUse = NULL, - limits = limits, - seed = seed[1], - nsim = nsim, - toInclude = toInclude, - msg = msg - ) + if(!simWithTheta){ + # set theta as nsim rows drawn from prior + thisPrior <- private$getSimPrior( + i = 1, + poppar = poppar, + split = split, + postToUse = NULL, + limits = limits, + seed = seed[1], + nsim = nsim, + toInclude = toInclude, + msg = msg + ) + } else { + thisPrior <- poppar$popPoints %>% mutate(prob = 1/n()) + } self$data <- private$getSim(thisPrior, template, mod, noise2, msg = msg) } @@ -1589,6 +1629,17 @@ PM_sim <- R6::R6Class( class(simout$data) <- c("PM_sim_data", "list") # ensure class is correct } self$data <- simout$data + } else if (type == "theta") { + theta_list <- list( + obs = NA, + amt = NA, + parValues = simout, + totalSets = NA, + totalMeans = NA, + totalCov = NA + ) + class(theta_list) <- c("PM_sim_data", "list") # ensure class is correct + self$data <- theta_list } return(self) }, # end populate @@ -2479,7 +2530,7 @@ print.summary.PM_sim <- function(x, ...) { } - +########### Internal functions ##################################################### # generate random samples from multivariate, multimodal normal distribution @@ -2573,4 +2624,12 @@ return(list( total_nsim = total_nsim )) +} + +check_valid_theta <- function(df){ + theta_df <- character(0) + if (!"prob" %in% names(df)) {theta_df <- c(theta_df, "Missing {.arg prob} column.")} + if (!all(purrr::map_lgl(df, ~ all(is.numeric(.x))))) {theta_df <- c(theta_df, "All columns must be numeric.")} + if (!all(purrr::map_lgl(df, ~all(!is.na(.x))))) {theta_df <- c("No missing values permitted.")} # check for NAs + return(theta_df) } \ No newline at end of file From 29e1550b632a3564f5adc7e35bca550c8ba64b34 Mon Sep 17 00:00:00 2001 From: Michael Neely Date: Tue, 25 Nov 2025 16:17:50 -0800 Subject: [PATCH 2/3] feat: simulate from theta --- R/PM_model.R | 4 ++-- R/PM_sim.R | 10 +++++----- R/PMoptions.R | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/R/PM_model.R b/R/PM_model.R index 2ed018d99..8b1c445d0 100644 --- a/R/PM_model.R +++ b/R/PM_model.R @@ -411,7 +411,7 @@ PM_model <- R6::R6Class( } self$arg_list <- private$R6fromFile(x) # read file and populate fields } else if (is.list(x)) { # x is a list in R - self$arg_list <- utils::modifyList(self$arg_list, x$arg_list) + self$arg_list <- utils::modifyList(self$arg_list, x) self$arg_list$x <- NULL } else if (inherits(x, "PM_model")) { # x is a PM_model object @@ -734,7 +734,7 @@ PM_model <- R6::R6Class( self$compile() } } else { # default is to compile - browser() + self$compile() } }, diff --git a/R/PM_sim.R b/R/PM_sim.R index 4c96d3081..56afa87e3 100755 --- a/R/PM_sim.R +++ b/R/PM_sim.R @@ -573,7 +573,7 @@ PM_sim <- R6::R6Class( } # CASE 1 - poppar is PM_result - + if (inherits(poppar, "PM_result")) { case <- 1 final <- poppar$final$data # PM_final_data @@ -680,6 +680,7 @@ PM_sim <- R6::R6Class( } } else { + # CASE 9 - last option, poppar is filename if (file.exists(poppar)) { if (grepl("rds$", poppar, perl = TRUE)) { # poppar is rds filename @@ -802,7 +803,7 @@ PM_sim <- R6::R6Class( } # end if !is.null(covariate) # finally, call the simulator, which updates self$data - + private$SIMrun( poppar = final, limits = limits, model = model, data = data, split = split, @@ -1321,7 +1322,6 @@ PM_sim <- R6::R6Class( # CALL SIMULATOR ---------------------------------------------------------------- - template <- PM_data$new(template, quiet = TRUE) mod <- PM_model$new(arg_list, quiet = TRUE) @@ -1404,9 +1404,9 @@ PM_sim <- R6::R6Class( msg = msg ) } else { - thisPrior <- poppar$popPoints %>% mutate(prob = 1/n()) + + thisPrior <- list(thetas = as.data.frame(poppar$popPoints)) } - self$data <- private$getSim(thisPrior, template, mod, noise2, msg = msg) } diff --git a/R/PMoptions.R b/R/PMoptions.R index d9140cb37..5e71043ea 100755 --- a/R/PMoptions.R +++ b/R/PMoptions.R @@ -294,7 +294,7 @@ setPMoptions <- function(launch.app = TRUE) { # Launch the app without trying to launch another browser if(launch.app){ - shiny::runApp(app, launch.browser = FALSE) + shiny::runApp(app, launch.browser = TRUE) } return(invisible(NULL)) From eed84ec9b5c1b33edd8da9f2461965223d7b3398 Mon Sep 17 00:00:00 2001 From: Michael Neely Date: Wed, 26 Nov 2025 11:35:15 -0800 Subject: [PATCH 3/3] feat: simulate using theta; remove template options --- DESCRIPTION | 2 +- R/PM_model.R | 16 +++++++-------- R/PMoptions.R | 40 ++++++++++++++++++------------------- R/extendr-wrappers.R | 8 ++++---- inst/options/PMoptions.json | 3 +-- man/compile_model.Rd | 10 +++++----- src/rust/src/lib.rs | 18 ++++++++--------- 7 files changed, 47 insertions(+), 50 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9c7728bf1..c60c5b796 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,7 +15,7 @@ Authors@R: c( person("David", "Bayard", email = "", role = "ctb"), person("Robert", "Leary", email = "", role = "ctb") ) -Version: 3.0.3 +Version: 3.0.5 URL: https://lapkb.github.io/Pmetrics_rust/ BugReports: https://github.com/LAPKB/Pmetrics_rust/issues SystemRequirements: Cargo (>= 1.82) (Rust's package manager), rustc diff --git a/R/PM_model.R b/R/PM_model.R index 8b1c445d0..80f1a0968 100644 --- a/R/PM_model.R +++ b/R/PM_model.R @@ -1363,21 +1363,20 @@ PM_model <- R6::R6Class( return(invisible(NULL)) } - model_path <- file.path(tempdir(), "model.rs") - private$write_model_to_rust(model_path) + model_rs <- private$create_rs() output_path <- tempfile(pattern = "model_", fileext = ".pmx") cli::cli_inform(c("i" = "Compiling model...")) # path inside Pmetrics package - template_path <- getPMoptions("model_template_path") + template_path <- if (Sys.getenv("env") == "Development") { temporary_path() } else { system.file(package = "Pmetrics")} if (file.access(template_path, 0) == -1 | file.access(template_path, 2) == -1){ cli::cli_abort(c("x" = "Template path {.path {template_path}} does not exist or is not writable.", "i" = "Please set the template path with {.fn setPMoptions} (choose {.emph Compile Options}), to an existing, writable folder." )) } - cat("Using template path:", template_path, "\n") + if (Sys.getenv("env") == "Development") {cat("Using template path:", template_path, "\n")} tryCatch( { - compile_model(model_path, output_path, private$get_primary(), template_path = template_path, kind = tolower(self$model_list$type)) + compile_model(model_rs, template_path, output_path, private$get_primary(), kind = tolower(self$model_list$type)) self$binary_path <- output_path }, error = function(e) { @@ -1387,7 +1386,6 @@ PM_model <- R6::R6Class( } ) - file.remove(model_path) # remove temporary model file return(invisible(self)) }, #' @description @@ -1578,7 +1576,7 @@ PM_model <- R6::R6Class( return(arg_list) }, # end R6fromFile - write_model_to_rust = function(file_path = "main.rs") { + create_rs = function() { # Check if model_list is not NULL if (is.null(self$model_list)) { cli::cli_abort(c("x" = "Model list is empty.", "i" = "Please provide a valid model list.")) @@ -1603,8 +1601,8 @@ PM_model <- R6::R6Class( # Replace placeholders in the base string with actual values from model_list base <- placeholders %>% purrr::reduce(\(x, y) stringr::str_replace(x, stringr::str_c("<", y, ">"), as.character(self$model_list[[y]])), .init = base) - # Write the model to a file - writeLines(base, file_path) + + return(base) }, from_file = function(file_path) { self$model_list <- private$makeR6model(model_filename) diff --git a/R/PMoptions.R b/R/PMoptions.R index 5e71043ea..70a642391 100755 --- a/R/PMoptions.R +++ b/R/PMoptions.R @@ -104,7 +104,7 @@ setPMoptions <- function(launch.app = TRUE) { title = "Pmetrics Options", tags$details( - tags$summary("\u1F4C1 Data File Reading"), + tags$summary("\U0001F4C1 Data File Reading"), selectInput("sep", "Field separator", choices = c(Comma = ",", Semicolon = ";", Tab = "\t"), selected = ","), @@ -115,31 +115,31 @@ setPMoptions <- function(launch.app = TRUE) { ), # Formatting options tags$details( - tags$summary("\u1F4CF Formatting Options"), + tags$summary("\U0001F4CF Formatting Options"), numericInput("digits", "Number of digits to display", value = 3, min = 0, max = 10, step = 1) ), - #C ompile options - tags$details( - tags$summary("\u2699\uFE0F Compile Options"), - markdown("Default Rust model template path is in Pmetrics package installation folder. Change if you have write permission errors."), - tags$div( - style = "display: flex; align-items: flex-start; gap: 8px;", - textAreaInput("model_template_path", NULL, value = system.file(package = "Pmetrics"), autoresize = TRUE), - actionButton("reset_model_template", "Reset to default", class = "btn-secondary") - ), - conditionalPanel( - condition = "input.show == false",selectInput("backend", "Default backend", - choices = c("Rust" = "rust"), - selected = "rust"), - markdown("*Rust is the only backend currently supported by Pmetrics.*") - ) - ), + # #C ompile options + # tags$details( + # tags$summary("\u2699\uFE0F Compile Options"), + # markdown("Default Rust model template path is in Pmetrics package installation folder. Change if you have write permission errors."), + # tags$div( + # style = "display: flex; align-items: flex-start; gap: 8px;", + # textAreaInput("model_template_path", NULL, value = system.file(package = "Pmetrics"), autoresize = TRUE), + # actionButton("reset_model_template", "Reset to default", class = "btn-secondary") + # ), + # conditionalPanel( + # condition = "input.show == false",selectInput("backend", "Default backend", + # choices = c("Rust" = "rust"), + # selected = "rust"), + # markdown("*Rust is the only backend currently supported by Pmetrics.*") + # ) + # ), tags$details( - tags$summary("\u1F4CA Prediction Error Metrics"), + tags$summary("\U0001F4CA Prediction Error Metrics"), br(), checkboxInput("show_metrics", "Display error metrics on obs-pred plots with linear regression", TRUE), selectInput("bias_method", "Bias Method", @@ -173,7 +173,7 @@ setPMoptions <- function(launch.app = TRUE) { ), tags$details( - tags$summary("\u1F4DD Report Generation"), + tags$summary("\U0001F4C8 Report Generation"), selectInput("report_template", "Default report template", choices = c("plotly", "ggplot2"), selected = "plotly") diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R index c7c3b3552..7da4b06c5 100755 --- a/R/extendr-wrappers.R +++ b/R/extendr-wrappers.R @@ -29,14 +29,14 @@ simulate_one <- function(data_path, model_path, spp, kind) .Call(wrap__simulate_ simulate_all <- function(data_path, model_path, theta, kind) .Call(wrap__simulate_all, data_path, model_path, theta, kind) #' Compiles the text representation of a model into a binary file. -#' @param model_path Path to the model file. -#' @param output_path Path to save the compiled model. +#' @param model_rs Text of the rust model. +#' @param template_path Path to the template directory where Rust model is compiled. +#' @param output_path Path to save the compiled model and execute during $fit(). #' @param params List of model parameters. -#' @param template_path Path to the template directory. #' @param kind Kind of model, which can either be "ODE" or "Analytical". #' @return Result of the compilation process. #' @export -compile_model <- function(model_path, output_path, params, template_path, kind) .Call(wrap__compile_model, model_path, output_path, params, template_path, kind) +compile_model <- function(model_rs, template_path, output_path, params, kind) .Call(wrap__compile_model, model_rs, template_path, output_path, params, kind) #' Dummy function to cache compilation artifacts. #' @param template_path Path to the template directory. diff --git a/inst/options/PMoptions.json b/inst/options/PMoptions.json index 9772c9d49..e86a815a2 100755 --- a/inst/options/PMoptions.json +++ b/inst/options/PMoptions.json @@ -7,6 +7,5 @@ "imp_method": "percent_rmbawse", "ic_method": "aic", "report_template": "plotly", - "backend": "rust", - "model_template_path": "" + "backend": "rust" } diff --git a/man/compile_model.Rd b/man/compile_model.Rd index d0727715a..b1ef51a64 100644 --- a/man/compile_model.Rd +++ b/man/compile_model.Rd @@ -4,16 +4,16 @@ \alias{compile_model} \title{Compiles the text representation of a model into a binary file.} \usage{ -compile_model(model_path, output_path, params, template_path, kind) +compile_model(model_rs, template_path, output_path, params, kind) } \arguments{ -\item{model_path}{Path to the model file.} +\item{model_rs}{Text of the rust model.} -\item{output_path}{Path to save the compiled model.} +\item{template_path}{Path to the template directory where Rust model is compiled.} -\item{params}{List of model parameters.} +\item{output_path}{Path to save the compiled model and execute during $fit().} -\item{template_path}{Path to the template directory.} +\item{params}{List of model parameters.} \item{kind}{Kind of model, which can either be "ODE" or "Analytical".} } diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs index d280ceff1..632c61ce9 100755 --- a/src/rust/src/lib.rs +++ b/src/rust/src/lib.rs @@ -6,7 +6,7 @@ mod simulation; use anyhow::Result; use extendr_api::prelude::*; -use pmcore::prelude::{data::read_pmetrics, pharmsol::exa::build, Analytical, ODE}; +use pmcore::prelude::{Analytical, ODE, data::read_pmetrics, pharmsol::{build::temp_path, exa::build}}; use simulation::SimulationRow; use std::process::Command; use tracing_subscriber::layer::SubscriberExt; @@ -179,27 +179,27 @@ fn parse_theta(matrix: RMatrix) -> Vec> { } /// Compiles the text representation of a model into a binary file. -/// @param model_path Path to the model file. -/// @param output_path Path to save the compiled model. +/// @param model_rs Text of the rust model. +/// @param template_path Path to the template directory where Rust model is compiled. +/// @param output_path Path to save the compiled model and execute during $fit(). /// @param params List of model parameters. -/// @param template_path Path to the template directory. /// @param kind Kind of model, which can either be "ODE" or "Analytical". /// @return Result of the compilation process. /// @export #[extendr] fn compile_model( - model_path: &str, + model_rs: &str, + template_path: &str, output_path: &str, params: Strings, - template_path: &str, kind: &str, ) -> Result<()> { let params: Vec = params.iter().map(|x| x.to_string()).collect(); - let model_txt = std::fs::read_to_string(model_path).expect("Failed to read model file"); + let model_rs = model_rs.to_string(); let template_path = std::path::PathBuf::from(template_path); match kind { "ode" => build::compile::( - model_txt, + model_rs, Some(output_path.into()), params.to_vec(), template_path, @@ -208,7 +208,7 @@ fn compile_model( }, )?, "analytical" => build::compile::( - model_txt, + model_rs, Some(output_path.into()), params.to_vec(), template_path,