#Step 1 Source files and load data
source('BuildEmulator/BuildEmulator.R')
source("BuildEmulator/DannyDevelopment.R")
source("BuildEmulator/Rotation.R")
source("HistoryMatching/HistoryMatchingFast.R")

load("Demonstrations/TnSReal.RData")
load("Demonstrations/TSensembleNW1.RData")
plot(RealWorld[,2],-RealWorld[,1],type='l',lwd=2,col=2,xlab="Temperature",ylab="Depth")
depths <- RealWorld[,1]
for(j in 1:400){
  lines(tData[j,23:52],-depths[-31],col=8,lwd=0.5)
}
lines(RealWorld[,2],-RealWorld[,1],lwd=2,col=2)

#Step 2 Extract the fields you want
tData <- tData[1:100,] # only use 100 runs so faster to fit emulators
ORCA2Data <- t(as.matrix(tData[,23:52]))
dim(ORCA2Data)

#Step 3: Centre the data and find the Singular vectors (EOFs)
ORCAsvd <- CentreAndBasis(ORCA2Data)
plot(ORCAsvd$tBasis[,1],-depths[-31],type='l')

#Step 4: Set a discrepancy variance and check the performance of basis reconstructions
Disc <- diag(rep(0.1^2,30))
DiscInv <- GetInverse(Disc) 
attributes(DiscInv)
ObsDat <- RealWorld[-31,2] - ORCAsvd$EnsembleMean
par(mfrow=c(1,2), mar = c(4,4,2,4))
vSVD <- VarMSEplot(ORCAsvd, ObsDat, weightinv = DiscInv, ylim = c(0,20.5), qmax = NULL)
abline(v = which(vSVD[,2] > 0.99)[1])

#Step 5: Rotate the basis for optimal calibration
rotOrca <- RotateBasis(ORCAsvd, ObsDat, kmax = 4, weightinv = DiscInv, v = c(0.3,0.15,0.1,0.1,0.1), vtot = 0.99, MaxTime = 10)
save(rotOrca,file="Demonstrations/rotOrca.RData")

load("Demonstrations/rotOrca.RData")
vROT <- VarMSEplot(rotOrca, ObsDat, weightinv = DiscInv, ylim = c(0,20.5), qmax = NULL)
qROT <- which(vROT[,2] > 0.99)[1]
abline(v = qROT)
qROT

#Step 6: Emulate the basis coefficients and tune emulators!
tDataORCA <- GetEmulatableDataWeighted(Design = tData[,1:20], EnsembleData = rotOrca, HowManyBasisVectors = qROT, weightinv = DiscInv)
tData = tDataORCA
StanEmsO <- InitialBasisEmulators(tData, HowManyEmulators=qROT)
names(StanEmsO[[1]])
StanEmsO[[1]]$StanModel <- NULL
StanEmsO[[2]]$StanModel <- NULL
StanEmsO[[3]]$StanModel <- NULL
save(StanEmsO,file="Demonstrations/StanEmsO.RData")

load("Demonstrations/StanEmsO.RData")
tLOOs1 <- LOO.plot(StanEmulator = StanEmsO[[1]], ParamNames = StanEmsO[[1]]$Names)
tLOOs2 <- LOO.plot(StanEmulator = StanEmsO[[2]], ParamNames = StanEmsO[[2]]$Names)
tLOOs3 <- LOO.plot(StanEmulator = StanEmsO[[3]], ParamNames = StanEmsO[[3]]$Names)

#Step 7: History Matching
#### We can now history match large fields without needing to do large matrix inversions for every x
#### Instead, need a) reconstruction error (fixed for all x) and b) implausibility on basis, with correct projection
###NOTE IF WE WILL RULE OUT ALL OF SPACE, THERE IS A DISCREPANCY SCALING STEP TO AVOID THIS, SEE SPATIALDEMO AND USE CAUTION!

# Predict and history match
FieldHM <- PredictAndHM(rotOrca, ObsDat, StanEmsO, tData, ns = 1000, Error = 0*diag(30), Disc = Disc, weightinv = DiscInv)
summary(FieldHM$impl)

# Plot best runs
par(mfrow=c(1,1), mar = c(4,4,2,2))
plot(RealWorld[,2],-RealWorld[,1],type='l',lwd=2,col=2,xlab="Temperature",ylab="Depth")
for(j in 1:100){
  lines(ORCA2Data[,j],-depths[-31],col=8,lwd=0.5)
}
lines(RealWorld[,2],-RealWorld[,1],type='l',lwd=2,col=2)

inds <- which(FieldHM$impl <= quantile(FieldHM$impl, probs = 0.01)) # plot best 1%
for(k in inds){
  anss <- Reconstruct(FieldHM$Expectation[k,],rotOrca$tBasis[,1:3])+rotOrca$EnsembleMean
  lines(anss,-depths[-31],col=3,lwd=1.5)
}
# Add best
ans <- Reconstruct(FieldHM$Expectation[which.min(FieldHM$impl),],rotOrca$tBasis[,1:3])+rotOrca$EnsembleMean
lines(ans,-depths[-31],col=4,lwd=1.5)


# Create density plot
source("HistoryMatching/HistoryMatching.R")
source("HistoryMatching/impLayoutplot.R")
# Generate a larger design for this (note - NROY space is <1%, should run for a lot longer to get a better plot)
FieldHM2 <- PredictAndHM(rotOrca, ObsDat, StanEmsO, tData, ns = 10000, Error = 0*diag(30), Disc = Disc, weightinv = DiscInv)
ImpData <- cbind(FieldHM2$Design, FieldHM2$impl)
ImpList <- CreateImpList(whichVars = 1:5, VarNames = colnames(tData)[10:14], ImpData, nEms=qROT, Resolution=c(15,15), whichMax=1, Cutoff=FieldHM2$bound)
imp.layoutm11(ImpList,VarNames = colnames(tData)[10:14],VariableDensity=TRUE,newPDF=FALSE,the.title=NULL,newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL)