twd <- getwd()
load("Demonstrations/DemonstrationsSCMLES.RData")

par(mfrow=c(1,1))
plot(days, SCMens[1,],type='l',col=8,xlab="Time (days)",ylab="Potential Temp Average (th_bl)",ylim=range(SCMens))
for(k in 2:length(SCMens[, 1])){
  lines(days,SCMens[k, ],col=8)
}
for(k in 1:length(LESens[, 1])){
  lines(days,LESens[k,],col=3)
}
source('BuildEmulator/BuildEmulator.R')
source('HistoryMatching/HistoryMatching.R')
### EMULATE SINGLE METRIC SCM-LES
tData <- data.frame(cbind(scaledParams, SCMens[, 72]))
names(tData) <- c(names(scaledParams), "thblsLAST")
cands <- names(tData)[-5]

myem.lm <- EMULATE.lm(Response=names(tData)[5], tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 2)
# Test if my Stan files work with constant mean function
#lin.model <- lm(thblsLAST ~ 1, data = tData)
#myem.lm$linModel <- lin.model
myem.gp <- EMULATE.gpstan(meanResponse=myem.lm, tData=tData, additionalVariables=names(tData)[1:3], FastVersion=TRUE)
### DIAGNOSTICS (LEAVE-ONE-OUT VALIDATIOn)
tLOOs <- LOO.plot(StanEmulator = myem.gp, ParamNames = names(tData)[1:3])
### HISTORY MATCHING SINGLE METRIC SCM-LES
LESruns <- LESens[,72]
Xpred <- 2*randomLHS(10000,3) - 1
Xp <- as.data.frame(Xpred)
names(Xp) <- names(tData)[1:3]
tObs <- LESruns[1]
tObsErr <- sd(LESruns)^2
Disc <- 0
Timps <- UniImplausibilityStan(NewData=Xp, Emulator=myem.gp, Discrepancy=Disc, Obs=tObs, ObsErr=tObsErr, is.GP=NULL, 
                               FastVersion = TRUE, StatVersion = TRUE)
ImpData <- cbind(Xp,Timps)
Cutoff <- 3
VarNames <- names(Xp)
ImpList <- CreateImpList(whichVars = 1:3, VarNames = VarNames, ImpData = ImpData,Resolution = c(15,15), whichMax=1)
imp.layoutm11(ImpList,VarNames,VariableDensity=FALSE,newPDF=FALSE,the.title="InputSpace.png",newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL)

# EMULATE TIME-SERIES
scmData <- ExtractCentredScaledDataAndBasis(OriginalField=t(SCMens), scaling=1)
scmData$tBasis
par(mfrow=c(1, 3), mar=c(4, 4, 2, 1))
plot(1:72,scmData$tBasis[,1],type='l',col=1, xlab = 'days')
plot(1:72,scmData$tBasis[,2],type='l',col=1, xlab = 'days')
plot(1:72,scmData$tBasis[,3],type='l',col=1, xlab = 'days')
tweight <- var(LESens)
VarMSEplot(scmData, LESens[1,]-scmData$EnsembleMean, weight = diag(diag(tweight)))
par(mfrow=c(1,2),mar=c(4,4,4,4))
VarMSEplot(scmData, LESens[1,]-scmData$EnsembleMean, weight = diag(0.001,72),ylim=c(1,5))
scmRotated <- RotateIterate(k=30, DataBasis=scmData, obs=LESens[1,]-scmData$EnsembleMean, weight = diag(0.001,72))
VarMSEplot(scmRotated, LESens[1,]-scmData$EnsembleMean, weight = diag(0.001,72),ylim=c(1,5),min.line = FALSE)
tData <- GetEmulatableData(Design=scaledParams[, 1:3], EnsembleData=scmRotated, HowManyBasisVectors=4, Noise=TRUE)
EmList <- InitialBasisEmulators(tData=tData, HowManyEmulators=4, additionalVariables=list(names(tData)[1:3],names(tData)[1:3],names(tData)[1:3],names(tData)[1:3]))
tLOOs1 <- LOO.plot(StanEmulator = EmList[[1]], ParamNames=names(tData)[1:3])
tLOOs1 <- LOO.plot(StanEmulator = EmList[[2]], ParamNames=names(tData)[1:3])
tLOOs1 <- LOO.plot(StanEmulator = EmList[[3]], ParamNames=names(tData)[1:3])
Xpred <- 2*randomLHS(10000,3) - 1
Xp <- as.data.frame(Xpred)
names(Xp) <- names(tData)[1:3]
EMpreds <- lapply(1:length(EmList), function(k) EMULATOR.gpstan.multicore(Xp, Emulator=EmList[[k]],FastVersion = TRUE,batches = 500))
tObs <- LESens[1,] - scmData$EnsembleMean
tObsErr <- diag(diag(var(LESens)))
tDisc <- tObsErr
tDisc[1:25,1:25] <- 10*tDisc[1:25,1:25]
TMimps <- BasisImplausibility(Basis=scmRotated$tBasis[,1:4], EmPreds=EMpreds, Obs=tObs, Error=tObsErr, Disc=tDisc)

tBound <- qchisq(0.995,length(tObs))
sum(TMimps>=tBound)
par(mfrow=c(1, 1), mar=c(4, 4, 2, 1))
plot(days,SCMens[1,],type='l',col=8,xlab="Time (days)",ylab="Potential Temp Average (th_bl)",ylim=range(SCMens))
for(k in 2:dim(SCMens)[1]){
  lines(days,SCMens[k,],col=8)
}
for(k in 1:length(LESens[,1])){
  lines(days,LESens[k,],col=3)
}
addReconLine <- function(index){
  ex <- lapply(EMpreds, function(e) e$Expectation[index])
  Recon(ex,scmRotated$tBasis[,1:4]) + scmData$EnsembleMean
  lines(days,Recon(ex,scmRotated$tBasis[,1:4]) + scmData$EnsembleMean,col=4)
}
addReconLine(which.min(TMimps))