setwd("~/R/ExeterUQ") source("BuildEmulator/DannyDevelopment.R") source("BuildEmulator/JamesNewDevelopment.R") #### Load and format the data #### #setwd("~/PhD/Canada/Wave1/MSLP") #nc000 <- nc_open(file="psl_Amon_DevAM4-4_dma_1-000_r1i1p1_200301-200812.nc") #ourDataPSL <- GetDataAndDesign(DataDir="", prefix="psl", variable="psl") #ourDataPSL$FieldData <- ourDataPSL$FieldData/100 # to MB #JJAmean <- ExtractData(OriginalField = ourDataPSL$FieldData, which.months=c(6,7,8), num.years=5) #JJAmean$Design <- ourDataPSL$Design #lat <- ncvar_get(nc000,"lat") #lon <- ncvar_get(nc000,"lon") #lon[lon>180] <- lon[lon>180] - 360 #JJAmean$lon <- lon #JJAmean$lat <- lat #save(JJAmean, file= "Demonstrations/JJAmean.RData") load("Demonstrations/JJAmean.RData") JJAdataPSL <- CentreAndBasis(JJAmean$Data, scaling=1) lat <- JJAmean$lat lon <- JJAmean$lon plotBasis(Basis=JJAdataPSL$tBasis, 9, lat, lon, pretty(range(JJAdataPSL$tBasis[,1:9]),65)) #### Load observations #### #tPSLFiles <- c("psl_Amon_ERA-Interim_reanalysis_r1i1p1_197901-201510.nc", "psl_Amon_MERRA_reanalysis_r1i1p1_197901-201511.nc", "psl_Amon_NCEP2_reanalysis_r1i1p1_197901-201512.nc") #tObsJJAPSL <- GetObservations(DataFiles=tPSLFiles, which.obs="psl", EnsembleData = JJAdataPSL, permanent.scaling=100) # Already centred load("Demonstrations/tObsJJAPSL.RData") #### Specify some (non-diagonal) discrepancy matrix - SLOW TO SPECIFY, INVERT ETC. FOR l=8192 (1000 OK) #### #LonLatGrid <- expand.grid(lon, lat) #Distances <- rdist.earth(LonLatGrid) #Disc <- exp(-(Distances/20)^2) # We need the inverse of this to calculate the projections in this - SLOW #DiscInv <- GetInverse(Disc) # So just specify a diagonal instead Disc <- diag(runif(8192, 0.1, 3)) # random variances DiscInv <- GetInverse(Disc) # GetInverse also assigns attribute to indicate whether diagonal, so don't need to perform this slow check every time we call other functions (Reconstruction error, CalcScores) - and still fast if have a diagonal attributes(DiscInv) par(mfrow=c(1,2), mar = c(4,4,2,4)) # Set qmax != NULL so don't calculate R_W for every truncated basis (e.g. if l is high) vSVD <- VarMSEplot(JJAdataPSL, tObsJJAPSL[[1]], weightinv = DiscInv, ylim = c(0,60), qmax = 60) abline(v = which(vSVD[,2] > 0.9)[1]) #### Rotate (SLOW if non-diagonal) #### #JJArotPSL <- RotateBasis(JJAdataPSL, tObsJJAPSL[[1]], kmax = 2, weightinv = DiscInv, v = c(0.2,0.15,0.1,0.1,0.1), vtot = 0.9, MaxTime = 60) load("Demonstrations/JJArotPSL.RData") vROT <- VarMSEplot(JJArotPSL, tObsJJAPSL[[1]], weightinv = DiscInv, ylim = c(0,60), qmax = 60) abline(v = which(vROT[,2] > 0.9)[1]) qROT <- which(vROT[,2] > 0.9)[1] #### Emulate #### tDataPSL <- GetEmulatableDataWeighted(Design = JJAmean$Design, EnsembleData = JJArotPSL, HowManyBasisVectors = qROT, weightinv = DiscInv) #emsPSL <- BasisEmulators(tDataPSL, qROT, maxdf = 12) load("Demonstrations/emsPSL.RData")#loads existing emulators. tData = tDataPSL tsts <- BuildStanEmulator(Response="C1", tData, cands=names(tData)[1:14], TryFouriers = TRUE, maxOrder = 2, maxdf = 6) tLOOs <- LOO.plot(StanEmulator = tsts, ParamNames = tsts$Names) StanEms <- InitialBasisEmulators(tData, HowManyEmulators=qROT) tLOOs1 <- LOO.plot(StanEmulator = StanEms[[1]], ParamNames = StanEms[[1]]$Names) tLOOs2 <- LOO.plot(StanEmulator = StanEms[[2]], ParamNames = StanEms[[2]]$Names) tLOOs3 <- LOO.plot(StanEmulator = StanEms[[3]], ParamNames = StanEms[[3]]$Names) tLOOs4 <- LOO.plot(StanEmulator = StanEms[[4]], ParamNames = StanEms[[4]]$Names) tLOOs5 <- LOO.plot(StanEmulator = StanEms[[5]], ParamNames = StanEms[[5]]$Names) tLOOs6 <- LOO.plot(StanEmulator = StanEms[[6]], ParamNames = StanEms[[6]]$Names) tLOOs7 <- LOO.plot(StanEmulator = StanEms[[7]], ParamNames = StanEms[[7]]$Names) tLOOs8 <- LOO.plot(StanEmulator = StanEms[[8]], ParamNames = StanEms[[8]]$Names) tLOOs9 <- LOO.plot(StanEmulator = StanEms[[9]], ParamNames = StanEms[[9]]$Names) tLOOs10 <- LOO.plot(StanEmulator = StanEms[[10]], ParamNames = StanEms[[10]]$Names) tLOOs11 <- LOO.plot(StanEmulator = StanEms[[11]], ParamNames = StanEms[[11]]$Names) tLOOs12 <- LOO.plot(StanEmulator = StanEms[[12]], ParamNames = StanEms[[12]]$Names) tLOOs13 <- LOO.plot(StanEmulator = StanEms[[13]], ParamNames = StanEms[[13]]$Names) #Insert Nonstationary Workflow here. #Assume can just use this emulator list for now. #### Set discrepancy to avoid ruling out all of space #### DiscMultiplier <- SetDiscrepancy(JJArotPSL$tBasis, q = qROT, obs = tObsJJAPSL[[1]], weightinv = DiscInv, level = 0.5) DiscMultiplier <- DiscMultiplier * (JJArotPSL$scaling)^2 DiscScaled <- DiscMultiplier*Disc DiscInvScaled <- (1/DiscMultiplier) * DiscInv # inverse just scales as well VarMSEplot(JJArotPSL, tObsJJAPSL[[1]], weightinv = DiscInvScaled, ylim = c(0,150), qmax = NULL) #### Set a bound so can match using the coefficients #### ImplDesign <- 2*as.data.frame(randomLHS(20, 13)) - 1 colnames(ImplDesign) <- colnames(JJAmean$Design) EmOutput1 <- lapply(1:length(emsPSL), function(e) pred.em(emsPSL[[e]], ImplDesign)) Expectation1 <- matrix(0, nrow = dim(ImplDesign)[1], ncol = length(emsPSL)) Variance1 <- matrix(0, nrow = dim(ImplDesign)[1], ncol = length(emsPSL)) for (j in 1:length(emsPSL)){ Expectation1[,j] <- EmOutput1[[j]]$post.m Variance1[,j] <- EmOutput1[[j]]$post.cov } #Stan version EMULATOR.gpstan(ImplDesign,StanEms[[1]]) EmOutput1 <- lapply(1:length(StanEms), function(e) EMULATOR.gpstan(ImplDesign,StanEms[[e]],FastVersion = TRUE)) Expectation1 <- matrix(0, nrow = dim(ImplDesign)[1], ncol = length(StanEms)) Variance1 <- matrix(0, nrow = dim(ImplDesign)[1], ncol = length(StanEms)) for (j in 1:length(StanEms)){ Expectation1[,j] <- EmOutput1[[j]]$Expectation Variance1[,j] <- EmOutput1[[j]]$Variance } 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(8150,9200), ylim = c(0,800)) abline(v = qchisq(0.995, 8192)) # bound on the field abline(h = qchisq(0.995, length(emsPSL)), lty = 2) # where the chi-squared bound is, based on the number of basis vectors abline(h = T_c) # new bound # Instead of calculating the implausibilities, can load these in and just fit the model #ImplModelInputs <- list(ImplDesign=ImplDesign, Expectation1=Expectation1, Variance1=Variance1, Field1=ImplModel$SampleIf, Coeff1=ImplModel$SampleIc) #save(ImplModelInputs, file = "Demonstrations/ImplModelInputs.RData") load("Demonstrations/ImplModelInputs.RData") ImplModel <- FitImplModel(FieldImpl = ImplModelInputs$Field1, CoeffImpl = ImplModelInputs$Coeff1, l = 8192) #save(ImplModel, file = "Demonstrations/ImplModel.RData") # Or just load the model load("Demonstrations/ImplModel.RData") BoundSamples <- extract(ImplModel, pars = "T_c") stan_dens(ImplModel, pars = "T_c") ModelSummary <- summary(ImplModel, pars = c("sigma_sq", "beta0", "beta1", "alpha0", "alpha1", "T_c"), probs = c(0.995))$summary T_c <- ModelSummary[rownames(ModelSummary) == "T_c", colnames(ModelSummary) == "99.5%"] # Plotting bound, data par(mfrow=c(1,1)) plot(ImplModelInputs$Field1, ImplModelInputs$Coeff1, pch=18, xlim = c(8150,9200), ylim = c(0,800)) abline(v = qchisq(0.995, 8192)) # bound on the field abline(h = qchisq(0.995, length(emsPSL)), lty = 2) # where the chi-squared bound is, based on the number of basis vectors abline(h = T_c) # new bound #### Sample implausiblities #### CoeffHM <- DesignHM <- NULL tempDesign1000 <- 2*as.data.frame(randomLHS(1000, 13)) - 1 colnames(tempDesign1000) <- colnames(JJAmean$Design) EmOutput1000 <- lapply(1:length(emsPSL), function(e) pred.em(emsPSL[[e]], tempDesign1000)) Expectation1000 <- Variance1000 <- matrix(0, nrow = dim(tempDesign1000)[1], ncol = length(emsPSL)) for (j in 1:length(emsPSL)){ Expectation1000[,j] <- EmOutput1000[[j]]$post.m Variance1000[,j] <- EmOutput1000[[j]]$post.cov } CoeffHM <- HistoryMatch(JJArotPSL, type="scoresMV", tObsJJAPSL[[1]], Expectation1000, Variance1000, Error = 0*diag(8192), Disc = DiscScaled, ModelBound = FALSE, weightinv = DiscInvScaled)$impl sum(CoeffHM < T_c) / length(CoeffHM) summary(CoeffHM) which.min(CoeffHM) tcoefstest <- t(t(Expectation1000[755,])) tfield <- Reconstruct(Expectation1000[755,],JJArotPSL$tBasis[,1:13]) plot.EnsembleMemberReconstruction(coefficients=tcoefstest, EnsembleData=JJArotPSL, which.plot=1, OriginalScale=TRUE) CoeffHM <- DesignHM <- NULL tempDesign1000 <- 2*as.data.frame(randomLHS(1000, 13)) - 1 colnames(tempDesign1000) <- colnames(JJAmean$Design) EmOutput1000 <- lapply(1:length(StanEms), function(e) EMULATOR.gpstan(tempDesign1000,StanEms[[e]], FastVersion = TRUE)) Expectation1000 <- Variance1000 <- matrix(0, nrow = dim(tempDesign1000)[1], ncol = length(StanEms)) for (j in 1:length(StanEms)){ Expectation1000[,j] <- EmOutput1000[[j]]$Expectation Variance1000[,j] <- EmOutput1000[[j]]$Variance } CoeffHM <- HistoryMatch(JJArotPSL, type="scoresMV", tObsJJAPSL[[1]], Expectation1000, Variance1000, Error = 0*diag(8192), Disc = DiscScaled, ModelBound = FALSE, weightinv = DiscInvScaled)$impl sum(CoeffHM < T_c) / length(CoeffHM) summary(CoeffHM) which.min(CoeffHM) tcoefstest <- t(Expectation1000[755,]) tfield <- Reconstruct(Expectation1000[755,],JJArotPSL$tBasis[,1:13]) plot.EnsembleMemberReconstruction(coefficients=tcoefstest, EnsembleData=JJArotPSL, which.plot=1, OriginalScale=TRUE) plot.DataEmulator1 <- function(x, Emulator, ObsField, EnsembleData, names.data, which.obs = 1, AnomCols=NULL, AnomBreaks=NULL, add.contours=FALSE, Anomaly.Only=FALSE, tcontours=seq(from=-40,by=2,to=40), contour.zero=TRUE, ...){ colnames(x) <- names.data EmOutput <- lapply(1:length(Emulator), function(e) unlist(EMULATOR.gpstan(x[1,], Emulator[[e]],FastVersion = TRUE))) Coef.exp <- unlist(lapply(EmOutput, function(e) e[1])) Coef.var <- unlist(lapply(EmOutput, function(e) e[2])) pred <- Reconstruct(coefficients=Coef.exp, basis=EnsembleData$tBasis[,1:length(Emulator)]) pred <- pred*EnsembleData$scaling + EnsembleData$EnsembleMean ObsToScale <- ObsField[[which.obs]]*EnsembleData$scaling + EnsembleData$EnsembleMean tanom <- pred - ObsToScale par(mar=c(3,3,4,2)) tdims <- c(length(lon), length(lat)) if(!Anomaly.Only){ par(mfrow=c(2,1), mar=c(3,3,2,2)) image.plot(lon[order(lon)],lat, matrix(pred,nrow=tdims[1],ncol=tdims[2])[order(lon),], ...) map("world",add=T,wrap=TRUE,interior=FALSE) } tanom1 <- tanom tanom1[tanom1>AnomBreaks[length(AnomBreaks)]] <- AnomBreaks[length(AnomBreaks)] tanom1[tanom1AnomBreaks[length(AnomBreaks)]] <- AnomBreaks[length(AnomBreaks)] tanom1[tanom1AnomBreaks[length(AnomBreaks)]] <- AnomBreaks[length(AnomBreaks)] tanom0[tanom0