From 78e69880302d1a0080724cab7746c58b5240eacd Mon Sep 17 00:00:00 2001
From: pqiao29 <pqiao@student.unimelb.edu.au>
Date: Wed, 1 Apr 2020 15:17:07 +1100
Subject: [PATCH] finish clean up arrivals.FUN
---
R/Chu_arrivals.R | 230 +++++++++++++++++++++++++++++++
R/FUN_arrivals.R | 72 ++++++++++
R/crosscheck.seiqhrf.R | 4 +-
R/mod_vital.R | 241 ---------------------------------
man/arrivals.FUN.Rd | 8 +-
man/departures.FUN.Rd | 2 +-
man/infection.FUN.Rd | 2 +-
man/recovery.FUN.Rd | 2 +-
tests/testthat/test-arrivals.R | 26 ++++
9 files changed, 338 insertions(+), 249 deletions(-)
create mode 100644 R/Chu_arrivals.R
create mode 100644 R/FUN_arrivals.R
delete mode 100644 R/mod_vital.R
create mode 100644 tests/testthat/test-arrivals.R
diff --git a/R/Chu_arrivals.R b/R/Chu_arrivals.R
new file mode 100644
index 0000000..db04277
--- /dev/null
+++ b/R/Chu_arrivals.R
@@ -0,0 +1,230 @@
+arrivals.seiqhrf.icm <- function(dat, at, seed = NULL) {
+
+ if(!is.null(seed)) set.seed(seed)
+
+ # Conditions --------------------------------------------------------------
+ if (dat$param$vital == FALSE) {
+ return(dat)
+ }
+
+ # Variables ---------------------------------------------------------------
+ a.rate <- dat$param$a.rate
+ a.rate.g2 <- dat$param$a.rate.g2
+ a.rand <- dat$control$a.rand
+ groups <- dat$param$groups
+ nOld <- dat$epi$num[at - 1]
+
+ # checking params, should be in control.icm or params.icm eventually
+ type <- dat$control$type
+ nsteps <- dat$control$nsteps
+
+ if (!(length(a.rate) == 1 || length(a.rate == nsteps))) {
+ stop("Length of a.rate must be 1 or the value of nsteps")
+ }
+ if (!is.null(a.rate.g2) &&
+ !(length(a.rate.g2) == 1 || length(a.rate.g2 == nsteps))) {
+ stop("Length of a.rate.g2 must be 1 or the value of nsteps")
+ }
+
+ a.prop.e <- dat$param$a.prop.e
+ if (!(length(a.prop.e) == 1 || length(a.prop.e == nsteps))) {
+ stop("Length of a.prop.e must be 1 or the value of nsteps")
+ }
+ a.prop.i <- dat$param$a.prop.i
+ if (!(length(a.prop.i) == 1 || length(a.prop.i == nsteps))) {
+ stop("Length of a.prop.i must be 1 or the value of nsteps")
+ }
+ a.prop.q <- dat$param$a.prop.q
+ if (!(length(a.prop.q) == 1 || length(a.prop.q == nsteps))) {
+ stop("Length of a.prop.q must be 1 or the value of nsteps")
+ }
+
+ a.prop.e.g2 <- dat$param$a.prop.e.g2
+ if (!is.null(a.prop.e.g2) &&
+ !(length(a.prop.e.g2) == 1 || length(a.prop.e.g2 == nsteps))) {
+ stop("Length of a.prop.e.g2 must be 1 or the value of nsteps")
+ }
+ a.prop.i.g2 <- dat$param$a.prop.i.g2
+ if (!is.null(a.prop.i.g2) &&
+ !(length(a.prop.i.g2) == 1 || length(a.prop.i.g2 == nsteps))) {
+ stop("Length of a.prop.i.g2 must be 1 or the value of nsteps")
+ }
+ a.prop.q.g2 <- dat$param$a.prop.q.g2
+ if (!is.null(a.prop.q.g2) &&
+ !(length(a.prop.q.g2) == 1 || length(a.prop.q.g2 == nsteps))) {
+ stop("Length of a.prop.q.g2 must be 1 or the value of nsteps")
+ }
+
+ # Process -----------------------------------------------------------------
+ nArrivals <- nArrivals.g2 <- 0
+
+ if (groups == 1) {
+ if (a.rand == TRUE) {
+ nArrivals <- sum(rbinom(nOld, 1, a.rate))
+ }
+ if (a.rand == FALSE) {
+ nArrivals <- round(nOld * a.rate)
+ }
+ }
+ if (groups == 2) {
+ nOldG2 <- dat$epi$num.g2[at - 1]
+ if (a.rand == TRUE) {
+ if (is.na(a.rate.g2)) {
+ nArrivals <- sum(rbinom(nOld, 1, a.rate))
+ nArrivals.g2 <- sum(rbinom(nOld, 1, a.rate))
+ } else {
+ nArrivals <- sum(rbinom(nOld, 1, a.rate))
+ nArrivals.g2 <- sum(rbinom(nOldG2, 1, a.rate.g2))
+ }
+ }
+ if (a.rand == FALSE) {
+ if (is.na(a.rate.g2)) {
+ nArrivals <- round(nOld * a.rate)
+ nArrivals.g2 <- round(nOld * a.rate)
+ } else {
+ nArrivals <- round(nOld * a.rate)
+ nArrivals.g2 <- round(nOldG2 * a.rate.g2)
+ }
+ }
+ }
+
+
+ ## Set attributes
+ totArrivals <- 0
+ totArrivals.g2 <- 0
+
+ # partition arrivals into compartments
+ if (length(a.prop.e) > 1) {
+ nArrivals.e <- round(nArrivals*(a.prop.e[at]))
+ totArrivals <- totArrivals + nArrivals.e
+ if (!is.null(a.prop.e.g2)) {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ } else {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ }
+ } else {
+ nArrivals.e <- round(nArrivals*a.prop.e)
+ totArrivals <- totArrivals + nArrivals.e
+ if (!is.null(a.prop.e.g2)) {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ } else {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ }
+ }
+
+ if (length(a.prop.i) > 1) {
+ nArrivals.i <- round(nArrivals*(a.prop.i[at]))
+ totArrivals <- totArrivals + nArrivals.i
+ if (!is.null(a.prop.i.g2)) {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ } else {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ }
+ } else {
+ nArrivals.i <- round(nArrivals*a.prop.i)
+ totArrivals <- totArrivals + nArrivals.i
+ if (!is.null(a.prop.i.g2)) {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ } else {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ }
+ }
+
+ if (length(a.prop.q) > 1) {
+ nArrivals.q <- round(nArrivals*(a.prop.q[at]))
+ totArrivals <- totArrivals + nArrivals.q
+ if (!is.null(a.prop.q.g2)) {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ } else {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ }
+ } else {
+ nArrivals.q <- round(nArrivals*a.prop.q)
+ totArrivals <- totArrivals + nArrivals.q
+ if (!is.null(a.prop.q.g2)) {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ } else {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ }
+ }
+
+ # debug
+ # print("totArrivals:")
+ # print(totArrivals)
+ # print("totArrivals.g2:")
+ # print(totArrivals.g2)
+ # print("----")
+ #
+ # group 1
+ 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),
+ rep("s", totArrivals - nArrivals.e - nArrivals.i - 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))
+ dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals))
+ dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals))
+ dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals))
+
+ # group 2
+ if (length(totArrivals.g2) > 0) {
+ dat$attr$active <- c(dat$attr$active, rep(1, totArrivals.g2))
+ dat$attr$group <- c(dat$attr$group, rep(2, totArrivals.g2))
+ dat$attr$status <- c(dat$attr$status,
+ rep("e", nArrivals.e.g2),
+ rep("i", nArrivals.i.g2),
+ rep("q", nArrivals.q.g2),
+ rep("s", totArrivals.g2 - nArrivals.e.g2 -
+ nArrivals.i.g2 - nArrivals.q.g2))
+ dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals.g2))
+ dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals.g2))
+ dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals.g2))
+ dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals.g2))
+ dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals.g2))
+ dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals.g2))
+ }
+
+ # 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 (length(totArrivals.g2) > 0 && dat$param$groups == 2) {
+ if (at == 2) {
+ dat$epi$a.flow.g2 <- c(0, totArrivals.g2)
+ dat$epi$a.e.flow.g2 <- c(0, nArrivals.e.g2)
+ dat$epi$a.i.flow.g2 <- c(0, nArrivals.i.g2)
+ dat$epi$a.q.flow.g2 <- c(0, nArrivals.q.g2)
+ } else {
+ dat$epi$a.flow.g2[at] <- totArrivals.g2
+ dat$epi$a.e.flow.g2[at] <- c(0, nArrivals.e.g2)
+ dat$epi$a.i.flow.g2[at] <- c(0, nArrivals.i.g2)
+ dat$epi$a.q.flow.g2[at] <- c(0, nArrivals.q.g2)
+ }
+ }
+
+ return(dat)
+}
\ No newline at end of file
diff --git a/R/FUN_arrivals.R b/R/FUN_arrivals.R
new file mode 100644
index 0000000..5488e3e
--- /dev/null
+++ b/R/FUN_arrivals.R
@@ -0,0 +1,72 @@
+#' Arrivals function
+#'
+#' Function to handel background demographics for the SEIQHRF model.
+#' Specifically arrivals (births and immigration). Uses the original EpiModel
+#' code currently. A replacement that implements modelling the arrival of
+#' infected individuals is under development -- but for now, all arrivals
+#' go into the S compartment..
+#'
+#' @param dat Input data needed.
+#' @param at time point
+#' @param seed random seed for checking consistency
+#'
+#' @return Updated dat
+#' @export
+arrivals.FUN <- function(dat, at, seed = NULL) {
+
+ if(!is.null(seed)) set.seed(seed)
+ # 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
+ a.rate.g2 <- dat$param$a.rate.g2
+ a.prop.e.g2 <- dat$param$a.prop.e.g2
+ a.prop.i.g2 <- dat$param$a.prop.i.g2
+ a.prop.q.g2 <- dat$param$a.prop.q.g2
+
+
+ # 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$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))
+ dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals))
+ 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)
+ 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)
+}
\ No newline at end of file
diff --git a/R/crosscheck.seiqhrf.R b/R/crosscheck.seiqhrf.R
index 2101c43..5f4a726 100644
--- a/R/crosscheck.seiqhrf.R
+++ b/R/crosscheck.seiqhrf.R
@@ -65,8 +65,8 @@ crosscheck.seiqhrf.icm <- function(param, init, control) {
}
}
- ## Check param (modified from infection.FUN) -------------------------------------
- check_list <- c("act.rate.i", "inf.prob.i")
+ ## Check param (modified from infection.FUN and arrivals.FUN) -------------------------------------
+ check_list <- c("act.rate.i", "inf.prob.i", "a.rate", "a.prop.e", "a.prop.i", "a.prop.q")
if(control$type %in% c("SEIQHR", "SEIQHRF")) check_list<- c(check_list, "quar.rate", "disch.rate", "act.rate.e", "inf.prob.e", "act.rate.q", "inf.prob.q")
check_idx <- which(names(param) %in% check_list)
diff --git a/R/mod_vital.R b/R/mod_vital.R
deleted file mode 100644
index 43d7f7b..0000000
--- a/R/mod_vital.R
+++ /dev/null
@@ -1,241 +0,0 @@
-#' Arrivals function
-#'
-#' Function to handel background demographics for the SEIQHRF model.
-#' Specifically arrivals (births and immigration). Uses the original EpiModel
-#' code currently. A replacement that implements modelling the arrival of
-#' infected individuals is under development -- but for now, all arrivals
-#' go into the S compartment..
-#'
-#' @param dat Input data needed.
-#' @param at Not sure???
-#'
-#' @return Updated dat
-#' @export
-arrivals.FUN <- function(dat, at) {
-
- # Conditions --------------------------------------------------------------
- if (dat$param$vital == FALSE) {
- return(dat)
- }
-
- # Variables ---------------------------------------------------------------
- a.rate <- dat$param$a.rate
- a.rate.g2 <- dat$param$a.rate.g2
- a.rand <- dat$control$a.rand
- groups <- dat$param$groups
- nOld <- dat$epi$num[at - 1]
-
- # checking params, should be in control.icm or params.icm eventually
- type <- dat$control$type
- nsteps <- dat$control$nsteps
-
- if (!(length(a.rate) == 1 || length(a.rate == nsteps))) {
- stop("Length of a.rate must be 1 or the value of nsteps")
- }
- if (!is.null(a.rate.g2) &&
- !(length(a.rate.g2) == 1 || length(a.rate.g2 == nsteps))) {
- stop("Length of a.rate.g2 must be 1 or the value of nsteps")
- }
-
- a.prop.e <- dat$param$a.prop.e
- if (!(length(a.prop.e) == 1 || length(a.prop.e == nsteps))) {
- stop("Length of a.prop.e must be 1 or the value of nsteps")
- }
- a.prop.i <- dat$param$a.prop.i
- if (!(length(a.prop.i) == 1 || length(a.prop.i == nsteps))) {
- stop("Length of a.prop.i must be 1 or the value of nsteps")
- }
- a.prop.q <- dat$param$a.prop.q
- if (!(length(a.prop.q) == 1 || length(a.prop.q == nsteps))) {
- stop("Length of a.prop.q must be 1 or the value of nsteps")
- }
-
- a.prop.e.g2 <- dat$param$a.prop.e.g2
- if (!is.null(a.prop.e.g2) &&
- !(length(a.prop.e.g2) == 1 || length(a.prop.e.g2 == nsteps))) {
- stop("Length of a.prop.e.g2 must be 1 or the value of nsteps")
- }
- a.prop.i.g2 <- dat$param$a.prop.i.g2
- if (!is.null(a.prop.i.g2) &&
- !(length(a.prop.i.g2) == 1 || length(a.prop.i.g2 == nsteps))) {
- stop("Length of a.prop.i.g2 must be 1 or the value of nsteps")
- }
- a.prop.q.g2 <- dat$param$a.prop.q.g2
- if (!is.null(a.prop.q.g2) &&
- !(length(a.prop.q.g2) == 1 || length(a.prop.q.g2 == nsteps))) {
- stop("Length of a.prop.q.g2 must be 1 or the value of nsteps")
- }
-
- # Process -----------------------------------------------------------------
- nArrivals <- nArrivals.g2 <- 0
-
- if (groups == 1) {
- if (a.rand == TRUE) {
- nArrivals <- sum(rbinom(nOld, 1, a.rate))
- }
- if (a.rand == FALSE) {
- nArrivals <- round(nOld * a.rate)
- }
- }
- if (groups == 2) {
- nOldG2 <- dat$epi$num.g2[at - 1]
- if (a.rand == TRUE) {
- if (is.na(a.rate.g2)) {
- nArrivals <- sum(rbinom(nOld, 1, a.rate))
- nArrivals.g2 <- sum(rbinom(nOld, 1, a.rate))
- } else {
- nArrivals <- sum(rbinom(nOld, 1, a.rate))
- nArrivals.g2 <- sum(rbinom(nOldG2, 1, a.rate.g2))
- }
- }
- if (a.rand == FALSE) {
- if (is.na(a.rate.g2)) {
- nArrivals <- round(nOld * a.rate)
- nArrivals.g2 <- round(nOld * a.rate)
- } else {
- nArrivals <- round(nOld * a.rate)
- nArrivals.g2 <- round(nOldG2 * a.rate.g2)
- }
- }
- }
-
-
- ## Set attributes
- totArrivals <- 0
- totArrivals.g2 <- 0
-
- # partition arrivals into compartments
- if (length(a.prop.e) > 1) {
- nArrivals.e <- round(nArrivals*(a.prop.e[at]))
- totArrivals <- totArrivals + nArrivals.e
- if (!is.null(a.prop.e.g2)) {
- nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at]))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
- } else {
- nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at]))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
- }
- } else {
- nArrivals.e <- round(nArrivals*a.prop.e)
- totArrivals <- totArrivals + nArrivals.e
- if (!is.null(a.prop.e.g2)) {
- nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
- } else {
- nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
- }
- }
-
- if (length(a.prop.i) > 1) {
- nArrivals.i <- round(nArrivals*(a.prop.i[at]))
- totArrivals <- totArrivals + nArrivals.i
- if (!is.null(a.prop.i.g2)) {
- nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at]))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
- } else {
- nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at]))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
- }
- } else {
- nArrivals.i <- round(nArrivals*a.prop.i)
- totArrivals <- totArrivals + nArrivals.i
- if (!is.null(a.prop.i.g2)) {
- nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
- } else {
- nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
- }
- }
-
- if (length(a.prop.q) > 1) {
- nArrivals.q <- round(nArrivals*(a.prop.q[at]))
- totArrivals <- totArrivals + nArrivals.q
- if (!is.null(a.prop.q.g2)) {
- nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at]))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
- } else {
- nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at]))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
- }
- } else {
- nArrivals.q <- round(nArrivals*a.prop.q)
- totArrivals <- totArrivals + nArrivals.q
- if (!is.null(a.prop.q.g2)) {
- nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
- } else {
- nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2))
- totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
- }
- }
-
- # debug
- print("totArrivals:")
- print(totArrivals)
- print("totArrivals.g2:")
- print(totArrivals.g2)
- print("----")
-
- # group 1
- 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),
- rep("s", totArrivals - nArrivals.e - nArrivals.i - 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))
- dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals))
- dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals))
- dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals))
-
- # group 2
- if (length(totArrivals.g2) > 0) {
- dat$attr$active <- c(dat$attr$active, rep(1, totArrivals.g2))
- dat$attr$group <- c(dat$attr$group, rep(2, totArrivals.g2))
- dat$attr$status <- c(dat$attr$status,
- rep("e", nArrivals.e.g2),
- rep("i", nArrivals.i.g2),
- rep("q", nArrivals.q.g2),
- rep("s", totArrivals.g2 - nArrivals.e.g2 -
- nArrivals.i.g2 - nArrivals.q.g2))
- dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals.g2))
- dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals.g2))
- dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals.g2))
- dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals.g2))
- dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals.g2))
- dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals.g2))
- }
-
- # 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 (length(totArrivals.g2) > 0 && dat$param$groups == 2) {
- if (at == 2) {
- dat$epi$a.flow.g2 <- c(0, totArrivals.g2)
- dat$epi$a.e.flow.g2 <- c(0, nArrivals.e.g2)
- dat$epi$a.i.flow.g2 <- c(0, nArrivals.i.g2)
- dat$epi$a.q.flow.g2 <- c(0, nArrivals.q.g2)
- } else {
- dat$epi$a.flow.g2[at] <- totArrivals.g2
- dat$epi$a.e.flow.g2[at] <- c(0, nArrivals.e.g2)
- dat$epi$a.i.flow.g2[at] <- c(0, nArrivals.i.g2)
- dat$epi$a.q.flow.g2[at] <- c(0, nArrivals.q.g2)
- }
- }
-
- return(dat)
-}
diff --git a/man/arrivals.FUN.Rd b/man/arrivals.FUN.Rd
index 05bc137..c65ecb9 100644
--- a/man/arrivals.FUN.Rd
+++ b/man/arrivals.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/FUN_arrivals.R
\name{arrivals.FUN}
\alias{arrivals.FUN}
\title{Arrivals function}
\usage{
-arrivals.FUN(dat, at)
+arrivals.FUN(dat, at, seed = NULL)
}
\arguments{
\item{dat}{Input data needed.}
-\item{at}{Not sure???}
+\item{at}{time point}
+
+\item{seed}{random seed for checking consistency}
}
\value{
Updated dat
diff --git a/man/departures.FUN.Rd b/man/departures.FUN.Rd
index 2231f95..4a8e952 100644
--- a/man/departures.FUN.Rd
+++ b/man/departures.FUN.Rd
@@ -1,5 +1,5 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/departure.FUN.R
+% Please edit documentation in R/FUN_departures.R
\name{departures.FUN}
\alias{departures.FUN}
\title{Departures function}
diff --git a/man/infection.FUN.Rd b/man/infection.FUN.Rd
index 9be730e..22fdd15 100644
--- a/man/infection.FUN.Rd
+++ b/man/infection.FUN.Rd
@@ -1,5 +1,5 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/infection.FUN.R
+% Please edit documentation in R/FUN_infection.R
\name{infection.FUN}
\alias{infection.FUN}
\title{Infection function}
diff --git a/man/recovery.FUN.Rd b/man/recovery.FUN.Rd
index 47e6150..7c5e74e 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/recovery.FUN.R
+% Please edit documentation in R/FUN_recovery.R
\name{recovery.FUN}
\alias{recovery.FUN}
\title{Progress icm}
diff --git a/tests/testthat/test-arrivals.R b/tests/testthat/test-arrivals.R
new file mode 100644
index 0000000..523f7fb
--- /dev/null
+++ b/tests/testthat/test-arrivals.R
@@ -0,0 +1,26 @@
+test_that("Identical output as Churches' original function: arrivals.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))
+ dat <- do.call(departures.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(arrivals.FUN, list(dat, at, seed))
+ dat2 <- do.call(arrivals.seiqhrf.icm, list(dat, at, seed))
+ comp[i] <- identical(dat1, dat2)
+ i <- i + 1
+ rm(.Random.seed)
+ }
+
+ expect_equal(sum(comp), No_seeds)
+})
+
+
+
+
--
GitLab