ARCHISSUR (version 0.0.1)

archissur: Run ARCHISSUR algorithm

Description

archissur adaptively enriches the Gaussian Process Classifier (GPC) using a learning criterion to achieve a precise approximation of the feasible area contour. This is done by iteratively adding the best learning point that minimizes future uncertainty over the feasible domain, following the Stepwise Uncertainty Reduction strategy (SUR).

Usage

archissur(
  design.init = NULL,
  cst.init = NULL,
  model = NULL,
  cst_function,
  lower = NULL,
  upper = NULL,
  n.ite = 10,
  seed = NULL,
  nb.integration = NULL,
  plot_2D_pn = FALSE,
  batchsize = 1,
  n_update = 1,
  gpc.options = NULL,
  optimcontrol = NULL,
  verbose = 1
)

Value

A list containing:

Xf

DoE

f

binary observations corresponding to the class labels.

alpha

a scalar representing the Vorob'ev threshold.

vorob_expect

Vorob’ev expectation.

vorob_dev

current Vorob’ev deviation.

model

an object of class gpcm containing the GPC model at last iteration.

An '.Rds' file model_gp.Rds containing an object of class gpcm corresponding to current GPC model. The class 1 probability map. Plots are available in '2D_plots' directory.

Arguments

design.init

optional matrix representing the initial design of experiments (DoE). If not provided, you must provide a model of type gpcm.

cst.init

optional vector of binary observations {0,1} corresponding to the initial class labels. If not provided, it will be calculated using cst_function.

model

optional object of type gpcm to start archissur. If not provided, you must provide an initial DoE design.init.

cst_function

constraint function with binary outputs {0,1} to be learn.

lower

inputs lower bound of design.init.

upper

inputs upper bound of design.init.

n.ite

number of iterations of archissur.

seed

to fix the seed.

nb.integration

number of integration points. Default is d*1000.

plot_2D_pn

if TRUE and d = 2, plot class 1 probability map and learning points. Plots are available in '2D_plots' directory. Default is FALSE.

batchsize

number of points to be learned at each iteration. Default is 1.

n_update

number of iterations between hyperparameter updates by likelihood maximization. Default is 1.

gpc.options

list with GPC model options: covariance kernel, noise variance, number of initial points for MLE optimization, standardization of inputs, constrained latent GP mean sign.If NULL, default options are list(normalize = T, multistart = 1, covtype = "matern5_2", MeanTransform=NULL). Default is NULL.

optimcontrol

optional list of control parameters for enrichment criterion optimization. Default is NULL. The field "method" defines which optimization method is used: it can be either "multistart" (default) for an optimization with L-BFGS-B multistart, or "discrete" for an optimization over a specified discrete set, or "genoud" for an optimization using the genoud algorithm. (See details).

verbose

Level of verbosity for printing information during iterations. 0: No printing. 1: Print iteration number and best point found. 2: Print iteration number, best point found, criterion value, and model hyperparameters. Default is 1.

Details

If the field "method" is set to "genoud", one can set some parameters of this algorithm: pop.size (default: 50d), max.generations (10d), wait.generations (2), BFGSburnin (2) and the mutations P1, P2, up to P9 (see genoud). Numbers into brackets are the default values. If the field method is set to "discrete", one can set the field optim.points: p * d matrix corresponding to the p points where the criterion will be evaluated. If nothing is specified, 100*d points are chosen randomly. Finally, one can control the field optim.option in order to decide how to optimize the sampling criterion. If optim.option is set to 2 (default), batchsize sequential optimizations in dimension d are performed to find the optimum. If optim.option is set to 1, only one optimization in dimension batchsize*d is performed. This option is only available with "genoud". This option might provide more global and accurate solutions, but is a lot more expensive.

References

Menz, M., Munoz-Zuniga, M., Sinoquet, D. Estimation of simulation failure set with active learning based on Gaussian Process classifiers and random set theory (2023). https://hal.science/hal-03848238.

Bachoc, F., Helbert, C. & Picheny, V. Gaussian process optimization with failures: classification and convergence proof. J Glob Optim 78, 483–506 (2020). tools:::Rd_expr_doi("10.1007/s10898-020-00920-0").

Examples

# \donttest{
#-------------------------------------------------------------------
#------------------------- archissur -------------------------------
#-------------------------------------------------------------------

#  20-points DoE, and the corresponding response
d <- 2
nb_PX <- 20
x <- matrix(c(0.205293785978832, 0.0159983370750337,
              0.684774733109666, 0.125251417595962,
              0.787208786290006, 0.700475706055049,
              0.480507717105934, 0.359730889653793,
              0.543665267336735, 0.565974761807069,
              0.303412043992361, 0.471502352650857,
              0.839505250127309, 0.504914690245002,
              0.573294917143728, 0.784444726564573,
              0.291681289223421, 0.255053812451938,
              0.87233450888786, 0.947168337730927,
              0.648262257638515, 0.973264712407035,
              0.421877310273815, 0.0686662506387988,
              0.190976166753807, 0.810964668176754,
              0.918527262507395, 0.161973686467513,
              0.0188128700859558, 0.43522031347403,
              0.99902788789426, 0.655561821513544,
              0.741113863862512, 0.321050086076934,
              0.112003007565305, 0.616551317575545,
              0.383511473487687, 0.886611679106771,
              0.0749211435982952, 0.205805968972305),
            byrow = TRUE, ncol = d)

require(DiceKriging)
cst_function <- function(z){
  fx <- apply(z, 1, branin)
  f <- ifelse(fx < 14, 0, 1)
  return(f)}

## constraint function
s <- cst_function(x)

# archissur parameters

design.init <- x
cst.init <- s
n.ite <- 2
n_update <- 5
lower <- rep(0,d)
upper <- rep(1,d)
### GPC model options
gpc.options <- list()
gpc.options$noise.var <- 1e-6
gpc.options$multistart <- 1

res <-  archissur(design.init = design.init, cst.init = cst.init,
                  cst_function = cst_function, lower = lower, upper = upper,
                  n.ite = n.ite, n_update = n_update,  gpc.options = gpc.options)
unlink('model_gp.Rds')
# }