| 1 | twd <- getwd() |
|---|
| 2 | |
|---|
| 3 | source('BuildEmulator/BuildEmulator.R') |
|---|
| 4 | source('BuildEmulator/BuildEmulatorNSt.R') |
|---|
| 5 | source('HistoryMatching/HistoryMatching.R') |
|---|
| 6 | #============================================= |
|---|
| 7 | # Read results of a series of SCM simulations |
|---|
| 8 | #============================================= |
|---|
| 9 | |
|---|
| 10 | load(file="Wave1_SCM.Rdata") |
|---|
| 11 | print(nparam) |
|---|
| 12 | |
|---|
| 13 | # PRESENTATION STOP ***************************************************************************** |
|---|
| 14 | |
|---|
| 15 | #source('StanEmulateCodeR.R') |
|---|
| 16 | cands <- names(tData)[-(nparam+1)] # Eliminate column which contains the model result |
|---|
| 17 | metric <- names(tData[nparam+1]) |
|---|
| 18 | cands |
|---|
| 19 | |
|---|
| 20 | |
|---|
| 21 | #==================================================================== |
|---|
| 22 | # Scatter plot |
|---|
| 23 | #==================================================================== |
|---|
| 24 | # Ploting for each parameter, the metric as a fuction of parameter values: |
|---|
| 25 | # directly from the SCM outputs |
|---|
| 26 | # Organizing sub-windows as a fuction of the number of plots |
|---|
| 27 | if(nparam<2){ par(mfrow = c(1, 1), mar=c(4, 4, 1, 1)) |
|---|
| 28 | } else if(nparam <=4){ par(mfrow = c(2, 2), mar=c(4, 4, 1, 1)) |
|---|
| 29 | } else if(nparam<=9){ par(mfrow = c(3, 3), mar=c(4, 4, 1, 1)) |
|---|
| 30 | } else if(nparam<=16){ par(mfrow = c(4, 4), mar=c(4, 4, 1, 1)) |
|---|
| 31 | } else if(nparam<=25){ par(mfrow = c(5,5), mar=c(4, 4, 1, 1)) } |
|---|
| 32 | for ( iparam in 1:nparam ) { |
|---|
| 33 | plot(tData[,iparam],tData[,nparam+1],col=2,xlab=cands[iparam],ylab=metric) |
|---|
| 34 | } |
|---|
| 35 | |
|---|
| 36 | #======================================================== |
|---|
| 37 | # fitting metrics (colimn nparam+1 in SCMdata) with a linear model "lm" : |
|---|
| 38 | #======================================================== |
|---|
| 39 | # maxOrder is the max number of fourier modes |
|---|
| 40 | # maxdf : max number of degree of freedom = number of fuctions retained in the linear model |
|---|
| 41 | myem.lm <- EMULATE.lm(Response=names(tData)[nparam+1], tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 20) |
|---|
| 42 | |
|---|
| 43 | #==================== |
|---|
| 44 | # Building emulator : |
|---|
| 45 | #==================== |
|---|
| 46 | # variables 1:nparam correspond to input parameters |
|---|
| 47 | # FastVersion = TRUE. Acceletes the following but less reliable. |
|---|
| 48 | # nugget : white noise. Variance, square of standard deviation with same unit as SCMdata. |
|---|
| 49 | # additionalVariables : |
|---|
| 50 | # The first linear fit may have kept only input parameters x1 and x2 for instance |
|---|
| 51 | # It can be decided or not to add part or all the other variables |
|---|
| 52 | # They are then specified in the additionalVariables |
|---|
| 53 | # if adding x3 to x1 and x2, puting names(tData)[3] or names(tData)[1:3] is equivalent. |
|---|
| 54 | myem.gp <- EMULATE.gpstan(meanResponse=myem.lm, nugget=0.0001, tData=tData, additionalVariables=names(tData)[1:nparam], FastVersion = TRUE) |
|---|
| 55 | # PRESENTATION STOP ******************************************************* |
|---|
| 56 | |
|---|
| 57 | #DIAGNOSTICS : verifying that the original SCMdata are well reproduced when eliminating the point from the emulator |
|---|
| 58 | tLOOs <- LOO.plot(StanEmulator = myem.gp, ParamNames=names(tData)[1:nparam]) |
|---|
| 59 | # PRESENTATION STOP ******************************************************* |
|---|
| 60 | |
|---|
| 61 | |
|---|
| 62 | # Generating 10000 samples varying nparam input parameters |
|---|
| 63 | Xp <- as.data.frame(2*randomLHS(10000, nparam)-1) |
|---|
| 64 | names(Xp) <- names(tData)[1:nparam] |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | ######################## TESTE JUSQUE LA ############################################# |
|---|
| 68 | # Reading observations (here LES) |
|---|
| 69 | load(file="Wave1_LES.Rdata") |
|---|
| 70 | print(paste("tObs",tObs," TobsErr",tObsErr)) |
|---|
| 71 | Disc <- 0 |
|---|
| 72 | |
|---|
| 73 | # Computing implausibility for the 10000 members |
|---|
| 74 | Timps <- UniImplausibilityStan(NewData=Xp, Emulator=myem.gp, Discrepancy=Disc, Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE) |
|---|
| 75 | ImpData <- cbind(Xp,Timps) |
|---|
| 76 | Cutoff <- 3 |
|---|
| 77 | VarNames <- names(Xp) |
|---|
| 78 | # PRESENTATION STOP ******************************************************* |
|---|
| 79 | |
|---|
| 80 | |
|---|
| 81 | ImpList <- CreateImpList(whichVars = 1:nparam, VarNames = VarNames, ImpData = ImpData,Resolution = c(15,15), whichMax=1) |
|---|
| 82 | imp.layoutm11(ImpList,VarNames,VariableDensity=FALSE,newPDF=FALSE,the.title="InputSpace.png",newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL) |
|---|
| 83 | # number of plausibile members divided by total size -> fractin of space retained after wave 1 |
|---|
| 84 | sum(Timps<3)/10000 |
|---|
| 85 | # PRESENTATION STOP ******************************************************* |
|---|
| 86 | |
|---|
| 87 | # Vector NROY |
|---|
| 88 | NROY1 <- which(Timps <= 3) |
|---|
| 89 | # restric Xp to the NROY space -> XP2 |
|---|
| 90 | Xp2 <- Xp[NROY1, ] |
|---|
| 91 | |
|---|
| 92 | # Creation of Wave 2 sample |
|---|
| 93 | design.wave2 <- sample(nrow(Xp2), 80, rep = F) |
|---|
| 94 | design.wave2 |
|---|
| 95 | Wave2 <- Xp2[design.wave2, ] |
|---|
| 96 | save(Wave2,file="Wave2.RData") |
|---|
| 97 | head(Wave2) |
|---|
| 98 | Wave2[,1] |
|---|
| 99 | |
|---|
