Nothing
# Copyright (C) 2013 - 2023 Metrum Research Group
#
# This file is part of mrgsolve.
#
# mrgsolve is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# mrgsolve is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with mrgsolve. If not, see <http://www.gnu.org/licenses/>.
#' @include modspec.R
NULL
#' Read a model specification file
#'
#' \code{mread} reads and parses the \code{mrgsolve} model specification file,
#' builds the model, and returns a model object for simulation.
#' \code{mread_cache} does the same, but caches the compilation result for
#' later use.
#'
#'
#' @param model model name
#' @param project location of the model specification file an any
#' headers to be included; see also the discussion about model; this argument
#' can be set via \code{options()}
#' library under details as well as the \code{\link{modlib}} help topic
#' @param file the full file name (with extension, but without path)
#' where the model is specified
#' @param soloc the directory location where the model shared object is built
#' and stored; see details; this argument can be set via \code{options()};
#' if the directory does not exist, `mread` will attempt to create it.
#' @param code a character string with model specification code to be
#' used instead of a model file
#' @param ignore.stdout passed to system call for compiling model
#' @param raw if TRUE, return a list of raw output
#' @param compile logical; if \code{TRUE}, the model will be built
#' @param check.bounds check boundaries of parameter list
#' @param audit check the model specification file for errors
#' @param warn logical; if \code{TRUE}, print warning messages that may arise
#' @param udll use unique name for shared object
#' @param quiet don't print messages when compiling
#' @param preclean logical; if \code{TRUE}, compilation artifacts are
#' cleaned up first
#' @param recover if \code{TRUE}, a list of build will be returned in case
#' the model shared object fails to compile; use this option to and
#' the returned object to collect information assist in debugging
#' @param capture a character vector or comma-separated string of additional
#' model variables to capture; these variables will be added to the capture
#' list for the current call to \code{\link{mread}} only
#' @param ... passed to \code{\link[mrgsolve]{update}}; also arguments passed
#' to mread from \code{\link{mread_cache}}.
#'
#' @details
#' The \code{model} argument is required. For typical use,
#' the \code{file} argument is omitted and the value
#' for \code{file} is generated from the value for \code{model}.
#' To determine the source file name, \code{mrgsolve} will look for
#' a file extension in \code{model}. A file extension is
#' assumed when it finds a period followed by one to three alpha-numeric
#' characters at the end of the string (e.g. \code{mymodel.txt} but not
#' \code{my.model}). If no file extension is found, the extension \code{.cpp}
#' is assumed (e.g. \code{file} is \code{<model-name>.cpp}). If a file
#' extension is found, \code{file} is \code{<model-name>}.
#'
#' Best practice is to avoid using \code{.} in \code{model} unless
#' you are using \code{model} to point to the model specification
#' file name. Otherwise, use \code{\link{mread_file}}.
#'
#' Use the \code{soloc} argument to specify a directory location for building
#' the model. This is the location where the model shared object will be
#' stored on disk. The default is a temporary directory, so compilation
#' artifacts are lost when R restarts when the default is used. Changing
#' \code{soloc} to a persistent directory location will preserve those
#' artifacts across R restarts. Also, if simulation from a single model is
#' being done in separate processes on separate compute nodes, it might be
#' necessary to store these compilation artifacts in a local directory
#' to make them accessible to the different nodes. If the \code{soloc}
#' directory does not exist, `mread` will attempt to create it.
#'
#' Similarly, using \code{mread_cache} will cache results in the temporary
#' directory and the cache cannot be accessed after the R process is
#' restarted.
#'
#' @section Model Library:
#'
#' \code{mrgsolve} comes bundled with several precoded PK, PK/PD, and
#' other systems models that are accessible via the \code{mread} interface.
#'
#' Models available in the library include:
#'
#' \itemize{
#' \item PK models: \code{pk1cmt}, \code{pk2cmt}, \code{pk3cmt},
#' \code{pk1}, \code{pk2}, \code{popex}, \code{tmdd}
#' \item PKPD models: \code{irm1}, \code{irm2}, \code{irm3}, \code{irm4},
#' \code{emax}, \code{effect}
#' \item Other models: \code{viral1}, \code{viral2}
#' }
#'
#' When the library model is accessed, \code{mrgsolve} will compile and load
#' the model as you would for any other model. It is only necessary to
#' reference the correct model name and point the \code{project} argument
#' to the \code{mrgsolve} model library location via \code{\link{modlib}}.
#'
#' For more details, see \code{\link{modlib_pk}}, \code{\link{modlib_pkpd}},
#' \code{\link{modlib_tmdd}}, \code{\link{modlib_viral}}, and
#' \code{\link{modlib_details}} for more information about the state
#' variables and parameters in each model.
#'
#' @examples
#'
#' \dontrun{
#' code <- '
#' $PARAM CL = 1, VC = 5
#' $CMT CENT
#' $ODE dxdt_CENT = -(CL/VC)*CENT;
#' '
#'
#' mod <- mcode("ex_mread", code)
#'
#' mod
#'
#' mod %>% init(CENT=1000) %>% mrgsim %>% plot
#'
#'
#' mod <- mread("irm3", modlib())
#'
#' mod
#'
#' # if the model is in the file mymodel.cpp
#' mod <- mread("mymodel")
#'
#' # if the model is in the file mymodel.txt
#' mod <- mread(file = "mymodel.txt")
#'
#' or
#'
#' mod <- mread_file("mymodel.txt")
#'
#'
#' }
#'
#' @seealso \code{\link{mcode}}, \code{\link{mcode_cache}}
#'
#' @export
mread <- function(model, project = getOption("mrgsolve.project", getwd()),
code = NULL, file = NULL,
udll = TRUE, ignore.stdout = TRUE,
raw = FALSE, compile = TRUE, audit = TRUE,
quiet = getOption("mrgsolve_mread_quiet", FALSE),
check.bounds = FALSE, warn = TRUE,
soloc = getOption("mrgsolve.soloc", tempdir()),
capture = NULL,
preclean = FALSE, recover = FALSE, ...) {
if(!identical(basename(model), as.character(model))) {
project <- dirname(model)
model <- basename(model)
}
quiet <- isTRUE(quiet)
warn <- isTRUE(warn) & (!quiet)
if(missing(model) & !missing(file)) {
model <- file
}
build <- new_build(
file = file, model = model, project = project,
soloc = soloc, code = code, preclean = preclean,
udll = udll, recover = recover
)
model <- build[["model"]]
## Read the model spec and parse:
is_rmd <- grepl("\\.[rR]md$", build[["modfile"]])[1]
if(is_rmd) {
spec <- modelparse_rmd(readLines(build[["modfile"]],warn = FALSE))
} else {
spec <- modelparse(readLines(build[["modfile"]], warn = FALSE))
}
## Block name aliases
incoming_block_names <- names(spec)
names(spec) <- toupper(names(spec))
names(spec) <- gsub("DES", "ODE", names(spec), fixed = TRUE)
names(spec) <- gsub("POST", "TABLE", names(spec), fixed = TRUE)
names(spec) <- gsub("ERROR", "TABLE", names(spec), fixed = TRUE)
names(spec) <- gsub("^PK$", "MAIN", names(spec), fixed = FALSE)
## Expand partial matches
index <- pmatch(names(spec), block_list, duplicates.ok = TRUE)
names(spec) <- ifelse(is.na(index), names(spec), block_list[index])
## Do a check on what we found in the spec
check_spec_contents(names(spec), warn = warn,...)
## Pull out the settings and ENV now
## We might be passing parse settings in here ...
SET <- tolist(dump_opts(spec[["SET"]]))
ENV <- eval_ENV_block(spec[["ENV"]],build[["project"]])
if("SET" %in% names(spec)) spec[["SET"]] <- ""
if("ENV" %in% names(spec)) spec[["ENV"]] <- ""
# mread parse env ----
mread.env <- parse_env(
spec,
incoming_block_names,
build,
ENV
)
# Find globals ----
# The main sections that need R processing:
# NOTE: this must happen prior to handling the various blocks
# `captured` items are injected into `$CAPTURE`; this can be changed but
# for now leave this here
spec <- move_global2(
spec,
mread.env,
build
)
# Find cpp objects with dot syntax; saved to mread.env$cpp_dot ----
find_cpp_dot(spec, mread.env)
# Parse blocks ----
# Each block gets assigned a class to dispatch the handler function
# Also, we use a position attribute so we know
# where we are when we're handling the block
specClass <- paste0("spec", names(spec))
for(i in seq_along(spec)) {
spec[[i]] <- structure(
.Data = spec[[i]],
class = specClass[i],
pos = i
)
}
# Call the handler for each block
spec <- lapply(spec, handle_spec_block, env = mread.env)
# collect -----
# TODO: move this to the plugin handler
plugin <- get_plugins(spec[["PLUGIN"]], mread.env)
param <- as.list(do.call("c",unname(mread.env[["param"]])))
fixed <- as.list(do.call("c",unname(mread.env[["fixed"]])))
init <- as.list(do.call("c",unname(mread.env[["init"]])))
ode <- do.call("c", unname(mread.env[["ode"]]))
namespace <- do.call("c", mread.env[["namespace"]])
omega <- omat(do.call("c", nonull.list(mread.env[["omega"]])))
sigma <- smat(do.call("c", nonull.list(mread.env[["sigma"]])))
if(isTRUE(SET[["collapse_omega"]])) {
omega <- collapse_matrix(omega)
}
if(isTRUE(SET[["collapse_sigma"]])) {
sigma <- collapse_matrix(sigma)
}
annot_list_maybe <- nonull.list(mread.env$annot)
if (!length(annot_list_maybe)) {
annot <- tibble()
} else {
annot <- bind_rows(annot_list_maybe)
}
# capture ----
capture_more <- capture
capture <- unlist(nonull.list(mread.env[["capture"]]))
capture <- .ren.create(capture)
capture <- .ren.sanitize(capture)
annot <- capture_param(annot, .ren.new(capture))
# Collect potential multiples
subr <- collect_subr(spec)
table <- unlist(spec[names(spec)=="TABLE"], use.names = FALSE)
spec[["ODE"]] <- unlist(spec[names(spec)=="ODE"], use.names = FALSE)
# TODO: deprecate audit argument
mread.env[["audit_dadt"]] <-
isTRUE(audit) && isTRUE(mread.env[["audit_dadt"]])
# Look for compartments we're dosing into: F/ALAG/D/R
# and add them to CMTN
dosing <- dosing_cmts(spec[["MAIN"]], names(init))
SET[["CMTN"]] <- c(spec[["CMTN"]], dosing)
# new model object ----
x <- new(
"mrgmod",
model = model,
soloc = build[["soloc"]],
package = build[["package"]],
project = build[["project"]],
fixed = fixed,
advan = subr[["advan"]],
trans = subr[["trans"]],
omega = omega,
sigma = sigma,
param = as.param(param),
init = as.init(init),
funs = funs_create(model),
envir = ENV,
capture = .ren.chr(capture),
plugin = names(plugin),
modfile = basename(build[["modfile"]])
)
# updates -----
# First update with what we found in the model specification file
x <- update(x, data = SET, open = TRUE, strict = FALSE)
# Next, update with what the user passed in as arguments
args <- list(...)
x <- update(x, data = args, open = TRUE, strict = FALSE)
x <- store_annot(x, annot)
# Arguments in $SET that will be passed to mrgsim
x <- set_simargs(x, SET)
# Modify SS compartments
x <- set_ss_cmt(x, SET[["ss_cmt"]])
# model r defs ----
# These are the various #define statements
# that go at the top of the .cpp.cpp file
rd <- generate_rdefs(
pars = Pars(x),
cmt = Cmt(x),
ode_func(x),
main_func(x),
table_func(x),
config_func(x),
model = model(x),
omats = omat(x),
smats = smat(x),
set = SET,
check.bounds = check.bounds,
plugin = plugin,
dbsyms = args[["dbsyms"]]
)
# Handle plugins ----
# Virtual compartments
if("VCMT" %in% names(spec)) {
what <- which(names(spec)=="VCMT")
vcmt <- unique(names(unlist(mread.env[["init"]][what])))
spec[["ODE"]] <- c(spec[["ODE"]], paste0("dxdt_", vcmt, "=0;"))
}
# Handle nm-vars plugin
if("nm-vars" %in% names(plugin)) {
nmv <- find_nm_vars(spec)
dfs <- generate_nmdefs(nmv)
rd <- c(rd, dfs)
plugin[["nm-vars"]][["nm-def"]] <- dfs
build[["nm-vars"]] <- nmv
audit_nm_vars(
spec,
param = param,
init = init,
build = build,
nmv = nmv,
env = mread.env
)
}
# autodec
if("autodec" %in% names(plugin)) {
auto_blocks <- c("PREAMBLE", "MAIN", "PRED", "ODE", "TABLE")
auto_skip <- cvec_cs(ENV[["MRGSOLVE_AUTODEC_SKIP"]])
autov <- autodec_vars(spec, blocks = auto_blocks)
autov <- autodec_clean(
autov,
rdefs = rd,
build = build,
skip = auto_skip
)
autodec_nm_vars(autov, mread.env)
autodec_save(autov, build, mread.env)
mread.env[["autodec"]] <- autodec_namespace(build, mread.env)
}
# Rcpp
if("Rcpp" %in% names(plugin)) {
spec <- global_rcpp(spec)
}
# mrgx
if("mrgx" %in% names(plugin)) {
toglob <- wrap_namespace("Rcpp::Environment _env;", NULL)
topream <- "_env = mrgx::get_envir(self);"
spec[["PREAMBLE"]] <- c(topream, spec[["PREAMBLE"]])
spec[["GLOBAL"]] <- c(toglob, spec[["GLOBAL"]])
}
# more captures ----
# Process @etas 1:n -----
x <- capture_etas(x, mread.env)
# Process capture passed into mread -----
if(is.character(capture_more)) {
valid_capture <- get_valid_capture(
param = param, omega = omega, sigma = sigma, build = build,
mread.env = mread.env
)
if(identical(capture_more, "(everything)")) {
capture_more <- valid_capture[valid_capture != "."]
}
capture_vars <- .ren.create(capture_more)
capture_vars <- .ren.sanitize(capture_vars)
if(!all(capture_vars[["old"]] %in% valid_capture)) {
bad <- setdiff(capture_vars[["old"]], valid_capture)
for(b in bad) {
msg <- glue(" - item `{b}` does not exist in model `{build$model}`")
message(msg)
}
stop(
"all requested `capture` variables must exist in the model",
call. = FALSE
)
}
x <- update_capture(x, .ren.chr(capture_vars))
build$preclean <- TRUE
}
# Check mod ----
check_pkmodel(x, subr, spec)
check_globals(mread.env[["move_global"]], Cmt(x))
check_cpp_dot(mread.env, x)
check_sim_eta_eps_n(x, spec)
# This must come after audit
# TODO: harmonize with other audit process
if(!has_name("ODE", spec)) {
spec[["ODE"]] <- "DXDTZERO();"
} else {
if(mread.env[["audit_dadt"]] && !mread.env[["using_nm-vars"]]) {
spec[["ODE"]] <- c(spec[["ODE"]], ode)
audit_spec(x, spec, warn = warn)
}
}
# set up shlib ----
# lock some of this down so we can check order later
x@code <- readLines(build[["modfile"]], warn=FALSE)
x@shlib[["covariates"]] <- mread.env[["covariates"]]
x@shlib[["param_tag"]] <- mread.env[["param_tag"]]
x@shlib[["cpp_variables"]] <- build$cpp_variables
inc <- spec[["INCLUDE"]]
if(is.null(inc)) inc <- character(0)
x@shlib[["include"]] <- inc
x@shlib[["nm_import"]] <- mread.env[["nm_import"]]
x@shlib[["source"]] <- file.path(build[["soloc"]],build[["compfile"]])
x@shlib[["md5"]] <- build[["md5"]]
# build----
# In soloc directory
cwd <- getwd()
setwd(build[["soloc"]])
to_restore <- set_up_env(
x = plugin,
clink = c(project(x), SET[["clink"]]),
CXX = ENV[["PKG_CXXFLAGS"]]
)
on.exit({
setwd(cwd)
do_restore(to_restore)
})
incl <- function(x) paste0('#include "', x, '"')
header_file <- paste0(build[["model"]], "-mread-header.h")
dbs <- NULL
if(isTRUE(args[["dbsyms"]])) {
dbs <- debug_symbols(names(init(x)))
}
cat(
paste0("// Source MD5: ", build[["md5"]]),
"\n// PLUGINS:",
plugin_code(plugin),
## This should get moved to rd
"\n// FIXED:",
fixed_parameters(fixed, SET[["fixed_type"]]),
"\n// NAMESPACES:",
namespace,
"\n// MODEL HEADER FILES:",
incl("mrgsolv.h"),
incl("modelheader.h"),
"\n// INCLUDE databox functions:",
incl("databox_cpp.h"),
"\n// USING plugins:",
plugin_using(plugin),
"\n// INCLUDES:",
form_includes(spec[["INCLUDE"]]),
"\n// GLOBAL CODE BLOCK:",
"// GLOBAL VARS FROM BLOCKS & TYPEDEFS:",
"// DECLARED BY USER",
mread.env[["global"]],
"// DECLARED VIA AUTODEC",
mread.env[["autodec"]],
"\n// GLOBAL START USER CODE:",
spec[["GLOBAL"]],
"\n// DEFS:",
rd,
"",
sep="\n", file = header_file)
## Write the model code to temporary file
temp_write <- tempfile()
def.con <- file(temp_write, open="w")
cat(
paste0("// Source MD5: ", build[["md5"]], "\n"),
incl(header_file),
"\n// PREAMBLE CODE BLOCK:",
"__BEGIN_config__",
spec[["PREAMBLE"]],
"__END_config__",
"\n// MAIN CODE BLOCK:",
"__BEGIN_main__",
dbs[["cmt"]],
spec[["MAIN"]],
advtr(x@advan,x@trans),
"__END_main__",
"\n// DIFFERENTIAL EQUATIONS:",
"__BEGIN_ode__",
dbs[["ode"]],
spec[["ODE"]],
"__END_ode__",
"\n// TABLE CODE BLOCK:",
"__BEGIN_table__",
dbs[["cmt"]],
table,
spec[["PRED"]],
write_capture(names(x@capture)),
"__END_table__",
"",
sep="\n", file=def.con)
close(def.con)
## this gets written in soloc
#write_build_env(build)
write_win_def(x)
same <- check_and_copy(
from = temp_write,
to = build[["compfile"]]
)
if(!compile) return(x)
if(ignore.stdout & !quiet) {
message("Building ", model(x), " ... ", appendLF = FALSE)
}
# Wait at least 2 sec since last compile
safe_wait(x)
cleanso(x)
# Compile the model
# The shared object is model-mread-source.cpp
out <- suppressWarnings(build_exec(build))
comp_success <- out[["status"]]==0 & file.exists(build[["compout"]])
if(!comp_success) {
if(ignore.stdout) message("error.\n", appendLF = FALSE)
return(build_failed(out, build, x, spec, ignore.stdout))
}
if(ignore.stdout) {
if(!quiet) message("done.\n", appendLF = FALSE)
} else {
out <- build_output_cleanup(out, build)
cat(out[["stdout"]], sep = "\n")
}
# Rename the shared object to unique name
# e.g model2340239403.so
z <- file.copy(build[["compout"]], sodll(x))
dyn.load(sodll(x))
stopifnot(all_loaded(x))
x <- compiled(x,TRUE)
return(x)
}
#' @rdname mread
#' @export
mread_cache <- function(model = NULL,
project = getOption("mrgsolve.project", getwd()),
file = paste0(model, ".cpp"),
code = NULL,
soloc = getOption("mrgsolve.soloc", tempdir()),
quiet = FALSE,
preclean = FALSE,
capture = NULL, ...) {
if (!identical(basename(model), as.character(model))) {
project <- dirname(model)
model <- basename(model)
}
if(is.character(capture)) preclean <- TRUE
build <- new_build(file, model, project, soloc, code, preclean)
cache_file <- file.path(build$soloc, "mrgmod_cache.RDS")
## If the cache file doesn't exist, build and return
te <- file.exists(cache_file) & !preclean
t0 <- t1 <- t2 <- t3 <- t4 <- FALSE
if(te) {
x <- readRDS(cache_file)
t0 <- is.mrgmod(x)
if(t1 <- is.character(x@shlib$md5)) {
t2 <- x@shlib$md5 == build$md5
}
t3 <- file.exists(sodll(x))
t4 <- length(x@shlib[["include"]])==0
}
if(all(t0,t1,t2,t3,t4,te)) {
if(!quiet) message("Loading model from cache.")
loadso(x)
return(update(x,...,strict=FALSE))
}
x <- mread(
build$model, project, soloc=soloc, quiet=quiet,
file = basename(build$modfile), capture = capture, ...
)
saveRDS(x,file=cache_file)
return(x)
}
#' @export
#' @rdname mread
mread_file <- function(file, ...) {
model <- tools::file_path_sans_ext(file)
mread(model = model, file = file, ...)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.