## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup_lies, warning=FALSE------------------------------------------------
library(EcoEnsemble)
library(mgcv)
library(ggplot2)
library(cowplot)
set.seed(1234)

## ----configure_priors---------------------------------------------------------
d <- 3 # The number of variables.

priors <- EnsemblePrior(d,
                        ind_st_params = IndSTPrior(AR_params=c(10,2)), 
                        ind_lt_params = IndLTPrior("beta",
                                                   cor_params = list(matrix(5, d, d), 
                                                                     matrix(1, d, d))),
                        sha_st_params = ShaSTPrior(AR_params=c(10,2)),
)

## ----sampling priors, eval = FALSE--------------------------------------------
#  M <- 4 # The number of models we are considering.
#  prior_density <- prior_ensemble_model(priors, M = M)
#  ex.fit <- rstan::extract(prior_density$samples) # Extracted samples

## ----truth_best_guesses, eval = FALSE-----------------------------------------
#  Time <- 50
#  
#  true_par <-which.max(ex.fit$lp__)
#  
#  latent <- matrix(NA, Time, (M+2)*d)
#  
#  #Priors on initial values
#  latent[1, 1:d] <- rnorm(d, 0, 1)
#  latent[1, -(1:d)] <- t(chol(ex.fit$SIGMA_init[true_par,-(1:d) ,-(1:d)]))
#                                            %*% rnorm((M+1)*d, 0, 1)
#  
#  #Find all the latent values of the dynamic linear model
#  SIGMA <- ex.fit$SIGMA[true_par,,]
#  SIGMA_chol <- t(chol(SIGMA))
#  for (t in 2:Time) {
#    latent[t,] <- ex.fit$AR_params[true_par,] * latent[t-1,] + SIGMA_chol
#                                            %*% rnorm((M+2)*d, 0, 1)
#  }
#  
#  #The truth is simply the first d values
#  truth_latent <- latent[,1:d]
#  
#  #The best guesses are sums of the truth and discrepancies
#  simulator_best_guesses <- array(NA,dim=c(Time,d,M))
#  for(i in 1:M){
#    simulator_best_guesses[,,i] <- t(
#      t(latent[,1:d] + latent[,(d+1):(2*d)] + latent[,(1+i) * d + (1:d)]) +
#        ex.fit$ind_lt[true_par,i,] + ex.fit$sha_lt[true_par,] )
#  }
#  

## ----load_datasets, echo = FALSE----------------------------------------------
#Here I will load in the data
M <- 4
Time <- 50
load("data/SyntheticData_PriorSamples.RData")

## ----plot_data, fig.dim = c(7, 4)---------------------------------------------
plot(truth_latent[,1], ylim=c(-2.5, 3), pch = 19)
lines(simulator_best_guesses[,1,1], col = 2)
lines(simulator_best_guesses[,1,2], col = 3)
lines(simulator_best_guesses[,1,3], col = 4)
lines(simulator_best_guesses[,1,4], col = 5)


## ----add_noise----------------------------------------------------------------
Times_obs <- round(Time * 0.8)
obs <- matrix(NA,Times_obs,d)
for(i in 1:d){
  g1 <- gam(y~s(x,k = 15),data=data.frame(x=1:Times_obs,y = truth_latent[1:Times_obs,i]))
  obs[,i] <- predict(g1)
}

obs.cov <- cov(obs - truth_latent[1:Times_obs,])


model.cov <- array(0,dim=c(M,d,d))
models_output <- array(NA,dim=c(M,Time,d))
for (j in 1:M){
  for(i in 1:d){
    g1 <- gam(y~s(x, k = 15),data=data.frame(x=1:Time,y=simulator_best_guesses[,i,j]))
    models_output[j,,i] <- predict(g1)
  }
  model.cov[j,,] <- cov(models_output[j,,] - simulator_best_guesses[,,j])
}


## ----plot_obs, fig.dim = c(7, 6)----------------------------------------------
#Observations and model outputs
plot(obs[,1], ylim=c(-2.5, 3), xlim = c(1,Time), pch = 19)
lines(models_output[1,,1], col=2, lwd = 2)
lines(models_output[2,,1], col=3, lwd = 2)
lines(models_output[3,,1], col=4, lwd = 2)
lines(models_output[4,,1], col=5, lwd = 2)

## -----------------------------------------------------------------------------
#Create the data frames that we'll use for EcoEnsemble
val_obs <- data.frame(obs); cov_obs <- obs.cov
val_model_1 <- data.frame(models_output[1,,]); cov_model_1 <- model.cov[1,,]
val_model_2 <- data.frame(models_output[2,,]); cov_model_2 <- model.cov[2,,]
val_model_3 <- data.frame(models_output[3,,]); cov_model_3 <- model.cov[3,,]
val_model_4 <- data.frame(models_output[4,,]); cov_model_4 <- model.cov[4,,]

#Set the dimnames to ensure EcoEnsemble can identify the information.
SPECIES_NAMES <- c("Species 1", "Species 2", "Species 3")
dimnames(val_obs) <- list(paste(1:Times_obs), SPECIES_NAMES)
dimnames(val_model_1) <- list(paste(1:Time), SPECIES_NAMES)
dimnames(val_model_2) <- list(paste(1:Time), SPECIES_NAMES)
dimnames(val_model_3) <- list(paste(1:Time), SPECIES_NAMES)
dimnames(val_model_4) <- list(paste(1:Time), SPECIES_NAMES)

dimnames(cov_obs) <- list(SPECIES_NAMES, SPECIES_NAMES)
dimnames(cov_model_1) <- list(SPECIES_NAMES, SPECIES_NAMES)
dimnames(cov_model_2) <- list(SPECIES_NAMES, SPECIES_NAMES)
dimnames(cov_model_3) <- list(SPECIES_NAMES, SPECIES_NAMES)
dimnames(cov_model_4) <- list(SPECIES_NAMES, SPECIES_NAMES)


## ----fitting, eval = FALSE----------------------------------------------------
#  fit <- fit_ensemble_model(observations = list(val_obs, cov_obs),
#                            simulators = list(list(val_model_1, cov_model_1, "Model 1"),
#                                              list(val_model_2, cov_model_2, "Model 2"),
#                                              list(val_model_3, cov_model_3, "Model 3"),
#                                              list(val_model_4, cov_model_4, "Model 4")),
#                            priors = EnsemblePrior(d),
#                            control = list(adapt_delta = 0.9))
#  samples <- generate_sample(fit)

## ----results, eval = FALSE----------------------------------------------------
#  var_index <- 1
#  
#  df <- data.frame("Year" = paste(1:Time),
#                   "Ensemble" = apply(samples@mle[, var_index, ], 1, median),
#                   "Lower"  = apply(samples@samples[, var_index, ], 1, quantile, 0.1),
#                   "Upper"  = apply(samples@samples[, var_index, ], 1, quantile, 0.9),
#                   "Actual" = truth_latent[,var_index])
#  df$Year <- as.numeric(df$Year)
#  df <-  reshape2::melt(df, id.vars=c("Year", "Lower", "Upper"), variable.name="Simulator")
#  
#  # We only want uncertainty bands for the Ensemble values, else we set zero width bands
#  df[df$Simulator == "Actual", c("Lower", "Upper")] <- df[df$Simulator != "Actual", "value"]
#  
#  ggplot(df, aes(x=`Year`, y=`value`)) +
#    geom_line(aes(group=`Simulator`,colour=`Simulator`)) +
#    geom_ribbon(aes(ymin=`Lower`, ymax =`Upper`, fill = `Simulator`), alpha=0.2)
#  

## ----plot_truth, echo = FALSE, fig.dim = c(7, 4)------------------------------
knitr::include_graphics("data/p_truth.png")

## ----eval = FALSE-------------------------------------------------------------
#  fit_M1 <- fit_ensemble_model(observations = list(val_obs, cov_obs),
#                               simulators = list(list(val_model_1, cov_model_1, "Model 1")),
#                               priors = EnsemblePrior(d),
#                               control = list(adapt_delta = 0.9))
#  samples_M1 <- generate_sample(fit_M1)
#  
#  fit_M2 <- fit_ensemble_model(observations = list(val_obs, cov_obs),
#                               simulators = list(list(val_model_1, cov_model_1, "Model 1"),
#                                              list(val_model_2, cov_model_2, "Model 2")),
#                            priors = EnsemblePrior(d),
#                            control = list(adapt_delta = 0.9))
#  samples_M2 <- generate_sample(fit_M2)
#  
#  fit_M3 <- fit_ensemble_model(observations = list(val_obs, cov_obs),
#                            simulators = list(list(val_model_1, cov_model_1, "Model 1"),
#                                              list(val_model_2, cov_model_2, "Model 2"),
#                                              list(val_model_3, cov_model_3, "Model 3")),
#                            priors = EnsemblePrior(d),
#                            control = list(adapt_delta = 0.9))
#  samples_M3 <- generate_sample(fit_M3)
#  
#  plot_grid(
#    plot(samples_M1, variable=1) + ggtitle("1 Model")  + theme(legend.position = "none"),
#    plot(samples_M2, variable=1) + ggtitle("2 Models") + theme(legend.position = "none"),
#    plot(samples_M3, variable=1) + ggtitle("3 Models") + theme(legend.position = "none"),
#    plot(samples   , variable=1) + ggtitle("4 Models") + theme(legend.position = "none"))

## ----model_comparison, echo = FALSE, fig.dim = c(7, 6)------------------------
knitr::include_graphics("data/p_NumberOfModels.png")

## ----quantify_variance, eval = FALSE------------------------------------------
#  def.par <- par(no.readonly=TRUE) #old pars
#  par(mfrow = c(d, 1))
#  legend_ys <- c(0.4, 0.18, 0.12)
#  for(i in 1:d){
#    plot.ts(apply(samples_M1@samples[40:50,i,],1, var),
#            main = paste0("Species ", i), ylab="Variance" )
#    lines(apply(samples_M2@samples[40:50,i,],1, var), col = 2)
#    lines(apply(samples_M3@samples[40:50,i,],1, var), col = 3)
#    lines(apply(samples@samples[40:50,i,],1, var), col = 4)
#    legend(1, legend_ys[i], legend = c("1 Model", "2 Models", "3 Models", "4 Models"),
#            col = 1:4, lty = 1)
#  }
#  par(def.par)

## ----plot_variances, echo=FALSE, fig.dim = c(7,6)-----------------------------
load("data/SyntheticData_AddingModelsData.RData")
def.par <- par(no.readonly=TRUE) #old pars
par(mfrow = c(d, 1), mar = c(4, 4, 1.4, 1.4))
legend_ys <- c(0.4, 0.18, 0.12)
for(i in 1:d){
  plot.ts(var_m1[, i], 
          main = paste0("Species ", i), ylab="Variance" )
  lines(var_m2[, i], col = 2)
  lines(var_m3[, i], col = 3)
  lines(var_m4[, i], col = 4)
  legend(1, legend_ys[i], legend = c("1 Model", "2 Models", "3 Models", "4 Models"),
          col = 1:4, lty = 1)
}
par(def.par)