twd <- getwd()
source('BuildEmulator/BuildEmulator.R')
source('BuildEmulator/BuildEmulatorNSt.R')
source('HistoryMatching/HistoryMatching.R')
load("Demonstrations/globalCanAM4W1.RData")
head(tData)
cands <- names(tData)[2:14]
Noise <- rnorm(dim(tData)[1],0,0.5)
tData <- cbind(tData,Noise)
cands <- c(cands, "Noise")
myem.lm.canTOA <- EMULATE.lm(Response="BALX", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 20)
myem.gp.canTOA = EMULATE.gpstan(meanResponse=myem.lm.canTOA, nugget=0.0001, tData=tData, additionalVariables=NULL,CreatePredictObject = FALSE,FastVersion = TRUE)
tLOOs <- stanLOOplot(StanEmulator = myem.gp.canTOA, ParamNames=myem.gp.canTOA$Names)

# Standard Errors to check for nonstationarity

std.err <- CalStError(tLOOs, tData$BALX)
par(mfrow = c(3, 3), mar = c(4, 4, 2, 1))
for(i in 1:dim(tData[, myem.gp.canTOA$Names])[2]) {
  plot(tData[, myem.gp.canTOA$Names][, i], std.err, pch = 20, xlab = myem.gp.canTOA$Names[i], 
       ylab = 'std error', ylim = c(-2, 2))
}
MixtureData <- data.frame(cbind(tData$CSIGMA, std.err))
names(MixtureData) <- c('CSIGMA', 'std.err')
# bayesian mixture in stan
myem.mixture <- MIXTURE.design(y.std = MixtureData[, 2], X =matrix(MixtureData[, 1], ncol = 1), L=2)
LOOMixture <- data.frame(cbind(MixtureData, myem.mixture$MixtureMat))
names(LOOMixture) <- c(names(LOOMixture), 'red', 'blue')
ggplot(LOOMixture, aes(x = CSIGMA, y = std.err, col = rgb(red = red, blue = blue, green = 0))) +
  geom_point(size = 2) + scale_color_identity()
myem.gp.nst <- EMULATE.gpstanNSt(myem.lm.canTOA, nugget = 0.0001, tData = tData, additionalVariables=NULL, 
                                 FastVersion = TRUE, mixtureComp = myem.mixture)
LOOs.nst <- stanLOOplotNSt(myem.gp.nst, ParamNames = myem.gp.canTOA$Names)

par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
ValidPlot(tLOOs,tData[, myem.gp.canTOA$Names], tData$BALX,  axis = 7, heading = "st-GP ", xrange = c(-1, 1), ParamNames = myem.gp.canTOA$Names, 
          interval = c(-15, 20))
ValidPlot(LOOs.nst,tData[, myem.gp.canTOA$Names], tData$BALX,  axis = 7, heading = "nst-GP ", xrange = c(-1, 1), ParamNames = myem.gp.canTOA$Names, 
          interval = c(range(tLOOs, LOOs.nst)))


myem.lm.canPCP <- EMULATE.lm(Response="PCP", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 20)
myem.gp.canPCP = EMULATE.gpstan(meanResponse=myem.lm.canPCP, nugget=0.0001, tData=tData, additionalVariables=NULL,CreatePredictObject = FALSE,FastVersion = TRUE)
tLOOs <- stanLOOplot(StanEmulator = myem.gp.canPCP, ParamNames=myem.gp.canPCP$Names)