source: src/htune_EOF.R @ 49

#Step 1 Source files and load data
source('BuildEmulator/BuildEmulator.R')
source("BuildEmulator/DannyDevelopment.R")
source("BuildEmulator/JamesNewDevelopment.R")
library("ncdf4") # to manipulate ncdf
#load("Demonstrations/TnSReal.RData")

#########################
# First settings

case_name="AYOTTE"
subcase_name="24SC"
variable="theta" # qv ou rneb # Don't use "var" which is un R function !!!

WAVEN=1

source('htune_case_setup.R')
file = paste("WAVE",WAVEN,"/Wave",WAVEN,".RData",sep="")
load(file)

casesu <-case_setup(case_name)
NLES=casesu[1]
TimeLES=casesu[2]
TimeSCM=casesu[3]

NRUNS=dim(wave_param_US)[1]
nparam=dim(wave_param_US)[2]-1

#########################
# Get vertical grid to project everything
print("Get vertical grid")

file=paste("WAVE",WAVEN,"/",case_name,"/",subcase_name,"/SCM_",WAVEN,"-001.nc",sep="")
print(file)
nc =  nc_open(file)
zf = ncvar_get(nc,"zf")
zfSCM = zf[,1]
nlev = length(zfSCM)
nc_close(nc)

# Get vertical grid of LES
file = paste("LES/",case_name,"/",subcase_name,"/LES",0,".nc",sep="")
nc =  nc_open(file)
zfLES = ncvar_get(nc,"zf")
nc_close(nc)

# Extract the SCM vertical below the top of the LES
zgrid = rep(0,nlev)
ii = 1
maxzLES = max(zfLES)
minzLES = min(zfLES)
for (k in 1:nlev) {
  if ( zfSCM[k] <= maxzLES & zfSCM[k] >= minzLES ) {
    zgrid[ii] = zfSCM[k]
    ii = ii +1
  }
}
zgrid=zgrid[1:ii-1]
nlevout = length(zgrid)[1]

print("Chosen vertical grid:")
print(zgrid)

#########################
# Interpolate LES data on the chosen vertical grid
print("Interpolate LES data on the chosen vertical grid")

dataLES = matrix(0,ncol=NLES+2,nrow=nlevout)
dataLES[,1] = zgrid

for ( k in 0:(NLES-1) ) {
  file = paste("LES/",case_name,"/",subcase_name,"/LES",k,".nc",sep="")
  nc =  nc_open(file)
  zf = ncvar_get(nc,"zf")
  data = ncvar_get(nc,variable)
  nc_close(nc)

  tempFun <- approxfun(zf,data[,TimeLES],yleft=data[1,TimeLES])
  dataLES[,k+2] <- tempFun(zgrid)  
  
  if ( k == 0 ) {
    # Check the interpolation
    plot(x=dataLES[,k+2],y=dataLES[,1],type='l')
    lines(x=data[,TimeLES],y=zf,col="red") 
  }
}

# Compute observation uncertainty
for ( k in 1:nlevout ) {
  dataLES[k,NLES+2] = sd(dataLES[k,2:NLES+1])^2
}

save(dataLES,file="dataLES.RData")
load("dataLES.RData")

#########################
# Interpolate SCM data on the chosen vertical grid
print("Interpolate SCM data on the chosen vertical grid")

dataSCM = matrix(0,ncol=NRUNS+1,nrow=nlevout)
dataSCM[,1] = zgrid

for ( k in 1:NRUNS ) {
  i = sprintf("%03i",k)
  file=paste("WAVE",WAVEN,"/",case_name,"/",subcase_name,"/SCM_",WAVEN,"-",i,".nc",sep="")
  nc =  nc_open(file)
  zf = ncvar_get(nc,"zf")
  data = ncvar_get(nc,variable)
  nc_close(nc)
  nlev = dim(zf)[1]
  
  # yright, yleft to constrain what happen at the edge. Not clear yet how to do it properly
  tempFun <- approxfun(zf[,TimeSCM],data[,TimeSCM],yright=data[1,TimeSCM],yleft=data[nlev,TimeSCM])
  dataSCM[,k+1] <- tempFun(zgrid)  
  
  if ( k == 0 ) {
    # Check the interpolation
    plot(x=dataSCM[,k+2],y=dataSCM[,1],type='l')
    lines(x=data[,TimeSCM],y=zf,col="red") 
  }
  
}

save(dataSCM,file="dataSCM.RData")
load("dataSCM.RData")

# For memory yet
#nt=dim(data)[2]
#dataout = matrix(0,ncol=nt,nrow=nlevout)
#for ( it in 1:nt){
#  tempFun <- approxfun(zf,data[,it],yleft=data[1,it])
#  dataout[,it] <- tempFun(zgrid)  
#}

#load("Demonstrations/TSensembleNW1.RData")
plot(x=dataLES[,2],y=dataLES[,1],type='l',lwd=2,col=2,xlab=variable,ylab="Altitude")
depths <- dataLES[,1]
for(j in 1:NRUNS){
  lines(x=dataSCM[,j],y=depths,col=8,lwd=0.5)
}
lines(x=dataLES[,2],y=dataLES[,1],lwd=2,col=2)
#Step 2 Extract the fields you want
SCMdata <- dataSCM[,-1]
#ORCA2Data <- t(as.matrix(tData[,23:52]))
#dim(ORCA2Data)

#Step 3: Centre the data and find the Singular vectors (EOFs)
SCMsvd <- CentreAndBasis(SCMdata)
par(mfrow=c(3,3))
for ( i in 1:9 ) {
  plot(x=SCMsvd$tBasis[,i],y=dataSCM[,1],type='l')
}

#Step 4: Set a discrepancy variance and check the performance of basis reconstructions
# To be thought about
tmp = rep(0.01,nlevout)
#for ( i in 1:nlevout ) {
#  if ( dataSCM[i,1] >= 2500. ) {tmp[i] = 1}
#}

#for ( i in 1:nlevout ) {
#  if ( dataSCM[i,1] >= 2000. ) {tmp[i] = 0.1}
#}

#Disc <- diag(0.01,nlevout)
Disc <- diag(tmp)
DiscInv <- GetInverse(Disc) 
attributes(DiscInv)
ObsDat <- dataLES[,2] - SCMsvd$EnsembleMean
par(mfrow=c(1,2), mar = c(4,4,2,4))
vSVD <- VarMSEplot(SCMsvd, ObsDat, weightinv = DiscInv, ylim = c(0,80), qmax = NULL)
abline(v = which(vSVD[,2] > 0.9)[1])

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

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

#Step 6: Emulate the basis coefficients and tune emulators!
tDataSCM <- GetEmulatableDataWeighted(Design = wave_param_US[,-1], EnsembleData = rotSCM, HowManyBasisVectors = qROT, weightinv = DiscInv)
tData = tDataSCM
StanEmsO <- InitialBasisEmulators(tData, HowManyEmulators=qROT)
# If it seems not working, relaunch source('BuildEmulator/BuildEmulator.R')
# A fast check that it worked
names(StanEmsO[[1]])

# To remove heavy data, especially if you want to save the emulators
#for ( i in 1:qROT ) {
#  StanEmsO[[i]]$StanModel <- NULL
#}
#save(StanEmsO,file="Demonstrations/StanEmsO.RData")

#load("Demonstrations/StanEmsO.RData")
for ( i in 1:qROT ) {
  tLOOs1 <- LOO.plot(StanEmulator = StanEmsO[[i]], ParamNames = StanEmsO[[i]]$Names)
}

#Step 7: History Matching preliminaries. If your emulator is good enough to history match, you don't always want to do it on the whole field as that inverts massive matrices everywhere in parameter space! We now model the relationship between field and coefficient implausibilities.
###NOTE IF WE WILL RULE OUT ALL OF SPACE, THERE IS A DISCREPANCY SCALING STEP TO AVOID THIS, SEE SPATIALDEMO AND USE CAUTION!
ImplDesign <- 2*as.data.frame(randomLHS(nparam, nparam)) - 1
colnames(ImplDesign) <- colnames(tData)[1:nparam]
EmOutput1 <- lapply(1:length(StanEmsO), function(e) EMULATOR.gpstan(ImplDesign,StanEmsO[[e]],FastVersion = TRUE))
Expectation1 <- matrix(0, nrow = dim(ImplDesign)[1], ncol = length(StanEmsO))
Variance1 <- matrix(0, nrow = dim(ImplDesign)[1], ncol = length(StanEmsO))
for (j in 1:length(StanEmsO)){
  Expectation1[,j] <- EmOutput1[[j]]$Expectation
  Variance1[,j] <- EmOutput1[[j]]$Variance
}

#ImplModel <- HMCoeffBound(JJArotSAT, tObsJJASAT[[1]], Expectation1, Variance1, Error = 0*diag(8192), Disc = DiscScaled, weightinv = DiscInvScaled) # Error set as 0 matrix as combined into Discrepancy
#With raw discrepancy (in case emulator uncertainty is large)
ImplModel <- HMCoeffBound(rotSCM, ObsDat, Expectation1, Variance1, Error = 0*diag(nlevout), Disc = Disc, weightinv = DiscInv)
#save(ImplModel, file="Demonstrations/ImplModel.RData")

#load("Demonstrations/ImplModel.RData")
T_c <- ImplModel$Bound
stan_dens(ImplModel$model, pars = "T_c")
par(mfrow=c(1,1))
plot(ImplModel$SampleIf, ImplModel$SampleIc, pch=18, xlim = c(0,max(ImplModel$SampleIf)), ylim = c(0,max(ImplModel$SampleIc)))
abline(v = qchisq(0.995, 30)) # bound on the field
abline(h = qchisq(0.995, length(StanEmsO)), lty = 2) # where the chi-squared bound is, based on the number of basis vectors
abline(h = T_c)

#Final step: History matching in the coefficient space. 
CoeffHM <- DesignHM <- NULL
tempDesign1000 <- 2*as.data.frame(randomLHS(1000, 20)) - 1
colnames(tempDesign1000) <- colnames(tData)[1:20]
EmOutput1000 <- lapply(1:length(StanEmsO), function(e) EMULATOR.gpstan(tempDesign1000,StanEmsO[[e]], FastVersion = TRUE))
Expectation1000 <- Variance1000 <- matrix(0, nrow = dim(tempDesign1000)[1], ncol = length(StanEmsO))
for (j in 1:length(StanEmsO)){
  Expectation1000[,j] <- EmOutput1000[[j]]$Expectation
  Variance1000[,j] <- EmOutput1000[[j]]$Variance
}
CoeffHM <- HistoryMatch(rotSCM, type="scoresMV", ObsDat, Expectation1000, Variance1000, Error = 0*diag(nlevout), Disc = Disc, ModelBound = FALSE, weightinv = DiscInv)$impl
sum(CoeffHM < T_c) / length(CoeffHM)
summary(CoeffHM)

ans <- Reconstruct(Expectation1000[which.min(CoeffHM),],rotSCM$tBasis[,1:qROT])+rotSCM$EnsembleMean
plot(x=dataLES[,2],y=dataLES[,1],type='l',lwd=2,col=2,xlab=variable,ylab="Altitude")
for(j in 2:NRUNS+1){
  lines(x=dataSCM[,j],y=dataSCM[,1],col=8,lwd=0.5)
}
lines(x=dataLES[,2],y=dataLES[,1],type='l',lwd=2,col=2)
for(k in 1:10){
anss <- Reconstruct(Expectation1000[k,],rotSCM$tBasis[,1:qROT])+rotSCM$EnsembleMean
lines(x=anss,y=dataSCM[,1],col=3,lwd=1.5)
}
lines(x=ans,y=dataSCM[,1],col=4,lwd=1.5)