| 1 | twd <- getwd() |
|---|
| 2 | source('htune_convert.R') |
|---|
| 3 | |
|---|
| 4 | #============================================= |
|---|
| 5 | # default values for optional arguments |
|---|
| 6 | #============================================= |
|---|
| 7 | WAVEN=1 |
|---|
| 8 | tau=0 |
|---|
| 9 | cutoff=3 |
|---|
| 10 | sample_size=1000000 |
|---|
| 11 | sample_size=40000000 |
|---|
| 12 | sample_size=5000000 |
|---|
| 13 | sample_size_next_design=80 |
|---|
| 14 | |
|---|
| 15 | #============================================= |
|---|
| 16 | # Reading line args |
|---|
| 17 | #============================================= |
|---|
| 18 | args = commandArgs(trailingOnly=TRUE) |
|---|
| 19 | if (length(args)%%2!=0) { print(paste("Bad number of arguments",length(args),"in htune_EmulatingMultiMetric.R")); q("no",1) ; } |
|---|
| 20 | if (length(args)>0) { |
|---|
| 21 | for (iarg in seq(1,length(args),by=2)) { |
|---|
| 22 | if (args[iarg]=="-wave") { WAVEN=as.numeric(args[iarg+1]) } |
|---|
| 23 | else if (args[iarg]=="-tau") {tau=as.numeric(args[iarg+1]) } |
|---|
| 24 | else if (args[iarg]=="-cutoff") {cutoff=as.numeric(args[iarg+1]) } |
|---|
| 25 | else { print(paste("Bad argument",args[iarg],"in htune_EmulatingMultiMetric.R")) ; q("no",1) ; } # quit with error code 1 |
|---|
| 26 | } |
|---|
| 27 | } |
|---|
| 28 | |
|---|
| 29 | #=============================================== |
|---|
| 30 | # Reading results of a series of SCM simulations |
|---|
| 31 | #=============================================== |
|---|
| 32 | print(paste("Arguments : -wave",WAVEN,"-cutoff",cutoff,"-tau",tau,sep=" ")) |
|---|
| 33 | source(paste("WAVE", WAVEN, "/ModelParam.R", sep="")) |
|---|
| 34 | load(file=paste("WAVE",WAVEN,"/Wave",WAVEN,"_SCM.Rdata",sep="")) |
|---|
| 35 | load(file=paste("WAVE",WAVEN,"/Wave",WAVEN,"_REF.Rdata",sep="")) |
|---|
| 36 | |
|---|
| 37 | if ( nmetrique == 1) { tau=0 } |
|---|
| 38 | |
|---|
| 39 | cands <- names(tData)[1:(nparam+1)] # Selects columns with the parameters |
|---|
| 40 | metric <- names(tData)[(nparam+2):(nparam+1+nmetrique)] # select columns with the metrics |
|---|
| 41 | print("History Matching for parameters: ") |
|---|
| 42 | print(cands[1:(nparam+1)]) |
|---|
| 43 | print(c("nmetrique,tau=",nmetrique,tau)) |
|---|
| 44 | |
|---|
| 45 | #==================================================================== |
|---|
| 46 | # Scatter plot |
|---|
| 47 | #==================================================================== |
|---|
| 48 | # Ploting for each parameter, the metric as a fuction of parameter values: |
|---|
| 49 | # directly from the SCM outputs |
|---|
| 50 | # Organizing sub-windows as a fuction of the number of plots |
|---|
| 51 | # For ploting with real prameters as axis |
|---|
| 52 | |
|---|
| 53 | file = paste("WAVE",WAVEN,"/Plots_Metrics.pdf",sep="") |
|---|
| 54 | pdf(file=file) |
|---|
| 55 | |
|---|
| 56 | param_SCM <- DesignConvert(tData[,1:nparam]) |
|---|
| 57 | logs=array(1:nparam) ; for ( iparam in 1:nparam ) { logs[iparam]<-"" } |
|---|
| 58 | logs[which.logs] <- "x" |
|---|
| 59 | if(nparam*nmetrique<2){ par(mfrow = c(1, 1), mar=c(4, 4, 1, 1)) |
|---|
| 60 | } else if(nparam*nmetrique <=4){ par(mfrow = c(2, 2), mar=c(4, 4, 1, 1)) |
|---|
| 61 | } else if(nparam*nmetrique<=9){ par(mfrow = c(3, 3), mar=c(4, 4, 1, 1)) |
|---|
| 62 | } else if(nparam*nmetrique<=16){ par(mfrow = c(4, 4), mar=c(4, 4, 1, 1)) |
|---|
| 63 | } else if(nparam*nmetrique<=25){ par(mfrow = c(5,5), mar=c(4, 4, 1, 1)) } |
|---|
| 64 | for ( iparam in 1:nparam ) { |
|---|
| 65 | for (j in 1:nmetrique) { |
|---|
| 66 | ymin=tObs[j]-sqrt(tObsErr[j]) |
|---|
| 67 | ymax=tObs[j]+sqrt(tObsErr[j]) |
|---|
| 68 | plot(param_SCM[,iparam],tData[,nparam+j+1],col=2,xlab=cands[iparam],ylab=metric[j],log=logs[iparam],ylim=range(tData[,nparam+j+1],ymin,ymax)) |
|---|
| 69 | abline(v=param.defaults[iparam],col='blue',lty=2) |
|---|
| 70 | abline(h=c(tObs[j],ymin,ymax ), col = 'blue', lty =2, lwd=c(1,3,3)) |
|---|
| 71 | } |
|---|
| 72 | } |
|---|
| 73 | dev.off() |
|---|
| 74 | |
|---|
| 75 | print(paste("Created figure file: ", file, sep="")) |
|---|
| 76 | |
|---|
| 77 | |
|---|
| 78 | ################################################################# |
|---|
| 79 | # Loading emulator for previous waves if it exists |
|---|
| 80 | ################################################################# |
|---|
| 81 | source('BuildEmulator/BuildEmulator.R') |
|---|
| 82 | options(mc.cores = parallel::detectCores()) |
|---|
| 83 | EMULATOR.LIST <- list() |
|---|
| 84 | tau_vec = c() |
|---|
| 85 | cutoff_vec = c() |
|---|
| 86 | |
|---|
| 87 | if(file.exists("EMULATOR_LIST_MULT_METRIC.RData")) { load("EMULATOR_LIST_MULT_METRIC.RData") } |
|---|
| 88 | |
|---|
| 89 | if ( length(EMULATOR.LIST) == WAVEN ) { |
|---|
| 90 | # We shall ovrewride the cutoff or tau for the last wave |
|---|
| 91 | tau_vec[WAVEN]=tau |
|---|
| 92 | cutoff_vec[WAVEN]=cutoff |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | } else { |
|---|
| 96 | ################################################################# |
|---|
| 97 | # Building emulator for the new WAVE (if not already done) |
|---|
| 98 | ################################################################# |
|---|
| 99 | # variables 1:nparam correspond to input parameters |
|---|
| 100 | # FastVersion = TRUE. Acceletes the following but less reliable. |
|---|
| 101 | # nugget : white noise. Variance, square of standard deviation with same unit as SCMdata. |
|---|
| 102 | # additionalVariables : |
|---|
| 103 | # The first linear fit may have kept only input parameters x1 and x2 for instance |
|---|
| 104 | # It can be decided or not to add part or all the other variables |
|---|
| 105 | # They are then specified in the additionalVariables |
|---|
| 106 | # if adding x3 to x1 and x2, puting names(tData)[3] or names(tData)[1:3] is equivalent. |
|---|
| 107 | # !!! CAUTION !!! InitialBasisEmulators has linear fit +gpstan |
|---|
| 108 | #======================================================== |
|---|
| 109 | # First fitting metrics (columns nparam+2:nparam+1+nmetrique of tData from file SCM.Rdata) with a linear model "lm" : |
|---|
| 110 | #======================================================== |
|---|
| 111 | # maxOrder is the max number of fourier modes |
|---|
| 112 | # maxdf : max number of degree of freedom = number of fuctions retained in the linear model |
|---|
| 113 | # This should be adapted as a function of the number of simulations |
|---|
| 114 | ################################################################# |
|---|
| 115 | |
|---|
| 116 | listnew <-lapply(1:nmetrique,function(k) names(tData[1:nparam])) # Est ce que ça fait vraiment ce qu'on veut ? |
|---|
| 117 | myem.lm <-InitialBasisEmulators(tData=tData,HowManyEmulators=nmetrique,additionalVariables = listnew) |
|---|
| 118 | for(i in 1:nmetrique) myem.lm[[i]]$StanModel = NULL |
|---|
| 119 | file = paste("WAVE",WAVEN,"/Plots_LOO.pdf",sep="") |
|---|
| 120 | pdf(file=file) |
|---|
| 121 | #DIAGNOSTICS : verifying that the original SCMdata are well reproduced when eliminating the point from the emulator |
|---|
| 122 | for (i in 1:nmetrique) { |
|---|
| 123 | tLOOs1 <- LOO.plot(StanEmulator = myem.lm[[i]], ParamNames=names(tData)[1:nparam]) |
|---|
| 124 | } |
|---|
| 125 | dev.off() |
|---|
| 126 | |
|---|
| 127 | EMULATOR.LIST[[WAVEN]] = myem.lm |
|---|
| 128 | tau_vec[WAVEN] = tau |
|---|
| 129 | #SAVE EMULATOR FOR HISTORY MATCHING |
|---|
| 130 | cutoff_vec[WAVEN] = cutoff |
|---|
| 131 | save(EMULATOR.LIST, tau_vec, cutoff_vec, file = "EMULATOR_LIST_MULT_METRIC.RData") |
|---|
| 132 | print(paste("A list of emulators has been saved under: ","EMULATOR_LIST_MULT_METRIC.RData",sep="")) |
|---|
| 133 | } |
|---|
| 134 | |
|---|
| 135 | |
|---|
| 136 | ####################################################################################################### |
|---|
| 137 | # Starting history matching or iterative refocusing |
|---|
| 138 | ####################################################################################################### |
|---|
| 139 | |
|---|
| 140 | source('HistoryMatching/HistoryMatching.R') |
|---|
| 141 | # Generating 10000 samples varying nparam input parameters |
|---|
| 142 | Xp <- as.data.frame(2*randomLHS(sample_size, nparam)-1) |
|---|
| 143 | names(Xp) <- names(tData)[1:nparam] |
|---|
| 144 | Disc = rep(0, nmetrique) |
|---|
| 145 | |
|---|
| 146 | totonew<-DesignantiConvert1D(param.defaults) |
|---|
| 147 | param.defaults.norm=rep(0,nparam) |
|---|
| 148 | nbatches=500 |
|---|
| 149 | for (i in 1:nparam) {param.defaults.norm[i]=totonew[i,]} |
|---|
| 150 | |
|---|
| 151 | if(WAVEN == 1) { |
|---|
| 152 | Timps <- ManyImplausibilitiesStan(NewData=Xp, Emulator=EMULATOR.LIST[[1]], Discrepancy=Disc, |
|---|
| 153 | Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE, |
|---|
| 154 | multicore=(ceiling(dim(Xp)[1]/nbatches)>1), batches=nbatches) # Multicore only if Xp is large enough |
|---|
| 155 | ImpData_wave1 = cbind(Xp, Timps) |
|---|
| 156 | VarNames <- names(Xp) |
|---|
| 157 | valmax = tau_vec[1] + 1 |
|---|
| 158 | ImpListM1 = CreateImpList(whichVars = 1:nparam, VarNames=VarNames, ImpData=ImpData_wave1, nEms=length(EMULATOR.LIST[[1]]), whichMax=valmax) |
|---|
| 159 | imp.layoutm11(ImpListM1,VarNames,VariableDensity=FALSE,newPDF=TRUE,the.title=paste("InputSpace_wave",WAVEN,".pdf",sep=""),newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=matrix(param.defaults.norm,ncol=nparam)) |
|---|
| 160 | |
|---|
| 161 | NROY1 <- which(rowSums(Timps <= cutoff_vec[1]) >=length(EMULATOR.LIST[[1]]) - tau_vec[1]) |
|---|
| 162 | TMimpls_wave1 <- apply(Timps, 1,MaxImp,whichMax=valmax) |
|---|
| 163 | XpNext <- Xp[NROY1, ] |
|---|
| 164 | mtext(paste("Remaining space:",length(NROY1)/dim(Xp)[1],sep=""), side=1) |
|---|
| 165 | |
|---|
| 166 | |
|---|
| 167 | # number of plausible members divided by total size -> fraction of space retained after wave N |
|---|
| 168 | print(paste("Remaining space after wave 1: ",length(NROY1)/dim(Xp)[1],sep="")) |
|---|
| 169 | cat(length(NROY1)/dim(Xp)[1],file="Remaining_space_after_wave_1.txt", sep="") |
|---|
| 170 | print("que fait on la") |
|---|
| 171 | |
|---|
| 172 | } else { |
|---|
| 173 | XpNext = Xp |
|---|
| 174 | NROY.list = list() |
|---|
| 175 | Impl.list = list() |
|---|
| 176 | for(i in 1:(length(EMULATOR.LIST))) { |
|---|
| 177 | Timps = ManyImplausibilitiesStan(NewData=XpNext, Emulator=EMULATOR.LIST[[i]], Discrepancy=Disc, |
|---|
| 178 | Obs=tObs, ObsErr=tObsErr, is.GP=NULL,FastVersion = TRUE, |
|---|
| 179 | multicore=(ceiling(dim(XpNext)[1]/nbatches)>1), batches=nbatches) # Multicore only if XpNext is large enough |
|---|
| 180 | valmax = tau_vec[i] + 1 |
|---|
| 181 | print(c("cutoff",i,cutoff_vec[i])) |
|---|
| 182 | Impl.list[[i]] = matrix(apply(Timps, 1,MaxImp,whichMax=valmax), ncol = 1) |
|---|
| 183 | NROY.list[[i]] = which(rowSums(Timps <= cutoff_vec[i]) >= length(EMULATOR.LIST[[i]]) - tau_vec[i]) |
|---|
| 184 | XpNext = XpNext[NROY.list[[i]], ] |
|---|
| 185 | print(paste("Remaining space after wave",i,": ",length(NROY.list[[i]])/dim(Xp)[1],sep="")) |
|---|
| 186 | cat(length(NROY.list[[i]])/dim(Xp)[1],file=paste("Remaining_space_after_wave_",i,".txt", sep=""),sep="") |
|---|
| 187 | } |
|---|
| 188 | Impl.list[[length(EMULATOR.LIST)]] = Timps |
|---|
| 189 | ImpData <- ImpDataWaveM(Xp, NROY.list, Impl.list) |
|---|
| 190 | VarNames <- names(Xp) |
|---|
| 191 | ImpList <- CreateImpListWaveM(whichVars = 1:nparam, VarNames=VarNames, ImpData = ImpData, |
|---|
| 192 | nEms=length(EMULATOR.LIST[[length(EMULATOR.LIST)]]), Resolution=c(15,15), whichMax=valmax) |
|---|
| 193 | imp.layoutm11(ImpList,VarNames,VariableDensity=FALSE,newPDF=TRUE,the.title=paste("InputSpace_wave",WAVEN,".pdf",sep=""),newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=matrix(param.defaults.norm,ncol=nparam)) |
|---|
| 194 | mtext(paste("Remaining space:",length(NROY.list[[length(EMULATOR.LIST)]])/dim(Xp)[1],sep=""), side=1) |
|---|
| 195 | } |
|---|
| 196 | |
|---|
| 197 | print(paste("Created figure file: ", "InputSpace_wave",WAVEN,".pdf", sep="")) |
|---|
| 198 | |
|---|
| 199 | |
|---|
| 200 | ################################################################################ |
|---|
| 201 | # Creation of Wave WAVEN+1 sample |
|---|
| 202 | ################################################################################ |
|---|
| 203 | |
|---|
| 204 | if (nrow(XpNext)>=sample_size_next_design) { samplesz=sample_size_next_design } else { samplesz=nrow(XpNext) |
|---|
| 205 | print(paste("Final Sample size reduced to ",samplesz)) |
|---|
| 206 | } |
|---|
| 207 | |
|---|
| 208 | design.waveNext <- sample(nrow(XpNext), samplesz, rep = F) |
|---|
| 209 | WaveNext <- XpNext[design.waveNext, ] |
|---|
| 210 | save(WaveNext,file=paste("Wave",strtoi(WAVEN, base = 0L)+1,".RData",sep="")) |
|---|
| 211 | print(paste("Next wave design has been saved under: ","Wave",strtoi(WAVEN, base = 0L)+1,".RData",sep="")) |
|---|
| 212 | #} else { print("The NROY space is not large enough to allow resampling for next wave") } |
|---|
