add flare.inf.num, flare.inf.point

DESCRIPTION

@@ -7,7 +7,7 @@ Authors@R:
     c(person("Puxue", "Qiao", role = c("aut", "cre"),
       email = "pqiao@svi.edu.au", 
       comment = c(ORCID = "0000-0002-1521-9086")),
-      person("Christina", "Azodi", role = c("aut", "cre"),
+      person("Christina", "Azodi", role = c("aut"),
       email = "cazodi@svi.edu.au",
       comment = c(ORCID = "0000-0002-6097-606X")),
       person("Ruqian", "Lyu", role = c("aut"),

R/MOD_arrivals.R

@@ -18,31 +18,38 @@ arrivals.FUN <- function(dat, at, seed = NULL) {
   # Conditions --------------------------------------------------------------
   if (!dat$param$vital) return(dat)
   
-  # Variables ---------------------------------------------------------------
-  a.rate <- dat$param$a.rate
-  a.prop.e <- dat$param$a.prop.e
-  a.prop.i <- dat$param$a.prop.i
-  a.prop.q <- dat$param$a.prop.q
-  a.rand <- dat$control$a.rand
-  groups <- dat$param$groups
-  nOld <- dat$epi$num[at - 1]
-  type <- dat$control$type
-  nsteps <- dat$control$nsteps
   
-  # Process: partition arrivals into compartments -----------------------------------------------------------------
-  nArrivals <- ifelse(a.rand, sum(stats::rbinom(nOld, 1, a.rate)), round(nOld * a.rate))
   
-  nArrivals.e <- round(nArrivals * ifelse(length(a.prop.e) > 1, a.prop.e[at], a.prop.e))
-  nArrivals.i <- round(nArrivals * ifelse(length(a.prop.i) > 1, a.prop.i[at], a.prop.i))
-  nArrivals.q <- round(nArrivals * ifelse(length(a.prop.q) > 1, a.prop.q[at], a.prop.q))
-  totArrivals <- nArrivals.e + nArrivals.i + nArrivals.q
+  flare.inf.point <- dat$param$flare.inf.point
+  flare.inf.num <- dat$param$flare.inf.num
+  
+  if(!(at %in% flare.inf.point)){
+    # Variables ---------------------------------------------------------------
+    a.rate <- dat$param$a.rate
+    a.prop.e <- dat$param$a.prop.e
+    a.prop.i <- dat$param$a.prop.i
+    a.prop.q <- dat$param$a.prop.q
+    a.rand <- dat$control$a.rand
+    nOld <- dat$epi$num[at - 1]
+    # Process: partition arrivals into compartments -----------------------------------------------------------------
+    nArrivals <- ifelse(a.rand, sum(stats::rbinom(nOld, 1, a.rate)), round(nOld * a.rate))
+    nArrivals.e <- round(nArrivals * ifelse(length(a.prop.e) > 1, a.prop.e[at], a.prop.e))
+    nArrivals.i <- round(nArrivals * ifelse(length(a.prop.i) > 1, a.prop.i[at], a.prop.i))
+    nArrivals.q <- round(nArrivals * ifelse(length(a.prop.q) > 1, a.prop.q[at], a.prop.q))
+    totArrivals <- nArrivals.e + nArrivals.i + nArrivals.q
+    dat$attr$status <- c(dat$attr$status,
+                         rep("e", nArrivals.e),
+                         rep("i", nArrivals.i),
+                         rep("q", nArrivals.q))
+    
+  }else{
+    nArrivals.i <- totArrivals <- flare.inf.num[which(flare.inf.point == at)]
+    dat$attr$status <- c(dat$attr$status, rep("i", nArrivals.i) )
+    nArrivals.e <- nArrivals.q <- 0
+  }
   
   dat$attr$active <- c(dat$attr$active, rep(1, totArrivals))
   dat$attr$group <- c(dat$attr$group, rep(1, totArrivals))
-  dat$attr$status <- c(dat$attr$status,
-                       rep("e", nArrivals.e),
-                       rep("i", nArrivals.i),
-                       rep("q", nArrivals.q))
   dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals))
   dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals))
   dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals))
@@ -50,6 +57,7 @@ arrivals.FUN <- function(dat, at, seed = NULL) {
   dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals))
   dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals))
   
+  
   # Output ------------------------------------------------------------------
   if (at == 2) {
     dat$epi$a.flow <- c(0, totArrivals)

R/PREP_crosscheck.R

@@ -56,6 +56,11 @@ crosscheck_seiqhrf <- function(param, init, control) {
         }
     }
     
+        ## Check flare parameters 
+        if(length(param$flare.inf.point) != length(param$flare.inf.num)) stop("flare.inf.point and flare.inf.num must have the same length!")
+        if(max(param$flare.inf.point) > control$nsteps) stop("flare.inf.point can not exceed the total time points: nsteps!")
+        if(min(param$flare.inf.point) <= 2) stop("flare.inf.point should all be larger than 2, please change the initial i.num if flare time point needs to be smaller than 2")
+    
     ## In-place assignment to update param and control
     on.exit(assign("param", param, pos = parent.frame()))
     on.exit(assign("control", control, pos = parent.frame()), add = TRUE)

R/PREP_param.R

@@ -103,6 +103,8 @@
 #'     fat.rand for details.
 #' @param vital Enables demographics, that is, arrivals and
 #'     departures, to and from the simulated population.
+#' @param flare.inf.point (Either a vector or a scalar) Time points where there is a sudden arrival of I's. 
+#' @param flare.inf.num (same dimension as flare.inf.point) The number of I's arriving at the specified time points in flare.inf.point.
 #' @param a.rate Background demographic arrival rate. Currently all
 #'     arrivals go into the S compartment, the default is approximately the
 #'     daily birth rate for Australia. Will be extended to cover immigration in
@@ -190,6 +192,8 @@ param_seiqhrf <- function(inf.prob.e = 0.02,
                           fat.rate.overcap = 1/25,
                           fat.tcoeff = 0.5,
                           vital = TRUE,
+                          flare.inf.point = NULL, 
+                          flare.inf.num = NULL, 
                           a.rate = (10.5/365)/1000,
                           a.prop.e = 0.01,
                           a.prop.i = 0.001,

man/param_seiqhrf.Rd

@@ -34,6 +34,8 @@ param_seiqhrf(
   fat.rate.overcap = 1/25,
   fat.tcoeff = 0.5,
   vital = TRUE,
+  flare.inf.point = NULL,
+  flare.inf.num = NULL,
   a.rate = (10.5/365)/1000,
   a.prop.e = 0.01,
   a.prop.i = 0.001,
@@ -178,6 +180,10 @@ fat.rand for details.}
 \item{vital}{Enables demographics, that is, arrivals and
 departures, to and from the simulated population.}
 
+\item{flare.inf.point}{(Either a vector or a scalar) Time points where there is a sudden arrival of I's.}
+
+\item{flare.inf.num}{(same dimension as flare.inf.point) The number of I's arriving at the specified time points in flare.inf.point.}
+
 \item{a.rate}{Background demographic arrival rate. Currently all
 arrivals go into the S compartment, the default is approximately the
 daily birth rate for Australia. Will be extended to cover immigration in