Merge pull request #159 from dd-harp/dev

xde_setup_cohort bug fixes

R/adult-GeRM.R

@@ -67,7 +67,7 @@ dMYZdt.GeRM_ode <- function(t, y, pars, s) {
     with(pars$MYZpar[[s]],{
 
       Omega <- make_Omega(g, sigma, calK, nPatches)
-      Upsilon <- expm(-Omega*eip)
+      Upsilon <- expm::expm(-Omega*eip)
 
       dMdt <- Lambda - (Omega %*% M)
       dGdt <-  f*(M - G) - nu*G  - (Omega %*% G)
@@ -121,7 +121,7 @@ dMYZdt.GeRM_dde <- function(t, y, pars, s){
       dGdt <-  f*(M - G) - nu*G  - (Omega %*% G)
       dYdt <- f*q*kappa*(M - Y) - (Omega %*% Y)
       dZdt <- Upsilon %*% (fqkappa_eip * (M_eip - Y_eip)) - (Omega %*% Z)
-      dUdt <- as.vector(((1-dEIPdt(t,pars))*Omega_eip - Omega) %*% Upsilon)
+      dUdt <- as.vector(((1-dEIPdt(t,EIPmod))*Omega_eip - Omega) %*% Upsilon)
 
       return(c(dMdt, dGdt, dYdt, dZdt, dUdt, f*q*kappa, g, sigma))
     })

R/adult-RM.R

@@ -63,6 +63,9 @@ dMYZdt.RM_ode <- function(t, y, pars, s) {
     Z <- y[Z_ix]
 
     with(pars$MYZpar[[s]],{
+      Omega <- make_Omega(g, sigma, calK, nPatches)
+      Upsilon <- expm::expm(-Omega*eip)
+
       dM <- Lambda - (Omega %*% M)
       dP <- f*(M - P) - (Omega %*% P)
       dY <- f*q*kappa*(M - Y) - (Omega %*% Y)
@@ -350,7 +353,7 @@ make_parameters_MYZ_RM <- function(pars, g, sigma, f, q, nu, eggsPerBatch, eip,
 
   Omega   <- make_Omega(g, sigma, calK, pars$nPatches)
   MYZpar$Omega <- Omega
-  MYZpar$Upsilon <- expm(-Omega*eip)
+  MYZpar$Upsilon <- expm::expm(-Omega*eip)
   MYZpar$EIPmod <- setup_eip_static(eip=eip)
   MYZpar$eip <- eip
   MYZpar$calK <- calK

R/compute_terms.R

@@ -66,8 +66,8 @@ compute_terms.xde <- function(varslist, deout, pars, s, i) {
 compute_terms.cohort <- function(varslist, deout, pars, s, i) {
   time = deout[,1]
   d1 = length(time)
-
-  eir = matrix(pars$F_eir(time, pars), d1, pars$Hpar[[i]]$nStrata)
+  eir = as.matrix(with(pars$EIRpar, sapply(time, pars$F_eir, bday=bday, scale=scale)))
+  eir = shapeIt(eir, d1, pars$Hpar[[i]]$nStrata)
   ni = compute_NI(deout, pars, i)
   pr = F_pr(varslist, pars, i)
   return(list(time=time,eir=eir,pr=pr,ni=ni))

R/diffeqn.R

@@ -138,13 +138,13 @@ xDE_diffeqn_mosy <- function(t, y, pars) {
 #' @param a age of a cohort
 #' @param y state vector
 #' @param pars a [list]
-#' @param F_eir a trace function that returns the eir
+#' @param F_eir a trace function that returns the eir as a function of time
 #' @return a [list] containing the vector of all state derivatives
 #' @export
 xDE_diffeqn_cohort <- function(a, y, pars, F_eir) {
 
   # EIR: entomological inoculation rate trace
-  pars$EIR[[1]] <- F_eir(a, pars)*pars$BFpar$relativeBitingRate[[1]][[1]]
+  pars$EIR[[1]] <- with(pars$EIRpar, F_eir(a, bday, scale))*pars$BFpar$relativeBitingRate[[1]][[1]]
 
   # FoI: force of infection
   pars <- Exposure(a, y, pars)

R/plot-terms.R

@@ -101,7 +101,7 @@ xde_plot_PR = function(pars, i=1, clrs="black", llty=1, stable=FALSE, add_axes=T
               ylab = "Prevalence", xlab = "Time")
   }
 
-  xde_lines_PR(tm, vars$terms$PR[[i]], pars$Hpar[[i]]$nStrata, clrs, llty)
+  xde_lines_PR(tm, vars$terms$pr[[i]], pars$Hpar[[i]]$nStrata, clrs, llty)
 }
 
 #' Add lines for the prevalence / parasite rate (PR) from a model of human infection and immunity

R/setup.R

@@ -301,15 +301,16 @@ xde_setup_human = function(modelName = "unnamed",
 
 #' @title Set up a model for xde_diffeqn_cohort
 #' @param F_eir is a function F_eir(t, pars) that returns the daily FoI
+#' @param bday the birthday of a cohort
+#' @param scale the birthday of a cohort
 #' @param modelName is a name for the model (arbitrary)
 #' @param Xname is a character string defining a X model
 #' @param HPop is the number of humans in each patch
 #' @param searchB is a vector of search weights for blood feeding
-#' @param residence is a vector that describes the patch where each human stratum lives
 #' @param Xopts a list to configure the X model
 #' @return a [list]
 #' @export
-xde_setup_cohort = function(F_eir,
+xde_setup_cohort = function(F_eir, bday=0, scale=1,
                            modelName = "unnamed",
 
                            # Dynamical Components
@@ -318,7 +319,6 @@ xde_setup_cohort = function(F_eir,
                            # Model Structure
                            HPop=1000,
                            searchB = 1,
-                           residence = 1,
 
                            # Human Strata / Options
                            Xopts = list()
@@ -329,18 +329,26 @@ xde_setup_cohort = function(F_eir,
   class(pars$xde) <- "cohort"
   class(pars$compute) = "cohort"
 
+  pars$nVectors = 1
+  pars$nHosts = 1
+  pars$nPatches = 1
+
   pars$modelName = modelName
   pars$Xname = Xname
 
   pars$F_eir = F_eir
+  pars$EIRpar = list()
+  pars$EIRpar$bday = bday
+  pars$EIRpar$scale = scale
 
   # Structure
   nStrata = length(HPop)
-  pars$nPatches = as.integer(nStrata)
-  pars$nHosts = 1
+  residence = rep(1, nStrata)
 
   pars = setup_Hpar_static(pars, 1, HPop)
   pars = setup_BloodFeeding(pars, 1, 1, list(), residence, searchB, NULL)
+  pars$BFpar$TimeSpent[[1]] = make_TimeSpent_athome(1, residence)
+  pars = make_TaR(0, pars, 1, 1)
 
   # Dynamics
   pars = setup_Xpar(Xname, pars, 1, Xopts)

_pkgdown.yml

@@ -809,6 +809,10 @@ reference:
   - compute_terms.human
   - compute_terms.na
   - compute_terms_steady
+  - compute_fqZ
+  - compute_fqZ_ix
+  - compute_fqM
+  - compute_fqM_ix
   - compute_NI
   - compute_NI_ix
   - compute_vars_full

vignettes/human_hmoi.Rmd

@@ -92,7 +92,7 @@ m20 <- 1.5
 h <- r2*m20
 m10 <- h/r1
 
-EIR <- h/b
+EIR <- rep(h/b, 3) 
 ```
 
 ```{r}
@@ -101,12 +101,13 @@ params$nStrata = nStrata
 params$nVectors = 1 
 params$nHosts = 1
 params$nPatches = 1 
-params$eir = rep(EIR, 3)
 params$EIR = list()
 
-F_eir = function(t, pars){
-  pars$EIR[[1]] = pars$eir
-}
+fF_eir = function(EIR){
+  EIR = as.vector(EIR) 
+  return(function(t, bday=0, scale=1){EIR})
+} 
+F_eir = fF_eir(EIR)
 
 params = make_parameters_demography_null(pars = params, H=H) 
 
@@ -152,17 +153,7 @@ Hpop = c(100, 500, 250)
 ```
 
 ```{r}
-h/b
-```
-
-
-```{r}
-fF_eir1 = function(EIR){return(function(t, pars){EIR})} 
-F_eir1 = fF_eir1(EIR)
-```
-
-```{r}
-xde_setup_cohort(F_eir1, Xname="hMoI", HPop=Hpop, Xopts = Xo) ->test_hMoI
+xde_setup_cohort(F_eir, Xname="hMoI", HPop=Hpop, Xopts = Xo) ->test_hMoI
 ```
 
 ```{r}

vignettes/human_sip.Rmd

@@ -91,9 +91,11 @@ params$nHosts = 1
 params$eir = EIR
 params$EIR = list()
 
-F_eir = function(t, pars){
-  pars$EIR[[1]] = params$eir
-}
+fF_eir = function(EIR){
+  EIR = as.vector(EIR) 
+  return(function(t, bday=0, scale=1){EIR})
+} 
+F_eir = fF_eir(EIR)
 
 params = make_parameters_demography_null(pars = params, H=H) 
 
@@ -138,15 +140,9 @@ Xo = list(b=0.55, c=0.15, r=1/200,
 Hpop = c(100, 500, 250)
 ```
 
-To use `xde_setup_cohort` we must set up a function `F_eir` returns the value of EIR we computed above:  
 
 ```{r}
-fF_eir1 = function(EIR){return(function(t, pars){EIR})} 
-F_eir1 = fF_eir1(EIR)
-```
-
-```{r}
-xde_setup_cohort(F_eir1, Xname="SIP", HPop=Hpop, Xopts = Xo) -> test_SIP
+xde_setup_cohort(F_eir, Xname="SIP", HPop=Hpop, Xopts = Xo) -> test_SIP
 ```
 
 ```{r}

vignettes/human_sis.Rmd

@@ -102,10 +102,11 @@ params$nHosts = 1
 params$eir = EIR
 params$EIR = list() 
 
-params$eir = as.vector(EIR)
-F_eir = function(t, pars){
-  pars$EIR[[1]] = pars$eir
-}
+fF_eir = function(EIR){
+  EIR = as.vector(EIR) 
+  return(function(t, bday=0, scale=1){EIR})
+} 
+F_eir = fF_eir(EIR)
 
 params = make_parameters_demography_null(pars = params, H=H) 
 
@@ -153,15 +154,8 @@ Hpop = c(100, 500, 250)
 ```
 
 
-To use `xde_setup_cohort` we must set up a function `F_eir` returns the value of EIR we computed above:  
-
-```{r}
-fF_eir1 = function(EIR){return(function(t, pars){EIR})} 
-F_eir1 = fF_eir1(EIR)
-```
-
 ```{r}
-xde_setup_cohort(F_eir1, Xname="SIS", HPop=Hpop, Xopts = Xo) -> test_SIS
+xde_setup_cohort(F_eir, Xname="SIS", HPop=Hpop, Xopts = Xo) -> test_SIS
 ```
 
 ```{r}