###################################
#LMDZ Case
twd <- getwd()
source('BuildEmulator/BuildEmulator.R')
source('BuildEmulator/BuildEmulatorNSt.R')
source('HistoryMatching/HistoryMatching.R')
load(file="Demonstrations/lmdzAWave1.RData")
head(tData)
cands <- names(tData)[-which(names(tData)=="SCMdata")]
myem.lm.lmdza1 <- EMULATE.lm(Response="SCMdata", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 20)
myem.gp.lmdza1 = EMULATE.gpstan(meanResponse=myem.lm.lmdza1, tData=tData, additionalVariables=NULL,FastVersion = TRUE)
tLOOs <- LOO.plot(StanEmulator = myem.gp.lmdza1, ParamNames = myem.gp.lmdza1$Names)

#*****NEW***NEW***NEW***#TO CHANGE NUGGET (OR OTHER PRIORS) 
tPriors = gpstan.default.params
#example of changing the nugget, other priors can be tuned
tPriors$nugget <- 0.0002
myem.gp.lmdza1 = EMULATE.gpstan(meanResponse=myem.lm.lmdza1, tData=tData, additionalVariables=NULL,FastVersion = TRUE, prior.params = tPriors)
tLOOs <- LOO.plot(StanEmulator = myem.gp.lmdza1, ParamNames = myem.gp.lmdza1$Names)

load(file="Demonstrations/LESens.RData")
tObs <- LESens[1,48]
tObsErr <- sd(LESens[,48])^2
Disc <- 0

load(file="Demonstrations/lmdzAWave2.RData")
head(tData)
cands <- names(tData)[-which(names(tData)=="SCMdata")]
myem.lm.lmdza2 <- EMULATE.lm(Response="SCMdata", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 20)
myem.gp.lmdza2 = EMULATE.gpstan(meanResponse=myem.lm.lmdza2, nugget=0.0001, tData=tData, additionalVariables=NULL, FastVersion = TRUE)
tLOOs <- LOO.plot(StanEmulator = myem.gp.lmdza2, ParamNames = myem.gp.lmdza2$Names)
Xpred <- 2*randomLHS(10000,6) - 1
Xp <- as.data.frame(Xpred)
names(Xp) <- names(tData)[1:6]
VarNames <- names(Xp)
Timps1 <- UniImplausibilityStan(NewData=Xp, Emulator=myem.gp.lmdza1, Discrepancy=Disc, Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE)
ImpData1 <- cbind(Xp,Timps1)
Xpred2 <- 2*randomLHS(10000,6) - 1
Xp.1 <- as.data.frame(Xpred)
names(Xp.1) <- names(tData)[1:6]
Timps12 <- UniImplausibilityStan(NewData=Xp.1, Emulator=myem.gp.lmdza1, Discrepancy=Disc, Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE)
ImpData12 <- cbind(Xp.1,Timps12)
colnames(ImpData12) <- colnames(ImpData1)
ImpData1 <- rbind(ImpData1,ImpData12)
Timps1 <- c(Timps1,Timps12)
NROY1 <- which(Timps1 <= 3)
Xp <- rbind(Xp,Xp.1)
Xp2 <- Xp[NROY1, ]
dim(Xp2)
Timps2 <- UniImplausibilityStan(NewData=Xp2, Emulator=myem.gp.lmdza2, Discrepancy=Disc, Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE)
ImpData2 <- cbind(Xp2,Timps2)
NROY2 <- which(Timps2 <= 3)
ImpDataM2 <- ImpDataWaveM(Design = Xp, NROY.list = list(NROY1, NROY2), Impl.list = list(Timps1, Timps2))
zz2 <- CreateImpListWaveM(whichVars = 1:6, VarNames= VarNames, ImpData=ImpDataM2, nEms=1, Resolution=c(15,15), whichMax=1)
imp.layoutm11(zz2,VarNames,VariableDensity=FALSE,newPDF=FALSE,the.title="InputSpace.png",newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL)

# CNRM Case
load(file="Demonstrations/cnrmWave1.RData")
head(tData)
cands <- names(tData)[-which(names(tData)=="SCMdata")]
myem.lm.cnrm1 <- EMULATE.lm(Response="SCMdata", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 20)

myem.gp.cnrm1 = EMULATE.gpstan(meanResponse=myem.lm.cnrm1, nugget=0.0001, tData=tData, additionalVariables=NULL,FastVersion = TRUE)
tLOOs <- LOO.plot(StanEmulator = myem.gp.cnrm1, ParamNames = myem.gp.cnrm1$Names)

load(file="Demonstrations/LESens.RData")
tObs <- LESens[1,52]
tObsErr <- sd(LESens[,52])^2
Disc <- 0
Xpred <- 2*randomLHS(10000,13) - 1
Xp <- as.data.frame(Xpred)
names(Xp) <- names(tData)[1:13]
VarNames <- names(Xp)
Timps1 <- UniImplausibilityStan(NewData=Xp, Emulator=myem.gp.cnrm1, Discrepancy=Disc, Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE)
ImpData1 <- cbind(Xp,Timps1)
ImpList <- CreateImpList(whichVars = 1:13, VarNames = VarNames, ImpData = ImpData1,Resolution = c(15,15), whichMax=1)

imp.layoutm11(ImpList,cands[1:5],VariableDensity=FALSE,newPDF=FALSE,the.title="InputSpace.png",newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL)

std.err <- CalStError(LOOs = tLOOs, true.y = tData$SCMdata)
par(mfrow = c(4, 4), mar = c(4, 4, 2, 1))
for(i in 1:(length(cands)-1)) plot(tData[, cands[i]], std.err, pch = 20, xlab = cands[i], ylab = 'std. error')

myem.mixture <- MIXTURE.design(y.std = std.err, X = tData[, cands[-14]], L=2)
mixture.cnrm <- data.frame(cbind(std.err, myem.mixture$MixtureMat, tData[, cands[-14]]))
names(mixture.cnrm) <- c('std.err', 'red', 'blue', cands[-10])
head(mixture.cnrm)
ggplot(mixture.cnrm, aes(x = ALFX, y = std.err, col = rgb(red = red, blue = blue, green = 0), ymin = -3.5, ymax = 3.5)) + geom_point(size = 2) + scale_color_identity()
myem.gp.cnrm.nst1 <- EMULATE.gpstanNSt(meanResponse=myem.lm.cnrm1, tData=tData, additionalVariables= NULL,
                                       FastVersion = TRUE, mixtureComp = myem.mixture)
tLOOs.nst <- LOO.plot.NSt(StanEmulator = myem.gp.cnrm.nst1, ParamNames=cands[-14])

# This code needs to be fixed.
#Timps1.nst <- UniImplausibilityStan(NewData=Xp, Emulator = myem.gp.cnrm.nst1, Discrepancy = Disc, 
#                                    Obs = tObs, ObsErr = tObsErr, is.GP=NULL, FastVersion=TRUE, StatVersion = FALSE, 
#                                    mixtureComp = myem.mixture)
#ImpData1.nst <- cbind(Xp,Timps1.nst)
#ImpList.nst <- CreateImpList(whichVars = 1:13, VarNames = VarNames, ImpData = ImpData1.nst, Resolution = c(15,15), whichMax=1)
#imp.layoutm11(ImpList.nst,cands[1:5],VariableDensity=FALSE,newPDF=FALSE,the.title="InputSpace.png",newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL)