Commit 4a2b2f56 authored by pqiao29's avatar pqiao29
Browse files

fix *.flow track, push version to 0.3.3

parent 0f90bc95
Pipeline #8776 passed with stages
in 7 minutes and 33 seconds
Package: sirplus
Type: Package
Title: Using stochastic individual compartment models (ICMs) to simulate COVID-19 spread
Version: 0.3.2
Version: 0.3.3
Date: 2020-07-07
Authors@R:
c(person("Puxue", "Qiao", role = c("aut", "cre"),
......@@ -57,7 +57,7 @@ Suggests:
devtools,
covr
biocViews: Epidemiology, Software
RoxygenNote: 7.1.0
RoxygenNote: 7.1.1
VignetteBuilder: knitr
Encoding: UTF-8
Language: en-GB
......@@ -57,7 +57,7 @@ seiqhrf <- function(init, control, param, priors = NULL, seed = NULL) {
# Timestep loop
for (at in 2:nsteps) {
## Infection
## Infection: S -> E
if (!is.null(control[["infection.FUN"]])) {
dat <- do.call(control[["infection.FUN"]], list(dat, at, seed))
}
......
......@@ -13,7 +13,7 @@ update_status <- function(rate, rand, active, status, label, state, at, expTime
vecProg <- which(stats::rbinom(nElig, 1, rate) == 1)
if (length(vecProg) > 0) {
at_idx <- idsElig[vecProg]
smp_sz <- length(at_idx)
smp_sz <- length(vecProg)
status[at_idx] <- state
}
}else{
......
......@@ -57,16 +57,18 @@ arrivals.FUN <- function(dat, at, seed = NULL) {
# Output ------------------------------------------------------------------
if (at == 2) {
dat$epi$a.flow <- c(0, totArrivals)
dat$epi$a.e.flow <- c(0, nArrivals.e)
dat$epi$a.i.flow <- c(0, nArrivals.i)
dat$epi$a.q.flow <- c(0, nArrivals.q)
} else {
dat$epi$a.flow[at] <- totArrivals
dat$epi$a.e.flow[at] <- c(0, nArrivals.e)
dat$epi$a.i.flow[at] <- c(0, nArrivals.i)
dat$epi$a.q.flow[at] <- c(0, nArrivals.q)
if(dat$control$track.arrival){
if (at == 2) {
dat$epi$a.flow <- c(0, totArrivals)
dat$epi$a.e.flow <- c(0, nArrivals.e)
dat$epi$a.i.flow <- c(0, nArrivals.i)
dat$epi$a.q.flow <- c(0, nArrivals.q)
} else {
dat$epi$a.flow[at] <- totArrivals
dat$epi$a.e.flow[at] <- c(0, nArrivals.e)
dat$epi$a.i.flow[at] <- c(0, nArrivals.i)
dat$epi$a.q.flow[at] <- c(0, nArrivals.q)
}
}
return(dat)
......
......@@ -28,42 +28,54 @@ departures.FUN <- function(dat, at, seed = NULL) {
nDepartures <- res[[1]]
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0
if (at == 2) dat$epi$ds.flow <- c(0, nDepartures) else dat$epi$ds.flow[at] <- nDepartures
if(dat$control$track.departure){
if (at == 2) dat$epi$sd.flow <- c(0, nDepartures) else dat$epi$sd.flow[at] <- nDepartures
}
# Exposed Departures ---------------------------------------------------------
res <- update_active(rate, rand, active, status, label = "e")
nDepartures <- res[[1]]
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0
if (at == 2) dat$epi$de.flow <- c(0, nDepartures) else dat$epi$de.flow[at] <- nDepartures
if(dat$control$track.departure){
if (at == 2) dat$epi$ed.flow <- c(0, nDepartures) else dat$epi$ed.flow[at] <- nDepartures
}
# Infected Departures ---------------------------------------------------------
res <- update_active(rate, rand, active, status, label = "i")
nDepartures <- res[[1]]
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0
if (at == 2) dat$epi$di.flow <- c(0, nDepartures) else dat$epi$di.flow[at] <- nDepartures
if(dat$control$track.departure){
if (at == 2) dat$epi$id.flow <- c(0, nDepartures) else dat$epi$id.flow[at] <- nDepartures
}
# Quarantined Departures ---------------------------------------------------------
res <- update_active(rate, rand, active, status, label = "q")
nDepartures <- res[[1]]
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0
if (at == 2) dat$epi$dq.flow <- c(0, nDepartures) else dat$epi$dq.flow[at] <- nDepartures
if(dat$control$track.departure){
if (at == 2) dat$epi$qd.flow <- c(0, nDepartures) else dat$epi$qd.flow[at] <- nDepartures
}
# Hospitalised Departures ---------------------------------------------------------
res <- update_active(rate, rand, active, status, label = "h")
nDepartures <- res[[1]]
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0
if (at == 2) dat$epi$dh.flow <- c(0, nDepartures) else dat$epi$dh.flow[at] <- nDepartures
if(dat$control$track.departure){
if (at == 2) dat$epi$hd.flow <- c(0, nDepartures) else dat$epi$hd.flow[at] <- nDepartures
}
# Recovered Departures --------------------------------------------------------
res <- update_active(rate, rand, active, status, label = "r")
nDepartures <- res[[1]]
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0
if (at == 2) dat$epi$dr.flow <- c(0, nDepartures) else dat$epi$dr.flow[at] <- nDepartures
if(dat$control$track.departure){
if (at == 2) dat$epi$rd.flow <- c(0, nDepartures) else dat$epi$rd.flow[at] <- nDepartures
}
# return --------------------------------------------------------
return(dat)
......
......@@ -56,14 +56,12 @@ infection.FUN <- function(dat, at, seed = NULL){
## Extract param
act.rate.i <- dat$param$act.rate.i
inf.prob.i <- dat$param$inf.prob.i
if(type %in% c("SEIQHR", "SEIQHRF")){
quar.rate <- dat$param$quar.rate
disch.rate <- dat$param$disch.rate
act.rate.e <- dat$param$act.rate.e
act.rate.q <- dat$param$act.rate.q
inf.prob.e <- dat$param$inf.prob.e
inf.prob.q <- dat$param$inf.prob.q
}
quar.rate <- dat$param$quar.rate
disch.rate <- dat$param$disch.rate
act.rate.e <- dat$param$act.rate.e
act.rate.q <- dat$param$act.rate.q
inf.prob.e <- dat$param$inf.prob.e
inf.prob.q <- dat$param$inf.prob.q
# Transmission from infected (a)
## Expected acts
......@@ -175,11 +173,12 @@ infection.FUN <- function(dat, at, seed = NULL){
}
## Output
tmp <- ifelse(type %in% c("SEIQHR", "SEIQHRF"), nExp.i + nExp.q, nExp.i)
if (at == 2) {
dat$epi$se.flow <- c(0, tmp)
} else {
dat$epi$se.flow[at] <- tmp
if(dat$control$track.flow){
if (at == 2) {
dat$epi$se.flow <- c(0, nExp.i + nExp.q)
} else {
dat$epi$se.flow[at] <- nExp.i + nExp.q
}
}
dat
......
......@@ -18,9 +18,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
if(!is.null(seed)) set.seed(seed)
# Conditions --------------------------------------------------------------------------
if (!(dat$control$type %in% c("SIR", "SIS", "SEIR", "SEIQHR", "SEIQHRF"))) {
return(dat)
}
if (dat$control$type != "SEIQHRF") return(dat)
# Variables ---------------------------------------------------------------------------
active <- dat$attr$active
......@@ -30,7 +28,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
type <- dat$control$type
expTime <- dat$attr$expTime
# ------------------------- progress from exposed to infectious --------------------------
# ------------------------- progress from exposed to infectious (E -> I)--------------------------
res <- update_status(rate = dat$param$prog.rate,
rand = dat$control$prog.rand,
active, dat$attr$status,
......@@ -49,7 +47,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
if (type %in% c("SEIQHR", "SEIQHRF")) {
# ----- quarantine -------
# ----- quarantine (I -> Q)-------
quar.rate <- dat$param$quar.rate
rate <- ifelse(length(quar.rate) > 1, quar.rate[at], quar.rate)
......@@ -67,7 +65,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
status <- res[[3]]
dat$attr$status <- status
# ----- need to be hospitalised -------
# ----- need to be hospitalised (I, Q -> H)-------
rate <- rep(dat$param$hosp.rate,
sum(dat$attr$active == 1 & (dat$attr$status == "i" | dat$attr$status == "q")))
......@@ -86,8 +84,8 @@ recovery.FUN <- function(dat, at, seed = NULL) {
status <- res[[3]]
dat$attr$status <- status
# ------------------------- discharge from need to be hospitalised ---------------------------
recovState <- ifelse(type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF"), "r", "s")
# ------------------------- discharge from need to be hospitalised (H -> R)---------------------------
recovState <- "r"
disch.rate <- dat$param$disch.rate
dcrate <- ifelse(length(disch.rate) > 1, disch.rate[at], disch.rate)
......@@ -108,7 +106,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
#cat("finished discharge \n")
}
# --------------------------------------- recover -------------------------------------------------
# --------------------------------------- recover (I, Q -> R)-------------------------------------------------
res <- update_status(rate = dat$param$rec.rate,
rand = dat$control$rec.rand,
active,
......@@ -125,7 +123,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
status <- res[[3]]
dat$attr$status <- status
# --------------------------------------- recover (E-> R)-------------------------------------------------
res <- update_status(rate = dat$param$arec.rate,
rand = dat$control$arec.rand,
active,
......@@ -137,13 +135,13 @@ recovery.FUN <- function(dat, at, seed = NULL) {
prog.dist.scale = dat$param$arec.dist.scale,
prog.dist.shape = dat$param$arec.dist.shape)
nRecov <- nRecov + res[[1]]
nRecov2 <- res[[1]]
dat$attr$recovTime[res[[2]]] <- at
status <- res[[3]]
dat$attr$status <- status
if (type %in% c("SEIQHRF")) {
# -------------------------------------- case fatality ----------------------------------------
# -------------------------------------- case fatality (H -> F)----------------------------------------
## obtain rates
idsElig <- which(active == 1 & status =="h")
if(length(idsElig) > 0){
......@@ -179,42 +177,24 @@ recovery.FUN <- function(dat, at, seed = NULL) {
#cat("finished fatality \n")
}
# Output ------------------------------------------------------------------
outName_a <- ifelse(type %in% c("SIR", "SEIR"), "ir.flow", "is.flow")
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
outName_b <- "ei.flow"
}
if (type %in% c("SEIQHR", "SEIQHRF")) {
outName_c <- "iq.flow"
outName_d <- "iq2h.flow"
}
if (type %in% c("SEIQHRF")) {
outName_e <- "hf.flow"
}
## Summary statistics
if (at == 2) {
dat$epi[[outName_a[1]]] <- c(0, nRecov)
if (type %in% c("SEIR", "SEIQHR")) {
dat$epi[[outName_b[1]]] <- c(0, nProg)
}
if (type %in% c("SEIQHR", "SEIQHRF")) {
dat$epi[[outName_c[1]]] <- c(0, nQuar)
dat$epi[[outName_d[1]]] <- c(0, nHosp)
}
if (type %in% c("SEIQHRF")) {
dat$epi[[outName_e[1]]] <- c(0, nFat)
}
} else {
dat$epi[[outName_a[1]]][at] <- nRecov
if (type %in% c("SEIR", "SEIQHR")) {
dat$epi[[outName_b[1]]][at] <- nProg
}
if (type %in% c("SEIQHR", "SEIQHRF")) {
dat$epi[[outName_c[1]]][at] <- nQuar
dat$epi[[outName_d[1]]][at] <- nHosp
}
if (type %in% c("SEIQHRF")) {
dat$epi[[outName_e[1]]][at] <- nFat
# Output Summary statistics------------------------------------------------------------------
if(dat$control$track.flow){
if (at == 2) {
dat$epi[["iq2r.flow"[1]]] <- c(0, nRecov)
dat$epi[["er.flow"[1]]] <- c(0, nRecov2)
dat$epi[["ei.flow"[1]]] <- c(0, nProg)
dat$epi[["iq.flow"[1]]] <- c(0, nQuar)
dat$epi[["iq2h.flow"[1]]] <- c(0, nHosp)
dat$epi[["hf.flow"[1]]] <- c(0, nFat)
dat$epi[["hr.flow"[1]]] <- c(0, nDisch)
} else {
dat$epi[["iq2r.flow"[1]]][at] <- nRecov
dat$epi[["er.flow"[1]]][at] <- nRecov2
dat$epi[["ei.flow"[1]]][at] <- nProg
dat$epi[["iq.flow"[1]]][at] <- nQuar
dat$epi[["iq2h.flow"[1]]][at] <- nHosp
dat$epi[["hf.flow"[1]]][at] <- nFat
dat$epi[["hr.flow"[1]]][at] <- nDisch
}
}
......
......@@ -69,6 +69,10 @@
#' this to \code{FALSE} is recommended when running models with new modules
#' specified.
#' @param ncores Number of physical CPU cores used for parallel computation.
#' @param track.flows If TRUE, function \code{\link{seiqhrf}} returns the number of changes in each procedure transition (See model description in vignette).
#' @param track.arrival If TRUE, function \code{\link{seiqhrf}} returns the number of departures of all compartments
#' in every time step in each simulation.
#' @param track.departure Similar to track.departure.
#'
#' @section New Modules:
#' Base ICM models use a set of module functions that specify
......@@ -89,7 +93,7 @@ control_seiqhrf <- function(nsteps = 366, nsims = 8,
disch.rand = TRUE, rec.rand = FALSE, arec.rand = TRUE,
fat.rand = TRUE, a.rand = TRUE, d.rand = TRUE,
verbose = FALSE, verbose.int = 0, skip.check = FALSE,
ncores = 4) {
ncores = 4, track.flows = TRUE, track.departure = FALSE, track.arrival = FALSE) {
# Get arguments
p <- list()
......
......@@ -19,7 +19,10 @@ control_seiqhrf(
verbose = FALSE,
verbose.int = 0,
skip.check = FALSE,
ncores = 4
ncores = 4,
track.flows = TRUE,
track.departure = FALSE,
track.arrival = FALSE
)
}
\arguments{
......@@ -101,6 +104,13 @@ this to \code{FALSE} is recommended when running models with new modules
specified.}
\item{ncores}{Number of physical CPU cores used for parallel computation.}
\item{track.flows}{If TRUE, function \code{\link{seiqhrf}} returns the number of changes in each procedure transition (See model description in vignette).}
\item{track.departure}{Similar to track.departure.}
\item{track.arrival}{If TRUE, function \code{\link{seiqhrf}} returns the number of departures of all compartments
in every time step in each simulation.}
}
\value{
A \code{control.seiqhrf} object containing required control parameters and module functions.
......
......@@ -34,7 +34,7 @@
</button>
<span class="navbar-brand">
<a class="navbar-link" href="index.html">sirplus</a>
<span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.3.2</span>
<span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.3.3</span>
</span>
</div>
......
......@@ -74,7 +74,7 @@ model types with additional compartments." />
</button>
<span class="navbar-brand">
<a class="navbar-link" href="../index.html">sirplus</a>
<span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.3.2</span>
<span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.3.3</span>
</span>
</div>
......@@ -146,7 +146,10 @@ model types with additional compartments.</p>
<span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>verbose.int</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>skip.check</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>ncores</span> <span class='kw'>=</span> <span class='fl'>4</span>
<span class='kw'>ncores</span> <span class='kw'>=</span> <span class='fl'>4</span>,
<span class='kw'>track.flows</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>,
<span class='kw'>track.departure</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>track.arrival</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>
)</pre>
<h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a>Arguments</h2>
......@@ -261,6 +264,19 @@ specified.</p></td>
<th>ncores</th>
<td><p>Number of physical CPU cores used for parallel computation.</p></td>
</tr>
<tr>
<th>track.flows</th>
<td><p>If TRUE, function <code><a href='seiqhrf.html'>seiqhrf</a></code> returns the number of changes in each procedure transition (See model description in vignette).</p></td>
</tr>
<tr>
<th>track.departure</th>
<td><p>Similar to track.departure.</p></td>
</tr>
<tr>
<th>track.arrival</th>
<td><p>If TRUE, function <code><a href='seiqhrf.html'>seiqhrf</a></code> returns the number of departures of all compartments
in every time step in each simulation.</p></td>
</tr>
</table>
<h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2>
......
......@@ -17,7 +17,7 @@ test_that("Identical output as Churches' original function: arrivals.FUN", {
for(seed in seed_list){
dat1 <- do.call(arrivals.FUN, list(dat, at, seed))
dat2 <- do.call(arrivals.seiqhrf.icm, list(dat, at, seed))
comp[i] <- identical(dat1, dat2)
comp[i] <- identical(dat1[1:4], dat2[1:4])
i <- i + 1
}
......
......@@ -16,7 +16,7 @@ test_that("Identical output as Churches' original function: departures.FUN", {
for(seed in seed_list){
dat1 <- do.call(departures.FUN, list(dat, at, seed))
dat2 <- do.call(departures.seiqhrf.icm, list(dat, at, seed))
comp[i] <- identical(dat1, dat2)
comp[i] <- identical(dat1[1:4], dat2[1:4])
i <- i + 1
}
......
......@@ -11,7 +11,7 @@ test_that("Identical output as Churches' original function: infection.FUN", {
dat1 <- do.call(infection.seiqhrf.icm, list(dat, at = 2, seed))
dat2 <- do.call(infection.FUN, list(dat, at = 2, seed))
comp[seed] <- identical(dat1, dat2)
comp[seed] <- identical(dat1[1:4], dat2[1:4])
}
......
......@@ -15,7 +15,7 @@ test_that("Identical output as Churches' original function: recovery.FUN", {
for(seed in seed_list){
dat1 <- do.call(recovery.FUN, list(dat, at, seed))
dat2 <- do.call(progress.seiqhrf.icm, list(dat, at, seed))
comp[i] <- identical(dat1, dat2)
comp[i] <- identical(dat1[1:4], dat2[1:4])
i <- i + 1
}
......
......@@ -9,7 +9,7 @@ test_that("seiqhrf and simulate_seiqhrf produce identical output", {
rec.rand = TRUE ### has to be fixed TRUE for comparison
No_seeds <- 10
seed_list <- sample(1:1000, No_seeds)
seed_list <- sample(1:100, No_seeds)
comp <- rep(NA, No_seeds)
i <- 1
for(seed in seed_list){
......@@ -26,7 +26,7 @@ test_that("seiqhrf and simulate_seiqhrf produce identical output", {
control <- control_seiqhrf(nsteps = nsteps, nsims = nsims, rec.rand = rec.rand)
sirplus_res <- seiqhrf(init, control, param, NULL, seed)
comp[i] <- identical(Churhes_res[3:4], sirplus_res[3:4]) # Due to $usr.specified in control and param, can only compare "epi" and "times"
comp[i] <- identical(Churhes_res[4], sirplus_res[4]) # Due to $usr.specified in control and param, 0.3.3 changed $epi, so can only compare times"
i <- i + 1
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment