diff --git a/.gitattributes b/.gitattributes
index 5ed691d..9c02166 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1,5 +1,5 @@
 eric-figs/design/ filter=lfs diff=lfs merge=lfs -text
 eric-figs/output/data/ filter=lfs diff=lfs merge=lfs -text
 eric-figs/packrat/bundles/* filter=lfs diff=lfs merge=lfs -text
-*.nb.html filter=lfs diff=lfs merge=lfs -text
-*.csv filter=lfs diff=lfs merge=lfs -text
+eric-figs/*.nb.html filter=lfs diff=lfs merge=lfs -text
+eric-figs/*.csv filter=lfs diff=lfs merge=lfs -text
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..3537ae5
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,39 @@
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+/*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# Packrat
+eric-figs/packrat/lib*/
+eric-figs/packrat/src
+tobyfigs/packrat/lib*/
+tobyfigs/packrat/src
diff --git a/simulation/.dockerignore b/simulation/.dockerignore
new file mode 100644
index 0000000..38d9849
--- /dev/null
+++ b/simulation/.dockerignore
@@ -0,0 +1,4 @@
+.git
+./data/*
+./src/*
+./work*/*
diff --git a/simulation/.gitignore b/simulation/.gitignore
new file mode 100644
index 0000000..daa8298
--- /dev/null
+++ b/simulation/.gitignore
@@ -0,0 +1,4 @@
+*~
+.Rhistory
+work*
+spaero_0.1.0*
diff --git a/simulation/Dockerfile b/simulation/Dockerfile
new file mode 100644
index 0000000..f34e5e0
--- /dev/null
+++ b/simulation/Dockerfile
@@ -0,0 +1,48 @@
+FROM r-base
+MAINTAINER Eamon O'Dea <[last name without apostrophe]35@gmail.com>
+
+RUN apt-get update && apt-get install -y -q --no-install-recommends \
+  openssh-server \
+  r-cran-coda \
+  r-cran-curl \
+  r-cran-desolve \
+  r-cran-doparallel \
+  r-cran-evaluate \
+  r-cran-foreach \
+  r-cran-formatr \
+  r-cran-ggplot2 \
+  r-cran-httr \
+  r-cran-igraph \
+  r-cran-iterators \
+  r-cran-jsonlite \
+  r-cran-lazyeval \
+  r-cran-memoise \
+  r-cran-mime \
+  r-cran-mvtnorm \
+  r-cran-nloptr \
+  r-cran-openssl \
+  r-cran-praise \
+  r-cran-rcpp \
+  r-cran-r6 \
+  r-cran-testthat \
+  r-cran-tgp \
+  r-cran-xtable \
+  r-cran-yaml
+
+RUN install2.r --error \
+  earlywarnings \
+  knitr \
+  lintr \
+  np \
+  pomp \
+  rmarkdown
+
+COPY ./spaero_0.2.0.1.tar.gz /home/docker/spaero_0.2.0.1.tar.gz
+RUN /usr/bin/R CMD check --no-build-vignettes --no-manual /home/docker/spaero_*.tar.gz > /home/docker/check.so \
+&& /usr/bin/R CMD INSTALL --build /home/docker/spaero_*.tar.gz
+
+RUN mkdir /var/run/sshd && echo 'docker:screencast' | chpasswd \
+  && sed 's@session\s*required\s*pam_loginuid.so@session optional pam_loginuid.so@g' -i /etc/pam.d/sshd
+EXPOSE 22
+
+CMD ["/usr/sbin/sshd", "-D"]
diff --git a/simulation/README.md b/simulation/README.md
new file mode 100644
index 0000000..df70c8e
--- /dev/null
+++ b/simulation/README.md
@@ -0,0 +1,18 @@
+The ``simulation`` subdirectory contains:
+1) simulation results used to generate the figures in the paper contained in the zip archive
+``work-2017-08-29_17-29-29-light.zip``. 
+2) All scripts used to generate the simulation results in the directory ``src``.
+3) Docker files (see below)
+4) A copy of the R package spaero
+5) Some additional non-essential but related scripts in ``src``.
+
+To facilitate reproducibility, the scripts were run within a Docker container based on an image that is available on Docker Hub. The file ``build-image`` was used to create the image based on
+the ``Dockerfile``. To reproduce the results on systems with Docker installed, run the bash script ``run-scripts`` the ``simulation`` subdirectory. This script will:
+
+1) Pull the Docker container down from Docker Hub and create a time-stamped child directory 
+2) Run the Makefile in ``src`` to reproduce the full simulation results in the time-stamped child directory
+
+Reproducing the full output is computationally intensive and not recommended on a desktop computer. The full output comes to about 30GB and includes many intermediate results which are not needed to recreate the figures. The script ``src/make-light.sh`` is included to generate a lighter compressed zip version of the results and only includes the data necessary for figures.
+
+It is also possible to run the scripts without Docker, although one might have to consult the Docker image to determine the correct software versions to use. One important dependency, the spaero
+R package for calculating the moving window statistics, is included in this subdirectory. The right sequence in which the scripts must be run to produce the results can be seen in ``src/Makefile``.
diff --git a/simulation/build-image b/simulation/build-image
new file mode 100755
index 0000000..a6ec7e6
--- /dev/null
+++ b/simulation/build-image
@@ -0,0 +1 @@
+docker build -t eamon/2016feasibility .
diff --git a/simulation/run-scripts b/simulation/run-scripts
new file mode 100755
index 0000000..35e48c0
--- /dev/null
+++ b/simulation/run-scripts
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+set -e
+
+topd=$(pwd)
+tstamp=$(date +%Y-%m-%d_%H-%M-%S)
+wkdir="${topd}/work-${tstamp}"
+ddir="${topd}/data"
+sdir="${topd}/src"
+hdir="/home/docker"
+
+mkdir ${wkdir}
+chmod 757 ${wkdir}
+docker run -v ${ddir}:${hdir}/data:ro \
+       -v ${sdir}:${hdir}/src:ro \
+       -v ${wkdir}:${hdir}/work \
+       --user docker -w ${hdir}/work \
+       eamon/2016feasibility:v20170829 /bin/bash -c "make -f ../src/Makefile"
+chown -R $(id -u):$(id -g) ${wkdir}
diff --git a/simulation/spaero_0.2.0.1.tar.gz b/simulation/spaero_0.2.0.1.tar.gz
new file mode 100644
index 0000000..6c5cd90
Binary files /dev/null and b/simulation/spaero_0.2.0.1.tar.gz differ
diff --git a/simulation/src/Makefile b/simulation/src/Makefile
new file mode 100644
index 0000000..320bead
--- /dev/null
+++ b/simulation/src/Makefile
@@ -0,0 +1,16 @@
+RS = '/usr/bin/Rscript'
+WD = /home/docker/work
+SD = /home/docker/src
+DD = /home/docker/data
+VPATH = $(SD):$(WD)
+
+all: checkpoint-03.rda
+
+checkpoint-03.rda: analyze_observations.R checkpoint-02.rda
+	$< >$(<F).so 2>$(<F).se
+
+checkpoint-02.rda: simulate_observation.R checkpoint-01.rda
+	$< >$(<F).so 2>$(<F).se
+
+checkpoint-01.rda: simulate_cases.R
+	$< >$(<F).so 2>$(<F).se
diff --git a/simulation/src/analyze_observations.R b/simulation/src/analyze_observations.R
new file mode 100755
index 0000000..f65f4d1
--- /dev/null
+++ b/simulation/src/analyze_observations.R
@@ -0,0 +1,94 @@
+#!/usr/bin/Rscript
+
+library(foreach)
+
+doParallel::registerDoParallel(cores=5)
+RNGkind("L'Ecuyer-CMRG")
+set.seed(1234)
+parallel::mc.reset.stream()
+
+load("checkpoint-02.rda")
+
+levs <- list()
+levs$bandwidth <- c(35, 100)
+levs$lag <- 1
+analysis_des_mat <- do.call(expand.grid, levs)
+
+calc_auc <- function(predictions, is_null){
+  test <- !is.na(predictions)
+  predictions <- predictions[test]
+  is_null <- is_null[test]
+  r <- rank(predictions)
+  r1 <- sum(r[!is_null])
+  n1 <- sum(!is_null)
+  n2 <- sum(is_null)
+  (r1 - n1 * (n1 + 1) / 2) / (n1 * n2)
+}
+
+analyze_observations <- function(bandwidth, lag, obs){
+  get_tau <- function(stat_ts, time){
+    stats::cor(stat_ts, time, use="pairwise.complete.obs", method="kendall")
+  }
+  get_ews_cor <- function(reports, bw=bandwidth, l=lag){
+      stats <- spaero::get_stats(reports,
+                                 center_bandwidth=bw, stat_bandwidth=bw,
+                                 center_trend="local_constant",
+                                 stat_trend="local_constant",
+                                 center_kernel="uniform",
+                                 stat_kernel="uniform", lag=l)$stats
+      taus <- sapply(stats, get_tau, time=seq_along(reports))
+      list(stats=stats, taus=taus)
+  }
+  loop_manip <- function(x){
+    loop_omult <- function(x) lapply(x, get_ews_cor)
+    sc <- lapply(x, loop_omult)
+    taus <- lapply(sc, function(x) sapply(x, "[[", "taus"))
+    taus <- do.call(cbind, taus)
+    list(sc=sc, taus=taus)
+  }
+  sct <- lapply(obs, loop_manip)
+  is_null <- c(rep(FALSE, ncol(sct$test$taus)), rep(TRUE, ncol(sct$null$taus)))
+  preds <- cbind(sct$test$taus, sct$null$taus)
+
+  df <- data.frame(is_null, t(preds), row.names=1:ncol(preds))
+  samp_stat <- function(x, w) {
+    apply(x[w, -1], 2, calc_auc, is_null=df$is_null[w])
+  }
+  bs <- boot::boot(data=df, statistic=samp_stat, R=300)
+  bssd <- apply(bs$t, 2, sd, na.rm=TRUE)
+  list(sct=sct, auc=bs$t0, auc_stderr=bssd)
+}
+
+analyzed_observations <-
+  foreach (i=seq(1, nrow(process_des_mat)),
+           .options.multicore=list(set.seed=TRUE, preschedule=FALSE)) %:%
+    foreach (j=seq(1, nrow(observation_des_mat)),
+             .options.multicore=list(set.seed=TRUE, preschedule=FALSE)) %:%
+      foreach (m=seq(1, nrow(analysis_des_mat)),
+               .options.multicore=list(set.seed=TRUE, preschedule=FALSE)) %dopar% {
+                 do.call(analyze_observations,
+                         c(analysis_des_mat[m, ],
+                           list(obs=simulated_observations[[i]][[j]])))
+    }
+warnings()
+
+res <- list()
+n <- 1
+for (i in seq(1, nrow(process_des_mat))){
+  pvars <- process_des_mat[i, ]
+  for (j in seq(1, nrow(observation_des_mat))){
+    ovars <- observation_des_mat[j, ]
+    for (m in seq(1, nrow(analysis_des_mat))){
+      auc <- analyzed_observations[[i]][[j]][[m]]$auc
+      names(auc) <- paste("AUC", names(auc), sep="_")
+      aucse <- analyzed_observations[[i]][[j]][[m]]$auc_stderr
+      names(aucse) <- paste(names(auc), "stderr", sep="_")
+      avars <- analysis_des_mat[m, ]
+      res[[n]] <- data.frame(c(pvars, ovars, avars, auc, aucse))
+      n <- n + 1
+    }
+  }
+}
+res <- do.call(rbind, res)
+
+save.image(file="checkpoint-03.rda")
diff --git a/simulation/src/extract-observations.R b/simulation/src/extract-observations.R
new file mode 100755
index 0000000..0b9e2e9
--- /dev/null
+++ b/simulation/src/extract-observations.R
@@ -0,0 +1,43 @@
+#!/usr/bin/Rscript
+
+load("checkpoint-02.rda")
+set.seed(1)
+
+res <- list()
+for (i in seq(1, nrow(process_des_mat))){
+  pvars <- process_des_mat[i, ]
+  jseq <- sample.int(nrow(observation_des_mat), 5) ## to reduce the size of the data
+##  for (j in seq(1, nrow(observation_des_mat))){
+  for (j in jseq){
+    ovars <- observation_des_mat[j, ]
+    sim_obs <- simulated_observations[[i]][[j]]
+    stopifnot(length(sim_obs$null) == 1L) ## To simplify code, we
+                                          ## assume that there is only
+                                          ## one realization of the
+                                          ## observation model per
+                                          ## realization of the
+                                          ## process model.
+    stopifnot(length(sim_obs$test) == 1L)
+    tdata <- sim_obs$test[[1L]]
+    ndata <- sim_obs$null[[1L]]
+    make_df <- function(dat, repid, is_emergence){
+      t <- seq_along(dat)
+      data.frame(c(pvars, ovars), process_replicate_id = repid,
+                 is_emergence = is_emergence,
+                 cbind(time_index = t, reported_cases = dat))
+    }
+    tdf <- Map(make_df, dat = tdata, repid = seq_along(tdata),
+               is_emergence = TRUE)
+    tdn <- Map(make_df, dat = ndata, repid = seq_along(ndata),
+               is_emergence = FALSE)
+    res <- c(res, tdf, tdn)
+  }
+}
+
+reports_ts <- do.call(rbind, res)
+## Delete potentially confusing column
+reports_ts$scenario <- NULL
+rownames(reports_ts) <- NULL
+
+write.csv(reports_ts, file="reports_ts.csv")
+save(reports_ts, file = "reports_ts.RData")
diff --git a/simulation/src/get-cases.R b/simulation/src/get-cases.R
new file mode 100755
index 0000000..8e94818
--- /dev/null
+++ b/simulation/src/get-cases.R
@@ -0,0 +1,27 @@
+#!/usr/bin/Rscript
+
+load("checkpoint-01.rda")
+
+## We'll delete observations that do not occur at multiples of
+## ``snapshots`` to save space because they are not used in any
+## further analysis
+
+test_snap <- function(x, period){
+  mod <- x %% period
+  abs(mod) < .Machine$double.eps^.5
+}
+
+get_snaps <- function(df) {
+  tms <- df$time
+  snapshots <- c(7, 30)
+  tests <- lapply(snapshots, function(per) test_snap(tms, period = per))
+  is_snap <- colSums(do.call(rbind, tests)) > 0
+  df[is_snap, ]
+}
+
+for(i in seq_len(nrow(process_des_mat))){
+  simulated_procs[[i]]$test <- get_snaps(simulated_procs[[i]]$test)
+  simulated_procs[[i]]$null <- get_snaps(simulated_procs[[i]]$null)
+}
+
+save(simulated_procs, process_des_mat, file="cases.RData")
diff --git a/simulation/src/get-res.R b/simulation/src/get-res.R
new file mode 100755
index 0000000..2e19244
--- /dev/null
+++ b/simulation/src/get-res.R
@@ -0,0 +1,3 @@
+#!/usr/bin/env Rscript
+load("checkpoint-03.rda")
+save(res, file = "res.RData")
diff --git a/simulation/src/make-light.sh b/simulation/src/make-light.sh
new file mode 100755
index 0000000..1658e81
--- /dev/null
+++ b/simulation/src/make-light.sh
@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+
+cd $outputdir
+../src/get-res.R
+../src/get-cases.R
+cd ..
+lightname="${outputdir}-light"
+zip -r $lightname ${outputdir}/ --exclude \*.rda
diff --git a/simulation/src/simulate_cases.R b/simulation/src/simulate_cases.R
new file mode 100755
index 0000000..e125eb5
--- /dev/null
+++ b/simulation/src/simulate_cases.R
@@ -0,0 +1,58 @@
+#!/usr/bin/Rscript
+
+library(foreach)
+
+doParallel::registerDoParallel(cores=30)
+RNGkind("L'Ecuyer-CMRG")
+set.seed(3234)
+parallel::mc.reset.stream()
+
+levs <- list()
+levs$external_forcing <- c(1 / 7)
+levs$host_lifetime <- c(70 * 365)
+
+
+levs$infectious_days <- c(30, 7)
+levs$observation_days <- c(20) * 365
+levs$population_size <- 10^c(6)
+levs$process_reps <- 1000
+
+levs$scenario <- c("emergence")
+process_des_mat <- do.call(expand.grid, levs)
+
+sample_process <- function(external_forcing, host_lifetime, infectious_days,
+                           observation_days, population_size, process_reps,
+                           scenario){
+    stopifnot(scenario == "emergence")
+
+    times <- seq(0, observation_days)
+    params <- c(gamma=1 / infectious_days, mu=1 / host_lifetime,
+                d=1 / host_lifetime, eta=external_forcing / population_size,
+                beta=0, rho=0.1, S_0=1, I_0=0, R_0=0, N_0=population_size)
+    covnul <- data.frame(gamma_t=c(0, 0), mu_t=c(0, 0), d_t=c(0, 0),
+                         eta_t=c(0, 0), beta_t=c(0, 0),
+                         time=c(0, observation_days))
+    beta_final <- (params["gamma"] + params["d"]) / population_size
+    covalt <- data.frame(gamma_t=c(0, 0), mu_t=c(0, 0), d_t=c(0, 0),
+                         eta_t=c(0, 0), beta_t=c(0, beta_final),
+                         time=c(0, observation_days))
+
+    simtest <- spaero::create_simulator(times=times, params=params, covar=covalt)
+    simnull <- spaero::create_simulator(times=times, params=params, covar=covnul)
+
+    ret <- list()
+    do_sim <- function(obj, nsim=process_reps){
+        cols_to_delete <- c("reports", "gamma_t", "mu_t", "d_t", "eta_t")
+        ret <- pomp::simulate(obj, nsim=nsim, as.data.frame=TRUE)
+        ret[, !colnames(ret) %in% cols_to_delete]
+    }
+    ret$test <- do_sim(simtest)
+    ret$null <- do_sim(simnull)
+    ret
+}
+
+simulated_procs <- foreach(i=seq(1, nrow(process_des_mat)),
+                           .options.multicore=list(set.seed=TRUE)) %dopar%
+  do.call(sample_process, process_des_mat[i, ])
+
+save.image(file="checkpoint-01.rda")
diff --git a/simulation/src/simulate_observation.R b/simulation/src/simulate_observation.R
new file mode 100755
index 0000000..0ea25c4
--- /dev/null
+++ b/simulation/src/simulate_observation.R
@@ -0,0 +1,47 @@
+#!/usr/bin/Rscript
+
+library(foreach)
+
+doParallel::registerDoParallel(cores=30)
+RNGkind("L'Ecuyer-CMRG")
+set.seed(3222)
+parallel::mc.reset.stream()
+
+load("checkpoint-01.rda")
+
+levs <- list()
+levs$aggregation_days <- c(7, 30)
+levs$neg_bin_k <- c(0.01, 1, 100)
+levs$observation_multiplier <- 1
+levs$reporting_prob <- 2^(seq(-8, 0, len=21))
+observation_des_mat <- do.call(expand.grid, levs)
+
+sample_observation <- function(aggregation_days, neg_bin_k,
+                               observation_multiplier, reporting_prob,
+                               procs){
+    splt <- function(x) {
+        split(x, x$sim)
+    }
+    aggr <- function(x, p=reporting_prob, m=observation_multiplier,
+                     k=neg_bin_k) {
+        tots <- aggregate.ts(x$cases, nfrequency=1 / aggregation_days)
+        mu <- tots * p
+        n <- length(mu)
+        rnbinom(n=n, mu=mu, size=k)
+    }
+    transform <- function(x, n=observation_multiplier){
+        sims <- splt(x)
+        replicate(n=n, lapply(sims, aggr), simplify=FALSE)
+    }
+    lapply(procs, transform)
+}
+
+simulated_observations <-
+  foreach (i=seq(1, nrow(process_des_mat)),
+           .options.multicore=list(set.seed=TRUE)) %:%
+    foreach (j=seq(1, nrow(observation_des_mat))) %dopar% {
+      do.call(sample_observation,
+              c(observation_des_mat[j, ], list(simulated_procs[[i]])))
+    }
+
+save.image(file="checkpoint-02.rda")
diff --git a/simulation/src/test-auc.R b/simulation/src/test-auc.R
new file mode 100644
index 0000000..d0ea529
--- /dev/null
+++ b/simulation/src/test-auc.R
@@ -0,0 +1,22 @@
+set.seed(234)
+
+calc_auc <- function(predictions, is_null){
+    r <- rank(predictions)
+    r1 <- sum(r[!is_null])
+    n1 <- sum(!is_null)
+    n2 <- sum(is_null)
+    (r1 - n1 * (n1 + 1) / 2) / (n1 * n2)
+}
+
+compare_calcs <- function(){
+  n <- runif(1, 10, 1000)
+  predictions <- rnorm(n)
+  is_null <- rep(TRUE, n)
+  is_null[seq(1, n / 2)] <- FALSE
+  us <- calc_auc(predictions, is_null)
+  pred <- ROCR::prediction(predictions, !is_null)
+  them <- ROCR::performance(pred, "auc")@y.values[[1]]
+  stopifnot(all.equal(us, them))
+}
+
+invisible(replicate(100, compare_calcs()))
diff --git a/simulation/work-2017-08-29_17-29-29-light.zip b/simulation/work-2017-08-29_17-29-29-light.zip
new file mode 100644
index 0000000..2412afe
Binary files /dev/null and b/simulation/work-2017-08-29_17-29-29-light.zip differ
diff --git a/tobyfigs/.Rprofile b/tobyfigs/.Rprofile
new file mode 100644
index 0000000..22ec147
--- /dev/null
+++ b/tobyfigs/.Rprofile
@@ -0,0 +1,3 @@
+#### -- Packrat Autoloader (version 0.4.8-1) -- ####
+source("packrat/init.R")
+#### -- End Packrat Autoloader -- ####
diff --git a/tobyfigs/.gitignore b/tobyfigs/.gitignore
new file mode 100644
index 0000000..d300c12
--- /dev/null
+++ b/tobyfigs/.gitignore
@@ -0,0 +1,2 @@
+packrat/lib*/
+packrat/src/*
diff --git a/tobyfigs/cases.RData b/tobyfigs/cases.RData
new file mode 100644
index 0000000..75c3afb
Binary files /dev/null and b/tobyfigs/cases.RData differ
diff --git a/tobyfigs/heat-plot.R b/tobyfigs/heat-plot.R
new file mode 100644
index 0000000..dc8e718
--- /dev/null
+++ b/tobyfigs/heat-plot.R
@@ -0,0 +1,214 @@
+#Script for heat plot 
+# 3-#2017-07-05 -- Toby
+
+library(colorspace)
+library(lattice)
+library(fields)
+library(dplyr)
+
+
+## Eric's AUC gradient ##
+# Given error of ~0.5, set number of levels to 10.  For a higher precision, use nlevels = 20 
+nlevels = 100  # number of bins = number of colors
+# Color scale with {colorspace}
+AUC.colors <- diverge_hcl(
+  n=nlevels,  # number of colors
+  h = c(45, 225),  # Hues (low, hi)
+  c = 100, # fixed Chroma or Chroma Range (edges, center)
+  l = c(90, 10),  # Lightness range (edges, center) 
+  power = 1  # exponent
+)
+#Gray colour used in paper:
+gray_colour <- rgb(0,0,0,.25)
+dl <- seq(0,1,0.01)
+
+## AUC data ##
+load("res.RData")
+res <- res[which(res$reporting_prob > 0.01),]
+
+
+heat_map_plot <- function(bw){
+    
+  ##Filter for bandwidth Choose either 35 or 100
+  res <- filter(res, bandwidth == bw)
+  
+  ## Results for 7 day recovery
+  res.7.wkly <- res[ which(res$infectious_days==7 
+                           & res$aggregation_days == 7), ]
+  
+  res.7.mon <- res[ which(res$infectious_days==7 
+                          & res$aggregation_days == 30), ]
+  
+  ## Results for 30 day recovery
+  res.30.wkly <- res[ which(res$infectious_days==30 
+                            & res$aggregation_days == 7), ]
+  
+  res.30.mon <- res[ which(res$infectious_days==30 
+                           & res$aggregation_days == 30), ]
+  
+  ## function to organize data for heat plots ##
+  org_data <- function(dat){
+    binLevs <- 3
+    ewsLevs <- ncol(dat[,14:23])
+    repLevs <- length(unique(dat$reporting_prob))
+    
+    tmpList <- list() # For each EWS
+    # Make each level of dispersion a new column
+    for(i in 1:ewsLevs){
+      cl <- 13 + i
+      tmp <- dat[, c(1:13, cl)]
+      tmp=filter(tmp, neg_bin_k!=.1)
+      tmp=filter(tmp, neg_bin_k!=10)
+      tmp <- arrange(tmp, neg_bin_k)
+      
+      tmp <- matrix(tmp[,14], ncol=repLevs, byrow = T)
+      row.names(tmp)=c("100",  "1" , ".01")
+      
+      colnames(tmp) <- round(unique(dat$reporting_prob), 5)
+      tmpList[[i]] <- tmp
+    }
+    
+    # format for heat plot
+    newDat <- do.call(rbind, tmpList)
+    newDat <- t(newDat)
+    
+    kurt=newDat[,c(28:30)]
+    skew=newDat[,c(25:27)]
+    cov =newDat[,c(22:24)]
+    ind =newDat[,c(19:21)]
+    men =newDat[,c(16:18)]
+    dec =newDat[,c(13:15)]
+    acf =newDat[,c(10:12)]
+    acv =newDat[,c(7:9)]
+    vdif =newDat[,c(4:6)]
+    var =newDat[,c(1:3)]
+    ord = cbind(men, var, vdif,
+                ind, acv, acf,
+                dec, cov, skew,
+                kurt)
+    return(ord)
+  }
+  
+  yLabs <- c("Mean", "Variance", "1st Diff. Var.", 
+             "Ind. of Dis.", "Autocovar.", "Autocorr.",
+             "Decay time", "Coeff. Var", "Skewness",
+             "Kurtosis")
+  
+  # For plotting
+  nRepProb <- length(unique(res$reporting_prob))
+  xVals <- seq(0,nRepProb-1)/(nRepProb-1)  # seq(0,20)/20
+  xLabVals <- xVals[seq(1, length(xVals), 6)]
+  xLabs <- round(unique(res$reporting_prob), 2)
+  xLabsNew <- xLabs[seq(1, length(xLabs), 6)]
+  yVals <- seq(2.5,50, by=5)/50
+  
+  addLines <- function(clr="white", lwd=1.5){
+    abline(h=18/200, col=clr, lwd=lwd)
+    abline(h=39/200, col=clr, lwd=lwd)
+    abline(h=59/200, col=clr, lwd=lwd)
+    abline(h=79.5/200, col=clr, lwd=lwd)
+    abline(h=100/200, col=clr, lwd=lwd)
+    abline(h=120/200, col=clr, lwd=lwd)
+    abline(h=141/200, col=clr, lwd=lwd)
+    abline(h=161/200, col=clr, lwd=lwd)
+    abline(h=181.5/200, col=clr, lwd=lwd)
+  }
+  
+  dispLabs3 <- function(){
+    #dispersion labels
+    axis(4, at = seq(46.5,50, length.out = 3)/50, 
+         labels = rep("", 3), tck=-.1, col.ticks = gray_colour,col = "white", pos = 1.05)
+    mtext(side=4, "less dispersed", line=.95, las=1, cex=.6, at=1.05, xpd=T)
+    mtext(side=4, "100", line=3, las=1, cex=.4, at=1, xpd=T)
+    mtext(side=4, " 1 ", line=3.1, las=1, cex=.4, at=.965, xpd=T)
+    mtext(side=4, "0.01", line=3, las=1, cex=.4, at=.93, xpd=T)
+    mtext(side=4, "more dispersed", line=.91, las=1, cex=.6, at=.88, xpd=T)
+  }
+  
+  addLines3 <- function(clr="white", lwd=1.5){
+    abline(h=17/200, col=clr, lwd=lwd)
+    abline(h=38/200, col=clr, lwd=lwd)
+    abline(h=59/200, col=clr, lwd=lwd)
+    abline(h=79.5/200, col=clr, lwd=lwd)
+    abline(h=100/200, col=clr, lwd=lwd)
+    abline(h=120/200, col=clr, lwd=lwd)
+    abline(h=141/200, col=clr, lwd=lwd)
+    abline(h=162/200, col=clr, lwd=lwd)
+    abline(h=182/200, col=clr, lwd=lwd)
+  }
+  
+  # Data for heat plots
+  H5 <- org_data(res.7.wkly)
+  H6 <- org_data(res.7.mon)
+  
+  H7 <- org_data(res.30.wkly)
+  H8 <- org_data(res.30.mon)
+  
+  
+  pdf(paste("heat-plot",bw,".pdf", sep=""), width=4.52,height=5.19)
+  # separate heat plots
+  layout(matrix(c(1, 2,
+                  3, 4,
+                  5, 5), nrow=3, byrow=TRUE),
+         heights = c(1,1,.25))
+  
+  par(oma=c(3,2,2,2), mar=c(1,5,3,3))
+  
+  image(H5, col=AUC.colors, xlab="", 
+        ylab="", axes=F, main="",
+        add.expr= abline(h=5/50, col="white"))
+  axis(3, at=xLabVals, labels=xLabsNew, col.ticks = gray_colour, col = "white", pos = 1.03)
+  axis(3, at=xVals, labels=FALSE, tck=-0.05, col.ticks = gray_colour, col = "white", pos = 1.03)
+  axis(2, at=yVals, labels=yLabs, las=1, col.ticks = gray_colour, col = "white", pos = -0.05)
+  mtext(side=3, text="Aggregated weekly", line=3)
+  mtext(side=2, expression(paste(gamma, "=1/7")), line=5.5)
+  addLines3()
+  dispLabs3()
+  
+  par(mar=c(1,4,3,4))
+  image(H6, col=AUC.colors, xlab="", 
+        ylab="", axes=F)
+  axis(3, at=xLabVals, labels=xLabsNew, col.ticks = gray_colour, col = "white", pos = 1.03)
+  axis(3, at=xVals, labels=FALSE, tck=-0.05, col.ticks = gray_colour, col = "white", pos = 1.03)
+  mtext(side=3, text="Aggregated monthly", line=3)
+  axis(4, at=yVals, labels=yLabs, las=1, col.ticks = gray_colour, col = "white", pos =  1.05)
+  addLines3()
+  axis(2, at = seq(46.5,50,length.out = 3)/50, 
+       labels = rep("", 3), tck=-.1,col.ticks = gray_colour, col = "white", pos = -0.05)
+  
+  
+  par(mar=c(4,5,0,3))
+  image(H7, col=AUC.colors, xlab="", 
+        ylab="", axes=F)
+  axis(1, at=xLabVals, labels=xLabsNew, col.ticks = gray_colour, col = "white", pos = -0.03)
+  axis(1, at=xVals, labels=FALSE, tck=-0.05, col.ticks = gray_colour, col = "white", pos = -0.03)
+  axis(2, at=yVals, labels=yLabs, las=1, col.ticks = gray_colour, col = "white",pos = -0.05)
+  mtext(side=2, expression(paste(gamma, "=1/30")), line=5.5)
+  mtext(side=1, "Reporting probability", line=3)
+  addLines3()
+  dispLabs3()
+  
+  
+  par(mar=c(4,4,0,4))
+  image(H8, col=AUC.colors, xlab="", 
+        ylab="", axes=F)
+  axis(1, at=xLabVals, labels=xLabsNew, col.ticks = gray_colour, col = "white", pos = -0.03)
+  axis(1, at=xVals, labels=FALSE, tck=-0.05, col.ticks = gray_colour, col = "white", pos = -0.03)
+  axis(4, at=yVals, labels=yLabs, las=1, col.ticks = gray_colour, col = "white", pos = 1.05)
+  mtext(side=1, "Reporting probability", line=3)
+  addLines3()
+  axis(2, at = seq(46.5,50, length.out = 3)/50, 
+       labels = rep("", 3), tck=-.1, col.ticks = gray_colour, col = "white", pos = -0.05)
+  
+  
+  par(mar=c(1.2,7,1.2,7)) 
+  legend_data <- matrix(dl, nrow = length(dl), ncol = 1)
+  image(legend_data, col=AUC.colors, xlab="", 
+              ylab="", axes=F, useRaster = T)
+  axis(1, at=seq(0.0,1,0.2), labels=seq(0.0,1,0.2), col.ticks = gray_colour, col = "white", pos = -1.5)
+  mtext(side=1, "AUC", line=3)
+  dev.off()
+}
+
+heat_map_plot(35)
+heat_map_plot(100)
diff --git a/tobyfigs/heat-plot100.pdf b/tobyfigs/heat-plot100.pdf
new file mode 100644
index 0000000..c772e6d
Binary files /dev/null and b/tobyfigs/heat-plot100.pdf differ
diff --git a/tobyfigs/heat-plot35.pdf b/tobyfigs/heat-plot35.pdf
new file mode 100644
index 0000000..19fe0b7
Binary files /dev/null and b/tobyfigs/heat-plot35.pdf differ
diff --git a/tobyfigs/high-low-plot.R b/tobyfigs/high-low-plot.R
new file mode 100644
index 0000000..22db66f
--- /dev/null
+++ b/tobyfigs/high-low-plot.R
@@ -0,0 +1,102 @@
+library(ggplot2)
+library(ggthemes) ##
+library(tidyverse)
+#library(plyr)
+library(reshape2)
+library(colorspace)
+
+## Eric's AUC gradient 
+# Color scale with {colorspace}
+nlevels = 100
+AUC_colors <- diverge_hcl(
+  n=nlevels,  # number of colors
+  h = c(45, 225),  # Hues (low, hi)
+  c = 100, # fixed Chroma or Chroma Range (edges, center)
+  l = c(90, 10),  # Lightness range (edges, center) 
+  power = 1  # exponent
+)
+names(AUC_colors) <- seq(0,1,1/(nlevels-1))
+
+## Load and select data to be displayed
+load("res.RData")
+res <- res[which(res$bandwidth==35),]
+res <- res[which(res$reporting_prob > 0.04),]
+rep_prob_lvls <- levels(as.factor(res$reporting_prob))[c(12,1)]
+neg_bin_k_lvls <- levels(as.factor(res$neg_bin_k))[c(3,1)]
+inf_days_lvls <- levels(as.factor(res$infectious_days))[1]
+agg_days_lvls <- levels(as.factor(res$aggregation_days))[1]
+
+hl_data <- res %>%
+  select_(.dots = names(.)[c(3,8,9,11,14:33)]) %>%
+  filter(neg_bin_k %in% neg_bin_k_lvls) %>%
+  filter(infectious_days == inf_days_lvls) %>%
+  filter(aggregation_days == agg_days_lvls) %>%
+  filter(reporting_prob %in% rep_prob_lvls) %>%
+  melt(id.vars = names(.)[1:4], value.name = "AUC") %>%
+  mutate(variable = as.factor(variable))
+
+nrows <- nrow(hl_data)/2
+
+hl_data <- hl_data %>%
+  mutate(err = lead(variable,nrows), AUC_err = lead(AUC,nrows)) %>%
+  filter(!is.na(err)) %>%
+  mutate(reporting_prob =  factor(reporting_prob, levels(as.factor(.$reporting_prob))[c(2,1)])) %>%
+  mutate(neg_bin_k = as.factor(neg_bin_k)) %>%
+  droplevels()
+new_levels <-  c("Mean", "Variance", "1st Diff. Var.", 
+                 "Index of Dis.", "Autocovar.", "Autocorr.",
+                 "Decay time", "Coeff. Var", "Skewness",
+                 "Kurtosis")
+hl_data$variable <- factor(hl_data$variable,levels(hl_data$variable)[c(6,1,2,7,3,4,5,8,9,10)])
+levels(hl_data$variable) <- new_levels
+levels(hl_data$reporting_prob) <- c("High reporting","Low reporting")
+levels(hl_data$neg_bin_k) <- c("High overdispersion", "Low overdispersion")
+
+#Define plotting function
+hl_plot <- function(hl_data){
+  ggplot(hl_data) + 
+    geom_bar(aes(x=variable,y=AUC-0.5, fill = AUC, color = AUC),stat="identity"  ) + 
+    facet_grid(reporting_prob~neg_bin_k) +
+    geom_rangeframe(colour ="black") +
+    scale_fill_gradientn(limits = c(0,1),colours=AUC_colors) +
+    scale_color_gradientn(limits = c(0,1),colours=AUC_colors) +
+    scale_y_continuous(name = "AUC",labels=c("0.0","0.25","0.5","0.75","1.0"))
+}
+
+
+#par(fg = "#FF6AD5", bg = "#FFFFFF", col.lab = "#94D0FF",
+#    col.main = "#C774E8", col.sub	= "#94D0FF",  col.axis = "#94D0FF",
+#    font.lab = 3)
+
+text_color <- "black"
+background_color <- "white"
+font_chosen <- "Helvetica"
+
+hl_plot(hl_data) +#theme_par() +
+  theme(text = element_text(color=text_color, size = 18, family=font_chosen),
+        title = element_text(size = 14),
+        line = element_line(color=text_color),
+        rect = element_rect(color="black"),
+        axis.title.x=element_blank(),
+        axis.title.y = element_text(color=text_color,face = "plain", 
+                                    family=font_chosen, size = 14),
+        axis.text.x = element_text(color = text_color, angle = 45,
+                                   vjust = 1, hjust=1),#family=font_chosen, size = 14),
+        axis.text.y = element_text(color = text_color, family=font_chosen),
+        axis.ticks = element_line(color=rgb(0,0,0,.25),line),
+        axis.line = element_blank(),
+        legend.text = element_text(color = text_color, family=font_chosen),
+        legend.title = element_text(color = text_color, family = font_chosen),
+        panel.grid.major = element_blank(),
+        panel.grid.major.y = element_line(color = rgb(0,0,0,.25)),
+        panel.grid.minor = element_blank(),
+        panel.border = element_blank(),
+        panel.background = element_rect(color=NA,fill=NA),
+        panel.spacing = unit(2, "lines"),
+        strip.text.y= element_text(color = text_color,angle = 270, family=font_chosen),
+        strip.text.x= element_text(color = text_color, family=font_chosen),
+        strip.background = element_rect(fill=NA),
+        plot.background = element_rect(color=NA, fill=background_color),
+        plot.margin = unit(c(1,1,1,1), "cm"))
+
+ggsave("high-low.pdf",width =1.2*8.64,height=1.2*4.20)
diff --git a/tobyfigs/high-low.pdf b/tobyfigs/high-low.pdf
new file mode 100644
index 0000000..7f1b192
Binary files /dev/null and b/tobyfigs/high-low.pdf differ
diff --git a/tobyfigs/packrat/init.R b/tobyfigs/packrat/init.R
new file mode 100644
index 0000000..a768be5
--- /dev/null
+++ b/tobyfigs/packrat/init.R
@@ -0,0 +1,217 @@
+local({
+
+  ## Helper function to get the path to the library directory for a
+  ## given packrat project.
+  getPackratLibDir <- function(projDir = NULL) {
+    path <- file.path("packrat", "lib", R.version$platform, getRversion())
+
+    if (!is.null(projDir)) {
+
+      ## Strip trailing slashes if necessary
+      projDir <- sub("/+$", "", projDir)
+
+      ## Only prepend path if different from current working dir
+      if (!identical(normalizePath(projDir), normalizePath(getwd())))
+        path <- file.path(projDir, path)
+    }
+
+    path
+  }
+
+  ## Ensure that we set the packrat library directory relative to the
+  ## project directory. Normally, this should be the working directory,
+  ## but we also use '.rs.getProjectDirectory()' if necessary (e.g. we're
+  ## rebuilding a project while within a separate directory)
+  libDir <- if (exists(".rs.getProjectDirectory"))
+    getPackratLibDir(.rs.getProjectDirectory())
+  else
+    getPackratLibDir()
+
+  ## Unload packrat in case it's loaded -- this ensures packrat _must_ be
+  ## loaded from the private library. Note that `requireNamespace` will
+  ## succeed if the package is already loaded, regardless of lib.loc!
+  if ("packrat" %in% loadedNamespaces())
+    try(unloadNamespace("packrat"), silent = TRUE)
+
+  if (suppressWarnings(requireNamespace("packrat", quietly = TRUE, lib.loc = libDir))) {
+
+    # Check 'print.banner.on.startup' -- when NA and RStudio, don't print
+    print.banner <- packrat::get_opts("print.banner.on.startup")
+    if (print.banner == "auto" && is.na(Sys.getenv("RSTUDIO", unset = NA))) {
+      print.banner <- TRUE
+    } else {
+      print.banner <- FALSE
+    }
+    return(packrat::on(print.banner = print.banner))
+  }
+
+  ## Escape hatch to allow RStudio to handle bootstrapping. This
+  ## enables RStudio to provide print output when automagically
+  ## restoring a project from a bundle on load.
+  if (!is.na(Sys.getenv("RSTUDIO", unset = NA)) &&
+      is.na(Sys.getenv("RSTUDIO_PACKRAT_BOOTSTRAP", unset = NA))) {
+    Sys.setenv("RSTUDIO_PACKRAT_BOOTSTRAP" = "1")
+    setHook("rstudio.sessionInit", function(...) {
+      # Ensure that, on sourcing 'packrat/init.R', we are
+      # within the project root directory
+      if (exists(".rs.getProjectDirectory")) {
+        owd <- getwd()
+        setwd(.rs.getProjectDirectory())
+        on.exit(setwd(owd), add = TRUE)
+      }
+      source("packrat/init.R")
+    })
+    return(invisible(NULL))
+  }
+
+  ## Bootstrapping -- only performed in interactive contexts,
+  ## or when explicitly asked for on the command line
+  if (interactive() || "--bootstrap-packrat" %in% commandArgs(TRUE)) {
+
+    message("Packrat is not installed in the local library -- ",
+            "attempting to bootstrap an installation...")
+
+    ## We need utils for the following to succeed -- there are calls to functions
+    ## in 'restore' that are contained within utils. utils gets loaded at the
+    ## end of start-up anyhow, so this should be fine
+    library("utils", character.only = TRUE)
+
+    ## Install packrat into local project library
+    packratSrcPath <- list.files(full.names = TRUE,
+                                 file.path("packrat", "src", "packrat")
+    )
+
+    ## No packrat tarballs available locally -- try some other means of installation
+    if (!length(packratSrcPath)) {
+
+      message("> No source tarball of packrat available locally")
+
+      ## There are no packrat sources available -- try using a version of
+      ## packrat installed in the user library to bootstrap
+      if (requireNamespace("packrat", quietly = TRUE) && packageVersion("packrat") >= "0.2.0.99") {
+        message("> Using user-library packrat (",
+                packageVersion("packrat"),
+                ") to bootstrap this project")
+      }
+
+      ## Couldn't find a user-local packrat -- try finding and using devtools
+      ## to install
+      else if (requireNamespace("devtools", quietly = TRUE)) {
+        message("> Attempting to use devtools::install_github to install ",
+                "a temporary version of packrat")
+        library(stats) ## for setNames
+        devtools::install_github("rstudio/packrat")
+      }
+
+      ## Try downloading packrat from CRAN if available
+      else if ("packrat" %in% rownames(available.packages())) {
+        message("> Installing packrat from CRAN")
+        install.packages("packrat")
+      }
+
+      ## Fail -- couldn't find an appropriate means of installing packrat
+      else {
+        stop("Could not automatically bootstrap packrat -- try running ",
+             "\"'install.packages('devtools'); devtools::install_github('rstudio/packrat')\"",
+             "and restarting R to bootstrap packrat.")
+      }
+
+      # Restore the project, unload the temporary packrat, and load the private packrat
+      packrat::restore(prompt = FALSE, restart = TRUE)
+
+      ## This code path only reached if we didn't restart earlier
+      unloadNamespace("packrat")
+      requireNamespace("packrat", lib.loc = libDir, quietly = TRUE)
+      return(packrat::on())
+
+    }
+
+    ## Multiple packrat tarballs available locally -- try to choose one
+    ## TODO: read lock file and infer most appropriate from there; low priority because
+    ## after bootstrapping packrat a restore should do the right thing
+    if (length(packratSrcPath) > 1) {
+      warning("Multiple versions of packrat available in the source directory;",
+              "using packrat source:\n- ", shQuote(packratSrcPath))
+      packratSrcPath <- packratSrcPath[[1]]
+    }
+
+
+    lib <- file.path("packrat", "lib", R.version$platform, getRversion())
+    if (!file.exists(lib)) {
+      dir.create(lib, recursive = TRUE)
+    }
+    lib <- normalizePath(lib, winslash = "/")
+
+    message("> Installing packrat into project private library:")
+    message("- ", shQuote(lib))
+
+    surround <- function(x, with) {
+      if (!length(x)) return(character())
+      paste0(with, x, with)
+    }
+
+    ## The following is performed because a regular install.packages call can fail
+    peq <- function(x, y) paste(x, y, sep = " = ")
+    installArgs <- c(
+      peq("pkgs", surround(packratSrcPath, with = "'")),
+      peq("lib", surround(lib, with = "'")),
+      peq("repos", "NULL"),
+      peq("type", surround("source", with = "'"))
+    )
+    installCmd <- paste(sep = "",
+                        "utils::install.packages(",
+                        paste(installArgs, collapse = ", "),
+                        ")")
+
+    fullCmd <- paste(
+      surround(file.path(R.home("bin"), "R"), with = "\""),
+      "--vanilla",
+      "--slave",
+      "-e",
+      surround(installCmd, with = "\"")
+    )
+    system(fullCmd)
+
+    ## Tag the installed packrat so we know it's managed by packrat
+    ## TODO: should this be taking information from the lockfile? this is a bit awkward
+    ## because we're taking an un-annotated packrat source tarball and simply assuming it's now
+    ## an 'installed from source' version
+
+    ## -- InstallAgent -- ##
+    installAgent <- 'InstallAgent: packrat 0.4.8-1'
+
+    ## -- InstallSource -- ##
+    installSource <- 'InstallSource: source'
+
+    packratDescPath <- file.path(lib, "packrat", "DESCRIPTION")
+    DESCRIPTION <- readLines(packratDescPath)
+    DESCRIPTION <- c(DESCRIPTION, installAgent, installSource)
+    cat(DESCRIPTION, file = packratDescPath, sep = "\n")
+
+    # Otherwise, continue on as normal
+    message("> Attaching packrat")
+    library("packrat", character.only = TRUE, lib.loc = lib)
+
+    message("> Restoring library")
+    restore(restart = FALSE)
+
+    # If the environment allows us to restart, do so with a call to restore
+    restart <- getOption("restart")
+    if (!is.null(restart)) {
+      message("> Packrat bootstrap successfully completed. ",
+              "Restarting R and entering packrat mode...")
+      return(restart())
+    }
+
+    # Callers (source-erers) can define this hidden variable to make sure we don't enter packrat mode
+    # Primarily useful for testing
+    if (!exists(".__DONT_ENTER_PACKRAT_MODE__.") && interactive()) {
+      message("> Packrat bootstrap successfully completed. Entering packrat mode...")
+      packrat::on()
+    }
+
+    Sys.unsetenv("RSTUDIO_PACKRAT_BOOTSTRAP")
+
+  }
+
+})
diff --git a/tobyfigs/packrat/packrat.lock b/tobyfigs/packrat/packrat.lock
new file mode 100644
index 0000000..0455458
--- /dev/null
+++ b/tobyfigs/packrat/packrat.lock
@@ -0,0 +1,422 @@
+PackratFormat: 1.4
+PackratVersion: 0.4.8.1
+RVersion: 3.4.2
+Repos: CRAN=http://cran.us.r-project.org
+
+Package: BH
+Source: CRAN
+Version: 1.62.0-1
+Hash: 14dfb3e8ffe20996118306ff4de1fab2
+
+Package: Formula
+Source: CRAN
+Version: 1.2-2
+Hash: 31c4658e6ae879eee45fcf7b6d109464
+
+Package: Hmisc
+Source: CRAN
+Version: 4.0-3
+Hash: 7e93c14752acdc238561bbdf56ad9ac2
+Requires: Formula, acepack, base64enc, data.table, ggplot2, gridExtra,
+    gtable, htmlTable, htmltools, latticeExtra, viridis
+
+Package: R6
+Source: CRAN
+Version: 2.2.1
+Hash: 530f0b839551f96ec991ce4f93156ee1
+
+Package: RColorBrewer
+Source: CRAN
+Version: 1.1-2
+Hash: c0d56cd15034f395874c870141870c25
+
+Package: Rcpp
+Source: CRAN
+Version: 0.12.13
+Hash: b9d8d89ff10ca2d8bfdaccdf83045ba1
+
+Package: acepack
+Source: CRAN
+Version: 1.4.1
+Hash: 5ed3b04bc1fa1a934188fc6022d3ba7c
+
+Package: assertthat
+Source: CRAN
+Version: 0.2.0
+Hash: e8805df54c65ac96d50235c44a82615c
+
+Package: backports
+Source: CRAN
+Version: 1.1.0
+Hash: 3f754426ca2b7dbbf63951fe5075373b
+
+Package: base64enc
+Source: CRAN
+Version: 0.1-3
+Hash: c590d29e555926af053055e23ee79efb
+
+Package: bindr
+Source: CRAN
+Version: 0.1
+Hash: e3a02070cf705d3ad1c5af1635a515a3
+
+Package: bindrcpp
+Source: CRAN
+Version: 0.2
+Hash: e011aef1e6b2f8ff15d833f56abcf0ba
+Requires: Rcpp, bindr, plogr
+
+Package: broom
+Source: CRAN
+Version: 0.4.2
+Hash: 7ebcffa46afb467e3f3c5687946f6e1a
+Requires: dplyr, plyr, psych, reshape2, stringr, tidyr
+
+Package: cellranger
+Source: CRAN
+Version: 1.1.0
+Hash: 4e1ef4d099b0c5fd531a3938cf4624bd
+Requires: rematch, tibble
+
+Package: checkmate
+Source: CRAN
+Version: 1.8.4
+Hash: fa853c0684e55ece9fc2de1bdcde9230
+Requires: backports
+
+Package: colorspace
+Source: CRAN
+Version: 1.3-2
+Hash: 0bf8618b585fa98eb23414cd3ab95118
+
+Package: curl
+Source: CRAN
+Version: 2.6
+Hash: 8162b82ca4809c0d63c30aedbd7348e0
+
+Package: data.table
+Source: CRAN
+Version: 1.10.4
+Hash: 2c7c47ccb40d01c7b214c8eff818108d
+
+Package: dichromat
+Source: CRAN
+Version: 2.0-0
+Hash: 08eed0c80510af29bb15f840ccfe37ce
+
+Package: digest
+Source: CRAN
+Version: 0.6.12
+Hash: e53fb8c58673df868183697e39a6a4d6
+
+Package: dplyr
+Source: CRAN
+Version: 0.7.4
+Hash: c028f5b3c279138df4ab90743ad68dc0
+Requires: BH, R6, Rcpp, assertthat, bindrcpp, glue, magrittr,
+    pkgconfig, plogr, rlang, tibble
+
+Package: evaluate
+Source: CRAN
+Version: 0.10.1
+Hash: 54d95f4ec6d0300100413ed0127d89ae
+Requires: stringr
+
+Package: fields
+Source: CRAN
+Version: 8.15
+Hash: 0236a31d894043dc56f152847ef9884b
+Requires: maps, spam
+
+Package: forcats
+Source: CRAN
+Version: 0.2.0
+Hash: e5a3b0b96a39f5581467b0c6366f7408
+Requires: magrittr, tibble
+
+Package: ggplot2
+Source: CRAN
+Version: 2.2.1
+Hash: 46e5cb78836848aa44655e577433f54b
+Requires: digest, gtable, lazyeval, plyr, reshape2, scales, tibble
+
+Package: ggthemes
+Source: CRAN
+Version: 3.4.0
+Hash: e9f85417507c5352fc05de654bb6f1ab
+Requires: assertthat, colorspace, ggplot2, scales
+
+Package: glue
+Source: CRAN
+Version: 1.1.1
+Hash: dfd5a27768175ae51d08dc6beba1ef11
+
+Package: gridExtra
+Source: CRAN
+Version: 2.3
+Hash: fa977bc1aab5588a08123b10ceb1ad3d
+Requires: gtable
+
+Package: gtable
+Source: CRAN
+Version: 0.2.0
+Hash: cd78381a9d3fea966ac39bd0daaf5554
+
+Package: haven
+Source: CRAN
+Version: 1.0.0
+Hash: 73ad91459a9d754299b3cbd2157a5f16
+Requires: BH, Rcpp, hms, readr, tibble
+
+Package: highr
+Source: CRAN
+Version: 0.6
+Hash: aa3d5b7912b5fed4b546ed5cd2a1760b
+
+Package: hms
+Source: CRAN
+Version: 0.3
+Hash: 3fca8a1c97e6cfb297fe3f4690f82c58
+
+Package: htmlTable
+Source: CRAN
+Version: 1.9
+Hash: b2a584fec4ff8c6206b01a55bed18e55
+Requires: checkmate, htmlwidgets, knitr, magrittr, stringr
+
+Package: htmltools
+Source: CRAN
+Version: 0.3.6
+Hash: 34e5c725d3ae30a76516792edcc16f24
+Requires: Rcpp, digest
+
+Package: htmlwidgets
+Source: CRAN
+Version: 0.9
+Hash: 7514f6ea9f3bef6d9b6c945095275302
+Requires: htmltools, jsonlite, yaml
+
+Package: httr
+Source: CRAN
+Version: 1.2.1
+Hash: 7de1f8f760441881804af7c1ff324340
+Requires: R6, curl, jsonlite, mime, openssl
+
+Package: jsonlite
+Source: CRAN
+Version: 1.5
+Hash: 9c51936d8dd00b2f1d4fe9d10499694c
+
+Package: knitr
+Source: CRAN
+Version: 1.16
+Hash: 3b8dc00d51027c6d041d56bc92136452
+Requires: digest, evaluate, highr, markdown, stringr, yaml
+
+Package: labeling
+Source: CRAN
+Version: 0.3
+Hash: ecf589b42cd284b03a4beb9665482d3e
+
+Package: latticeExtra
+Source: CRAN
+Version: 0.6-28
+Hash: efdf3a10a129129f7b8f4ae5a0d7ec3d
+Requires: RColorBrewer
+
+Package: lazyeval
+Source: CRAN
+Version: 0.2.0
+Hash: 3d6e7608e65bbf5cb170dab1e3c9ed8b
+
+Package: lubridate
+Source: CRAN
+Version: 1.6.0
+Hash: b90f4cbefe0b3c545dd68b22c66a8a12
+Requires: stringr
+
+Package: magrittr
+Source: CRAN
+Version: 1.5
+Hash: bdc4d48c3135e8f3b399536ddf160df4
+
+Package: maps
+Source: CRAN
+Version: 3.1.1
+Hash: 4a0a1b1764bfe2c41c6fc21790408caf
+
+Package: markdown
+Source: CRAN
+Version: 0.8
+Hash: 045d7c594d503b41f1c28946d076c8aa
+Requires: mime
+
+Package: mime
+Source: CRAN
+Version: 0.5
+Hash: 463550cf44fb6f0a2359368f42eebe62
+
+Package: mnormt
+Source: CRAN
+Version: 1.5-5
+Hash: d0d5efbb1fb26d2dc5f9394c223084b5
+
+Package: modelr
+Source: CRAN
+Version: 0.1.0
+Hash: 7c9848bf4d734f38b8ce91022d8de949
+Requires: broom, dplyr, lazyeval, magrittr, purrr, tibble, tidyr
+
+Package: munsell
+Source: CRAN
+Version: 0.4.3
+Hash: f96d896947fcaf9b6d0074002e9f4f9d
+Requires: colorspace
+
+Package: openssl
+Source: CRAN
+Version: 0.9.6
+Hash: 5f4711e142a44655dfea4d64fcf2f641
+
+Package: packrat
+Source: CRAN
+Version: 0.4.8-1
+Hash: 6ad605ba7b4b476d84be6632393f5765
+
+Package: pkgconfig
+Source: CRAN
+Version: 2.0.1
+Hash: 0dda4a2654a22b36a715c2b0b6fbacac
+
+Package: plogr
+Source: CRAN
+Version: 0.1-1
+Hash: fb19215402e2d9f1c7f803dcaa806fc2
+
+Package: plyr
+Source: CRAN
+Version: 1.8.4
+Hash: 9598fbff0cce13bee89ae3a50230d93f
+Requires: Rcpp
+
+Package: psych
+Source: CRAN
+Version: 1.7.5
+Hash: 0c076a96de916d0d26d866e83909d961
+Requires: mnormt
+
+Package: purrr
+Source: CRAN
+Version: 0.2.3
+Hash: d741f262ed1b74543e33c5426a497996
+Requires: magrittr, rlang, tibble
+
+Package: readr
+Source: CRAN
+Version: 1.1.1
+Hash: 9f1ac7e796ff384e3ce9e082fdb2d17c
+Requires: BH, R6, Rcpp, hms, tibble
+
+Package: readxl
+Source: CRAN
+Version: 1.0.0
+Hash: 006771e12b2ae14b7dd6f77dfb70d802
+Requires: Rcpp, cellranger, tibble
+
+Package: rematch
+Source: CRAN
+Version: 1.0.1
+Hash: ad4faf59e7611117ff165817074c50c7
+
+Package: reshape2
+Source: CRAN
+Version: 1.4.2
+Hash: cca35d2c477d99d8bc53908939eb659e
+Requires: Rcpp, plyr, stringr
+
+Package: rlang
+Source: CRAN
+Version: 0.1.2
+Hash: c17b638e3282aa222a47112d7a60bc66
+
+Package: rvest
+Source: CRAN
+Version: 0.3.2
+Hash: c69f7526520bad66fd2111ebe8b1364b
+Requires: httr, magrittr, selectr, xml2
+
+Package: scales
+Source: CRAN
+Version: 0.5.0
+Hash: 1e70c4e1dcc1c51847edc9f4bc11a81e
+Requires: R6, RColorBrewer, Rcpp, dichromat, labeling, munsell, plyr,
+    viridisLite
+
+Package: selectr
+Source: CRAN
+Version: 0.3-1
+Hash: 367275e3dcdd208339e131c7a41bec56
+Requires: stringr
+
+Package: spaero
+Source: CRAN
+Version: 0.2.0
+Hash: 3e128e5aa0bed3af97cdbb3147e4887f
+
+Package: spam
+Source: CRAN
+Version: 1.4-0
+Hash: a82fb29f79c3d556fe2584db0000bfbc
+
+Package: stringi
+Source: CRAN
+Version: 1.1.5
+Hash: b6308e49357a0b475f433599e0d8b5eb
+
+Package: stringr
+Source: CRAN
+Version: 1.2.0
+Hash: 25a86d7f410513ebb7c0bc6a5e16bdc3
+Requires: magrittr, stringi
+
+Package: tibble
+Source: CRAN
+Version: 1.3.3
+Hash: e34ba28025602bd2341af5b8845156c9
+Requires: Rcpp, rlang
+
+Package: tidyr
+Source: CRAN
+Version: 0.6.3
+Hash: e99b95cbb31017cebaab8ba3ddabcb67
+Requires: Rcpp, dplyr, lazyeval, magrittr, stringi, tibble
+
+Package: tidyverse
+Source: CRAN
+Version: 1.1.1
+Hash: 72d5fada870c90b835bbdfc281283c99
+Requires: broom, dplyr, forcats, ggplot2, haven, hms, httr, jsonlite,
+    lubridate, magrittr, modelr, purrr, readr, readxl, rvest, stringr,
+    tibble, tidyr, xml2
+
+Package: viridis
+Source: CRAN
+Version: 0.4.0
+Hash: 5bdac1bcf74a10a7a96f82191f498ab7
+Requires: ggplot2, gridExtra, viridisLite
+
+Package: viridisLite
+Source: CRAN
+Version: 0.2.0
+Hash: 10f0c25af3dc84eaae10f5854f47efdb
+
+Package: xml2
+Source: CRAN
+Version: 1.1.1
+Hash: f312e4beff105290588bf035a69428b0
+Requires: BH, Rcpp
+
+Package: yaml
+Source: CRAN
+Version: 2.1.14
+Hash: c81230c3a7d9ba20607ad6b4331173d1
diff --git a/tobyfigs/packrat/packrat.opts b/tobyfigs/packrat/packrat.opts
new file mode 100644
index 0000000..7f6839c
--- /dev/null
+++ b/tobyfigs/packrat/packrat.opts
@@ -0,0 +1,15 @@
+auto.snapshot: TRUE
+use.cache: FALSE
+print.banner.on.startup: auto
+vcs.ignore.lib: TRUE
+vcs.ignore.src: FALSE
+external.packages: 
+local.repos: 
+load.external.packages.on.startup: TRUE
+ignored.packages: 
+quiet.package.installation: TRUE
+snapshot.recommended.packages: FALSE
+snapshot.fields:
+    Imports
+    Depends
+    LinkingTo
diff --git a/tobyfigs/res.RData b/tobyfigs/res.RData
new file mode 100644
index 0000000..b1eec87
Binary files /dev/null and b/tobyfigs/res.RData differ
diff --git a/tobyfigs/snapshots-auc-calculation.R b/tobyfigs/snapshots-auc-calculation.R
new file mode 100644
index 0000000..2932d40
--- /dev/null
+++ b/tobyfigs/snapshots-auc-calculation.R
@@ -0,0 +1,123 @@
+
+library(tidyverse)
+library(reshape2)
+
+##Used to calculate AUC directly from snapshots data.
+
+
+#load simulated cases
+load("cases.RData")
+
+#Calculate EWS for times series data using spaero
+get_ews <- function(emerge){
+
+  ## Wrapper for spaero::get_stats, sets function variables and returns only stats as a data frame.
+  #spaero parameters are hardcoded in here
+  ews <- function(x){
+    return(
+      as.data.frame(spaero::get_stats(x, center_trend ="local_constant",
+                                      center_kernel = "uniform",
+                                      center_bandwidth =35,
+                                      stat_trend = "local_constant",
+                                      stat_kernel = "uniform",
+                                      stat_bandwidth=35,
+                                      lag = 1)$stats)
+    )
+  }
+
+  emerge %>%
+    group_by(sim) %>%
+    do(ews = ews(.$I)) -> a
+  bind_rows(a$ews) -> b
+  emerge_ews <- cbind(emerge,b)
+}
+
+#calculate taus
+get_taus <- function(df){
+  df <- melt(df, id.vars = c("time", "S", "I", "R", "N", "cases", "beta_t", "sim"))
+  df %>%
+    group_by(variable, sim) %>%
+    summarise(tau = stats::cor(value,time, use="pairwise.complete.obs", method="kendall")) -> taus
+  return(taus)
+}
+
+#calculate AUC
+get_auc <- function(df){
+  calc_auc <- function(predictions, is_null){
+    test <- !is.na(predictions)
+    predictions <- predictions[test]
+    is_null <- is_null[test]
+    r <- rank(predictions)
+    r1 <- sum(r[!is_null])
+    n1 <- sum(!is_null)
+    n2 <- sum(is_null)
+    (r1 - n1 * (n1 + 1) / 2) / (n1 * n2)
+  }
+  df %>%
+    group_by(variable) %>%
+    summarise(AUC = calc_auc(tau, isnull)) -> out
+
+  new_levels <-  c("Mean", "Variance", "Var. 1st Diff",
+                   "Index of Dis.", "Autocovar.", "Autocorr.",
+                   "Decay time", "Coeff. Var", "Skewness",
+                   "Kurtosis")
+  out$variable <- factor(out$variable,levels(out$variable)[c(6,1,2,7,3,4,5,8,9,10)])
+  #names(new_levels) <- levels(out$variable)
+  levels(out$variable) <- new_levels
+
+
+
+  return(out)
+}
+
+
+
+
+#Simulates 1000 replicates, N=1e6, 20yrs of data, rng seed = 1
+
+calculate_auc <- function(infectious_period){
+
+  if (infectious_period == "-weekly"){
+    i <- which(dplyr::near(process_des_mat$infectious_days, 7))
+  } else {
+    i <- which(dplyr::near(process_des_mat$infectious_days, 30))
+  }
+  emerge_data <- simulated_procs[[i]]$test
+  null_data <- simulated_procs[[i]]$null
+
+  calculate_auc_agg <- function(aggregation_period){
+    edf <- emerge_data
+    ndf <- null_data
+    is_snap <- function(x, period){
+      mod <- x %% period
+      dplyr::near(mod, 0)
+    }
+    if(aggregation_period == "-monthly"){
+      per <- 30
+    } else {
+      per <- 7
+    }
+    edf <- edf[is_snap(edf$time, per), ]
+    ndf <- ndf[is_snap(ndf$time, per), ]
+
+    emerge_ews <- get_ews(edf)
+    not_ews <- get_ews(ndf)
+    taus_test <- get_taus(emerge_ews)
+    taus_null <- get_taus(not_ews)
+    taus_test$isnull <- FALSE
+    taus_null$isnull <- TRUE
+    taus <- rbind(taus_test, taus_null)
+
+    df <- get_auc(taus)
+    df$`Infectious period` <- as.factor(infectious_period)
+    df$`Aggregation period` <- as.factor(aggregation_period)
+    return(df)
+  }
+
+  df <- do.call(rbind,lapply(c("-weekly","-monthly"),calculate_auc_agg))
+  return(df)
+}
+
+auc_data <- do.call(rbind,lapply(c("-weekly", "-monthly"),calculate_auc))
+
+write_csv(auc_data,"snapshots-auc.csv")
diff --git a/tobyfigs/snapshots-auc-plot.R b/tobyfigs/snapshots-auc-plot.R
new file mode 100644
index 0000000..478f3c6
--- /dev/null
+++ b/tobyfigs/snapshots-auc-plot.R
@@ -0,0 +1,84 @@
+library(ggplot2)
+library(ggthemes) ##
+library(tidyverse)
+#library(plyr)
+library(reshape2)
+library(colorspace)
+
+auc_data <- read.csv("snapshots-auc.csv")
+auc_data$`Infectious period` <- as.factor(auc_data$Infectious.period)
+levels(auc_data$`Infectious period`) <- c(expression(paste(gamma, "=1/30")),expression(paste(gamma, "=1/7")))
+
+auc_data$`Aggregation period` <- as.factor(auc_data$Aggregation.period)
+levels(auc_data$`Aggregation period`) <- c("Monthly~snapshots","Weekly~snapshots")
+
+auc_data$variable <- factor(auc_data$variable, levels = levels(auc_data$variable)[c(7,10,9,5,2,1,4,3,8,6)])
+
+
+auc_plot <- function(df){
+  
+  ## Eric's AUC gradient ##
+  # Given error of ~0.5, set number of levels to 10.  For a higher precision, use nlevels = 20 
+  nlevels = 100  # number of bins = number of colors
+  # Color scale with {colorspace}
+  AUC_colors <- diverge_hcl(
+    n=nlevels,  # number of colors
+    h = c(45, 225),  # Hues (low, hi)
+    c = 100, # fixed Chroma or Chroma Range (edges, center)
+    l = c(90, 10),  # Lightness range (edges, center) 
+    power = 1  # exponent
+  )
+  names(AUC_colors) <- seq(0,1,1/(nlevels-1))
+  
+  
+  ggplot(df) + 
+    # geom_pointrange(aes(x=variable,y=AUC, ymin=AUC-AUC_err, ymax= AUC+AUC_err, color=variable)) +
+    geom_bar(aes(x=variable,y=AUC-0.5, fill = AUC, color = AUC),stat="identity"  ) + #, fill = "#088E7C") +
+    facet_grid(`Infectious period`~`Aggregation period`,labeller = label_parsed)+
+    geom_rangeframe(colour ="black") +
+    scale_fill_gradientn(limits = c(0,1),colours=AUC_colors) +
+    scale_color_gradientn(limits = c(0,1),colours=AUC_colors) +
+    scale_y_continuous(name = "AUC",limits = c(-0.5,0.5),labels=c("0.0","0.25","0.5","0.75","1.0"))
+  #scale_x_continuous(name = "")+#,labels=c("Variance","Variance convexity","Autocovariance",
+  #                                      "Autocorrelation","Decay time", "Mean",
+  #                                      "Index of Dispersion", "Coefficient of variation",
+  #                                      "Skewness", "Kurtosis"))+
+  # theme_minimal() +
+  #theme(axis.text.x = element_text(angle = 45, hjust = 1))
+}
+
+par(font.lab = 3)
+text_color <- "black"
+background_color <- "white"
+font_chosen <- "Helvetica"
+auc_plot(auc_data) +
+  theme(text = element_text(color=text_color, size = 18, family=font_chosen),
+        title = element_text(size = 14),
+        line = element_line(color=text_color),
+        rect = element_rect(color="black"),
+        axis.title.x=element_blank(),
+        axis.title.y = element_text(color=text_color,face = "plain", 
+                                    family=font_chosen, size = 14),
+        axis.text.x = element_text(color = text_color, angle = 45,
+                                   vjust = 1, hjust=1),#family=font_chosen, size = 14),
+        axis.text.y = element_text(color = text_color, family=font_chosen),
+        axis.ticks = element_line(color=rgb(0,0,0,.25),line),
+        axis.line = element_blank(),
+        legend.text = element_text(color = text_color, family=font_chosen),
+        legend.title = element_text(color = text_color, family = font_chosen),
+        panel.grid.major = element_blank(),
+        panel.grid.major.y = element_line(color = rgb(0,0,0,.25)),
+        panel.grid.minor = element_blank(),
+        panel.border = element_blank(),
+        panel.background = element_rect(color=NA,fill=NA),
+        panel.spacing = unit(2, "lines"),
+        strip.text.y= element_text(color = text_color,angle = 270, family=font_chosen),
+        strip.text.x= element_text(color = text_color, family=font_chosen),
+        strip.background = element_rect(fill=NA),
+        plot.background = element_rect(color=NA, fill=background_color),
+        plot.margin = unit(c(1,1,1,1), "cm"))
+
+
+ggsave("snapshots-auc.pdf",width =1.2*8.64,height=1.2*4.20)
+
+
diff --git a/tobyfigs/snapshots-auc.csv b/tobyfigs/snapshots-auc.csv
new file mode 100644
index 0000000..5371124
--- /dev/null
+++ b/tobyfigs/snapshots-auc.csv
@@ -0,0 +1,41 @@
+variable,AUC,Infectious period,Aggregation period
+Variance,1,-weekly,-weekly
+Var. 1st Diff,0.97501,-weekly,-weekly
+Autocovar.,1,-weekly,-weekly
+Autocorr.,0.999677,-weekly,-weekly
+Decay time,0.999677,-weekly,-weekly
+Mean,1,-weekly,-weekly
+Index of Dis.,0.999993,-weekly,-weekly
+Coeff. Var,0.179537,-weekly,-weekly
+Skewness,0.489022,-weekly,-weekly
+Kurtosis,0.444602,-weekly,-weekly
+Variance,0.999798,-weekly,-monthly
+Var. 1st Diff,0.9969355,-weekly,-monthly
+Autocovar.,0.9749975,-weekly,-monthly
+Autocorr.,0.9436765,-weekly,-monthly
+Decay time,0.946374245472837,-weekly,-monthly
+Mean,0.999709,-weekly,-monthly
+Index of Dis.,0.9959355,-weekly,-monthly
+Coeff. Var,0.513333,-weekly,-monthly
+Skewness,0.613984,-weekly,-monthly
+Kurtosis,0.603714,-weekly,-monthly
+Variance,0.999835,-monthly,-weekly
+Var. 1st Diff,0.912403,-monthly,-weekly
+Autocovar.,0.999702,-monthly,-weekly
+Autocorr.,0.9534735,-monthly,-weekly
+Decay time,0.9534735,-monthly,-weekly
+Mean,0.999998,-monthly,-weekly
+Index of Dis.,0.97325,-monthly,-weekly
+Coeff. Var,0.052677,-monthly,-weekly
+Skewness,0.4368965,-monthly,-weekly
+Kurtosis,0.344703,-monthly,-weekly
+Variance,0.999623,-monthly,-monthly
+Var. 1st Diff,0.987731,-monthly,-monthly
+Autocovar.,0.99791,-monthly,-monthly
+Autocorr.,0.9629065,-monthly,-monthly
+Decay time,0.9628425,-monthly,-monthly
+Mean,0.99998,-monthly,-monthly
+Index of Dis.,0.9887905,-monthly,-monthly
+Coeff. Var,0.469341,-monthly,-monthly
+Skewness,0.5989215,-monthly,-monthly
+Kurtosis,0.5566645,-monthly,-monthly
diff --git a/tobyfigs/snapshots-auc.pdf b/tobyfigs/snapshots-auc.pdf
new file mode 100644
index 0000000..341d1f9
Binary files /dev/null and b/tobyfigs/snapshots-auc.pdf differ