context("ode")
test_that("ode interface", {
skip_if_not_installed("deSolve")
tt <- seq(0, 1, length.out = 200)
m1 <- run_lorenz_deSolve(tt)
dimnames(m1) <- NULL
m2 <- run_lorenz_dde(tt)
expect_equal(m2, m1, tolerance = 1e-6)
m3 <- run_lorenz_dde(tt, return_minimal = TRUE)
expect_equal(t(m3), m1[-1, -1, drop = FALSE], tolerance = 1e-6)
expect_identical(m3, t(m2[-1, -1]))
})
test_that("integer time input", {
tt <- seq(0, 10, by = 1)
expect_identical(run_lorenz_dde(as.integer(tt)),
run_lorenz_dde(tt))
})
test_that("zero time difference", {
tt <- c(0, 0, 1, 1, 2, 2, 2, 3, 3, 3)
res1 <- run_lorenz_dde(tt, return_initial = TRUE)
res2 <- run_lorenz_dde(unique(tt), return_initial = TRUE)
expect_identical(res1, res2[tapply(tt, tt), ])
})
test_that("ode, 873 stepper", {
tt <- seq(0, 1, length.out = 200)
m5 <- run_lorenz_dde(tt)
m8 <- run_lorenz_dde(tt, method = "dopri853")
expect_equal(m5, m8, tolerance = 1e-7)
expect_false(identical(m5, m8))
})
test_that("dense output", {
skip_if_not_installed("deSolve")
tt <- seq(0, 1, length.out = 200)
m1 <- run_lorenz_deSolve(tt)
dimnames(m1) <- NULL
for (method in dopri_methods()) {
m2 <- run_lorenz_dde(c(0, max(tt)), n_history = 1000L, method = method)
## single row of output:
expect_equal(dim(m2), c(2, 4))
h2 <- attr(m2, "history")
expect_is(h2, "dopri_history")
expect_equal(nrow(h2),
if (method == "dopri5") 17L else 26L) # c(5, 8) * 3 + 2
expect_identical(attr(h2, "n"), 3L)
expect_output(print(h2), "<dopri_history>", fixed = TRUE)
m3 <- dopri_interpolate(h2, tt)
m4 <- run_lorenz_dde(tt, method = method)
expect_equal(m3, m4[, -1], tolerance = 1e-10)
expect_equal(m3, m1[, -1], tolerance = 1e-6)
expect_equal(dopri_interpolate(m2, tt), m3)
## Check row output:
m5 <- run_lorenz_dde(tt, n_history = 1000L, return_by_column = FALSE,
method = method)
expect_identical(attr(m5, "history"), attr(m2, "history"))
attr(m5, "history") <- NULL
expect_identical(m5, t(m4))
expect_error(dopri_interpolate(h2, tt - 1),
"Time falls outside of range of known history")
expect_error(dopri_interpolate(h2, tt + 1),
"Time falls outside of range of known history")
expect_error(dopri_interpolate(h2[-1, ], tt),
"Corrupt history object")
h2 <- h2[-1, ]
attr(h2, "n") <- 3L
class(h2) <- "dopri_history"
expect_error(dopri_interpolate(h2, tt),
"Corrupt history object: incorrect number of rows")
}
})
test_that("output", {
tt <- seq(0, 1, length.out = 200)
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
res1 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz")
res2 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
return_output_with_y = FALSE,
n_out = 2L, output = "lorenz_output")
expect_equal(names(attributes(res1)), "dim")
output <- attr(res2, "output")
expect_equal(dim(output), c(nrow(res1), 2L))
expect_equal(output[, 1L], pmin(res1[, 2], res1[, 3], res1[, 4]))
expect_equal(output[, 2L], pmax(res1[, 2], res1[, 3], res1[, 4]))
attr(res2, "output") <- NULL
expect_identical(res1, res2)
res3 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output",
return_output_with_y = FALSE, return_by_column = FALSE)
expect_equal(attr(res3, "output"), t(output))
attr(res3, "output") <- NULL
expect_identical(res3, t(res2))
res4 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output",
return_by_column = FALSE, return_output_with_y = TRUE)
expect_null(attr(res4, "output"))
expect_equal(dim(res4), c(6, length(tt)))
expect_equal(res4[1:4, ], res3, check.attributes = FALSE)
expect_equal(res4[5:6, ], t(output), check.attributes = FALSE)
res5 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output",
return_by_column = FALSE, return_output_with_y = TRUE,
return_time = FALSE, return_initial = FALSE)
expect_null(attr(res5, "output"))
expect_equal(dim(res5), c(5, length(tt) - 1L))
expect_equal(res5[1:3, ], res3[-1, -1], check.attributes = FALSE)
expect_equal(res5[4:5, ], t(output[-1, ]), check.attributes = FALSE)
})
test_that("keep initial", {
tt <- seq(0, 1, length.out = 200)
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
res1 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
return_initial = TRUE)
res2 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
return_initial = FALSE)
expect_equal(nrow(res1), length(tt))
expect_identical(res1[1, -1], y0)
expect_identical(res1[-1, ], res2)
res3 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output", return_initial = TRUE,
return_output_with_y = FALSE)
res4 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output", return_initial = FALSE,
return_output_with_y = FALSE)
expect_equal(nrow(res3), length(tt))
expect_identical(res3[1, -1], y0)
out3 <- attr(res3, "output")
out4 <- attr(res4, "output")
expect_equal(nrow(out3), length(tt))
expect_equal(out3[1, ], range(y0)) # based on definition of output function
expect_equal(out3[-1, ], out4)
output <- attr(res3, "output")
attr(res3, "output") <- attr(res4, "output") <- NULL
expect_equal(nrow(res3), length(tt))
expect_identical(res3[1, -1], y0)
expect_identical(res3[-1, ], res4)
res5 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output", return_initial = TRUE,
return_output_with_y = TRUE)
res6 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output", return_initial = FALSE,
return_output_with_y = TRUE)
expect_null(attr(res5, "output"))
expect_null(attr(res6, "output"))
expect_equal(res5[, 1:4], res3)
expect_equal(res6[, 1:4], res4)
expect_equal(res5[, 5:6], output)
expect_equal(res6[, 5:6], output[-1, ])
})
test_that("R interface", {
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
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]])
}
lorenz_output <- function(t, y, p) {
c(min(y), max(y))
}
tt <- seq(0, 1, length.out = 200)
res1 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz")
res2 <- dopri(y0, tt, lorenz, p)
expect_equal(res1, res2, tolerance = 1e-10)
res3 <- dopri(y0, tt, "lorenz", p, dllname = "dde_lorenz",
n_out = 2L, output = "lorenz_output")
res4 <- dopri(y0, tt, lorenz, p, n_out = 2L, output = lorenz_output)
expect_equal(res3, res4, tolerance = 1e-10)
})
test_that("critical times", {
target <- function(t, y, p) {
if (t <= 1) {
y
} else {
-5 * y
}
}
tt <- seq(0, 2, length.out = 200)
res1 <- dopri(1, tt, target, numeric(0), return_statistics = TRUE)
res2 <- dopri(1, tt, target, numeric(0), tcrit = 1, return_statistics = TRUE)
res3 <- dopri(1, tt, target, numeric(0), tcrit = c(-1, 1),
return_statistics = TRUE)
s1 <- attr(res1, "statistics")
s2 <- attr(res2, "statistics")
s3 <- attr(res3, "statistics")
expect_lt(s2[["n_step"]], s1[["n_step"]])
expect_equal(s2, s3)
expect_identical(res2, res3)
expect_is(attr(res1, "step_size"), "numeric")
## I also need to check a pathology here:
res4 <- dopri(1, tt, target, numeric(0), tcrit = c(tt[[1]], 1),
return_statistics = TRUE)
expect_equal(res4, res3)
## Bunch of pathalogical times:
res5 <- dopri(1, tt, target, numeric(0), tcrit = rep(tt[[1]], 3),
return_statistics = TRUE)
expect_equal(res5, res1)
## Pile up some times in the middle:
res6 <- dopri(1, tt, target, numeric(0), tcrit = rep(1, 3),
return_statistics = TRUE)
expect_equal(res4, res3)
})
test_that("names", {
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
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]])
}
lorenz_output <- function(t, y, p) {
c(min(y), max(y))
}
tt <- seq(0, 1, length.out = 200)
nms <- letters[1:3]
cmp <- list(NULL, c("time", nms))
expect_null(dimnames(dopri(y0, tt, lorenz, p)))
expect_equal(dimnames(dopri(y0, tt, lorenz, p, ynames = nms)), cmp)
expect_equal(dimnames(dopri(setNames(y0, nms), tt, lorenz, p)), cmp)
expect_null(dimnames(dopri(setNames(y0, nms), tt, lorenz, p, ynames = FALSE)))
expect_error(dopri(y0, tt, lorenz, p, ynames = nms[1]),
"ynames must be the same length as y")
expect_error(dopri(y0, tt, lorenz, p, ynames = 1),
"Invalid value for ynames")
## Similar for output names:
onms <- LETTERS[1:2]
ocmp <- list(NULL, onms)
f <- function(..., return_output_with_y = FALSE) {
dopri(y0, tt, lorenz, p, n_out = 2L, output = lorenz_output,
return_output_with_y = return_output_with_y, ...)
}
expect_null(dimnames(attr(f(), "output")))
expect_null(dimnames(attr(f(outnames = NULL), "output")))
expect_equal(dimnames(attr(f(outnames = onms), "output")), ocmp)
expect_error(f(outnames = nms), "outnames must have length n_out")
expect_error(f(outnames = 1), "Invalid value for outnames")
## Check both together:
res <- f(ynames = nms, outnames = onms)
expect_equal(dimnames(res), cmp)
expect_equal(dimnames(attr(res, "output")), ocmp)
## And check when output is combined
m <- f(ynames = nms, outnames = NULL, return_output_with_y = TRUE)
expect_equal(colnames(m), c("time", nms, rep("", length(onms))))
m <- f(ynames = NULL, outnames = onms, return_output_with_y = TRUE)
expect_equal(colnames(m), c("time", rep("", length(nms)), onms))
m <- f(ynames = nms, outnames = onms, return_output_with_y = TRUE)
expect_equal(colnames(m), c("time", nms, onms))
## And check time
expect_null(dimnames(f(return_time = TRUE)))
expect_equal(colnames(f(return_time = TRUE, ynames = nms)),
c("time", nms))
expect_equal(colnames(f(return_time = TRUE, return_output_with_y = TRUE,
ynames = nms)),
c("time", nms, rep("", length(onms))))
expect_equal(colnames(f(return_time = TRUE, return_output_with_y = TRUE,
ynames = nms, outnames = onms)),
c("time", nms, onms))
expect_equal(colnames(f(return_time = TRUE, return_output_with_y = TRUE,
ynames = NULL, outnames = onms)),
c("time", rep("", length(nms)), onms))
expect_null(colnames(f(return_time = TRUE, return_output_with_y = TRUE,
ynames = NULL, outnames = NULL)))
})
test_that("return_time", {
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
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]])
}
tt <- seq(0, 1, length.out = 200)
res1 <- dopri(y0, tt, lorenz, p)
res2 <- dopri(y0, tt, lorenz, p, return_time = FALSE)
expect_equal(ncol(res1), 4)
expect_equal(ncol(res2), 3)
expect_equal(res1[, 1L], tt)
expect_equal(res1[, -1L], res2)
res3 <- dopri(y0, tt, lorenz, p, return_initial = FALSE)
## include the first time:
expect_equal(res3[, 1], tt[-1])
## with names:
nms <- letters[1:3]
res <- dopri(y0, tt, lorenz, p, ynames = nms)
expect_equal(colnames(res), c("time", nms))
res <- dopri(y0, tt, lorenz, p, ynames = nms, return_by_column = FALSE)
expect_equal(rownames(res), c("time", nms))
})
test_that("minimal mode", {
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
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]])
}
tt <- seq(0, 1, length.out = 200)
names(y0) <- letters[1:3]
res1 <- dopri(y0, tt, lorenz, p)
res2 <- dopri(y0, tt, lorenz, p, return_minimal = TRUE)
expect_equal(dim(res1), c(length(tt), 4))
expect_equal(dim(res2), c(3, length(tt) - 1L))
expect_equal(dimnames(res1), list(NULL, c("time", names(y0))))
expect_null(dimnames(res2))
expect_equal(res1[, 1], tt)
expect_equal(t(unname(res1[-1, -1])), res2)
})
test_that("step tuning", {
tt <- seq(0, 1, length.out = 200)
T_IDX <- 16L
m0 <- run_lorenz_dde(tt, n_history = 500, return_statistics = TRUE)
t0 <- attr(m0, "history")[T_IDX, ]
s0 <- attr(m0, "statistics")
hmax <- max(diff(t0)) / 5
m1 <- run_lorenz_dde(tt, step_size_max = hmax, n_history = 500)
t1 <- attr(m1, "history")[T_IDX, ]
expect_gt(length(t1), length(t0))
expect_equal(max(diff(t1)), hmax)
h0 <- diff(t0[1:2])
m2 <- run_lorenz_dde(tt, n_history = 500, return_statistics = TRUE,
step_size_initial = h0)
t2 <- attr(m0, "history")[T_IDX, ]
expect_equal(t0, t2)
expect_equal(attr(m2, "statistics")[["n_eval"]],
s0[["n_eval"]] - 1L)
expect_error(run_lorenz_dde(tt, step_max_n = length(t0) - 20),
"too many steps")
expect_error(run_lorenz_dde(tt, step_size_min = min(diff(t0)) * 1.1),
"step size too small")
})
test_that("allow step size to get small", {
tt <- seq(0, 1, length.out = 200)
T_IDX <- 16L
m0 <- run_lorenz_dde(tt, n_history = 500, return_statistics = TRUE)
t0 <- attr(m0, "history")[T_IDX, ]
s0 <- attr(m0, "statistics")
step_size_min <- 0.01
m1 <- run_lorenz_dde(tt,
step_size_min = step_size_min,
step_size_min_allow = TRUE,
n_history = 500,
return_statistics = TRUE)
t1 <- attr(m1, "history")[T_IDX, ]
## This is really annoying, as we might be off by a _very_ small
## amount (~1e-18) due to general floating point horrors.
expect_lt(min(diff(t1) - step_size_min), 1e-16)
expect_equal(min(diff(t1)), step_size_min)
s1 <- attr(m1, "statistics")
expect_lt(s1[["n_eval"]], s0[["n_eval"]])
})
test_that("integrate function with no absolute error", {
deriv <- function(t, y, p) {
1
}
expect_equal(drop(dopri(0, c(0, 1), deriv, 0, return_minimal = TRUE)), 1)
})
test_that("Native Symbol interface", {
p <- c(10, 28, 8 / 3)
y <- c(10, 1, 1)
times <- seq(0, 10, length.out = 101)
func <- getNativeSymbolInfo("lorenz", PACKAGE = "dde_lorenz")
y <- dopri(y, times, func, p)
expect_equal(y, run_lorenz_dde(times, tol = 1e-6))
})
test_that("NULL pointer safety", {
p <- c(10, 28, 8 / 3)
y <- c(10, 1, 1)
times <- seq(0, 10, length.out = 101)
func <- getNativeSymbolInfo("lorenz", PACKAGE = "dde_lorenz")
func <- make_null_pointer(func)
expect_error(dopri(y, times, func, p), "null pointer")
expect_error(dopri(y, times, "lorenz", p, output = func, n_out = 2L),
"null pointer")
})
test_that("always return history when asked", {
has_history <- function(x) {
"history" %in% names(attributes(x))
}
has_statistics <- function(x) {
"statistics" %in% names(attributes(x))
}
growth <- function(t, y, p) {
y * p
}
y0 <- runif(1)
r <- runif(1)
tt <- seq(0, 2, length.out = 101)
expect_false(has_history(dopri(y0, tt, growth, r)))
expect_true(has_history(dopri(y0, tt, growth, r, n_history = 100)))
expect_true(has_history(
dopri(y0, tt, growth, r, n_history = 100, return_minimal = TRUE)))
expect_true(has_history(
dopri(y0, tt, growth, r, n_history = 100, return_by_column = TRUE)))
expect_true(has_history(
dopri(y0, tt, growth, r, n_history = 100, return_initial = TRUE)))
expect_true(has_history(
dopri(y0, tt, growth, r, n_history = 100, return_time = TRUE)))
expect_false(has_statistics(dopri(y0, tt, growth, r)))
expect_true(has_statistics(dopri(y0, tt, growth, r,
return_statistics = TRUE)))
expect_true(has_statistics(
dopri(y0, tt, growth, r, return_statistics = TRUE,
return_minimal = TRUE)))
expect_true(has_statistics(
dopri(y0, tt, growth, r, return_statistics = TRUE,
return_by_column = TRUE)))
expect_true(has_statistics(
dopri(y0, tt, growth, r, return_statistics = TRUE,
return_initial = TRUE)))
expect_true(has_statistics(
dopri(y0, tt, growth, r, return_statistics = TRUE,
return_time = TRUE)))
})
test_that("negative time", {
skip_if_not_installed("deSolve")
growth <- function(t, y, p) {
y * p
}
set.seed(1)
y0 <- runif(4)
r <- runif(4)
## True solution is:
##
## exp(r * t) * y0
##
## which we can compute easily enough foe this problem:
tt <- seq(0, -2, length.out = 11)
real <- t(y0 * outer(r, tt, function(t, r) exp(r * t)))
real_fwd <- t(y0 * outer(-r, -tt, function(t, r) exp(r * t)))
expect_identical(real, real_fwd)
## Then solve the system using deSolve (sanity check for the above)
cmp <- deSolve::ode(y0, tt, function(...) list(growth(...)), r)
cmp_fwd <- deSolve::ode(y0, -tt, function(...) list(growth(...)), -r)
expect_equal(unname(cmp[, -1]), real, 5e-6)
expect_equal(cmp[, -1], cmp_fwd[, -1], tolerance = 1e-16)
## and dde:
res <- dopri(y0, tt, growth, r, n_history = 100)
res_fwd <- dopri(y0, -tt, growth, -r, n_history = 100)
expect_equal(res[, -1], res_fwd[, -1], tolerance = 1e-16)
## Interestingly, we do stop at the same points here, so the error
## calculation and underlying stepping is probably OK.
t_idx <- length(y0) * 5 + 1
h1 <- attr(res, "history")
h2 <- attr(res_fwd, "history")
expect_equal(-h1[t_idx, ], h2[t_idx, ], tolerance = 1e-16)
expect_equal(-h1[t_idx + 1, ], h2[t_idx + 1, ], tolerance = 1e-16)
expect_equal(h1[1:20, ], h2[1:20, ], tolerance = 1e-16)
})
test_that("negative time with tcrit", {
skip_if_not_installed("deSolve")
target1 <- function(t, y, p) {
if (t <= 1) y else -5 * y
}
target2 <- function(t, y, p) {
if (t >= -1) -y else 5 * y
}
## Here's our system. It's very silly again.
tt <- seq(0, 2, length.out = 200)
y0 <- 1
y1 <- deSolve::lsoda(y0, tt, function(...) list(target1(...)), numeric())
y2 <- deSolve::lsoda(y0, -tt, function(...) list(target2(...)), numeric())
expect_equal(y1[, 2], y2[, 2], tolerance = 1e-16)
res1a <- dopri(y0, tt, target1, numeric(0), return_statistics = TRUE)
res2a <- dopri(y0, -tt, target2, numeric(0), return_statistics = TRUE)
## Not terribly accurate because of the nasty discontinuity
expect_equal(res1a[, 2], y1[, 2], tolerance = 1e-4)
## As above, the dopri solutions should agree well
expect_equal(res1a[, 2], res2a[, 2], tolerance = 1e-16)
## Add the critical time in:
res1b <- dopri(y0, tt, target1, numeric(0), tcrit = 1,
return_statistics = TRUE)
res2b <- dopri(y0, -tt, target2, numeric(0), tcrit = -1,
return_statistics = TRUE)
expect_equal(res1b[, 2], res2b[, 2], tolerance = 1e-16)
s1a <- attr(res1a, "statistics")
s2a <- attr(res2a, "statistics")
s1b <- attr(res1b, "statistics")
s2b <- attr(res2b, "statistics")
expect_equal(s1a, s2a)
expect_equal(s1b, s2b)
expect_lt(s2b[["n_step"]], s2a[["n_step"]])
})
test_that("time validation", {
expect_error(run_lorenz_dde(numeric(0)), "At least two times must be given")
expect_error(run_lorenz_dde(0), "At least two times must be given")
expect_error(run_lorenz_dde(c(0, 0)),
"Beginning and end times are the same")
expect_error(run_lorenz_dde(c(0, 2, 1)),
"Times have inconsistent sign")
})
test_that("non-real input", {
tt <- seq(0, 10, length.out = 101)
cmp <- run_lorenz_dde(tt)
p <- c(10, 28, 8 / 3)
y <- c(10, 1, 1)
res <- dopri(y, tt, "lorenz", p, atol = 1e-7, rtol = 1e-7,
parms_are_real = FALSE,
dllname = "dde_lorenz2")
expect_identical(res, cmp)
})
test_that("externalptr input", {
tt <- seq(0, 10, length.out = 101)
cmp <- run_lorenz_dde(tt)
ptr <- .Call("lorenz_init", c(10, 28, 8 / 3), PACKAGE = "dde_lorenz3")
expect_is(ptr, "externalptr")
y <- c(10, 1, 1)
res <- dopri(y, tt, "lorenz", ptr, atol = 1e-7, rtol = 1e-7,
parms_are_real = FALSE,
dllname = "dde_lorenz3")
expect_identical(res, cmp)
})
test_that("grow history", {
tt <- seq(0, 10, length.out = 101)
res <- run_lorenz_dde(tt, n_history = 5, grow_history = TRUE,
return_history = TRUE, return_statistics = TRUE)
s <- attr(res, "statistics")
h <- attr(res, "history")
expect_equal(s[["n_accept"]], ncol(h))
cmp <- run_lorenz_dde(tt, n_history = s[["n_step"]] * 2,
return_history = TRUE, return_statistics = TRUE)
expect_equal(h, attr(cmp, "history"))
expect_equal(s, attr(cmp, "statistics"))
})
test_that("initial derivative validation", {
deriv <- function(t, y, p) {
dy <- rep(1, length(y))
if (length(p) > 0L) {
dy[p] <- NA_real_
}
dy
}
tt <- c(0, 1)
y0 <- rep(1, 10)
expect_error(dopri(y0, tt, deriv, 1),
"non-finite derivative at initial time for element 1")
expect_error(dopri(y0, tt, deriv, 4:9),
"non-finite derivative at initial time for element 4")
})
test_that("verbose mode prints trace", {
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
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]])
}
tt <- seq(0, 1, length.out = 200)
names(y0) <- letters[1:3]
expect_silent(ans1 <- dopri(y0, tt, lorenz, p))
expect_output(ans2 <- dopri(y0, tt, lorenz, p, verbose = TRUE),
".step. t: .+, h: .+\n")
expect_output(ans3 <- dopri(y0, tt, lorenz, p, verbose = VERBOSE_EVAL),
".eval. t: .+, h: .+\n")
expect_identical(ans1, ans2)
expect_identical(ans1, ans3)
})
test_that("verbose with callback", {
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
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]])
}
tt <- seq(0, 1, length.out = 200)
names(y0) <- letters[1:3]
steps <- list()
evals <- list()
callback <- function(t, h, y) {
if (is.na(h)) {
evals <<- c(evals, list(list(t = t, h = h, y = y)))
} else {
steps <<- c(steps, list(list(t = t, h = h, y = y)))
}
}
## Callback successfully returns all steps - we could have done
## something more interesting here, and typically this would involve
## printing not collection.
ans1 <- dopri(y0, tt, lorenz, p, verbose = TRUE, callback = callback)
n <- length(steps)
expect_gt(n, 0)
expect_equal(steps[[1]]$t, 0)
expect_equal(steps[[1]]$y, unname(y0))
expect_equal(steps[[n]]$t + steps[[n]]$h, 1)
expect_equal(length(evals), 0)
ans2 <- dopri(y0, tt, lorenz, p, verbose = VERBOSE_EVAL, callback = callback)
expect_identical(steps[1:n], steps[-(1:n)])
m <- length(evals)
expect_gt(m, 0)
expect_equal(evals[[1]]$t, 0)
expect_equal(evals[[1]]$y, unname(y0))
expect_true(all(is.na(vapply(evals, function(x) x$h, numeric(1)))))
## Time must have changed
expect_gt(evals[[2]]$t, evals[[1]]$t)
## Output must have changed
expect_true(all(evals[[2]]$y != evals[[1]]$y))
})
test_that("check verbose argument", {
expect_identical(dopri_verbose(FALSE), VERBOSE_QUIET)
expect_identical(dopri_verbose(0.0), VERBOSE_QUIET)
expect_identical(dopri_verbose(0L), VERBOSE_QUIET)
expect_identical(dopri_verbose(TRUE), VERBOSE_STEP)
expect_identical(dopri_verbose(1.0), VERBOSE_STEP)
expect_identical(dopri_verbose(1L), VERBOSE_STEP)
expect_identical(dopri_verbose(2.0), VERBOSE_EVAL)
expect_identical(dopri_verbose(2L), VERBOSE_EVAL)
expect_error(dopri_verbose(-1), "Invalid value for verbose")
expect_error(dopri_verbose(1.5), "verbose must be integer")
expect_error(dopri_verbose(c(1, 2)), "verbose must be a scalar")
expect_error(dopri_verbose(integer(0)), "verbose must be a scalar")
expect_error(dopri_verbose(NULL), "verbose must be a scalar")
})
test_that("check callback argument", {
expect_null(dopri_callback(NULL))
f <- function(a, b, c) NULL
res <- dopri_callback(f)
expect_identical(f, res[[1]])
expect_is(res[[2]], "environment")
expect_equal(parent.env(res[[2]]), environment(f))
expect_error(dopri_callback(TRUE),
"Expected a function for 'callback'")
expect_error(dopri_callback(function() 1),
"Expected a function with 3 arguments for 'callback'")
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.