library(ggplot2)

t = seq(1,40)
realcases = 100*exp(t/4)
realincrement = diff(c(0, realcases))


testseekers = 10e3 + realincrement*4
maxtests    = 20000

## now assume that you test *up to* 20k people. if more people are
## seeking tests, you test a random subset of the seekers
## getting a binomial count of positives for the given frequency

ntests = rep(0,NROW(t));
ntests[1] = 10e3;
confinc = rep(0,NROW(t));
confinc[1] = 100;
for(i in 2:(NROW(t)-1)){
    if(testseekers[i] < maxtests){
        confinc[i] = realincrement[i]
        ntests[i] = testseekers[i]
    }
    else if(testseekers[i] > maxtests){
        confinc[i] = min(realincrement[i],rbinom(1, maxtests, realincrement[i]/testseekers[i]))
        ntests[i] = maxtests
    }
}

confinc = head(confinc, - 1)/1000
ntests  = head(ntests, -1)/1000

# cumconf = cumsum(confinc)
# cumtests = cumsum(ntests)

# ggplot(data.frame(t=t,conf=cumconf,nt=cumtests,real=realcases))+geom_line(aes(t,cumconf),color="blue")  + geom_line(aes(t,nt),color="green")+ geom_line(aes(t,real),color="red") +coord_cartesian(xlim=c(0,35),ylim=c(0,400000));
# ggplot(data.frame(t=t,conf=cumconf,nt=cumtests,real=realcases))+geom_line(aes(t,log(cumconf)),color="blue") + geom_line(aes(t,log(nt)),color="green")+ geom_line(aes(t,log(real)),color="red") +coord_cartesian(xlim=c(0,30),ylim=c(0,log(400000)));


par(mfrow = c(2, 2))
plot(ntests,
     ylim = c(0, max(ntests)),
     type = "b",
     col  = "Blue",
     ylab = "# per day (Thousands)",
     panel.first = grid())
lines(confinc, type = "b", col = "Red")
legend("topleft",
       legend = c("Tests", "Cases"),
       col    = c("Blue", "Red"),
       lwd    = 2)

plot(confinc/ntests,
     type = "b",
     col  = "Blue",
     ylab = "% Positive",
     panel.first = grid())

plot(ntests, confinc,
     type = "p",
     col  = "Blue",
     xlab = "# Tests (Thousands)",
     ylab = "# Cases (Thousands)",
     panel.first = grid())