tryna again

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

DESCRIPTION

@@ -0,0 +1,21 @@
+Package: ramp.falciparum
+Title: A Probabilistic Synthesis of Plasmodium falciparum Infection and Immunity
+Version: 0.1.0
+Authors@R: c(
+      person("John M.", "Henry", email = "[email protected]", role = c("aut","cre"), comment = c(ORCID = "0000-0001-8684-2003")),
+       person("David L.", "Smith", email = "[email protected]", role = c("aut"), comment = c(ORCID = "0000-0003-4367-3849"))
+    )
+Description: A Probabilistic Synthesis of Plasmodium falciparum Infection and Immunity, a part of RAMP     
+License: GPL (>= 3)
+Encoding: UTF-8
+RoxygenNote: 7.3.1
+URL: https://dd-harp.github.io/ramp.falciparum/
+LazyData: true
+Imports: 
+  stats,
+  deSolve
+Suggests: 
+  knitr, 
+  graphics, 
+  viridisLite
+VignetteBuilder: knitr

NAMESPACE

@@ -0,0 +1,119 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(Fmu,W)
+S3method(Fmu,base)
+S3method(Fmu,chronic)
+S3method(Wda,delta)
+S3method(Wda,none)
+S3method(ageFoI,flat)
+S3method(ageFoI,type2)
+S3method(d_Omega,beta)
+S3method(d_counts,nb)
+S3method(d_counts,pois)
+S3method(d_detect,nb)
+S3method(d_detect,pois)
+S3method(log10RBC,static)
+S3method(p_Omega,beta)
+S3method(p_counts,nb)
+S3method(p_counts,pois)
+S3method(q_Omega,beta)
+S3method(r_Omega,beta)
+S3method(seasonalFoI,exp)
+S3method(seasonalFoI,flat)
+S3method(seasonalFoI,sin)
+S3method(sigma_mu,abc)
+S3method(trendFoI,flat)
+export(FQ)
+export(Fmu)
+export(FoI)
+export(MoIDistPlot)
+export(Wda)
+export(ageFoI)
+export(cdfConvolve2)
+export(cdfConvolve2a)
+export(cdfConvolve2b)
+export(dAoI)
+export(dAoIda)
+export(dAoY)
+export(dAoYN)
+export(dAoYda)
+export(dDensityPaConvolve2)
+export(dMoIda)
+export(d_Omega)
+export(d_alpha2density)
+export(d_clone_counts)
+export(d_clone_density)
+export(d_counts)
+export(d_detect)
+export(d_detect_mesh)
+export(d_moi_count)
+export(d_nz_counts_log)
+export(d_nz_counts_log_binned)
+export(d_parasite_counts)
+export(d_parasite_density)
+export(d_parasite_detect)
+export(d_parasite_detect_moi)
+export(dbeta1)
+export(dmda)
+export(dpda)
+export(log10RBC)
+export(meanMoI)
+export(mean_clone_counts)
+export(mean_parasite_counts)
+export(momentAoI)
+export(momentAoY)
+export(moments_clone_density)
+export(moments_parasite_density)
+export(nzPois)
+export(pAoI)
+export(pAoY)
+export(pAoYN)
+export(pAoY_long)
+export(p_Omega)
+export(p_alpha2density)
+export(p_clone_counts)
+export(p_counts)
+export(p_nz_counts_log)
+export(p_nz_counts_log_binned)
+export(p_parasite_counts)
+export(p_parasite_density)
+export(par_Fmu_W)
+export(par_Fmu_base)
+export(par_Fmu_chronic)
+export(par_Omega_beta)
+export(par_Wda_delta)
+export(par_Wda_none)
+export(par_expSeason)
+export(par_flatAge)
+export(par_flatSeason)
+export(par_flatTrend)
+export(par_lRBC_static)
+export(par_nb)
+export(par_pois)
+export(par_sigma_abc)
+export(par_sinSeason)
+export(par_type2Age)
+export(parasite_density)
+export(pbeta1)
+export(q_Omega)
+export(q_alpha2density)
+export(qbeta1)
+export(rAoI)
+export(rAoY)
+export(rAoYN)
+export(rRda)
+export(r_Omega)
+export(r_alpha2density)
+export(r_parasite_density)
+export(rbeta1)
+export(sFmu)
+export(seasonalFoI)
+export(sigma_mu)
+export(solveMMinfty)
+export(solve_dAoI)
+export(solve_dAoYda)
+export(solve_dm)
+export(solve_dpda)
+export(trendFoI)
+export(truePR)
+export(zda)

R/Fmu-base.R

@@ -0,0 +1,39 @@
+
+#' Compute expected log10 parasite densities, `mu`, as a function of the age of infection `alpha`
+#' @description Compute mean log10 parasites for a model with no immunity
+#'
+#' @inheritParams Fmu
+#'
+#' @return mean log10 parasite densities
+#' @export
+#'
+Fmu.base = function(alpha, W, par_Fmu){with(par_Fmu,{
+  B = tildel + (tildeb-tildel)*exp(-(Sa*(alpha-peak)))
+  ix = which(alpha<peak)
+  if(length(ix>0)) B[ix] = tildel+(tildeb-tildel)*alpha[ix]/peak
+  #  ix = which(alpha<liver)
+  #  if(length(ix>0)) B[ix] = sqrt(-1)
+  B
+})}
+
+#' Set up parameters for [Fmu.base]
+#'
+#' @param peak The age of infection (in days) when parasite densities peak
+#' @param liver The age of infection (in days) when parasites emerge from the liver
+#' @param tildeb The maximum expected log10 parasite densities
+#' @param tildel The minimum expected log10 parasite densities
+#' @param Sa The decline in mu with respect to alpha
+#'
+#' @return a [list]
+#' @export
+#'
+par_Fmu_base = function(peak=20, liver=7, tildeb=10.3, tildel=2, Sa=0.0033){
+  par = list()
+  class(par) <- "base"
+  par$peak=peak
+  par$liver=liver
+  par$tildeb=tildeb
+  par$tildel=tildel
+  par$Sa=Sa
+  par
+}

R/Fmu-chronic.R

@@ -0,0 +1,29 @@
+#' Compute expected log10 parasite densities, `mu`, as a function of the age of infection `alpha`
+#' @description Compute mean log10 parasites for a model with no immunity
+#'
+#' @inheritParams Fmu
+#'
+#' @return mean log10 parasite densities
+#' @export
+#'
+Fmu.chronic = function(alpha, W, par_Fmu){with(par_Fmu,{
+  B = tildel + (tildeb-tildel)*exp(-(Sa*(alpha)))
+})}
+
+#' Set up parameters for [Fmu.chronic]
+#'
+#' @param tildeb The maximum expected log10 parasite densities
+#' @param tildel The minimum expected log10 parasite densities
+#' @param Sa The decline in mu with respect to alpha
+#'
+#' @return a [list]
+#' @export
+#'
+par_Fmu_chronic = function(tildeb=10.3, tildel=2, Sa=0.0033){
+  par = list()
+  class(par) <- "chronic"
+  par$tildeb=tildeb
+  par$tildel=tildel
+  par$Sa=Sa
+  par
+}

R/Fmu-immune.R

@@ -0,0 +1,41 @@
+
+#' Compute expected log10 parasite densities, `mu`, as a function of the age of infection `alpha`
+#' @description Compute mean log10 parasites for a model with no immunity
+#'
+#' @inheritParams Fmu
+#'
+#' @return mean log10 parasite densities
+#' @export
+#'
+Fmu.W = function(alpha, W, par_Fmu){with(par_Fmu,{
+  B = tildel + (tildeb-tildel)*exp(-Sa*(alpha-peak)-Sw*W)
+  ix = which(alpha<=peak)
+  if(length(ix>0)) B[ix] = tildel+(tildeb-tildel)*exp(-Sw*W)*alpha[ix]/peak
+  # ix = which(alpha<=liver)
+  #  if(length(ix>0)) B[ix] = sqrt(-1)
+  B
+})}
+
+#' Set up parameters for [Fmu.W]
+#'
+#' @param peak The age of infection (in days) when parasite densities peak
+#' @param liver The age of infection (in days) when parasites emerge from the liver
+#' @param tildeb The maximum expected log10 parasite densities
+#' @param tildel The minimum expected log10 parasite densities
+#' @param Sa The decline in mu with respect to alpha
+#' @param Sw The decline in mu with respect to immunity
+#'
+#' @return mean log10 parasite densities
+#' @export
+#'
+par_Fmu_W = function(peak=20, liver=7, tildeb=10.3, tildel=2, Sa=0.0033, Sw=0.001){
+  par = list()
+  class(par) <- "W"
+  par$peak=peak
+  par$liver=liver
+  par$tildeb=tildeb
+  par$tildel=tildel
+  par$Sa=Sa
+  par$Sw=Sw
+  par
+}

R/Fmu-interface.R

@@ -0,0 +1,27 @@
+#' Compute mean, expected parasite densities `mu` as a function of the age of infection `alpha`
+#'
+#' @param alpha the age of a parasite infection
+#' @param W the immune tracking variables
+#' @param par_Fmu a [list] that defines a model
+#'
+#' @return mean log10 parasite densities, a [numeric] vector of length(alpha)
+#' @export
+#'
+Fmu = function(alpha, W, par_Fmu){
+  UseMethod("Fmu", par_Fmu)
+}
+
+#' Compute mean, expected parasite densities `mu` for multiple values of the age of infection `alpha`
+#'
+#' @param alpha a parasite's age of infection
+#' @param W immune tracking variables
+#' @param par_Fmu a [list] that defines a model
+#'
+#' @return mean log10 parasite densities, a [numeric] vector of length(alpha)
+#' @export
+#'
+sFmu = function(alpha, W, par_Fmu){
+  sapply(alpha, Fmu, W=W, par_Fmu=par_Fmu)
+}
+
+

R/MMinfinity.R

@@ -0,0 +1,81 @@
+
+#' Derivatives for the queuing model \eqn{M/M/\infty}
+#'
+#' @description
+#'
+#' This queuing model \eqn{M/M/\infty} tracks the MoI in a cohort of humans
+#' as it ages. It assumes a time- and age-dependent hazard rate for infection,
+#' called the force of infection (FoI, \eqn{h_\tau(a)}). Infections do not affect
+#' each other, and each one clears independently at the rate \eqn{r}.
+#'
+#' Let \eqn{\zeta_i} the fraction of the population with
+#' MoI = i, then
+#' \deqn{\frac{d\zeta_0}{da}= -h_\tau(a) \zeta_0 + r \zeta_1}
+#' and for \eqn{i\geq 1}
+#' \deqn{\frac{d\zeta_i}{da}= h_\tau(a) \left( \zeta_{i-1} - \zeta_i \right)  - ri \zeta_i + r(i+1)\zeta_{i+1}}
+#'
+#' This function computes the derivatives in a form that can be used by [deSolve::ode].
+#'
+#' @param a the host age
+#' @param M the state variables
+#' @param pars the parameters
+#' @param FoIpar \eqn{h_\tau(a)}, a [list] formatted to compute [FoI]
+#'
+#' @return the derivatives as a [list]
+#' @seealso [solveMMinfty]
+#' @export
+#'
+dMoIda = function(a, M, pars, FoIpar){with(as.list(c(M,pars)),{
+  foi = h*FoI(a, FoIpar, tau)
+  i = 1:N
+  m = i-1
+  dM = 0*M-(foi + r*m)*M
+  dM[-N] = dM[-N] + r*m[-1]*M[-1]
+  dM[-1] = dM[-1] + foi*M[-N]
+  list(c(dM))
+})}
+
+#' Solve the queuing model \eqn{M/M/\infty}
+#'
+#' @description
+#'
+#' A wrapper to solve the queuing model \eqn{M/M/\infty} (see [dMoIda]).
+#'
+#' The function automatically sets the maximum MoI to be computed, and it
+#' sets initial conditions. The equations are solved using [deSolve::ode] and returned at
+#' regular intervals dt from age 0 up to Amax (in days).
+#'
+#' @param h the force of infection
+#' @param FoIpar \eqn{h_\tau(a)}, a [list] formatted to compute [FoI]
+#' @param r the clearance rate for a simple infection
+#' @param tau the cohort birthday
+#' @param Amax The maximum runtime (in days)
+#' @param dt The output frequency (in days)
+#'
+#' @return a [list] with the orbits by name
+#' @seealso [dMoIda]
+#' @export
+#'
+solveMMinfty = function(h, FoIpar, r=1/200, tau=0, Amax=730, dt=1){
+  tms = seq(0, Amax, by = dt)
+  N = round(max(4*h/r,20))
+  prms = c(h=h,r=r,N=N,tau=tau)
+  inits = rep(0,N)
+  inits[1]=1
+  out = deSolve::ode(inits, times=tms, dMoIda, prms, FoIpar=FoIpar)
+  time = out[,1]; moi = out[,-1]
+  m = moi %*% c(0:(N-1))
+  list(time=time, moi=moi, m=m)
+}
+
+#' Plot the output of [solveMMinfty]
+#'
+#' @param moi the mean moi
+#' @param t the time
+#' @param clr1 the color
+#' @export
+MoIDistPlot = function(moi, t, clr1 = "red"){
+  N = dim(moi)[2]-2
+  mm = 1:N -1
+  plot(mm, moi[t,1:N +1], type="h", xlab = "MoI", ylab = expression(M[tau](a)), main = paste ("Age = ", t, "Days"))
+}