#' Create new ODE model
#'
#' @param model model name from model library
#' @param code C++ code specifying ODE system
#' @param pk_code C++ code called at any event
#' @param dose_code C++ code called at dose event only
#' @param file file containing C++ code
#' @param func R function to be used with deSolve library
#' @param state_init vector of state init
#' @param parameters list or vector of parameter values
#' @param reparametrization list of parameters with definitions that reparametrize the linear PK model to a 1-, 2- o4 3-compartment PK with standardized parametrization.
#' @param mixture for mixture models, provide a list of the parameter associated with the mixture and it's possible values and probabilities (of the first value), e.g. `list(CL = list(value = c(10, 20), probability = 0.3)`.
#' @param units list or vector of parameter units
#' @param size size of state vector for model. Size will be extracted automatically from supplied code, use this argument to override.
#' @param lagtime lag time
#' @param obs list with "scale": character string with definition for scale, e.g. "V" or "V*(WT/70)". If NULL, scale defaults to 1., and "cmt" the observation compartment
#' @param dose specify default dose compartment, e.g. list(cmt = 1)
#' @param covariates specify covariates, either as a character vector or a list. if specified as list, it allows use of timevarying covariates (see `new_covariate()` function for more info)
#' @param declare_variables declare variables
#' @param fixed parameters that should not have iiv added.
#' @param iiv inter-individual variability, can optionally be added to library
#' @param iov inter-occasion variability, can optionally be added to library
#' @param development Information about the model development population, can optionally be added to library
#' @param omega_matrix variance-covariance matrix for inter-individual variability, can optionally be added to library
#' @param ruv residual variability, can optionally be added to library
#' @param ltbs log-transform both sides. Not used in simulations, only for fitting (sets attribute `ltbs`).
#' @param misc a list of miscellaneous model metadata
#' @param cmt_mapping list indicating which administration routes apply to which compartments. Example: `list("oral" = 1, "infusion" = 2)`
#' @param int_step_size step size for integrator. Can be pre-specified for model, to override default for `sim_ode()`
#' @param default_parameters population or specific patient values, can optionally be added to library
#' @param cpp_show_code show generated C++ code
#' @param package package name when saving as package
#' @param test_file optional test file to be included with package
#' @param install install package after compilation?
#' @param folder base folder name to create package in
#' @param lib_location install into folder (`--library` argument)
#' @param verbose show more output
#' @param as_is use C-code as-is, don't substitute line-endings or shift indices
#' @param nonmem add NONMEM code as attribute to model object
#' @param comments comments for model
#' @param version number of library
#' @param quiet passed on to `system2` as setting for stderr and stdout; how to
#' output cmd line output. Default (`""`) is R console, NULL or FALSE discards.
#' TRUE captures the output and saves as a file.
#' @param definition optional, filename for the JSON file the full definition
#' for the model. The definition file will be stored as `definition.json` in the
#' resulting package.
#' @export
#' @return If package name is NULL, returns the model object. Otherwise has no
#' return value.
new_ode_model <- function (model = NULL,
code = NULL,
pk_code = NULL,
dose_code = NULL,
file = NULL,
func = NULL,
state_init = NULL,
parameters = NULL,
reparametrization = NULL,
mixture = NULL,
units = NULL,
size = NULL,
lagtime = NULL,
obs = list("cmt" = 1, "scale" = 1),
dose = list("cmt" = 1),
covariates = NULL,
declare_variables = NULL,
iiv = NULL,
iov = NULL,
development = NULL,
omega_matrix = NULL,
ruv = NULL,
ltbs = NULL,
misc = NULL,
cmt_mapping = NULL,
int_step_size = NULL,
default_parameters = NULL,
fixed = NULL,
cpp_show_code = FALSE,
package = NULL,
test_file = NULL,
install = TRUE,
folder = NULL,
lib_location = NULL,
verbose = FALSE,
as_is = FALSE,
nonmem = NULL,
comments = NULL,
version = "0.1.0",
quiet = "",
definition = NULL
) {
if (is.null(model) & is.null(code) & is.null(file) & is.null(func)) {
stop(paste0("Either a model name (from the PKPDsim library), ODE code, an R function, or a file containing code for the ODE system have to be supplied to this function. The following models are available:\n ", model_library()))
}
if (!is.null(func)) { # R function supplied, use deSolve approach
if(!inherits(func, "function")) {
stop("The object specified in 'func' is not a valid R function.")
}
sim_out <- func
attr(sim_out, "cpp") <- FALSE
} else {
if (!is.null(model)) {
mod_def <- model_library(model)
code <- mod_def$code
obs <- mod_def$obs
dose <- mod_def$dose
if(!is.null(mod_def$size)) {
size <- mod_def$size
} else {
size <- get_ode_model_size(mod_def$code)
}
}
code_init_text <- ""
code <- gsub("(#).*?\\n", "\n", code) # remove comments
if(inherits(state_init, "character")) {
code <- gsub("dadt", "dAdt", code)
code <- gsub("DADT", "dAdt", code)
code <- gsub("^\\n", "", code)
code <- gsub("^(\\s)*", "", code)
state_init <- gsub("(#).*?\\n", "\n", state_init) # remove comments
state_init <- gsub("\\n", ";\n", state_init) # make sure each line has a line ending
code_init_text <- paste0(state_init, ";\n")
}
if(is.null(size) && !is.null(code)) {
size <- get_ode_model_size(code)
if(size == 0) {
stop("Sorry, no ODE system detected. Please specify ODEs using 'dAdt[...] = ...'")
}
}
code_init_text <- shift_state_indices(code_init_text, -1)
if(is.null(size)) {
size <- get_ode_model_size(code)
}
## IIV
if(!is.null(omega_matrix)) {
if(!is_positive_definite(omega_matrix)) {
stop("Specified omega matrix is not positive definite.")
}
}
## IOV
use_iov <- FALSE
if(!is.null(iov)) {
check_iov_specification(iov, code, pk_code)
# add parameters
for(i in rev(seq(iov$cv))) {
txt <- paste0(" kappa_", names(iov$cv)[i], " = 1e-6;\ \n")
if(length(grep(paste0("kappa_", names(iov$cv)[i]), code)) > 0) {
for(j in 1:iov$n_bins) {
txt_occ <- paste0(" if(t >= iov_bin[",(j-1),"] && t < iov_bin[",j,"]) { kappa_", names(iov$cv)[i], " = kappa_", names(iov$cv)[i], "_", j, " + 1e-6; } \ \n")
txt <- paste0(txt, txt_occ)
par_tmp <- paste0("kappa_", names(iov$cv)[i], "_", j)
if(! par_tmp %in% parameters) {
parameters <- c(parameters, par_tmp)
}
default_parameters[[par_tmp]] <- 0
}
code <- paste0(txt, code)
}
txt <- paste0(" kappa_", names(iov$cv)[i], " = 1e-6;\ \n")
if(length(grep(paste0("kappa_", names(iov$cv)[i]), pk_code)) > 0) {
for(j in 1:iov$n_bins) {
txt_occ <- paste0(" if(times[i] >= iov_bin[",(j-1),"] && times[i] < iov_bin[",j,"]) { kappa_", names(iov$cv)[i], " = kappa_", names(iov$cv)[i], "_", j, "; } \ \n")
txt <- paste0(txt, txt_occ)
par_tmp <- paste0("kappa_", names(iov$cv)[i], "_", j)
if(! par_tmp %in% parameters) {
parameters <- c(parameters, par_tmp)
}
default_parameters[[par_tmp]] <- 0
}
pk_code <- paste0(txt, pk_code)
}
var_tmp <- paste0("kappa_", names(iov$cv)[i])
if(! var_tmp %in% declare_variables) {
declare_variables <- c(declare_variables, var_tmp)
}
}
} else {
iov <- list(n_bins = 1) # dummy
}
if(is.null(parameters)) {
comb_code <- code
if(!is.null(pk_code)) {
comb_code <- paste(code, pk_code, sep = "\n")
}
parameters <- get_parameters_from_code(comb_code, state_init, declare_variables)
} else {
if(inherits(parameters, "list")) {
parameters <- names(parameters)
}
}
check_mixture_model(mixture, parameters)
obs <- check_obs_input(obs)
if(is.null(dose$cmt)) { dose$cmt <- 1 }
cov_names <- NULL
if(!is.null(covariates)) {
if(inherits(covariates, "character")) {
cov_names <- covariates
}
if(inherits(covariates, "list")) {
cov_names <- names(covariates)
}
}
variables <- declare_variables
if(!is.null(cov_names) && !is.null(variables)) {
if(length(intersect(cov_names, variables))) {
stop("Covariates and specified variables cannot overlap.")
}
}
if(!is.null(cov_names)) {
declare_variables <- c(declare_variables,
cov_names,
paste0(cov_names, "_0"),
paste0("gr_", cov_names),
paste0("t_prv_", cov_names))
}
code_orig <- code
code <- cleanup_code(code)
pk_code <- cleanup_code(pk_code)
dose_code <- cleanup_code(dose_code)
compile <- TRUE
if(!is.null(package)) { # don't compile if saving as library
compile = FALSE
}
cmp <- compile_sim_cpp(
code = code,
pk_code = pk_code,
dose_code = dose_code,
size = size,
p = parameters,
cpp_show_code = cpp_show_code,
code_init = code_init_text,
state_init = state_init,
declare_variables = declare_variables,
variables = variables,
covariates = covariates,
obs = obs,
dose = dose,
verbose = verbose,
compile = compile,
as_is = as_is,
iov = iov
)
reqd <- parameters
if(length(grep("cov_", reqd)) > 0) {
reqd <- reqd[-grep("cov_", reqd)]
}
if(!is.null(declare_variables)) {
reqd <- reqd[!reqd %in% declare_variables]
}
if(!is.null(cov_names)) {
reqd <- reqd[!reqd %in% cov_names]
}
if(exists("sim_wrapper_cpp", envir = globalenv())) {
if(compile) {
sim_out <- get("sim_wrapper_cpp")
}
rm("sim_wrapper_cpp", envir = globalenv())
}
if(compile) {
attr(sim_out, "code") <- code
if(!is.null(pk_code)) {
attr(sim_out, "pk_code") <- pk_code
}
attr(sim_out, "parameters") <- reqd
attr(sim_out, "mixture") <- mixture
attr(sim_out, "reparametrization") <- reparametrization
attr(sim_out, "covariates") <- cov_names
attr(sim_out, "variables") <- variables
attr(sim_out, "fixed") <- fixed
attr(sim_out, "cpp") <- TRUE
attr(sim_out, "size") <- size
attr(sim_out, "obs") <- obs
attr(sim_out, "dose") <- dose
attr(sim_out, "lagtime") <- lagtime
attr(sim_out, "ltbs") <- ltbs
attr(sim_out, "misc") <- misc
attr(sim_out, "cmt_mapping") <- cmt_mapping
attr(sim_out, "iov") <- iov
attr(sim_out, "development") <- development
attr(sim_out, "comments") <- paste0("\n", as.character(paste0(paste0(" - ", comments), collapse = "\n")))
if(!is.null(int_step_size)) {
attr(sim_out, "int_step_size") <- int_step_size
}
class(sim_out) <- c("PKPDsim", class(sim_out))
return(sim_out)
}
## save to library, if requested:
if(!is.null(package)) {
if(is.null(folder)) {
folder <- tempdir()
}
## Copy template files
new_folder <- file.path(folder, package)
templ_folder <- file.path(system.file(package="PKPDsim"), "template")
templ_folder <- file.path(templ_folder, dir(templ_folder))
if(!file.exists(new_folder)) {
dir.create(new_folder)
}
file.copy(
from = templ_folder,
to = new_folder,
overwrite = TRUE,
recursive = TRUE,
copy.mode = FALSE
)
## copy definition JSON, if model created from definition file.
if(!is.null(definition) && file.exists(definition)) {
dir.create(file.path(new_folder, "inst"))
copy_result <- file.copy(
from = definition,
to = file.path(new_folder, "inst", "definition.json5")
)
if(verbose) {
if(isTRUE(copy_result)) {
message("* Adding model definition to R package.")
} else {
message("* Could not add model definition to R package.")
}
}
}
## Write new source file
fileConn <- file(file.path(new_folder, "src", "sim_wrapper_cpp.cpp"))
writeLines(cmp$cpp, fileConn)
close(fileConn)
## Replace module name and other info
if(is.null(pk_code)) { pk_code <- "" }
if(is.null(dose_code)) { dose_code <- "" }
lagtime <- vector_to_R_code(lagtime)
dose$bioav <- bioavailability_to_R_code(dose$bioav)
if(is.null(size)) { size <- "1" }
if(is.null(ltbs)) { ltbs <- FALSE }
if(is.null(state_init)) {
state_init <- "NULL"
} else {
state_init <- add_quotes(state_init)
}
if(is.null(nonmem) || length(nonmem) == 0) { nonmem <- "NULL" }
if(is.null(int_step_size)) { int_step_size <- "NULL" }
pars <- vector_to_R_code(reqd)
covs <- vector_to_R_code(cov_names)
fixed <- vector_to_R_code(fixed)
vars <- vector_to_R_code(variables)
repl <- matrix(c("\\[MODULE\\]", package,
"\\[N_COMP\\]", size,
"\\[OBS_COMP\\]", obs$cmt,
"\\[DOSE_COMP\\]", dose$cmt,
"\\[OBS_SCALE\\]", obs$scale,
"\\[OBS_VARIABLE\\]", paste0(deparse(obs$variable), collapse = ""),
"\\[DOSE_BIOAV\\]", dose$bioav,
"\\[DOSE_DUPLICATE\\]", paste0(deparse(dose$duplicate), collapse = ""),
"\\[CODE\\]", code,
"\\[PK_CODE\\]", pk_code,
"\\[DOSE_CODE\\]", dose_code,
"\\[STATE_INIT\\]", state_init,
"\\[PARS\\]", pars,
"\\[REPARAM\\]", paste0(deparse(reparametrization), collapse = ""),
"\\[MIXTURE\\]", paste0(deparse(mixture), collapse = ""),
"\\[VARS\\]", vars,
"\\[COVS\\]", covs,
"\\[FIXED\\]", fixed,
"\\[LAGTIME\\]", lagtime,
"\\[DURATION_SCALE\\]", deparse(dose$duration_scale),
"\\[USE_IOV\\]", as.character(use_iov),
"\\[IOV\\]", PKPDsim::print_list(iov, FALSE),
"\\[DEVELOPMENT\\]", print_list(development, FALSE),
"\\[LTBS\\]", as.character(ltbs),
"\\[MISC\\]", paste0(deparse(misc), collapse = ""),
"\\[CMT_MAPPING\\]", paste0(deparse(cmt_mapping), collapse = ""),
"\\[INT_STEP_SIZE\\]", as.character(int_step_size),
"\\[COMMENTS\\]", paste0("\n", as.character(paste0(paste0(" - ", comments), collapse = "\n"))),
"\\[NONMEM\\]", as.character(nonmem)
), ncol=2, byrow=TRUE)
if(verbose) {
print(repl)
}
search_replace_in_file(paste0(new_folder, "/R/model.R"), repl[,1], repl[,2])
if(is.null(iiv)) { iiv <- "" } else {
iiv <- PKPDsim::print_list(iiv, FALSE)
}
if(is.null(iov)) { iov <- "" } else {
iov <- PKPDsim::print_list(iov, FALSE)
}
if(is.null(ruv)) { ruv <- "" } else {
ruv <- PKPDsim::print_list(ruv, FALSE)
}
if(is.null(omega_matrix)) {
omega_matrix <- ""
} else {
omega_matrix <- paste0("c(", paste(omega_matrix, collapse = ", "), ")")
}
if(is.null(default_parameters)) { default_parameters <- "" } else {
default_parameters <- PKPDsim::print_list(default_parameters, FALSE)
}
if(is.null(units)) { units <- "''" } else {
units <- PKPDsim::print_list(units, wrapper = TRUE)
}
search_replace_in_file(file.path(new_folder, "R", "iiv.R"), "\\[IIV\\]", iiv)
search_replace_in_file(file.path(new_folder, "R", "iov.R"), "\\[IOV\\]", iov)
search_replace_in_file(file.path(new_folder, "R", "omega_matrix.R"), "\\[OMEGA_MATRIX\\]", omega_matrix)
search_replace_in_file(file.path(new_folder, "R", "parameters.R"), c("\\[PARAMETERS\\]", "\\[UNITS\\]"), c(default_parameters, units))
search_replace_in_file(file.path(new_folder, "R", "fixed.R"), "\\[FIXED\\]", fixed)
search_replace_in_file(file.path(new_folder, "R", "ruv.R"), "\\[RUV\\]", ruv)
search_replace_in_file(file.path(new_folder, "DESCRIPTION"), "\\[MODULE\\]", package)
search_replace_in_file(file.path(new_folder, "DESCRIPTION"), "\\[VERSION\\]", version)
search_replace_in_file(file.path(new_folder, "NAMESPACE"), "\\[MODULE\\]", package)
search_replace_in_file(file.path(new_folder, "man", "modulename-package.Rd"), "\\[MODULE\\]", package)
file.rename(
file.path(new_folder, "man", "modulename-package.Rd"),
file.path(new_folder, "man", paste0(package, ".Rd"))
)
## copy test file into package
if(!is.null(test_file)) {
if(file.exists(test_file[1])) {
t_dir <- file.path(new_folder, "tests")
if(!file.exists(t_dir)) dir.create(t_dir)
file.copy(test_file[1], file.path(t_dir, paste0(package, ".R")))
} else {
warning("Specified test file not found.")
}
}
## Compile / build / install
if(file.exists(file.path(new_folder, "R", "RcppExports.R"))) {
file.remove(file.path(new_folder, "R", "RcppExports.R"))
}
if(file.exists(file.path(new_folder, "src", "RcppExports.cpp"))) {
file.remove(file.path(new_folder, "src", "RcppExports.cpp"))
}
Rcpp::compileAttributes(new_folder)
cmd <- file.path(Sys.getenv("R_HOME"), "bin", "R")
if (install) { # install into R
lib_location_arg <- ""
if(!is.null(lib_location)) {
lib_location_arg <- paste0("--library=", lib_location)
}
# Run R CMD INSTALL with appropriate settings
args <- c("CMD", "INSTALL",
lib_location_arg,
"--no-docs",
"--no-demo",
"--no-help",
"--no-multiarch",
"--with-keep.source",
"--pkglock",
"--no-staged-install",
shQuote(normalizePath(file.path(folder, package)))
)
system2(cmd, args, stdout = quiet, stderr = quiet)
} else {
# Run R CMD BUILD to zip file
args <- c("CMD", "build", shQuote(normalizePath(file.path(folder, package))))
system2(cmd, args, stdout = quiet, stderr = quiet)
pkg_file <- paste0(package, "_", version, ".tar.gz")
if(file.exists(pkg_file)) {
message(paste0("Package built in: ", file.path(getwd(), pkg_file)))
} else {
message("Package not created.")
}
}
}
return()
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.