Fixed trap algorithms

Trap algorithm vignette added

Author: smitdave

Diff for: NAMESPACE

@@ -106,6 +106,7 @@ export(make_tiles)
 export(net2patches)
 export(outline_points_b)
 export(outline_points_q)
+export(outline_points_s)
 export(plotDDcmf)
 export(plotDDpmf)
 export(plot_Kbb)

Diff for: R/adult-BQS.R

@@ -7,7 +7,7 @@
 #' @export
 adult_dynamics.BQS = function(t, model){
   with(model,{
-    with(c(Mpars, Mvars, terms), {
+    with(c(Mpar, Mvars, terms), {
       eggs_t = ova*psiQ*Q
       # compute variables
       B_t = Mbb %*% B + Mbq %*% Q + Mbs %*% S + Mbl %*% Lambda
@@ -89,10 +89,10 @@ init_adult_model.BQS = function(model, M0_opts){
 #' @return model, a compound [list]
 #' @export
 init_adult_model_BQS = function(model, opts=list(), B0=10, Q0=10, S0=10){with(opts,{
-  model$Mvars$B = c(B0, model$nb, 1)
-  model$Mvars$Q = c(Q0, model$nq, 1)
-  model$Mvars$S = c(S0, model$ns, 1)
-  model$Mvars$eggs = c(0, model$nq, 1)
+  model$Mvars$B = rep(B0, model$nb)[1:model$nb]
+  model$Mvars$Q = rep(Q0, model$nq)[1:model$nq]
+  model$Mvars$S = rep(S0, model$ns)[1:model$ns]
+  model$Mvars$eggs = rep(0, model$nq)[1:model$nq]
   return(model)
 })}
 
@@ -181,17 +181,17 @@ setup_bionomics_BQS = function(model, opts=list(),
                                sigb=0.1, sigq=0.1, sigf=0.1,
                                sigL=0.1, ova=20){
   with(opts,{
-    model$Mpar$setup$pB = pB
-    model$Mpar$setup$pQ = pQ
-    model$Mpar$setup$pS = pS
-    model$Mpar$setup$psiB = psiB
-    model$Mpar$setup$psiQ = psiQ
-    model$Mpar$setup$psiS = psiS
-    model$Mpar$setup$sigb = sigb
-    model$Mpar$setup$sigq = sigq
-    model$Mpar$setup$sigf = sigf
-    model$Mpar$setup$sigL = sigL
-    model$Mpar$setup$ova = ova
+    model$Mpar$pB = pB
+    model$Mpar$pQ = pQ
+    model$Mpar$pS = pS
+    model$Mpar$psiB = psiB
+    model$Mpar$psiQ = psiQ
+    model$Mpar$psiS = psiS
+    model$Mpar$sigb = sigb
+    model$Mpar$sigq = sigq
+    model$Mpar$sigf = sigf
+    model$Mpar$sigL = sigL
+    model$Mpar$ova = ova
   return(model)
 })}
 

Diff for: R/adult-interface.R

@@ -1,4 +1,4 @@
-#' Adult dynamics
+#' Update States for an Adult Mosquito Model
 #'
 #' @param model a model defined as a compound [list]
 #'
@@ -9,7 +9,7 @@ adult_dynamics = function(t, model){
   UseMethod("adult_dynamics", model$Mpar)
 }
 
-#' Save state variables
+#' Store State Variables for an Adult Mosquito Model
 #'
 #' @param t the current time
 #' @param model a model defined as a compound [list]
@@ -31,7 +31,7 @@ init_states_M = function(model){
 }
 
 
-#' Set initial values for
+#' Set Initial Values for **M** Module
 #'
 #' @param model a model defined as a compound [list]
 #' @param M0_opts a list of options to overwrite defaults
@@ -52,7 +52,7 @@ compute_diffs_M = function(model){
   UseMethod("compute_diffs_M", model$Mpar)
 }
 
-#' Setup an adult model
+#' Setup Adult Mosquito Module
 #'
 #' @param model a model defined as a compound [list]
 #' @param b a point set defining blood feeding sites

Diff for: R/compute_K_matrices.R

@@ -10,67 +10,106 @@ compute_Kqb = function(model, Tmax){
   UseMethod("compute_Kqb", model$Mpar)
 }
 
-#' Compute dispersal to lay eggs within one feeding cycle
+#' Compute net dispersal matrix to blood feed within one feeding cycle
 #'
 #' @param model a model defined as a compound [list]
 #' @param Tmax the number of time steps
 #'
 #' @return the model, a compound [list]
 #' @export
-compute_Kqb.BQS = function(model, Tmax=200){with(model,with(Mpar,{
-  nb = dim(Psi_bb)[2]
-  nq = dim(Psi_qq)[2]
-  ns = dim(Psi_ss)[2]
-
-  Mbb = pB*(1-sigb)*(1-psiB)*Psi_bb
-  success = pB*(1-sigb)*psiB*Psi_bb
-  Mbq = pB*psiB*Psi_bq
-  Mbs = pB*sigb*psiB*Psi_bs
-  Mqb = pQ*(1-sigf)*psiQ*Psi_qb
-  Mqq = pQ*(1-sigq)*(1-psiQ)*Psi_qq
-  Mqs = pQ*(sigf*psiQ + sigq*(1-psiQ))*Psi_qs
-  Msb = pS*psiS*Psi_sb
-  Msq = 0*t(Mqs)
-  Mss = (1-psiS)*Psi_ss
-
-  M1 = cbind(Mbb, Mqb, Msb, 0*Mbb)
-  M2 = cbind(0*Mbq, Mqq, Msq, 0*Mbq)
-  M3 = cbind(Mbs, Mqs, Mss, 0*Mbs)
-  M4 = cbind(success, 0*Mqb, 0*Msb, diag(1,nb))
-  M = rbind(M1, M2, M3, M4)
+compute_Kqb.BQ = function(model, Tmax=100){with(model, with(Mpar,{
 
-  cohort = Psi_qb%*%diag(pQ,nq)
-  Cyes = diag(psiB, nb)%*%cohort
-  Cno = diag(1-psiB, nb)%*%cohort
+  cohort = Psi_qb %*% diag(pB, nb)
 
-  Kt = rbind(Cno, 0*Mqq, 0*Mqs, Cyes)
+  notpsiQ = diag(1-psiQ, nq)
+  psiQ = diag(psiQ, nq)
+  pQ = diag(pQ, nq)
 
-  for(i in 1:Tmax) Kt = M%*%Kt
-  Kt[-c(1:(nb+nq+ns)),]
+  Kbq = 0*cohort  # Success
+  Bt = cohort # Failure
 
-  model$KGV$Kqb = Kt[nb+nq+ns+c(1:nb),]
+  Kqb = psiQ %*% cohort  # Success
+  Qt = notpsiQ %*% cohort # Failure
+  for(i in 1:Tmax){
+    Kqb = Kqb + psiQ %*% Qt
+    Qt = Psi_qq %*% pQ %*% notpsiQ %*% Qt
+  }
+  model$KGV$Kqb = Kqb
   return(model)
 }))}
 
-#' Compute net dispersal matrix to blood feed within one feeding cycle
+
+#' Compute dispersal to lay eggs within one feeding cycle
 #'
 #' @param model a model defined as a compound [list]
 #' @param Tmax the number of time steps
 #'
 #' @return the model, a compound [list]
 #' @export
-compute_Kqb.BQ = function(model, Tmax=100){with(model, with(Mpar,{
-  cohort = Psi_qb %*% diag(pB, nb)
-  Kqb = diag(psiQ, nq) %*% cohort  # Success
-  Qt = diag(1-psiQ, nq) %*% cohort # Failure
-  for(i in 1:Tmax){
-    Kqb = Kqb + diag(pQ*psiQ, nq) %*% Psi_qq %*% Qt
-    Qt =diag(pQ*(1-psiQ), nq) %*% Psi_qq %*%Qt
+compute_Kqb.BQS = function(model, Tmax=200){with(model,with(Mpar,{
+
+  nb = dim(b)[1]
+  nq = dim(q)[1]
+  ns = dim(s)[1]
+
+  sigF = diag(sigf, nq)
+  nsigF = diag(1-sigf, nq)
+  npsiB = diag(1-psiB, nb)
+  psiB = diag(psiB, nb)
+  pB = diag(pB, nb)
+  nsigB = diag(1-sigb, nb)
+  sigB = diag(sigb, nb)
+
+  npsiQ = diag(1-psiQ, nq)
+  psiQ = diag(psiQ, nq)
+  pQ = diag(pQ, nq)
+  nsigQ = diag(1-sigq, nq)
+  sigQ = diag(sigq, nq)
+
+  npsiS = diag(1-psiS, ns)
+  psiS = diag(psiS, ns)
+  pS = diag(pS, ns)
+
+  Mbb = Psi_bb %*% nsigB %*% pB %*% npsiB
+  Mbq = 0*Psi_bq
+  trap = psiQ
+  Mbs = Psi_bs %*% psiS %*% pS
+
+  Mqb = Psi_qb %*% psiB %*% pB
+  Mqq = Psi_qq %*% nsigQ %*% pQ %*% npsiQ
+  Mqs = 0*t(Psi_sq)
+
+  Msb = Psi_sb %*% sigB %*% pB %*% npsiB
+  Msq = Psi_sq %*% sigQ %*% pQ %*% npsiQ
+  Mss = Psi_ss %*% npsiS %*% pS
+
+  M1 = cbind(Mbb, Mbq, Mbs, 0*Mbq)
+  M2 = cbind(Mqb, Mqq, Mqs, 0*Mqq)
+  M3 = cbind(Msb, Msq, Mss, 0*Msq)
+  M4 = cbind(0*Mqb, trap, 0*Mqs, diag(1,nq))
+  M = rbind(M1, M2, M3, M4)
+
+  B0 = pB
+  Q0 = 0*(Psi_qb %*% pB)
+  S0 = 0*(Psi_sb %*% pB)
+  T0 = 0*Q0
+
+  nix = nb + nq + ns + c(1:nq)
+  BSQT = rbind(B0, S0, Q0, T0)
+
+  looking = 1
+  while(looking > 1e-6){
+    BSQT = M%*%BSQT
+    looking = sum(BSQT[-nix,-nix])
   }
+
+  nix = 1:(nb+nq+ns)
+  Kqb = BSQT[-nix,]
   model$KGV$Kqb = Kqb
   return(model)
 }))}
 
+
 #' Compute dispersal matrix to blood feed within one feeding cycle
 #'
 #' @param model a model defined as a compound [list]
@@ -90,13 +129,18 @@ compute_Kbq = function(model, Tmax){
 #' @return the model, a compound [list]
 #' @export
 compute_Kbq.BQ = function(model, Tmax=100){with(model, with(Mpar,{
+
+  notpsiB = diag(1-psiB, nb)
+  psiB = diag(psiB, nb)
+  pB = diag(pB, nb)
+
   cohort = Psi_bq %*% diag(pQ, nq)
-  Kbq = diag(psiB, nb) %*% cohort  # Success
-  Bt = diag(1-psiB, nb) %*% cohort # Failure
+  Kbq = 0*cohort  # Success
+  Bt = cohort # Failure
 
   for(i in 1:Tmax){
-    Kbq = Kbq + diag(pB*psiB, nb) %*% Psi_bb %*% Bt
-    Bt = diag(pB*(1-psiB), nb) %*% Psi_bb %*% Bt
+    Kbq = Kbq + psiB %*% Bt
+    Bt = Psi_bb %*% pB %*% notpsiB %*% Bt
   }
   model$KGV$Kbq = Kbq
   return(model)
@@ -116,31 +160,43 @@ compute_Kbq.BQS = function(model, Tmax = 200){with(model,with(Mpar,{
   nq = dim(q)[1]
   ns = dim(s)[1]
 
-  Mbb = pB*(1-sigb)*(1-psiB)*Psi_bb
-  success = pB*(1-sigb)*psiB*Psi_bb
-  Mbq = pB*psiB*Psi_bq
-  Mbs = pB*sigb*psiB*Psi_bs
-  Mqb = pQ*(1-sigf)*psiQ*Psi_qb
-  Mqq = pQ*(1-sigq)*(1-psiQ)*Psi_qq
-  Mqs = pQ*(sigf*psiQ + sigq*(1-psiQ))*Psi_qs
-  Msb = pS*psiS*Psi_sb
-  Msq = 0*t(Mqs)
-  Mss = (1-psiS)*Psi_ss
-
-  M1 = cbind(Mbb, 0*Mqb, Msb, 0*Mqb)
-  M2 = cbind(Mbq, Mqq, Msq, 0*Mqq)
-  M3 = cbind(Mbs, Mqs, Mss, 0*Mqs)
-  M4 = cbind(success, diag(1,nq), 0*Msq, diag(1,nq))
-  M = rbind(M1, M2, M3, M4)
-
-  cohort = Psi_bq%*%diag(pB,nb)
-  Cyes = diag(psiQ, nq)%*%cohort
-  Cno = diag(1-psiQ, nq)%*%cohort
-
-  Kt = rbind(0*Mbb, Cno, 0*Mbs, Cyes)
-  for(i in 1:Tmax) Kt = M%*%Kt
-  Kt[-c(1:(nb+nq+ns)),]
-  model$KGV$Kbq = Kt[nb+nq+ns+c(1:nq),]
+  sigF = diag(sigf, nq)
+  nsigF = diag(1-sigf, nq)
+  npsiB = diag(1-psiB, nb)
+  psiB = diag(psiB, nb)
+  pB = diag(pB, nb)
+  nsigB = diag(1-sigb, nb)
+  sigB = diag(sigb, nb)
+  npsiS = diag(1-psiS, ns)
+  psiS = diag(psiS, ns)
+  pS = diag(pS, ns)
+
+  cohort = diag(pQ,nq)
+  B = Psi_bq %*% nsigF %*% diag(pQ,nq)
+  S = Psi_sq %*% sigF %*% diag(pQ,nq)
+
+  Mbb = Psi_bb %*% nsigB %*% pB %*% npsiB
+  trap = psiB
+  Msb = Psi_sb %*% sigB %*% pB %*% npsiB
+  Mbs = Psi_bs %*% pS %*% psiS
+  Mss = Psi_ss %*% pS %*% npsiS
+
+  M1 = cbind(Mbb, Mbs, 0*Mbb)
+  M2 = cbind(Msb, Mss, 0*Msb)
+  M3 = cbind(trap, 0*Psi_bs, diag(1,nb))
+  M = rbind(M1, M2, M3)
+
+  nix = nb + ns + c(1:nb)
+  TBS = rbind(B, S, 0*B)
+
+  looking = 1
+  while(looking > 1e-6) {
+    TBS = M%*%TBS
+    looking = sum(TBS[-nix,-nix])
+  }
+  nix = 1:(nb+ns)
+  Kbq = TBS[-nix,]
+  model$KGV$Kbq = Kbq
   return(model)
 }))}