Figure8.Rmd

---
title: "Figure 8"
subtitle: "A Probabilistic Synthesis of Malaria Epidemiology: Exposure, Infection, Parasite Densities, and Detection"
date: "`r format(Sys.time(), '%B %d, %Y')`"
author: John M. Henry, Austin R. Carter, Sean L. Wu, and David L. Smith
output:
  html_document 
---

***

[Home](Memory.html) | 
[Fig 4](Figure4.html) |
[Fig 5](Figure5.html) |
[Fig 6](Figure6.html) |
[Fig 7](Figure7.html) |
[Fig 8](Figure8.html) |
[Fig 9](Figure9.html) |
[Fig 10](Figure10.html)

***

# {.tabset}

## $\odot$

![Figure8.png](Figure8.png)

## Setup 

```{r}
library(ramp.falciparum)
library(viridisLite)
library(knitr)
```

```{r}
foiP3 = list(hbar = 1, 
             agePar = par_type2Age(), 
             seasonPar = par_sinSeason(), 
             trendPar = par_flatTrend())
```

## Figure 8 

```{r}
aa = seq(5, 5*365, by = 5)
```

This runs the code and saves it into a file.

```{r, eval=F}
meanP = sapply(aa, moments_clone_density, FoIpar=foiP3, hhat=5/365)
meanB = sapply(aa, moments_parasite_density, FoIpar=foiP3, hhat=5/365)
write.table(data.frame(P=meanP, B=meanB), "PB.txt")
```

To save time, we simply read it in.  

```{r}
PB = read.table("PB.txt")
meanP = PB$P
meanB = PB$B
PB1 = read.table("PB1.txt")
meanP1 = PB1$P
meanB1 = PB1$B
PB2 = read.table("PB2.txt")
meanP2 = PB2$P
meanB2 = PB2$B
```


```{r}
pFmu = par_Fmu_base()
solve_dAoYda(5/365, foiP3, Amax=5*365, dt=5) -> hybrid
tm = hybrid$time
approxP = Fmu(hybrid$x, 0, pFmu)
approxB = Fmu(hybrid$y, 0, pFmu)

plot(aa, meanB, type = "l", ylim = c(0, 13), col = "darkred", 
     xlab = "a - Cohort Age", 
     ylab = expression(list(xi, paste(log[10], " Parasite Densities"))))
#lines(tm, approxB, col = "salmon3", lwd=2)
lines(aa, meanP, type = "l", ylim = c(0, 13), col = "darkblue")
lines(tm, approxP, col = "cyan4")
mtext("Expected Densities Exactly vs. Hybrid Model Predictions", 3, 1, at=365)
```

## Remake

```{r purl Figure 8, eval=F}
print("Making Figure 8")
purl("Figure8.Rmd", "Figure8.R")
```


```{r source Figure 8, eval=F}
png("Figure8.png", height=360, width=540)
source("Figure8.R")
invisible(dev.off(dev.cur()))
```