Testing custom state object

library(ppmcmc)

We can create a model state object in R that contains the state of each model element (compartment) at a single time point.

# Age groups
na <- 5
# Heterogeneity groups
nh <- 5
# New state
s <- state$new()
s$S <- matrix(0L, ncol = na, nrow = nh)
s$I <- matrix(0L, ncol = na, nrow = nh)
s$R <- matrix(0L, ncol = na, nrow = nh)

The state can have custom methods that under the hood call c++ functionality:

s$I
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    0    0    0    0
#> [2,]    0    0    0    0    0
#> [3,]    0    0    0    0    0
#> [4,]    0    0    0    0    0
#> [5,]    0    0    0    0    0
s$update()
s$I
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    0    0    0    0
#> [2,]    1    0    0    0    0
#> [3,]    1    0    0    0    0
#> [4,]    1    0    0    0    0
#> [5,]    1    0    0    0    0

This maintains passing by reference, so is fast:

s$I
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    0    0    0    0
#> [2,]    1    0    0    0    0
#> [3,]    1    0    0    0    0
#> [4,]    1    0    0    0    0
#> [5,]    1    0    0    0    0
system.time({
  for(i in 1:10000){
    s$update()
  }
})
#>    user  system elapsed 
#>   0.193   0.030   0.223
s$I
#>       [,1] [,2] [,3] [,4] [,5]
#> [1,] 10001    0    0    0    0
#> [2,] 10001    0    0    0    0
#> [3,] 10001    0    0    0    0
#> [4,] 10001    0    0    0    0
#> [5,] 10001    0    0    0    0

For the pMCMC, we will need to be able to copy particles. To do this we need to force a copy of the state in memory, which can be specified with a custom deep clone method.

# Without a deep clone, updating s, also changes s2
s2 <- s$clone()
s2$I
#>       [,1] [,2] [,3] [,4] [,5]
#> [1,] 10001    0    0    0    0
#> [2,] 10001    0    0    0    0
#> [3,] 10001    0    0    0    0
#> [4,] 10001    0    0    0    0
#> [5,] 10001    0    0    0    0
s$update()
s2$I
#>       [,1] [,2] [,3] [,4] [,5]
#> [1,] 10002    0    0    0    0
#> [2,] 10002    0    0    0    0
#> [3,] 10002    0    0    0    0
#> [4,] 10002    0    0    0    0
#> [5,] 10002    0    0    0    0

# With a deep clone it doesn't
s3 <- s$clone(deep = TRUE)
s3$I
#>       [,1] [,2] [,3] [,4] [,5]
#> [1,] 10002    0    0    0    0
#> [2,] 10002    0    0    0    0
#> [3,] 10002    0    0    0    0
#> [4,] 10002    0    0    0    0
#> [5,] 10002    0    0    0    0
s$update()
s3$I
#>       [,1] [,2] [,3] [,4] [,5]
#> [1,] 10002    0    0    0    0
#> [2,] 10002    0    0    0    0
#> [3,] 10002    0    0    0    0
#> [4,] 10002    0    0    0    0
#> [5,] 10002    0    0    0    0