| Title: | Compute and Visualize Profile Extrema Functions |
|---|---|
| Description: | Computes profile extrema functions for arbitrary functions. If the function is expensive-to-evaluate it computes profile extrema by emulating the function with a Gaussian process (using package 'DiceKriging'). In this case uncertainty quantification on the profile extrema can also be computed. The different plotting functions for profile extrema give the user a tool to better locate excursion sets. |
| Authors: | Dario Azzimonti [aut, cre, cph]
|
| Maintainer: | Dario Azzimonti <[email protected]> |
| License: | GPL-3 |
| Version: | 0.2.1 |
| Built: | 2024-11-07 08:22:37 UTC |
| Source: | https://github.com/dazzimonti/profextrema |
Evaluate profile extrema over other variables with approximations at few values
approxMaxMin(f, fprime = NULL, d, opts = NULL)approxMaxMin(f, fprime = NULL, d, opts = NULL)
f |
the function to be evaluated |
fprime |
derivative of the function |
d |
dimension of the input domain |
opts |
a list containing the options for this function and the subfunctions getMax, getMin or getMaxMinMC, see documentation of getMax, getMin, getMaxMinMC for details. The options only for approxMaxMin are
|
a list of two data frames (min, max) of the evaluations of and
for each i at the design Design. By default Design is a 100 equally spaced points for each dimension. It can be changed by defining it in options$Design
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute the coordinate profile extrema with full optimization on 2d example # Define the function g=function(x){ return(-branin(x)) } # Define the gradient gprime = function(x){ x1 = x[1]*15-5 x2 = x[2]*15 f1prime = (15*25)/(4*pi^4)*x1^3 - (15*75)/(2*pi^3)*x1^2 + (80*15)/(pi^2)*x1 - (5*15)/(pi^2)*x2*x1 + 10*15/pi*x2 - 60*15/pi-10*15* (1 - 1/(8*pi))*sin(x1) f2prime = 2*15*(x2-5/(4*pi^2)*x1^2 +5/pi*x1-6) return(matrix(c(-f1prime,-f2prime),nrow=1)) } # generic approximation options init_des<-lhs::maximinLHS(15,2) options_approx<- list(multistart=4,heavyReturn=TRUE,initDesign=init_des,fullDesignSize=100) # 1order approximation options_approx$smoother<-"1order" coordProf_approx_1order<-approxMaxMin(f = g,fprime = gprime,d=2,opts = options_approx) # quantile regression options_approx$smoother<-"quantSpline" coordProf_approx_quantReg<-approxMaxMin(f = g,fprime = gprime,d=2,opts = options_approx) # Consider threshold=-10 threshold<- -10 # obtain the points where the profiles take the threshold value pp_change<-getChangePoints(threshold = threshold,allRes = coordProf_approx_quantReg) # evaluate g at a grid and plot the image x<-seq(0,1,,100) grid<-expand.grid(x,x) g_evals<- apply(X = grid,MARGIN = 1,FUN = g) image(x = x,y = x,z = matrix(g_evals,nrow = 100),col = grey.colors(20)) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, nlevels = 20) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, levels = threshold,col=2) abline(h = pp_change$neverEx$`-10`[[2]],col="darkgreen",lwd=2) abline(v = pp_change$neverEx$`-10`[[1]],col="darkgreen",lwd=2) # Plot the coordinate profiles and a threshold plotMaxMin(allRes = coordProf_approx_1order,threshold = threshold,changes = TRUE) plotMaxMin(allRes = coordProf_approx_quantReg,threshold = threshold,changes = TRUE)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute the coordinate profile extrema with full optimization on 2d example # Define the function g=function(x){ return(-branin(x)) } # Define the gradient gprime = function(x){ x1 = x[1]*15-5 x2 = x[2]*15 f1prime = (15*25)/(4*pi^4)*x1^3 - (15*75)/(2*pi^3)*x1^2 + (80*15)/(pi^2)*x1 - (5*15)/(pi^2)*x2*x1 + 10*15/pi*x2 - 60*15/pi-10*15* (1 - 1/(8*pi))*sin(x1) f2prime = 2*15*(x2-5/(4*pi^2)*x1^2 +5/pi*x1-6) return(matrix(c(-f1prime,-f2prime),nrow=1)) } # generic approximation options init_des<-lhs::maximinLHS(15,2) options_approx<- list(multistart=4,heavyReturn=TRUE,initDesign=init_des,fullDesignSize=100) # 1order approximation options_approx$smoother<-"1order" coordProf_approx_1order<-approxMaxMin(f = g,fprime = gprime,d=2,opts = options_approx) # quantile regression options_approx$smoother<-"quantSpline" coordProf_approx_quantReg<-approxMaxMin(f = g,fprime = gprime,d=2,opts = options_approx) # Consider threshold=-10 threshold<- -10 # obtain the points where the profiles take the threshold value pp_change<-getChangePoints(threshold = threshold,allRes = coordProf_approx_quantReg) # evaluate g at a grid and plot the image x<-seq(0,1,,100) grid<-expand.grid(x,x) g_evals<- apply(X = grid,MARGIN = 1,FUN = g) image(x = x,y = x,z = matrix(g_evals,nrow = 100),col = grey.colors(20)) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, nlevels = 20) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, levels = threshold,col=2) abline(h = pp_change$neverEx$`-10`[[2]],col="darkgreen",lwd=2) abline(v = pp_change$neverEx$`-10`[[1]],col="darkgreen",lwd=2) # Plot the coordinate profiles and a threshold plotMaxMin(allRes = coordProf_approx_1order,threshold = threshold,changes = TRUE) plotMaxMin(allRes = coordProf_approx_quantReg,threshold = threshold,changes = TRUE)
Evaluate profile extrema for a set of Psi with approximations at few values
approxProfileExtrema(f, fprime = NULL, d, allPsi, opts = NULL)approxProfileExtrema(f, fprime = NULL, d, allPsi, opts = NULL)
f |
the function to be evaluated |
fprime |
derivative of the function |
d |
dimension of the input domain |
allPsi |
a list containing the matrices Psi (dim |
opts |
a list containing the options for this function and the subfunctions getProfileSup_optim, getProfileInf_optim or getProfileExtrema. The options only for approxProfileExtrema are
|
a list of two data frames (min, max) of the evaluations of and
for each i at the design Design. By default Design is a 100 equally spaced points for each dimension. It can be changed by defining it in options$Design
Dario Azzimonti
# Compute the oblique profile extrema with approximate optimization on 2d example # Define the function testF <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(sin(crossprod(v1,x)*params[1]+params[2])+cos(crossprod(v2,x)*params[3]+params[4])-1.5) } testFprime <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(matrix(c(params[1]*v1[1]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[1]*sin(crossprod(v2,x)*params[3]+params[4]), params[1]*v1[2]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[2]*sin(crossprod(v2,x)*params[3]+params[4])),ncol=1)) } # Define the main directions of the function theta=pi/6 pparams<-c(1,0,10,0) vv1<-c(cos(theta),sin(theta)) vv2<-c(cos(theta+pi/2),sin(theta+pi/2)) # Define optimizer friendly function f <-function(x){ return(testF(x,pparams,vv1,vv2)) } fprime <- function(x){ return(testFprime(x,pparams,vv1,vv2)) } # Define list of directions where to evaluate the profile extrema all_Psi <- list(Psi1=vv1,Psi2=vv2) # Evaluate profile extrema along directions of all_Psi allOblique<-approxProfileExtrema(f=f,fprime = fprime,d = 2,allPsi = all_Psi, opts = list(plts=FALSE,heavyReturn=TRUE)) # Consider threshold=0 threshold <- 0 # Plot oblique profile extrema functions plotMaxMin(allOblique,allOblique$Design,threshold = threshold) ## Since the example is two dimensional we can visualize the regions excluded by the profile extrema # evaluate the function at a grid for plots inDes<-seq(0,1,,100) inputs<-expand.grid(inDes,inDes) outs<-apply(X = inputs,MARGIN = 1,function(x){return(testF(x,pparams,v1=vv1,v2=vv2))}) # obtain the points where the profiles take the threshold value cccObl<-getChangePoints(threshold = threshold,allRes = allOblique,Design = allOblique$Design) # visualize the functions and the regions excluded image(inDes,inDes,matrix(outs,ncol=100),col=grey.colors(20),main="Example and oblique profiles") contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,nlevels = 20) contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,levels = c(threshold),col=4,lwd=1.5) plotOblique(cccObl$alwaysEx$`0`[[1]],all_Psi[[1]],col=3) plotOblique(cccObl$alwaysEx$`0`[[2]],all_Psi[[2]],col=3) plotOblique(cccObl$neverEx$`0`[[1]],all_Psi[[1]],col=2) plotOblique(cccObl$neverEx$`0`[[2]],all_Psi[[2]],col=2)# Compute the oblique profile extrema with approximate optimization on 2d example # Define the function testF <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(sin(crossprod(v1,x)*params[1]+params[2])+cos(crossprod(v2,x)*params[3]+params[4])-1.5) } testFprime <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(matrix(c(params[1]*v1[1]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[1]*sin(crossprod(v2,x)*params[3]+params[4]), params[1]*v1[2]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[2]*sin(crossprod(v2,x)*params[3]+params[4])),ncol=1)) } # Define the main directions of the function theta=pi/6 pparams<-c(1,0,10,0) vv1<-c(cos(theta),sin(theta)) vv2<-c(cos(theta+pi/2),sin(theta+pi/2)) # Define optimizer friendly function f <-function(x){ return(testF(x,pparams,vv1,vv2)) } fprime <- function(x){ return(testFprime(x,pparams,vv1,vv2)) } # Define list of directions where to evaluate the profile extrema all_Psi <- list(Psi1=vv1,Psi2=vv2) # Evaluate profile extrema along directions of all_Psi allOblique<-approxProfileExtrema(f=f,fprime = fprime,d = 2,allPsi = all_Psi, opts = list(plts=FALSE,heavyReturn=TRUE)) # Consider threshold=0 threshold <- 0 # Plot oblique profile extrema functions plotMaxMin(allOblique,allOblique$Design,threshold = threshold) ## Since the example is two dimensional we can visualize the regions excluded by the profile extrema # evaluate the function at a grid for plots inDes<-seq(0,1,,100) inputs<-expand.grid(inDes,inDes) outs<-apply(X = inputs,MARGIN = 1,function(x){return(testF(x,pparams,v1=vv1,v2=vv2))}) # obtain the points where the profiles take the threshold value cccObl<-getChangePoints(threshold = threshold,allRes = allOblique,Design = allOblique$Design) # visualize the functions and the regions excluded image(inDes,inDes,matrix(outs,ncol=100),col=grey.colors(20),main="Example and oblique profiles") contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,nlevels = 20) contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,levels = c(threshold),col=4,lwd=1.5) plotOblique(cccObl$alwaysEx$`0`[[1]],all_Psi[[1]],col=3) plotOblique(cccObl$alwaysEx$`0`[[2]],all_Psi[[2]],col=3) plotOblique(cccObl$neverEx$`0`[[1]],all_Psi[[1]],col=2) plotOblique(cccObl$neverEx$`0`[[2]],all_Psi[[2]],col=2)
The function bound_profiles computes the upper and lower bounds for the profile extrema quantiles of a Gaussian process model.
bound_profiles( objectUQ, mean_var_delta = NULL, beta = 0.0124, alpha = 0.025, allPsi = NULL, options_approx = NULL, options_full_sims = NULL )bound_profiles( objectUQ, mean_var_delta = NULL, beta = 0.0124, alpha = 0.025, allPsi = NULL, options_approx = NULL, options_full_sims = NULL )
objectUQ |
an object returned by coordProf_UQ or the object saved in |
mean_var_delta |
the profile extrema functions at |
beta |
the level of confidence for the approximate simulations |
alpha |
the level of confidence for the bound |
allPsi |
optional list of matrices (dim |
options_approx |
an optional list of options for approxMaxMin (or approxProfileExtrema if |
options_full_sims |
an optional list of options for getAllMaxMin (or getProfileExtrema if |
a list containing
bound: a list containing the upper/lower bound for profile sup and inf
approx: a list containing the upper/lower approximate quantiles for profile sup and inf
Dario Azzimonti
The function cleanProfileResults cleans a profile extrema object to partially redo some computations.
cleanProfileResults(object, level = 1)cleanProfileResults(object, level = 1)
object |
a list containing profile extrema results. |
level |
an integer 1-4 denoting how much it should be removed from |
returns object with the deleted parts as selected by level. In particular
1: keep only profMean_full.
2: keep profMean_full and profMean_approx. Remove all UQ results.
3: keep profMean_full and profMean_approx and the pilot points. Remove all UQ simulations.
4: Remove only the bound computations.
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g=function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 # Compute coordinate profiles on the posterior mean # Increase multistart and size of designs for more precise results options_full<-list(multistart=2,heavyReturn=TRUE, Design = replicate(2,seq(0,1,,50))) init_des<-lhs::maximinLHS(12,2) options_approx<- list(multistart=2,heavyReturn=TRUE,initDesign=init_des,fullDesignSize=50) cProfilesMean<-coordinateProfiles(object=kmModel,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=FALSE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) # If we want to run again the computation of approximate coordinate profiles # we delete that result and run again the coordinate profiles function cProfiles_full <- cleanProfileResults(cProfilesMean,level=1) ## Not run: # Coordinate profiles with UQ with approximate profiles plot_options<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) cProfilesUQ<-coordinateProfiles(object=cProfilesMean,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=TRUE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) # If we would like to remove all UQ results cProf_noUQ <- cleanProfileResults(cProfilesUQ,level=2) # If we would like to remove the simulations but keep the pilot points cProf_noSims <- cleanProfileResults(cProfilesUQ,level=3) # the line above is useful, for example, if we need a more accurate UQ. In that case # we obtain more simulations with the same pilot points and then combine the results. ## End(Not run)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g=function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 # Compute coordinate profiles on the posterior mean # Increase multistart and size of designs for more precise results options_full<-list(multistart=2,heavyReturn=TRUE, Design = replicate(2,seq(0,1,,50))) init_des<-lhs::maximinLHS(12,2) options_approx<- list(multistart=2,heavyReturn=TRUE,initDesign=init_des,fullDesignSize=50) cProfilesMean<-coordinateProfiles(object=kmModel,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=FALSE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) # If we want to run again the computation of approximate coordinate profiles # we delete that result and run again the coordinate profiles function cProfiles_full <- cleanProfileResults(cProfilesMean,level=1) ## Not run: # Coordinate profiles with UQ with approximate profiles plot_options<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) cProfilesUQ<-coordinateProfiles(object=cProfilesMean,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=TRUE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) # If we would like to remove all UQ results cProf_noUQ <- cleanProfileResults(cProfilesUQ,level=2) # If we would like to remove the simulations but keep the pilot points cProf_noSims <- cleanProfileResults(cProfilesUQ,level=3) # the line above is useful, for example, if we need a more accurate UQ. In that case # we obtain more simulations with the same pilot points and then combine the results. ## End(Not run)
A dataset containing the results of a numerical simulation conducted with the MARS model (Lazure and Dumas, 2008) for coastal flooding. The numerical model was adapted to the Boucholeurs area (French Atlantic coast), close to La Rochelle, and validated with data from the 2010 Xynthia storm event. See Azzimonti et al. (2017+) and Rohmer et al. (2018) for more details.
coastal_floodingcoastal_flooding
A data frame with 200 rows and 6 variables:
High tide level in meters;
Surge peak amplitude in meters;
Phase difference between high tide and surge peak;
Duration of the increasing part of the surge temporal signal (assumed to be triangular);
Duration of the decreasing part of the surge temporal signal (assumed to be triangular);
Flooded area in m^2.
The data frame contains 5 input variables: Tide, Surge, phi, t-, t+ detailing
the offshore forcing conditions for the model. All input variables are normalized
in . The response is Area, the area flooded in m^2.
Azzimonti, D., Ginsbourger, D., Rohmer, J. and Idier, D. (2017+). Profile extrema for visualizing and quantifying uncertainties on excursion regions. Application to coastal flooding. arXiv:1710.00688.
Rohmer, J., Idier, D., Paris, F., Pedreros, R., and Louisor, J. (2018). Casting light on forcing and breaching scenarios that lead to marine inundation: Combining numerical simulations with a random-forest classification approach. Environmental Modelling & Software, 104:64-80.
Lazure, P. and Dumas, F. (2008). An external-internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS). Advances in Water Resources, 31:233-250.
# Define inputs inputs<-data.frame(coastal_flooding[,-6]) colnames(inputs)<-colnames(coastal_flooding[,-6]) colnames(inputs)[4:5]<-c("tPlus","tMinus") # put response in areaFlooded variable areaFlooded<-data.frame(coastal_flooding[,6]) colnames(areaFlooded)<-colnames(coastal_flooding)[6] response = sqrt(areaFlooded) model <- km(formula=~Tide+Surge+I(phi^2)+tMinus+tPlus, design = inputs,response = response,covtype="matern3_2") # Fix threshold threshold<-sqrt(c(1.2e6,1.9e6,3.1e6,6.5e6)) # use the coordinateProfile function ## set up plot options options_plots <- list(save=FALSE, folderPlots = "./" , titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) # set up full profiles options options_full<-list(multistart=15,heavyReturn=TRUE) # set up approximation options d <- model@d init_des<-lhs::maximinLHS(5*d , d ) options_approx<- list(multistart=2,heavyReturn=TRUE, initDesign=init_des, fullDesignSize=100, smoother="quantSpline") # run the coordinate profile extrema on the mean CF_CoordProf_mean<- coordinateProfiles(object = model, threshold = threshold, uq_computations = FALSE, options_approx = options_approx, plot_level=3, plot_options= options_plots, return_level=3, options_full=options_full) ## Not run: ## UQ computations might require a long time # set up simulation options ## reduce nsims and batchsize for faster/less accurate UQ nsims=200 opts_sims<-list(algorithm="B", lower=rep(0,d ), upper=rep(1,d ), batchsize=150, optimcontrol=list(method="genoud", pop.size=100,print.level=0), integcontrol = list(distrib="sobol",n.points=2000),nsim=nsims) opts_sims$integration.param <- integration_design(opts_sims$integcontrol, d , opts_sims$lower, opts_sims$upper, model,threshold) opts_sims$integration.param$alpha <- 0.5 # run UQ computations CF_CoordProf_UQ<- coordinateProfiles(object = CF_CoordProf_mean, threshold = threshold, uq_computations = TRUE, options_approx = options_approx, plot_level=3, plot_options= options_plots, return_level=3, options_sims=opts_sims,options_full=options_full, options_bound = list(beta=0.024,alpha=0.05)) ## End(Not run)# Define inputs inputs<-data.frame(coastal_flooding[,-6]) colnames(inputs)<-colnames(coastal_flooding[,-6]) colnames(inputs)[4:5]<-c("tPlus","tMinus") # put response in areaFlooded variable areaFlooded<-data.frame(coastal_flooding[,6]) colnames(areaFlooded)<-colnames(coastal_flooding)[6] response = sqrt(areaFlooded) model <- km(formula=~Tide+Surge+I(phi^2)+tMinus+tPlus, design = inputs,response = response,covtype="matern3_2") # Fix threshold threshold<-sqrt(c(1.2e6,1.9e6,3.1e6,6.5e6)) # use the coordinateProfile function ## set up plot options options_plots <- list(save=FALSE, folderPlots = "./" , titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) # set up full profiles options options_full<-list(multistart=15,heavyReturn=TRUE) # set up approximation options d <- model@d init_des<-lhs::maximinLHS(5*d , d ) options_approx<- list(multistart=2,heavyReturn=TRUE, initDesign=init_des, fullDesignSize=100, smoother="quantSpline") # run the coordinate profile extrema on the mean CF_CoordProf_mean<- coordinateProfiles(object = model, threshold = threshold, uq_computations = FALSE, options_approx = options_approx, plot_level=3, plot_options= options_plots, return_level=3, options_full=options_full) ## Not run: ## UQ computations might require a long time # set up simulation options ## reduce nsims and batchsize for faster/less accurate UQ nsims=200 opts_sims<-list(algorithm="B", lower=rep(0,d ), upper=rep(1,d ), batchsize=150, optimcontrol=list(method="genoud", pop.size=100,print.level=0), integcontrol = list(distrib="sobol",n.points=2000),nsim=nsims) opts_sims$integration.param <- integration_design(opts_sims$integcontrol, d , opts_sims$lower, opts_sims$upper, model,threshold) opts_sims$integration.param$alpha <- 0.5 # run UQ computations CF_CoordProf_UQ<- coordinateProfiles(object = CF_CoordProf_mean, threshold = threshold, uq_computations = TRUE, options_approx = options_approx, plot_level=3, plot_options= options_plots, return_level=3, options_sims=opts_sims,options_full=options_full, options_bound = list(beta=0.024,alpha=0.05)) ## End(Not run)
The function coordinateProfiles computes the profile extrema functions for the posterior mean of a Gaussian process and its confidence bounds
coordinateProfiles( object, threshold, options_full = NULL, options_approx = NULL, uq_computations = FALSE, plot_level = 0, plot_options = NULL, CI_const = NULL, return_level = 1, ... )coordinateProfiles( object, threshold, options_full = NULL, options_approx = NULL, uq_computations = FALSE, plot_level = 0, plot_options = NULL, CI_const = NULL, return_level = 1, ... )
object |
either a km model or a list containing partial results. If |
threshold |
the threshold of interest |
options_full |
an optional list of options for getAllMaxMin, see getAllMaxMin for details. |
options_approx |
an optional list of options for approxMaxMin, see approxMaxMin for details. |
uq_computations |
boolean, if TRUE the uq computations for the profile mean are computed. |
plot_level |
an integer to select the plots to return (0=no plots, 1=basic plots, 2= all plots) |
plot_options |
an optional list of parameters for plots. See setPlotOptions for currently available options. |
CI_const |
an optional vector containing the constants for the CI. If not NULL, then profiles extrema for |
return_level |
an integer to select the amount of details returned |
... |
additional parameters to be passed to coordProf_UQ. |
If return_level=1 a list containing
profMean_full: the results of getAllMaxMin for the posterior mean
profMean_approx: the results of approxMaxMin for the posterior mean
res_UQ: the results of coordProf_UQ for the posterior mean
if return_level=2 the same list as above but also including
abs_err: the vector of maximum absolute approximation errors for the profile inf /sup on posterior mean for the chosen approximation
times: a list containing
full: computational time for the full computation of profile extrema
approx: computational time for the approximate computation of profile extrema
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g=function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 # Compute coordinate profiles on the posterior mean # Increase multistart and size of designs for more precise results options_full<-list(multistart=2,heavyReturn=TRUE, Design = replicate(2,seq(0,1,,50))) init_des<-lhs::maximinLHS(12,2) options_approx<- list(multistart=2,heavyReturn=TRUE,initDesign=init_des,fullDesignSize=50) cProfilesMean<-coordinateProfiles(object=kmModel,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=FALSE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) ## Not run: # Coordinate profiles with UQ with approximate profiles plot_options<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) cProfilesUQ<-coordinateProfiles(object=cProfilesMean,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=TRUE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) # Coordinate profiles with UQ with fully optim profiles options_full_sims<-list(multistart=4,heavyReturn=TRUE, Design = replicate(2,seq(0,1,,60))) cProfilesUQ<-coordinateProfiles(object=cProfilesMean,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=TRUE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2, options_full_sims=options_full_sims) ## End(Not run)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g=function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 # Compute coordinate profiles on the posterior mean # Increase multistart and size of designs for more precise results options_full<-list(multistart=2,heavyReturn=TRUE, Design = replicate(2,seq(0,1,,50))) init_des<-lhs::maximinLHS(12,2) options_approx<- list(multistart=2,heavyReturn=TRUE,initDesign=init_des,fullDesignSize=50) cProfilesMean<-coordinateProfiles(object=kmModel,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=FALSE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) ## Not run: # Coordinate profiles with UQ with approximate profiles plot_options<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) cProfilesUQ<-coordinateProfiles(object=cProfilesMean,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=TRUE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2) # Coordinate profiles with UQ with fully optim profiles options_full_sims<-list(multistart=4,heavyReturn=TRUE, Design = replicate(2,seq(0,1,,60))) cProfilesUQ<-coordinateProfiles(object=cProfilesMean,threshold=threshold,options_full=options_full, options_approx=options_approx,uq_computations=TRUE, plot_level=3,plot_options=NULL,CI_const=NULL,return_level=2, options_full_sims=options_full_sims) ## End(Not run)
The function coordProf_UQ computes the profile extrema functions for posterior realizations of a Gaussian process and its confidence bounds
coordProf_UQ( object, threshold, allResMean = NULL, quantiles_uq = c(0.05, 0.95), options_approx = NULL, options_full_sims = NULL, options_sims = NULL, options_bound = NULL, plot_level = 0, plot_options = NULL, return_level = 1 )coordProf_UQ( object, threshold, allResMean = NULL, quantiles_uq = c(0.05, 0.95), options_approx = NULL, options_full_sims = NULL, options_sims = NULL, options_bound = NULL, plot_level = 0, plot_options = NULL, return_level = 1 )
object |
either a km model or a list containing partial results. If |
threshold |
the threshold of interest |
allResMean |
a list resulting from |
quantiles_uq |
a vector containing the quantiles to be computed |
options_approx |
an optional list of options for approxMaxMin, see approxMaxMin for details. |
options_full_sims |
an optional list of options for getAllMaxMin, see getAllMaxMin for details. If NULL the full computations are not excuted. NOTE: this computations might be very expensive! |
options_sims |
an optional list of options for the posterior simulations.
|
options_bound |
an optional list containing |
plot_level |
an integer to select the plots to return (0=no plots, 1=basic plots, 2= all plots) |
plot_options |
an optional list of parameters for plots. See setPlotOptions for currently available options. |
return_level |
an integer to select the amount of details returned |
If return_level=1 a list containing
profSups: an array dxfullDesignSizexnsims containing the profile sup for each coordinate for each realization.
profInfs: an array dxfullDesignSizexnsims containing the profile inf for each coordinate for each realization.
prof_quantiles_approx: a list containing the quantiles (levels set by quantiles_uq) of the profile extrema functions.
if return_level=2 the same list as above but also including more: a list containing
times: a list containing
tSpts: computational time for selecting pilot points.
tApprox1ord: vector containing the computational time required for profile extrema computation for each realization
simuls: a matrix containing the value of the field simulated at the pilot points
sPts: the pilot points
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g<-function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 d<-2 # Compute coordinate profiles UQ starting from GP model # define simulation options options_sims<-list(algorithm="B", lower=rep(0,d), upper=rep(1,d), batchsize=80, optimcontrol = list(method="genoud",pop.size=100,print.level=0), integcontrol = list(distrib="sobol",n.points=1000), nsim=150) # define 1 order approximation options init_des<-lhs::maximinLHS(15,d) options_approx<- list(multistart=4,heavyReturn=TRUE, initDesign=init_des,fullDesignSize=100, smoother="1order") # define plot options options_plots<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) ## Not run: # profile UQ on approximate coordinate profiles cProfiles_UQ<-coordProf_UQ(object = kmModel,threshold = threshold,allResMean = NULL, quantiles_uq = c(0.05,0.95),options_approx = options_approx, options_full_sims = NULL,options_sims = options_sims, options_bound = NULL,plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim coordinate profiles options_full_sims<-list(multistart=4,heavyReturn=TRUE) cProfiles_UQ_full<-coordProf_UQ(object = cProfiles_UQ,threshold = threshold,allResMean = NULL, quantiles_uq = c(0.05,0.95),options_approx = options_approx, options_full_sims = options_full_sims,options_sims = options_sims, options_bound = NULL,plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim coordinate profiles with bound cProfiles_UQ_full_bound<-coordProf_UQ(object = cProfiles_UQ_full,threshold = threshold, allResMean = NULL, quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = options_full_sims, options_sims = options_sims, options_bound = list(beta=0.024,alpha=0.05), plot_level = 3, plot_options = options_plots, return_level = 3) ## End(Not run)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g<-function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 d<-2 # Compute coordinate profiles UQ starting from GP model # define simulation options options_sims<-list(algorithm="B", lower=rep(0,d), upper=rep(1,d), batchsize=80, optimcontrol = list(method="genoud",pop.size=100,print.level=0), integcontrol = list(distrib="sobol",n.points=1000), nsim=150) # define 1 order approximation options init_des<-lhs::maximinLHS(15,d) options_approx<- list(multistart=4,heavyReturn=TRUE, initDesign=init_des,fullDesignSize=100, smoother="1order") # define plot options options_plots<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) ## Not run: # profile UQ on approximate coordinate profiles cProfiles_UQ<-coordProf_UQ(object = kmModel,threshold = threshold,allResMean = NULL, quantiles_uq = c(0.05,0.95),options_approx = options_approx, options_full_sims = NULL,options_sims = options_sims, options_bound = NULL,plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim coordinate profiles options_full_sims<-list(multistart=4,heavyReturn=TRUE) cProfiles_UQ_full<-coordProf_UQ(object = cProfiles_UQ,threshold = threshold,allResMean = NULL, quantiles_uq = c(0.05,0.95),options_approx = options_approx, options_full_sims = options_full_sims,options_sims = options_sims, options_bound = NULL,plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim coordinate profiles with bound cProfiles_UQ_full_bound<-coordProf_UQ(object = cProfiles_UQ_full,threshold = threshold, allResMean = NULL, quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = options_full_sims, options_sims = options_sims, options_bound = list(beta=0.024,alpha=0.05), plot_level = 3, plot_options = options_plots, return_level = 3) ## End(Not run)
Evaluate coordinate profile extrema with full optimization.
getAllMaxMin(f, fprime = NULL, d, options = NULL)getAllMaxMin(f, fprime = NULL, d, options = NULL)
f |
the function to be evaluated |
fprime |
derivative of the function |
d |
dimension of the input domain |
options |
a list containing the options for this function and the subfunctions getMax, getMin see documentation of getMax, getMin for details. The options only for getAllMaxMin are
|
a list of two data frames (min, max) of the evaluations of and
for each i at the design Design. By default Design is a 100 equally spaced points for each dimension. It can be changed by defining it in options$Design
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute the coordinate profile extrema with full optimization on 2d example # Define the function g=function(x){ return(-branin(x)) } # Define the gradient gprime = function(x){ x1 = x[1]*15-5 x2 = x[2]*15 f1prime = (15*25)/(4*pi^4)*x1^3 - (15*75)/(2*pi^3)*x1^2 + (80*15)/(pi^2)*x1 - (5*15)/(pi^2)*x2*x1 + 10*15/pi*x2 - 60*15/pi-10*15* (1 - 1/(8*pi))*sin(x1) f2prime = 2*15*(x2-5/(4*pi^2)*x1^2 +5/pi*x1-6) return(c(-f1prime,-f2prime)) } # set up dimension coordProf<-getAllMaxMin(f = g,fprime = gprime,d=2,options = list(multistart=4,heavyReturn=TRUE)) # Consider threshold=-10 threshold<- -10 # obtain the points where the profiles take the threshold value pp_change<-getChangePoints(threshold = threshold,allRes = coordProf) # evaluate g at a grid and plot the image x<-seq(0,1,,100) grid<-expand.grid(x,x) g_evals<- apply(X = grid,MARGIN = 1,FUN = g) image(x = x,y = x,z = matrix(g_evals,nrow = 100),col = grey.colors(20)) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, nlevels = 20) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, levels = threshold,col=2) abline(h = pp_change$neverEx$`-10`[[2]],col="darkgreen",lwd=2) abline(v = pp_change$neverEx$`-10`[[1]],col="darkgreen",lwd=2) # Plot the coordinate profiles and a threshold plotMaxMin(allRes = coordProf,threshold = threshold,changes = TRUE)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute the coordinate profile extrema with full optimization on 2d example # Define the function g=function(x){ return(-branin(x)) } # Define the gradient gprime = function(x){ x1 = x[1]*15-5 x2 = x[2]*15 f1prime = (15*25)/(4*pi^4)*x1^3 - (15*75)/(2*pi^3)*x1^2 + (80*15)/(pi^2)*x1 - (5*15)/(pi^2)*x2*x1 + 10*15/pi*x2 - 60*15/pi-10*15* (1 - 1/(8*pi))*sin(x1) f2prime = 2*15*(x2-5/(4*pi^2)*x1^2 +5/pi*x1-6) return(c(-f1prime,-f2prime)) } # set up dimension coordProf<-getAllMaxMin(f = g,fprime = gprime,d=2,options = list(multistart=4,heavyReturn=TRUE)) # Consider threshold=-10 threshold<- -10 # obtain the points where the profiles take the threshold value pp_change<-getChangePoints(threshold = threshold,allRes = coordProf) # evaluate g at a grid and plot the image x<-seq(0,1,,100) grid<-expand.grid(x,x) g_evals<- apply(X = grid,MARGIN = 1,FUN = g) image(x = x,y = x,z = matrix(g_evals,nrow = 100),col = grey.colors(20)) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, nlevels = 20) contour(x=x,y=x,z=matrix(g_evals,nrow = 100), add=TRUE, levels = threshold,col=2) abline(h = pp_change$neverEx$`-10`[[2]],col="darkgreen",lwd=2) abline(v = pp_change$neverEx$`-10`[[1]],col="darkgreen",lwd=2) # Plot the coordinate profiles and a threshold plotMaxMin(allRes = coordProf,threshold = threshold,changes = TRUE)
Obtain the points where the coordinate profile extrema functions cross the threshold
getChangePoints(threshold, Design = NULL, allRes)getChangePoints(threshold, Design = NULL, allRes)
threshold |
if not null plots the level |
Design |
a d dimensional design corresponding to the points |
allRes |
list containing the list |
returns a list containing two lists with d components where
alwaysEx: each component is a numerical vector indicating the points where threshold;
neverEx: each component is a numerical vector indicating the points where threshold.
Dario Azzimonti
Obtain points close in one specific dimension
getClosePoints(x, allPts, whichDim)getClosePoints(x, allPts, whichDim)
x |
one dimensional point |
allPts |
dataframe containing a list of d dimensional points |
whichDim |
integer defining the dimension of x |
the index in allPts (row number) of the closest point in allPts to x along the whichDim dimension
Dario Azzimonti
Compute coordinate profile sup functions
getMax(x, f, fprime, coord, d, options = NULL)getMax(x, f, fprime, coord, d, options = NULL)
x |
one dimensional point where the function is to be evaluated |
f |
function to be optimized (takes a vector y of dimension d and returns a real number) |
fprime |
derivative of f (same format) |
coord |
integer selecting the dimension that is fixed, the other ones are optimized over |
d |
dimension of the input for f |
options |
a list containing the options to be passed to optim:
|
a real value corresponding to
Dario Azzimonti
Compute coordinate profile extrema with MC
getMaxMinMC(x, f, fprime, coord, d, options = NULL)getMaxMinMC(x, f, fprime, coord, d, options = NULL)
x |
one dimensional point where the function is to be evaluated |
f |
function to be optimized (takes a vector y of dimension d and returns a real number) |
fprime |
derivative of f (same format) |
coord |
integer selecting the dimension that is fixed, the other ones are optimized over |
d |
dimension of the input for f |
options |
a list containing the options to be passed to the MC optimizer:
|
a real value corresponding to
Dario Azzimonti
Compute coordinate profile inf functions
getMin(x, f, fprime, coord, d, options = NULL)getMin(x, f, fprime, coord, d, options = NULL)
x |
one dimensional point where the function is to be evaluated |
f |
function to be optimized (takes a vector y of dimension d and returns a real number) |
fprime |
derivative of f (same format) |
coord |
integer selecting the dimension that is fixed, the other ones are optimized over |
d |
dimension of the input for f |
options |
a list containing the options to be passed to optim:
|
a real value corresponding to
Dario Azzimonti
Computes the proportion of observations in the excursion set from true function evaluations,
binned by the grid determined with xBins, yBins.
getPointProportion(pp, xBins, yBins, whichAbove, plt = FALSE)getPointProportion(pp, xBins, yBins, whichAbove, plt = FALSE)
pp |
a matrix of dimension nPts x 2 with the true points locations in 2 dimensions. |
xBins |
numerical vector with the ordered breaks of the grid along the x axis |
yBins |
numerical vector with the ordered breaks of the grid along the y axis |
whichAbove |
boolean vector of dimension nPts, selects the points above |
plt |
if not |
a list containing above, the counts of points in excursion, full the counts per cell of all points,
freq, the relative frequence.
Dario Azzimonti
Evaluate profile extrema for a set of matrices allPsi with full optimization.
getProfileExtrema(f, fprime = NULL, d, allPsi, opts = NULL)getProfileExtrema(f, fprime = NULL, d, allPsi, opts = NULL)
f |
the function to be evaluated |
fprime |
derivative of the function |
d |
dimension of the input domain |
allPsi |
a list containing the matrices Psi (dim |
opts |
a list containing the options for this function and the subfunctions getProfileSup_optim, getProfileInf_optim. The options only for getProfileExtrema are
|
a list of two data frames (min, max) of the evaluations of and
discretized over 50 equally spaced points for each dimension for each Psi in allPsi. This number can be changed by defining it in options$discretization.
Dario Azzimonti
# Compute the oblique profile extrema with full optimization on 2d example # Define the function testF <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(sin(crossprod(v1,x)*params[1]+params[2])+cos(crossprod(v2,x)*params[3]+params[4])-1.5) } testFprime <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(matrix(c(params[1]*v1[1]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[1]*sin(crossprod(v2,x)*params[3]+params[4]), params[1]*v1[2]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[2]*sin(crossprod(v2,x)*params[3]+params[4])),ncol=1)) } # Define the main directions of the function theta=pi/6 pparams<-c(1,0,10,0) vv1<-c(cos(theta),sin(theta)) vv2<-c(cos(theta+pi/2),sin(theta+pi/2)) # Define optimizer friendly function f <-function(x){ return(testF(x,pparams,vv1,vv2)) } fprime <- function(x){ return(testFprime(x,pparams,vv1,vv2)) } # Define list of directions where to evaluate the profile extrema all_Psi <- list(Psi1=vv1,Psi2=vv2) # Evaluate profile extrema along directions of all_Psi allOblique<-getProfileExtrema(f=f,fprime = fprime,d = 2,allPsi = all_Psi, opts = list(plts=FALSE,discretization=100,multistart=8)) # Consider threshold=0 threshold <- 0 # Plot oblique profile extrema functions plotMaxMin(allOblique,allOblique$Design,threshold = threshold) ## Since the example is two dimensional we can visualize the regions excluded by the profile extrema # evaluate the function at a grid for plots inDes<-seq(0,1,,100) inputs<-expand.grid(inDes,inDes) outs<-apply(X = inputs,MARGIN = 1,function(x){return(testF(x,pparams,v1=vv1,v2=vv2))}) # obtain the points where the profiles take the threshold value cccObl<-getChangePoints(threshold = threshold,allRes = allOblique,Design = allOblique$Design) # visualize the functions and the regions excluded image(inDes,inDes,matrix(outs,ncol=100),col=grey.colors(20),main="Example and oblique profiles") contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,nlevels = 20) contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,levels = c(threshold),col=4,lwd=1.5) plotOblique(cccObl$alwaysEx$`0`[[1]],all_Psi[[1]],col=3) plotOblique(cccObl$alwaysEx$`0`[[2]],all_Psi[[2]],col=3) plotOblique(cccObl$neverEx$`0`[[1]],all_Psi[[1]],col=2) plotOblique(cccObl$neverEx$`0`[[2]],all_Psi[[2]],col=2)# Compute the oblique profile extrema with full optimization on 2d example # Define the function testF <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(sin(crossprod(v1,x)*params[1]+params[2])+cos(crossprod(v2,x)*params[3]+params[4])-1.5) } testFprime <- function(x,params,v1=c(1,0),v2=c(0,1)){ return(matrix(c(params[1]*v1[1]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[1]*sin(crossprod(v2,x)*params[3]+params[4]), params[1]*v1[2]*cos(crossprod(v1,x)*params[1]+params[2])- params[3]*v2[2]*sin(crossprod(v2,x)*params[3]+params[4])),ncol=1)) } # Define the main directions of the function theta=pi/6 pparams<-c(1,0,10,0) vv1<-c(cos(theta),sin(theta)) vv2<-c(cos(theta+pi/2),sin(theta+pi/2)) # Define optimizer friendly function f <-function(x){ return(testF(x,pparams,vv1,vv2)) } fprime <- function(x){ return(testFprime(x,pparams,vv1,vv2)) } # Define list of directions where to evaluate the profile extrema all_Psi <- list(Psi1=vv1,Psi2=vv2) # Evaluate profile extrema along directions of all_Psi allOblique<-getProfileExtrema(f=f,fprime = fprime,d = 2,allPsi = all_Psi, opts = list(plts=FALSE,discretization=100,multistart=8)) # Consider threshold=0 threshold <- 0 # Plot oblique profile extrema functions plotMaxMin(allOblique,allOblique$Design,threshold = threshold) ## Since the example is two dimensional we can visualize the regions excluded by the profile extrema # evaluate the function at a grid for plots inDes<-seq(0,1,,100) inputs<-expand.grid(inDes,inDes) outs<-apply(X = inputs,MARGIN = 1,function(x){return(testF(x,pparams,v1=vv1,v2=vv2))}) # obtain the points where the profiles take the threshold value cccObl<-getChangePoints(threshold = threshold,allRes = allOblique,Design = allOblique$Design) # visualize the functions and the regions excluded image(inDes,inDes,matrix(outs,ncol=100),col=grey.colors(20),main="Example and oblique profiles") contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,nlevels = 20) contour(inDes,inDes,matrix(outs,ncol=100),add=TRUE,levels = c(threshold),col=4,lwd=1.5) plotOblique(cccObl$alwaysEx$`0`[[1]],all_Psi[[1]],col=3) plotOblique(cccObl$alwaysEx$`0`[[2]],all_Psi[[2]],col=3) plotOblique(cccObl$neverEx$`0`[[1]],all_Psi[[1]],col=2) plotOblique(cccObl$neverEx$`0`[[2]],all_Psi[[2]],col=2)
Compute profile inf function for an arbitrary matrix Psi with with the L-BFGS-B algorithm of optim. Here the linear equality constraint is eliminated by using the Null space of Psi.
getProfileInf_optim(eta, Psi, f, fprime, d, options = NULL)getProfileInf_optim(eta, Psi, f, fprime, d, options = NULL)
eta |
|
Psi |
projection matrix of dimension |
f |
function to be optimized (takes a vector y of dimension d and returns a real number) |
fprime |
derivative of f (same format, returning a |
d |
dimension of the input for f |
options |
a list containing the options to be passed to optim:
|
a real value corresponding to
Dario Azzimonti
getProfileSup_optim, plotMaxMin
Compute profile sup function for an arbitrary matrix Psi with the L-BFGS-B algorithm of optim.
getProfileSup_optim(eta, Psi, f, fprime, d, options = NULL)getProfileSup_optim(eta, Psi, f, fprime, d, options = NULL)
eta |
|
Psi |
projection matrix of dimensions |
f |
function to be optimized (takes a vector y of dimension d and returns a real number) |
fprime |
derivative of f (same format, returning a |
d |
dimension of the input for f |
options |
a list containing the options to be passed to optim:
|
a real value corresponding to
Dario Azzimonti
getProfileInf_optim, plotMaxMin
Auxiliary function for getChangePoints
getSegments(y)getSegments(y)
y |
a vector |
plots the sup and inf of the function for each dimension. If threshold is not NULL
Dario Azzimonti
The function grad_mean_Delta_T computes the gradient for the mean function of the difference process at x.
grad_mean_Delta_T(x, kmModel, simupoints, T.mat, F.mat)grad_mean_Delta_T(x, kmModel, simupoints, T.mat, F.mat)
x |
a matrix |
kmModel |
the km model of the Gaussian process |
simupoints |
the matrix |
T.mat |
the upper triangular factor of the Choleski decomposition of the covariance matrix of |
F.mat |
the evaluation of the trend function at |
the value of the gradient for the mean function at x for the difference process .
Dario Azzimonti
The function grad_var_Delta_T computes the gradient for the variance function of the difference process at x.
grad_var_Delta_T(x, kmModel, simupoints, T.mat, F.mat)grad_var_Delta_T(x, kmModel, simupoints, T.mat, F.mat)
x |
a matrix |
kmModel |
the km model of the Gaussian process |
simupoints |
the matrix |
T.mat |
the upper triangular factor of the Choleski decomposition of the covariance matrix of |
F.mat |
the evaluation of the trend function at |
the value of the gradient for the variance function at x for the difference process .
Dario Azzimonti
Computes the gradient of the posterior mean and variance of the kriging model in object at the points newdata.
gradKm_dnewdata( object, newdata, type, se.compute = TRUE, light.return = FALSE, bias.correct = FALSE )gradKm_dnewdata( object, newdata, type, se.compute = TRUE, light.return = FALSE, bias.correct = FALSE )
object |
a km object |
newdata |
a vector, matrix or data frame containing the points where to perform predictions. |
type |
a character corresponding to the type of kriging family ( |
se.compute |
an optional boolean indicating whether to compute the posterior variance or not. Default is TRUE. |
light.return |
an optional boolean indicating whether to return additional variables. Default is FALSE. |
bias.correct |
an optional boolean to correct bias in the UK variance. Default is FALSE. |
Returns a list containing
mean: the gradient of the posterior mean at newdata.
trend: the gradient of the trend at newdata.
s2: the gradient of the posterior variance at newdata.
Dario Azzimonti
Compute first order approximation of function from evaluations and gradient
kGradSmooth(newPoints, profPoints, profEvals, profGradient)kGradSmooth(newPoints, profPoints, profEvals, profGradient)
newPoints |
vector of points where to approximate the function |
profPoints |
locations where the function was evaluated |
profEvals |
value of the evaluation at profPoints |
profGradient |
value of the gradient at profPoints |
approximated values of the function at newPoints
Dario Azzimonti
The function mean_Delta_T computes the mean function of the difference process at x.
mean_Delta_T(x, kmModel, simupoints, T.mat, F.mat)mean_Delta_T(x, kmModel, simupoints, T.mat, F.mat)
x |
a matrix |
kmModel |
the km model of the Gaussian process |
simupoints |
the matrix |
T.mat |
the upper triangular factor of the Choleski decomposition of the covariance matrix of |
F.mat |
the evaluation of the trend function at |
the value of the mean function at x for the difference process .
Dario Azzimonti
The function obliqueProf_UQ computes the profile extrema functions for posterior realizations of a Gaussian process and its confidence bounds
obliqueProf_UQ( object, allPsi, threshold, allResMean = NULL, quantiles_uq = c(0.05, 0.95), options_approx = NULL, options_full_sims = NULL, options_sims = NULL, options_bound = NULL, plot_level = 0, plot_options = NULL, return_level = 1 )obliqueProf_UQ( object, allPsi, threshold, allResMean = NULL, quantiles_uq = c(0.05, 0.95), options_approx = NULL, options_full_sims = NULL, options_sims = NULL, options_bound = NULL, plot_level = 0, plot_options = NULL, return_level = 1 )
object |
either a km model or a list containing partial results. If |
allPsi |
a list containing the matrices Psi (dim |
threshold |
the threshold of interest |
allResMean |
a list resulting from |
quantiles_uq |
a vector containing the quantiles to be computed |
options_approx |
an optional list of options for approxProfileExtrema, see approxProfileExtrema for details. |
options_full_sims |
an optional list of options for getProfileExtrema, see getProfileExtrema for details. If NULL the full computations are not executed. NOTE: this computations might be very expensive! |
options_sims |
an optional list of options for the posterior simulations.
|
options_bound |
an optional list containing |
plot_level |
an integer to select the plots to return (0=no plots, 1=basic plots, 2= all plots) |
plot_options |
an optional list of parameters for plots. See setPlotOptions for currently available options. |
return_level |
an integer to select the amount of details returned |
If return_level=1 a list containing
profSups: an array dxfullDesignSizexnsims containing the profile sup for each coordinate for each realization.
profInfs: an array dxfullDesignSizexnsims containing the profile inf for each coordinate for each realization.
prof_quantiles_approx: a list containing the quantiles (levels set by quantiles_uq) of the profile extrema functions.
if return_level=2 the same list as above but also including more: a list containing
times: a list containing
tSpts: computational time for selecting pilot points.
tApprox1ord: vector containing the computational time required for profile extrema computation for each realization
simuls: a matrix containing the value of the field simulated at the pilot points
sPts: the pilot points
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g<-function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 d<-2 # Compute oblique profiles UQ starting from GP model # define simulation options options_sims<-list(algorithm="B", lower=rep(0,d), upper=rep(1,d), batchsize=80, optimcontrol = list(method="genoud",pop.size=100,print.level=0), integcontrol = list(distrib="sobol",n.points=1000), nsim=150) # define approximation options options_approx<- list(multistart=4,heavyReturn=TRUE, initDesign=NULL,fullDesignSize=100, smoother=NULL) # define plot options options_plots<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) # Define the oblique directions # (for theta=0 it is equal to coordinateProfiles) theta=pi/4 allPsi = list(Psi1=matrix(c(cos(theta),sin(theta)),ncol=2), Psi2=matrix(c(cos(theta+pi/2),sin(theta+pi/2)),ncol=2)) ## Not run: # here we reduce the number of simulations to speed up the example # a higher number should be used options_sims$nsim <- 50 # profile UQ on approximate oblique profiles oProfiles_UQ<-obliqueProf_UQ(object = kmModel,threshold = threshold,allPsi=allPsi, allResMean = NULL,quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = NULL, options_sims = options_sims,options_bound = NULL, plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim oblique profiles options_full_sims<-list(multistart=4,heavyReturn=TRUE) oProfiles_UQ_full<- obliqueProf_UQ(object = oProfiles_UQ,threshold = threshold,allPsi=allPsi, allResMean = NULL,quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = options_full_sims, options_sims = options_sims,options_bound = NULL, plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim oblique profiles with bound oProfiles_UQ_full_bound<-obliqueProf_UQ(object = oProfiles_UQ_full,threshold = threshold, allPsi=allPsi, allResMean = NULL, quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = options_full_sims, options_sims = options_sims, options_bound = list(beta=0.024,alpha=0.05), plot_level = 3, plot_options = options_plots, return_level = 3) ## End(Not run)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g<-function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 d<-2 # Compute oblique profiles UQ starting from GP model # define simulation options options_sims<-list(algorithm="B", lower=rep(0,d), upper=rep(1,d), batchsize=80, optimcontrol = list(method="genoud",pop.size=100,print.level=0), integcontrol = list(distrib="sobol",n.points=1000), nsim=150) # define approximation options options_approx<- list(multistart=4,heavyReturn=TRUE, initDesign=NULL,fullDesignSize=100, smoother=NULL) # define plot options options_plots<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) # Define the oblique directions # (for theta=0 it is equal to coordinateProfiles) theta=pi/4 allPsi = list(Psi1=matrix(c(cos(theta),sin(theta)),ncol=2), Psi2=matrix(c(cos(theta+pi/2),sin(theta+pi/2)),ncol=2)) ## Not run: # here we reduce the number of simulations to speed up the example # a higher number should be used options_sims$nsim <- 50 # profile UQ on approximate oblique profiles oProfiles_UQ<-obliqueProf_UQ(object = kmModel,threshold = threshold,allPsi=allPsi, allResMean = NULL,quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = NULL, options_sims = options_sims,options_bound = NULL, plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim oblique profiles options_full_sims<-list(multistart=4,heavyReturn=TRUE) oProfiles_UQ_full<- obliqueProf_UQ(object = oProfiles_UQ,threshold = threshold,allPsi=allPsi, allResMean = NULL,quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = options_full_sims, options_sims = options_sims,options_bound = NULL, plot_level = 3, plot_options = options_plots,return_level = 3) # profile UQ on full optim oblique profiles with bound oProfiles_UQ_full_bound<-obliqueProf_UQ(object = oProfiles_UQ_full,threshold = threshold, allPsi=allPsi, allResMean = NULL, quantiles_uq = c(0.05,0.95), options_approx = options_approx, options_full_sims = options_full_sims, options_sims = options_sims, options_bound = list(beta=0.024,alpha=0.05), plot_level = 3, plot_options = options_plots, return_level = 3) ## End(Not run)
The function obliqueProfiles computes the (oblique) profile extrema functions for the posterior mean of a Gaussian process and its confidence bounds
obliqueProfiles( object, allPsi, threshold, options_full = NULL, options_approx = NULL, uq_computations = FALSE, plot_level = 0, plot_options = NULL, CI_const = NULL, return_level = 1, ... )obliqueProfiles( object, allPsi, threshold, options_full = NULL, options_approx = NULL, uq_computations = FALSE, plot_level = 0, plot_options = NULL, CI_const = NULL, return_level = 1, ... )
object |
either a km model or a list containing partial results. If |
allPsi |
a list containing the matrices Psi (dim |
threshold |
the threshold of interest |
options_full |
an optional list of options for getProfileExtrema, see getProfileExtrema for details. |
options_approx |
an optional list of options for approxProfileExtrema, see approxProfileExtrema for details. |
uq_computations |
boolean, if TRUE the uq computations for the profile mean are computed. |
plot_level |
an integer to select the plots to return (0=no plots, 1=basic plots, 2= all plots) |
plot_options |
an optional list of parameters for plots. See setPlotOptions for currently available options. |
CI_const |
an optional vector containing the constants for the CI. If not NULL, then profiles extrema for |
return_level |
an integer to select the amount of details returned |
... |
additional parameters to be passed to obliqueProf_UQ. |
If return_level=1 a list containing
profMean_full: the results of getProfileExtrema for the posterior mean
profMean_approx: the results of approxProfileExtrema for the posterior mean
res_UQ: the results of obliqueProf_UQ for the posterior mean
if return_level=2 the same list as above but also including
abs_err: the vector of maximum absolute approximation errors for the profile inf /sup on posterior mean for the chosen approximation
times: a list containing
full: computational time for the full computation of profile extrema
approx: computational time for the approximate computation of profile extrema
Dario Azzimonti
if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g=function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 # Compute oblique profiles on the posterior mean # (for theta=0 it is equal to coordinateProfiles) options_full<-list(multistart=4,heavyReturn=TRUE,discretization=100) options_approx<- list(multistart=4,heavyReturn=TRUE,initDesign=NULL,fullDesignSize=100) theta=pi/4 allPsi = list(Psi1=matrix(c(cos(theta),sin(theta)),ncol=2), Psi2=matrix(c(cos(theta+pi/2),sin(theta+pi/2)),ncol=2)) ## Not run: profMeans<-obliqueProfiles(object = kmModel,allPsi = allPsi,threshold = threshold, options_full = options_full,options_approx = options_approx, uq_computations = FALSE,plot_level = 3,plot_options = NULL, CI_const = NULL,return_level = 2) # Approximate oblique profiles with UQ plot_options<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) options_sims<-list(nsim=150) obProfUQ<-obliqueProfiles(object=profMeans,threshold=threshold,allPsi = allPsi, options_full=options_full, options_approx=options_approx, uq_computations=TRUE, plot_level=3,plot_options=NULL, CI_const=NULL,return_level=2,options_sims=options_sims) ## End(Not run)if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } # Compute a kriging model from 50 evaluations of the Branin function # Define the function g=function(x){ return(-branin(x)) } gp_des<-lhs::maximinLHS(20,2) reals<-apply(gp_des,1,g) kmModel<-km(design = gp_des,response = reals,covtype = "matern3_2") threshold=-10 # Compute oblique profiles on the posterior mean # (for theta=0 it is equal to coordinateProfiles) options_full<-list(multistart=4,heavyReturn=TRUE,discretization=100) options_approx<- list(multistart=4,heavyReturn=TRUE,initDesign=NULL,fullDesignSize=100) theta=pi/4 allPsi = list(Psi1=matrix(c(cos(theta),sin(theta)),ncol=2), Psi2=matrix(c(cos(theta+pi/2),sin(theta+pi/2)),ncol=2)) ## Not run: profMeans<-obliqueProfiles(object = kmModel,allPsi = allPsi,threshold = threshold, options_full = options_full,options_approx = options_approx, uq_computations = FALSE,plot_level = 3,plot_options = NULL, CI_const = NULL,return_level = 2) # Approximate oblique profiles with UQ plot_options<-list(save=FALSE, titleProf = "Coordinate profiles", title2d = "Posterior mean",qq_fill=TRUE) options_sims<-list(nsim=150) obProfUQ<-obliqueProfiles(object=profMeans,threshold=threshold,allPsi = allPsi, options_full=options_full, options_approx=options_approx, uq_computations=TRUE, plot_level=3,plot_options=NULL, CI_const=NULL,return_level=2,options_sims=options_sims) ## End(Not run)
Function to plot the univariate profile extrema functions with UQ
plot_univariate_profiles_UQ( objectUQ, plot_options, nsims, threshold, nameFile = "prof_UQ", quantiles_uq = c(0.05, 0.95), profMean = NULL, typeProf = "approx" )plot_univariate_profiles_UQ( objectUQ, plot_options, nsims, threshold, nameFile = "prof_UQ", quantiles_uq = c(0.05, 0.95), profMean = NULL, typeProf = "approx" )
objectUQ |
an object returned by coordProf_UQ or the object saved in |
plot_options |
a list containing the same elements as the one passed to coordinateProfiles |
nsims |
number of simulations |
threshold |
threshold of interest |
nameFile |
the central name of the plot file |
quantiles_uq |
a vector containing the quantiles to be computed |
profMean |
the profile coordinate extrema functions for the mean. It is saved in |
typeProf |
a string to choose with type of profile extrema for simulations to plot
|
Plots either to the default graphical device or to pdf (according to the options passed in plot_options)
Plot bivariate profiles, for dimension up to 6.
plotBivariateProfiles( bivProf, allPsi, Design = NULL, threshold = NULL, whichIQR = NULL, plot_options = NULL, ... )plotBivariateProfiles( bivProf, allPsi, Design = NULL, threshold = NULL, whichIQR = NULL, plot_options = NULL, ... )
bivProf |
list returned by |
allPsi |
a list containing the matrices Psi (dim |
Design |
a matrix of dimension |
threshold |
if not |
whichIQR |
which quantiles to use for the inter-quantile range plot. If |
plot_options |
list as returned by |
... |
additional parameters to be passed to the plot function |
plots the 2d maps of the profile sup and inf of the function for each Psi in allPsi. If threshold is not NULL also contours the threshold level.
Dario Azzimonti
Plot coordinate profiles, for dimension up to 6.
plotMaxMin( allRes, Design = NULL, threshold = NULL, changes = FALSE, trueEvals = NULL, ... )plotMaxMin( allRes, Design = NULL, threshold = NULL, changes = FALSE, trueEvals = NULL, ... )
allRes |
list containing the list |
Design |
a d dimensional design corresponding to the points |
threshold |
if not |
changes |
boolean, if not |
trueEvals |
if not |
... |
additional parameters to be passed to the plot function |
plots the sup and inf of the function for each dimension. If threshold is not NULL
Dario Azzimonti
Auxiliary function for 2d plotting of excluded regions
plotOblique(changePoints, direction, ...)plotOblique(changePoints, direction, ...)
changePoints |
Numerical vector with the change points (usually if |
direction |
The Psi vector used for the direction |
... |
parameters to be passed to abline |
adds to the current plot the lines s.t. direction^T = changePoints[i] for all i
Dario Azzimonti
Plots the bivariate profiles stored in allRes for each Psi in allPsi.
plotOneBivProfile( allRes, allPsi, Design = NULL, threshold = NULL, trueEvals = NULL, main_addendum = "", ... )plotOneBivProfile( allRes, allPsi, Design = NULL, threshold = NULL, trueEvals = NULL, main_addendum = "", ... )
allRes |
list containing the list |
allPsi |
a list containing the matrices Psi (dim |
Design |
a matrix of dimension |
threshold |
if not |
trueEvals |
if not |
main_addendum |
additional string to add to image title. Default is empty string. |
... |
additional parameters to be passed to the plot function |
plots the 2d maps of the profile sup and inf in allRes for each Psi in allPsi. If threshold is not NULL also contours the threshold level.
Dario Azzimonti
plotBivariateProfiles
The function prof_mean_var_Delta computes the profile extrema functions for the mean and variance functions of the difference process at x.
prof_mean_var_Delta( kmModel, simupoints, allPsi = NULL, options_full_sims = NULL, options_approx = NULL, F.mat = NULL, T.mat = NULL )prof_mean_var_Delta( kmModel, simupoints, allPsi = NULL, options_full_sims = NULL, options_approx = NULL, F.mat = NULL, T.mat = NULL )
kmModel |
the km model of the Gaussian process |
simupoints |
the matrix |
allPsi |
optional list of matrices (dim |
options_full_sims |
an optional list of options for getAllMaxMin(or approxProfileExtrema if |
options_approx |
an optional list of options for approxMaxMin (or approxProfileExtrema if |
F.mat |
the evaluation of the trend function at |
T.mat |
the upper triangular factor of the Choleski decomposition of the covariance matrix of |
the profile extrema functions at options_approx$design for the mean and variance function of the difference process .
Dario Azzimonti
Computation and plots of profile extrema functions. The package main functions are:
coordinateProfiles: Given a km objects computes the coordinate profile extrema function for the posterior mean and its quantiles.
coordProf_UQ: UQ part of coordinateProfiles.
obliqueProfiles: Given a km objects computes the profile extrema functions for a generic list of matrices Psi for the posterior mean and its quantiles.
obliqueProf_UQ: The UQ part of obliqueProfiles.
getAllMaxMin: computes coordinate profile extrema with full optimization for a deterministic function.
approxMaxMin: approximates coordinate profile extrema for a deterministic function.
getProfileExtrema: computes profile extrema given a list of matrices Psi for a deterministic function.
approxProfileExtrema: approximates profile extrema given a list of matrices Psi for a deterministic function.
plot_univariate_profiles_UQ: plots for the results of coordProf_UQ or obliqueProf_UQ. Note that this function only works for univariate profiles.
plotBivariateProfiles: plots the bivariate maps results of a call to obliqueProfiles with a two dimensional projection matrix Psi.
plotMaxMin: simple plotting function for univariate profile extrema.
plotOneBivProfile: simple plotting function for bivariate profile extrema.
Package: profExtrema
Type: Package
Version: 0.2.1
Date: 2020-03-20
This work was supported in part the Hasler Foundation, grant number 16065 and by the Swiss National Science Foundation, grant number 167199. The author warmly thanks David Ginsbourger, Jérémy Rohmer and Déborah Idier for fruitful discussions and accurate, thought provoking suggestions.
Dario Azzimonti ([email protected]) .
Azzimonti, D., Bect, J., Chevalier, C., and Ginsbourger, D. (2016). Quantifying uncertainties on excursion sets under a Gaussian random field prior. SIAM/ASA Journal on Uncertainty Quantification, 4(1):850–874.
Azzimonti, D., Ginsbourger, D., Rohmer, J. and Idier, D. (2017+). Profile extrema for visualizing and quantifying uncertainties on excursion regions. Application to coastal flooding. arXiv:1710.00688.
Chevalier, C. (2013). Fast uncertainty reduction strategies relying on Gaussian process models. PhD thesis, University of Bern.
Chevalier, C., Picheny, V., Ginsbourger, D. (2014). An efficient and user-friendly implementation of batch-sequential inversion strategies based on kriging. Computational Statistics & Data Analysis, 71: 1021-1034.
Johnson, S. G. The NLopt nonlinear-optimization package, http://ab-initio.mit.edu/nlopt
Koenker, R. (2017). quantreg: Quantile Regression. R package version 5.33.
Nocedal, J. and Wright, S. J. (2006). Numerical Optimization, second edition. Springer- Verlag, New York.
Neuwirth, E. (2014). RColorBrewer: ColorBrewer Palettes. R package version 1.1-2.
Roustant, O., Ginsbourger, D., Deville, Y. (2012). DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization. Journal of Statistical Software, 51(1): 1-55.
Function to set-up plot options for plot_univariate_profiles_UQ, plotBivariateProfiles, coordinateProfiles, coordProf_UQ, obliqueProfiles and obliqueProf_UQ.
setPlotOptions(plot_options = NULL, d, num_T, kmModel = NULL)setPlotOptions(plot_options = NULL, d, num_T, kmModel = NULL)
plot_options |
the list of plot options to set-up |
d |
number of coordinates |
num_T |
number of thresholds of interest |
kmModel |
a km model, used to obtain the coordinates names. |
the properly set-up list containing the following fields
save: boolean, if TRUE saves the plots in folderPlots
folderPlots: a string containing the destination folder for plots, if save==TRUE default is ./
ylim: a matrix coordx2 containing the ylim for each coordinate, if NULL in plot_options this is left NULL and automatically set at the plot time.
titleProf: a string containing the title for the coordinate profile plots, default is "Coordinate profiles"
title2d: a string containing the title for the 2d plots (if the input is 2d), default is "Posterior mean"
design: a matrix where is the input dimension and is the size of the discretization for plots at each dimension
coord_names: a -vector of characters naming the dimensions. If NULL and kmModel not NULL then it is the names of kmModel@X otherwise x_1,...,x_d
id_save: a string to be added to the plot file names, useful for serial computations on HPC, left as in plot_options.
qq_fill: if TRUE it fills the region between the first 2 quantiles in quantiles_uq and between the upper and lower bound in objectUQ$bound$bound, if NULL, it is set as FALSE.
bound_cols: a vector of two strings containing the names of the colors for upper and lower bound plots.
qq_fill_colors: a list containing the colors for qq_fill: approx for 2 quantiles, bound_min for bounds on the profile inf, bound_max for profile sup. Initialized only if qq_fill==TRUE.
col_CCPthresh_nev: Color palette of dimension num_T for the colors of the vertical lines delimiting the intersections between the profiles sup and the thresholds
col_CCPthresh_alw: Color palette of dimension num_T for the colors of the vertical lines delimiting the intersections between the profiles inf and the thresholds
col_thresh: Color palette of dimension num_T for the colors of the thresholds
fun_evals: integer denoting the level of plot for the true evaluations.
0: default, no plots for true evaluations;
1: plot the true evaluations as points in 2d plots, no true evaluation plots in 1d;
2: plot true evaluations, in 2d with different color for values above threshold;
3: plot true evaluations, in 2d plots in color, with background of the image colored as proportion of points inside excursion;
if all the fields are already filled then returns plot_options
Dario Azzimonti
The function var_Delta_T computes the gradient for the variance function of the difference process at x.
var_Delta_T(x, kmModel, simupoints, T.mat, F.mat)var_Delta_T(x, kmModel, simupoints, T.mat, F.mat)
x |
a matrix |
kmModel |
the km model of the Gaussian process |
simupoints |
the matrix |
T.mat |
the upper triangular factor of the Choleski decomposition of the covariance matrix of |
F.mat |
the evaluation of the trend function at |
the value of the variance function at x for the difference process .
Dario Azzimonti