Merge pull request #141 from mrc-ide/mrc-5834

Allow dt > 1

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: dust2
 Title: Next Generation dust
-Version: 0.3.13
+Version: 0.3.14
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "rich.fitzjohn@gmail.com"),
              person("Imperial College of Science, Technology and Medicine",

R/filter-support.R

@@ -53,12 +53,18 @@ check_dt <- function(dt, allow_null = FALSE, name = deparse(substitute(dt)),
   if (dt <= 0) {
     cli::cli_abort("Expected '{name}' to be greater than 0")
   }
-  if (dt > 1) {
-    cli::cli_abort("Expected '{name}' to be at most 1")
-  }
-  if (!rlang::is_integerish(1 / dt)) {
-    cli::cli_abort("Expected '{name}' to be the inverse of an integer",
-                   arg = "dt", call = call)
+  if (dt >= 1) {
+    if (!rlang::is_integerish(dt) || is.infinite(dt)) {
+      cli::cli_abort(
+        "Expected '{name}' to be an integer, if greater than  1",
+        arg = "dt", call = call)
+    }
+  } else {
+    if (!rlang::is_integerish(1 / dt)) {
+      cli::cli_abort(
+        "Expected '{name}' to be the inverse of an integer, if less than 1",
+        arg = "dt", call = call)
+    }
   }
   dt
 }

R/interface.R

@@ -700,7 +700,7 @@ check_time_sequence <- function(time, time_control,
   }
 
   dt <- time_control$dt
-  if (!is.null(dt)) {
+  if ((!is.null(dt)) && (dt <= 1)) {
     #rem <- time %% dt # The problem is 10 %% 0.25 = 0, but 10 %% 0.1 = 0.1
     rem <- fmod(time, dt)
     err <- abs(rem) > sqrt(.Machine$double.eps)

inst/include/dust2/r/helpers.hpp

@@ -236,11 +236,15 @@ inline double check_dt(cpp11::list r_time_control, bool enabled = true,
     cpp11::stop("Expected 'dt' to be greater than 0");
   }
   if (dt > 1) {
-    cpp11::stop("Expected 'dt' to be at most 1");
+    if (!is_integer_like(dt, eps)) {
+      cpp11::stop("Expected 'dt' to be an integer, if above 1");
+    }
   }
-  const auto inv_dt = 1 / dt;
-  if (!is_integer_like(inv_dt, eps)) {
-    cpp11::stop("Expected 'dt' to be the inverse of an integer");
+  if (dt < 1) {
+    const auto inv_dt = 1 / dt;
+    if (!is_integer_like(inv_dt, eps)) {
+      cpp11::stop("Expected 'dt' to be the inverse of an integer, if below 1");
+    }
   }
   return dt;
 }

tests/testthat/test-filter-support.R

@@ -1,15 +1,18 @@
 test_that("can validate that 'dt' is reasonable", {
   expect_no_error(check_dt(1))
   expect_no_error(check_dt(1 / 5))
+  expect_no_error(check_dt(2))
   expect_no_error(check_dt(NULL, allow_null = TRUE))
   expect_error(check_dt(-1, name = "dt"),
                "Expected 'dt' to be greater than 0")
-  expect_error(check_dt(2, name = "dt"),
-               "Expected 'dt' to be at most 1")
+  expect_error(check_dt(2.5, name = "dt"),
+               "Expected 'dt' to be an integer, if greater than 1")
   expect_error(check_dt(1 / 3.5, name = "dt"),
                "Expected 'dt' to be the inverse of an integer")
   expect_error(check_dt(NULL, name = "dt"),
                "'dt' must be a scalar")
+  expect_error(check_dt(Inf, name = "dt"),
+               "Expected 'dt' to be an integer, if greater than 1")
 })
 
 
@@ -23,3 +26,21 @@ test_that("check index", {
   expect_error(check_index(idx),
                "All elements of 'idx' must be at least 1")
 })
+
+test_that("can apply stochastic updates by setting dt", {
+  sys <- dust_system_create(malaria(), list(), n_particles = 10, dt = 10)
+  dust_system_set_state_initial(sys)
+  t <- seq(0, 20)
+  y <- dust_system_simulate(sys, t)
+  i_beta <- dust_unpack_index(sys)$beta
+  beta <- y[i_beta, , ]
+
+  ## Equivalence to a deterinistic model:
+  sys0 <- dust_system_create(malaria(), list(), n_particles = 10)
+  for (i in seq_len(length(t) - 1)) {
+    dust_system_set_state(sys0, y[, , i])
+    dust_system_run_to_time(sys0, t[[i + 1]])
+    y0 <- dust_system_state(sys0)
+    #expect_equal(y0[-i_beta, ], y[-i_beta, , i + 1])
+  }
+})

tests/testthat/test-walk.R

@@ -85,7 +85,7 @@ test_that("provided dt is reasonable", {
     "Expected 'dt' to be greater than 0")
   expect_error(
     dust_system_create(walk(), pars, n_particles = 10, dt = 1.5),
-    "Expected 'dt' to be at most 1")
+    "Expected 'dt' to be an integer, if greater than 1")
   expect_error(
     dust_system_create(walk(), pars, n_particles = 10, dt = sqrt(2) / 2),
     "Expected 'dt' to be the inverse of an integer")
@@ -111,6 +111,18 @@ test_that("time aligns to grid when dt is less than 1", {
     fixed = TRUE)
 })
 
+test_that("time aligns to grid when dt is more than 1", {
+  pars <- list(sd = 1, random_initial = TRUE)
+  obj <- dust_system_create(walk(), pars, n_particles = 10, dt = 2.0,
+                            time = 0.0)
+  expect_equal(dust_system_time(obj), 0.0)
+
+  expect_error(
+    dust_system_create(walk(), pars, n_particles = 10, dt = 2.0,
+                            time = 1.0),
+    "'time' must be a multiple of 'dt'")
+})
+
 
 ## Not a test at all of the system, but of the error handling in
 ## helpers.hpp