From c4b21ba15a52966c1246e00614d7e393d2f47022 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Thu, 3 Jul 2025 13:29:49 -0700 Subject: [PATCH 01/20] first working version --- CRAN-SUBMISSION | 3 +++ R/create_event_table.R | 7 +++--- R/sim.R | 6 ++--- inst/cpp/sim.cpp | 57 +++++++++++++++++++++++++++++++++++++++++- 4 files changed, 64 insertions(+), 9 deletions(-) create mode 100644 CRAN-SUBMISSION diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION new file mode 100644 index 00000000..54b4befd --- /dev/null +++ b/CRAN-SUBMISSION @@ -0,0 +1,3 @@ +Version: 1.4.1 +Date: 2025-04-16 16:44:08 UTC +SHA: 9419161ff5c39cefefeda5695f91f8c65bb99ab9 diff --git a/R/create_event_table.R b/R/create_event_table.R index 6567d0ba..4deb4fd5 100644 --- a/R/create_event_table.R +++ b/R/create_event_table.R @@ -152,10 +152,9 @@ create_event_table <- function( design <- dos[order(dos$t, -dos$dose),] if(!is.null(t_obs) && length(t_obs) != 0) { # make sure observation times are in dataset t_obs <- round(t_obs, 6) - t_diff <- setdiff(t_obs, design$t) - if(length(t_diff) > 0) { - design[(length(design[,1])+1) : (length(design[,1])+length(t_diff)),] <- cbind( - t = t_diff, + if(length(t_obs) > 0) { + design[(length(design[,1])+1) : (length(design[,1])+length(t_obs)),] <- cbind( + t = t_obs, dose = 0, type = 0, dum = 0, diff --git a/R/sim.R b/R/sim.R index f7dd0c50..91a87788 100644 --- a/R/sim.R +++ b/R/sim.R @@ -152,9 +152,6 @@ sim <- function (ode = NULL, lagtime <- attr(ode, "lagtime") } } - if(!is.null(lagtime)) { - regimen <- apply_lagtime(regimen, lagtime, parameters, attr(ode, "cmt_mapping")) - } if(!is.null(attr(ode, "dose")$duration_scale)) { regimen <- apply_duration_scale( regimen, @@ -507,11 +504,12 @@ sim <- function (ode = NULL, #################### Main call to ODE solver / analytical eq solver ####################### if(!is.null(ode)) { - tmp <- ode(A_init, design_i, p_i, iov_bins, int_step_size) + tmp <- ode(A_init, design_i, p_i, iov_bins, lagtime, int_step_size) } else { tmp <- analytical_eqn_wrapper(analytical, design_i, p_i) } ##################################################################### + browser() tmp$y <- matrix(unlist(tmp$y), nrow = length(tmp$time), byrow = TRUE) tmp <- as.data.frame(tmp) diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index ce8487a9..1d7b1d21 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -55,6 +55,59 @@ int first_dose_index(std::vector time, std::vector evid) { return(index); } +List apply_lagtime(List design, NumericVector lagtime) { + + List new_design = clone(design); + std::vector times = as >(new_design["t"]); + std::vector evid = as >(new_design["evid"]); + std::vector cmt = as >(new_design["dose_cmt"]); + + // Apply lagtime to dose events (evid == 1) + for(int i = 0; i < times.size(); i++) { + if(evid[i] == 1) { + times[i] += lagtime[cmt[i]-1]; + } + } + + // Create sorted index according to "t" + std::vector indices(times.size()); + std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ... + std::sort(indices.begin(), indices.end(), [×](size_t i1, size_t i2) { + return times[i1] < times[i2]; + }); + + // Reorder all elements in `t` in new_design + std::vector sorted_times(times.size()); + for (size_t i = 0; i < indices.size(); i++) { + sorted_times[i] = times[indices[i]]; + } + new_design["t"] = sorted_times; + + // Sort all other vectors in the design object + for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) { + if (new_design.containsElementNamed(key)) { + SEXP vec = new_design[key]; + if (TYPEOF(vec) == REALSXP) { + std::vector old_vec = as >(vec); + std::vector new_vec(old_vec.size()); + for (size_t i = 0; i < indices.size(); i++) { + new_vec[i] = old_vec[indices[i]]; + } + new_design[key] = new_vec; + } else if (TYPEOF(vec) == INTSXP) { + std::vector old_vec = as >(vec); + std::vector new_vec(old_vec.size()); + for (size_t i = 0; i < indices.size(); i++) { + new_vec[i] = old_vec[indices[i]]; + } + new_design[key] = new_vec; + } + } + } + + return new_design; +} + void set_covariates(int i) { // insert covariates for integration period } @@ -64,7 +117,7 @@ void pk_code (int i, std::vector times, std::vector doses, doubl } // [[Rcpp::export]] -List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_bins, double step_size) { +List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVector iov_bins, NumericVector lagtime, double step_size) { std::vector t; std::vector y; // insert observation variable definition @@ -72,6 +125,8 @@ List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_ std::vector times, doses, dummy, rates; std::vector dose_cmt, dose_type, evid, obs_type, y_type; // insert variable definitions + List design = apply_lagtime(input_design, lagtime); + times = as >(design["t"]); doses = as >(design["dose"]); evid = as >(design["evid"]); From fde2c24462532ab416e7f6c27ad49ab5c37f4936 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 4 Jul 2025 00:35:53 -0700 Subject: [PATCH 02/20] handle lagtime conversion to numeric --- R/sim.R | 1 - inst/cpp/sim.cpp | 31 +++++++++++++++++++++++++++---- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/R/sim.R b/R/sim.R index 91a87788..ec4fae8f 100644 --- a/R/sim.R +++ b/R/sim.R @@ -509,7 +509,6 @@ sim <- function (ode = NULL, tmp <- analytical_eqn_wrapper(analytical, design_i, p_i) } ##################################################################### - browser() tmp$y <- matrix(unlist(tmp$y), nrow = length(tmp$time), byrow = TRUE) tmp <- as.data.frame(tmp) diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index 1d7b1d21..bbf86c1c 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -55,6 +55,27 @@ int first_dose_index(std::vector time, std::vector evid) { return(index); } +NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) { + NumericVector lagtime_numeric; + if (TYPEOF(lagtime) == REALSXP) { + lagtime_numeric = as(lagtime); + } else if (TYPEOF(lagtime) == STRSXP) { + CharacterVector lagtime_char = as(lagtime); + lagtime_numeric = NumericVector(lagtime_char.size()); + for (int i = 0; i < lagtime_char.size(); i++) { + String param_name = lagtime_char[i]; + if (parameters.containsElementNamed(param_name.get_cstring())) { + lagtime_numeric[i] = as(parameters[param_name]); + } else { + lagtime_numeric[i] = 0.0; // default value + } + } + } else { + stop("lagtime must be either numeric or character vector"); + } + return(lagtime_numeric); +} + List apply_lagtime(List design, NumericVector lagtime) { List new_design = clone(design); @@ -117,7 +138,7 @@ void pk_code (int i, std::vector times, std::vector doses, doubl } // [[Rcpp::export]] -List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVector iov_bins, NumericVector lagtime, double step_size) { +List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVector iov_bins, SEXP lagtime, double step_size) { std::vector t; std::vector y; // insert observation variable definition @@ -125,8 +146,11 @@ List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVecto std::vector times, doses, dummy, rates; std::vector dose_cmt, dose_type, evid, obs_type, y_type; // insert variable definitions - List design = apply_lagtime(input_design, lagtime); + // Handle lagtime parameter - can be numeric or character + NumericVector lagtime_numeric = lagtime_to_numeric(lagtime, par); + List design = apply_lagtime(input_design, lagtime_numeric); + times = as >(design["t"]); doses = as >(design["dose"]); evid = as >(design["evid"]); @@ -151,8 +175,7 @@ List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVecto // call ode() once to pre-calculate any initial variables // insert A dAdt state_init set_covariates(0); - t_prv_dose = times[0]; - + // Main call to ODE solver, initialize any variables in ode code ode(A_dum, dAdt_dum, 0); From 672f56175a32e71e731993f3698896c6981b46a4 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 4 Jul 2025 01:02:51 -0700 Subject: [PATCH 03/20] properly convert lagtime + separate functions off --- R/compile_sim_cpp.R | 30 ++++++++++++++- inst/cpp/lagtime.utils.cpp | 75 +++++++++++++++++++++++++++++++++++++ inst/cpp/lagtime.utils.h | 10 +++++ inst/cpp/sim.cpp | 76 +------------------------------------- 4 files changed, 116 insertions(+), 75 deletions(-) create mode 100644 inst/cpp/lagtime.utils.cpp create mode 100644 inst/cpp/lagtime.utils.h diff --git a/R/compile_sim_cpp.R b/R/compile_sim_cpp.R index 417d75c1..6ee5c4c1 100644 --- a/R/compile_sim_cpp.R +++ b/R/compile_sim_cpp.R @@ -270,12 +270,40 @@ compile_sim_cpp <- function( if(length(grep("-w", flg)) == 0) { Sys.setenv("PKG_CXXFLAGS" = paste(flg, "-w")) } + if(compile) { - Rcpp::sourceCpp(code=sim_func, rebuild = TRUE, env = globalenv(), verbose = verbose, showOutput = verbose) + main_cpp <- copy_cpp_files_to_tempfolder(sim_func) + Rcpp::sourceCpp(file = main_cpp, rebuild = TRUE, env = globalenv(), verbose = verbose, showOutput = verbose) Sys.setenv("PKG_CXXFLAGS" = flg) } + return(list( ode = ode_def_cpp, cpp = sim_func )) } + +## Copy cpp and header files to temp folder (first remove existing files) +copy_cpp_files_to_tempfolder <- function(sim_func) { + tmp_folder <- tempdir() + existing_files <- c( + dir(tmp_folder, pattern = "\\.cpp$"), + dir(tmp_folder, pattern = "\\.h$") + ) + for(f in existing_files) { + unlink(file.path(tmp_folder, f)) + } + main_cpp <- file.path(tmp_folder, "sim.cpp") + writeLines(sim_func, main_cpp) + files_to_copy <- setdiff(c( + dir(system.file(package="PKPDsim", "cpp"), pattern = "\\.cpp$"), + dir(system.file(package="PKPDsim", "cpp"), pattern = "\\.h$") + ), "sim.cpp") + for(f in files_to_copy) { + file.copy( + system.file(package="PKPDsim", "cpp", f), + file.path(tmp_folder, f) + ) + } + main_cpp +} \ No newline at end of file diff --git a/inst/cpp/lagtime.utils.cpp b/inst/cpp/lagtime.utils.cpp new file mode 100644 index 00000000..b2d311a4 --- /dev/null +++ b/inst/cpp/lagtime.utils.cpp @@ -0,0 +1,75 @@ +#include "lagtime.utils.h" + +NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) { + NumericVector lagtime_numeric; + if (TYPEOF(lagtime) == REALSXP) { + lagtime_numeric = as(lagtime); + } else if (TYPEOF(lagtime) == STRSXP) { + CharacterVector lagtime_char = as(lagtime); + lagtime_numeric = NumericVector(lagtime_char.size()); + for (int i = 0; i < lagtime_char.size(); i++) { + String param_name = lagtime_char[i]; + if (parameters.containsElementNamed(param_name.get_cstring())) { + lagtime_numeric[i] = as(parameters[param_name]); + } else { + lagtime_numeric[i] = 0.0; // default value + } + } + } else { + stop("lagtime must be either numeric or character vector"); + } + return(lagtime_numeric); +} + +List apply_lagtime(List design, NumericVector lagtime) { + + List new_design = clone(design); + std::vector times = as >(new_design["t"]); + std::vector evid = as >(new_design["evid"]); + std::vector cmt = as >(new_design["dose_cmt"]); + + // Apply lagtime to dose events (evid == 1) + for(int i = 0; i < times.size(); i++) { + if(evid[i] == 1) { + times[i] += lagtime[cmt[i]-1]; + } + } + + // Create sorted index according to "t" + std::vector indices(times.size()); + std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ... + std::sort(indices.begin(), indices.end(), [×](size_t i1, size_t i2) { + return times[i1] < times[i2]; + }); + + // Reorder all elements in `t` in new_design + std::vector sorted_times(times.size()); + for (size_t i = 0; i < indices.size(); i++) { + sorted_times[i] = times[indices[i]]; + } + new_design["t"] = sorted_times; + + // Sort all other vectors in the design object + for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) { + if (new_design.containsElementNamed(key)) { + SEXP vec = new_design[key]; + if (TYPEOF(vec) == REALSXP) { + std::vector old_vec = as >(vec); + std::vector new_vec(old_vec.size()); + for (size_t i = 0; i < indices.size(); i++) { + new_vec[i] = old_vec[indices[i]]; + } + new_design[key] = new_vec; + } else if (TYPEOF(vec) == INTSXP) { + std::vector old_vec = as >(vec); + std::vector new_vec(old_vec.size()); + for (size_t i = 0; i < indices.size(); i++) { + new_vec[i] = old_vec[indices[i]]; + } + new_design[key] = new_vec; + } + } + } + + return new_design; +} \ No newline at end of file diff --git a/inst/cpp/lagtime.utils.h b/inst/cpp/lagtime.utils.h new file mode 100644 index 00000000..93c41ac5 --- /dev/null +++ b/inst/cpp/lagtime.utils.h @@ -0,0 +1,10 @@ +#ifndef LAGTIME_UTILS_H +#define LAGTIME_UTILS_H + +#include +using namespace Rcpp; + +NumericVector lagtime_to_numeric(SEXP lagtime, List parameters); +List apply_lagtime(List design, NumericVector lagtime); + +#endif \ No newline at end of file diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index bbf86c1c..b81e0f1b 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -2,6 +2,8 @@ runge_kutta4 < state_type > stepr; // runge_kutta_cash_karp54 < state_type > stepr; +#include "lagtime.utils.h" + struct push_back_solution { std::vector< state_type >& m_states; @@ -55,80 +57,6 @@ int first_dose_index(std::vector time, std::vector evid) { return(index); } -NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) { - NumericVector lagtime_numeric; - if (TYPEOF(lagtime) == REALSXP) { - lagtime_numeric = as(lagtime); - } else if (TYPEOF(lagtime) == STRSXP) { - CharacterVector lagtime_char = as(lagtime); - lagtime_numeric = NumericVector(lagtime_char.size()); - for (int i = 0; i < lagtime_char.size(); i++) { - String param_name = lagtime_char[i]; - if (parameters.containsElementNamed(param_name.get_cstring())) { - lagtime_numeric[i] = as(parameters[param_name]); - } else { - lagtime_numeric[i] = 0.0; // default value - } - } - } else { - stop("lagtime must be either numeric or character vector"); - } - return(lagtime_numeric); -} - -List apply_lagtime(List design, NumericVector lagtime) { - - List new_design = clone(design); - std::vector times = as >(new_design["t"]); - std::vector evid = as >(new_design["evid"]); - std::vector cmt = as >(new_design["dose_cmt"]); - - // Apply lagtime to dose events (evid == 1) - for(int i = 0; i < times.size(); i++) { - if(evid[i] == 1) { - times[i] += lagtime[cmt[i]-1]; - } - } - - // Create sorted index according to "t" - std::vector indices(times.size()); - std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ... - std::sort(indices.begin(), indices.end(), [×](size_t i1, size_t i2) { - return times[i1] < times[i2]; - }); - - // Reorder all elements in `t` in new_design - std::vector sorted_times(times.size()); - for (size_t i = 0; i < indices.size(); i++) { - sorted_times[i] = times[indices[i]]; - } - new_design["t"] = sorted_times; - - // Sort all other vectors in the design object - for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) { - if (new_design.containsElementNamed(key)) { - SEXP vec = new_design[key]; - if (TYPEOF(vec) == REALSXP) { - std::vector old_vec = as >(vec); - std::vector new_vec(old_vec.size()); - for (size_t i = 0; i < indices.size(); i++) { - new_vec[i] = old_vec[indices[i]]; - } - new_design[key] = new_vec; - } else if (TYPEOF(vec) == INTSXP) { - std::vector old_vec = as >(vec); - std::vector new_vec(old_vec.size()); - for (size_t i = 0; i < indices.size(); i++) { - new_vec[i] = old_vec[indices[i]]; - } - new_design[key] = new_vec; - } - } - } - - return new_design; -} - void set_covariates(int i) { // insert covariates for integration period } From 128be21c291576b0b81b3b90da15a45e75a3c090 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 4 Jul 2025 02:20:13 -0700 Subject: [PATCH 04/20] misc fix --- R/compile_sim_cpp.R | 25 ++++++++++++++++++++----- R/sim.R | 2 +- R/sim_core.R | 1 + 3 files changed, 22 insertions(+), 6 deletions(-) diff --git a/R/compile_sim_cpp.R b/R/compile_sim_cpp.R index 6ee5c4c1..13b551eb 100644 --- a/R/compile_sim_cpp.R +++ b/R/compile_sim_cpp.R @@ -272,11 +272,19 @@ compile_sim_cpp <- function( } if(compile) { - main_cpp <- copy_cpp_files_to_tempfolder(sim_func) - Rcpp::sourceCpp(file = main_cpp, rebuild = TRUE, env = globalenv(), verbose = verbose, showOutput = verbose) + res <- copy_cpp_files_to_tempfolder(sim_func) + withr::with_dir(res$folder, { + Rcpp::sourceCpp( + file = res$main_cpp, + rebuild = TRUE, + env = globalenv(), + verbose = verbose, + showOutput = verbose + ) + }) Sys.setenv("PKG_CXXFLAGS" = flg) } - + return(list( ode = ode_def_cpp, cpp = sim_func @@ -288,12 +296,15 @@ copy_cpp_files_to_tempfolder <- function(sim_func) { tmp_folder <- tempdir() existing_files <- c( dir(tmp_folder, pattern = "\\.cpp$"), - dir(tmp_folder, pattern = "\\.h$") + dir(tmp_folder, pattern = "\\.h$"), + dir(tmp_folder, pattern = "\\.o$") ) for(f in existing_files) { unlink(file.path(tmp_folder, f)) } main_cpp <- file.path(tmp_folder, "sim.cpp") + if(file.exists(main_cpp)) + unlink(main_cpp) writeLines(sim_func, main_cpp) files_to_copy <- setdiff(c( dir(system.file(package="PKPDsim", "cpp"), pattern = "\\.cpp$"), @@ -305,5 +316,9 @@ copy_cpp_files_to_tempfolder <- function(sim_func) { file.path(tmp_folder, f) ) } - main_cpp + list( + main_cpp = main_cpp, + folder = tmp_folder + ) + } \ No newline at end of file diff --git a/R/sim.R b/R/sim.R index ec4fae8f..c627c1ff 100644 --- a/R/sim.R +++ b/R/sim.R @@ -93,7 +93,7 @@ sim <- function (ode = NULL, n_ind = 1, event_table = NULL, regimen = NULL, - lagtime = NULL, + lagtime = c(0), covariates = NULL, covariates_table = NULL, covariates_implementation = list(), diff --git a/R/sim_core.R b/R/sim_core.R index 44fe7bd8..e8c124b1 100644 --- a/R/sim_core.R +++ b/R/sim_core.R @@ -16,6 +16,7 @@ sim_core <- function( design = sim_object$design, par = sim_object$p, iov_bins = sim_object$iov_bins, + lagtime = c(0), step_size = sim_object$int_step_size) out <- data.frame(t = tmp$time, y = tmp$obs, obs_type = tmp$obs_type) if(duplicate_t_obs) { From d919a6b4d7863f800fd9e6973a13d4e29dce9f62 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 4 Jul 2025 02:29:11 -0700 Subject: [PATCH 05/20] rollback split off functions --- R/compile_sim_cpp.R | 17 +++++----- inst/cpp/sim.cpp | 76 +++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 81 insertions(+), 12 deletions(-) diff --git a/R/compile_sim_cpp.R b/R/compile_sim_cpp.R index 13b551eb..fb623108 100644 --- a/R/compile_sim_cpp.R +++ b/R/compile_sim_cpp.R @@ -272,16 +272,13 @@ compile_sim_cpp <- function( } if(compile) { - res <- copy_cpp_files_to_tempfolder(sim_func) - withr::with_dir(res$folder, { - Rcpp::sourceCpp( - file = res$main_cpp, - rebuild = TRUE, - env = globalenv(), - verbose = verbose, - showOutput = verbose - ) - }) + Rcpp::sourceCpp( + code = sim_func, + rebuild = TRUE, + env = globalenv(), + verbose = verbose, + showOutput = verbose + ) Sys.setenv("PKG_CXXFLAGS" = flg) } diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index b81e0f1b..60635cb4 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -2,8 +2,6 @@ runge_kutta4 < state_type > stepr; // runge_kutta_cash_karp54 < state_type > stepr; -#include "lagtime.utils.h" - struct push_back_solution { std::vector< state_type >& m_states; @@ -65,6 +63,80 @@ void pk_code (int i, std::vector times, std::vector doses, doubl // insert custom pk event code } +NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) { + NumericVector lagtime_numeric; + if (TYPEOF(lagtime) == REALSXP) { + lagtime_numeric = as(lagtime); + } else if (TYPEOF(lagtime) == STRSXP) { + CharacterVector lagtime_char = as(lagtime); + lagtime_numeric = NumericVector(lagtime_char.size()); + for (int i = 0; i < lagtime_char.size(); i++) { + String param_name = lagtime_char[i]; + if (parameters.containsElementNamed(param_name.get_cstring())) { + lagtime_numeric[i] = as(parameters[param_name]); + } else { + lagtime_numeric[i] = 0.0; // default value + } + } + } else { + stop("lagtime must be either numeric or character vector"); + } + return(lagtime_numeric); +} + +List apply_lagtime(List design, NumericVector lagtime) { + + List new_design = clone(design); + std::vector times = as >(new_design["t"]); + std::vector evid = as >(new_design["evid"]); + std::vector cmt = as >(new_design["dose_cmt"]); + + // Apply lagtime to dose events (evid == 1) + for(int i = 0; i < times.size(); i++) { + if(evid[i] == 1) { + times[i] += lagtime[cmt[i]-1]; + } + } + + // Create sorted index according to "t" + std::vector indices(times.size()); + std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ... + std::sort(indices.begin(), indices.end(), [×](size_t i1, size_t i2) { + return times[i1] < times[i2]; + }); + + // Reorder all elements in `t` in new_design + std::vector sorted_times(times.size()); + for (size_t i = 0; i < indices.size(); i++) { + sorted_times[i] = times[indices[i]]; + } + new_design["t"] = sorted_times; + + // Sort all other vectors in the design object + for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) { + if (new_design.containsElementNamed(key)) { + SEXP vec = new_design[key]; + if (TYPEOF(vec) == REALSXP) { + std::vector old_vec = as >(vec); + std::vector new_vec(old_vec.size()); + for (size_t i = 0; i < indices.size(); i++) { + new_vec[i] = old_vec[indices[i]]; + } + new_design[key] = new_vec; + } else if (TYPEOF(vec) == INTSXP) { + std::vector old_vec = as >(vec); + std::vector new_vec(old_vec.size()); + for (size_t i = 0; i < indices.size(); i++) { + new_vec[i] = old_vec[indices[i]]; + } + new_design[key] = new_vec; + } + } + } + + return new_design; +} + // [[Rcpp::export]] List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVector iov_bins, SEXP lagtime, double step_size) { std::vector t; From 218784dde2f0b84edf4e226f03e4960cce9c8b61 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 4 Jul 2025 03:41:08 -0700 Subject: [PATCH 06/20] misc fixes --- R/sim.R | 3 +++ inst/cpp/sim.cpp | 18 ++++++++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/R/sim.R b/R/sim.R index c627c1ff..ae4b59fb 100644 --- a/R/sim.R +++ b/R/sim.R @@ -231,6 +231,9 @@ sim <- function (ode = NULL, } else { size <- attr(analytical, "size") } + if(is.null(lagtime)) { + lagtime <- rep(0, size) # needs to have at least 1 zero value, cannot be NULL when passed to cpp func + } if(is.null(ode) && is.null(analytical)) { stop("Please specify at least the required arguments 'ode' or 'analytical' for simulations.") } diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index 60635cb4..7b8001f5 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -84,13 +84,27 @@ NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) { return(lagtime_numeric); } -List apply_lagtime(List design, NumericVector lagtime) { +List apply_lagtime(List design, NumericVector tlag, int n_comp) { List new_design = clone(design); std::vector times = as >(new_design["t"]); std::vector evid = as >(new_design["evid"]); std::vector cmt = as >(new_design["dose_cmt"]); + NumericVector lagtime; + if(tlag.size() < n_comp) { // fill in with zeroes, if needed + int current_size = tlag.size(); + lagtime = NumericVector(n_comp); + for(int i = 0; i < current_size; i++) { + lagtime[i] = tlag[i]; + } + for(int i = current_size; i < n_comp; i++) { + lagtime[i] = 0.0; + } + } else { + lagtime = tlag; + } + // Apply lagtime to dose events (evid == 1) for(int i = 0; i < times.size(); i++) { if(evid[i] == 1) { @@ -149,7 +163,7 @@ List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVecto // Handle lagtime parameter - can be numeric or character NumericVector lagtime_numeric = lagtime_to_numeric(lagtime, par); - List design = apply_lagtime(input_design, lagtime_numeric); + List design = apply_lagtime(input_design, lagtime_numeric, n_comp); times = as >(design["t"]); doses = as >(design["dose"]); From 09ec2b66aaf545a8c7bf173b5d49611b9c80a5b4 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 02:30:34 -0700 Subject: [PATCH 07/20] filter out dose event times afterwards --- R/sim.R | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/R/sim.R b/R/sim.R index ae4b59fb..e8de7a4d 100644 --- a/R/sim.R +++ b/R/sim.R @@ -602,6 +602,12 @@ sim <- function (ode = NULL, all_names <- unique(c(par_names, cov_names, var_names)) all_names <- intersect(all_names, names(comb)) # only cols that appear in data + ## remove dose-times from regimen, leave ones that were also in t_obs + # comb |> dplyr::filter(t >= 24 & t <= 27 & comp == 2) + dose_times <- design_i[design_i$evid > 0, ]$t + # comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$comp, comb$t, sep="_"))),] + comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, sep="_"))),] + if(!extra_t_obs) { ## include the observations at which a bolus dose is added into the output object too comb <- comb[!duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, sep="_")),] From 7c8c95b8b9b12b3992b9d0646173f46bd93a0e74 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 03:59:29 -0700 Subject: [PATCH 08/20] clean up + additional fixes --- R/apply_lagtime.R | 40 --------------- R/sim.R | 13 +++-- inst/cpp/sim.cpp | 2 +- tests/testthat/test_apply_lagtime.R | 60 ---------------------- tests/testthat/test_is_positive_definite.R | 4 +- tests/testthat/test_new_ode_model.R | 5 +- 6 files changed, 15 insertions(+), 109 deletions(-) delete mode 100644 R/apply_lagtime.R delete mode 100644 tests/testthat/test_apply_lagtime.R diff --git a/R/apply_lagtime.R b/R/apply_lagtime.R deleted file mode 100644 index 85718ccf..00000000 --- a/R/apply_lagtime.R +++ /dev/null @@ -1,40 +0,0 @@ -#' Apply lagtime to a regimen -#' -#' @param regimen PKPDsim regimen -#' @param lagtime lagtime object, either single value / parameter name or vector of values/parameter names for all compartments. -#' @param parameters parameter list, required if parameters are specified. -#' @param cmt_mapping map of administration types to compartments, e.g. `list("oral" = 1, "infusion" = 2, "bolus" = 2)`. -#' -#' @export -#' @return Original regimen with lagtime added to dose times -apply_lagtime <- function( - regimen, - lagtime, - parameters, - cmt_mapping = NULL -) { - if(class(lagtime) %in% c("numeric", "integer")) { - if(length(lagtime) == 1) { - regimen$dose_times <- regimen$dose_times + lagtime - } else { - regimen$dose_times <- regimen$dose_times + lagtime[regimen$cmt] - } - } - if(class(lagtime) %in% c("character")) { - if(length(lagtime) == 1) { - regimen$dose_times <- regimen$dose_times + parameters[[lagtime]] - } else { - if(is.null(regimen$cmt)) { - if(!is.null(cmt_mapping)) { - regimen$cmt <- as.numeric(cmt_mapping[regimen$type]) - } else { - regimen$cmt <- rep(1, length(regimen$dose_times)) - } - } - par_tmp <- parameters - par_tmp[["0"]] <- 0 - regimen$dose_times <- regimen$dose_times + as.numeric(unlist(par_tmp[lagtime[regimen$cmt]])) - } - } - return(regimen) -} diff --git a/R/sim.R b/R/sim.R index e8de7a4d..94b2478d 100644 --- a/R/sim.R +++ b/R/sim.R @@ -93,7 +93,7 @@ sim <- function (ode = NULL, n_ind = 1, event_table = NULL, regimen = NULL, - lagtime = c(0), + lagtime = NULL, covariates = NULL, covariates_table = NULL, covariates_implementation = list(), @@ -152,6 +152,9 @@ sim <- function (ode = NULL, lagtime <- attr(ode, "lagtime") } } + if(is.null(lagtime)) { + lagtime <- c(0) + } if(!is.null(attr(ode, "dose")$duration_scale)) { regimen <- apply_duration_scale( regimen, @@ -602,10 +605,10 @@ sim <- function (ode = NULL, all_names <- unique(c(par_names, cov_names, var_names)) all_names <- intersect(all_names, names(comb)) # only cols that appear in data - ## remove dose-times from regimen, leave ones that were also in t_obs - # comb |> dplyr::filter(t >= 24 & t <= 27 & comp == 2) - dose_times <- design_i[design_i$evid > 0, ]$t - # comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$comp, comb$t, sep="_"))),] + ## remove dose-times from regimen + ## but leave ones that were also in t_obs + ## also leave extra obs when bolus, because this may be needed (controllec below using `extra_t_obs`) + dose_times <- design_i[design_i$evid > 0 & design_i$type != 0, ]$t comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, sep="_"))),] if(!extra_t_obs) { diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index 7b8001f5..2c9604e8 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -105,7 +105,7 @@ List apply_lagtime(List design, NumericVector tlag, int n_comp) { lagtime = tlag; } - // Apply lagtime to dose events (evid == 1) + // Apply lagtime to dose events (evid == 1) for(int i = 0; i < times.size(); i++) { if(evid[i] == 1) { times[i] += lagtime[cmt[i]-1]; diff --git a/tests/testthat/test_apply_lagtime.R b/tests/testthat/test_apply_lagtime.R deleted file mode 100644 index 5ff1c599..00000000 --- a/tests/testthat/test_apply_lagtime.R +++ /dev/null @@ -1,60 +0,0 @@ -reg1 <- new_regimen( - amt = 1000, - n = 12, - interval = 12, - type = "oral" -) - -test_that("lagtime applied to regimens using parameters", { - lag1 <- apply_lagtime( - regimen = reg1, - lagtime = "TLAG", - parameters = list(CL = 5, V = 50, TLAG = .5) - ) - expect_true("regimen" %in% class(lag1)) - expect_equal(lag1$dose_times, reg1$dose_times + 0.5) -}) - -test_that("lagtime applied to regimens using lagtime arg", { - lag2 <- apply_lagtime(regimen = reg1, lagtime = .75) - expect_true("regimen" %in% class(lag2)) - expect_equal(lag2$dose_times, reg1$dose_times + 0.75) -}) - -test_that("lagtime applied when multiple compartments and no compartment specified in regimen", { - lag3 <- apply_lagtime( - regimen = reg1, - lagtime = c("TLAG", 0, 0), - parameters = list(CL = 5, V = 50, TLAG = .5) - ) - expect_true("regimen" %in% class(lag3)) - expect_equal(lag3$dose_times, reg1$dose_times + 0.5) -}) - -test_that("lagtime applied when multiple compartments and no compartment specified in regimen but cmt_mapping available (oral)", { - lag4a <- apply_lagtime( - regimen = reg1, - lagtime = c("TLAG", 0, 0), - parameters = list(CL = 5, V = 50, TLAG = .5), - cmt_mapping = list(oral = 1, infusion = 2, bolus = 2) - ) - expect_true("regimen" %in% class(lag4a)) - expect_equal(lag4a$dose_times, reg1$dose_times + 0.5) -}) - -test_that("lagtime applied when multiple compartments and no compartment specified in regimen but cmt_mapping available (infusion): dose times should not be altered", { - reg2 <- new_regimen( - amt = 1000, - n = 12, - interval = 12, - type = "infusion" - ) - lag4b <- apply_lagtime( - regimen = reg2, - lagtime = c("TLAG", 0, 0), - parameters = list(CL = 5, V = 50, TLAG = .5), - cmt_mapping = list(oral = 1, infusion = 2, bolus = 2) - ) - expect_true("regimen" %in% class(lag4b)) - expect_equal(lag4b$dose_times, reg2$dose_times) -}) diff --git a/tests/testthat/test_is_positive_definite.R b/tests/testthat/test_is_positive_definite.R index 4ccc82be..8dd334cd 100644 --- a/tests/testthat/test_is_positive_definite.R +++ b/tests/testthat/test_is_positive_definite.R @@ -105,7 +105,9 @@ test_that("Matrix where numerical solver issues result in very tiny complex comp ), nrow = 16 ) - is_x_positive_def <- is_positive_definite(x) + suppressWarnings( + is_x_positive_def <- is_positive_definite(x) + ) expect_true(is_x_positive_def) }) diff --git a/tests/testthat/test_new_ode_model.R b/tests/testthat/test_new_ode_model.R index 5b5141a6..b51b3ec4 100644 --- a/tests/testthat/test_new_ode_model.R +++ b/tests/testthat/test_new_ode_model.R @@ -11,7 +11,7 @@ test_that("Double absorption models work appropriately", { bioav = c(0.3, 0.6, 0.9) ), lagtime = c(0, 1.5), - cpp = F + cpp = T ) reg <- new_regimen(amt = 100, n = 12, interval = 12, type = "oral") res <- sim( @@ -25,7 +25,8 @@ test_that("Double absorption models work appropriately", { ), only_obs = F, output_include = list(variables=TRUE), - t_obs = seq(0, 24, .1) + t_obs = seq(0, 24, .1), + return_design = F ) expect_equal(attr(mod, "dose")$duplicate, 2) expect_equal(round(res[res$comp == 1 & res$t == 0.1, ]$y, 2), 27.15) # abs comp 1 From c420315930efd63360d8bbb7c2236e972411a765 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 05:22:53 -0700 Subject: [PATCH 09/20] fixes for all current tests --- NAMESPACE | 1 - R/sim.R | 4 +- R/sim_core.R | 17 ++++--- inst/cpp/sim.cpp | 20 ++++----- man/apply_lagtime.Rd | 23 ---------- tests/testthat/setup.R | 10 +++++ tests/testthat/test_sim.R | 75 +++++++++++++++---------------- tests/testthat/test_sim_lagtime.R | 28 +++++++----- 8 files changed, 85 insertions(+), 93 deletions(-) delete mode 100644 man/apply_lagtime.Rd diff --git a/NAMESPACE b/NAMESPACE index 145be003..92cbb2de 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,7 +15,6 @@ export(advan_create_data) export(advan_parse_output) export(advan_process_infusion_doses) export(apply_duration_scale) -export(apply_lagtime) export(available_default_literature_models) export(calc_auc_analytic) export(calc_ss_analytic) diff --git a/R/sim.R b/R/sim.R index 94b2478d..f49c8b2e 100644 --- a/R/sim.R +++ b/R/sim.R @@ -608,8 +608,8 @@ sim <- function (ode = NULL, ## remove dose-times from regimen ## but leave ones that were also in t_obs ## also leave extra obs when bolus, because this may be needed (controllec below using `extra_t_obs`) - dose_times <- design_i[design_i$evid > 0 & design_i$type != 0, ]$t - comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, sep="_"))),] + dose_times <- design_i[design_i$evid > 0, ]$t + comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, comb$y, sep="_"))),] if(!extra_t_obs) { ## include the observations at which a bolus dose is added into the output object too diff --git a/R/sim_core.R b/R/sim_core.R index e8c124b1..23c0df5e 100644 --- a/R/sim_core.R +++ b/R/sim_core.R @@ -11,13 +11,16 @@ sim_core <- function( sim_object = NULL, ode, duplicate_t_obs = FALSE, - t_init = 0) { - tmp <- ode(A = sim_object$A_init, - design = sim_object$design, - par = sim_object$p, - iov_bins = sim_object$iov_bins, - lagtime = c(0), - step_size = sim_object$int_step_size) + t_init = 0 +) { + tmp <- ode( + A = sim_object$A_init, + design = sim_object$design, + par = sim_object$p, + iov_bins = sim_object$iov_bins, + lagtime = c(0), + step_size = sim_object$int_step_size + ) out <- data.frame(t = tmp$time, y = tmp$obs, obs_type = tmp$obs_type) if(duplicate_t_obs) { # use match to ensure that duplicates in t_obs is possible diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index 2c9604e8..3adb0659 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -152,7 +152,7 @@ List apply_lagtime(List design, NumericVector tlag, int n_comp) { } // [[Rcpp::export]] -List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVector iov_bins, SEXP lagtime, double step_size) { +List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_bins, SEXP lagtime, double step_size) { std::vector t; std::vector y; // insert observation variable definition @@ -163,16 +163,16 @@ List sim_wrapper_cpp (NumericVector A, List input_design, List par, NumericVecto // Handle lagtime parameter - can be numeric or character NumericVector lagtime_numeric = lagtime_to_numeric(lagtime, par); - List design = apply_lagtime(input_design, lagtime_numeric, n_comp); + List events = apply_lagtime(design, lagtime_numeric, n_comp); - times = as >(design["t"]); - doses = as >(design["dose"]); - evid = as >(design["evid"]); - dummy = as >(design["dum"]); - rates = as >(design["rate"]); - dose_cmt = as >(design["dose_cmt"]); - dose_type = as >(design["type"]); - obs_type = as >(design["obs_type"]); + times = as >(events["t"]); + doses = as >(events["dose"]); + evid = as >(events["evid"]); + dummy = as >(events["dum"]); + rates = as >(events["rate"]); + dose_cmt = as >(events["dose_cmt"]); + dose_type = as >(events["type"]); + obs_type = as >(events["obs_type"]); int len = times.size(); int start; memset(rate, 0, sizeof(rate)); diff --git a/man/apply_lagtime.Rd b/man/apply_lagtime.Rd deleted file mode 100644 index 3718dc5c..00000000 --- a/man/apply_lagtime.Rd +++ /dev/null @@ -1,23 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/apply_lagtime.R -\name{apply_lagtime} -\alias{apply_lagtime} -\title{Apply lagtime to a regimen} -\usage{ -apply_lagtime(regimen, lagtime, parameters, cmt_mapping = NULL) -} -\arguments{ -\item{regimen}{PKPDsim regimen} - -\item{lagtime}{lagtime object, either single value / parameter name or vector of values/parameter names for all compartments.} - -\item{parameters}{parameter list, required if parameters are specified.} - -\item{cmt_mapping}{map of administration types to compartments, e.g. \code{list("oral" = 1, "infusion" = 2, "bolus" = 2)}.} -} -\value{ -Original regimen with lagtime added to dose times -} -\description{ -Apply lagtime to a regimen -} diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 9dc6bb2c..d4e2f51c 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -1,6 +1,16 @@ mod_1cmt_iv <- new_ode_model("pk_1cmt_iv") mod_2cmt_iv <- new_ode_model("pk_2cmt_iv") mod_1cmt_oral <- new_ode_model("pk_1cmt_oral") +mod_1cmt_oral_lagtime <- new_ode_model( + code = " + dAdt[0] = -KA * A[0] + dAdt[1] = +KA * A[0] -(CL/V) * A[1] + ", + lagtime = c("TLAG", 0), + obs = list(cmt = 2, scale = "V"), + dose = list(cmt = 1, bioav = 1), + parameters = pars +) oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence factor code = " if(t<168.0) { diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R index 0e3892b7..1be57a6d 100644 --- a/tests/testthat/test_sim.R +++ b/tests/testthat/test_sim.R @@ -20,22 +20,25 @@ test_that("return_event_table=TRUE returns an appropriate event table", { ) expect_equal( evtab1, - data.frame( - t = c(0, 2, 12, 14, 24, 26, 48, 48), - dose = c(100, 0, 100, 0, 100, 0, 0, 0), - type = c(1, 1, 1, 1, 1, 1, 0, 0), - dum = c(0, 1, 0, 1, 0, 1, 0, 0), - dose_cmt = c(1, 1, 1, 1, 1, 1, 0, 0), - t_inf = c(2, 0, 2, 0, 2, 0, 0, 0), - evid = c(1, 2, 1, 2, 1, 2, 0, 0), - bioav = c(1, 0, 1, 0, 1, 0, 0, 0), - rate = c(50,-50, 50,-50, 50,-50, 0, 0), - obs_type = c(0, 1, 0, 0, 0, 0, 1, 1) - ) - ) + structure(list( + t = c(0, 2, 2, 12, 14, 24, 26, 48, 48), + dose = c(100, 0, 0, 100, 0, 100, 0, 0, 0), + type = c(1, 0, 1, 1, 1, 1, 1, 0, 0), + dum = c(0, 0, 1, 0, 1, 0, 1, 0, 0), + dose_cmt = c(1, 0, 1, 1, 1, 1, 1, 0, 0), + t_inf = c(2, 0, 0, 2, 0, 2, 0, 0, 0), + evid = c(1, 0, 2, 1, 2, 1, 2, 0, 0), + bioav = c(1, 0, 0, 1, 0, 1, 0, 0, 0), + rate = c(50, 0, -50, 50, -50, 50, -50, 0, 0), + obs_type = c(0, 1, 1, 0, 0, 0, 0, 1, 1) + ), + row.names = c(1L, 3L, 2L, 4L, 5L, 6L, 7L, 8L, 9L), + class = "data.frame" + )) }) test_that("return_event_table=TRUE returns an appropriate event table with covariate", { + mod_1cmt_iv <- new_ode_model("pk_1cmt_iv") covs <- list(CRCL = new_covariate(value = c(70, 80), t = c(0, 24)), WT = new_covariate(70)) evtab2 <- sim_ode( mod_1cmt_iv, @@ -47,30 +50,27 @@ test_that("return_event_table=TRUE returns an appropriate event table with covar ) expect_equal( evtab2, - structure( - list( - t = c(0, 0, 2, 12, 14, 24, 24, 26, 48), - dose = c(0, 100, 0, 100, 0, 0, 100, 0, 0), - type = c(0, 1, 1, 1, 1, 0, 1, 1, 0), - dum = c(0, 0, 1, 0, 1, 0, 0, 1, 0), - dose_cmt = c(0, 1, 1, 1, 1, 0, 1, 1, 0), - t_inf = c(0, 2, 0, 2, 0, 0, 2, 0, 0), - evid = c(2, 1, 2, 1, 2, 2, 1, 2, 0), - bioav = c(1, 1, 0, 1, 0, 1, 1, 0, 0), - rate = c(0, 50,-50, 50,-50, 0, 50,-50, 0), - cov_CRCL = c(70, 70, 70, 70, 70, 80, 80, 80, 80), - cov_t_CRCL = c(0, 0, 0, 0, 0, 24, 24, 24, 24), - gradients_CRCL = c(0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0, 0, 0, 0), - cov_WT = c(70, 70, 70, 70, 70, 70, 70, 70, 70), - cov_t_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0), - gradients_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0), - obs_type = c(0, 0, 1, 0, 0, 0, 0, 0, 1) - ), - row.names = c(2L, 1L, 3L, 4L, 5L, - 6L, 7L, 8L, 9L), - class = "data.frame" - ) - ) + structure(list( + t = c(0, 0, 2, 2, 12, 14, 24, 24, 26, 48), + dose = c(0, 100, 0, 0, 100, 0, 0, 100, 0, 0), + type = c(0, 1, 0, 1, 1, 1, 0, 1, 1, 0), + dum = c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0), + dose_cmt = c(0, 1, 0, 1, 1, 1, 0, 1, 1, 0), + t_inf = c(0, 2, 0, 0, 2, 0, 0, 2, 0, 0), + evid = c(2, 1, 0, 2, 1, 2, 2, 1, 2, 0), + bioav = c(1, 1, 0, 0, 1, 0, 1, 1, 0, 0), + rate = c(0, 50, 0, -50, 50, -50, 0, 50, -50, 0), + cov_CRCL = c(70, 70, 70, 70, 70, 70, 80, 80, 80, 80), + cov_t_CRCL = c(0, 0, 0, 0, 0, 0, 24, 24, 24, 24), + gradients_CRCL = c(0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0.416666666666667, 0, 0, 0, 0), + cov_WT = c(70, 70, 70, 70, 70, 70, 70, 70, 70, 70), + cov_t_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), + gradients_WT = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), + obs_type = c(0, 0, 1, 1, 0, 0, 0, 0, 0, 1) + ), + row.names = c(2L, 1L, 4L, 3L, 5L, 6L, 8L, 7L, 9L, 10L), + class = "data.frame" + )) }) test_that("sim works properly for a model where bioavailability is dependent on dose", { @@ -183,7 +183,6 @@ test_that("covariates and doses are shifted correctly when t_init != 0", { }) - test_that("covariates_table and doses are shifted correctly when t_init != 0", { cov_table <- data.frame( id = c(1, 1), diff --git a/tests/testthat/test_sim_lagtime.R b/tests/testthat/test_sim_lagtime.R index 91951509..77a8c5f7 100644 --- a/tests/testthat/test_sim_lagtime.R +++ b/tests/testthat/test_sim_lagtime.R @@ -2,21 +2,25 @@ test_that("dose dump after lagtime in correct order in output data", { skip_on_cran() reg <- new_regimen(amt = 500, n = 4, interval = 12, type = 'oral') pars <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83) - mod <- new_ode_model( - code = " - dAdt[0] = -KA * A[0] - dAdt[1] = +KA * A[0] -(CL/V) * A[1] - ", - lagtime = c("TLAG", 0), - obs = list(cmt = 2, scale = "V"), - dose = list(cmt = 1, bioav = 1), - parameters = pars - ) dat <- sim_ode( - ode = mod, + ode = mod_1cmt_oral_lagtime, regimen = reg, parameters = pars, only_obs = FALSE ) - expect_equal(round(dat[dat$t == 12.83 & dat$comp == 1,]$y, 1), c(1.2, 501.2)) + ## Change after RXR-2394: time of TLAG not in dataset anymore unless requested by user in t_obs + ## Before change: expect_equal(round(dat[dat$t == 12.83 & dat$comp == 1,]$y, 1), c(1.2, 501.2)) + ## After change: + expect_equal(nrow(dat[dat$t == 12.83,]), 0) + ## When grid requested by user, lagtime should be visible + dat <- sim_ode( + ode = mod_1cmt_oral_lagtime, + regimen = reg, + parameters = pars, + t_obs = seq(0, 1, 0.01), + only_obs = TRUE + ) + tmp <- dat[dat$t >= 0.82 & dat$t <= 0.85, ] + expect_equal(tmp$t, c(0.82, 0.83, 0.84, 0.85)) + expect_equal(round(tmp$y, 2), c(0, 0, 0.05, 0.10)) }) From 5899fa683c402a158d1dbbe54a87e698f1461a01 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 05:32:31 -0700 Subject: [PATCH 10/20] fix setup.R --- tests/testthat/setup.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index d4e2f51c..906863f0 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -9,7 +9,7 @@ mod_1cmt_oral_lagtime <- new_ode_model( lagtime = c("TLAG", 0), obs = list(cmt = 2, scale = "V"), dose = list(cmt = 1, bioav = 1), - parameters = pars + parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83) ) oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence factor code = " From 90e003b0fe1503215336e31168c42ac81962ee37 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 05:55:40 -0700 Subject: [PATCH 11/20] update sim_core to support lagtime + added test --- R/sim_core.R | 14 ++++++----- man/sim_core.Rd | 16 +++++++++---- tests/testthat/test_sim_core.R | 44 ++++++++++++++++++++++++++++++++++ 3 files changed, 64 insertions(+), 10 deletions(-) diff --git a/R/sim_core.R b/R/sim_core.R index 23c0df5e..6196940d 100644 --- a/R/sim_core.R +++ b/R/sim_core.R @@ -1,24 +1,26 @@ #' Only core function of the simulation function, always just returns observations. #' Mostly useful for estimations / optimal design. Has no checks (for speed)! #' +#' @inheritParams sim #' @param sim_object list with design and simulation parameters -#' @param ode ode -#' @param duplicate_t_obs allow duplicate t_obs in output? E.g. for optimal design calculations when t_obs = c(0,1,2,2,3). Default is FALSE. -#' @param t_init time of initialization of the ODE system. Usually 0. +#' #' @export -#' @return Data frame with simulation results +#' +#' @return data.frame with simulation results +#' sim_core <- function( sim_object = NULL, ode, duplicate_t_obs = FALSE, - t_init = 0 + t_init = 0, + lagtime = c(0) ) { tmp <- ode( A = sim_object$A_init, design = sim_object$design, par = sim_object$p, iov_bins = sim_object$iov_bins, - lagtime = c(0), + lagtime = lagtime, step_size = sim_object$int_step_size ) out <- data.frame(t = tmp$time, y = tmp$obs, obs_type = tmp$obs_type) diff --git a/man/sim_core.Rd b/man/sim_core.Rd index 502bdb81..6e2d35e6 100644 --- a/man/sim_core.Rd +++ b/man/sim_core.Rd @@ -5,19 +5,27 @@ \title{Only core function of the simulation function, always just returns observations. Mostly useful for estimations / optimal design. Has no checks (for speed)!} \usage{ -sim_core(sim_object = NULL, ode, duplicate_t_obs = FALSE, t_init = 0) +sim_core( + sim_object = NULL, + ode, + duplicate_t_obs = FALSE, + t_init = 0, + lagtime = c(0) +) } \arguments{ \item{sim_object}{list with design and simulation parameters} -\item{ode}{ode} +\item{ode}{function describing the ODE system} \item{duplicate_t_obs}{allow duplicate t_obs in output? E.g. for optimal design calculations when t_obs = c(0,1,2,2,3). Default is FALSE.} -\item{t_init}{time of initialization of the ODE system. Usually 0.} +\item{t_init}{initialization time before first dose, default 0.} + +\item{lagtime}{either a value (numeric) or a parameter (character) or NULL.} } \value{ -Data frame with simulation results +data.frame with simulation results } \description{ Only core function of the simulation function, always just returns observations. diff --git a/tests/testthat/test_sim_core.R b/tests/testthat/test_sim_core.R index ce09b429..1e18d2a0 100644 --- a/tests/testthat/test_sim_core.R +++ b/tests/testthat/test_sim_core.R @@ -16,3 +16,47 @@ test_that("sim core works", { } expect_equal(f1(), f2()) }) + +test_that("sim core works for absorption model with lagtime", { + reg <- new_regimen(amt = 100, n = 5, interval = 12, t_inf = 1, type = "infusion") + par <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83) + + ## have to be explicit about t_obs with sim_core! + f_ref <- function() { + res <- sim( + ode = mod_1cmt_oral_lagtime, + regimen = reg, + parameters = par, + only_obs = TRUE, + t_obs=c(0:24) + )$y + return(res) + } + f_core_nolag <- function() { + obj <- sim( + ode = mod_1cmt_oral_lagtime, + regimen = reg, + parameters = par, + only_obs = TRUE, + t_obs=c(0:24), + return_design=TRUE + ) + sim_core(obj, ode = mod_1cmt_oral_lagtime)$y + } + f_core <- function() { + obj <- sim( + ode = mod_1cmt_oral_lagtime, + regimen = reg, + parameters = par, + only_obs = TRUE, + t_obs=c(0:24), + return_design=TRUE + ) + sim_core( + obj, + ode = mod_1cmt_oral_lagtime, + lagtime = c(0.83, 0) + )$y ## Lagtime parameter needed! + } + expect_equal(f_ref(), f_core()) +}) From 4a7468c93ad283bbadb3f6ab050ffd88b7a5a7a2 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 06:48:55 -0700 Subject: [PATCH 12/20] split off parse_lagtime funciton --- R/parse_lagtime.R | 20 +++++++++++ R/sim.R | 9 +---- man/parse_lagtime.Rd | 19 ++++++++++ tests/testthat/test_parse_lagtime.R | 56 +++++++++++++++++++++++++++++ 4 files changed, 96 insertions(+), 8 deletions(-) create mode 100644 R/parse_lagtime.R create mode 100644 man/parse_lagtime.Rd create mode 100644 tests/testthat/test_parse_lagtime.R diff --git a/R/parse_lagtime.R b/R/parse_lagtime.R new file mode 100644 index 00000000..2da62b34 --- /dev/null +++ b/R/parse_lagtime.R @@ -0,0 +1,20 @@ +#' Parse lagtime specified to main sim() function +#' +#' @inheritParams sim +#' +#' @returns a vector of character or numeric values +#' +parse_lagtime <- function( + lagtime, + ode +) { + if(!is.null(attr(ode, "lagtime")) && attr(ode, "lagtime")[1] != "undefined" && attr(ode, "lagtime")[1] != "NULL") { + if(is.null(lagtime)) { # only override from metadata if not specified by user + lagtime <- attr(ode, "lagtime") + } + } + if(is.null(lagtime)) { + lagtime <- c(0) + } + lagtime +} \ No newline at end of file diff --git a/R/sim.R b/R/sim.R index f49c8b2e..569276ef 100644 --- a/R/sim.R +++ b/R/sim.R @@ -147,14 +147,7 @@ sim <- function (ode = NULL, regimen <- merge_regimen(list(regimen, regimen_dupl)) } } - if(!is.null(attr(ode, "lagtime")) && attr(ode, "lagtime")[1] != "undefined" && attr(ode, "lagtime")[1] != "NULL") { - if(is.null(lagtime)) { # only override from metadata if not specified by user - lagtime <- attr(ode, "lagtime") - } - } - if(is.null(lagtime)) { - lagtime <- c(0) - } + lagtime <- parse_lagtime(lagtime, ode) if(!is.null(attr(ode, "dose")$duration_scale)) { regimen <- apply_duration_scale( regimen, diff --git a/man/parse_lagtime.Rd b/man/parse_lagtime.Rd new file mode 100644 index 00000000..658ec85d --- /dev/null +++ b/man/parse_lagtime.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_lagtime.R +\name{parse_lagtime} +\alias{parse_lagtime} +\title{Parse lagtime specified to main sim() function} +\usage{ +parse_lagtime(lagtime, ode) +} +\arguments{ +\item{lagtime}{either a value (numeric) or a parameter (character) or NULL.} + +\item{ode}{function describing the ODE system} +} +\value{ +a vector of character or numeric values +} +\description{ +Parse lagtime specified to main sim() function +} diff --git a/tests/testthat/test_parse_lagtime.R b/tests/testthat/test_parse_lagtime.R new file mode 100644 index 00000000..ecf4997b --- /dev/null +++ b/tests/testthat/test_parse_lagtime.R @@ -0,0 +1,56 @@ +test_that("parse_lagtime works correctly", { + # Test 1: lagtime provided by user, no ode attribute + ode1 <- list() + result1 <- parse_lagtime(lagtime = 5, ode = ode1) + expect_equal(result1, 5) + + # Test 2: no lagtime provided, no ode attribute + ode2 <- list() + result2 <- parse_lagtime(lagtime = NULL, ode = ode2) + expect_equal(result2, c(0)) + + # Test 3: no lagtime provided, ode has lagtime attribute + ode3 <- list() + attr(ode3, "lagtime") <- c(2.5) + result3 <- parse_lagtime(lagtime = NULL, ode = ode3) + expect_equal(result3, c(2.5)) + + # Test 4: lagtime provided by user, ode has lagtime attribute (user takes precedence) + ode4 <- list() + attr(ode4, "lagtime") <- c(3) + result4 <- parse_lagtime(lagtime = 7, ode = ode4) + expect_equal(result4, 7) + + # Test 5: ode lagtime attribute is "undefined" + ode5 <- list() + attr(ode5, "lagtime") <- "undefined" + result5 <- parse_lagtime(lagtime = NULL, ode = ode5) + expect_equal(result5, c(0)) + + # Test 6: ode lagtime attribute is "NULL" + ode6 <- list() + attr(ode6, "lagtime") <- "NULL" + result6 <- parse_lagtime(lagtime = NULL, ode = ode6) + expect_equal(result6, c(0)) + + # Test 7: ode lagtime attribute is NULL + ode7 <- list() + attr(ode7, "lagtime") <- NULL + result7 <- parse_lagtime(lagtime = NULL, ode = ode7) + expect_equal(result7, c(0)) + + # Test 8: vector lagtime inputs + ode8 <- list() + result8 <- parse_lagtime(lagtime = c(1, 2, 3), ode = ode8) + expect_equal(result8, c(1, 2, 3)) + + # Test 9: vector lagtime from ode attribute + ode9 <- list() + attr(ode9, "lagtime") <- c(0.5, 1.5, 2.5) + result9 <- parse_lagtime(lagtime = NULL, ode = ode9) + expect_equal(result9, c(0.5, 1.5, 2.5)) + + # Test 10: character vector lagtime inputs + result10 <- parse_lagtime(lagtime = c("TLAG", 0, 0), ode = list()) + expect_equal(result10, c("TLAG", 0, 0)) +}) From c967795ff71ae36c39043a2093824942a97979a1 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Tue, 8 Jul 2025 01:16:45 -0700 Subject: [PATCH 13/20] rm cran file --- CRAN-SUBMISSION | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 CRAN-SUBMISSION diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION deleted file mode 100644 index 54b4befd..00000000 --- a/CRAN-SUBMISSION +++ /dev/null @@ -1,3 +0,0 @@ -Version: 1.4.1 -Date: 2025-04-16 16:44:08 UTC -SHA: 9419161ff5c39cefefeda5695f91f8c65bb99ab9 From 62bcd77991026917434f7076b01c395797454a1f Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Tue, 8 Jul 2025 01:33:03 -0700 Subject: [PATCH 14/20] cleanup + bump version --- DESCRIPTION | 2 +- R/compile_sim_cpp.R | 32 ------------ inst/cpp/lagtime.utils.cpp | 75 ----------------------------- inst/cpp/lagtime.utils.h | 10 ---- inst/cpp/sim.cpp | 3 +- tests/testthat/test_new_ode_model.R | 5 +- tests/testthat/test_sim.R | 1 - 7 files changed, 5 insertions(+), 123 deletions(-) delete mode 100644 inst/cpp/lagtime.utils.cpp delete mode 100644 inst/cpp/lagtime.utils.h diff --git a/DESCRIPTION b/DESCRIPTION index f26ac97d..fbac24ab 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: PKPDsim Type: Package Title: Tools for Performing Pharmacokinetic-Pharmacodynamic Simulations -Version: 1.4.1 +Version: 1.5.0 Date: 2025-04-09 Authors@R: c( person("Ron", "Keizer", email = "ron@insight-rx.com", role = c("aut", "cre")), diff --git a/R/compile_sim_cpp.R b/R/compile_sim_cpp.R index fb623108..32393448 100644 --- a/R/compile_sim_cpp.R +++ b/R/compile_sim_cpp.R @@ -287,35 +287,3 @@ compile_sim_cpp <- function( cpp = sim_func )) } - -## Copy cpp and header files to temp folder (first remove existing files) -copy_cpp_files_to_tempfolder <- function(sim_func) { - tmp_folder <- tempdir() - existing_files <- c( - dir(tmp_folder, pattern = "\\.cpp$"), - dir(tmp_folder, pattern = "\\.h$"), - dir(tmp_folder, pattern = "\\.o$") - ) - for(f in existing_files) { - unlink(file.path(tmp_folder, f)) - } - main_cpp <- file.path(tmp_folder, "sim.cpp") - if(file.exists(main_cpp)) - unlink(main_cpp) - writeLines(sim_func, main_cpp) - files_to_copy <- setdiff(c( - dir(system.file(package="PKPDsim", "cpp"), pattern = "\\.cpp$"), - dir(system.file(package="PKPDsim", "cpp"), pattern = "\\.h$") - ), "sim.cpp") - for(f in files_to_copy) { - file.copy( - system.file(package="PKPDsim", "cpp", f), - file.path(tmp_folder, f) - ) - } - list( - main_cpp = main_cpp, - folder = tmp_folder - ) - -} \ No newline at end of file diff --git a/inst/cpp/lagtime.utils.cpp b/inst/cpp/lagtime.utils.cpp deleted file mode 100644 index b2d311a4..00000000 --- a/inst/cpp/lagtime.utils.cpp +++ /dev/null @@ -1,75 +0,0 @@ -#include "lagtime.utils.h" - -NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) { - NumericVector lagtime_numeric; - if (TYPEOF(lagtime) == REALSXP) { - lagtime_numeric = as(lagtime); - } else if (TYPEOF(lagtime) == STRSXP) { - CharacterVector lagtime_char = as(lagtime); - lagtime_numeric = NumericVector(lagtime_char.size()); - for (int i = 0; i < lagtime_char.size(); i++) { - String param_name = lagtime_char[i]; - if (parameters.containsElementNamed(param_name.get_cstring())) { - lagtime_numeric[i] = as(parameters[param_name]); - } else { - lagtime_numeric[i] = 0.0; // default value - } - } - } else { - stop("lagtime must be either numeric or character vector"); - } - return(lagtime_numeric); -} - -List apply_lagtime(List design, NumericVector lagtime) { - - List new_design = clone(design); - std::vector times = as >(new_design["t"]); - std::vector evid = as >(new_design["evid"]); - std::vector cmt = as >(new_design["dose_cmt"]); - - // Apply lagtime to dose events (evid == 1) - for(int i = 0; i < times.size(); i++) { - if(evid[i] == 1) { - times[i] += lagtime[cmt[i]-1]; - } - } - - // Create sorted index according to "t" - std::vector indices(times.size()); - std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ... - std::sort(indices.begin(), indices.end(), [×](size_t i1, size_t i2) { - return times[i1] < times[i2]; - }); - - // Reorder all elements in `t` in new_design - std::vector sorted_times(times.size()); - for (size_t i = 0; i < indices.size(); i++) { - sorted_times[i] = times[indices[i]]; - } - new_design["t"] = sorted_times; - - // Sort all other vectors in the design object - for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) { - if (new_design.containsElementNamed(key)) { - SEXP vec = new_design[key]; - if (TYPEOF(vec) == REALSXP) { - std::vector old_vec = as >(vec); - std::vector new_vec(old_vec.size()); - for (size_t i = 0; i < indices.size(); i++) { - new_vec[i] = old_vec[indices[i]]; - } - new_design[key] = new_vec; - } else if (TYPEOF(vec) == INTSXP) { - std::vector old_vec = as >(vec); - std::vector new_vec(old_vec.size()); - for (size_t i = 0; i < indices.size(); i++) { - new_vec[i] = old_vec[indices[i]]; - } - new_design[key] = new_vec; - } - } - } - - return new_design; -} \ No newline at end of file diff --git a/inst/cpp/lagtime.utils.h b/inst/cpp/lagtime.utils.h deleted file mode 100644 index 93c41ac5..00000000 --- a/inst/cpp/lagtime.utils.h +++ /dev/null @@ -1,10 +0,0 @@ -#ifndef LAGTIME_UTILS_H -#define LAGTIME_UTILS_H - -#include -using namespace Rcpp; - -NumericVector lagtime_to_numeric(SEXP lagtime, List parameters); -List apply_lagtime(List design, NumericVector lagtime); - -#endif \ No newline at end of file diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index 3adb0659..30dbdf12 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -189,7 +189,8 @@ List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_ // call ode() once to pre-calculate any initial variables // insert A dAdt state_init set_covariates(0); - + t_prv_dose = times[0]; + // Main call to ODE solver, initialize any variables in ode code ode(A_dum, dAdt_dum, 0); diff --git a/tests/testthat/test_new_ode_model.R b/tests/testthat/test_new_ode_model.R index b51b3ec4..5b5141a6 100644 --- a/tests/testthat/test_new_ode_model.R +++ b/tests/testthat/test_new_ode_model.R @@ -11,7 +11,7 @@ test_that("Double absorption models work appropriately", { bioav = c(0.3, 0.6, 0.9) ), lagtime = c(0, 1.5), - cpp = T + cpp = F ) reg <- new_regimen(amt = 100, n = 12, interval = 12, type = "oral") res <- sim( @@ -25,8 +25,7 @@ test_that("Double absorption models work appropriately", { ), only_obs = F, output_include = list(variables=TRUE), - t_obs = seq(0, 24, .1), - return_design = F + t_obs = seq(0, 24, .1) ) expect_equal(attr(mod, "dose")$duplicate, 2) expect_equal(round(res[res$comp == 1 & res$t == 0.1, ]$y, 2), 27.15) # abs comp 1 diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R index 1be57a6d..4b488bec 100644 --- a/tests/testthat/test_sim.R +++ b/tests/testthat/test_sim.R @@ -38,7 +38,6 @@ test_that("return_event_table=TRUE returns an appropriate event table", { }) test_that("return_event_table=TRUE returns an appropriate event table with covariate", { - mod_1cmt_iv <- new_ode_model("pk_1cmt_iv") covs <- list(CRCL = new_covariate(value = c(70, 80), t = c(0, 24)), WT = new_covariate(70)) evtab2 <- sim_ode( mod_1cmt_iv, From e7579fa565a24c07584c6e9767209e78192dd948 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Tue, 8 Jul 2025 01:33:47 -0700 Subject: [PATCH 15/20] update date --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index fbac24ab..4ff62243 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: PKPDsim Type: Package Title: Tools for Performing Pharmacokinetic-Pharmacodynamic Simulations Version: 1.5.0 -Date: 2025-04-09 +Date: 2025-07-08 Authors@R: c( person("Ron", "Keizer", email = "ron@insight-rx.com", role = c("aut", "cre")), person("Jasmine", "Hughes", email = "jasmine@insight-rx.com", role = "aut"), From 1fc77444161533939c159f332cb4cccc53d3604a Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Wed, 9 Jul 2025 02:22:35 -0700 Subject: [PATCH 16/20] minor fix for models with lagtime on infusion --- inst/cpp/sim.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/inst/cpp/sim.cpp b/inst/cpp/sim.cpp index 30dbdf12..66ae6848 100644 --- a/inst/cpp/sim.cpp +++ b/inst/cpp/sim.cpp @@ -89,6 +89,7 @@ List apply_lagtime(List design, NumericVector tlag, int n_comp) { List new_design = clone(design); std::vector times = as >(new_design["t"]); std::vector evid = as >(new_design["evid"]); + std::vector rate = as >(new_design["rate"]); std::vector cmt = as >(new_design["dose_cmt"]); NumericVector lagtime; @@ -107,7 +108,7 @@ List apply_lagtime(List design, NumericVector tlag, int n_comp) { // Apply lagtime to dose events (evid == 1) for(int i = 0; i < times.size(); i++) { - if(evid[i] == 1) { + if(evid[i] == 1 | (evid[i] == 2 & rate[i] != 0)) { // dose events or infusion stop events times[i] += lagtime[cmt[i]-1]; } } From dafa5c3dcd25b92fac60be8fa606939410c9613d Mon Sep 17 00:00:00 2001 From: roninsightrx Date: Wed, 9 Jul 2025 04:03:48 -0700 Subject: [PATCH 17/20] Update R/parse_lagtime.R Co-authored-by: Jasmine Hughes <43552465+jasmineirx@users.noreply.github.com> --- R/parse_lagtime.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/parse_lagtime.R b/R/parse_lagtime.R index 2da62b34..5a000435 100644 --- a/R/parse_lagtime.R +++ b/R/parse_lagtime.R @@ -8,13 +8,13 @@ parse_lagtime <- function( lagtime, ode ) { - if(!is.null(attr(ode, "lagtime")) && attr(ode, "lagtime")[1] != "undefined" && attr(ode, "lagtime")[1] != "NULL") { - if(is.null(lagtime)) { # only override from metadata if not specified by user - lagtime <- attr(ode, "lagtime") - } + lagtime_ode <- attr(ode, "lagtime") + # override from ode if not specified by user and defined in ode + if(is.null(lagtime) && !is.null(lagtime_ode) && !lagtime_ode %in% c("NULL", "undefined")) { + lagtime <- lagtime_ode } if(is.null(lagtime)) { - lagtime <- c(0) + lagtime <- 0 } lagtime } \ No newline at end of file From e79239bb7a1becd46829ca55519626f458332380 Mon Sep 17 00:00:00 2001 From: roninsightrx Date: Wed, 9 Jul 2025 04:03:58 -0700 Subject: [PATCH 18/20] Update R/sim.R Co-authored-by: Jasmine Hughes <43552465+jasmineirx@users.noreply.github.com> --- R/sim.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/sim.R b/R/sim.R index 569276ef..df662a88 100644 --- a/R/sim.R +++ b/R/sim.R @@ -600,7 +600,7 @@ sim <- function (ode = NULL, ## remove dose-times from regimen ## but leave ones that were also in t_obs - ## also leave extra obs when bolus, because this may be needed (controllec below using `extra_t_obs`) + ## also leave extra obs when bolus, because this may be needed (controlled below using `extra_t_obs`) dose_times <- design_i[design_i$evid > 0, ]$t comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, comb$y, sep="_"))),] From 513ab136f07875aa6b434d6ece15761f1feee3aa Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Wed, 9 Jul 2025 04:16:03 -0700 Subject: [PATCH 19/20] misc updates from peer review --- R/parse_lagtime.R | 9 +- R/sim.R | 2 +- tests/testthat/test_parse_lagtime.R | 129 +++++++++++++++++----------- tests/testthat/test_sim_core.R | 11 --- 4 files changed, 87 insertions(+), 64 deletions(-) diff --git a/R/parse_lagtime.R b/R/parse_lagtime.R index 5a000435..b36cd443 100644 --- a/R/parse_lagtime.R +++ b/R/parse_lagtime.R @@ -6,7 +6,8 @@ #' parse_lagtime <- function( lagtime, - ode + ode, + parameters ) { lagtime_ode <- attr(ode, "lagtime") # override from ode if not specified by user and defined in ode @@ -16,5 +17,11 @@ parse_lagtime <- function( if(is.null(lagtime)) { lagtime <- 0 } + if(inherits(lagtime, "character")) { + idx <- grep("[a-zA-Z]", lagtime) # only pick character names, not "0" + if(! all(lagtime[idx] %in% names(parameters))) { + warning("Lagtime parameter(s) not found. Please check model and parameters.") + } + } lagtime } \ No newline at end of file diff --git a/R/sim.R b/R/sim.R index df662a88..8524f6a1 100644 --- a/R/sim.R +++ b/R/sim.R @@ -147,7 +147,7 @@ sim <- function (ode = NULL, regimen <- merge_regimen(list(regimen, regimen_dupl)) } } - lagtime <- parse_lagtime(lagtime, ode) + lagtime <- parse_lagtime(lagtime, ode, parameters) if(!is.null(attr(ode, "dose")$duration_scale)) { regimen <- apply_duration_scale( regimen, diff --git a/tests/testthat/test_parse_lagtime.R b/tests/testthat/test_parse_lagtime.R index ecf4997b..43dfb03b 100644 --- a/tests/testthat/test_parse_lagtime.R +++ b/tests/testthat/test_parse_lagtime.R @@ -1,56 +1,83 @@ -test_that("parse_lagtime works correctly", { - # Test 1: lagtime provided by user, no ode attribute +test_that("parse_lagtime handles NULL lagtime input", { + # Mock ODE object without lagtime attribute + ode <- list() + parameters <- list(CL = 5, V = 50) + + result <- parse_lagtime(NULL, ode, parameters) + expect_equal(result, 0) +}) + +test_that("parse_lagtime uses ODE lagtime when user lagtime is NULL", { + # Mock ODE object with lagtime attribute + ode <- list() + attr(ode, "lagtime") <- "TLAG" + parameters <- list(CL = 5, V = 50, TLAG = 0.6) + + result <- parse_lagtime(NULL, ode, parameters) + expect_equal(result, "TLAG") +}) + +test_that("parse_lagtime ignores ODE lagtime when user provides lagtime", { + # Mock ODE object with lagtime attribute + ode <- list() + attr(ode, "lagtime") <- "TLAG" + parameters <- list(CL = 5, V = 50, TLAG = 0.6) + + result <- parse_lagtime(1.5, ode, parameters) + expect_equal(result, 1.5) +}) + +test_that("parse_lagtime handles character lagtime with valid parameters", { + ode <- list() + parameters <- list(CL = 5, V = 50, TLAG = 0.6) + lagtime <- c("TLAG", "0", "0") + + result <- parse_lagtime(lagtime, ode, parameters) + expect_equal(result, c("TLAG", "0", "0")) +}) + +test_that("parse_lagtime warns when character lagtime parameters not found", { + ode <- list() + parameters <- list(CL = 5, V = 50) + lagtime <- c("TLAG", "0", "0") + + expect_warning( + parse_lagtime(lagtime, ode, parameters), + "Lagtime parameter\\(s\\) not found. Please check model and parameters." + ) +}) + +test_that("parse_lagtime handles numeric lagtime", { + ode <- list() + parameters <- list(CL = 5, V = 50) + lagtime <- 1.5 + + result <- parse_lagtime(lagtime, ode, parameters) + expect_equal(result, 1.5) +}) + +test_that("parse_lagtime ignores NULL and undefined ODE lagtime", { + # Test NULL lagtime in ODE ode1 <- list() - result1 <- parse_lagtime(lagtime = 5, ode = ode1) - expect_equal(result1, 5) + attr(ode1, "lagtime") <- "NULL" + parameters <- list(CL = 5, V = 50) + + result1 <- parse_lagtime(NULL, ode1, parameters) + expect_equal(result1, 0) - # Test 2: no lagtime provided, no ode attribute + # Test undefined lagtime in ODE ode2 <- list() - result2 <- parse_lagtime(lagtime = NULL, ode = ode2) - expect_equal(result2, c(0)) - - # Test 3: no lagtime provided, ode has lagtime attribute - ode3 <- list() - attr(ode3, "lagtime") <- c(2.5) - result3 <- parse_lagtime(lagtime = NULL, ode = ode3) - expect_equal(result3, c(2.5)) - - # Test 4: lagtime provided by user, ode has lagtime attribute (user takes precedence) - ode4 <- list() - attr(ode4, "lagtime") <- c(3) - result4 <- parse_lagtime(lagtime = 7, ode = ode4) - expect_equal(result4, 7) - - # Test 5: ode lagtime attribute is "undefined" - ode5 <- list() - attr(ode5, "lagtime") <- "undefined" - result5 <- parse_lagtime(lagtime = NULL, ode = ode5) - expect_equal(result5, c(0)) - - # Test 6: ode lagtime attribute is "NULL" - ode6 <- list() - attr(ode6, "lagtime") <- "NULL" - result6 <- parse_lagtime(lagtime = NULL, ode = ode6) - expect_equal(result6, c(0)) - - # Test 7: ode lagtime attribute is NULL - ode7 <- list() - attr(ode7, "lagtime") <- NULL - result7 <- parse_lagtime(lagtime = NULL, ode = ode7) - expect_equal(result7, c(0)) - - # Test 8: vector lagtime inputs - ode8 <- list() - result8 <- parse_lagtime(lagtime = c(1, 2, 3), ode = ode8) - expect_equal(result8, c(1, 2, 3)) - - # Test 9: vector lagtime from ode attribute - ode9 <- list() - attr(ode9, "lagtime") <- c(0.5, 1.5, 2.5) - result9 <- parse_lagtime(lagtime = NULL, ode = ode9) - expect_equal(result9, c(0.5, 1.5, 2.5)) + attr(ode2, "lagtime") <- "undefined" + + result2 <- parse_lagtime(NULL, ode2, parameters) + expect_equal(result2, 0) +}) - # Test 10: character vector lagtime inputs - result10 <- parse_lagtime(lagtime = c("TLAG", 0, 0), ode = list()) - expect_equal(result10, c("TLAG", 0, 0)) +test_that("parse_lagtime handles mixed character and numeric lagtime", { + ode <- list() + parameters <- list(CL = 5, V = 50, TLAG1 = 0.6, TLAG2 = 1.2) + lagtime <- c("TLAG1", "0", "TLAG2") + + result <- parse_lagtime(lagtime, ode, parameters) + expect_equal(result, c("TLAG1", "0", "TLAG2")) }) diff --git a/tests/testthat/test_sim_core.R b/tests/testthat/test_sim_core.R index 1e18d2a0..fd81ca48 100644 --- a/tests/testthat/test_sim_core.R +++ b/tests/testthat/test_sim_core.R @@ -32,17 +32,6 @@ test_that("sim core works for absorption model with lagtime", { )$y return(res) } - f_core_nolag <- function() { - obj <- sim( - ode = mod_1cmt_oral_lagtime, - regimen = reg, - parameters = par, - only_obs = TRUE, - t_obs=c(0:24), - return_design=TRUE - ) - sim_core(obj, ode = mod_1cmt_oral_lagtime)$y - } f_core <- function() { obj <- sim( ode = mod_1cmt_oral_lagtime, From a8654b32acacd3b2db882df7bc308292282d00a1 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Wed, 9 Jul 2025 04:40:22 -0700 Subject: [PATCH 20/20] minor fix --- R/parse_lagtime.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/parse_lagtime.R b/R/parse_lagtime.R index b36cd443..d87e88f2 100644 --- a/R/parse_lagtime.R +++ b/R/parse_lagtime.R @@ -11,7 +11,7 @@ parse_lagtime <- function( ) { lagtime_ode <- attr(ode, "lagtime") # override from ode if not specified by user and defined in ode - if(is.null(lagtime) && !is.null(lagtime_ode) && !lagtime_ode %in% c("NULL", "undefined")) { + if(is.null(lagtime) && !is.null(lagtime_ode) && lagtime_ode[1] != "NULL" && lagtime_ode[1] != "undefined") { lagtime <- lagtime_ode } if(is.null(lagtime)) {