dgpsi: Interface to 'dgpsi' for Deep and Linked Gaussian Process Emulations

Description

Interface to the 'python' package 'dgpsi' for Gaussian process, deep Gaussian process, and linked deep Gaussian process emulations of computer models and networks using stochastic imputation (SI). The implementations follow Ming & Guillas (2021) doi:10.1137/20M1323771 and Ming, Williamson, & Guillas (2023) doi:10.1080/00401706.2022.2124311 and Ming & Williamson (2023) doi:10.48550/arXiv.2306.01212. To get started with the package, see https://mingdeyu.github.io/dgpsi-R/.

Author(s)

Maintainer: Deyu Ming deyu.ming.16@ucl.ac.uk [copyright holder]

Authors:

• Daniel Williamson d.williamson@exeter.ac.uk

See Also

Useful links:

• https://github.com/mingdeyu/dgpsi-R

• https://mingdeyu.github.io/dgpsi-R/

• Report bugs at https://github.com/mingdeyu/dgpsi-R/issues

Locate the next design point(s) for a (D)GP emulator or a bundle of (D)GP emulators using Active Learning MacKay (ALM)

Description

This function searches from a candidate set to locate the next design point(s) to be added to a (D)GP emulator or a bundle of (D)GP emulators using the Active Learning MacKay (ALM) criterion (see the reference below).

Usage

alm(object, ...)

## S3 method for class 'gp'
alm(
  object,
  x_cand = NULL,
  n_start = 20,
  batch_size = 1,
  M = 50,
  workers = 1,
  limits = NULL,
  int = FALSE,
  ...
)

## S3 method for class 'dgp'
alm(
  object,
  x_cand = NULL,
  n_start = 20,
  batch_size = 1,
  M = 50,
  workers = 1,
  limits = NULL,
  int = FALSE,
  aggregate = NULL,
  ...
)

## S3 method for class 'bundle'
alm(
  object,
  x_cand = NULL,
  n_start = 20,
  batch_size = 1,
  M = 50,
  workers = 1,
  limits = NULL,
  int = FALSE,
  aggregate = NULL,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

1. If x_cand is not NULL:

   • When object is an instance of the gp class, a vector of length batch_size is returned, containing the positions (row numbers) of the next design points from x_cand.

   • When object is an instance of the dgp class, a vector of length batch_size * D is returned, containing the positions (row numbers) of the next design points from x_cand to be added to the DGP emulator.

     • D is the number of output dimensions of the DGP emulator if no likelihood layer is included.

     • For a DGP emulator with a Hetero or NegBin likelihood layer, D = 2.

     • For a DGP emulator with a Categorical likelihood layer, D = 1 for binary output or D = K for multi-class output with K classes.

   • When object is an instance of the bundle class, a matrix is returned with batch_size rows and a column for each emulator in the bundle, containing the positions (row numbers) of the next design points from x_cand for individual emulators.

2. If x_cand is NULL:

   • When object is an instance of the gp class, a matrix with batch_size rows is returned, giving the next design points to be evaluated.

   • When object is an instance of the dgp class, a matrix with batch_size * D rows is returned, where:

     • D is the number of output dimensions of the DGP emulator if no likelihood layer is included.

     • For a DGP emulator with a Hetero or NegBin likelihood layer, D = 2.

     • For a DGP emulator with a Categorical likelihood layer, D = 1 for binary output or D = K for multi-class output with K classes.

   • When object is an instance of the bundle class, a list is returned with a length equal to the number of emulators in the bundle. Each element of the list is a matrix with batch_size rows, where each row represents a design point to be added to the corresponding emulator.

Note

The first column of the matrix supplied to the first argument of aggregate must correspond to the first output dimension of the DGP emulator if object is an instance of the dgp class, and so on for subsequent columns and dimensions. If object is an instance of the bundle class, the first column must correspond to the first emulator in the bundle, and so on for subsequent columns and emulators.

References

MacKay, D. J. (1992). Information-based objective functions for active data selection. Neural Computation, 4(4), 590-604.

Examples

## Not run: 

# load packages and the Python env
library(lhs)
library(dgpsi)

# construct a 1D non-stationary function
f <- function(x) {
 sin(30*((2*x-1)/2-0.4)^5)*cos(20*((2*x-1)/2-0.4))
}

# generate the initial design
X <- maximinLHS(10,1)
Y <- f(X)

# training a 2-layered DGP emulator with the global connection off
m <- dgp(X, Y, connect = F)

# specify the input range
lim <- c(0,1)

# locate the next design point using ALM
X_new <- alm(m, limits = lim)

# obtain the corresponding output at the located design point
Y_new <- f(X_new)

# combine the new input-output pair to the existing data
X <- rbind(X, X_new)
Y <- rbind(Y, Y_new)

# update the DGP emulator with the new input and output data and refit
m <- update(m, X, Y, refit = TRUE)

# plot the LOO validation
plot(m)

## End(Not run)

Combine layers

Description

This function is deprecated and will be removed in the next release, as it is simply a wrapper for the list() function. To construct linked (D)GP structures, please use the updated lgp() function, which provides a simpler and more efficient approach to building (D)GP emulators.

Usage

combine(...)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

A list defining a linked (D)GP structure to be passed to struc of lgp().

Continue training a DGP emulator

Description

This function implements additional training iterations for a DGP emulator.

Usage

continue(
  object,
  N = NULL,
  cores = 1,
  ess_burn = 10,
  verb = TRUE,
  burnin = NULL,
  B = NULL
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object.

Note

• One can also use this function to fit an untrained DGP emulator constructed by dgp() with training = FALSE.

• The following slots:

  • loo and oos created by validate(); and

  • results created by predict() in object will be removed and not contained in the returned object.

Examples

## Not run: 

# See dgp() for an example.

## End(Not run)

Restore the serialized emulator

Description

This function restores the serialized emulator created by serialize().

Usage

deserialize(object)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

The S3 class of a GP emulator, a DGP emulator, a linked (D)GP emulator, or a bundle of (D)GP emulators.

Note

See the Note section in serialize().

Examples

## Not run: 

library(future)
library(future.apply)
library(dgpsi)

# model
f <- function(x) {
 (sin(7.5*x)+1)/2
}

# training data
X <- seq(0, 1, length = 10)
Y <- sapply(X, f)

# train a DGP emulator
m <- dgp(X, Y, name = "matern2.5")

# testing input data
X_dgp <- seq(0, 1, length = 100)

# serialize the DGP emulator
m_serialized <- serialize(m)

# start a multi-session with three cores for parallel predictions
plan(multisession, workers = 3)

# perform parallel predictions
results <- future_lapply(1:length(X_dgp), function(i) {
  m_deserialized <- deserialize(m_serialized)
  mean_i <- predict(m_deserialized, X_dgp[i])$results$mean
}, future.seed = TRUE)

# reset the future plan to sequential
plan(sequential)

# combine mean predictions
pred_mean <- do.call(rbind, results)

## End(Not run)

Sequential design of a (D)GP emulator or a bundle of (D)GP emulators

Description

This function implements sequential design and active learning for a (D)GP emulator or a bundle of (D)GP emulators, supporting an array of popular methods as well as user-specified approaches. It can also be used as a wrapper for Bayesian optimization methods.

Usage

design(
  object,
  N,
  x_cand,
  y_cand,
  n_sample,
  n_cand,
  limits,
  f,
  reps,
  freq,
  x_test,
  y_test,
  reset,
  target,
  method,
  batch_size,
  eval,
  verb,
  autosave,
  new_wave,
  M_val,
  cores,
  ...
)

## S3 method for class 'gp'
design(
  object,
  N,
  x_cand = NULL,
  y_cand = NULL,
  n_sample = 200,
  n_cand = lifecycle::deprecated(),
  limits = NULL,
  f = NULL,
  reps = 1,
  freq = c(1, 1),
  x_test = NULL,
  y_test = NULL,
  reset = FALSE,
  target = NULL,
  method = vigf,
  batch_size = 1,
  eval = NULL,
  verb = TRUE,
  autosave = list(),
  new_wave = TRUE,
  M_val = 50,
  cores = 1,
  ...
)

## S3 method for class 'dgp'
design(
  object,
  N,
  x_cand = NULL,
  y_cand = NULL,
  n_sample = 200,
  n_cand = lifecycle::deprecated(),
  limits = NULL,
  f = NULL,
  reps = 1,
  freq = c(1, 1),
  x_test = NULL,
  y_test = NULL,
  reset = FALSE,
  target = NULL,
  method = vigf,
  batch_size = 1,
  eval = NULL,
  verb = TRUE,
  autosave = list(),
  new_wave = TRUE,
  M_val = 50,
  cores = 1,
  train_N = NULL,
  refit_cores = 1,
  pruning = TRUE,
  control = list(),
  ...
)

## S3 method for class 'bundle'
design(
  object,
  N,
  x_cand = NULL,
  y_cand = NULL,
  n_sample = 200,
  n_cand = lifecycle::deprecated(),
  limits = NULL,
  f = NULL,
  reps = 1,
  freq = c(1, 1),
  x_test = NULL,
  y_test = NULL,
  reset = FALSE,
  target = NULL,
  method = vigf,
  batch_size = 1,
  eval = NULL,
  verb = TRUE,
  autosave = list(),
  new_wave = TRUE,
  M_val = 50,
  cores = 1,
  train_N = NULL,
  refit_cores = 1,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object is returned with a slot called design that contains:

• S slots, named ⁠wave1, wave2,..., waveS⁠, that contain information of S waves of sequential design that have been applied to the emulator. Each slot contains the following elements:

  • N, an integer that gives the numbers of iterations implemented in the corresponding wave;

  • rmse, a matrix providing the evaluation metric values for emulators constructed during the corresponding wave, when eval = NULL. Each row of the matrix represents an iteration.

    • for an object of class gp, the matrix contains a single column of RMSE values.

    • for an object of class dgp without a categorical likelihood, each row contains mean/median squared errors corresponding to different output dimensions.

    • for an object of class dgp with a categorical likelihood, the matrix contains a single column of log-loss values.

    • for an object of class bundle, each row contains either mean/median squared errors or log-loss values for the emulators in the bundle.

  • metric: a matrix providing the values of custom evaluation metrics, as computed by the user-supplied eval function, for emulators constructed during the corresponding wave.

  • freq, an integer that gives the frequency that the emulator validations are implemented during the corresponding wave.

  • enrichment, a vector of size N that gives the number of new design points added after each step of the sequential design (if object is an instance of the gp or dgp class), or a matrix that gives the number of new design points added to emulators in a bundle after each step of the sequential design (if object is an instance of the bundle class).

  If target is not NULL, the following additional elements are also included:

  • target: the target evaluating metric computed by the eval or built-in function to stop the sequential design.

  • reached: indicates whether the target was reached at the end of the sequential design:

    • a bool if object is an instance of the gp or dgp class.

    • a vector of bools if object is an instance of the bundle class, with its length determined as follows:

      • equal to the number of emulators in the bundle when eval = NULL.

      • equal to the length of the output from eval when a custom eval function is provided.

• a slot called type that gives the type of validation:

  • either LOO ('loo') or OOS ('oos') if eval = NULL. See validate() for more information about LOO and OOS.

  • 'customized' if a customized R function is provided to eval.

• two slots called x_test and y_test that contain the data points for the OOS validation if the type slot is 'oos'.

• If y_cand = NULL and x_cand is supplied, and there are NAs returned from the supplied f during the sequential design, a slot called exclusion is included that records the located design positions that produced NAs via f. The sequential design will use this information to avoid re-visiting the same locations in later runs of design().

See Note section below for further information.

Note

• Validation of an emulator is forced after the final step of a sequential design even if N is not a multiple of the second element in freq.

• Any loo or oos slot that already exists in object will be cleaned, and a new slot called loo or oos will be created in the returned object depending on whether x_test and y_test are provided. The new slot gives the validation information of the emulator constructed in the final step of the sequential design. See validate() for more information about the slots loo and oos.

• If object has previously been used by design() for sequential design, the information of the current wave of the sequential design will replace those of old waves and be contained in the returned object, unless

  • the validation type (LOO or OOS depending on whether x_test and y_test are supplied or not) of the current wave of the sequential design is the same as the validation types (shown in the type of the design slot of object) in previous waves, and if the validation type is OOS, x_test and y_test in the current wave must also be identical to those in the previous waves;

  • both the current and previous waves of the sequential design supply customized evaluation functions to eval. Users need to ensure the customized evaluation functions are consistent among different waves. Otherwise, the trace plot of RMSEs produced by draw() will show values of different evaluation metrics in different waves.

  For the above two cases, the information of the current wave of the sequential design will be added to the design slot of the returned object under the name waveS.

• If object is an instance of the gp class and eval = NULL, the matrix in the rmse slot is single-columned. If object is an instance of the dgp or bundle class and eval = NULL, the matrix in the rmse slot can have multiple columns that correspond to different output dimensions or different emulators in the bundle.

• If object is an instance of the gp class and eval = NULL, target needs to be a single value giving the RMSE threshold. If object is an instance of the dgp or bundle class and eval = NULL, target can be a vector of values that gives the thresholds of evaluating metrics for different output dimensions or different emulators. If a single value is provided, it will be used as the threshold for all output dimensions (if object is an instance of the dgp) or all emulators (if object is an instance of the bundle). If a customized function is supplied to eval and target is given as a vector, the user needs to ensure that the length of target is equal to that of the output from eval.

• When defining f, it is important to ensure that:

  • the column order of the first argument of f is consistent with the training input used for the emulator;

  • the column order of the output matrix of f is consistent with the order of emulator output dimensions (if object is an instance of the dgp class), or the order of emulators placed in object (if object is an instance of the bundle class).

• The output matrix produced by f may include NAs. This is especially beneficial as it allows the sequential design process to continue without interruption, even if errors or NA outputs are encountered from f at certain input locations identified by the sequential design. Users should ensure that any errors within f are handled by appropriately returning NAs.

• When defining eval, the output metric needs to be positive if draw() is used with log = T. And one needs to ensure that a lower metric value indicates a better emulation performance if target is set.

Examples

## Not run: 

# load packages and the Python env
library(lhs)
library(dgpsi)

# construct a 2D non-stationary function that takes a matrix as the input
f <- function(x) {
  sin(1/((0.7*x[,1,drop=F]+0.3)*(0.7*x[,2,drop=F]+0.3)))
}

# generate the initial design
X <- maximinLHS(5,2)
Y <- f(X)

# generate the validation data
validate_x <- maximinLHS(30,2)
validate_y <- f(validate_x)

# training a 2-layered DGP emulator with the initial design
m <- dgp(X, Y)

# specify the ranges of the input dimensions
lim_1 <- c(0, 1)
lim_2 <- c(0, 1)
lim <- rbind(lim_1, lim_2)

# 1st wave of the sequential design with 10 steps
m <- design(m, N=10, limits = lim, f = f, x_test = validate_x, y_test = validate_y)

# 2nd wave of the sequential design with 10 steps
m <- design(m, N=10, limits = lim, f = f, x_test = validate_x, y_test = validate_y)

# 3rd wave of the sequential design with 10 steps
m <- design(m, N=10, limits = lim, f = f, x_test = validate_x, y_test = validate_y)

# draw the design created by the sequential design
draw(m,'design')

# inspect the trace of RMSEs during the sequential design
draw(m,'rmse')

# reduce the number of imputations for faster OOS
m_faster <- set_imp(m, 5)

# plot the OOS validation with the faster DGP emulator
plot(m_faster, x_test = validate_x, y_test = validate_y)

## End(Not run)

Deep Gaussian process emulator construction

Description

This function builds and trains a DGP emulator.

Usage

dgp(
  X,
  Y,
  depth = 2,
  node = ncol(X),
  name = "sexp",
  lengthscale = 1,
  bounds = NULL,
  prior = "ga",
  share = TRUE,
  nugget_est = FALSE,
  nugget = NULL,
  scale_est = TRUE,
  scale = 1,
  connect = TRUE,
  likelihood = NULL,
  training = TRUE,
  verb = TRUE,
  check_rep = TRUE,
  vecchia = FALSE,
  M = 25,
  ord = NULL,
  N = ifelse(vecchia, 200, 500),
  cores = 1,
  blocked_gibbs = TRUE,
  ess_burn = 10,
  burnin = NULL,
  B = 10,
  internal_input_idx = NULL,
  linked_idx = NULL,
  id = NULL
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An S3 class named dgp that contains five slots:

• id: A number or character string assigned through the id argument.

• data: a list that contains two elements: X and Y which are the training input and output data respectively.

• specs: a list that contains

  1. L (i.e., the number of layers in the DGP hierarchy) sub-lists named ⁠layer1, layer2,..., layerL⁠. Each sub-list contains D (i.e., the number of GP/likelihood nodes in the corresponding layer) sub-lists named ⁠node1, node2,..., nodeD⁠. If a sub-list corresponds to a likelihood node, it contains one element called type that gives the name (Hetero, Poisson, NegBin, or Categorical) of the likelihood node. If a sub-list corresponds to a GP node, it contains four elements:

     • kernel: the type of the kernel function used for the GP node.

     • lengthscales: a vector of lengthscales in the kernel function.

     • scale: the variance value in the kernel function.

     • nugget: the nugget value in the kernel function.

  2. internal_dims: the column indices of X that correspond to the linked emulators in the preceding layers of a linked system. The slot will be removed in the next release.

  3. external_dims: the column indices of X that correspond to global inputs to the linked system of emulators. It is shown as FALSE if internal_input_idx = NULL. The slot will be removed in the next release.

  4. linked_idx: the value passed to argument linked_idx. It is shown as FALSE if the argument linked_idx is NULL. The slot will be removed in the next release.

  5. seed: the random seed generated to produce imputations. This information is stored for reproducibility when the DGP emulator (that was saved by write() with the light option light = TRUE) is loaded back to R by read().

  6. B: the number of imputations used to generate the emulator.

  7. vecchia: whether the Vecchia approximation is used for the GP emulator training.

  8. M: the size of the conditioning set for the Vecchia approximation in the DGP emulator training. M is generated only when vecchia = TRUE.

• constructor_obj: a 'python' object that stores the information of the constructed DGP emulator.

• container_obj: a 'python' object that stores the information for the linked emulation.

• emulator_obj: a 'python' object that stores the information for the predictions from the DGP emulator.

The returned dgp object can be used by

• predict() for DGP predictions.

• continue() for additional DGP training iterations.

• validate() for LOO and OOS validations.

• plot() for validation plots.

• lgp() for linked (D)GP emulator constructions.

• window() for model parameter trimming.

• summary() to summarize the trained DGP emulator.

• write() to save the DGP emulator to a .pkl file.

• set_imp() to change the number of imputations.

• design() for sequential design.

• update() to update the DGP emulator with new inputs and outputs.

• alm(), mice(), and vigf() to locate next design points.

Note

Any R vector detected in X and Y will be treated as a column vector and automatically converted into a single-column R matrix. Thus, if X is a single data point with multiple dimensions, it must be given as a matrix.

Examples

## Not run: 

# load the package and the Python env
library(dgpsi)

# construct a step function
f <- function(x) {
  if (x < 0.5) return(-1)
  if (x >= 0.5) return(1)
  }

# generate training data
X <- seq(0, 1, length = 10)
Y <- sapply(X, f)

# set a random seed
set_seed(999)

# training a DGP emulator
m <- dgp(X, Y)

# continue for further training iterations
m <- continue(m)

# summarizing
summary(m)

# trace plot
trace_plot(m)

# trim the traces of model parameters
m <- window(m, 800)

# LOO cross validation
m <- validate(m)
plot(m)

# prediction
test_x <- seq(0, 1, length = 200)
m <- predict(m, x = test_x)

# OOS validation
validate_x <- sample(test_x, 10)
validate_y <- sapply(validate_x, f)
plot(m, validate_x, validate_y)

# write and read the constructed emulator
write(m, 'step_dgp')
m <- read('step_dgp')

## End(Not run)

Validation and diagnostic plots for a sequential design

Description

This function draws diagnostic and validation plots for a sequential design of a (D)GP emulator or a bundle of (D)GP emulators.

Usage

draw(object, ...)

## S3 method for class 'gp'
draw(object, type = "rmse", log = FALSE, ...)

## S3 method for class 'dgp'
draw(object, type = "rmse", log = FALSE, ...)

## S3 method for class 'bundle'
draw(object, type = "rmse", log = FALSE, emulator = NULL, ...)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

A patchwork object.

Examples

## Not run: 

# See design() for an example.

## End(Not run)

Get the number of threads

Description

This function gets the number of threads used for parallel computations involved in the package.

Usage

get_thread_num()

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

the number of threads.

Gaussian process emulator construction

Description

This function builds and trains a GP emulator.

Usage

gp(
  X,
  Y,
  name = "sexp",
  lengthscale = rep(0.1, ncol(X)),
  bounds = NULL,
  prior = "ref",
  nugget_est = FALSE,
  nugget = ifelse(nugget_est, 0.01, 1e-08),
  scale_est = TRUE,
  scale = 1,
  training = TRUE,
  verb = TRUE,
  vecchia = FALSE,
  M = 25,
  ord = NULL,
  internal_input_idx = NULL,
  linked_idx = NULL,
  id = NULL
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An S3 class named gp that contains five slots:

• id: A number or character string assigned through the id argument.

• data: a list that contains two elements: X and Y which are the training input and output data respectively.

• specs: a list that contains seven elements:

  1. kernel: the type of the kernel function used. Either "sexp" for squared exponential kernel or "matern2.5" for Matérn-2.5 kernel.

  2. lengthscales: a vector of lengthscales in the kernel function.

  3. scale: the variance value in the kernel function.

  4. nugget: the nugget value in the kernel function.

  5. internal_dims: the column indices of X that correspond to the linked emulators in the preceding layers of a linked system. The slot will be removed in the next release.

  6. external_dims: the column indices of X that correspond to global inputs to the linked system of emulators. It is shown as FALSE if internal_input_idx = NULL. The slot will be removed in the next release.

  7. linked_idx: the value passed to argument linked_idx. It is shown as FALSE if the argument linked_idx is NULL. The slot will be removed in the next release.

  8. vecchia: whether the Vecchia approximation is used for the GP emulator training.

  9. M: the size of the conditioning set for the Vecchia approximation in the GP emulator training.

• constructor_obj: a 'python' object that stores the information of the constructed GP emulator.

• container_obj: a 'python' object that stores the information for the linked emulation.

• emulator_obj: a 'python' object that stores the information for the predictions from the GP emulator.

The returned gp object can be used by

• predict() for GP predictions.

• validate() for LOO and OOS validations.

• plot() for validation plots.

• lgp() for linked (D)GP emulator constructions.

• summary() to summarize the trained GP emulator.

• write() to save the GP emulator to a .pkl file.

• design() for sequential designs.

• update() to update the GP emulator with new inputs and outputs.

• alm(), mice(), and vigf() to locate next design points.

Note

Any R vector detected in X and Y will be treated as a column vector and automatically converted into a single-column R matrix. Thus, if X is a single data point with multiple dimensions, it must be given as a matrix.

References

• Gu, M. (2019). Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection. Bayesian Analysis, 14(3), 857-885.

• Katzfuss, M., Guinness, J., & Lawrence, E. (2022). Scaled Vecchia approximation for fast computer-model emulation. SIAM/ASA Journal on Uncertainty Quantification, 10(2), 537-554.

Examples

## Not run: 
# load the package and the Python env
library(dgpsi)

# construct a step function
f <- function(x) {
   if (x < 0.5) return(-1)
   if (x >= 0.5) return(1)
  }

# generate training data
X <- seq(0, 1, length = 10)
Y <- sapply(X, f)

# training
m <- gp(X, Y)

# summarizing
summary(m)

# LOO cross validation
m <- validate(m)
plot(m)

# prediction
test_x <- seq(0, 1, length = 200)
m <- predict(m, x = test_x)

# OOS validation
validate_x <- sample(test_x, 10)
validate_y <- sapply(validate_x, f)
plot(m, validate_x, validate_y)

# write and read the constructed emulator
write(m, 'step_gp')
m <- read('step_gp')

## End(Not run)

'python' environment initialization

Description

This function initializes the 'python' environment for the package.

Usage

init_py(
  py_ver = NULL,
  dgpsi_ver = NULL,
  reinstall = FALSE,
  uninstall = FALSE,
  verb = TRUE
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

No return value, called to install required 'python' environment.

Examples

## Not run: 

# See gp(), dgp(), or lgp() for an example.

## End(Not run)

Linked (D)GP emulator construction

Description

This function constructs a linked (D)GP emulator for a model chain or network.

Usage

lgp(struc, emulators = NULL, B = 10, activate = TRUE, verb = TRUE, id = NULL)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An S3 class named lgp that contains three slots:

• id: A number or character string assigned through the id argument.

• constructor_obj: a list of 'python' objects that stores the information of the constructed linked emulator.

• emulator_obj, a 'python' object that stores the information for predictions from the linked emulator.

• specs: a list that contains

  1. seed: the random seed generated to produce the imputations. This information is stored for reproducibility when the linked (D)GP emulator (that was saved by write() with the light option light = TRUE) is loaded back to R by read().

  2. B: the number of imputations used to generate the linked (D)GP emulator.

  If struc is a data frame, specs also includes:

  1. metadata: a data frame providing configuration details for each emulator in the linked system, with following columns:

     • Emulator: the ID of an emulator.

     • Layer: the layer in the linked system where the emulator is positioned. A lower Layer number indicates a position closer to the input, with layer numbering increasing as you move away from the input.

     • Pos_in_Layer: the position of the emulator within its layer. A lower Pos_in_Layer number indicates a position higher up in that layer.

     • Total_Input_Dims: the total number of input dimensions of the emulator.

     • Total_Output_Dims: the total number of output dimensions of the emulator.

  2. struc: The linked system structure, as supplied by struc.

The returned lgp object can be used by

• predict() for linked (D)GP predictions.

• validate() for OOS validation.

• plot() for validation plots.

• summary() to summarize the constructed linked (D)GP emulator.

• write() to save the linked (D)GP emulator to a .pkl file.

Examples

## Not run: 

# load the package and the Python env
library(dgpsi)

# model 1
f1 <- function(x) {
  (sin(7.5*x)+1)/2
}
# model 2
f2 <- function(x) {
  2/3*sin(2*(2*x - 1))+4/3*exp(-30*(2*(2*x-1))^2)-1/3
}
# linked model
f12 <- function(x) {
  f2(f1(x))
}

# training data for Model 1
X1 <- seq(0, 1, length = 9)
Y1 <- sapply(X1, f1)
# training data for Model 2
X2 <- seq(0, 1, length = 13)
Y2 <- sapply(X2, f2)

# emulation of model 1
m1 <- gp(X1, Y1, name = "matern2.5", id = "emulator1")
# emulation of model 2
m2 <- dgp(X2, Y2, depth = 2, name = "matern2.5", id = "emulator2")

struc <- data.frame(From_Emulator = c("Global", "emulator1"),
                    To_Emulator = c("emulator1", "emulator2"),
                    From_Output = c(1, 1),
                    To_Input = c(1, 1))
emulators <- list(m1, m2)

# construct the linked emulator for visual inspection
m_link <- lgp(struc, emulators, activate = FALSE)

# visual inspection
summary(m_link)

# build the linked emulator for prediction
m_link <- lgp(struc, emulators, activate = TRUE)
test_x <- seq(0, 1, length = 300)
m_link <- predict(m_link, x = test_x)

# OOS validation
validate_x <- sample(test_x, 20)
validate_y <- sapply(validate_x, f12)
plot(m_link, validate_x, validate_y, style = 2)

# write and read the constructed linked emulator
write(m_link, 'linked_emulator')
m_link <- read('linked_emulator')

## End(Not run)

Locate the next design point for a (D)GP emulator or a bundle of (D)GP emulators using MICE

Description

This function searches from a candidate set to locate the next design point(s) to be added to a (D)GP emulator or a bundle of (D)GP emulators using the Mutual Information for Computer Experiments (MICE), see the reference below.

Usage

mice(object, ...)

## S3 method for class 'gp'
mice(
  object,
  x_cand = NULL,
  n_cand = 200,
  batch_size = 1,
  M = 50,
  nugget_s = 1e-06,
  workers = 1,
  limits = NULL,
  int = FALSE,
  ...
)

## S3 method for class 'dgp'
mice(
  object,
  x_cand = NULL,
  n_cand = 200,
  batch_size = 1,
  M = 50,
  nugget_s = 1e-06,
  workers = 1,
  limits = NULL,
  int = FALSE,
  aggregate = NULL,
  ...
)

## S3 method for class 'bundle'
mice(
  object,
  x_cand = NULL,
  n_cand = 200,
  batch_size = 1,
  M = 50,
  nugget_s = 1e-06,
  workers = 1,
  limits = NULL,
  int = FALSE,
  aggregate = NULL,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

1. If x_cand is not NULL:

   • When object is an instance of the gp class, a vector of length batch_size is returned, containing the positions (row numbers) of the next design points from x_cand.

   • When object is an instance of the dgp class, a vector of length batch_size * D is returned, containing the positions (row numbers) of the next design points from x_cand to be added to the DGP emulator.

     • D is the number of output dimensions of the DGP emulator if no likelihood layer is included.

     • For a DGP emulator with a Hetero or NegBin likelihood layer, D = 2.

     • For a DGP emulator with a Categorical likelihood layer, D = 1 for binary output or D = K for multi-class output with K classes.

   • When object is an instance of the bundle class, a matrix is returned with batch_size rows and a column for each emulator in the bundle, containing the positions (row numbers) of the next design points from x_cand for individual emulators.

2. If x_cand is NULL:

   • When object is an instance of the gp class, a matrix with batch_size rows is returned, giving the next design points to be evaluated.

   • When object is an instance of the dgp class, a matrix with batch_size * D rows is returned, where:

     • D is the number of output dimensions of the DGP emulator if no likelihood layer is included.

     • For a DGP emulator with a Hetero or NegBin likelihood layer, D = 2.

     • For a DGP emulator with a Categorical likelihood layer, D = 1 for binary output or D = K for multi-class output with K classes.

   • When object is an instance of the bundle class, a list is returned with a length equal to the number of emulators in the bundle. Each element of the list is a matrix with batch_size rows, where each row represents a design point to be added to the corresponding emulator.

Note

The first column of the matrix supplied to the first argument of aggregate must correspond to the first output dimension of the DGP emulator if object is an instance of the dgp class, and so on for subsequent columns and dimensions. If object is an instance of the bundle class, the first column must correspond to the first emulator in the bundle, and so on for subsequent columns and emulators.

References

Beck, J., & Guillas, S. (2016). Sequential design with mutual information for computer experiments (MICE): emulation of a tsunami model. SIAM/ASA Journal on Uncertainty Quantification, 4(1), 739-766.

Examples

## Not run: 

# load packages and the Python env
library(lhs)
library(dgpsi)

# construct a 1D non-stationary function
f <- function(x) {
 sin(30*((2*x-1)/2-0.4)^5)*cos(20*((2*x-1)/2-0.4))
}

# generate the initial design
X <- maximinLHS(10,1)
Y <- f(X)

# training a 2-layered DGP emulator with the global connection off
m <- dgp(X, Y, connect = F)

# generate a candidate set
x_cand <- maximinLHS(200,1)

# locate the next design point using MICE
next_point <- mice(m, x_cand = x_cand)
X_new <- x_cand[next_point,,drop = F]

# obtain the corresponding output at the located design point
Y_new <- f(X_new)

# combine the new input-output pair to the existing data
X <- rbind(X, X_new)
Y <- rbind(Y, Y_new)

# update the DGP emulator with the new input and output data and refit
m <- update(m, X, Y, refit = TRUE)

# plot the LOO validation
plot(m)

## End(Not run)

Calculate the predictive negative log-likelihood

Description

This function computes the predictive negative log-likelihood from a DGP emulator with a likelihood layer.

Usage

nllik(object, x, y)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object is returned with an additional slot named NLL that contains two elements. The first one, named meanNLL, is a scalar that gives the average negative predicted log-likelihood across all testing data points. The second one, named allNLL, is a vector that gives the negative predicted log-likelihood for each testing data point.

Pack GP and DGP emulators into a bundle

Description

This function packs GP emulators and DGP emulators into a bundle class for sequential designs if each emulator emulates one output dimension of the underlying simulator.

Usage

pack(..., id = NULL)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An S3 class named bundle to be used by design() for sequential designs. It has:

• a slot called id that is assigned through the id argument.

• N slots named ⁠emulator1,...,emulatorN⁠, each of which contains a GP or DGP emulator, where N is the number of emulators that are provided to the function.

• a slot called data which contains two elements X and Y. X contains N matrices named ⁠emulator1,...,emulatorN⁠ that are training input data for different emulators. Y contains N single-column matrices named ⁠emulator1,...,emulatorN⁠ that are training output data for different emulators.

Examples

## Not run: 

# load packages
library(lhs)
library(dgpsi)

# construct a function with a two-dimensional output
f <- function(x) {
 y1 = sin(30*((2*x-1)/2-0.4)^5)*cos(20*((2*x-1)/2-0.4))
 y2 = 1/3*sin(2*(2*x - 1))+2/3*exp(-30*(2*(2*x-1))^2)+1/3
 return(cbind(y1,y2))
}

# generate the initial design
X <- maximinLHS(10,1)
Y <- f(X)

# generate the validation data
validate_x <- maximinLHS(30,1)
validate_y <- f(validate_x)

# training a 2-layered DGP emulator with respect to each output with the global connection off
m1 <- dgp(X, Y[,1], connect = F)
m2 <- dgp(X, Y[,2], connect = F)

# specify the range of the input dimension
lim <- c(0, 1)

# pack emulators to form an emulator bundle
m <- pack(m1, m2)

# 1st wave of the sequential design with 10 iterations and the target RMSE of 0.01
m <- design(m, N = 10, limits = lim, f = f, x_test = validate_x, y_test = validate_y, target = 0.01)

# 2rd wave of the sequential design with additional 10 iterations and the same target
m <- design(m, N = 10, limits = lim, f = f, x_test = validate_x, y_test = validate_y, target = 0.01)

# draw sequential designs of the two packed emulators
draw(m, type = 'design')

# inspect the traces of RMSEs of the two packed emulators
draw(m, type = 'rmse')

# write and read the constructed emulator bundle
write(m, 'bundle_dgp')
m <- read('bundle_dgp')

# unpack the bundle into individual emulators
m_unpacked <- unpack(m)

# plot OOS validations of individual emulators
plot(m_unpacked[[1]], x_test = validate_x, y_test = validate_y[,1])
plot(m_unpacked[[2]], x_test = validate_x, y_test = validate_y[,2])

## End(Not run)

Validation plots of a constructed GP, DGP, or linked (D)GP emulator

Description

This function draws validation plots of a GP, DGP, or linked (D)GP emulator.

Usage

## S3 method for class 'dgp'
plot(
  x,
  x_test = NULL,
  y_test = NULL,
  dim = NULL,
  method = NULL,
  sample_size = 50,
  style = 1,
  min_max = TRUE,
  normalize = TRUE,
  color = "turbo",
  type = "points",
  verb = TRUE,
  M = 50,
  force = FALSE,
  cores = 1,
  ...
)

## S3 method for class 'lgp'
plot(
  x,
  x_test = NULL,
  y_test = NULL,
  dim = NULL,
  method = NULL,
  sample_size = 50,
  style = 1,
  min_max = TRUE,
  color = "turbo",
  type = "points",
  M = 50,
  verb = TRUE,
  force = FALSE,
  cores = 1,
  ...
)

## S3 method for class 'gp'
plot(
  x,
  x_test = NULL,
  y_test = NULL,
  dim = NULL,
  method = NULL,
  sample_size = 50,
  style = 1,
  min_max = TRUE,
  color = "turbo",
  type = "points",
  verb = TRUE,
  M = 50,
  force = FALSE,
  cores = 1,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

A patchwork object.

Note

• plot() calls validate() internally to obtain validation results for plotting. However, plot() will not export the emulator object with validation results. Instead, it only returns the plotting object. For small-scale validations (i.e., small training or testing data points), direct execution of plot() works well. However, for moderate- to large-scale validation, it is recommended to first run validate() to obtain and store validation results in the emulator object, and then supply the object to plot(). plot() checks the object's loo and oos slots prior to calling validate() and will not perform further calculation if the required information is already stored.

• plot() will only use stored OOS validation if x_test and y_test are identical to those used by validate() to produce the data contained in the object's oos slot, otherwise plot() will re-evaluate OOS validation before plotting.

• The returned patchwork::patchwork object contains the ggplot2::ggplot2 objects. One can modify the included individual ggplots by accessing them with double-bracket indexing. See https://patchwork.data-imaginist.com/ for further information.

Examples

## Not run: 

# See gp(), dgp(), or lgp() for an example.

## End(Not run)

Prediction from GP, DGP, or linked (D)GP emulators

Description

This function implements prediction from GP, DGP, or linked (D)GP emulators.

Usage

## S3 method for class 'dgp'
predict(
  object,
  x,
  method = NULL,
  mode = "label",
  full_layer = FALSE,
  sample_size = 50,
  M = 50,
  cores = 1,
  chunks = NULL,
  ...
)

## S3 method for class 'lgp'
predict(
  object,
  x,
  method = NULL,
  full_layer = FALSE,
  sample_size = 50,
  M = 50,
  cores = 1,
  chunks = NULL,
  ...
)

## S3 method for class 'gp'
predict(
  object,
  x,
  method = NULL,
  sample_size = 50,
  M = 50,
  cores = 1,
  chunks = NULL,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

• If object is an instance of the gp class:

  1. if method = "mean_var": an updated object is returned with an additional slot called results that contains two matrices named mean for the predictive means and var for the predictive variances. Each matrix has only one column with its rows corresponding to testing positions (i.e., rows of x).

  2. if method = "sampling": an updated object is returned with an additional slot called results that contains a matrix whose rows correspond to testing positions and columns correspond to sample_size number of samples drawn from the predictive distribution of GP.

• If object is an instance of the dgp class:

  1. if method = "mean_var" and full_layer = FALSE: an updated object is returned with an additional slot called results that contains two matrices named mean for the predictive means and var for the predictive variances respectively. Each matrix has its rows corresponding to testing positions and columns corresponding to DGP global output dimensions (i.e., the number of GP/likelihood nodes in the final layer).

  2. if method = "mean_var" and full_layer = TRUE: an updated object is returned with an additional slot called results that contains two sub-lists named mean for the predictive means and var for the predictive variances respectively. Each sub-list contains L (i.e., the number of layers) matrices named ⁠layer1, layer2,..., layerL⁠. Each matrix has its rows corresponding to testing positions and columns corresponding to output dimensions (i.e., the number of GP/likelihood nodes from the associated layer).

  3. if method = "sampling" and full_layer = FALSE: an updated object is returned with an additional slot called results that contains D (i.e., the number of GP/likelihood nodes in the final layer) matrices named ⁠output1, output2,..., outputD⁠. Each matrix has its rows corresponding to testing positions and columns corresponding to samples of size: B * sample_size, where B is the number of imputations specified in dgp().

  4. if method = "sampling" and full_layer = TRUE: an updated object is returned with an additional slot called results that contains L (i.e., the number of layers) sub-lists named ⁠layer1, layer2,..., layerL⁠. Each sub-list represents samples drawn from the GP/likelihood nodes in the corresponding layer, and contains D (i.e., the number of GP/likelihood nodes in the corresponding layer) matrices named ⁠output1, output2,..., outputD⁠. Each matrix gives samples of the output from one of D GP/likelihood nodes, and has its rows corresponding to testing positions and columns corresponding to samples of size: B * sample_size, where B is the number of imputations specified in dgp().

• If object is an instance of the dgp class with a categorical likelihood:

  1. if full_layer = FALSE and mode = "label": an updated object is returned with an additional slot called results that contains one matrix named label. The matrix has rows corresponding to testing positions and columns corresponding to sample labels of size: B * sample_size, where B is the number of imputations specified in dgp().

  2. if full_layer = FALSE and mode = "proba", an updated object is returned with an additional slot called results. This slot contains D matrices (where D is the number of classes in the training output), where each matrix gives probability samples for the corresponding class with its rows corresponding to testing positions and columns containing probabilities. The number of columns of each matrix is B * sample_size, where B is the number of imputations specified in the dgp() function.

  3. if method = "mean_var" and full_layer = TRUE: an updated object is returned with an additional slot called results that contains L (i.e., the number of layers) sub-lists named ⁠layer1, layer2,..., layerL⁠. Each of first L-1 sub-lists contains two matrices named mean for the predictive means and var for the predictive variances of the GP nodes in the associated layer. Rows of each matrix correspond to testing positions.

     • when mode = "label", the sub-list LayerL contains one matrix named label. The matrix has its rows corresponding to testing positions and columns corresponding to label samples of size: B * sample_size. B is the number of imputations specified in dgp().

     • when mode = "proba", the sub-list LayerL contains D matrices (where D is the number of classes in the training output), where each matrix gives probability samples for the corresponding class with its rows corresponding to testing positions and columns containing probabilities. The number of columns of each matrix is B * sample_size. B is the number of imputations specified in dgp().

  4. if method = "sampling" and full_layer = TRUE: an updated object is returned with an additional slot called results that contains L (i.e., the number of layers) sub-lists named ⁠layer1, layer2,..., layerL⁠. Each of first L-1 sub-lists represents samples drawn from the GP nodes in the corresponding layer, and contains D (i.e., the number of GP nodes in the corresponding layer) matrices named ⁠output1, output2,..., outputD⁠. Each matrix gives samples of the output from one of D GP nodes, and has its rows corresponding to testing positions and columns corresponding to samples of size: B * sample_size.

     • when mode = "label", the sub-list LayerL contains one matrix named label. The matrix has its rows corresponding to testing positions and columns corresponding to label samples of size: B * sample_size.

     • when mode = "proba", the sub-list LayerL contains D matrices (where D is the number of classes in the training output), where each matrix gives probability samples for the corresponding class with its rows corresponding to testing positions and columns containing probabilities. The number of columns of each matrix is B * sample_size.

     B is the number of imputations specified in dgp().

• If object is an instance of the lgp class:

  1. if method = "mean_var" and full_layer = FALSE: an updated object is returned with an additional slot called results that contains two sub-lists named mean for the predictive means and var for the predictive variances respectively. Each sub-list contains K (same number of emulators in the final layer of the system) matrices named using the IDs of the corresponding emulators in the final layer. Each matrix has rows corresponding to global testing positions and columns corresponding to output dimensions of the associated emulator in the final layer.

  2. if method = "mean_var" and full_layer = TRUE: an updated object is returned with an additional slot called results that contains two sub-lists named mean for the predictive means and var for the predictive variances respectively. Each sub-list contains L (i.e., the number of layers in the emulated system) components named ⁠layer1, layer2,..., layerL⁠. Each component represents a layer and contains K (same number of emulators in the corresponding layer of the system) matrices named using the IDs of the corresponding emulators in that layer. Each matrix has its rows corresponding to global testing positions and columns corresponding to output dimensions of the associated GP/DGP emulator in the corresponding layer.

  3. if method = "sampling" and full_layer = FALSE: an updated object is returned with an additional slot called results that contains K (same number of emulators in the final layer of the system) sub-lists named using the IDs of the corresponding emulators in the final layer. Each sub-list contains D matrices, named ⁠output1, output2,..., outputD⁠, that correspond to the output dimensions of the GP/DGP emulator. Each matrix has rows corresponding to testing positions and columns corresponding to samples of size: B * sample_size, where B is the number of imputations specified in lgp().

  4. if method = "sampling" and full_layer = TRUE: an updated object is returned with an additional slot called results that contains L (i.e., the number of layers of the emulated system) sub-lists named ⁠layer1, layer2,..., layerL⁠. Each sub-list represents a layer and contains K (same number of emulators in the corresponding layer of the system) components named using the IDs of the corresponding emulators in that layer. Each component contains D matrices, named ⁠output1, output2,..., outputD⁠, that correspond to the output dimensions of the GP/DGP emulator. Each matrix has its rows corresponding to testing positions and columns corresponding to samples of size: B * sample_size, where B is the number of imputations specified in lgp().

  If object is an instance of the lgp class created by lgp() without specifying the struc argument in data frame form, the IDs, that are used as names of sub-lists or matrices within results, will be replaced by emulator1, emulator2, and so on.

The results slot will also include:

• the value of M, which represents the size of the conditioning set for the Vecchia approximation, if used, in the emulator prediction.

• the value of sample_size if method = "sampling".

Examples

## Not run: 

# See gp(), dgp(), or lgp() for an example.

## End(Not run)

Static pruning of a DGP emulator

Description

This function implements static pruning for a DGP emulator.

Usage

prune(object, control = list(), verb = TRUE)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object that could be an instance of gp, dgp, or bundle (of GP emulators) class.

Note

• The function requires a DGP emulator that has been trained with a dataset comprising a minimum size equal to min_size in control. If the training dataset size is smaller than this, it is recommended that the design of the DGP emulator is enriched and its structure pruned dynamically using the design() function. Depending on the design of the DGP emulator, static pruning may not be accurate. It is thus recommended that dynamic pruning is implemented as a part of a sequential design via design().

• The following slots:

  • loo and oos created by validate(); and

  • results created by predict();

  in object will be removed and not contained in the returned object.

Examples

## Not run: 

# load the package and the Python env
library(dgpsi)

# construct the borehole function over a hypercube
f <- function(x){
  x[,1] <- (0.15 - 0.5) * x[,1] + 0.5
  x[,2] <- exp((log(50000) - log(100)) * x[,2] + log(100))
  x[,3] <- (115600 - 63070) *x[,3] + 63070
  x[,4] <- (1110 - 990) * x[,4] + 990
  x[,5] <- (116 - 63.1) * x[,5] + 63.1
  x[,6] <- (820 - 700) * x[,6] + 700
  x[,7] <- (1680 - 1120) * x[,7] + 1120
  x[,8] <- (12045 - 9855) * x[,8] + 9855
  y <- apply(x, 1, RobustGaSP::borehole)
}

# set a random seed
set_seed(999)

# generate training data
X <- maximinLHS(80, 8)
Y <- f(X)

# generate validation data
validate_x <- maximinLHS(500, 8)
validate_y <- f(validate_x)

# training a DGP emulator with anisotropic squared exponential kernels
m <- dgp(X, Y, share = F)

# OOS validation of the DGP emulator
plot(m, validate_x, validate_y)

# prune the emulator until no more GP nodes are removable
m <- prune(m)

# OOS validation of the resulting emulator
plot(m, validate_x, validate_y)

## End(Not run)

Load the stored emulator

Description

This function loads the .pkl file that stores the emulator.

Usage

read(pkl_file)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

The S3 class of a GP emulator, a DGP emulator, a linked (D)GP emulator, or a bundle of (D)GP emulators.

Examples

## Not run: 

# See gp(), dgp(), lgp(), or pack() for an example.

## End(Not run)

Serialize the constructed emulator

Description

This function serializes the constructed emulator.

Usage

serialize(object, light = TRUE)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

A serialized version of object.

Note

Since the constructed emulators are 'python' objects, they cannot be directly exported to other R processes for parallel processing in multi-session workers created through spawning. This function provides a solution by converting the emulators into serialized objects, which can be restored using deserialize() for multi-session processing. Note that in forking, serialization is generally not required.

Examples

## Not run: 

library(future)
library(future.apply)
library(dgpsi)

# model
f <- function(x) {
 (sin(7.5*x)+1)/2
}

# training data
X <- seq(0, 1, length = 10)
Y <- sapply(X, f)

# train a DGP emulator
m <- dgp(X, Y, name = "matern2.5")

# testing input data
X_dgp <- seq(0, 1, length = 100)

# serialize the DGP emulator
m_serialized <- serialize(m)

# start a multi-session with three cores for parallel predictions
plan(multisession, workers = 3)

# perform parallel predictions
results <- future_lapply(1:length(X_dgp), function(i) {
  m_deserialized <- deserialize(m_serialized)
  mean_i <- predict(m_deserialized, X_dgp[i])$results$mean
}, future.seed = TRUE)

# reset the future plan to sequential
plan(sequential)

# combine mean predictions
pred_mean <- do.call(rbind, results)

## End(Not run)

Set Emulator ID

Description

This function assigns a unique identifier to an emulator.

Usage

set_id(object, id)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

The updated object, with the assigned ID stored in its id slot.

Examples

## Not run: 

# See lgp() for an example.

## End(Not run)

Reset number of imputations for a DGP emulator

Description

This function resets the number of imputations for prediction from a DGP emulator.

Usage

set_imp(object, B = 5)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object with the information of B incorporated.

Note

• This function is useful when a DGP emulator has been trained and one wants to make faster predictions by decreasing the number of imputations without rebuilding the emulator.

• The following slots:

  • loo and oos created by validate(); and

  • results created by predict() in object will be removed and not contained in the returned object.

Examples

## Not run: 

# See design() for an example.

## End(Not run)

Set linked indices

Description

This function is deprecated and will be removed in the next release. The updated lgp() function now offers a simpler, more efficient way to specify linked information for (D)GP emulators.

Usage

set_linked_idx(object, idx)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object with the information of idx incorporated.

Note

This function is useful when different models are emulated by different teams. Each team can create their (D)GP emulator even without knowing how different emulators are connected together. When this information is available and different emulators are collected, the connection information between emulators can then be assigned to individual emulators with this function.

Examples

## Not run: 

# See lgp() for an example.

## End(Not run)

Random seed generator

Description

This function initializes a random number generator that sets the random seed in both R and Python to ensure reproducible results from the package.

Usage

set_seed(seed)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

No return value.

Examples

## Not run: 

# See dgp() for an example.

## End(Not run)

Set the number of threads

Description

This function sets the number of threads for parallel computations involved in the package.

Usage

set_thread_num(num)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

No return value.

Add or remove the Vecchia approximation

Description

This function adds or removes the Vecchia approximation from a GP, DGP or linked (D)GP emulator constructed by gp(), dgp() or lgp().

Usage

set_vecchia(object, vecchia = TRUE, M = 25, ord = NULL)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object with the Vecchia approximation either added or removed.

Note

This function is useful for quickly switching between Vecchia and non-Vecchia approximations for an existing emulator without the need to reconstruct the emulator. If the emulator was built without the Vecchia approximation, the function can add it, and if the emulator was built with the Vecchia approximation, the function can remove it. If the current state already matches the requested state, the emulator remains unchanged.

Summary of a constructed GP, DGP, or linked (D)GP emulator

Description

This function provides a summary of key information for a GP, DGP, or linked (D)GP emulator by generating either a table or an interactive plot of the emulator’s structure.

Usage

## S3 method for class 'gp'
summary(object, type = "plot", ...)

## S3 method for class 'dgp'
summary(object, type = "plot", ...)

## S3 method for class 'lgp'
summary(object, type = "plot", group_size = 1, ...)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

Either a summary table (returned as kableExtra object) or an interactive visualization (returned as a visNetwork object) of the emulator. The visualization is compatible with R Markdown documents and the RStudio Viewer. The summary table can be further customized by kableExtra::kableExtra package. The resulting visNetwork object can be saved as an HTML file using visNetwork::visSave() from the visNetwork::visNetwork package.

Examples

## Not run: 

# See gp(), dgp(), or lgp() for an example.

## End(Not run)

Trace plot for DGP hyperparameters

Description

This function draws trace plots for the hyperparameters of a chosen GP node in a DGP emulator.

Usage

trace_plot(object, layer = NULL, node = 1)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

A ggplot object.

Examples

## Not run: 

# See dgp() for an example.

## End(Not run)

Unpack a bundle of (D)GP emulators

Description

This function unpacks a bundle of (D)GP emulators safely so that any further manipulations of unpacked individual emulators will not impact those in the bundle.

Usage

unpack(object)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

A named list that contains individual emulators (named ⁠emulator1,...,emulatorS⁠) packed in object, where S is the number of emulators in object.

Examples

## Not run: 

# See pack() for an example.

## End(Not run)

Update a GP or DGP emulator

Description

This function updates the training input and output of a GP or DGP emulator with an option to refit the emulator.

Usage

update(object, X, Y, refit, reset, verb, ...)

## S3 method for class 'dgp'
update(
  object,
  X,
  Y,
  refit = TRUE,
  reset = FALSE,
  verb = TRUE,
  N = NULL,
  cores = 1,
  ess_burn = 10,
  B = NULL,
  ...
)

## S3 method for class 'gp'
update(object, X, Y, refit = TRUE, reset = FALSE, verb = TRUE, ...)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object.

Note

• The following slots:

  • loo and oos created by validate();

  • results created by predict(); and

  • design created by design()

  in object will be removed and not contained in the returned object.

Examples

## Not run: 

# See alm(), mice(), or vigf() for an example.

## End(Not run)

Validate a constructed GP, DGP, or linked (D)GP emulator

Description

This function calculates Leave-One-Out (LOO) cross validation or Out-Of-Sample (OOS) validation statistics for a constructed GP, DGP, or linked (D)GP emulator.

Usage

validate(
  object,
  x_test,
  y_test,
  method,
  sample_size,
  verb,
  M,
  force,
  cores,
  ...
)

## S3 method for class 'gp'
validate(
  object,
  x_test = NULL,
  y_test = NULL,
  method = NULL,
  sample_size = 50,
  verb = TRUE,
  M = 50,
  force = FALSE,
  cores = 1,
  ...
)

## S3 method for class 'dgp'
validate(
  object,
  x_test = NULL,
  y_test = NULL,
  method = NULL,
  sample_size = 50,
  verb = TRUE,
  M = 50,
  force = FALSE,
  cores = 1,
  ...
)

## S3 method for class 'lgp'
validate(
  object,
  x_test = NULL,
  y_test = NULL,
  method = NULL,
  sample_size = 50,
  verb = TRUE,
  M = 50,
  force = FALSE,
  cores = 1,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

• If object is an instance of the gp class, an updated object is returned with an additional slot called loo (for LOO cross validation) or oos (for OOS validation) that contains:

  • two slots called x_train (or x_test) and y_train (or y_test) that contain the validation data points for LOO (or OOS).

  • a column matrix called mean, if method = "mean_var", or median, if method = "sampling", that contains the predictive means or medians of the GP emulator at validation positions.

  • three column matrices called std, lower, and upper that contain the predictive standard deviations and credible intervals of the GP emulator at validation positions. If method = "mean_var", the upper and lower bounds of a credible interval are two standard deviations above and below the predictive mean. If method = "sampling", the upper and lower bounds of a credible interval are 2.5th and 97.5th percentiles.

  • a numeric value called rmse that contains the root mean/median squared error of the GP emulator.

  • a numeric value called nrmse that contains the (max-min) normalized root mean/median squared error of the GP emulator. The max-min normalization uses the maximum and minimum values of the validation outputs contained in y_train (or y_test).

  • an integer called M that contains the size of the conditioning set used for the Vecchia approximation, if used, for emulator validation.

  • an integer called sample_size that contains the number of samples used for validation if method = "sampling".

  The rows of matrices (mean, median, std, lower, and upper) correspond to the validation positions.

• If object is an instance of the dgp class, an updated object is returned with an additional slot called loo (for LOO cross validation) or oos (for OOS validation) that contains:

  • two slots called x_train (or x_test) and y_train (or y_test) that contain the validation data points for LOO (or OOS).

  • a matrix called mean, if method = "mean_var", or median, if method = "sampling", that contains the predictive means or medians of the DGP emulator at validation positions.

  • three matrices called std, lower, and upper that contain the predictive standard deviations and credible intervals of the DGP emulator at validation positions. If method = "mean_var", the upper and lower bounds of a credible interval are two standard deviations above and below the predictive mean. If method = "sampling", the upper and lower bounds of a credible interval are 2.5th and 97.5th percentiles.

  • a vector called rmse that contains the root mean/median squared errors of the DGP emulator across different output dimensions.

  • a vector called nrmse that contains the (max-min) normalized root mean/median squared errors of the DGP emulator across different output dimensions. The max-min normalization uses the maximum and minimum values of the validation outputs contained in y_train (or y_test).

  • an integer called M that contains size of the conditioning set used for the Vecchia approximation, if used, for emulator validation.

  • an integer called sample_size that contains the number of samples used for validation if method = "sampling".

  The rows and columns of matrices (mean, median, std, lower, and upper) correspond to the validation positions and DGP emulator output dimensions, respectively.

• If object is an instance of the dgp class with a categorical likelihood, an updated object is returned with an additional slot called loo (for LOO cross validation) or oos (for OOS validation) that contains:

  • two slots called x_train (or x_test) and y_train (or y_test) that contain the validation data points for LOO (or OOS).

  • a matrix called label that contains predictive samples of labels from the DGP emulator at validation positions. The matrix has its rows corresponding to validation positions and columns corresponding to samples of labels.

  • a list called probability that contains predictive samples of probabilities for each class from the DGP emulator at validation positions. The element in the list is a matrix that has its rows corresponding to validation positions and columns corresponding to samples of probabilities.

  • a scalar called log_loss that represents the average log loss of the predicted labels in the DGP emulator across all validation positions. Log loss measures the accuracy of probabilistic predictions, with lower values indicating better classification performance. log_loss ranges from 0 to positive infinity, where a value closer to 0 suggests more confident and accurate predictions.

  • an integer called M that contains size of the conditioning set used for the Vecchia approximation, if used, in emulator validation.

  • an integer called sample_size that contains the number of samples used for validation.

• If object is an instance of the lgp class, an updated object is returned with an additional slot called oos (for OOS validation) that contains:

  • two slots called x_test and y_test that contain the validation data points for OOS.

  • a list called mean, if method = "mean_var", or median, if method = "sampling", that contains the predictive means or medians of the linked (D)GP emulator at validation positions.

  • three lists called std, lower, and upper that contain the predictive standard deviations and credible intervals of the linked (D)GP emulator at validation positions. If method = "mean_var", the upper and lower bounds of a credible interval are two standard deviations above and below the predictive mean. If method = "sampling", the upper and lower bounds of a credible interval are 2.5th and 97.5th percentiles.

  • a list called rmse that contains the root mean/median squared errors of the linked (D)GP emulator.

  • a list called nrmse that contains the (max-min) normalized root mean/median squared errors of the linked (D)GP emulator. The max-min normalization uses the maximum and minimum values of the validation outputs contained in y_test.

  • an integer called M that contains size of the conditioning set used for the Vecchia approximation, if used, in emulator validation.

  • an integer called sample_size that contains the number of samples used for validation if method = "sampling".

  Each element in mean, median, std, lower, upper, rmse, and nrmse corresponds to a (D)GP emulator in the final layer of the linked (D)GP emulator.

Note

• When both x_test and y_test are NULL, LOO cross validation will be implemented. Otherwise, OOS validation will be implemented. LOO validation is only applicable to a GP or DGP emulator (i.e., object is an instance of the gp or dgp class). If a linked (D)GP emulator (i.e., object is an instance of the lgp class) is provided, x_test and y_test must also be provided for OOS validation.

Examples

## Not run: 

# See gp(), dgp(), or lgp() for an example.

## End(Not run)

Locate the next design point for a (D)GP emulator or a bundle of (D)GP emulators using VIGF

Description

This function searches from a candidate set to locate the next design point(s) to be added to a (D)GP emulator or a bundle of (D)GP emulators using the Variance of Improvement for Global Fit (VIGF). For VIGF on GP emulators, see the reference below.

Usage

vigf(object, ...)

## S3 method for class 'gp'
vigf(
  object,
  x_cand = NULL,
  n_start = 10,
  batch_size = 1,
  M = 50,
  workers = 1,
  limits = NULL,
  int = FALSE,
  ...
)

## S3 method for class 'dgp'
vigf(
  object,
  x_cand = NULL,
  n_start = 10,
  batch_size = 1,
  M = 50,
  workers = 1,
  limits = NULL,
  int = FALSE,
  aggregate = NULL,
  ...
)

## S3 method for class 'bundle'
vigf(
  object,
  x_cand = NULL,
  n_start = 10,
  batch_size = 1,
  M = 50,
  workers = 1,
  limits = NULL,
  int = FALSE,
  aggregate = NULL,
  ...
)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

1. If x_cand is not NULL:

   • When object is an instance of the gp class, a vector of length batch_size is returned, containing the positions (row numbers) of the next design points from x_cand.

   • When object is an instance of the dgp class, a vector of length batch_size * D is returned, containing the positions (row numbers) of the next design points from x_cand to be added to the DGP emulator.

     • D is the number of output dimensions of the DGP emulator if no likelihood layer is included.

     • For a DGP emulator with a Hetero or NegBin likelihood layer, D = 2.

     • For a DGP emulator with a Categorical likelihood layer, D = 1 for binary output or D = K for multi-class output with K classes.

   • When object is an instance of the bundle class, a matrix is returned with batch_size rows and a column for each emulator in the bundle, containing the positions (row numbers) of the next design points from x_cand for individual emulators.

2. If x_cand is NULL:

   • When object is an instance of the gp class, a matrix with batch_size rows is returned, giving the next design points to be evaluated.

   • When object is an instance of the dgp class, a matrix with batch_size * D rows is returned, where:

     • D is the number of output dimensions of the DGP emulator if no likelihood layer is included.

     • For a DGP emulator with a Hetero or NegBin likelihood layer, D = 2.

     • For a DGP emulator with a Categorical likelihood layer, D = 1 for binary output or D = K for multi-class output with K classes.

   • When object is an instance of the bundle class, a list is returned with a length equal to the number of emulators in the bundle. Each element of the list is a matrix with batch_size rows, where each row represents a design point to be added to the corresponding emulator.

Note

The first column of the matrix supplied to the first argument of aggregate must correspond to the first output dimension of the DGP emulator if object is an instance of the dgp class, and so on for subsequent columns and dimensions. If object is an instance of the bundle class, the first column must correspond to the first emulator in the bundle, and so on for subsequent columns and emulators.

References

Mohammadi, H., & Challenor, P. (2022). Sequential adaptive design for emulating costly computer codes. arXiv:2206.12113.

Examples

## Not run: 

# load packages and the Python env
library(lhs)
library(dgpsi)

# construct a 1D non-stationary function
f <- function(x) {
 sin(30*((2*x-1)/2-0.4)^5)*cos(20*((2*x-1)/2-0.4))
}

# generate the initial design
X <- maximinLHS(10,1)
Y <- f(X)

# training a 2-layered DGP emulator with the global connection off
m <- dgp(X, Y, connect = F)

# specify the input range
lim <- c(0,1)

# locate the next design point using VIGF
X_new <- vigf(m, limits = lim)

# obtain the corresponding output at the located design point
Y_new <- f(X_new)

# combine the new input-output pair to the existing data
X <- rbind(X, X_new)
Y <- rbind(Y, Y_new)

# update the DGP emulator with the new input and output data and refit
m <- update(m, X, Y, refit = TRUE)

# plot the LOO validation
plot(m)

## End(Not run)

Trim the sequence of hyperparameter estimates within a DGP emulator

Description

This function trims the sequence of hyperparameter estimates within a DGP emulator generated during training.

Usage

window(object, start, end = NULL, thin = 1)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An updated object with a trimmed sequence of hyperparameters.

Note

• This function is useful when a DGP emulator has been trained and one wants to trim the sequence of hyperparameters estimated and to use the trimmed sequence to generate point estimates of the DGP model parameters for prediction.

• The following slots:

  • loo and oos created by validate(); and

  • results created by predict() in object will be removed and not contained in the returned object.

Examples

## Not run: 

# See dgp() for an example.

## End(Not run)

Save the constructed emulator

Description

This function saves the constructed emulator to a .pkl file.

Usage

write(object, pkl_file, light = TRUE)

Arguments

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

No return value. object will be saved to a local .pkl file specified by pkl_file.

Note

Since emulators built from the package are 'python' objects, save() from R will not work as it would for R objects. If object was processed by set_vecchia() to add or remove the Vecchia approximation, light should be set to FALSE to ensure reproducibility after the saved emulator is reloaded by read().

Examples

## Not run: 

# See gp(), dgp(), lgp(), or pack() for an example.

## End(Not run)