multi-species transition bug fixes

R/Xmod-SIS.R

@@ -3,34 +3,46 @@
 #' @inheritParams split_stratum_by_biting
 #'
 #' @return pars a list
-#' @param i the stratum to split
+#' @param i the host species index
+#' @param j the stratum to split
 #' @param p the fraction that gets multiplied by `fac`
 #' @param fac a factor
 #' @return a [list]
 #' @export
-split_stratum_by_biting.SIS = function(pars, i, p, fac){
-  stopifnot(i <= pars$nStrata)
-  pars$nStrata = pars$nStrata + 1
+split_stratum_by_biting.SIS = function(pars, i, j, p, fac){
+  stopifnot(i <= pars$Hpar[[i]]$nStrata)
+  nj_ix = pars$Hpar[[i]]$nStrata + 1
+  pars$Hpar[[i]]$nStrata = nj_ix
+
 
   Xpar = pars$Xpar[[1]]
-  Xpar$b <- c(Xpar$b, Xpar$b[i])
-  Xpar$r <- c(Xpar$r, Xpar$r[i])
-  Xpar$c <- c(Xpar$c, Xpar$c[i])
-
-  Xinits = pars$inits[[1]]
-  Xinits$X0 <- c(Xinits$X0, Xinits$X0[i])
-  Xinits$X0[i] <- Xinits$X0[i]*p
-  Xinits$X0[pars$nStrata] <- Xinits$X0[i]*(1-p)
-
-  pars$Hpar$residence = c(pars$Hpar$residence, pars$Hpar$residence[i])
-  pars$Hpar$wts_f = c(pars$Hpar$wts_f, pars$Hpar$wts_f[i])
-  pars$Hpar$wts_f[pars$nStrata] = pars$Hpar$wts_f[i]*fac
-  pars$Hpar$wts_f[i] = pars$Hpar$wts_f[i]/fac
-  pars$Hpar$rbr = with(pars$Hpar, wts_f*sum(H)/sum(wts_f*H))
-  pars$Hpar$H = c(pars$Hpar$H, pars$Hpar$H[i])
-  pars$Hpar$H[pars$nStrata] = pars$Hpar$H[i]*p
-  pars$Hpar$H[i] = pars$Hpar$H[i]*(1-p)
-  pars$Hpar$TaR = cbind(pars$Hpar$TaR, pars$Hpar$TaR[,i])
+  Xpar$b <- c(Xpar$b, Xpar$b[j])
+  Xpar$r <- c(Xpar$r, Xpar$r[j])
+  Xpar$c <- c(Xpar$c, Xpar$c[j])
+  pars$Xpar[[1]] = Xpar
+
+  Xinits = pars$Xinits[[1]]
+  Xinits$S <- c(Xinits$S, p*Xinits$S[j])
+  Xinits$S[j] <- Xinits$S[j]*(1-p)
+  Xinits$I <- c(Xinits$I, p*Xinits$I[j])
+  Xinits$I[j] <- Xinits$I[j]*(1-p)
+  pars$Xinits[[1]] = Xinits
+
+  residence = pars$BFpar$residence[[i]]
+  pars$BFpar$residence[[i]] = c(residence, residence[j])
+  for(s in 1:pars$nVectors){
+    wts_f = pars$BFpar$searchWts[[i]][[s]]
+    pars$BFpar$searchWts[[i]][[s]] = c(wts_f, fac*wts_f[j])
+  }
+
+  H =  pars$Hpar[[i]]$H
+  pars$Hpar[[i]]$H = c(H, p*H[j])
+  pars$Hpar[[i]]$H[j] = H[j]*(1-p)
+
+  TimeSpent = pars$BFpar$TimeSpent[[i]]
+  TimeSpent = cbind(TimeSpent, TimeSpent[,j])
+
+  for(s in 1:pars$nVectors) pars = make_TaR(0, pars, i, 1)
 
   pars <- exDE::make_indices(pars)
 

R/plot-terms.R

@@ -2,44 +2,47 @@
 #' Plot EIR(t) *vs.* the PR(t)
 #'
 #' @param pars a [list] that defines an `exDE` model (*e.g.*,  generated by `exDE::xde_setup()`)
+#' @param i the host species index
 #' @param clrs a [character] vector of colors
 #' @param llty an [integer] that specifies `lty` for plotting
 #' @param stable a [logical] set to FALSE for `orbits` and TRUE for `stable_orbits`
 #' @param add_axes a [logical] to plot add_axes only if FALSE
 #'
 #' @export
-plot_eirVpr <- function(pars, clrs="black", llty=1, stable=FALSE, add_axes=TRUE){
+plot_eirVpr <- function(pars, i=1, clrs="black", llty=1, stable=FALSE, add_axes=TRUE){
   vars=with(pars$outputs,if(stable==TRUE){stable_orbits}else{orbits})
 
-  if(add_axes==TRUE)
-  with(vars$terms,{
-    aeir = 365*eir
-      plot(aeir, pr, type = "l", xaxt="n", lty = llty,
+  pr = vars$terms$pr[[i]]
+  eir = vars$terms$eir[[i]]
+  aeir = 365*eir
+  if(add_axes==TRUE){
+    plot(aeir, pr, type = "n", xaxt="n", lty = llty,
            xlab = "aEIR", ylab = "PR",
            xlim = range(0, aeir), ylim = c(0,1), col = clrs)
       axis(1, 10^c(-1:3), c(".1", "1", "10", "100","1000"))
-  })
-  lines_eirVpr(vars$terms, pars, clrs, llty)
+  }
+  lines_eirVpr(eir, pr, pars$Hpar[[i]]$nStrata, clrs, llty)
 }
 
 #' Add lines for the EIR(t) *vs.* the PR(t)
 #'
-#' @param terms a [list] with the outputs of `exDE::parse_deout()`
-#' @param pars a [list] that defines an `exDE` model (*e.g.*,  generated by `exDE::xde_setup()`)
+#' @param eir the daily EIR
+#' @param pr the parasite rate
+#' @param nStrata the number of population strata
 #' @param clrs a [character] vector of colors
 #' @param llty an [integer] (or integers) that specifies `lty` for plotting
 #'
 #' @export
-lines_eirVpr <- function(terms, pars, clrs= "black", llty = 1){
-  with(terms,{
+lines_eirVpr <- function(eir, pr, nStrata, clrs= "black", llty = 1){
     aeir = 365*eir
-    if(pars$nStrata==1) lines(aeir, pr, col=clrs[1], lty = llty[1])
-    if(pars$nStrata>1){
-      if(length(clrs)==1) clrs=rep(clrs, pars$nStrata)
-      if(length(llty)==1) llty=rep(llty, pars$nStrata)
-      for(i in 1:pars$nStrata)
+    if(nStrata==1)
+      lines(aeir, pr, col=clrs[1], lty = llty[1])
+    if(nStrata>1){
+      if(length(clrs)==1) clrs=rep(clrs, nStrata)
+      if(length(llty)==1) llty=rep(llty, nStrata)
+      for(i in 1:nStrata)
         lines(aeir[,i], pr[,i], col=clrs[i], lty = llty[i])
-  }})
+  }
 }
 
 #' Plot the eir-pr scaling relationship

R/scaling.R

@@ -8,38 +8,32 @@ xde_scaling_eir = function(model, N=25){
   stopifnot(model$nPatches == 1)
   stopifnot(class(model$xde) == "cohort")
 
-  # make sure that eir & ni are output
-  model <- exDE::make_indices(model)
+  scale_base = model$EIRpar$scale
 
   # Get the
-  F_eir_base = model$F_eir
-  wts_f = as.vector(model$BFpar$searchWts[[1]][[1]])
-  H = F_H(0, get_inits(model), model, 1)
-  scl_a = sum(wts_f*H)/sum(H)
-  scl_b = stats::integrate(F_eir_base, 0, 365, pars=model)$value
-  model$scl = 1/scl_a/scl_b
-  model$F_eir = function(t, model){with(model,aeir*scl*F_eir_base(t, model))}
+  scl = with(model$EIRpar, stats::integrate(model$F_eir, 0, 365, scale=scale)$value)
 
   eir = 10^seq(-1,3, length.out=N)
   pr = 0*eir
   ni = 0*eir
   scaling = list()
-
   for(i in 1:N){
-    model$aeir = eir[i]
+    model$EIRpar$scale = eir[i]/scl
     xde_tmp <- exDE::xde_stable_orbit(model)
     tmp <- xde_tmp$outputs$stable_orbit
     H = tmp$XH[[1]]$H
-    ni = tmp$terms$ni
-    tot_pr <- rowSums(as.matrix(tmp$terms$pr[[1]]*H))/rowSums(as.matrix(H))
-    mean_ni <- rowSums(as.matrix(ni*wts_f*H))/rowSums(as.matrix(wts_f*H))
+    wts_f = as.vector(model$BFpar$searchWts[[1]][[1]])
+    pr_t = tmp$terms$pr
+    ni_t = tmp$terms$ni
+    tot_pr <- rowSums(as.matrix(pr_t*H))/rowSums(as.matrix(H))
+    mean_ni <- rowSums(as.matrix(ni_t*wts_f*H))/rowSums(as.matrix(wts_f*H))
     scaling[[i]] = with(tmp$terms, list(aeir=365*eir, eir=eir, pr=tot_pr, ni=mean_ni, pr_t = pr, ni_t = ni))
     pr[i] = mean(tot_pr)
     ni[i] = mean(mean_ni)
   }
 
   model$outputs$eirpr <- list(aeir=eir, eir=eir/365, pr=pr, ni=ni, scaling=scaling)
-  model$F_eir = F_eir_base
+  model$EIRPar$scale = scale_base
 
   return(model)
 }

R/stratify.R

@@ -2,14 +2,15 @@
 #' Split a stratum into two strata, assigning
 #'
 #' @param pars a [list] defining a model
-#' @param i the stratum to split
+#' @param i the host species index
+#' @param j the stratum to split
 #' @param p the fraction in the higher exposure stratum
 #' @param fac the factor increase
 #'
 #' @return pars a list
 #' @export
-split_stratum_by_biting = function(pars, i, p, fac){
-  UseMethod("split_stratum_by_biting", pars$Xpar)
+split_stratum_by_biting = function(pars, i, j, p, fac){
+  UseMethod("split_stratum_by_biting", pars$Xpar[[i]])
 }