update README

Author: ChihLi

README.md

@@ -1,17 +1,24 @@
-This folder consists of the R code for the manuscript “*Stacking designs: designing multi-fidelity computer experiments with confidence*” by Sung, Ji, Tang, and Mak. 
+This instruction aims to reproduce the results in the paper “*Stacking
+designs: designing multi-fidelity computer experiments with confidence*”
+by [Sung, Ji, Mak, Wang, and Tang
+(2023+)](https://arxiv.org/abs/2211.00268), recently accepted by
+SIAM/ASA Journal on Uncertainty Quantification.
 
-The code folder reproduces the result in Section 5 of the manuscript. 
+The code folder reproduces the result in Section 5 of the manuscript.
 
-- `currin_example.R` reproduces the results of Section 5.1. 
-- `poisson_example.R` reproduces the results of Section 5.2. 
-- `blade_example.R` reproduces the results of Section 5.3. 
-	- Note: this code could take a long time because it runs very high fidelity simulations.
-	
+-   `currin_example.R` reproduces the results of Section 5.1.
+-   `poisson_example.R` reproduces the results of Section 5.2.
+-   `blade_example.R` reproduces the results of Section 5.3.
 
 The following R packages are required to run the code:
-  - `matlabr (>=1.5.2)`
-  - `randtoolbox (>=2.0.2)`
-  - `RColorBrewer (>=1.1-3)`
-  - `ggplot2 (>=3.3.6)`
-  - `gridExtra (>=2.3)`
-  - `plgp (>=1.1-11)`
+
+-   `matlabr (>=1.5.2)`
+-   `randtoolbox (>=2.0.2)`
+-   `RColorBrewer (>=1.1-3)`
+-   `ggplot2 (>=3.3.6)`
+-   `gridExtra (>=2.3)`
+-   `plgp (>=1.1-11)`
+
+`MATLAB` needs to be installed in order to run the finite element
+simulations. The MATLAB code can be accessed from the folder
+`Rmatlab_files`.

blade_example.R

@@ -38,9 +38,11 @@ log.fg <- TRUE
 
 set.seed(1)
 X.test <- maximinLHS(20,2)
-y.test <- f(X.test,l=5, MM, tt, return.time=TRUE)
-# load("2d_jet_blade_true.Rdata")
-# y.test <- true.save
+#y.test <- f(X.test,l=5, MM, tt, return.time=TRUE) 
+# running this high fidelity can take a long time
+# therefore, the outputs of y.test are stored as "2d_jet_blade_true.Rdata"
+load("2d_jet_blade_true.Rdata")
+y.test <- true.save
 
 ML.out <- stacking_design(cost, d, n.init, epsilon=5, alpha, 
                           tt, MM, Lmax, f, k, xnew, n.max, 
@@ -151,40 +153,7 @@ g4 <- g4 + theme(axis.text.y = element_text(size = 12))
 grid.arrange(g1, g2, g3, g4, ncol = 4)
 dev.off()
 
-
-# Figure 12 (no longer needed)
-# pdf("alpha_estimation_blade.pdf", width=6, height=2.5)
-# z.val <- c(ML.out$z[[2]][1:ML.out$n.save[[3]][2]], ML.out$z[[3]][1:ML.out$n.save[[3]][3]])
-# df_alpha <- data.frame(mesh_l=rep(1/tt^(2:3), ML.out$n.save[[3]][2:3]),
-#                        z=abs(z.val))
-# lm.fit <- lm(log(z) ~ log(mesh_l), data=df_alpha)
-# 
-# g1 <- ggplot(df_alpha, aes(x=log(mesh_l), y=log(z), group=log(mesh_l))) + 
-#   geom_boxplot(fill='#A4A4A4', color="black", width=0.25)+
-#   xlab(expression(log(italic(T^-l))))+ylab(bquote("log|" ~ f[italic(l)] ~ "(x) -" ~ f[italic(l)-1] ~ "(x)|"))+#ylab(bquote("log|" ~ z[italic(l)] ~  "|"))+
-#   theme_bw()+ggtitle(expression(italic(L)==3)) + theme(plot.title = element_text(hjust = 0.5))
-# 
-# g1 <- g1 + geom_abline(intercept = lm.fit$coefficients[1], 
-#                        slope = lm.fit$coefficients[2], color = "red", linetype = "dashed") 
-# 
-# z.val <- unlist(ML.out$z[2:4])
-# df_alpha <- data.frame(mesh_l=rep(1/tt^(2:4), ML.out$n[2:4]),
-#                        z=abs(z.val))
-# lm.fit <- lm(log(z) ~ log(mesh_l), data=df_alpha)
-# 
-# g2 <- ggplot(df_alpha, aes(x=log(mesh_l), y=log(z), group=log(mesh_l))) + 
-#   geom_boxplot(fill='#A4A4A4', color="black", width=0.25)+
-#   xlab(expression(log(italic(T^-l))))+ylab(bquote("log|" ~ f[italic(l)] ~ "(x) -" ~ f[italic(l)-1] ~ "(x)|"))+#ylab(bquote("log|" ~ z[italic(l)] ~  "|"))+
-#   theme_bw() + ggtitle(expression(italic(L)==4)) + theme(plot.title = element_text(hjust = 0.5))
-# 
-# g2 <- g2 + geom_abline(intercept = lm.fit$coefficients[1], 
-#                        slope = lm.fit$coefficients[2], color = "red", linetype = "dashed") 
-# 
-# grid.arrange(g1, g2, ncol = 2)
-# dev.off()
-
-
-# Figure 13
+# Figure 10
 pdf("blade_prediction.pdf", width=7.5, height = 3.3)
 par(mfrow=c(1,2))
 par(mar=c(2,1,0,0))

currin_example.R

@@ -170,36 +170,6 @@ grid.arrange(g1, g2, g3, g4, ncol = 4)
 dev.off()
 
 
-# Figure 5 (no longer needed)
-# pdf("alpha_estimation.pdf", width=6, height=2.5)
-# z.val <- c(ML.out$z[[2]][1:ML.out$n.save[[3]][2]], ML.out$z[[3]][1:ML.out$n.save[[3]][3]])
-# df_alpha <- data.frame(mesh_l=rep(1/tt^(2:3), ML.out$n.save[[3]][2:3]),
-#                        z=abs(z.val))
-# lm.fit <- lm(log(z) ~ log(mesh_l), data=df_alpha)
-# 
-# g1 <- ggplot(df_alpha, aes(x=log(mesh_l), y=log(z), group=log(mesh_l))) + 
-#   geom_boxplot(fill='#A4A4A4', color="black", width=0.25)+
-#   xlab(expression(log(italic(T^-l))))+ylab(bquote("log|" ~ f[italic(l)] ~ "(x) -" ~ f[italic(l)-1] ~ "(x)|"))+#ylab(bquote("log|" ~ z[italic(l)] ~  "|"))+
-#   theme_bw()+ggtitle(expression(italic(L)==3)) + theme(plot.title = element_text(hjust = 0.5))
-# 
-# g1 <- g1 + geom_abline(intercept = lm.fit$coefficients[1], 
-#                        slope = lm.fit$coefficients[2], color = "red", linetype = "dashed") 
-# 
-# z.val <- unlist(ML.out$z[2:4])
-# df_alpha <- data.frame(mesh_l=rep(1/tt^(2:4), ML.out$n[2:4]),
-#                        z=abs(z.val))
-# lm.fit <- lm(log(z) ~ log(mesh_l), data=df_alpha)
-# 
-# g2 <- ggplot(df_alpha, aes(x=log(mesh_l), y=log(z), group=log(mesh_l))) + 
-#   geom_boxplot(fill='#A4A4A4', color="black", width=0.25)+
-#   xlab(expression(log(italic(T^-l))))+ylab(bquote("log|" ~ f[italic(l)] ~ "(x) -" ~ f[italic(l)-1] ~ "(x)|"))+#ylab(bquote("log|" ~ z[italic(l)] ~  "|"))+
-#   theme_bw() + ggtitle(expression(italic(L)==4)) + theme(plot.title = element_text(hjust = 0.5))
-# 
-# g2 <- g2 + geom_abline(intercept = lm.fit$coefficients[1], 
-#                        slope = lm.fit$coefficients[2], color = "red", linetype = "dashed") 
-# 
-# grid.arrange(g1, g2, ncol = 2)
-# dev.off()
 
 
 
@@ -261,7 +231,7 @@ df2$meshsize <- factor(MM/tt^as.numeric(df2$meshsize))
 df1$meshsize <- factor(df1$meshsize, levels = c(8, 4, 2, 1, 0.5))
 df2$meshsize <- factor(df2$meshsize, levels = c(8, 4, 2, 1, 0.5))
 
-# Figure 6
+# Figure 5
 my.cols <- brewer.pal(6, "Blues")
 pdf(paste0("illustration epsilon.pdf"), width=11, height=2.8)
 g1 <- ggplot(data=df1, aes(x=factor(epsilon))) +