context("dde")
test_that("dde", {
skip_if_not_installed("deSolve")
## The pre-lag part of the integration
tt1 <- seq(0, 14, length.out = 101)
## Post lag, but the first bit
tt2 <- seq(0, 30, length.out = 301)
## The entire interesting phase
tt3 <- seq(0, 200, length.out = 301)
yy1 <- run_seir_deSolve(tt1)
yy2 <- run_seir_deSolve(tt2)
yy3 <- run_seir_deSolve(tt3)
for (method in dopri_methods()) {
## Before the lag; this one is easy:
expect_equal(run_seir_dde(tt1, method = method), yy1)
expect_equal(run_seir_dde(tt2, method = method), yy2)
## Entire interesting region
expect_equal(run_seir_dde(tt3, method = method), yy3, tolerance = 1e-6)
}
## Confirm that different integrators were actually run here:
y5 <- run_seir_dde(tt3, return_statistics = TRUE)
y8 <- run_seir_dde(tt3, method = "dopri853", return_statistics = TRUE)
expect_false(identical(y5, y8))
## The 853 stepper should take fewer steps (though in this case it's
## more ~10% more function evaluations because of more rejected
## steps).
s5 <- attr(y5, "statistics")
s8 <- attr(y8, "statistics")
expect_gt(s5[["n_step"]], s8[["n_step"]])
expect_lt(s5[["n_eval"]], s8[["n_eval"]])
## Run again with a critical time at the point the delay starts:
y5_2 <- run_seir_dde(tt3, return_statistics = TRUE, tcrit = 14)
y8_2 <- run_seir_dde(tt3, method = "dopri853", return_statistics = TRUE,
tcrit = 14)
s5_2 <- attr(y5_2, "statistics")
s8_2 <- attr(y8_2, "statistics")
expect_gt(s5_2[["n_step"]], s8_2[["n_step"]])
expect_gt(s5_2[["n_eval"]], s8_2[["n_eval"]])
expect_lt(s8_2[["n_reject"]], s8[["n_reject"]])
})
test_that("output", {
tt <- seq(0, 200, length.out = 301)
p <- numeric(0)
y0 <- c(1e7 - 1, 0, 1, 0)
res1 <- dopri(y0, tt, "seir", p, n_history = 1000L,
dllname = "dde_seir", return_history = FALSE)
expect_equal(names(attributes(res1)), "dim")
res2 <- dopri(y0, tt, "seir", p, n_history = 1000L,
n_out = 1L, output = "seir_output",
dllname = "dde_seir", return_history = FALSE,
return_output_with_y = FALSE)
output <- attr(res2, "output")
expect_equal(dim(output), c(length(tt), 1L))
expect_equal(drop(output), rowSums(res2[, -1]), tolerance = 1e-14)
attr(res2, "output") <- NULL
expect_identical(res1, res2)
## Corner case with the first output entry
res3 <- dopri(y0, tt, "seir", p, n_history = 1000L,
n_out = 1L, output = "seir_output",
dllname = "dde_seir", return_history = FALSE,
return_initial = FALSE, return_output_with_y = FALSE)
m3 <- res3[]
attr(m3, "output") <- NULL
expect_identical(m3, res1[-1, ])
out3 <- attr(res3, "output")
expect_identical(out3, output[-1, , drop = FALSE])
})
test_that("R interface", {
tt <- seq(0, 200, length.out = 301)
p <- numeric(0)
y0 <- c(1e7 - 1, 0, 1, 0)
res1 <- dopri(y0, tt, "seir", p, n_history = 1000L,
dllname = "dde_seir", return_history = FALSE)
seir <- function(t, y, p) {
b <- 1.0 / 10.0
N <- 1e7
beta <- 10.0
sigma <- 1.0 / 3.0
delta <- 1.0 / 21.0
lat_hum <- 14.0
Births <- N * b
surv <- exp(-b * lat_hum)
S <- y[[1L]]
E <- y[[2L]]
I <- y[[3L]]
R <- y[[4L]]
if (p == "nonnumeric") {
tau <- p
} else {
tau <- t - lat_hum
if (p == "as.integer" && abs(tau - as.integer(tau)) < 1e-6) {
tau <- as.integer(tau)
} else if (p == "nonscalar") {
tau <- rep(tau, 2)
}
}
if (p == "one") {
S_lag <- ylag(tau, 1L)
I_lag <- ylag(tau, 3L)
} else if (p == "idx") {
y_lag <- ylag(tau, c(1L, 3L))
S_lag <- y_lag[[1L]]
I_lag <- y_lag[[2L]]
} else {
y_lag <- ylag(tau)
S_lag <- y_lag[[1L]]
I_lag <- y_lag[[3L]]
}
new_inf <- beta * S * I / N
lag_inf <- beta * S_lag * I_lag * surv / N
c(Births - b * S - new_inf + delta * R,
new_inf - lag_inf - b * E,
lag_inf - (b + sigma) * I,
sigma * I - b * R - delta * R)
}
for (method in dopri_methods()) {
res2 <- dopri(y0, tt, seir, "one", n_history = 1000L,
return_history = FALSE, method = method)
res3 <- dopri(y0, tt, seir, "idx", n_history = 1000L,
return_history = FALSE, method = method)
res4 <- dopri(y0, tt, seir, "all", n_history = 1000L,
return_history = FALSE, method = method)
expect_equal(res2, res1,
tolerance = if (method == "dopri5") 1e-8 else 1e-5)
expect_identical(res3, res2)
expect_identical(res4, res2)
}
## Check some invalid input here too
expect_error(dopri(y0, tt, seir, "nonnumeric", n_history = 1000L,
return_history = FALSE, method = method),
"Expected a double")
expect_error(dopri(y0, tt, seir, "nonscalar", n_history = 1000L,
return_history = FALSE, method = method),
"Expected a scalar")
res <- dopri(y0, tt, seir, "as.integer", n_history = 1000L,
return_history = FALSE, method = method)
expect_equal(res, res2)
})
test_that("R interface with output", {
growth <- function(t, y, p) {
tau <- t - 2.0
c(y[[1L]], ylag(tau, 2L)) * p
}
output <- function(t, y, p) {
tau <- t - 2.0
ylag(tau, 1L)
}
tt <- seq(0, 20, length.out = 501)
p <- 0.1
y0 <- c(.1, .1)
for (method in dopri_methods()) {
if (method == "dopri5") {
res <- dopri(y0, tt, growth, p,
n_out = 1, output = output,
n_history = 1000L, return_output_with_y = FALSE,
atol = 1e-8, rtol = 1e-8,
tcrit = 2, method = method)
} else {
expect_error(dopri(y0, tt, growth, p,
n_out = 1, output = output,
n_history = 1000L, return_output_with_y = FALSE,
tcrit = 2, method = method),
"step size vanished")
## The fix here is to include a number of multiples of the delay
## time that causes the discontinutites in the solution and all
## the bits that the integration relies on; this should be dealt
## with using some sort of discontinuty pruning approach I think.
tcrit <- seq(2, 20, by = 2)
res <- dopri(y0, tt, growth, p,
n_out = 1, output = output,
n_history = 1000L, return_output_with_y = FALSE,
atol = 1e-8, rtol = 1e-8,
tcrit = tcrit, method = method)
}
tol <- 1e-7
i <- tt <= 2.0
## The first entry is easy:
expect_equal(res[, 2], y0[1] * exp(p * tt), tolerance = tol)
## The second entry is in two parts, and only the first part is easy
## (for me) to compute:
expect_equal(res[i, 3], y0[2] + y0[2] * p * tt[i], tolerance = tol)
## The output:
out <- drop(attr(res, "output"))
expect_equal(out[i], rep(y0[2], sum(i)))
expect_equal(out[!i], y0[1] * exp(p * (tt[!i] - 2.0)), tolerance = tol)
}
})
## There's a lot of really nasty bits here, so I'll build this test up
## in pieces. Eventually we want to show that running a DDE with
## negative time gives the same result as an equivalent DDE with
## positive time.
test_that("delay, negative time", {
growth0 <- function(t, y, p) {
y * p
}
output0 <- function(t, y, p) {
ylag(t - sign(p[1]) * 2)
}
set.seed(1)
tt <- seq(0, 5, length.out = 501)
r <- runif(2)
y0 <- runif(2)
ih <- seq_len(length(y0) * 5)
real_fwd <- t(y0 * exp(outer(-r, -tt)))
real_back <- t(y0 * exp(outer(r, tt)))
expect_equal(real_fwd, real_back)
## First, the simplest case; running with no output
res_fwd <- dopri(y0, tt, growth0, r, n_history = 1000L,
atol = 1e-8, rtol = 1e-8)
res_back <- dopri(y0, -tt, growth0, -r, n_history = 1000L,
atol = 1e-8, rtol = 1e-8)
expect_true(all.equal(res_fwd[, -1], real_fwd, check.attributes = FALSE))
expect_true(all.equal(res_back[, -1], real_back, check.attributes = FALSE))
expect_identical(attr(res_fwd, "history")[ih, ],
attr(res_back, "history")[ih, ])
res_fwd <- dopri(y0, tt, growth0, r, output = output0, n_out = 2L,
n_history = 1000L, return_output_with_y = FALSE,
atol = 1e-8, rtol = 1e-8)
res_back <- dopri(y0, -tt, growth0, -r, output = output0, n_out = 2L,
n_history = 1000L, return_output_with_y = FALSE,
atol = 1e-8, rtol = 1e-8)
## After this, the output agrees:
real_output <- matrix(rep(y0, length(tt)), length(y0))
i <- tt > 2
real_output[, i] <- y0 * exp(outer(r, tt[i] - 2))
real_output <- t(real_output)
out_fwd <- attr(res_fwd, "output")
out_back <- attr(res_back, "output")
expect_true(all.equal(out_fwd, real_output))
expect_equal(out_fwd, out_back)
})
test_that("failure to fetch history", {
tt <- seq(0, 30, length.out = 301)
expect_error(run_seir_dde(tt, n_history = 2L),
"Integration failure: did not find time")
})
test_that("ylag_vec_int", {
times <- seq(0, 20, length.out = 101)
p <- numeric(0)
y0 <- c(1e7 - 1, 0, 1, 0)
for (method in dopri_methods()) {
cmp <- dopri(y0, times, "seir", p, n_history = 1000L, method = method,
dllname = "dde_seir")
res <- dopri(y0, times, "seir", p, n_history = 1000L, method = method,
dllname = "dde_seir_int")
expect_equal(cmp, res, tolerance = 1e-8)
}
})
test_that("Zero lag time", {
## Oops; this is not good; delay is too short
growth <- function(t, y, p) {
c(y[[1L]], ylag(t, 2L)) * p
}
tt <- seq(0, 20, length.out = 501)
p <- 0.1
y0 <- c(.1, .1)
expect_error(dopri(y0, tt, growth, p, n_history = 1000L),
"did not find time in history")
})
## OK, most of the fundamentals are here now, but need to get the
## argument handling correct (I'm passing a NULL through where there
## should have been a TRUE/FALSE). Once this is basically working,
## I'll get this going for the difeq version too.
test_that("restart", {
tt <- seq(0, 200, length.out = 101)
tt1 <- tt[tt < 80]
tt2 <- tt[tt >= tt1[length(tt1)]]
cmp <- run_seir_dde(tt, return_minimal = TRUE)
cmp1 <- run_seir_dde(tt1, return_minimal = TRUE)
cmp2 <- cmp[, tt[-1] >= tt1[length(tt1)]]
res1 <- run_seir_dde(tt1, restartable = TRUE, return_minimal = TRUE)
expect_is(attr(res1, "ptr"), "externalptr")
expect_equal(res1, cmp1, check.attributes = FALSE)
## Then start work on the continuation
res2 <- dopri_continue(res1, tt2, return_initial = TRUE)
## NOTE: testthat will report incorrectly here on failure because of
## attributes (at least of 1.0.2)
expect_equal(res2, cmp2, check.attributes = FALSE)
expect_error(dopri_continue(res1, tt2, return_initial = TRUE),
"Incorrect initial time on integration restart")
})
test_that("restart and copy", {
tt <- seq(0, 200, length.out = 101)
tt1 <- tt[tt < 80]
tt2 <- tt[tt >= tt1[length(tt1)]]
cmp <- run_seir_dde(tt, return_minimal = TRUE)
cmp1 <- run_seir_dde(tt1, return_minimal = TRUE)
cmp2 <- cmp[, tt[-1] >= tt1[length(tt1)]]
res1 <- run_seir_dde(tt1, restartable = TRUE, return_minimal = TRUE)
expect_is(attr(res1, "ptr"), "externalptr")
res2 <- dopri_continue(res1, tt2, return_initial = TRUE, copy = TRUE)
res3 <- dopri_continue(res1, tt2, return_initial = TRUE, copy = TRUE)
expect_equal(res2, cmp2, check.attributes = FALSE)
expect_equal(res2, res3)
## Continue off and things should work OK here
res4 <- dopri_continue(res1, tt2, return_initial = TRUE, copy = FALSE)
expect_equal(res4, cmp2, check.attributes = FALSE)
## But we've modified the pointer so this will no longer work:
expect_error(dopri_continue(res1, tt2, return_initial = TRUE),
"Incorrect initial time on integration restart")
})
test_that("restart errors", {
tt <- seq(0, 200, length.out = 101)
tt1 <- tt[tt < 80]
tt2 <- tt[tt >= tt1[length(tt1)]]
cmp <- run_seir_dde(tt, return_minimal = TRUE)
cmp1 <- run_seir_dde(tt1, return_minimal = TRUE)
cmp2 <- cmp[, tt[-1] >= tt1[length(tt1)]]
res1 <- run_seir_dde(tt1, restartable = TRUE, return_minimal = TRUE)
expect_error(dopri_continue(res1, tt2[[1]], copy = TRUE),
"At least two times must be given")
expect_error(dopri_continue(res1, rev(tt1), copy = TRUE),
"Incorrect sign for the times")
y1 <- res1[, ncol(res1)]
expect_error(dopri_continue(res1, tt2, y1[-1], copy = TRUE),
"Incorrect size 'y' on integration restart")
expect_error(dopri_continue(res1, tt2, rep(y1, 2), copy = TRUE),
"Incorrect size 'y' on integration restart")
ptr <- attr(res1, "ptr")
## Make a NULL pointer:
attr(res1, "ptr") <- make_null_pointer(ptr)
expect_error(dopri_continue(res1, tt2, copy = TRUE),
"pointer has been freed")
## Delete the pointer:
attr(res1, "ptr") <- NULL
expect_error(dopri_continue(res1, tt2, copy = TRUE),
"Expected an external pointer")
## Something stupid as a pointer:
attr(res1, "ptr") <- res1
expect_error(dopri_continue(res1, tt2, copy = TRUE),
"Expected an external pointer")
})
test_that("change y on restart", {
growth <- function(t, y, p) {
y * p
}
output <- function(t, y, p) {
ylag(t - 2.0)
}
tt <- seq(0, 10, length.out = 101)
tc <- 4
i1 <- tt <= tc
i2 <- tt >= tc
tt1 <- tt[i1]
tt2 <- tt[i2]
y0 <- 1
r <- -0.5
j <- seq_along(y0)
## TODO: error is not very clear when output = NULL but n_out > 0
##
## TODO: error is not very clear when n_history = 0
##
## TODO: error if tcrit[1] == t[1]
res <- dopri(y0, tt, growth, r, n_out = length(y0), output = output,
n_history = 1000, return_history = FALSE,
return_time = FALSE, return_by_column = FALSE,
tcrit = tt2[[1]], atol = 1e-8, rtol = 1e-8)
res1 <- dopri(y0, tt1, growth, r, n_out = length(y0), output = output,
n_history = 1000, return_history = FALSE,
return_time = FALSE, return_by_column = FALSE,
restartable = TRUE)
y1 <- res1[j, ncol(res1)]
res2 <- dopri_continue(res1, tt2, y1, copy = TRUE)
## Check that options have been preserved:
expect_equal(nrow(res2), nrow(res1))
expect_equal(ncol(res2), length(tt2))
expect_null(attr(res2, "output"))
## TODO: not sure why this is very slightly off, but it's possibly
## due to a change in step size? I don't think that this should be
## the case though and would like to see this shrink.
expect_equal(res2[j, ], res[j, i2], tolerance = 5e-7)
expect_equal(res2[-j, ], res[-j, i2], tolerance = 5e-7)
## Change y on re-entry:
y2 <- y1 * 2
res3 <- dopri_continue(res1, tt2, y2, copy = TRUE, tcrit = tt2[[1]] + 2)
## TODO: this is not terrific accuracy either:
expect_equal(res3[j, ] / res2[j, ], rep(2, length(tt2)), tolerance = 1e-5)
k <- tt2 - tt2[[1]] < 2.0
expect_equal(res3[-j, k] / res2[-j, k], rep(1, sum(k)))
expect_equal(res3[-j, !k] / res2[-j, !k], rep(2, sum(!k)), tolerance = 1e-5)
})
## This indicates to me that we might still have a problem with the
## solver, but it might also be that the high order difficulty from
## the echo of the delay is causing trouble. If you push out tcrit
## values as multiples of 14 and drop tolerance to 1e-5 you can get as
## far as t = 56 before it gives up.
test_that("stiffness detection", {
tt <- seq(0, 200, length.out = 301)
yy <- run_seir_dde(tt)
expect_equal(run_seir_dde(tt, stiff_check = 1), yy)
expect_error(run_seir_dde(tt, method = "dopri853", stiff_check = 1),
"problem became stiff")
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.