tests/SimInf_model.R

## This file is part of SimInf, a framework for stochastic
## disease spread simulations.
##
## Copyright (C) 2015 Pavol Bauer
## Copyright (C) 2017 -- 2019 Robin Eriksson
## Copyright (C) 2015 -- 2019 Stefan Engblom
## Copyright (C) 2015 -- 2024 Stefan Widgren
##
## SimInf is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## SimInf is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <https://www.gnu.org/licenses/>.

library(SimInf)
library(Matrix)
library(tools)
source("util/check.R")

## Specify the number of threads to use.
set_num_threads(1)

## For debugging
sessionInfo()

## Initialize test data
S <- as(Matrix(c(-1,  0,  0,
                  1,  0,  0,
                  0, -1,  0,
                  0,  1,  0,
                  0,  0, -1,
                  0,  0,  1),
               nrow   = 6,
               ncol   = 3,
               byrow  = TRUE,
               dimnames = list(LETTERS[1:6], NULL)),
        "dgCMatrix")

Nn <- 6L

G <- as(Matrix(c(1, 0, 0,
                 0, 1, 0,
                 0, 0, 1),
               nrow   = 3,
               ncol   = 3,
               byrow  = TRUE,
               dimnames = list(c("A -> B", "C -> D", "E -> F"), NULL)),
        "dgCMatrix")

u0 <- structure(c(0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 2, 0, 2, 0, 2,
                  0, 3, 0, 3, 0, 3, 0, 4, 0, 4, 0, 4, 0, 5, 0, 5, 0, 5, 0),
                .Dim = c(6L, 6L))
storage.mode(u0) <- "integer"

U <- matrix(nrow = 0, ncol = 0)
storage.mode(U) <- "integer"

## Check valid_SimInf_model_object
m <- SISe(u0 = c(S = 10, I = 0), tspan = 1:10, phi = 0,
          upsilon = 0.1, gamma = 0.1, alpha = 1.0, beta_t1 = 0.1,
          beta_t2 = 0.1, beta_t3 = 0.1, beta_t4 = 0.1, end_t1  = 91,
          end_t2  = 182, end_t3  = 273, end_t4  = 365, epsilon = 0.1)
stopifnot(isTRUE(SimInf:::valid_SimInf_model_object(m)))

## Check valid_SimInf_model_object with invalid tspan.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@tspan <- integer(0)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Input time-span must be a double vector."))
m@tspan <- numeric(0)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Input time-span must be an increasing vector."))

## Check valid_SimInf_model_object with invalid u0.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
storage.mode(m@u0) <- "double"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Initial state 'u0' must be an integer matrix."))
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@u0[1, 1] <- -1L
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Initial state 'u0' has negative elements."))

## Check valid_SimInf_model_object with invalid U.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
storage.mode(m@U) <- "double"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Output state 'U' must be an integer matrix."))
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@U <- matrix(-1L)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Output state 'U' has negative elements."))

## Check valid_SimInf_model_object with invalid v0.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@v0 <- matrix(1L)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Initial model state 'v0' must be a double matrix."))
m@v0 <- matrix(1)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'v0' must have rownames."))
m@v0 <- matrix(c(1, 1), ncol = 2)
rownames(m@v0) <- "test"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "The number of nodes in 'u0' and 'v0' must match."))

## Check valid_SimInf_model_object with invalid V.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
storage.mode(m@V) <- "integer"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Output model state 'V' must be a double matrix."))

## Check valid_SimInf_model_object with invalid S.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@S@x <- m@S@x * 0.5
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'S' matrix must be an integer matrix."))

## Check valid_SimInf_model_object with different rownames in S and E.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1,
         events = data.frame(event = 1, node = 1, n = 1, time = 1, dest = 0,
                             proportion = 0, select = 1, shift = 0))
rownames(m@events@E) <- NULL
stopifnot(identical(
    SimInf:::valid_SimInf_model_object(m),
    "'S' and 'E' must have rownames matching the compartments."))
rownames(m@events@E) <- rownames(m@S)
rownames(m@events@E)[1] <- "Z"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'S' and 'E' must have identical compartments."))

## Check valid_SimInf_model_object with invalid G.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@G <- m@G[1, 1:2, drop = FALSE]
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "Wrong size of dependency graph."))

m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
rownames(m@G) <- NULL
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'G' must have rownames that specify transitions."))

m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
rownames(m@G)[1] <- ""
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'G' must have rownames that specify transitions."))

m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
rownames(m@G)[1] <- "A"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'G' rownames have invalid transitions."))

m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
rownames(m@G)[1] <- "Z -> beta*Z*I/(S+I+R) -> I"
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'G' and 'S' must have identical compartments."))

## Check valid_SimInf_model_object with invalid ldata.
m <- SIR(u0 = data.frame(S = 10, I = 0, R = 0),
         tspan = 1:10, beta = 0.1, gamma = 0.1)
m@ldata <- matrix(1L)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'ldata' matrix must be a double matrix."))
m@ldata <- matrix(c(1, 2), ncol = 2)
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "The number of nodes in 'u0' and 'ldata' must match."))

## Check valid_SimInf_model_object with invalid gdata.
m <- SISe(u0 = data.frame(S = 10, I = 0), tspan = 1:10, phi = 0,
          upsilon = 0.1, gamma = 0.1, alpha = 1.0, beta_t1 = 0.1,
          beta_t2 = 0.1, beta_t3 = 0.1, beta_t4 = 0.1, end_t1  = 91,
          end_t2  = 182, end_t3  = 273, end_t4  = 365, epsilon = 0.1)
m@gdata <- as.integer(m@gdata)
names(m@gdata) <- c("upsilon", "gamma", "alpha", "beta_t1",
                    "beta_t2", "beta_t3", "beta_t4", "epsilon")
stopifnot(identical(SimInf:::valid_SimInf_model_object(m),
                    "'gdata' must be a double vector."))

## Check v0
res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = U,
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = c(1, 2),
                                u0    = u0,
                                v0    = 1))
check_error(res, "'v0' must be a double matrix.")
res <- SimInf_model(G     = G,
                    S     = S,
                    U     = U,
                    ldata = matrix(rep(0, Nn), nrow = 1),
                    tspan = c(1, 2),
                    u0    = u0,
                    v0    = matrix(rep(1L, Nn), nrow = 1,
                                   dimnames = list("phi")))

## Check tspan
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = numeric(0),
                       u0    = u0))
check_error(res, "Input time-span must be an increasing vector.", FALSE)

res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(c(3, 2, 1)),
                       u0    = u0))
check_error(res, "Input time-span must be an increasing vector.", FALSE)

res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = c(1, NA, 3),
                       u0    = u0))
check_error(res, "Input time-span must be an increasing vector.", FALSE)

## Check that tspan can be a Date vector
res <- SimInf_model(G     = G,
                    S     = S,
                    U     = U,
                    ldata = matrix(rep(0, Nn), nrow = 1),
                    tspan = as.Date(c("2017-01-01",
                                      "2017-01-02",
                                      "2017-01-03")),
                    u0    = u0)
stopifnot(identical(res@tspan,
                    structure(c(1, 2, 3),
                              .Names = c("2017-01-01",
                                         "2017-01-02",
                                         "2017-01-03"))))

## Check u0
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0 * -1L))
check_error(res, "Initial state 'u0' has negative elements.", FALSE)

## Change storage mode of u0 to double.
## Should not raise error
u0_double <- u0
storage.mode(u0_double) <- "double"
SimInf_model(G     = G,
             S     = S,
             U     = U,
             ldata = matrix(rep(0, Nn), nrow = 1),
             tspan = as.numeric(1:10),
             u0    = u0_double)

## Change storage mode of u0 to double and change to non-integer values.
## Should raise error
u0_double <- u0
storage.mode(u0_double) <- "double"
u0_double <- 1.2 * u0_double
res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = U,
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = as.numeric(1:10),
                                u0    = u0_double))
check_error(res, "'u0' must be an integer matrix.")

## Change u0 to vector. Should raise an error.
res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = U,
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = as.numeric(1:10),
                                u0    = as.numeric(u0)))
check_error(res, "'u0' must be an integer matrix.")

## Check S
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S * 1.1,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'S' matrix must be an integer matrix.", FALSE)

## Check G
## Error: Wrong size of dependency graph
res <- assertError(new("SimInf_model",
                       G     = G[-1, ],
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "Wrong size of dependency graph.", FALSE)

## Check specication of transition
rownames(G) <- NULL
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'G' must have rownames that specify transitions.", FALSE)

rownames(G) <- c("", "  ", "E -> F")
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'G' must have rownames that specify transitions.", FALSE)

rownames(G) <- c("A -> B", "C -> D", "E ->")
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'G' rownames have invalid transitions.", FALSE)

rownames(G) <- c("A -> B", "C -> D", "E ->")
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'G' rownames have invalid transitions.", FALSE)

rownames(G) <- c("A -> B", "C -> D", "E -> G")
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'G' and 'S' must have identical compartments", FALSE)
rownames(G) <- c("A -> B", "C -> D", "E -> F")

## Check gdata
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = matrix(rep(0, Nn), nrow = 1),
                       gdata = 1L,
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "'gdata' must be a double vector.", FALSE)

## Check ldata
ldata <- matrix(rep(0, Nn), nrow = 1)
ldata <- ldata[, 1:3, drop = FALSE]

## Wrong size of ldata matrix
res <- assertError(new("SimInf_model",
                       G     = G,
                       S     = S,
                       U     = U,
                       ldata = ldata,
                       tspan = as.numeric(1:10),
                       u0    = u0))
check_error(res, "The number of nodes in 'u0' and 'ldata' must match.", FALSE)

## Check initial state
u0 <- data.frame(S_1 = c(0, 1, 2, 3, 4, 5),
                 I_1 = c(0, 0, 0, 0, 0, 0),
                 S_2 = c(0, 1, 2, 3, 4, 5),
                 I_2 = c(0, 0, 0, 0, 0, 0),
                 S_3 = c(0, 1, 2, 3, 4, 5),
                 I_3 = c(0, 0, 0, 0, 0, 0))

## 'u0' is NULL
res <- assertError(SimInf_model())
check_error(res, "'u0' must be an integer matrix.")

## Check first lines of show method without events
model <- SISe(u0      = data.frame(S = 99, I = 1),
              tspan   = seq_len(1000) - 1,
              events  = NULL,
              phi     = 1,
              upsilon = 0.017,
              gamma   = 0.1,
              alpha   = 1,
              beta_t1 = 0.19,
              beta_t2 = 0.085,
              beta_t3 = 0.075,
              beta_t4 = 0.185,
              end_t1  = 91,
              end_t2  = 182,
              end_t3  = 273,
              end_t4  = 365,
              epsilon = 0.000011)

show_expected <- c("Model: SISe",
                   "Number of nodes: 1",
                   "Number of transitions: 2",
                   "Number of scheduled events: 0",
                   "")

stopifnot(identical(show(model), model))
show_observed <- capture.output(show(model))

stopifnot(identical(show_observed[1:5], show_expected))

## Check summary method without events
summary(run(model))

## Check first lines of show method with events
u0 <- data.frame(S_1 = c(0, 1, 2, 3, 4, 5),
                 I_1 = c(0, 0, 0, 0, 0, 0),
                 S_2 = c(0, 1, 2, 3, 4, 5),
                 I_2 = c(0, 0, 0, 0, 0, 0),
                 S_3 = c(0, 1, 2, 3, 4, 5),
                 I_3 = c(0, 0, 0, 0, 0, 0))

events <- data.frame(
    event      = c(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3),
    time       = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
    node       = c(2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6),
    dest       = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
    n          = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5),
    proportion = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
    select     = c(4, 5, 6, 4, 5, 6, 4, 5, 6, 4, 5, 6, 4, 5, 6),
    shift      = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

model <- SISe3(u0        = u0,
               tspan     = 0:10,
               events    = events,
               phi       = rep(1, 6),
               upsilon_1 = 1,
               upsilon_2 = 1,
               upsilon_3 = 1,
               gamma_1   = 1,
               gamma_2   = 1,
               gamma_3   = 1,
               alpha     = 1,
               beta_t1   = 1,
               beta_t2   = 1,
               beta_t3   = 1,
               beta_t4   = 1,
               end_t1    = 91,
               end_t2    = 182,
               end_t3    = 273,
               end_t4    = 365,
               epsilon   = 1)

show_observed <- capture.output(show(model))

stopifnot(identical(
    show_observed[1:5],
    c("Model: SISe3",
      "Number of nodes: 6",
      "Number of transitions: 6",
      "Number of scheduled events: 15",
      "")))

stopifnot(identical(
    show_observed[22:28],
    c("Local data",
      "----------",
      " Parameter Value",
      " end_t1     91  ",
      " end_t2    182  ",
      " end_t3    273  ",
      " end_t4    365  ")))

## Check summary method with events
summary(run(model))

## Check U. Change storage mode of U to double.
## Should not raise error
U <- structure(c(
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 2L, 0L, 1L, 1L,
    1L, 1L, 2L, 1L, 3L, 0L, 2L, 1L, 2L, 2L, 0L, 4L, 1L, 3L, 2L, 3L,
    3L, 2L, 1L, 4L, 6L, 9L, 7L, 8L, 4L, 11L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 7L, 8L, 7L, 8L, 5L, 10L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 8L, 7L, 6L, 9L,
    8L, 7L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    5L, 10L, 4L, 11L, 6L, 9L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 7L, 8L, 5L, 10L, 7L, 8L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 11L, 5L, 10L, 3L, 12L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 8L, 7L, 5L,
    10L, 4L, 11L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 6L, 9L, 2L, 13L, 4L, 11L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 5L, 10L, 2L, 13L, 7L, 8L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 9L, 6L, 2L, 13L, 6L, 9L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
    .Dim = c(36L, 11L))

U_double <- U
storage.mode(U_double) <- "double"

SimInf_model(G     = G,
             S     = S,
             U     = U_double,
             ldata = matrix(rep(0, Nn), nrow = 1),
             tspan = as.numeric(1:10),
             u0    = u0)

## Check U. Change storage mode of U to double and change to non-integer values.
## Should raise error
U_double <- U
storage.mode(U_double) <- "double"
U_double <- U_double * 1.2

res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = U_double,
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = as.numeric(1:10),
                                u0    = u0))
check_error(res, "'U' must be an integer matrix.")

## Check U. Should not raise an error if U is an integer vector of length 0
SimInf_model(G     = G,
             S     = S,
             U     = integer(0),
             ldata = matrix(rep(0, Nn), nrow = 1),
             tspan = as.numeric(1:10),
             u0    = u0)

## Check U. Should raise error if U is an integer vector of length > 0
res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = c(1L),
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = as.numeric(1:10),
                                u0    = u0))
check_error(res, "'U' must be equal to a 0 x 0 matrix.")

## Check V. Change storage mode of V to double.
## Should not raise error
V <- structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                 1, 1, 1, 1, 1, 1, 1, 1, 1),
               .Dim = c(6L, 10L))

V_integer <- V
storage.mode(V) <- "integer"

SimInf_model(G     = G,
             S     = S,
             U     = U,
             V     = V_integer,
             ldata = matrix(rep(0, Nn), nrow = 1),
             tspan = as.numeric(1:10),
             u0    = u0)

## Check V. Change storage mode of V to character
## Should raise error
V_character <- V
storage.mode(V_character) <- "character"

res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = U,
                                V     = V_character,
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = as.numeric(1:10),
                                u0    = u0))
check_error(res, "'V' must be a double matrix.")

## Check V. Should raise error if V is a vector of length > 0
res <- assertError(SimInf_model(G     = G,
                                S     = S,
                                U     = U,
                                V     = 1,
                                ldata = matrix(rep(0, Nn), nrow = 1),
                                tspan = as.numeric(1:10),
                                u0    = u0))
check_error(res, "'V' must be equal to a 0 x 0 matrix.")

## Check V. Should not raise an error if V is an integer vector of length 0
SimInf_model(G     = G,
             S     = S,
             U     = U,
             V     = integer(0),
             ldata = matrix(rep(0, Nn), nrow = 1),
             tspan = as.numeric(1:10),
             u0    = u0)

## Check that plot raises an error if the model hasn't run
model <- SISe(u0      = data.frame(S = 99, I = 1),
              tspan   = seq_len(1000) - 1,
              events  = NULL,
              phi     = 1,
              upsilon = 0.017,
              gamma   = 0.1,
              alpha   = 1,
              beta_t1 = 0.19,
              beta_t2 = 0.085,
              beta_t3 = 0.075,
              beta_t4 = 0.185,
              end_t1  = 91,
              end_t2  = 182,
              end_t3  = 273,
              end_t4  = 365,
              epsilon = 0.000011)
res <- assertError(plot(model))
check_error(res, "Please run the model first, the trajectory is empty.")

## Check that the SimInf_model initialisation fails if the events
## argument is not either NULL or a data.frame
u0 <- data.frame(S_1 = c(0, 1, 2, 3, 4, 5),
                 I_1 = c(0, 0, 0, 0, 0, 0),
                 S_2 = c(0, 1, 2, 3, 4, 5),
                 I_2 = c(0, 0, 0, 0, 0, 0),
                 S_3 = c(0, 1, 2, 3, 4, 5),
                 I_3 = c(0, 0, 0, 0, 0, 0))

res <- assertError(model <- SISe3(u0        = u0,
                                  tspan     = 0:10,
                                  events    = "events",
                                  phi       = rep(1, 6),
                                  upsilon_1 = 1,
                                  upsilon_2 = 1,
                                  upsilon_3 = 1,
                                  gamma_1   = 1,
                                  gamma_2   = 1,
                                  gamma_3   = 1,
                                  alpha     = 1,
                                  beta_t1   = 1,
                                  beta_t2   = 1,
                                  beta_t3   = 1,
                                  beta_t4   = 1,
                                  end_t1    = 91,
                                  end_t2    = 182,
                                  end_t3    = 273,
                                  end_t4    = 365,
                                  epsilon   = 1))
check_error(res, "'events' must be NULL or a data.frame.")

## Check arguments to 'trajectory' method
u0 <- data.frame(S = 100:105, I = 1:6, R = rep(0, 6))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)
result <- run(model)
res <- assertError(trajectory(result, compartments = c("A", "S")))
check_error(res, "Non-existing compartment(s) in model: 'A'.")
res <- assertError(trajectory(result, index = c("A", "S")))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")
res <- assertError(trajectory(result, index = 3.4))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")
res <- assertError(trajectory(result, index = 0))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")
res <- assertError(trajectory(result, index = 10))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")

## Check arguments to 'prevalence' method
u0 <- data.frame(S = 100:105, I = 1:6, R = rep(0, 6))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)
result <- run(model)
res <- assertError(prevalence(result, A + S ~ S))
check_error(res, "Non-existing compartment(s) in model: 'A'.")
res <- assertError(prevalence(result, S ~ A + S))
check_error(res, "Non-existing compartment(s) in model: 'A'.")
res <- assertError(prevalence(result, I ~ S + I + R, i = c("A", "S")))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")
res <- assertError(prevalence(result, I ~ S + I + R, i = 3.4))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")
res <- assertError(prevalence(result, I ~ S + I + R, i = 0))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")
res <- assertError(prevalence(result, I ~ S + I + R, i = 10))
check_error(
    res,
    "The node index must be an integer > 0 and <= number of nodes.")

## Check 'gdata'
model <- SISe(u0 = data.frame(S = 10, I = 0), tspan = 1:10, phi = 0,
              upsilon = 1, gamma = 2, alpha = 3, beta_t1 = 4,
              beta_t2 = 5, beta_t3 = 6, beta_t4 = 7, end_t1  = 91,
              end_t2  = 182, end_t3  = 273, end_t4  = 365, epsilon = 8)

stopifnot(identical(gdata(model),
                    c(upsilon = 1, gamma = 2, alpha = 3, beta_t1 = 4,
                      beta_t2 = 5, beta_t3 = 6, beta_t4 = 7, epsilon = 8)))

gdata(model, "epsilon") <- 9

stopifnot(identical(gdata(model),
                    c(upsilon = 1, gamma = 2, alpha = 3, beta_t1 = 4,
                      beta_t2 = 5, beta_t3 = 6, beta_t4 = 7, epsilon = 9)))

res <- assertError(gdata(model) <- 6)
check_error(res, "Missing 'parameter' argument.")

res <- assertError(gdata(model, "epsilon") <- "6")
check_error(res, "'value' argument must be a numeric.")

res <- assertError(gdata(model, 5) <- 6)
check_error(res, "'parameter' argument must be a character.")

res <- assertError("gdata<-" (model, "epsilon"))
check_error(res, "Missing 'value' argument.")

## Check 'ldata'
model@ldata <- matrix(1, dimnames = list("test", NULL))
res <- assertError(ldata(model))
check_error(res, "Missing 'node' argument.")

res <- assertError(ldata(model, "0"))
check_error(res, "Invalid 'node' argument.")

res <- assertError(ldata(model, 0))
check_error(res, "Invalid 'node' argument.")

res <- assertError(ldata(model, c(0, 0)))
check_error(res, "Invalid 'node' argument.")

ldata(model, 1)

Try the SimInf package in your browser

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

SimInf documentation built on June 22, 2024, 9:27 a.m.