devtools::install_github("seroanalytics/serosolver")
library(serosolver)

library(serosolver)
library(plyr)
library(data.table)
library(ggplot2)

## Load in example parameter values and antigenic map
data(example_par_tab)
data(example_antigenic_map)

## Get all possible infection times
strain_isolation_times <- unique(example_antigenic_map$inf_times)

## Vector of strains that have titres (note only one representative strain per time)
sampled_viruses <- seq(min(strain_isolation_times), max(strain_isolation_times), by=2)

## Times at which serum samples can be taken
sampling_times <- 2010:2015

## Number of serum samples taken
n_samps <- 2

## Simulate some random attack rates
attack_rates <- runif(length(strain_isolation_times), 0.05, 0.15)

## Simulate a full serosurvey with these parameters
all_simulated_data <- simulate_data(par_tab=example_par_tab, group=1, n_indiv=50,
                                  strain_isolation_times=strain_isolation_times,
                                  measured_strains=sampled_viruses,
                                  sampling_times=2010:2015, nsamps=n_samps,
                                  antigenic_map=example_antigenic_map,
                                  age_min=10,age_max=75,
                                  attack_rates=attack_rates, repeats=2)

## Pull out the simulated titre data and infection histories
titre_dat <- all_simulated_data$data
ages <- all_simulated_data$ages
example_inf_hist <- all_simulated_data$infection_histories
example_titre_dat <- merge(titre_dat, ages)

## Run the MCMC
# This example uses prior version 2 (i.e. beta prior on phi with parameters alpha, beta)
# We have to remove the explicit specification of phi in the parameter table
par_tab <- example_par_tab[example_par_tab$names != "phi",]
res <- run_MCMC(par_tab, example_titre_dat, example_antigenic_map, 
                filename="test", version=2,
                mcmc_pars=c(adaptive_period=20000, iterations=80000, 
                            inf_propn=0.5,hist_sample_prob=0.5,
                            save_block=10000,thin=10,thin_hist=100))

## Read in the MCMC chains and plot posteriors
chain <- read.csv(res$chain_file)
inf_chain <- data.table::fread(res$history_file)
plot(coda::as.mcmc(chain[chain$sampno > 20000,c("mu","wane","lnlike")]))

# Plot model predicted titres for a subset of individuals
plot_infection_histories(chain = chain,infection_histories = inf_chain,
                         titre_dat = example_titre_dat,individuals=c(1:4),
                         antigenic_map=example_antigenic_map,par_tab=par_tab)