#
#  R code to demonstrate application of GP emulators.  Train emulator on
#  tropospheric ozone burden from the 250 training runs and use it to
#  predict the burden for the validation runs, comparing with CTM results.
#
#  Author: Oliver Wild, based on original code by Ed Ryan
#-----------------------------------------------------------------------

# Remove all current objects stored in R environment
rm(list = ls())
cat("\014")
graphics.off()

# Load required R packages (install on first use if necessary) 
#install.packages("lhs")
#install.packages("DiceKriging")
#install.packages("DiceOptim")
library(lhs)
library(DiceKriging)
library(DiceOptim)

# Read in normalised parameter values (removing 2 rows of headers)
params_ALL <-read.csv('parameter_values_normalised.csv',skip=2,header=FALSE)

# Select all 36 parameters from training runs (removing row label)
training <- as.matrix(params_ALL[1:250,2:37], rownames.force = NA)
X <- training

# Select parameters from validation runs
valid <- as.matrix(params_ALL[251:280,2:37], rownames.force = NA)

# Read in outputs (O3 burden only as an example, removing row labels)
outputs_ALL <-read.csv('model_budgets_all_runs.csv',skip=3,header=FALSE)
outputs <- as.matrix(outputs_ALL[1:250,2:2], rownames.force = NA)
out_val <- as.matrix(outputs_ALL[251:280,2:2], rownames.force = NA)

# Create empty matrices to store the emulator outputs
N=30    # Number of validation runs
V=1     # Number of variables (O3 burden only)
outputs_emul <- matrix(NA,nrow=N,ncol=V)
outputs_min  <- matrix(NA,nrow=N,ncol=V)
outputs_max  <- matrix(NA,nrow=N,ncol=V)
outputs_col  <- matrix(NA,nrow=N,ncol=V*4)

for (j in 1:V){
    print(j)
    y <- as.matrix(outputs[1:250,j], rownames.force = NA)
    m <- km(~ ., design = X, response = y, covtype = "matern5_2")
    outputs_emul[,j] <- matrix(t(predict.km(m, valid, "UK", se.compute = TRUE, checkNames = FALSE)$mean))
    outputs_min[,j] <- matrix(t(predict.km(m, valid, "UK", se.compute = TRUE, checkNames = FALSE)$lower95))
    outputs_max[,j] <- matrix(t(predict.km(m, valid, "UK", se.compute = TRUE, checkNames = FALSE)$upper95))
}

# Write emulator outputs to csv file with appropriate headers (model, emulator mean, lower95th, upper95th)
outputs_col[] <- cbind(out_val[1:N,1],outputs_emul[1:N,1],outputs_min[1:N,1],outputs_max[1:N,1])
colnames(outputs_col) <- c("Model_run","Emulated","Emulated_low95th","Emulated_high95th")
write.table(data.frame(outputs_col), file = paste("Outputs.csv",sep=""), sep=",",row.names=FALSE,col.names=TRUE)