## ***************************************************************************** ## PART 3_BaseFunGpFun ## ***************************************************************************** ## ----fgpm----------------------------------------------------------------------------- library("funGp") set.seed(100) n.tr <- 25 sIn <- expand.grid(x1 = seq(0, 1, length = sqrt(n.tr)), x2 = seq(0, 1, length = sqrt(n.tr))) fIn <- list(f1 = matrix(runif(n.tr * 10), ncol = 10), f2 = matrix(runif(n.tr * 22), ncol = 22)) sOut <- fgp_BB3(sIn, fIn, n.tr) m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut) ## ----retrieveCovPar, results = "markup"----------------------------------------------- c(m1@kern@varHyp, m1@kern@s_lsHyps, m1@kern@f_lsHyps) ## ----plotLOO, fig.show = "hide", fig.height = 5, crop = TRUE-------------------------- plot(m1) ## ----show, results = "markup"--------------------------------------------------------- summary(m1) ## ----pred1---------------------------------------------------------------------------- n.pr <- 100 sIn.pr <- as.matrix(expand.grid(x1 = seq(0, 1, length = sqrt(n.pr)), x2 = seq(0, 1, length = sqrt(n.pr)))) fIn.pr <- list(f1 = matrix(runif(n.pr * 10), ncol = 10), f2 = matrix(runif(n.pr * 22), ncol = 22)) ## ----pred2---------------------------------------------------------------------------- m1.preds <- predict(m1, sIn.pr = sIn.pr, fIn.pr = fIn.pr) t(data.frame(row.names = names(m1.preds), Length = lengths(m1.preds))) ## ----plotPred, fig.show = "hide", fig.height = 4, crop = TRUE------------------------- plot(m1.preds) ## ----pbOptim, crop = TRUE, fig.width = 5, fig.height = 3.5, include = FALSE----------- m1Hack <- m1 m1Hack@kern@f_lsHyps <- c(10^(-2), m1@kern@f_lsHyps[2]) predHack <- predict(m1Hack, sIn.pr = sIn.pr, fIn.pr = fIn.pr) plot(predHack, sortp.gpars = list(legends = FALSE, main = "", xlab = "")) m1Hack <- m1 m1Hack@kern@f_lsHyps <- c(10^(-1), m1@kern@f_lsHyps[2]) predHack <- predict(m1Hack, sIn.pr = sIn.pr, fIn.pr = fIn.pr) plot(predHack, sortp.gpars = list(legends = FALSE, main = "", xlab = "")) ## ----predOut, fig.show = "hide", fig.height = 6, crop = TRUE-------------------------- sOut.pr <- fgp_BB3(sIn.pr, fIn.pr, n.pr) plot(m1.preds, sOut.pr = sOut.pr) ## ----predOut2------------------------------------------------------------------------- m1.preds <- predict(m1, sIn.pr = sIn.pr, fIn.pr = fIn.pr, detail = "full") t(data.frame(row.names = names(m1.preds), Length = lengths(m1.preds))) ## ----simLight------------------------------------------------------------------------- n.sm <- 100 sIn.sm <- as.matrix(expand.grid(x1 = seq(0, 1, length = sqrt(n.sm)), x2 = seq(0, 1, length = sqrt(n.sm)))) fIn.sm <- list(f1 = matrix(runif(n.sm * 10), ncol = 10), f2 = matrix(runif(n.sm * 22), ncol = 22)) m1.sims_l <- simulate(m1, nsim = 10, sIn.sm = sIn.sm, fIn.sm = fIn.sm) ## ----plotSimLight, fig.show = "hide", fig.height = 4, crop = TRUE--------------------- plot(m1.sims_l) ## ----simDetails----------------------------------------------------------------------- m1.sims_f <- simulate(m1, nsim = 10, sIn.sm = sIn.sm, fIn.sm = fIn.sm, detail = "full") t(data.frame(row.names = names(m1.sims_f), Length = lengths(m1.sims_f))) ## ----plotSimFull, fig.show = "hide", fig.height = 4, crop = TRUE---------------------- plot(m1.sims_f) ## ----updateDelete--------------------------------------------------------------------- ind.dl <- sample(1:m1@n.tot, 2) m1up <- update(m1, ind.dl = ind.dl) ## ----updateAdd------------------------------------------------------------------------ n.nw <- 5 sIn.nw <- matrix(runif(n.nw * m1@ds), nrow = n.nw) fIn.nw <- list(f1 = matrix(runif(n.nw * 10), ncol = 10), f2 = matrix(runif(n.nw * 22), ncol = 22)) sOut.nw <- fgp_BB3(sIn.nw, fIn.nw, n.nw) m1up <- update(m1, sIn.nw = sIn.nw, fIn.nw = fIn.nw, sOut.nw = sOut.nw, trace = FALSE) ## ----updateSubst---------------------------------------------------------------------- n.sb <- 2 sIn.sb <- matrix(runif(n.sb * m1@ds), nrow = n.sb) fIn.sb <- list(f1 = matrix(runif(n.sb * 10), ncol = 10), f2 = matrix(runif(n.sb * 22), ncol = 22)) sOut.sb <- fgp_BB3(sIn.sb, fIn.sb, n.sb) ind.sb <- sample(1:(m1@n.tot), n.sb) m1up <- update(m1, sIn.sb = sIn.sb, fIn.sb = fIn.sb, sOut.sb = sOut.sb, ind.sb = ind.sb, trace = FALSE) ## ----updateSubstInStruc--------------------------------------------------------------- m1up1 <- update(m1, sIn.sb = sIn.sb, ind.sb = ind.sb, trace = FALSE) m1up2 <- update(m1, sOut.sb = sOut.sb, ind.sb = ind.sb, trace = FALSE) ## ----updateSubstParam----------------------------------------------------------------- var.sb <- 3 ls_s.sb <- c(2.44, 1.15) ls_f.sb <- c(5.83, 4.12) m1up <- update(m1, var.sb = var.sb, ls_s.sb = ls_s.sb, ls_f.sb = ls_f.sb, trace = FALSE) ## ----updateSubstParamStruct----------------------------------------------------------- m1up <- update(m1, var.sb = var.sb, trace = FALSE) m1up <- update(m1, ls_f.sb = ls_f.sb, trace = FALSE) ## ----updateReestim-------------------------------------------------------------------- m1up <- update(m1, var.re = TRUE, trace = FALSE) m1up <- update(m1, ls_s.re = TRUE, trace = FALSE, control.optim = list(trace = FALSE)) m1up <- update(m1, ls_s.re = TRUE, ls_f.re = TRUE, trace = FALSE, control.optim = list(trace = FALSE)) m1up <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE, trace = FALSE, control.optim = list(trace = FALSE)) ## ----updateReestimExtend-------------------------------------------------------------- m1up <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE, trace = FALSE, control.optim = list(trace = FALSE), extend = TRUE) ## ***************************************************************************** ## PART 4_Customize ## ***************************************************************************** ## ----kernel--------------------------------------------------------------------------- n.tr <- 25 sIn <- expand.grid(x1 = seq(0, 1, length = sqrt(n.tr)), x2 = seq(0, 1, length = sqrt(n.tr))) fIn <- list(f1 = matrix(runif(n.tr * 10), ncol = 10), f2 = matrix(runif(n.tr * 22), ncol = 22)) sOut <- fgp_BB3(sIn, fIn, n.tr) m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, kerType = "gauss", trace = FALSE, control.optim = list(trace = FALSE)) ## ----proj----------------------------------------------------------------------------- m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, f_basType = c("B-splines", "PCA"), trace = FALSE, control.optim = list(trace = FALSE)) ## ----projDim-------------------------------------------------------------------------- m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, f_pdims = c(0, 7), trace = FALSE, control.optim = list(trace = FALSE)) ## ----dist----------------------------------------------------------------------------- n.tr <- 25 sIn <- expand.grid(x1 = seq(0, 1, length = sqrt(n.tr)), x2 = seq(0, 1, length = sqrt(n.tr))) fIn <- list(f1 = matrix(runif(n.tr * 10), ncol = 10), f2 = matrix(runif(n.tr * 22), ncol = 22)) sOut <- fgp_BB3(sIn, fIn, n.tr) ## ----dist1---------------------------------------------------------------------------- m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, f_pdims = c(0, 5), f_disType = c("L2_byindex", "L2_bygroup"), trace = FALSE, control.optim = list(trace = FALSE)) ## ----dist2---------------------------------------------------------------------------- m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, f_pdims = c(0, 5), f_disType = "L2_bygroup", trace = FALSE, control.optim = list(trace = FALSE)) ## ----dist3---------------------------------------------------------------------------- m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, f_pdims = c(0, 5), f_disType = "L2_byindex", trace = FALSE, control.optim = list(trace = FALSE)) ## ***************************************************************************** ## PART 5_Automatic ## ***************************************************************************** ## ----BB7------------------------------------------------------------------------------ set.seed(100) n.tr <- 32 sIn <- expand.grid(x1 = seq(0, 1, length = n.tr^(1/5)), x2 = seq(0, 1, length = n.tr^(1/5)), x3 = seq(0, 1, length = n.tr^(1/5)), x4 = seq(0, 1, length = n.tr^(1/5)), x5 = seq(0, 1, length = n.tr^(1/5))) fIn <- list(f1 = matrix(runif(n.tr * 10), ncol = 10), f2 = matrix(runif(n.tr * 22), ncol = 22)) sOut <- fgp_BB7(sIn, fIn, n.tr) ## ----Xfgpm, results = "hide"---------------------------------------------------------- xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut) ## ----XfgpmModelPlot, fig.show = "hide", fig.width = 4.5, fig.height = 3.8, crop=TRUE---- plot(xm@model, main = "") m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, trace = FALSE, control.optim = list(trace = FALSE)) plot(m1, main = "") ## ----plotDiag, fig.show = "hide", fig.height = 6, fig.width = 8----------------------- plot(xm, which = "diag") ## ----plotXEvol, fig.show = "hide", fig.height = 5, fig.width = 8---------------------- plot(xm, which = "evol") ## ----XfgmpNIter, results = "hide"----------------------------------------------------- set.seed(100) xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, setup = list(n.iter = 30), trace = FALSE) ## ----XfgpmNIterPlotX, echo = FALSE, fig.show = "hide", fig.height = 6, fig.width = 8---- plot(xm, which = "diag") ## ----XfgpmVal15, results = "hide"----------------------------------------------------- ind.vl <- sample(seq_len(n.tr), 5) xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, ind.vl = ind.vl, trace = FALSE) ## ----XfgpmVal15Repl, results = "hide", eval = FALSE----------------------------------- ## ind.vl <- replicate(30, sample(seq_len(n.tr), 5)) ## xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, ind.vl = ind.vl) ## ----CustomHeur, results = "hide"----------------------------------------------------- mysup <- list(tao0 = .15, dop.s = 1.2, dop.f = 1.3, delta.f = 4, dispr.f = 1.1, rho.l = .2, u.gbest = TRUE, n.ibest = 2, rho.g = .08, n.iter = 30, n.pop = 12, q0 = .85) xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, setup = mysup, trace = FALSE) ## ----customConstr, results = "hide", eval = FALSE------------------------------------- ## myctr <- list(s_keepOn = c(1, 2), f_keepOn = c(2), ## f_disTypes = list("2" = c("L2_byindex")), ## f_fixDims = matrix(c(2, 4), ncol = 1), ## f_maxDims = matrix(c(1,5), ncol = 1), ## f_basTypes = list("1" = c("B-splines")), ## kerTypes = c("matern5_2", "gauss")) ## xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, ctraints = myctr) ## ----XfgpmSummary, results = "hide"--------------------------------------------------- summary(xm) ## ----timeStopping, results = "hide"--------------------------------------------------- mysup <- list(n.iter = 2000) xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, setup = mysup, time.lim = 60) ## ----XfgpmPoints, results = "hide"---------------------------------------------------- set.seed(100) n.tr <- 32 sIn <- expand.grid(x1 = seq(0, 1, length = n.tr^(1/5)), x2 = seq(0, 1, length = n.tr^(1/5)), x3 = seq(0, 1, length = n.tr^(1/5)), x4 = seq(0, 1, length = n.tr^(1/5)), x5 = seq(0, 1, length = n.tr^(1/5))) fIn <- list(f1 = matrix(runif(n.tr * 10), ncol = 10), f2 = matrix(runif(n.tr * 22), ncol = 22)) sOut <- fgp_BB7(sIn, fIn, n.tr) xm <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, trace = FALSE) ## ----XfgpmPoint2, results = "hide"---------------------------------------------------- n.tr <- 243 sIn <- expand.grid(x1 = seq(0, 1, length = n.tr^(1/5)), x2 = seq(0, 1, length = n.tr^(1/5)), x3 = seq(0, 1, length = n.tr^(1/5)), x4 = seq(0, 1, length = n.tr^(1/5)), x5 = seq(0, 1, length = n.tr^(1/5))) fIn <- list(f1 = matrix(runif(n.tr * 10), ncol = 10), f2 = matrix(runif(n.tr * 22), ncol = 22)) sOut <- fgp_BB7(sIn, fIn, n.tr) m1 <- eval(modelDef(xm, ind = 1, trace = FALSE, control.optim = list(trace = FALSE))) m2 <- eval(modelDef(xm, ind = 2, trace = FALSE, control.optim = list(trace = FALSE))) m3 <- eval(modelDef(xm, ind = 3, trace = FALSE, control.optim = list(trace = FALSE))) ## ----modStack------------------------------------------------------------------------- modStack <- list(m1, m2, m3) argStack <- xm@log.success@args[1:3] ## ----dataPred------------------------------------------------------------------------- n.pr <- 32 sIn.pr <- expand.grid(x1 = seq(0, 1, length = n.pr^(1/5)), x2 = seq(0, 1, length = n.pr^(1/5)), x3 = seq(0, 1, length = n.pr^(1/5)), x4 = seq(0, 1, length = n.pr^(1/5)), x5 = seq(0, 1, length = n.pr^(1/5))) fIn.pr <- list(f1 = matrix(runif(n.pr * 10), ncol = 10), f2 = matrix(runif(n.pr * 22), ncol = 22)) ## ----doPred--------------------------------------------------------------------------- preds <- do.call(cbind, Map(function(model, args) { active <- get_active_in(sIn = sIn.pr, fIn = fIn.pr, args) predict(model, sIn.pr = active$sIn.on, fIn.pr = active$fIn.on)$mean }, modStack, argStack)) ## ----plotPred3, fig.witdth = 10, fig.height = 5, crop = TRUE, include = FALSE--------- preds <- do.call(cbind, Map(function(model, args) { active <- get_active_in(sIn = sIn.pr, fIn = fIn.pr, args) predict(model, sIn.pr = active$sIn.on, fIn.pr = active$fIn.on)$mean }, modStack, argStack)) plot(1, xlim = c(1, nrow(preds)), ylim = range(preds), xaxt = "n", xlab = "Prediction point index", ylab = "Output", main = "Predictions with best 3 structural configurations") axis(1, 1:nrow(preds)) for (i in seq_len(n.pr)) { lines(rep(i, 2), range(preds[i, 1:3]), col = "grey35", lty = 3) } points(preds[,1], pch = 21, bg = "black") points(preds[,2], pch = 23, bg = "red") points(preds[,3], pch = 24, bg = "green") legend("bottomleft", legend = c("Model 1", "Model 2", "Model 3"), pch = c(21, 23, 24), pt.bg = c("black", "red", "green"), inset = c(.02, .08)) ## ***************************************************************************** ## PART 6_Parallel ## ***************************************************************************** ## ----fgpmNonParallel, results="hide"-------------------------------------------------- m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, n.starts = 10, control.optim = list(trace = FALSE)) ## ----fgpmParallel--------------------------------------------------------------------- cl <- parallel::makeCluster(3) m1 <- fgpm(sIn = sIn, fIn = fIn, sOut = sOut, n.starts = 10, par.clust = cl,control.optim = list(trace = FALSE)) parallel::stopCluster(cl) ## ----XfgpmParallel, results = "hide", message = FALSE, warning = FALSE---------------- cl <- parallel::makeCluster(3) xm.par <- fgpm_factory(sIn = sIn, fIn = fIn, sOut = sOut, par.clust = cl) parallel::stopCluster(cl) ## ***************************************************************************** ## PART 7_RealCase ## ***************************************************************************** ## The dataset "FloodingCase.rda", if reused, should be referenced as an ## extract of the CFMDG dataset whose reference is: Idier D., Rohmer J., Pedreros ## R., Le Roy S., Betancourt J. (2023). CFMDG: a Coastal Flood Modelling Dataset in ## Gavres (France) to support risk prevention and metamodels development ## (V1(15_02_2023)) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.7533336 ## ----FloodingY, fig.show = "hide", fig.height = 5, crop = TRUE------------------------ load("FloodingCase.rda") df <- length(Xf) n <- nrow(Xs) hist(Y, xlab = "Flooded Area (m^2)", main = "Model outputs") ## ----Xf, fig.show = "hide", fig.height = 10, crop = TRUE, echo = FALSE---------------- NOM <- c("Tide (m | msl)","Surge (m)","Hsx (m)","Hsy (m)","Tp (s)","U10x (m/s)","U10y (m/s)") par(mfrow = c(4, 2), mar = c(5, 5, 2, 2), cex.lab = 1.5, cex.axis = 1.5) for (i in 1:df) { matplot(seq(-3, 3, length.out = 37), t(Xf[[i]] + Xs[,i+1]), ty = "l", main = "", xlab = "Time (hr)", ylab = NOM[i]) } ## ----FloodingFit---------------------------------------------------------------------- set.seed(100) id.train <- sample(1:n, 100, replace = FALSE) id.test <- (1:n)[-id.train] Xs.train <- Xs[id.train, ] Xs.test <- Xs[id.test, ] Xf.train <- Xf.test <- list() for (i in 1:df){ Xf.train[[i]] <- Xf[[i]][id.train, ] Xf.test[[i]] <- Xf[[i]][id.test, ] } Y.train <- Y[id.train] Y.test <- Y[id.test] m1 <- fgpm(sIn = Xs.train, fIn = Xf.train, sOut = Y.train, trace = FALSE, control.optim = list(trace = FALSE)) ## ----FloodingQ2, fig.show = "hide", fig.height = 5, crop = TRUE----------------------- plot(m1) ## ----FloodingPredict, fig.show = "hide", fig.height = 5, crop = TRUE------------------ m1.preds <- predict(m1, sIn.pr = Xs.test, fIn.pr = Xf.test) plot(m1.preds) ## ----FloodingACOGP, fig.show = "hide", fig.height = 5, crop = TRUE-------------------- xm <- fgpm_factory(sIn = Xs.train, fIn = Xf.train, sOut = Y.train, trace = FALSE) plot(xm, which = "diag") summary(xm, n = 5) ## ----FloodingFitCheckOptimUn---------------------------------------------------------- m2 <- fgpm(sIn = Xs.train, fIn = Xf.train, sOut = Y.train) ## ----FloodingFitCheckOptimDeux-------------------------------------------------------- m2@convergence ## ----FloodingFitCheckOptimTrois------------------------------------------------------- m3 <- fgpm(sIn = Xs.train, fIn = Xf.train, sOut = Y.train, trace = TRUE, control.optim = list(trace = TRUE, maxit = 500)) ## ----FloodingFitCheckOptimQuatre------------------------------------------------------ m3@convergence ## ----FloodingFitCheckOptimCinq-------------------------------------------------------- m4 <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE, extend = FALSE, trace = FALSE) ## ----FloodingFitCheckOptimSix--------------------------------------------------------- m5 <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE, extend = TRUE, trace = FALSE) ## ***************************************************************************** ## PART 9b_ManualPrediction ## ***************************************************************************** ## ----eval = FALSE--------------------------------------------------------------------- ## y <- m1@sOut ## ----eval = FALSE--------------------------------------------------------------------- ## preds <- predict(m1, sIn.pr = sIn.pr, fIn.pr = fIn.pr, detail = "full") ## K.tp <- preds$K.tp ## K.pp <- preds$K.pp ## ----eval = FALSE--------------------------------------------------------------------- ## K.tt <- tcrossprod(m1@preMats$L) ## ----eval = FALSE--------------------------------------------------------------------- ## y.pr <- t(K.tp) %*% solve(K.tt) %*% y ## v.pr <- diag(K.pp - t(K.tp) %*% solve(K.tt) %*% K.tp) ## ll <- y.pr - 1.96 * sqrt(v.pr) ## ul <- y.pr + 1.96 * sqrt(v.pr)