Merge pull request #1 from ku-awdc/CppExpansion

Cpp expansion to numI, numR etc

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: IPDMR
 Title: Support Code for the Introduction to Practical Disease Modelling Course
-Version: 0.4.3-1
-Date: 2024-11-07
+Version: 0.5.0-1
+Date: 2024-11-08
 Description: A collection of functions and classes that are used to illustrate
   fundamental concepts of disease modelling as part of a physical course series
   taught jointly by the University of Copenhagen and University of Sydney. The

NAMESPACE

@@ -19,6 +19,7 @@ export(si_continuous)
 export(si_discrete)
 export(sir_det)
 export(sir_stoc)
+export(suggest_boxes)
 import(R6)
 import(Rcpp)
 import(dplyr)
@@ -33,6 +34,8 @@ importFrom(methods,new)
 importFrom(pbapply,pblapply)
 importFrom(rlang,.data)
 importFrom(rlang,set_names)
+importFrom(stats,dgamma)
+importFrom(stats,qgamma)
 importFrom(stats,quantile)
 importFrom(stats,rmultinom)
 importFrom(tibble,tibble)

R/BetweenGroupClass.R

@@ -117,7 +117,14 @@ BetweenGroupClass <- R6::R6Class(
 
       ## TODO: C++ function if private$.allcpp
 
-      trans_b <- colSums(private$.beta_matrix * sapply(private$.groups, \(x) x$I))
+      if(private$.trans_between == "frequency"){
+        multby <- sapply(private$.groups, \(x){ x$I / x$N })
+      }else if(private$.trans_between == "density"){
+        multby <- sapply(private$.groups, \(x) x$I)
+      }else{
+        stop("Unrecognised private$.trans_between")
+      }
+      trans_b <- colSums(private$.beta_matrix * multby)
       for(i in seq_along(private$.groups)){
         private$.groups[[i]]$trans_external <- trans_b[i]
         private$.groups[[i]]$update(d_time)
@@ -206,6 +213,7 @@ BetweenGroupClass <- R6::R6Class(
     .ngroups = numeric(),
     .groups = list(),
     .time = numeric(),
+    .trans_between = "density",
 
     .dummy=NULL
   ),
@@ -229,6 +237,13 @@ BetweenGroupClass <- R6::R6Class(
       private$.time
     },
 
+    #' @field transmission_between the between-group transmission type (frequency or density)
+    transmission_between = function(value){
+      if(missing(value)) return(private$.trans_between)
+      value <- match.arg(value, c("frequency","density"))
+      private$.trans_between <- value
+    },
+
     #' @field state a data frame of the current state of each group
     state = function(){
       lapply(private$.groups, function(x) x$state) |>

R/SEIRclass.R

@@ -169,6 +169,8 @@ SEIRclass <- R6::R6Class("SEIRclass",
         "beta/omega/gamma/delta = ", self$beta, "/", self$omega, "/", self$gamma, "/", self$delta, "\n\t",
         "vacc/repl/cull = ", self$vacc, "/", self$repl, "/", self$cull, "\n\t",
         "E compartments = ", private$.numE, "\n\t",
+        "I compartments = 1\n\t",
+        "R compartments = 1\n\t",
         "external transmission = ", private$.trans_external, "\n\t",
         "update type = ", private$.update_type, "\n\t",
         "transmission type = ", private$.transmission_type, "\n\t",

R/make_group.R

@@ -8,9 +8,10 @@
 #' implementations should be identical).
 #'
 #' @param update_type either stochastic or deterministic
-#' @param numE the number of sub-compartments desired for the E state
+#' @param numE the number of sub-compartments desired for the E state (0 or more)
+#' @param numI the number of sub-compartments desired for the I state (1 or more)
+#' @param numR the number of sub-compartments desired for the R state (0 or more)
 #' @param group_name an optional name for the group
-#' @param model_type the compartmental model representation desired (currently only SEIR is supported)
 #' @param implementation either C++ or R6
 #'
 #' @importFrom methods new
@@ -26,35 +27,33 @@
 #' stopifnot(unlist(r6res) - unlist(cppres) == 0L)
 #'
 #' @export
-make_group <- function(update_type = c("deterministic","stochastic"), numE = 3L, group_name=NA_character_, model_type = c("SEIR", "SIR", "SI"), implementation = c("C++","R6")){
+make_group <- function(update_type = c("deterministic","stochastic"), numE = 3L, numI = 1L, numR = 1L, group_name=NA_character_, implementation = c("C++","R6")){
 
   update_type <- match.arg(update_type)
   implementation <- toupper(implementation)
   implementation <- match.arg(implementation)
-  model_type <- toupper(model_type)
-  model_type <- match.arg(model_type)
 
-  stopifnot(model_type=="SEIR")
+  qassert(numE, "X1[0,)")
+  qassert(numI, "X1(0,)")
+  qassert(numR, "X1[0,)")
 
   if(implementation=="R6"){
 
+    if(numE==0L) stop("The R6 implementation is limited to numE>=1")
+    if(numI!=1L) stop("The R6 implementation is limited to numI==1")
+    if(numR!=1L) stop("The R6 implementation is limited to numR==1")
+
     model <- SEIRclass$new(update_type=update_type, numE=numE, group_name=group_name)
 
   }else if(implementation=="C++"){
 
-    if(update_type=="deterministic"){
-      if(numE==3L){
-        model <- new(SEIRdet3, 3L, group_name)
-      }else{
-        model <- new(SEIRdetN, numE, group_name)
-      }
-    }else{
-      if(numE==3L){
-        model <- new(SEIRstoc3, 3L, group_name)
-      }else{
-        model <- new(SEIRstocN, numE, group_name)
-      }
-    }
+    modname <- str_c("SEIR",
+      if(update_type=="deterministic") "det" else "stoc",
+      if(numE %in% c(0,1,3)) numE else "N",
+      if(numI %in% c(1)) numI else "N",
+      if(numR %in% c(0,1)) numR else "N"
+    )
+    model <- new(get(modname), numE, numI, numR, group_name)
 
   }else{
     stop("Unmatched implementation type")

R/model_wrappers.R

@@ -31,11 +31,47 @@ NULL
 
 #' @rdname group_models
 #' @export
-sir_det <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type="frequency", d_time=1, max_time=100){
+sir_det <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type=c("frequency","density"), d_time=1, max_time=100){
 
-  qassert(transmission_type, "S1")
+  transmission_type <- match.arg(transmission_type)
+
+  model <- make_group(update_type="deterministic", numE=0, numI=1, numR=1)
+  #model <- WithinGroupModel$new(model_type="sir", update_type="deterministic", transmission_type=transmission_type, d_time=d_time)
+
+  model$transmission_type <- transmission_type
+  model$S <- S
+  model$I <- I
+  model$R <- R
+  model$beta <- beta
+  model$gamma <- gamma
+  model$delta <- delta
+  out <- model$run(max_time, d_time)
+  rm(model)
+
+  class(out) <- c("ipdmr_dt", class(out))
+  attr(out, "plot_caption") <- str_c("deterministic; discrete; ", transmission_type)
+  return(out)
+
+  ### OLDER CODE
+
+
+  update_type <- "deterministic"
+  if(update_type=="deterministic"){
+
+
+  }else{
+
+    out |>
+      mutate(Iteration = 1L) |>
+      select(Iteration, everything()) ->
+      out
+
+    class(out) <- c("ipdmr_st", class(out))
+    attr(out,'iterations') <- 1L
+    attr(out, "plot_caption") <- str_c("stochastic; discrete; ", transmission_type)
+
+  }
 
-  model <- WithinGroupModel$new(model_type="sir", update_type="deterministic", transmission_type=transmission_type, d_time=d_time)
 
   model$S <- S
   model$I <- I
@@ -54,16 +90,11 @@ sir_det <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmissi
 #' @export
 seir_det <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_type=c("frequency","density"), d_time=1, max_time=100){
 
-  #qassert(transmission_type, "S1")
   transmission_type <- match.arg(transmission_type)
 
-  model <- make_group(update_type="deterministic", numE=numE, group_name=NA_character_, model_type="SEIR", implementation="C++")
+  model <- make_group(update_type="deterministic", numE=numE, numI=1, numR=1)
   model$transmission_type <- transmission_type
 
-  if(FALSE){
-    model <- WithinGroupModel$new(model_type="seir", update_type="deterministic", transmission_type=transmission_type, numE=numE, d_time=d_time)
-  }
-
   model$S <- S
   model$E <- E
   model$I <- I
@@ -77,7 +108,7 @@ seir_det <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0
   model$repl <- repl
   model$cull <- cull
 
-  model$run(ceiling(max_time/d_time), d_time) ->
+  model$run(max_time, d_time) ->
     output
 
   class(output) <- c("ipdmr_dt", class(output))
@@ -89,23 +120,27 @@ seir_det <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0
 
 #' @rdname group_models
 #' @export
-sir_stoc <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type="frequency", d_time=1, max_time=100, iterations=1L){
+sir_stoc <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type=c("frequency","density"), d_time=1, max_time=100, iterations=1L){
 
-  qassert(transmission_type, "S1")
+  transmission_type <- match.arg(transmission_type)
+  qassert(iterations, "X1")
 
-  pblapply(seq_len(iterations), \(i){
+  model <- make_group(update_type="stochastic", numE=0, numI=1, numR=1)
 
-    model <- WithinGroupModel$new(model_type="sir", update_type="stochastic", transmission_type=transmission_type, d_time=d_time)
+  model$transmission_type <- transmission_type
+  model$S <- S
+  model$I <- I
+  model$R <- R
+  model$beta <- beta
+  model$gamma <- gamma
+  model$delta <- delta
+  model$save()
 
-    model$S <- S
-    model$I <- I
-    model$R <- R
+  pblapply(seq_len(iterations), \(i){
 
-    model$beta <- beta
-    model$gamma <- gamma
-    model$delta <- delta
+    model$reset()
 
-    model$run(ceiling(max_time/d_time), d_time) |>
+    model$run(max_time, d_time) |>
       as_tibble() |>
       mutate(Iteration = i) |>
       select("Iteration", everything())
@@ -124,12 +159,12 @@ sir_stoc <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmiss
 
 #' @rdname group_models
 #' @export
-seir_stoc <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_type="frequency", d_time=1, max_time=100, iterations=1L){
+seir_stoc <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_type=c("frequency","density"), d_time=1, max_time=100, iterations=1L){
 
   #qassert(transmission_type, "S1")
   transmission_type <- match.arg(transmission_type)
 
-  model <- make_group(update_type="stochastic", numE=numE, group_name=NA_character_, model_type="SEIR", implementation="C++")
+  model <- make_group(update_type="stochastic", numE=numE, numI=1, numR=1)
   model$transmission_type <- transmission_type
   model$S <- S
   model$E <- E
@@ -149,7 +184,7 @@ seir_stoc <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=
     #model <- WithinGroupModel$new(model_type="seir", update_type="stochastic", transmission_type=transmission_type, num_E=num_E, d_time=d_time)
 
     model$reset()
-    model$run(ceiling(max_time/d_time), d_time) |>
+    model$run(max_time, d_time) |>
       as_tibble() |>
       mutate(Iteration = i) |>
       select("Iteration", everything())

R/multi_wrappers.R

@@ -22,6 +22,8 @@
 #' @param vacc the vaccination rate parameter per unit time (must be positive)
 #' @param repl the replacement rate parameter per unit time (must be positive)
 #' @param cull the targeted culling rate parameter per unit time (must be positive)
+#' @param transmission_within the within-group transmission type (frequency or density)
+#' @param transmission_between the between-group transmission type (frequency or density)
 #' @param iterations the number of iterations to run (stochastic models only)
 #' @param d_time the desired time step (delta time)
 #' @param max_time the desired maximum time point (must be greater than the time step)
@@ -37,9 +39,11 @@ NULL
 
 #' @rdname multi_models
 #' @export
-multi_seir_det <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, d_time=1, max_time=100){
+multi_seir_det <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_within=c("frequency","density"), transmission_between=c("density","frequency"), d_time=1, max_time=100){
 
-  output <- multi_wrapper("deterministic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, d_time=d_time, max_time=max_time, iterations=1L)
+  transmission_within <- match.arg(transmission_within)
+  transmission_between <- match.arg(transmission_between)
+  output <- multi_wrapper("deterministic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, transmission_within=transmission_within, transmission_between=transmission_between, d_time=d_time, max_time=max_time, iterations=1L)
 
   class(output) <- c("ipdmr_dm", class(output))
   attr(output, "plot_caption") <- str_c("deterministic; ", n_groups, " groups")
@@ -49,9 +53,11 @@ multi_seir_det <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L,
 
 #' @rdname multi_models
 #' @export
-multi_seir_stoc <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, d_time=1, max_time=100, iterations=1L){
+multi_seir_stoc <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_within=c("frequency","density"), transmission_between=c("density","frequency"), d_time=1, max_time=100, iterations=1L){
 
-  output <- multi_wrapper("stochastic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, d_time=d_time, max_time=max_time, iterations=iterations)
+  transmission_within <- match.arg(transmission_within)
+  transmission_between <- match.arg(transmission_between)
+  output <- multi_wrapper("stochastic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, transmission_within=transmission_within, transmission_between=transmission_between, d_time=d_time, max_time=max_time, iterations=iterations)
 
   class(output) <- c("ipdmr_sm", class(output))
   attr(output, "iterations") <- iterations
@@ -61,9 +67,11 @@ multi_seir_stoc <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L,
 }
 
 
-multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, d_time=1, max_time=100L, iterations=1L){
+multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_within=c("frequency","density"), transmission_between=c("density","frequency"), d_time=1, max_time=100L, iterations=1L){
 
   ## Check arguments:
+  transmission_within <- match.arg(transmission_within)
+  transmission_between <- match.arg(transmission_between)
   qassert(n_groups, "X1(0,)")
   assert_matrix(beta_matrix, "numeric", any.missing=FALSE, nrows=n_groups, ncols=n_groups)
   qassert(iterations, "X1(0,)")
@@ -73,6 +81,7 @@ multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=
   model_pars <- list(
     S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega,
     gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull,
+    transmission_type=transmission_within,
     group_number=seq_len(n_groups)
   )
 
@@ -89,9 +98,7 @@ multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=
     group_split() |>
     set_names(str_c("Group_", seq_len(n_groups) |> format() |> str_replace_all(" ", "0"))) |>
     lapply(\(x){
-      #md <- WithinGroupModel$new("seir", update_type, "frequency", numE=x$numE, d_time=d_time, group_number=x[["group_number"]])
-      md <- make_group(update_type=update_type, numE=x$numE, group_name=x[["group_number"]], model_type="SEIR", implementation="C++")
-      md$transmission_type <- "frequency"
+      md <- make_group(update_type=update_type, numE=x$numE, numI=1, numR=1, group_name=x[["group_number"]])
       x$numE <- NULL
       x$group_number <- NULL
       for(nm in names(x)){
@@ -105,6 +112,7 @@ multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=
   ##  Set up the model:
   model <- BetweenGroupClass$new(models)
   model$beta_matrix <- beta_matrix
+  model$transmission_between <- transmission_between
   model$save()
 
   ## Iterate:

R/si_discrete.R

@@ -28,6 +28,23 @@ si_discrete <- function(S=9, I=1, beta=0.05, transmission_type=c("frequency","de
   ntime <- ceiling(max_time/d_time)+1L
   type <- match.arg(transmission_type)
 
+  model <- make_group("deterministic", numE=0L, numI=1L, numR=0L)
+  model$transmission_type <- transmission_type
+  model$S <- S
+  model$I <- I
+  model$beta <- beta
+  model$gamma <- 0
+  model$repl <- 0
+  model$cull <- 0
+  output <- model$run(max_time, d_time)
+  rm(model)
+
+  class(output) <- c("ipdmr_dt", class(output))
+  attr(output, "plot_caption") <- str_c("discrete; ", type, "; d_time=", round(d_time,3))
+  return(output)
+
+  ### OLDER CODE BELOW HERE
+
   i_dens <- function(i) 1 - exp(-beta * i)
   i_freq <- function(i) 1 - exp(-beta * i/N)
   i_fun <- if(type=="density") i_dens else i_freq