#' check if system is ready for model run
check_model <- function(ip) {
if (!is(ip, "isopath")) stop ("can only check model ready for an isopath", call. = FALSE)
message(sprintf("Checking model for %d scenario(s)...", nrow(ip$parameters)))
tryCatch({
# do a sample calculation for each parameter set
ip$parameters %>%
group_by(!!!purrr::map(names(ip$parameters), rlang::sym)) %>%
do({
# safety checks - will through errors if needed
if (nrow(.) != 1) stop("there seem to be multiple identical run scenarios, please make sure each set of parameters is unique", call. = FALSE)
dx_var <- get_ode_matrix(ip, eval = T, param = .)
if (any(is.na(dx_var$`dx/dt`))) {
dx_text <- get_ode_matrix(ip, eval = F)$`dx/dt`
stop("some derivatives could not be computed, check the expressions: ",
paste(dx_var$x[is.na(dx_var$`dx/dt`)], "=", dx_text[is.na(dx_var$`dx/dt`)]) %>%
paste(collapse = ", "), call. = FALSE)
}
tibble()
})
},
# error handling
error = function(e) {
stop("encountered the following error during pre-check of the model run: <<< ", e$message, " >>>. If this is a missing object error, please make sure that all referenced columns are defined in the parameter set and that scenarios are unique. Currently available to the system are: ", names(ip$parameters) %>% paste(collapse = ", "), call. = FALSE)
})
return(invisible(ip))
}
#' compile the ODE expression for the isopath
get_ode_expression <- function(ip) {
exp_text <- paste0("list(",
(ip %>% get_ode_matrix() %>% mutate(exp = paste(x, "=", `dx/dt`)))$exp %>%
paste(collapse = ", "), ")")
interp(lazy(x), x = parse(text = exp_text, keep.source = F, n = NULL)[[1]])
}
#' get ode derivative
get_ode_function <- function(ip) {
# lazy expression for optimized derivative calculation
exp_lazy <- ip %>% get_ode_expression()
# vector of variable component pools
pools <- (ip %>% get_component_matrix() %>% filter(variable == T))$component
# derivatives function
function(t, y, p) {
# current ODE state
ode_state <- c(as.list(y), p)
# evaluate the expression for the variable dx/dt
dx_var <- lazy_eval(exp_lazy, data = ode_state)
# check if there are any missing values other computational problems
if ( any(sapply(dx_var, length) == 0) || any(is.nan(unlist(dx_var))) )
stop("some derivatives could not be computed: ",
dx_var[sapply(dx_var, length) == 0 | is.nan(unlist(dx_var))] %>% names() %>% paste(collapse = ", "),
" (state variables: ", paste(names(y), "=", y) %>% paste(collapse = ", "), ")",
call. = FALSE)
# make sure not running out of any component or running into infinity problems (with simple loop for speed)
for (pool in pools){
if (ode_state[[pool]] > 1e100)
stop("overflowing pool: ", pool, " (size: ", ode_state[[pool]], " | dx/dt: ", signif(dx_var[[pool]],3), ")",
". Check your system to make sure that none of the component pools overflows (reaches a size of 10^100). ",
"State variables: ", paste(names(y), "=", y) %>% paste(collapse = ", "),
call. = FALSE)
else if (ode_state[[pool]] + dx_var[[pool]] < 0)
stop("depleted the following pool: ", pool, " (size: ", ode_state[[pool]], " | dx/dt: ", signif(dx_var[[pool]],3), ")",
". Please make sure that none of the component pools runs out completely. ",
"State variables: ", paste(names(y), "=", y) %>% paste(collapse = ", "),
call. = FALSE)
}
# all dx/dt, default value 0 for all non modified variables (=parameters)
dx <- modifyList(setNames(as.list(rep(0, length(y))), names(y)), dx_var)
# convert to correct list format
return(list(unlist(dx)))
}
}
#' Run the reaction model
#'
#' Run it for a certain number of time steps.
#'
#' @param time_steps the number of time steps, can be a number or expression (referring to a parameter in the isopath)
#' @param ... additional parameters passed on to the \link{ode} solver
#' @param make_state_var vector of parameters that should be included as state variables (so they can be changed as part of any special events passed to \code{...}, see \link{ode} for detail on the \code{events} parameter). All variable components' pools and isotopic compositions are always inclduded as state parameters and don't need to be added explicitly here.
#' @export
run_model <- function(ip, time_steps, ..., make_state_var = c()) {
if (!is(ip, "isopath")) stop ("can only run model for an isopath", call. = FALSE)
# model ready checks
if (missing(time_steps))
stop("time_steps is required when calling 'run_model' (either a number or an expression referencing a parameter)")
check_model(ip)
# info
message(sprintf("Running model for %d scenario(s)...", nrow(ip$parameters)))
# derivatives function
func <- get_ode_function(ip)
# steps evaluation
steps_exp <- lazy(time_steps)
# state variables
state_vars <- c(get_ode_matrix(ip)$x, make_state_var) %>% unique()
constants <- setdiff(names(ip$parameters), state_vars)
# run each scenario by grouping by each row
result <-
ip$parameters %>%
group_by(!!!purrr::map(names(ip$parameters), rlang::sym)) %>%
do({
# time steps
times <- seq(0, lazy_eval(steps_exp, data = .), by = 1)
# attempt to solve ODE
sln <- tibble()
tryCatch({
sln <- ode(y = unlist(.[state_vars]), times = times, func = func,
parms = as.list(.[constants]), ...)
},
error = function(e) {
message("WARNING: encountered the following error while trying to solve one parameter set (skipping to the next set of parameters):\n<<< ", e$message, " >>>.\nParameters: ", paste(names(unlist(.)), "=", unlist(.)) %>% paste(collapse = ", "))
})
# output
as.data.frame(sln)
}) %>%
ungroup()
if (nrow(result) == 0) stop("None of the scenarios could be computed successfully", call. = FALSE)
return(select(result, !!!c("time",c(state_vars, constants))))
}
#' Run model to steady state
#'
#' Reports the time required to reach steady-state (in model time, with dt=1) and creates <variable>.t0 columns for each state variable. Uses \link{runsteady} internally for greatest flexibility although \link{stode} can work for some reaction systems and would be a little faster.
#'
#' @param ... additional parameters passed on to the steady state ODE (\link{runsteady}) solver, use \link{runsteady} parameter \code{verbose = TRUE} for additional output during computations, use parameters \code{rtol}, \code{atol} and \code{ctol} to adjust steady-state tolerances
#' @export
run_steady_state <- function(ip, ...) {
if (!is(ip, "isopath")) stop ("can only run model for an isopath", call. = FALSE)
# model ready checks
check_model(ip)
# info
message(sprintf("Looking for steady-state for %d scenario(s)...", nrow(ip$parameters)))
# ode function
func <- get_ode_function(ip)
# state variables
state_vars <- get_ode_matrix(ip)$x
state_vars_t0 <- paste0(state_vars, "_t0")
constants <- setdiff(names(ip$parameters), state_vars)
# run each scenario by grouping by each row
result <-
ip$parameters %>%
group_by(!!!purrr::map(names(ip$parameters), rlang::sym)) %>%
do({
sln <- tibble()
# attempt to solve ODE
tryCatch({
timing <-
system.time(
out <- runsteady(y = unlist(.[state_vars]), func = func, times = c(0, Inf),
parms = as.list(.[constants]), ...)
)
if (attributes(out)$steady) {
sln <- bind_cols(
tibble(time = attributes(out)$time),
rename(., !!!setNames(state_vars, state_vars_t0))[state_vars_t0], # t0 values
as_tibble(as.list(out$y)) # solutions
)
} else {
message("WARNING: did not reach steady state (elapsed time: ", signif(timing[['elapsed']], 5), "s) for one parameter set (skipping to the next set of parameters). Run with verbose = TRUE for additional detail.\nParameters: ", paste(names(unlist(.)), "=", unlist(.)) %>% paste(collapse = ", "))
}
},
error = function(e) {
message("WARNING: encountered the following error while trying to find steady state for one reaction system (skipping to the next set of parameters):\n<<< ", e$message, " >>>.\nParameters: ",
paste(names(unlist(.)), "=", unlist(.)) %>% paste(collapse = ", "))
})
# output
sln
}) %>%
ungroup()
if (nrow(result) == 0) stop("None of the scenarios could be run to steady-state.", call. = FALSE)
return(select(result, !!!c("time", state_vars, state_vars_t0, constants)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.