Merge pull request #154 from dd-harp/dev

Multi-Species mop up

R/adult-Gtrace.R

@@ -34,7 +34,7 @@ F_fqM.Gtrace <- function(t, y, pars, s) {
 #' @return a [numeric] vector of length `nPatches`
 #' @export
 F_eggs.Gtrace <- function(t, y, pars, s) {
-  with(pars$MYZpar[[s]], return(Gm*Gf(t, pars$MYZpar[[s]])))
+  with(pars$MYZpar[[s]], return(scale*Gf(t)))
 }
 
 #' @title Derivatives for aquatic stage mosquitoes
@@ -70,8 +70,8 @@ make_MYZpar_Gtrace = function(nPatches, MYZopts, Gm = 1, Gf=NULL){
     MYZpar <- list()
     class(MYZpar) <- "Gtrace"
 
-    MYZpar$Gm <- checkIt(Gm, nPatches)
-    if(is.null(Gf)) Gf = function(t, MYZpar){return(1)}
+    MYZpar$scale <- checkIt(Gm, nPatches)
+    if(is.null(Gf)) Gf = function(t){return(1)}
     MYZpar$Gf = Gf
 
     return(MYZpar)
@@ -110,14 +110,15 @@ parse_deout_MYZ.Gtrace <- function(deout, pars, s) {
 #' @param Gf a [function] of the form Gf(t, pars) that computes temporal fluctuations
 #' @return none
 #' @export
-make_parameters_MYZ_Gtrace <- function(pars, Gm, Gf) {
+make_parameters_MYZ_Gtrace <- function(pars, Gm, Gf=NULL) {
   stopifnot(is.numeric(Gm))
   MYZpar <- list()
   class(MYZpar) <- 'Gtrace'
   xde <- "trace"
   class(xde) <- "trace"
   MYZpar$xde <- xde
-  MYZpar$Gm <- Gm
+  MYZpar$scale <- checkIt(Gm, pars$nPatches)
+  if(is.null(Gf)) Gf = function(t){return(1)}
   MYZpar$Gf = Gf
   pars$MYZpar[[1]] <- MYZpar
   return(pars)

R/adult-Ztrace.R

@@ -19,7 +19,7 @@ MBionomics.Ztrace <- function(t, y, pars, s) {
 #' @return a [numeric] vector of length `nHabitats`
 #' @export
 F_fqZ.Ztrace <- function(t, y, pars, s) {
-  with(pars$MYZpar[[s]], return(f*q*Zm*Zf(t, pars)))
+  with(pars$MYZpar[[s]], return(f*q*scale*Zf(t)))
 }
 
 #' @title Blood feeding rate of the infective mosquito population
@@ -82,8 +82,8 @@ make_MYZpar_Ztrace = function(nPatches, MYZopts,
     MYZpar$f0 <- MYZpar$f
     MYZpar$q0 <- MYZpar$q
 
-    MYZpar$Zm <- checkIt(Zm, pars$nPatches)
-    if(is.null(Zf)) Zf = function(t, pars){return(0*t + 1)}
+    MYZpar$scale <- checkIt(Zm, pars$nPatches)
+    if(is.null(Zf)) Zf = function(t){return(1)}
     MYZpar$Zf <- Zf
 
     return(MYZpar)
@@ -125,18 +125,19 @@ make_indices_MYZ.Ztrace <- function(pars, s) {
 #' @param Zf a [function] of the form Zf(t, pars) that computes temporal fluctuations
 #' @return none
 #' @export
-make_parameters_MYZ_Ztrace <- function(pars, Zm, f, q, Zf) {
+make_parameters_MYZ_Ztrace <- function(pars, Zm, f, q, Zf=NULL) {
   stopifnot(is.numeric(Zm))
   MYZpar <- list()
   class(MYZpar) <- 'Ztrace'
   xde <- "trace"
   class(xde) <- "trace"
   MYZpar$xde <- xde
-  MYZpar$Zm <- Zm
   MYZpar$f0 <- f
   MYZpar$f <- f
   MYZpar$q0 <- q
   MYZpar$q <- q
+  MYZpar$scale <- checkIt(Zm, pars$nPatches)
+  if(is.null(Zf)) Zf = function(t){return(1)}
   MYZpar$Zf = Zf
   pars$MYZpar <- MYZpar
   return(pars)

R/aquatic-trace.R

@@ -15,7 +15,7 @@ LBionomics.trace <- function(t, y, pars, s) {
 #' @return a [numeric] vector of length `nHabitats`
 #' @export
 F_alpha.trace <- function(t, y, pars, s) {
-  pars$Lpar[[s]]$Lt(t)
+  with(pars$Lpar[[s]], scale*Lt(t))
 }
 
 #' @title Derivatives for aquatic stage mosquitoes
@@ -58,9 +58,8 @@ make_Lpar_trace = function(nHabitats, Lopts=list(), Lambda=1000, Lt = NULL){
   with(Lopts,{
     Lpar = list()
     class(Lpar) <- "trace"
-    Lpar$Lambda = checkIt(Lambda, nHabitats)
-    #if(is.null(Lt)) Lt = function(t, Lpar){Lpar$Lambda}
-    if(is.null(Lt)) Lt = function(t){Lambda}
+    Lpar$scale = checkIt(Lambda, nHabitats)
+    if(is.null(Lt)) Lt = function(t){1}
     Lpar$Lt = Lt
     return(Lpar)
 })}
@@ -93,8 +92,8 @@ make_parameters_L_trace <- function(pars, Lambda, Lt=NULL) {
   stopifnot(is.numeric(Lambda))
   Lpar <- list()
   class(Lpar) <- 'trace'
-  Lpar$Lambda = checkIt(Lambda, pars$nHabitats)
-  if(is.null(Lt)) Lt = function(t){Lambda}
+  Lpar$scale = checkIt(Lambda, pars$nHabitats)
+  if(is.null(Lt)) Lt = function(t){return(1)}
   Lpar$Lt = Lt
   pars$Lpar[[1]] <- Lpar
   return(pars)

R/diffeqn.R

@@ -37,12 +37,18 @@ xDE_diffeqn <- function(t, y, pars) {
   pars <- Exposure(t, y, pars)
 
   # state derivatives
-  dL <- c()
-  for(s in 1:pars$nVectors) dL <- c(dL, dLdt(t, y, pars, s))
-  dMYZ <- c()
-  for(s in 1:pars$nVectors) dMYZ <- c(dMYZ, dMYZdt(t, y, pars, s))
-  dX <- c()
-  for(i in 1:pars$nHosts) dX <- c(dX, dXdt(t, y, pars, i))
+  dL <- dLdt(t, y, pars, 1)
+  dMYZ <- dMYZdt(t, y, pars, 1)
+  if(pars$nVectors > 1)
+    for(s in 1:pars$nVectors){
+      dL <- c(dL, dLdt(t, y, pars, s))
+      dMYZ <- c(dMYZ, dMYZdt(t, y, pars, s))
+    }
+
+  dX <- dXdt(t, y, pars, 1)
+  if(pars$nHosts > 1)
+    for(i in 1:pars$nHosts)
+      dX <- c(dX, dXdt(t, y, pars, i))
 
   return(list(c(dL, dMYZ, dX)))
 }
@@ -75,8 +81,10 @@ xDE_diffeqn_human <- function(t, y, pars) {
   pars <- Exposure(t, y, pars)
 
   # state derivatives
-  dX <- c()
-  for(i in 1:pars$nHosts) dX <- c(dX, dXdt(t, y, pars, i))
+  dX <- dXdt(t, y, pars, 1)
+  if(pars$nHosts > 1)
+    for(i in 1:pars$nHosts)
+      dX <- c(dX, dXdt(t, y, pars, i))
 
   return(list(c(dX)))
 }
@@ -112,10 +120,13 @@ xDE_diffeqn_mosy <- function(t, y, pars) {
   pars <- Emergence(t, y, pars)
 
   # state derivatives
-  dL <- c()
-  for(s in 1:pars$nVectors) dL <- c(dL, dLdt(t, y, pars, s))
-  dM <- c()
-  for(s in 1:pars$nVectors) dM <- c(dM, dMYZdt(t, y, pars, s))
+  dL <- dLdt(t, y, pars, 1)
+  dM <- dMYZdt(t, y, pars, 1)
+  if (pars$nVectors > 1)
+    for(s in 2:pars$nVectors){
+      dL <- c(dL, dLdt(t, y, pars, s))
+      dM <- c(dM, dMYZdt(t, y, pars, s))
+    }
 
   return(list(c(dL, dM)))
 }
@@ -138,8 +149,10 @@ xDE_diffeqn_cohort <- function(a, y, pars, F_eir) {
   pars <- Exposure(a, y, pars)
 
   # state derivatives
-  dX <- c()
-  for(i in 1:pars$nHosts) dX <- c(dX, dXdt(t, y, pars, i))
+  dX <- dXdt(t, y, pars, 1)
+  if(pars$nHosts > 1)
+    for(i in 1:pars$nHosts)
+      dX <- c(dX, dXdt(t, y, pars, i))
 
   return(list(c(dX)))
 }
@@ -168,8 +181,10 @@ xDE_diffeqn_aquatic <- function(t, y, pars) {
   pars <- EggLaying(t, y, pars)
 
   # state derivatives
-  dL <- c()
-  for(s in 1:pars$nVectors) dL <- c(dL, dLdt(t, y, pars, s))
+  dL <- dLdt(t, y, pars, 1)
+  if(pars$nVectors > 1)
+    for(s in 1:pars$nVectors)
+      dL <- c(dL, dLdt(t, y, pars, s))
 
   return(list(c(dL)))
 }

R/human-trace.R

@@ -6,7 +6,7 @@
 #' @return a [numeric] vector of length `nStrata`
 #' @export
 F_X.trace <- function(t, y, pars, i) {
-  with(pars$Xpar[[i]], kappa)
+  with(pars$Xpar[[i]], scale*Kf(t))
 }
 
 #' @title Size of the human population
@@ -80,8 +80,8 @@ make_Xpar_trace = function(nPatches, Xopts=list(),
     Xpar = list()
     class(Xpar) <- "trace"
 
-    Xpar$kappa = checkIt(kappa, nPatches)
-    if(is.null(Kf)) Kf = function(t, y, pars){return(1)}
+    Xpar$scale = checkIt(kappa, nPatches)
+    if(is.null(Kf)) Kf = function(t){return(1)}
     Xpar$Kf = Kf
     return(Xpar)
 })}
@@ -128,7 +128,8 @@ parse_deout_X.trace <- function(deout, pars,i) {
 make_parameters_X_trace <- function(pars, kappa) {
   Xpar <- list()
   class(Xpar) <- c('trace')
-  Xpar$kappa <- kappa
+  Xpar$scale = checkIt(kappa, pars$nPatches)
+  Xpar$Kf = function(t){return(1)}
   pars$Xpar[[1]] <- Xpar
   return(pars)
 }

_pkgdown.yml

@@ -7,6 +7,10 @@ navbar:
     href: reference/index.html
   - text: Articles
     menu:
+    - text: Understanding exDE
+      href: articles/Understanding_exDE.html
+    - text: Modular Forms for Disease Dynamics
+      href: articles/modularity.html
     - text: 5-3-4 Spatial Example
       href: articles/ex_534.html
     - text: Spatial Metrics

vignettes/Understanding_exDE.Rmd

@@ -0,0 +1,77 @@
+---
+title: "Understanding exDE"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Understanding exDE}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+
+## Introduction 
+
+The goal of developing `exDE` was to lower the costs of building models  (*e.g.* time spent developing, coding, verifying, debugging, & *etc.*) that are realistic enough to support decisions affecting malaria policies. The operating assumption of the development process has been that robust decision support for malaria policies would put demands on model-building that differ from what is expected from a scientific publication. In some ways, it was the need for *realism* and the associated *computational complexity* that drove the development of several *individual-based models* (IBMs). While models developed as systems of differential equations can not replicate some advantages that come from being *individual-based,* we wanted software that could handle computational complexity just as well as IBMs. The advantage is that the resulting models would be much easier to analyze. 
+
+We wanted a framework that could *scale complexity,* to start with a simple model, and then progressively modifying the models to add realism including:
+
++ realistic human demography, including age structure, births & deaths, and migration. 
+
++ spatial heterogeneity and spatial dynamics, including human mobility and mosquito dispersal 
+
++ exogenous forcing by weather and other factors
+
++ multiple host and vector species or types 
+
+ If we think of a model as defining a *skill set,* or output variables that naturally represent a subset of all possible features, then the framework should make it easy to build models that have the right subset of features.  The models should be -- as Einstein suggested -- as simple as possible, but no simpler. While most people would agree with Einstein in principle, he provided no usable advice about how to do this. How would you know a model had the right level of complexity unless, at some point, a model had been developed that was clearly too complex? 
+
+To support decisions, model building must be *nimble.*  As malaria programs integrate sophisticated analytics into their decision-making, the conversation in a room can shift. Over time, programmatic priorities and needs can change. To keep up, the model builders *ought* to have the capability of building new models that could address the new concerns. While this might not be possible to do in real time, it should be possible to have a new model developed with preliminary results within a week or so.  
+
+With these goals in mind, we developed a mathematical framework that could make computation for dynamical systems **modular.** The mathematical basis for modularity was described in [Spatial Dynamics of Malaria Transmission](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010684){target="_blank"}. We have also developed a *vignette* that describes the [modular forms](modularity.html){target="_blank"} that we use to describe models implemented in `exDE`.
+
+## Structural Elements
+
+Each model has several structural parameters that describe the basic features of a model: 
+
++ `nPatches` -- the number of distinct *patches* in a model. 
+
+    - A model has *spatial dynamics* if `nPatches`$>1,$ making it necessary to 
+
++ Aqu
++ `nHabitats` -- the number of distinct *patches* in a model. 
+
+
+## Dynamical Components
+
+We identified three inseparable chunks of any model that would need to be internally coherent but that could be represented in several different ways. The chunks represented five distinct processes: 
+
++ Models for human ecology and parasite / pathogen infection dynamics would appear in one dynamical component. The part that computes the dynamics of infection and immunity was called $\cal X$, and the part that describes human demography was called $\cal H.$  These two components can't be separated in any easy way, so we call this chunk $\cal XH$.
+
+    - The derivatives for a model of this type are computed by a `S3` function `dXdt` 
+
+    - The parameters for the  model are in an object called `Xpar` and the model dispatches on `class(Xpar)`
+
+    - Since there could be multiple host species, the $i^{th}$ species is in  `pars$Xpar[[i]].`
+
+    - The function `dXdt` dispatches on `class(pars$Xpar[[i]])` 
+
++ Models for adult mosquito ecology and parasite / pathogen infection dynamics would appear in a second dynamical component. The part that computes the dynamics of infection was called $\cal YZ$, and the part that computes mosquito population dynamics was called $\cal M$. These two components can't be separated in any easy way, so we call this chunk $\cal MYZ$.   
+
++ Models for aquatic mosquito ecology were called $\cal L$. 
+
+## Interfaces 
+
+In developing a modular framework, we recognized the need to develop a rigorous yet flexible *interface* that would allow different dynamical components to interact. 
+To connect these two dynamical components, we developed two well-defined interfaces: 
+
++ A *model* for **blood feeding** and parasite / pathogen **transmission** by mosquitoes. This model *should* constrain mosquito blood feeding rates and the human fraction in sensible ways. 
+
++ A description of the locations of aquatic habitats and a *model* for **egg laying** and **emergence.**
+
+In a related vignette, a Ross-Maconald model is presented in its [modular form](modularity.html)

vignettes/adult_RM.Rmd

@@ -16,7 +16,7 @@ knitr::opts_chunk$set(
 
 The Ross-Macdonald adult mosquito model fulfills the generic interface of the adult mosquito component. 
 
-Here, we use a Ross-Macdonald model based on a model first published in 1982 by Joan Aron and Robert May (Aron JL, and May RM (1982). The population dynamics of malaria. In *The Population Dynamics of Infectious Diseases: Theory and Applications,* R. M. Anderson, ed. (Springer US), pp. 139–179). It includes state variables for total mosquito density $M$, infected mosquito density $Y$, and infectious mosquito density $Z$. Because of the interest in parity, the model has been extended to include a variable that tracks parous mosquitoes, $P$.
+Here, we use a Ross-Macdonald model based on a model first published in 1982 by Joan Aron and Robert May^[The population dynamics of malaria. In *The Population Dynamics of Infectious Diseases: Theory and Applications,* R. M. Anderson, ed. (Springer US), pp. 139–179. [online](https://link.springer.com/chapter/10.1007/978-1-4899-2901-3_5){target="_blank"}]. It includes state variables for total mosquito density $M$, infected mosquito density $Y$, and infectious mosquito density $Z$. Because of the interest in parity, the model has been extended to include a variable that tracks parous mosquitoes, $P$.
 
 # Differential Equations 
 
@@ -54,9 +54,7 @@ The system of ODEs is the same as above except for the equation giving the rate
 $$
 \dot{Z} = e^{-\Omega \tau} \cdot \mbox{diag}(fq\kappa) \cdot (M-Y) - \Omega \cdot Z
 $$
-The resulting set of equations is similar in spirit to the simple model presented in
-[Smith & McKenzie (2004)](https://malariajournal.biomedcentral.com/articles/10.1186/1475-2875-3-13) in that
-mortality and dispersal over the EIP is accounted for, but the time lag is not. While
+The resulting set of equations is similar in spirit to the simple model presented in  Smith & McKenzie (2004)^[Smith, D.L., Ellis McKenzie, F. Statics and dynamics of malaria infection in Anopheles mosquitoes. Malar J 3, 13 (2004). [online](https://doi.org/10.1186/1475-2875-3-13){target="_blank"}]. in that mortality and dispersal over the EIP is accounted for, but the time lag is not. While
 transient dynamics of the ODE model will not equal the DDE model, they have the same
 equilibrium values, and so for numerical work requiring finding equilibrium points, the faster
 ODE model can be safely substituted.
@@ -115,15 +113,19 @@ We can use the same equation for $P$ as above.
 
 # Example
 
+
 ```{r, message=FALSE, warning=FALSE}
 library(exDE)
-#devtools::load_all()
 library(expm)
 library(deSolve)
 library(data.table)
 library(ggplot2)
 ```
 
+```{r, echo=F, eval=F}
+#devtools::load_all()
+``` 
+
 Here we show an example of starting and solving a model at equilibrium. Please note that this only
 runs this adult mosquito model and that most users should read [our fully worked example](ex_534.html)
 to run a full simulation.

vignettes/aqua_basic.Rmd

@@ -55,6 +55,10 @@ library(data.table)
 library(ggplot2)
 ```
 
+
+```{r, echo=F}
+#devtools::load_all()
+``` 
 ## The long way 
 
 Here we run a simple example with 3 aquatic habitats at equilibrium. We use `exDE::make_parameters_L_basic` to

vignettes/aqua_trace.Rmd

@@ -24,6 +24,11 @@ library(deSolve)
 library(data.table)
 library(ggplot2)
 ```
+
+
+```{r, echo=F}
+#devtools::load_all()
+``` 
 
 This is the null model of aquatic mosquito dynamics; there are no endogenous dynamics
 and the model is simply specified by `Lambda`, the rate at which adult mosquitoes