Merge pull request #157 from dd-harp/dev

Multi-species upgrade, compute_var_, documentation

DESCRIPTION

@@ -23,6 +23,7 @@ Imports:
   deSolve,
   rootSolve,
   expm,
+  graphics,
   MASS
 Suggests:
     ggplot2,

NAMESPACE

@@ -344,6 +344,11 @@ export(compute_kappa)
 export(compute_local_frac)
 export(compute_terms)
 export(compute_terms_steady)
+export(compute_vars_aqua)
+export(compute_vars_cohort)
+export(compute_vars_full)
+export(compute_vars_human)
+export(compute_vars_mosy)
 export(dEIPdt)
 export(dHdt)
 export(dLdt)

R/compute.R → R/compute_terms.R

@@ -19,17 +19,41 @@ compute_terms <- function(varslist, deout, pars, s, i) {
 #' @export
 compute_terms.xde <- function(varslist, deout, pars, s, i) {
   time = deout[,1]
-  eir = c()
-  kappa = c()
+
+  eir = list()
+  if(pars$nHosts>0) for(i in 1:pars$nHosts)
+     eir[[i]] = matrix(0, nrow = length(time), ncol = pars$Hpar[[i]]$nStrata)
+
+  kappa = list()
+  fqZ = list()
+  fqM = list()
+  for(s in 1:pars$nVectors){
+     kappa[[s]] = matrix(0, nrow = length(time), ncol = pars$nPatches)
+     fqZ[[s]] = matrix(0, nrow = length(time), ncol = pars$nPatches)
+     fqM[[s]] = matrix(0, nrow = length(time), ncol = pars$nPatches)
+  }
+
   for (ix in 1:length(time)){
     yt = deout[ix,-1]
-    pars = Transmission(time[ix], yt, pars)
-    eir = c(eir, pars$EIR[[1]])
-    kappa = c(kappa, pars$kappa[[1]])
+    pars = compute_vars_full(time[ix], yt, pars)
+
+    for(i in 1:pars$nHosts)
+      eir[[i]][ix,] = pars$EIR[[i]]
+
+
+    for(s in 1:pars$nVectors){
+      kappa[[s]][ix,] = pars$kappa[[i]]
+      fqZ[[s]][ix,] = F_fqZ(time[ix], yt, pars, s)
+      fqM[[s]][ix,] = F_fqM(time[ix], yt, pars, s)
+    }
   }
-  ni = compute_NI(deout, pars, i)
-  fqZ = compute_fqZ(deout, pars, s)
-  pr = F_pr(varslist, pars, i)
+  pr = list()
+  ni = list()
+  for(i in 1:pars$nHosts){
+    ni[[i]] = compute_NI(deout, pars, i)
+    pr[[i]] = F_pr(varslist, pars, i)
+  }
+
   return(list(time=time,eir=eir,pr=pr,ni=ni,kappa=kappa,fqZ=fqZ))
 }
 

R/compute_variables.R

@@ -0,0 +1,146 @@
+
+#' @title Compute other variables at time t
+#' @description Compute everything but the derivatives for the generalized
+#' spatial differential equation model
+#' @param t current simulation time
+#' @param y state vector
+#' @param pars a [list]
+#' @return **pars** a [list]
+#' @export
+compute_vars_full <- function(t, y, pars) {
+
+  # set the values of exogenous forcing variables
+  pars <- Abiotic(t, pars)
+  pars <- Shock(t, pars)
+  pars <- Control(t, y, pars)
+  pars <- Behavior(t, y, pars)
+  pars <- Visitors(t, pars)
+  pars <- VectorControlEffects(t, y, pars)
+  pars <- Resources(t, y, pars)
+
+  # set and modify the baseline bionomic parameters
+  pars <- Bionomics(t, y, pars)
+  pars <- VectorControlEffectSizes(t, y, pars)
+
+  # egg laying: compute eta
+  pars <- EggLaying(t, y, pars)
+
+  # emergence: compute Lambda
+  pars <- Emergence(t, y, pars)
+
+  # compute beta, EIR, and kappa
+  pars <- Transmission(t, y, pars)
+
+  # compute the FoI
+  pars <- Exposure(t, y, pars)
+
+  return(pars)
+}
+
+#' @title Compute other variables at time t
+#' @description Compute everything but the derivatives for a human-only
+#' differential equation model
+#' @param t current simulation time
+#' @param y state vector
+#' @param pars a [list]
+#' @return **pars** a [list]
+#' @export
+compute_vars_human <- function(t, y, pars) {
+
+  # set the values of exogenous forcing variables
+  pars <- Abiotic(t, pars)
+  pars <- Shock(t,  pars)
+  pars <- Control(t, y, pars)
+  pars <- Behavior(t, y, pars)
+  pars <- Resources(t, y, pars)
+
+  # set and modify the baseline mosquito bionomic parameters
+  pars <- MBionomics(t, y, pars, 1)
+  pars <- VectorControlEffectSizes(t, y, pars)
+
+  # compute beta, EIR, and kappa
+  pars <- Transmission(t, y, pars)
+
+  # compute the FoI
+  pars <- Exposure(t, y, pars)
+
+  return(pars)
+}
+
+#' @title Compute other variables at time t
+#' @description Compute everything but the derivatives for a mosquito-only
+#' differential equation model
+#' @param t current simulation time
+#' @param y state vector
+#' @param pars a [list]
+#' the appropriate adult mosquito model
+#' @return **pars** a [list]
+#' @export
+compute_vars_mosy <- function(t, y, pars) {
+
+  # set the values of exogenous forcing variables
+  pars <- Abiotic(t, pars)
+  pars <- Shock(t, pars)
+  pars <- Control(t, y, pars)
+  pars <- Behavior(t, y, pars)
+  #pars <- Resources(t, y, pars)
+
+  # set baseline mosquito bionomic parameters
+  pars <- Bionomics(t, y, pars)
+  pars <- VectorControlEffectSizes(t, y, pars)
+
+  # egg laying: compute eta
+  pars <- EggLaying(t, y, pars)
+
+  # emergence: compute Lambda
+  pars <- Emergence(t, y, pars)
+
+  return(pars)
+}
+
+#' @title Differential equation models for human cohorts
+#' @description Compute everything but the derivatives for a human cohort
+#' differential equation model
+#' @param a age of a cohort
+#' @param y state vector
+#' @param pars a [list]
+#' @param F_eir a trace function that returns the eir
+#' @return **pars** a [list]
+#' @export
+compute_vars_cohort <- function(a, y, pars, F_eir) {
+
+  # EIR: entomological inoculation rate trace
+  pars$EIR[[1]] <- F_eir(a, pars)*pars$BFpar$relativeBitingRate[[1]][[1]]
+
+  # FoI: force of infection
+  pars <- Exposure(a, y, pars)
+
+  return(pars)
+}
+
+#' @title Differential equation models for aquatic mosquito populations
+#' @description Compute everything but the derivatives for an aquatic mosquito-only
+#' differential equation model
+#' @param t current simulation time
+#' @param y state vector
+#' @param pars a [list]
+#' @return **pars** a [list]
+#' @export
+compute_vars_aqua <- function(t, y, pars) {
+
+  # set the values of exogenous forcing variables
+  pars <- Abiotic(t, pars)
+  pars <- Shock(t, pars)
+  pars <- Control(t, y, pars)
+  pars <- HabitatDynamics(t, pars)
+
+  # modify baseline mosquito bionomic parameters
+  pars <- LBionomics(t, y, pars, 1)
+  pars <- VectorControlEffectSizes(t, y, pars)
+
+  # egg laying: compute eta
+
+  pars$eggs_laid[[1]] = F_eggs(t, y, pars, 1)
+
+  return(pars)
+}

R/exposure.R

@@ -1,8 +1,14 @@
-#' @title A model for exposure. The function `F_b` must be define
+#' @title Exposure and Infection
+#' @description This function translates the local daily entomological
+#' inoculation rate (dEIR), computed as one of the transmission terms
+#' into a measure of the daily force of infection (dFoI). The daily FoI
+#' is the sum of two terms: 1) a function [F_foi] computes the local dFoI;
+#' 2) a function [travel_malaria] that computes the FoI resulting from
+#' exposure while traveling.
 #' @param t the time
-#' @param y the variables
-#' @param pars a [list]
-#' @return [list]
+#' @param y the state variables
+#' @param pars the model object as a [list]
+#' @return the function modifies **pars** and returns it: the computed FoI are stored as `pars$FoI`
 #' @export
 Exposure <- function(t, y, pars){
 
@@ -13,18 +19,24 @@ Exposure <- function(t, y, pars){
   return(pars)
 }
 
-#' @title A model for exposure
-#' @description This method dispatches on the type of `pars$FOIpar`.
-#' @param eir the daily eir
+#' @title A model for daily FoI as a function of the daily EIR.
+#' @description This function compute the daily local FoI as a function
+#' of the daily EIR and effects of partial immunity.
+#' It is computed based on a probabilistic model of exposure and
+#' possibly including environmental heterogeneity. If a model of human / host
+#' infection has terms that describe partial immunity, *e.g.* affecting
+#' pre-erythrocytic immunity to malaria called by [F_b], those effects are implemented here.
+#' The method dispatches on the type of `pars$FOIpar`.
+#' @param eir the daily eir for each stratum
 #' @param b the probability of infection, per bite
-#' @param pars a [list]
-#' @return see help pages for specific methods
+#' @param pars the model object as a [list]
+#' @return the daily, local force of infection as a [numeric] vector
 #' @export
 F_foi <- function(eir, b, pars){
   UseMethod("F_foi", pars$FOIpar)
 }
 
-#' @title Null human population births
+#' @title A Poisson model for the daily local FoI as a function of the daily EIR.
 #' @description Implements [F_foi] for a Poisson model
 #' @inheritParams F_foi
 #' @return a [numeric] vector of length `nStrata`
@@ -33,7 +45,7 @@ F_foi.pois <- function(eir, b, pars){
   b*eir
 }
 
-#' @title Null human population births
+#' @title The daily FoI as a function of the daily EIR under a negative binomial model of exposure.
 #' @description Implements [F_foi] for a negative binomial model
 #' @inheritParams F_foi
 #' @return a [numeric] vector of length `nStrata`

R/human-SIP.R

@@ -272,15 +272,15 @@ xde_plot_X.SIP = function(pars, i=1, clrs=c("darkblue", "darkred", "darkgreen"),
 #' @export
 xde_lines_X_SIP = function(XH, pars, clrs=c("darkblue", "darkred", "darkgreen"), llty=1){
   with(XH,{
-    if(pars$Hpar[[i]]$nStrata==1) {
+    if(pars$Hpar[[1]]$nStrata==1) {
       lines(time, S, col=clrs[1], lty = llty[1])
       lines(time, I, col=clrs[2], lty = llty[1])
       lines(time, P, col=clrs[3], lty = llty[1])
     }
-    if(pars$Hpar[[i]]$nStrata>1){
+    if(pars$Hpar[[1]]$nStrata>1){
       if (length(clrs)==1) clrs=matrix(clrs, 3, pars$Hpar[[i]]$nStrata)
       if (length(llty)==1) llty=rep(llty, pars$Hpar[[i]]$nStrata)
-      for(i in 1:pars$Hpar[[i]]$nStrata){
+      for(i in 1:pars$Hpar[[1]]$nStrata){
         lines(time, S[,i], col=clrs[1,i], lty = llty[i])
         lines(time, I[,i], col=clrs[2,i], lty = llty[i])
         lines(time, P[,i], col=clrs[3,i], lty = llty[i])

R/human-SIS.R

@@ -230,29 +230,30 @@ xde_plot_X.SIS = function(pars, i=1, clrs=c("darkblue","darkred"), llty=1, stabl
          plot(time, 0*time, type = "n", ylim = c(0, max(H)),
               ylab = "# Infected", xlab = "Time"))
 
-  xde_lines_X_SIS(vars$XH[[i]], pars, clrs, llty)
+
+  xde_lines_X_SIS(vars$XH[[i]], pars$Hpar[[i]]$nStrata, clrs, llty)
 }
 
 
 #' Add lines for the density of infected individuals for the SIS model
 #'
 #' @param XH a list with the outputs of parse_deout_X_SIS
-#' @param pars a list that defines an `exDE` model (*e.g.*,  generated by `xde_setup()`)
+#' @param nStrata the number of population strata
 #' @param clrs a vector of colors
 #' @param llty an integer (or integers) to set the `lty` for plotting
 #'
 #' @export
-xde_lines_X_SIS = function(XH, pars, clrs=c("darkblue","darkred"), llty=1){
+xde_lines_X_SIS = function(XH, nStrata, clrs=c("darkblue","darkred"), llty=1){
   with(XH,{
-    if(pars$Hpar[[i]]$nStrata==1) {
+    if(nStrata==1) {
       lines(time, S, col=clrs[1], lty = llty[1])
       lines(time, I, col=clrs[2], lty = llty[1])
     }
-    if(pars$Hpar[[i]]$nStrata>1){
-      if (length(clrs)==2) clrs=matrix(clrs, 2, pars$Hpar[[i]]$nStrata)
-      if (length(llty)==1) llty=rep(llty, pars$Hpar[[i]]$nStrata)
+    if(nStrata>1){
+      if (length(clrs)==2) clrs=matrix(clrs, 2, nStrata)
+      if (length(llty)==1) llty=rep(llty, nStrata)
 
-      for(i in 1:pars$Hpar[[i]]$nStrata){
+      for(i in 1:nStrata){
         lines(time, S[,i], col=clrs[1,i], lty = llty[i])
         lines(time, I[,i], col=clrs[2,i], lty = llty[i])
       }

R/human-trace.R

@@ -96,8 +96,7 @@ make_Hpar_trace = function(nPatches, Hopts=list()){
     Hpar = list()
     Hpar$wts_f = rep(1, nPatches)
     Hpar$H = rep(1, nPatches)
-    Hpar$residence = c(1:nPatches)
-    Hpar$TaR = diag(1, nPatches)
+    Hpar$nStrata = nPatches
     return(Hpar)
   })}