From 7eb6ff398de2aef147b930a9ab5ab74abe171968 Mon Sep 17 00:00:00 2001
From: pqiao29 <pqiao@student.unimelb.edu.au>
Date: Wed, 1 Apr 2020 12:15:27 +1100
Subject: [PATCH] finish clean up departures.FUN
---
R/Chu_departures.R | 259 +++++++++++++++++++++++++++++
R/internal_update_status.R | 74 +++++++++
R/mod_vital.R | 270 -------------------------------
R/recovery.FUN.R | 52 ------
man/departures.FUN.Rd | 6 +-
man/recovery.FUN.Rd | 2 +-
tests/testthat/test-departures.R | 21 +++
7 files changed, 359 insertions(+), 325 deletions(-)
create mode 100644 R/Chu_departures.R
create mode 100644 R/internal_update_status.R
create mode 100644 tests/testthat/test-departures.R
diff --git a/R/Chu_departures.R b/R/Chu_departures.R
new file mode 100644
index 0000000..0e5acfc
--- /dev/null
+++ b/R/Chu_departures.R
@@ -0,0 +1,259 @@
+departures.seiqhrf.icm <- function(dat, at, seed = NULL) {
+
+ if(!is.null(seed)) set.seed(seed)
+
+ # Conditions --------------------------------------------------------------
+ if (dat$param$vital == FALSE) {
+ return(dat)
+ }
+
+
+ # Variables ---------------------------------------------------------------
+ groups <- dat$param$groups
+ group <- dat$attr$group
+
+
+ # Susceptible departures ------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "s")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$ds.rate, dat$param$ds.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$ds.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$ds.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$ds.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$ds.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ # Exposed Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "e")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$de.rate, dat$param$de.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$de.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$de.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$de.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$de.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+
+ # Infected Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "i")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$di.rate, dat$param$di.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$di.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$di.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$di.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$di.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ # Quarantined Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "q")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$dq.rate, dat$param$dq.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$dq.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$dq.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$dq.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$dq.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ # Hospitalised Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "h")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$dh.rate, dat$param$dh.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$dh.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$dh.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$dh.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$dh.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+
+ # Recovered Departures --------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "r")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$dr.rate, dat$param$dr.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$dr.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$dr.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$dr.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$dr.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ return(dat)
+}
\ No newline at end of file
diff --git a/R/internal_update_status.R b/R/internal_update_status.R
new file mode 100644
index 0000000..dfa4e69
--- /dev/null
+++ b/R/internal_update_status.R
@@ -0,0 +1,74 @@
+# Needed for recovery.FUN --------------------------------------------------------------
+update_status <- function(rate, rand, active, status, label, state, at, prog, expTime = NULL, prog.dist.scale = NULL, prog.dist.shape = NULL){
+
+ smp_sz <- 0
+ at_idx <- NULL
+
+ idsElig <- which(active == 1 & status %in% label)
+ nElig <- length(idsElig)
+
+ if (nElig > 0) {
+ gElig <- rep(1, nElig)
+
+ if (rand) {
+ vecProg <- which(stats::rbinom(nElig, 1, rate) == 1)
+ if (length(vecProg) > 0) {
+ at_idx <- idsElig[vecProg]
+ smp_sz <- length(at_idx)
+ status[at_idx] <- state
+ }
+ }else{
+ do_sample <- TRUE
+ if(!prog){
+ smp_sz <- min(round(sum(rates[gElig == 1])), sum(gElig == 1))
+ smp_prob <- NULL
+ }else{
+ vecTimeSinceExp <- at - expTime[idsElig]
+ vecTimeSinceExp[is.na(vecTimeSinceExp)] <- 0
+ gammaRatesElig <- stats::pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale)
+ smp_sz <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE))
+ smp_prob <- gammaRatesElig[gElig == 1]
+ if(smp_sz <= 0) do_sample <- FALSE
+ }
+
+ if(do_sample){
+ at_idx <- EpiModel::ssample(idsElig[gElig == 1], smp_sz, prob = smp_prob)
+ status[at_idx] <- state
+ }
+
+
+ }
+
+ }
+
+ list(smp_sz, at_idx, status)
+}
+
+# Needed for departures.FUN ------------------------------------------------------------
+update_active <- function(rate, rand, active, status, label){
+
+ smp_sz <- 0
+ act_idx <- NULL
+
+ idsElig <- which(active == 1 & status == label)
+ nElig <- length(idsElig)
+
+ if (nElig > 0) {
+
+ gElig <- rep(1, nElig)
+
+ if (rand) {
+ vec <- which(stats::rbinom(nElig, 1, rate) == 1)
+
+ if (length(vec) > 0) {
+ act_idx <- idsElig[vec]
+ smp_sz <- length(act_idx)
+ }
+ } else {
+ smp_sz <- min(round(sum(rate[gElig == 1])), sum(gElig == 1))
+ act_idx <- EpiModel::ssample(idsElig[gElig == 1], smp_sz)
+ }
+ }
+
+ return(list(smp_sz, act_idx))
+}
\ No newline at end of file
diff --git a/R/mod_vital.R b/R/mod_vital.R
index 7082e67..43d7f7b 100644
--- a/R/mod_vital.R
+++ b/R/mod_vital.R
@@ -1,273 +1,3 @@
-#' Departures function
-#'
-#' Function to handel background demographics for the SEIQHRF model.
-#' Specifically departures (deaths not due to the virus, and emigration).
-#'
-#' @param dat Merged input parameters.
-#' @param at Step number
-#'
-#' @return Updated dat
-#' @export
-departures.FUN <- function(dat, at) {
-
- # Conditions --------------------------------------------------------------
- if (dat$param$vital == FALSE) {
- return(dat)
- }
-
-
- # Variables ---------------------------------------------------------------
- groups <- dat$param$groups
- group <- dat$attr$group
-
-
- # Susceptible departures ------------------------------------------------------
- nDepartures <- nDeparturesG2 <- 0
- idsElig <- which(dat$attr$active == 1 & dat$attr$status == "s")
- nElig <- length(idsElig)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- c(dat$param$ds.rate, dat$param$ds.rate.g2)
- ratesElig <- rates[gElig]
-
- if (dat$control$d.rand == TRUE) {
- vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecDepartures) > 0) {
- idsDpt <- idsElig[vecDepartures]
- nDepartures <- sum(group[idsDpt] == 1)
- nDeparturesG2 <- sum(group[idsDpt] == 2)
- dat$attr$active[idsDpt] <- 0
- }
- } else {
- nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
- if (groups == 2) {
- nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
- dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
- }
- }
- }
-
- if (at == 2) {
- dat$epi$ds.flow <- c(0, nDepartures)
- if (groups == 2) {
- dat$epi$ds.flow.g2 <- c(0, nDeparturesG2)
- }
- } else {
- dat$epi$ds.flow[at] <- nDepartures
- if (groups == 2) {
- dat$epi$ds.flow.g2[at] <- nDeparturesG2
- }
- }
-
-# Exposed Departures ---------------------------------------------------------
- nDepartures <- nDeparturesG2 <- 0
- idsElig <- which(dat$attr$active == 1 & dat$attr$status == "e")
- nElig <- length(idsElig)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- c(dat$param$de.rate, dat$param$de.rate.g2)
- ratesElig <- rates[gElig]
-
- if (dat$control$d.rand == TRUE) {
- vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecDepartures) > 0) {
- idsDpt <- idsElig[vecDepartures]
- nDepartures <- sum(group[idsDpt] == 1)
- nDeparturesG2 <- sum(group[idsDpt] == 2)
- dat$attr$active[idsDpt] <- 0
- }
- } else {
- nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
- if (groups == 2) {
- nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
- dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
- }
- }
- }
-
- if (at == 2) {
- dat$epi$de.flow <- c(0, nDepartures)
- if (groups == 2) {
- dat$epi$de.flow.g2 <- c(0, nDeparturesG2)
- }
- } else {
- dat$epi$de.flow[at] <- nDepartures
- if (groups == 2) {
- dat$epi$de.flow.g2[at] <- nDeparturesG2
- }
- }
-
-
- # Infected Departures ---------------------------------------------------------
- nDepartures <- nDeparturesG2 <- 0
- idsElig <- which(dat$attr$active == 1 & dat$attr$status == "i")
- nElig <- length(idsElig)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- c(dat$param$di.rate, dat$param$di.rate.g2)
- ratesElig <- rates[gElig]
-
- if (dat$control$d.rand == TRUE) {
- vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecDepartures) > 0) {
- idsDpt <- idsElig[vecDepartures]
- nDepartures <- sum(group[idsDpt] == 1)
- nDeparturesG2 <- sum(group[idsDpt] == 2)
- dat$attr$active[idsDpt] <- 0
- }
- } else {
- nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
- if (groups == 2) {
- nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
- dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
- }
- }
- }
-
- if (at == 2) {
- dat$epi$di.flow <- c(0, nDepartures)
- if (groups == 2) {
- dat$epi$di.flow.g2 <- c(0, nDeparturesG2)
- }
- } else {
- dat$epi$di.flow[at] <- nDepartures
- if (groups == 2) {
- dat$epi$di.flow.g2[at] <- nDeparturesG2
- }
- }
-
- # Quarantined Departures ---------------------------------------------------------
- nDepartures <- nDeparturesG2 <- 0
- idsElig <- which(dat$attr$active == 1 & dat$attr$status == "q")
- nElig <- length(idsElig)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- c(dat$param$dq.rate, dat$param$dq.rate.g2)
- ratesElig <- rates[gElig]
-
- if (dat$control$d.rand == TRUE) {
- vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecDepartures) > 0) {
- idsDpt <- idsElig[vecDepartures]
- nDepartures <- sum(group[idsDpt] == 1)
- nDeparturesG2 <- sum(group[idsDpt] == 2)
- dat$attr$active[idsDpt] <- 0
- }
- } else {
- nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
- if (groups == 2) {
- nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
- dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
- }
- }
- }
-
- if (at == 2) {
- dat$epi$dq.flow <- c(0, nDepartures)
- if (groups == 2) {
- dat$epi$dq.flow.g2 <- c(0, nDeparturesG2)
- }
- } else {
- dat$epi$dq.flow[at] <- nDepartures
- if (groups == 2) {
- dat$epi$dq.flow.g2[at] <- nDeparturesG2
- }
- }
-
- # Hospitalised Departures ---------------------------------------------------------
- nDepartures <- nDeparturesG2 <- 0
- idsElig <- which(dat$attr$active == 1 & dat$attr$status == "h")
- nElig <- length(idsElig)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- c(dat$param$dh.rate, dat$param$dh.rate.g2)
- ratesElig <- rates[gElig]
-
- if (dat$control$d.rand == TRUE) {
- vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecDepartures) > 0) {
- idsDpt <- idsElig[vecDepartures]
- nDepartures <- sum(group[idsDpt] == 1)
- nDeparturesG2 <- sum(group[idsDpt] == 2)
- dat$attr$active[idsDpt] <- 0
- }
- } else {
- nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
- if (groups == 2) {
- nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
- dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
- }
- }
- }
-
- if (at == 2) {
- dat$epi$dh.flow <- c(0, nDepartures)
- if (groups == 2) {
- dat$epi$dh.flow.g2 <- c(0, nDeparturesG2)
- }
- } else {
- dat$epi$dh.flow[at] <- nDepartures
- if (groups == 2) {
- dat$epi$dh.flow.g2[at] <- nDeparturesG2
- }
- }
-
-
- # Recovered Departures --------------------------------------------------------
- nDepartures <- nDeparturesG2 <- 0
- idsElig <- which(dat$attr$active == 1 & dat$attr$status == "r")
- nElig <- length(idsElig)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- c(dat$param$dr.rate, dat$param$dr.rate.g2)
- ratesElig <- rates[gElig]
-
- if (dat$control$d.rand == TRUE) {
- vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecDepartures) > 0) {
- idsDpt <- idsElig[vecDepartures]
- nDepartures <- sum(group[idsDpt] == 1)
- nDeparturesG2 <- sum(group[idsDpt] == 2)
- dat$attr$active[idsDpt] <- 0
- }
- } else {
- nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
- if (groups == 2) {
- nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
- dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
- }
- }
- }
-
- if (at == 2) {
- dat$epi$dr.flow <- c(0, nDepartures)
- if (groups == 2) {
- dat$epi$dr.flow.g2 <- c(0, nDeparturesG2)
- }
- } else {
- dat$epi$dr.flow[at] <- nDepartures
- if (groups == 2) {
- dat$epi$dr.flow.g2[at] <- nDeparturesG2
- }
- }
-
- return(dat)
-}
-
-
-
#' Arrivals function
#'
#' Function to handel background demographics for the SEIQHRF model.
diff --git a/R/recovery.FUN.R b/R/recovery.FUN.R
index 37af740..7bf3df6 100644
--- a/R/recovery.FUN.R
+++ b/R/recovery.FUN.R
@@ -1,50 +1,3 @@
-# internal function --------------------------------------------------------------
-update_status <- function(rate, rand, active, status, label, state, at, prog, expTime = NULL, prog.dist.scale = NULL, prog.dist.shape = NULL){
-
- smp_sz <- 0
- at_idx <- NULL
-
- idsElig <- which(active == 1 & status %in% label)
- nElig <- length(idsElig)
-
- if (nElig > 0) {
- gElig <- rep(1, nElig)
-
- if (rand) {
- vecProg <- which(stats::rbinom(nElig, 1, rate) == 1)
- if (length(vecProg) > 0) {
- at_idx <- idsElig[vecProg]
- smp_sz <- length(at_idx)
- status[at_idx] <- state
- }
- }else{
- do_sample <- TRUE
- if(!prog){
- smp_sz <- min(round(sum(rates[gElig == 1])), sum(gElig == 1))
- smp_prob <- NULL
- }else{
- vecTimeSinceExp <- at - expTime[idsElig]
- vecTimeSinceExp[is.na(vecTimeSinceExp)] <- 0
- gammaRatesElig <- stats::pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale)
- smp_sz <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE))
- smp_prob <- gammaRatesElig[gElig == 1]
- if(smp_sz <= 0) do_sample <- FALSE
- }
-
- if(do_sample){
- at_idx <- EpiModel::ssample(idsElig[gElig == 1], smp_sz, prob = smp_prob)
- status[at_idx] <- state
- }
-
-
- }
-
- }
-
- list(smp_sz, at_idx, status)
-}
-
-
#' Progress icm
#'
#' Function to get progress of icms
@@ -73,14 +26,9 @@ recovery.FUN <- function(dat, at, seed = NULL) {
# Variables ---------------------------------------------------------------------------
active <- dat$attr$active
status <- dat$attr$status
-
groups <- dat$param$groups
group <- dat$attr$group
-
type <- dat$control$type
- recovState <- ifelse(type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF"), "r", "s")
- hospState <- "h"
- fatState <- "f"
# ------------------------- progress from exposed to infectious --------------------------
res <- update_status(rate = dat$param$prog.rate,
diff --git a/man/departures.FUN.Rd b/man/departures.FUN.Rd
index 6175e9f..2231f95 100644
--- a/man/departures.FUN.Rd
+++ b/man/departures.FUN.Rd
@@ -1,15 +1,17 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mod_vital.R
+% Please edit documentation in R/departure.FUN.R
\name{departures.FUN}
\alias{departures.FUN}
\title{Departures function}
\usage{
-departures.FUN(dat, at)
+departures.FUN(dat, at, seed = NULL)
}
\arguments{
\item{dat}{Merged input parameters.}
\item{at}{Step number}
+
+\item{seed}{random seed for checking consistency}
}
\value{
Updated dat
diff --git a/man/recovery.FUN.Rd b/man/recovery.FUN.Rd
index 0f81dc3..47e6150 100644
--- a/man/recovery.FUN.Rd
+++ b/man/recovery.FUN.Rd
@@ -1,5 +1,5 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mod_status.R
+% Please edit documentation in R/recovery.FUN.R
\name{recovery.FUN}
\alias{recovery.FUN}
\title{Progress icm}
diff --git a/tests/testthat/test-departures.R b/tests/testthat/test-departures.R
new file mode 100644
index 0000000..6aca443
--- /dev/null
+++ b/tests/testthat/test-departures.R
@@ -0,0 +1,21 @@
+test_that("Identical output as Churches' original function: departures.FUN", {
+
+ at <- 2
+ dat <- do.call(initialize.FUN, set_param())
+ dat <- do.call(infection.FUN, list(dat, at))
+ dat <- do.call(recovery.FUN, list(dat, at))
+
+ No_seeds <- 10
+ seed_list <- sample(1:1000, No_seeds)
+ comp <- rep(NA, No_seeds)
+ i <- 1
+ 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)
+ i <- i + 1
+ rm(.Random.seed)
+ }
+
+ expect_equal(sum(comp), No_seeds)
+})
\ No newline at end of file
--
GitLab