#
#  R code to demonstrate global sensitivity analysis.  Train emulator on
#  tropospheric ozone burden from the 250 training runs and use the fast99
#  routine to calculate the sensitivity indices for each parameter. The
#  indices are output as percentages for each contributing parameter.
#
#  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("sensitivity")
#install.packages("DiceKriging")
#install.packages("DiceOptim")
library(sensitivity)
library(DiceKriging)
library(DiceOptim)

# Inputs
N=10000

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

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

# 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)


SI <- matrix(-9999,nrow=36,ncol=1)
for (j in 1:1){
  y <- as.matrix(outputs[,j], rownames.force = NA)
  m <- km(~ ., design = X, response = y, covtype = "matern3_2")
  kriging.mean <- function(Xnew, m){
    predict.km(m, Xnew, "UK", se.compute = FALSE, checkNames = FALSE)$mean
  }
  temp <- fast99(model = kriging.mean, factors = 36, n = N, q = "qunif", q.arg = list(min = 0, max = 1), m = m)
  SI[,j] <- t(as.matrix(temp$D1/temp$V))
}

# Write emulator outputs to csv file as percentages
out = data.frame(SI=SI*100)
names(out)=c("O3 Burd")
write.csv(out, file = paste("SIs_O3burden.csv",sep=""), row.names=FALSE)