# R/dopri.R In dde: Solve Delay Differential Equations

#### Documented in dopridopri5dopri853dopri_continuedopri_interpolateylag

##' Integrate an ODE or DDE with dopri.
##'
##' Like \code{deSolve::lsoda}, this function has \emph{many}
##' arguments.  This is far from ideal, and I would welcome any
##' approach for simplifying it a bit.
##'
##' The options \code{return_by_column}, \code{return_initial},
##' \code{return_time}, \code{return_output_with_y} exist because
##' these options all carry out modifications of the data at the end
##' of solving the ODE and this can incur a small but measurable cost.
##' When solving an ODE repeatedly (e.g., in the context of an MCMC or
##' optimisation) it may be useful to do as little as possible.  For
##' simple problems this can save around 5-10\% of the total
##' computational time (especially the transpose).  The shorthand
##' option \code{return_minimal} will set all to \code{FALSE} when
##' used.
##'
##' @title Integrate ODE/DDE with dopri
##'
##' @param y Initial conditions for the integration
##'
##' @param times Times where output is needed.  Unlike \code{deSolve}
##'   we won't actually stop at these times, but instead interpolate
##'   back to get the result.
##'
##' @param func Function to integrate.  Can be an R function of
##'   arguments \code{t, y, parms}, returning a numeric vector, or it
##'   can be the name or address of a C function with arguments
##'   \code{size_t n, double t, const double *y, double *dydt, void *data}.
##'
##' @param parms Parameters to pass through to the derivatives.
##'
##' @param ... Dummy arguments - nothing is allowed here, but this
##'   means that all further arguments \emph{must} be specified by
##'   name (not order) so I can easily reorder them later on.
##'
##' @param n_out Number of "output" variables (not differential
##'   equation variables) to compute via the routine \code{output}.
##'
##' @param output The output routine; either an R function taking
##'   arguments \code{t, y, parms} or the name/address of a C function
##'   taking arguments \code{size_t n, double t, const double *y,
##'   size_t n_out, double *out, void *data}.
##'
##' @param rtol The per-step relative tolerance.  The total accuracy
##'   will be less than this.
##'
##' @param atol The per-step absolute tolerance.
##'
##' @param step_size_min The minimum step size.  The actual minimum
##'   used will be the largest of the absolute value of this
##'   \code{step_size_min} or \code{.Machine$double.eps}. If the ##' integration attempts to make a step smaller than this, it will ##' throw an error by default, stopping the integration (note that ##' this differs from the treatment of \code{hmin} in ##' \code{deSolve::lsoda}). See \code{allow_step_size_min} to change ##' this behaviour. ##' ##' @param step_size_max The largest step size. By default there is ##' no maximum step size (Inf) so the solver can take as large a ##' step as it wants to. If you have short-lived fluctuations in ##' your rhs that the solver may skip over by accident, then specify ##' a smaller maximum step size here (or use \code{tcrit} below). ##' ##' @param step_size_initial The initial step size. By default the ##' integrator will guess the step size automatically, but one can ##' be given here instead. ##' ##' @param step_max_n The maximum number of steps allowed. If the ##' solver takes more steps than this it will throw an error. Note ##' the number of evaluations of \code{func} will be about 6 times ##' the number of steps (or 11 times if using \code{method = ##' "dopri853"}). ##' ##' @param step_size_min_allow Logical, indicating if when a step size ##' is driven down to \code{step_size_min} we should allow it to ##' proceed. This is the behaviour in of \code{hmin} in ##' \code{deSolve::lsoda}. ##' ##' @param tcrit An optional vector of critical times that the solver ##' must stop at (rather than interpolating over). This can include ##' an end time that we can't go past, or points within the ##' integration that must be stopped at exactly (for example cases ##' where the derivatives change abruptly). Note that this differs ##' from the interpretation of this parameter in deSolve; there ##' \code{tcrit} is a single time that integration may not go past ##' -- with dde we never go past the final time, and this is just ##' for times that fall \emph{within} the range of times in ##' \code{times}. ##' ##' @param event_time Vector of times to fire events listed in ##' \code{event_function} at ##' ##' @param event_function Function to fire at events. For R models ##' (\code{func} is an R function and \code{dllname} is empty), this ##' must be either a single R function (same function for all ##' events) or a \code{list} of R functions. For C models, this ##' must be a singe C function (same requirements as \code{func} or ##' \code{output} or a list/vector of these as appropriate). ##' ##' @param method The integration method to use, as a string. The ##' supported methods are \code{"dopri5"} (5th order method with 4th ##' order dense output) and \code{"dopri853"} (8th order method with ##' 7th order output and embedded 5th and 3rd order schemes). ##' Alternatively, use the functions \code{dopri5} or ##' \code{dopri853} which simply sets this argument. ##' ##' @param stiff_check How often to check that the problem has become ##' stiff. If zero, then the problem is never checked, and if ##' positive then the problem is checked every \code{stiff_check} ##' accepted steps. The actual check is based off the algorithm in ##' Hairer's implementation of the solvers and may be overly strict, ##' especially for delay equations with the 853 method (in my ##' limited experience with it). ##' ##' @param verbose Be verbose, and print information about each step. ##' This may be useful for learning about models that misbehave. ##' Valid values are \code{TRUE} (enable debugging) or \code{FALSE} ##' (disable debugging) or use one of \code{dopri:::VERBOSE_QUIET}, ##' \code{dopri:::VERBOSE_STEP} or \code{VERBOSE:::VERBOSE_EVAL}. ##' If an R function is provided as the argument \code{callback} ##' then this function will also be called at each step or ##' evaluation (see below for details). ##' ##' @param callback Callback function that can be used to make verbose ##' output more useful. This can be used to return more information ##' about the evaluation as it proceeds, generally as information ##' printed to the screen. The function must accept arguments ##' \code{t}, \code{y} and \code{dydt}. See Details for further ##' information. ##' ##' @param n_history Number of history points to retain. This needs ##' to be greater than zero for delay differential equations to ##' work. Alternatively, this may be greater than zero to return ##' model outputs that can be inspected later. ##' ##' @param grow_history Logical indicating if history should be grown ##' during the simulation. If \code{FALSE} (the default) then when ##' history is used it is overwritten as needed (so only the most ##' recent \code{n_history} elements are saved. This may require ##' some tuning so that you have enough history to run your ##' simulation (i.e. to the longest delay) or an error will be ##' thrown when it underflows. The required history length will ##' vary with your delay sizes and with the timestep for dopri. If ##' \code{TRUE}, then history will grow as the buffer is exhausted. ##' The growth is geometric, so every time it reaches the end of the ##' buffer it will increase by a factor of about 1.6 (see the ##' \code{ring} documentation). This may consume more memory than ##' necessary, but may be useful where you don't want to care about ##' picking the history length carefully. ##' ##' @param return_history Logical indicating if history should be ##' returned alongside the output or discarded. By default, history ##' is retained if \code{n_history} is greater than 0, but that ##' might change (and may not be desirable unless you plan on ##' actually using it). ##' ##' @param dllname Name of the shared library (without extension) to ##' find the function \code{func} (and \code{output} if given) in ##' the case where \code{func} refers to compiled function. ##' ##' @param parms_are_real Logical, indicating if \code{parms} should ##' be treated as vector of doubles by \code{func} (when it is a ##' compiled function). If \code{TRUE} (the default), then ##' \code{REAL(parms)}, which is \code{double*} is passed through. ##' If \code{FALSE} then if \code{params} is an externalptr type ##' (\code{EXTPTRSXP}) we pass through the result of ##' \code{R_ExternalPtrAddr}, otherwise we pass \code{params} ##' through unmodified as a \code{SEXP}. In the last case, in your ##' target function you will need to include \code{<Rinternals.h>}, ##' cast to \code{SEXP} and then pull it apart using the R API (or ##' Rcpp). ##' ##' @param ynames Logical, indicating if the output should be named ##' following the names of the input vector \code{y}. ##' Alternatively, if \code{ynames} is a character vector of the ##' same length as \code{y}, these will be used as the output names. ##' ##' @param outnames An optional character vector, used when ##' \code{n_out} is greater than 0, to name the model output matrix. ##' ##' @param return_by_column Logical, indicating if the output should be ##' returned organised by column (rather than row). This incurs a ##' slight cost for transposing the matrices. If you can work with ##' matrices that are transposed relative to \code{deSolve}, then ##' set this to \code{FALSE}. ##' ##' @param return_initial Logical, indicating if the output should ##' include the initial conditions. Specifying \code{FALSE} avoids ##' binding this onto the output. ##' ##' @param return_time Logical, indicating if a row (or column if ##' \code{return_by_column} is \code{TRUE}) representing time is included. ##' If \code{FALSE}, this is not added. ##' ##' @param return_output_with_y Logical, indicating if the output ##' should be bound together with the returned matrix \code{y} (as ##' it is with \code{deSolve}). If \code{FALSE}, then output will ##' be returned as the attribute \code{output}. ##' ##' @param return_statistics Logical, indicating if statistics about ##' the run should be included. If \code{TRUE}, then an integer ##' vector containing the number of target evaluations, steps, ##' accepted steps and rejected steps is returned (the vector is ##' named). ##' ##' @param return_minimal Shorthand option - if set to \code{TRUE} ##' then it sets all of \code{return_by_column}, ##' \code{return_initial}, \code{return_time}, ##' \code{return_output_with_y} to \code{FALSE} ##' ##' @param restartable Logical, indicating if the problem should be ##' restartable. If \code{TRUE}, then the return value of an ##' integration can be passed to \code{dopri_restart} to continue ##' the integration after arbitrary changes to the state or the ##' parameters. Note that when using delay models, the integrator ##' is fairly naive about how abrupt changes in the state space are ##' dealt with, and may perform very badly with \code{method = ##' "dopri853"} which assumes a fairly smooth problem. Note that ##' this is really only useful for delay differential equations ##' where you want to keep the history but make changes to the ##' parameters or to the state vector while keeping the history of ##' the problem so far. ##' ##' @return At present the return value is transposed relative to ##' deSolve. This might change in future. ##' ##' @export ##' ##' @seealso \code{\link{dopri_interpolate}} which can be used to ##' efficiently sample from output of \code{dopri}, and the package ##' vignette which shows in more detail how to solve delay ##' differential equations and to use compiled objective functions. ##' ##' @section Verbose output and callbacks: ##' ##' Debugging a failed integration can be difficult, but \code{dopri} ##' provides a couple of tools to get more information about where a ##' failure might have occurred. Most simply, one can pass ##' \code{verbose = TRUE} which will print information about the ##' time and the step size at each point just before the step is ##' stated. Passing in \code{verbose = dde:::VERBOSE_EVAL} will ##' print information just before every evaluation of the target ##' function (there are several evaluations per step). ##' ##' However, this does not provide information about the state just ##' before failure. To get that, one must provide a \code{callback} ##' function - this is an R function that will be called just before ##' a step or evaluation (based on the value of the \code{verbose} ##' argument) in place of the default print. Define a callback ##' function with arguments \code{t}, \code{h} and \code{y} where ##' \code{t} is the time (beginning of a step or location of an ##' evaluation), \code{h} is the step size (or \code{NA} for an ##' evaluation) and \code{y} is the state at the point of the step ##' or evaluation. Your callback function can do anything - you can ##' print to the screen (using \code{cat} or \code{message}), you ##' can store results using a closure and \code{<<-} or you could ##' conditionally use a \code{browser()} call to debug ##' interactively. However, it is not possible for the callback to ##' affect the solution short of throwing an error and interrupting ##' it. See the Examples for an example of use. ##' ##' @examples ##' ##' # The lorenz attractor: ##' lorenz <- function(t, y, p) { ##' sigma <- p[[1L]] ##' R <- p[[2L]] ##' b <- p[[3L]] ##' c(sigma * (y[[2L]] - y[[1L]]), ##' R * y[[1L]] - y[[2L]] - y[[1L]] * y[[3L]], ##' -b * y[[3L]] + y[[1L]] * y[[2L]]) ##' } ##' ##' p <- c(10, 28, 8 / 3) ##' y0 <- c(10, 1, 1) ##' ##' tt <- seq(0, 100, length.out = 40000) ##' y <- dde::dopri(y0, tt, lorenz, p, return_time = FALSE) ##' plot(y[, c(1, 3)], type = "l", lwd = 0.5, col = "#00000066") ##' ##' # If we want to print progress as the integration progresses we can ##' # use the verbose argument: ##' y <- dde::dopri(y0, c(0, 0.1), lorenz, p, verbose = TRUE) ##' ##' # Or print the y values too using a callback: ##' callback <- function(t, h, y) { ##' message(sprintf("t: %f, h: %e, y: [%s]", t, h, ##' paste(format(y, 5), collapse = ", "))) ##' } ##' y <- dde::dopri(y0, c(0, 0.1), lorenz, p, verbose = TRUE, ##' callback = callback) dopri <- function(y, times, func, parms, ..., n_out = 0L, output = NULL, rtol = 1e-6, atol = 1e-6, step_size_min = 0, step_size_max = Inf, step_size_initial = 0, step_max_n = 100000L, step_size_min_allow = FALSE, tcrit = NULL, event_time = NULL, event_function = NULL, method = "dopri5", stiff_check = 0, verbose = FALSE, callback = NULL, n_history = 0, grow_history = FALSE, return_history = n_history > 0, dllname = "", parms_are_real = TRUE, ynames = names(y), outnames = NULL, return_by_column = TRUE, return_initial = TRUE, return_time = TRUE, return_output_with_y = TRUE, return_statistics = FALSE, restartable = FALSE, return_minimal = FALSE) { ## TODO: include "deSolve" mode where we do the transpose, add the ## time column too? DOTS <- list(...) if (length(DOTS) > 0L) { stop("Invalid dot arguments!") } if (return_minimal) { return_by_column <- FALSE return_initial <- FALSE return_time <- FALSE return_output_with_y <- FALSE ynames <- NULL } func <- find_function_address(func, dllname) is_r_target <- is.function(func) if (is_r_target) { parms_are_real <- FALSE parms <- list(func = func, parms = parms, rho = parent.frame()) func <- NULL if (nzchar(dllname)) { stop("dllname must not be given when using an R function for 'func'") } } use_853 <- match_value(method, dopri_methods()) == "dopri853" assert_scalar(rtol) assert_scalar(atol) assert_scalar(step_size_min) assert_scalar(step_size_max) assert_scalar(step_size_initial) assert_size(step_max_n) assert_scalar_logical(step_size_min_allow) assert_size(n_history) assert_scalar_logical(grow_history) assert_scalar_logical(return_history) assert_scalar_logical(parms_are_real) assert_scalar_logical(return_by_column) assert_scalar_logical(return_initial) assert_scalar_logical(return_statistics) assert_scalar_logical(return_time) assert_scalar_logical(return_output_with_y) assert_scalar_logical(restartable) assert_size(stiff_check) verbose <- dopri_verbose(verbose) callback <- dopri_callback(callback) ynames <- check_ynames(y, ynames) assert_size(n_out) outnames <- check_outnames(n_out, outnames) if (n_out > 0L) { output <- find_function_address(output, dllname) ## NOTE: The same-typedness of output/func is not really ## necessary, but I think it's simplest to think about if we ## enforce it. if (is_r_target) { if (!is.function(output)) { stop("output must be an R function") } parms[[DOPRI_IDX_OUTPUT]] <- output output <- NULL } else { if (is.function(output)) { stop("output must be a compiled function (name or address)") } } } else if (!is.null(output)) { stop("If 'output' is given, then n_out must be specified") } dat <- check_events(event_time, event_function, tcrit, dllname) tcrit <- dat$tcrit
## This is needed to allow things like tcrit = 1:5 (which is int)
if (!is.null(tcrit)) {
assert_numeric(tcrit)
tcrit <- as.numeric(tcrit)
}

is_event <- dat$is_event event <- dat$event

if (!is.null(event) && is_r_target) {
parms[[DOPRI_IDX_EVENT]] <- dat$event_r } ret <- .Call(Cdopri, y, as.numeric(times), func, parms, as.integer(n_out), output, parms_are_real, ## Tolerance: rtol, atol, ## Step control: step_size_min, step_size_max, step_size_initial, as.integer(step_max_n), step_size_min_allow, ## Critical and events as.numeric(tcrit), is_event, event, ## Other: use_853, as.integer(stiff_check), verbose, callback, ## Return information: as.integer(n_history), grow_history, return_history, return_initial, return_statistics, restartable) has_output <- n_out > 0L ret <- prepare_output(ret, times, ynames, outnames, has_output, return_by_column, return_initial, return_time, return_output_with_y, "time") if (restartable) { restart_data <- list(parms = parms, parms_are_real = parms_are_real, ynames = ynames, outnames = outnames, has_output = has_output, tcrit = tcrit, return_history = return_history, return_by_column = return_by_column, return_initial = return_initial, return_statistics = return_statistics, return_time = return_time, return_output_with_y = return_output_with_y) attr(ret, "restart_data") <- restart_data } ret } ##' @export ##' @rdname dopri dopri5 <- function(y, times, func, parms, ...) { dopri(y, times, func, parms, ..., method = "dopri5") } ##' @export ##' @rdname dopri dopri853 <- function(y, times, func, parms, ...) { dopri(y, times, func, parms, ..., method = "dopri853") } ##' @export ##' @rdname dopri ##' @param obj An object to continue from; this must be the results of ##' running an integration with the option \code{restartable = ##' TRUE}. Note that continuing a problem moves the pointer along ##' in time (unless \code{copy = TRUE}, and that the incoming time ##' (\code{times[[1]]}) must equal the previous time \emph{exactly}. ##' ##' @param copy Logical, indicating if the pointer should be copied ##' before continuing. If \code{TRUE}, this is non-destructive with ##' respect to the data in the original pointer so the problem can ##' be restarted multiple times. By default this is \code{FALSE} ##' because there is a (potentially very small) cost to this ##' operation. dopri_continue <- function(obj, times, y = NULL, ..., copy = FALSE, parms = NULL, tcrit = NULL, return_history = NULL, return_by_column = NULL, return_initial = NULL, return_statistics = NULL, return_time = NULL, return_output_with_y = NULL, restartable = NULL) { DOTS <- list(...) if (length(DOTS) > 0L) { stop("Invalid dot arguments!") } ptr <- attr(obj, "ptr", exact = TRUE) dat <- attr(obj, "restart_data", exact = TRUE) if (copy) { ptr <- .Call(Cdopri_copy, ptr) } if (is.null(tcrit)) { tcrit <- dat$tcrit
}
if (is.null(parms)) {
parms <- dat$parms } ## Process any given options, falling back on the previous values return_history <- logopt(return_history, dat$return_history)
return_by_column <- logopt(return_by_column, dat$return_by_column) return_initial <- logopt(return_initial, dat$return_initial)
return_statistics <- logopt(return_statistics, dat$return_statistics) return_time <- logopt(return_time, dat$return_time)
return_output_with_y <- logopt(return_output_with_y, dat$return_output_with_y) restartable <- logopt(restartable, TRUE) ret <- .Call(Cdopri_continue, ptr, y, as.numeric(times), parms, dat$parms_are_real, tcrit,
return_history, return_initial, return_statistics, restartable)

## TODO: make this work as restartable
prepare_output(ret, times, dat$ynames, dat$outnames, dat\$has_output,
return_by_column, return_initial, return_time,
return_output_with_y, "time")
}

##' Interpolate the Dormand-Prince output after an integration.  This
##' only interpolates the core integration variables and not any
##'
##' This decouples the integration of the equations and the generation
##' of output; it is not necessary for use of the package, but may
##' come in useful where you need to do (for example) root finding on
##' the time course of a problem, or generate minimal output in some
##' cases and interrogate the solution more deeply in others.  See the
##' examples and the package vignette for a full worked example.
##'
##' @title Interpolate Dormand-Prince output
##' @param h The interpolation history.  This can be the output
##'   running \code{dopri} with \code{return_history = TRUE}, or the
##'   history attribute of this object (retrievable with
##'   \code{attr(res, "history")}).
##' @param t The times at which interpolated output is required.
##'   These times must fall within the included history (i.e., the
##'   times that the original simulation was run) or an error will be
##'   thrown.
##' @export
##' @author Rich FitzJohn
##'
##' @examples
##' # Here is the Lorenz attractor implemented as an R function
##' lorenz <- function(t, y, p) {
##'   sigma <- p[[1L]]
##'   R <- p[[2L]]
##'   b <- p[[3L]]
##'   c(sigma * (y[[2L]] - y[[1L]]),
##'     R * y[[1L]] - y[[2L]] - y[[1L]] * y[[3L]],
##'     -b * y[[3L]] + y[[1L]] * y[[2L]])
##' }
##'
##' # Standard parameters and a reasonable starting point:
##' p <- c(10, 28, 8 / 3)
##' y0 <- c(10, 1, 1)
##'
##' # Run the integration for times [0, 50] and return minimal output,
##' # but *do* record and return history.
##' y <- dopri(y0, c(0, 50), lorenz, p,
##'            n_history = 5000, return_history = TRUE,
##'            return_time = FALSE, return_initial = FALSE,
##'            return_by_column = FALSE)
##'
##' # Very little output is returned (just 3 numbers being the final
##' # state of the system), but the "history" attribute is fairly
##' # large matrix of history information.  It is not printed though
##' # as its contents should not be relied on.  What does matter is
##' # the range of supported times printed (i.e., [0, 50]) and the
##' # number of entries (~2000).
##' y
##'
##' # Generate an interpolated set of variables using this; first for
##' # 1000 steps over the full range:
##' tt <- seq(0, 50, length.out = 1000)
##' yy <- dopri_interpolate(y, tt)
##' plot(yy[, c(1, 3)], type = "l")
##'
##' # Then for 50000
##' tt <- seq(0, 50, length.out = 50000)
##' yy <- dopri_interpolate(y, tt)
##' plot(yy[, c(1, 3)], type = "l")
dopri_interpolate <- function(h, t) {
if (!inherits(h, "dopri_history")) {
h <- attr(h, "history", exact = TRUE)
}

nd <- attr(h, "n") # number of equations
if (is.null(nd)) {
stop("Corrupt history object: 'n' is missing")
}
nh <- ncol(h)      # number of history entries
it <- nrow(h) - 1L # time index
ih <- nrow(h)      # step size index
order <- (nrow(h) - 2) / nd # order of integration
if (!(order == 5 || order == 8)) {
## This one should really never be triggered but acts as a
## safeguard against real weirdness with subsetting.
stop("Corrupt history object: incorrect number of rows")
}

tr <- c(h[it, 1L], h[it, nh] + h[ih, nh])
if (min(t) < tr[[1L]] || max(t) > tr[[2L]]) {
stop(sprintf("Time falls outside of range of known history [%s, %s]",
tr[[1L]], tr[[2L]]))
}

idx <- findInterval(t, h[it, ])
theta  <- (t - h[it, idx]) / h[ih, idx]
theta1 <- 1.0 - theta

ret <- matrix(NA_real_, length(t), nd)

## Next, we need to put history into the right shape.
history <- array(h[seq_len(order * nd), ], c(nd, order, nh))
## Then pulling apart into the 5 coefficients
h1 <- history[, 1L, ]
h2 <- history[, 2L, ]
h3 <- history[, 3L, ]
h4 <- history[, 4L, ]
h5 <- history[, 5L, ]
if (order == 5) {
for (i in seq_len(nd)) {
ret[, i] <- h1[i, idx] + theta *
(h2[i, idx] + theta1 *
(h3[i, idx] + theta *
(h4[i, idx] + theta1 *
h5[i, idx])))
}
} else { # order == 8
h6 <- history[, 6L, ]
h7 <- history[, 7L, ]
h8 <- history[, 8L, ]
for (i in seq_len(nd)) {
tmp <- h5[i, idx] + theta *
(h6[i, idx] + theta1 *
(h7[i, idx] + theta *
h8[i, idx]))
ret[, i] <- h1[i, idx] + theta *
(h2[i, idx] + theta1 *
(h3[i, idx] + theta *
(h4[i, idx] + theta1 *
tmp)))
}
}

ret
}

##' @export
print.dopri_history <- function(x, ...) {
nd <- attr(x, "n") # number of equations
nh <- ncol(x)
it <- nrow(x) - 1L
ih <- nrow(x)
order <- (nrow(x) - 2) / nd
cat("<dopri_history>\n")
cat(sprintf("  - equations: %d\n", nd))
cat(sprintf("  - time: [%s, %s]\n", x[it, 1L], x[it, nh] + x[ih, nh]))
cat(sprintf("  - entries: %d\n", nh))
cat(sprintf("  - order: %d\n", order))
}

##' @export
##' @rdname dopri
##'
##' @param t The time to access (not that this is not an offset,
##'   but the actual time; within your target function you'd write
##'   things like \code{tlag(t - 1)} to get 1 time unit ago.
##'
##' @param i index within the state vector \code{y} to return.  The
##'   index here is R-style base-1 indexing, so pass \code{1} in to
##'   access the first element.  This can be left \code{NULL} to
##'   return all the elements or a vector longer than one.
ylag <- function(t, i = NULL) {
.Call(Cylag, t, i)
}

dopri_methods <- function() {
c("dopri5", "dopri853")
}

dopri_verbose <- function(verbose) {
assert_scalar(verbose)
if (is.logical(verbose)) {
verbose <- if (verbose) VERBOSE_STEP else VERBOSE_QUIET
} else {
assert_integer(verbose)
valid <- c(VERBOSE_QUIET, VERBOSE_STEP, VERBOSE_EVAL)
if (!any(verbose == valid)) {
stop("Invalid value for verbose")
}
verbose <- as.integer(verbose)
}
verbose
}

dopri_callback <- function(callback) {
if (is.null(callback)) {
return(NULL)
}
if (!is.function(callback)) {
stop("Expected a function for 'callback'")
}
if (length(formals(callback)) != 3) {
stop("Expected a function with 3 arguments for 'callback'")
}
list(callback, new.env(parent = environment(callback)))
}

DOPRI_IDX_TARGET <- 1L
DOPRI_IDX_PARMS <- 2L
DOPRI_IDX_ENV <- 3L
DOPRI_IDX_OUTPUT <- 4L
DOPRI_IDX_EVENT <- 5L

## Must match up with enum dopri_verbose in dopri.h
VERBOSE_QUIET <- 0L
VERBOSE_STEP <- 1L
VERBOSE_EVAL <- 2L


## Try the dde package in your browser

Any scripts or data that you put into this service are public.

dde documentation built on May 29, 2024, 8:15 a.m.