Package {magp}

magp: Mapping-Based Additive Gaussian Process Models

Description

Fits additive Gaussian process models for experiments in which each component has both a quantitative level and a position in a sequence.

Details

For q components, the first q input columns contain quantitative levels and the next q columns contain sequence positions. Every sequence row must be a permutation of 1:q. Quantitative columns outside ⁠[0, 1]⁠ are scaled during fitting, and the same transformation is used for new data.

magp2d_fit() represents the sequence positions in two latent dimensions. magpfull_fit() uses q - 1 latent dimensions. Both fitted model classes support stats::predict() for point predictions and plug-in predictive uncertainty. The computational kernels for covariance matrices, analytical gradients, and cross-covariances are implemented in C++ with Rcpp.

magp_initial_design() constructs a space-filling quantitative-sequence design before responses are collected. It combines a Latin hypercube with sequence permutations generated randomly or improved with simulated annealing or space-filling threshold accepting. Joint alignment preserves both component designs.

magp_expected_improvement() evaluates improvement using latent predictive uncertainty. magp_next_point() searches quantitative bounds and sequence permutations for the next experiment, while magp_bayes_optimize() runs the sequential fitting and evaluation loop from completed experiments. magp_bayes_optimize_from_scratch() first generates and evaluates an initial design, then continues through the same sequential loop.

Author(s)

Maintainer: Tony Wang wangtony883@gmail.com [copyright holder]

Authors:

• Tony Wang wangtony883@gmail.com [copyright holder]

• Qian Xiao [copyright holder]

Other contributors:

• Yaping Wang [copyright holder]

• Abhyuday Mandal [copyright holder]

• Xinwei Deng [copyright holder]

See Also

Useful links:

• https://github.com/twskr/magp

• https://CRAN.R-project.org/package=magp

• Report bugs at https://github.com/twskr/magp/issues

Fit a MaGP model with a two-dimensional sequence map

Description

Fits an additive Gaussian process for data that combine component amounts with a component order. Sequence positions are represented by points in a compact two-dimensional latent map.

Usage

magp2d_fit(
  X,
  y = NULL,
  q = NULL,
  tau = 0.001,
  maxeval = 500,
  xtol_rel = 1e-05,
  lb_sigma = 10,
  ub_sigma = 1000,
  lb_theta = 0.5,
  ub_theta = 1000,
  lb_delta = -1,
  ub_delta = 1,
  seed = NULL,
  n_starts = 1,
  workers = 1
)

Arguments

Details

The covariance is a sum of component-specific terms. Each term combines the distance between two quantitative levels with the distance between their mapped sequence positions. The variance, correlation, and mapping parameters are estimated together with bounded optimization.

Quantitative columns outside ⁠[0, 1]⁠ are transformed with column-wise min-max scaling. Their training ranges are stored in the fitted object and reused for prediction. Columns already in ⁠[0, 1]⁠ are left unchanged.

When several starts are requested, each start receives a separate seed. Supplying seed makes the starts and the selected fit reproducible for both sequential and parallel execution. Start-level objective values, convergence codes, warnings, and errors are stored in fit$multistart$starts.

Value

An object of class magp2d.

Examples

train <- read.table(
  system.file("extdata", "example_train.txt", package = "magp"),
  header = TRUE
)
fit <- magp2d_fit(train, seed = 1, n_starts = 2)
fit

Calculate root-mean-squared prediction error

Description

Calculates the root-mean-squared error between predicted and observed values.

Usage

magp2d_rmse(pred, actual)

Arguments

Value

A single nonnegative numeric value.

Continue Bayesian optimization from completed experiments

Description

Use this function when initial experiments and their responses are already available. It fits a MaGP model, selects a new input with expected improvement, evaluates FUN, and adds the new result to the data. This process repeats for at most n_iter new evaluations.

Usage

magp_bayes_optimize(
  FUN,
  X,
  y = NULL,
  model = c("2d", "full"),
  direction = c("minimize", "maximize"),
  n_iter = 10L,
  xi = 0,
  stop_ei = 0,
  stop_patience = 3L,
  seed = NULL,
  fit_control = list(),
  acquisition_control = list(),
  objective_args = list(),
  verbose = TRUE
)

Arguments

Value

An object of class magp_bayes_opt. The final fitted model includes every completed evaluation.

How the loop works

The initial rows of X are treated as completed experiments and are not evaluated again. Each iteration performs four steps:

1. fit the selected MaGP model to all results collected so far;

2. use magp_next_point() to select an unobserved input;

3. call FUN at that input; and

4. add the response and refit the model.

FUN receives one named argument for each input column. For columns A, B, a, and b, for example, the call is equivalent to FUN(A = ..., B = ..., a = ..., b = ...).

Reading the result

The returned object contains:

• call: the function call;

• best_point: the input row with the best observed response;

• best_value: the best observed response;

• best_index: the row number of the best result in X and y;

• history: the initial and newly evaluated rows in evaluation order;

• model: the final fitted MaGP model;

• X and y: all inputs and responses used by the final model;

• acquisitions: details from each call to magp_next_point();

• mapping and direction: the model and optimization direction;

• iterations_requested and iterations_completed: the requested and completed numbers of new evaluations;

• stop_reason: why the loop ended;

• xi, stop_ei, and stop_patience: the acquisition and stopping settings; and

• fit_control and acquisition_control: the control lists used in the search.

References

Jones, D. R., Schonlau, M., and Welch, W. J. (1998). Efficient Global Optimization of Expensive Black-Box Functions. Journal of Global Optimization, 13, 455-492. doi:10.1023/A:1008306431147.

Examples

design <- magp_initial_design(n = 6, q = 3, seed = 1)
objective <- function(quantity_1, quantity_2, quantity_3,
                      sequence_1, sequence_2, sequence_3) {
  quantities <- c(quantity_1, quantity_2, quantity_3)
  positions <- c(sequence_1, sequence_2, sequence_3)
  -sum((quantities - c(0.2, 0.6, 0.8))^2) -
    0.02 * sum((positions - c(1, 3, 2))^2)
}
initial_y <- apply(design$design, 1L, function(row) {
  do.call(objective, as.list(row))
})
result <- magp_bayes_optimize(
  FUN = objective,
  X = design$design,
  y = initial_y,
  direction = "maximize",
  n_iter = 1,
  seed = 2,
  fit_control = list(maxeval = 100),
  acquisition_control = list(n_starts = 2, maxit = 20),
  verbose = FALSE
)
result$best_point
result$best_value
result$history

Start Bayesian optimization before any experiments have been run

Description

Use this function when no initial data are available. It creates an initial quantitative-sequence design, evaluates FUN at those runs, and then uses MaGP expected improvement to select later runs.

Usage

magp_bayes_optimize_from_scratch(
  FUN,
  n_initial,
  q,
  model = c("2d", "full"),
  direction = c("minimize", "maximize"),
  n_iter = 10L,
  xi = 0,
  stop_ei = 0,
  stop_patience = 3L,
  seed = NULL,
  design_control = list(),
  fit_control = list(),
  acquisition_control = list(),
  objective_args = list(),
  verbose = TRUE
)

Arguments

Value

An object of class magp_bayes_opt containing the initial design, all completed evaluations, and the final fitted model.

How to use this function

Define FUN, choose n_initial and q, and state whether the response should be minimized or maximized. The function then:

1. creates an initial design;

2. calls FUN once for each initial row;

3. fits the selected MaGP model; and

4. selects and evaluates as many as n_iter additional rows.

The generated columns are named quantity_1 through quantity_q and sequence_1 through sequence_q. Use magp_bayes_optimize() instead when initial experiments and responses already exist.

Reading the result

The result has the same main components as magp_bayes_optimize(), including best_point, best_value, history, the final model, and all recorded search settings. It also contains:

• initial_design: the generated quantitative-sequence design;

• initial_response: the response from each initial run;

• initial_evaluations: the number of initial runs; and

• design_control: the initial-design settings that were used; and

• started_from_initial_design: TRUE, marking that the workflow generated its own starting design.

Examples

objective <- function(quantity_1, quantity_2, quantity_3,
                      sequence_1, sequence_2, sequence_3) {
  quantities <- c(quantity_1, quantity_2, quantity_3)
  sequence <- c(sequence_1, sequence_2, sequence_3)
  -sum((quantities - c(0.2, 0.6, 0.8))^2) -
    0.01 * sum((sequence - c(1, 3, 2))^2)
}
result <- magp_bayes_optimize_from_scratch(
  FUN = objective,
  n_initial = 6,
  q = 3,
  direction = "maximize",
  n_iter = 1,
  seed = 1,
  design_control = list(
    sequence_maxit = 100,
    quantity_maxit = 100,
    alignment_maxit = 100
  ),
  fit_control = list(maxeval = 100),
  acquisition_control = list(n_starts = 2, maxit = 20),
  verbose = FALSE
)
result$initial_design$design
result$best_point
result$best_value
result$history

Score candidate experiments with expected improvement

Description

Calculates expected improvement for one or more candidate rows. Higher values indicate candidates that offer a better combination of predicted improvement and uncertainty.

Usage

magp_expected_improvement(
  object,
  newdata,
  direction = c("minimize", "maximize"),
  best = NULL,
  xi = 0
)

Arguments

Details

Expected improvement uses the model's latent predictive mean and standard error. A candidate can receive a high value because its predicted response is good, its uncertainty is large, or both. Rows already present in the training data receive a value of zero. Predictions with variance below 1e-8 also receive zero.

Value

A nonnegative numeric vector with one expected-improvement value per row of newdata.

References

Jones, D. R., Schonlau, M., and Welch, W. J. (1998). Efficient Global Optimization of Expensive Black-Box Functions. Journal of Global Optimization, 13, 455-492. doi:10.1023/A:1008306431147.

Examples

train <- read.table(
  system.file("extdata", "example_train.txt", package = "magp"),
  header = TRUE
)
test <- read.table(
  system.file("extdata", "example_test.txt", package = "magp"),
  header = TRUE
)
fit <- magp2d_fit(train, seed = 1)
magp_expected_improvement(
  fit,
  test[1:3, ],
  direction = "minimize"
)

Construct a quantitative-sequence initial design

Description

Builds an initial design in three steps. It first generates the sequence permutations, then constructs a maximin-style Latin hypercube for the quantitative levels, and finally aligns the two fixed designs by permuting whole quantitative rows. The final alignment preserves both the Latin hypercube and every sequence permutation.

Usage

magp_initial_design(
  n,
  q,
  pair_weight = 0.2,
  sequence_space_weight = 0.8,
  quantity_weight = 0.5,
  sequence_weight = 0.5,
  p = 15L,
  sequence_maxit = 10000L,
  quantity_maxit = 10000L,
  alignment_maxit = 10000L,
  sequence_temp = 0.1,
  quantity_temp = 0.01,
  alignment_temp = 0.001,
  tmax = 10L,
  initial_sequence = NULL,
  initial_quantity = NULL,
  seed = NULL,
  sequence_method = c("sann", "random", "sfta"),
  sfta_control = list()
)

Arguments

Details

The sequence method does not change how quantitative levels are generated. With a fixed seed, all three methods use the same quantitative design before the final row-alignment step. The alignment can change its row order but not its values or Latin-hypercube structure.

Value

An object of class magp_initial_design. Its design element is an n by 2*q matrix ready to use as the input to a MAGP fitting function. The first q columns contain quantitative levels and the last q columns contain sequence positions. The component searches and their criterion values are retained for inspection.

References

Xiao, Q., Wang, Y., Mandal, A., and Deng, X. (2024). Modeling and Active Learning for Experiments with Quantitative-Sequence Factors. Journal of the American Statistical Association. doi:10.1080/01621459.2022.2123335.

Examples

design <- magp_initial_design(
  n = 8,
  q = 4,
  sequence_maxit = 500,
  quantity_maxit = 500,
  alignment_maxit = 500,
  seed = 1
)
design$design
design$criteria

random <- magp_initial_design(
  8, 4, sequence_method = "random", seed = 1,
  quantity_maxit = 100, alignment_maxit = 100
)
sfta <- magp_initial_design(
  8, 4, sequence_method = "sfta", seed = 1, sequence_maxit = 200,
  quantity_maxit = 100, alignment_maxit = 100,
  sfta_control = list(nstarts = 2, ncalibrate = 50)
)
rbind(random = random$criteria, sfta = sfta$criteria)

Evaluate a complete quantitative-sequence initial design

Description

Combines Euclidean distance for the quantitative portion with Hamming distance for the sequence portion. Smaller values indicate better joint separation of the design runs.

Usage

magp_joint_criterion(
  quantity,
  sequence,
  quantity_weight = 0.5,
  sequence_weight = 0.5,
  p = 15L
)

Arguments

Value

One numeric criterion value. Smaller values are preferred.

References

Xiao, Q., Wang, Y., Mandal, A., and Deng, X. (2024). Modeling and Active Learning for Experiments with Quantitative-Sequence Factors. Journal of the American Statistical Association. doi:10.1080/01621459.2022.2123335.

Examples

quantity <- cbind(
  c(0.125, 0.375, 0.625, 0.875),
  c(0.625, 0.125, 0.875, 0.375),
  c(0.375, 0.875, 0.125, 0.625),
  c(0.875, 0.625, 0.375, 0.125)
)
sequence <- rbind(
  c(1, 2, 3, 4),
  c(3, 1, 4, 2),
  c(2, 4, 1, 3),
  c(4, 3, 2, 1)
)
magp_joint_criterion(quantity, sequence)

Select the next quantitative-sequence experiment

Description

Searches the allowed quantitative values and sequence permutations, then returns the unobserved point with the largest expected improvement.

Usage

magp_next_point(
  object,
  direction = c("minimize", "maximize"),
  xi = 0,
  best = NULL,
  lower = NULL,
  upper = NULL,
  sequences = NULL,
  max_sequences = 120L,
  n_starts = 5L,
  workers = 1L,
  maxit = 100L,
  factr = 1e+07,
  pgtol = 0,
  exclude_observed = TRUE,
  duplicate_tolerance = sqrt(.Machine$double.eps),
  seed = NULL
)

Arguments

Value

An object of class magp_next_point containing the selected point, its prediction, and search diagnostics.

Sequence search

If sequences is supplied, only those rows are searched. Otherwise, the function searches every permutation when there are no more than max_sequences. For a larger sequence space, it searches a reproducible sample of max_sequences permutations when seed is supplied.

Reading the result

The returned object contains:

• call: the function call;

• point: the selected input as a one-row data frame;

• expected_improvement: the score of the selected point;

• predicted_mean and predicted_standard_error: the MaGP prediction;

• direction, best_observed, and xi: the expected-improvement settings;

• bounds: the quantitative lower and upper bounds;

• sequences and sequence_source: the permutations searched and how they were obtained;

• selected_sequence and selected_start: the winning search indices;

• n_starts: the number of quantitative starts per sequence;

• workers_requested, workers_used, and execution: the parallel-search settings actually used; and

• diagnostics: the result of every sequence and starting-point search.

Examples

train <- read.table(
  system.file("extdata", "example_train.txt", package = "magp"),
  header = TRUE
)
fit <- magp2d_fit(train, seed = 1)
next_run <- magp_next_point(
  fit,
  direction = "minimize",
  n_starts = 2,
  maxit = 20,
  seed = 2
)
next_run$point
next_run$expected_improvement

Evaluate a quantitative Latin hypercube

Description

Measures the separation of design rows under Euclidean distance. The criterion is an inverse-distance p-norm, so smaller values indicate a more space-filling design.

Usage

magp_quantitative_criterion(quantity, p = 15L)

Arguments

Value

One numeric criterion value. Smaller values are preferred. A design with duplicate rows has value Inf.

Examples

quantity <- cbind(
  c(0.125, 0.375, 0.625, 0.875),
  c(0.625, 0.125, 0.875, 0.375)
)
magp_quantitative_criterion(quantity)

Construct the quantitative portion of an initial design

Description

Uses simulated annealing to improve a Latin hypercube under a maximin-style Euclidean-distance criterion. A candidate is created by swapping two values within one column, which preserves the Latin hypercube property.

Usage

magp_quantitative_design(
  n,
  q,
  p = 15L,
  maxit = 10000L,
  temp = 0.01,
  tmax = 10L,
  initial = NULL,
  seed = NULL
)

Arguments

Value

An object of class magp_quantitative_design, containing the optimized Latin hypercube, the starting design, criterion values, minimum distances, and search settings.

References

Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P. (1983). Optimization by Simulated Annealing. Science, 220, 671-680. doi:10.1126/science.220.4598.671.

Examples

design <- magp_quantitative_design(
  n = 8,
  q = 4,
  maxit = 500,
  seed = 1
)
design$quantity
design$criterion

Evaluate a sequence initial design

Description

Measures two properties of a sequence design: how evenly ordered adjacent component pairs are represented and how well separated the design rows are under Hamming distance. Smaller values indicate a better design.

Usage

magp_sequence_criterion(
  sequence,
  pair_weight = 0.2,
  space_weight = 0.8,
  p = 15L
)

Arguments

Details

Each row uses the same sequence format as the fitting functions. The value in column j is the position assigned to component j, and every row must be a permutation of 1:q.

Value

One numeric criterion value. Smaller values are preferred.

References

Xiao, Q., Wang, Y., Mandal, A., and Deng, X. (2024). Modeling and Active Learning for Experiments with Quantitative-Sequence Factors. Journal of the American Statistical Association. doi:10.1080/01621459.2022.2123335.

Examples

sequence <- rbind(
  c(1, 2, 3, 4),
  c(3, 1, 4, 2),
  c(2, 4, 1, 3),
  c(4, 3, 2, 1)
)
magp_sequence_criterion(sequence)

Construct the sequence portion of an initial design

Description

Constructs random sequences, or searches for a sequence design with balanced ordered adjacent pairs and well-separated rows using simulated annealing or space-filling threshold accepting (SFTA). A neighbor is generated by selecting one design row and exchanging two component positions, so every candidate remains a valid permutation.

Usage

magp_sequence_design(
  n,
  q,
  pair_weight = 0.2,
  space_weight = 0.8,
  p = 15L,
  maxit = 10000L,
  temp = 0.1,
  tmax = 10L,
  initial = NULL,
  seed = NULL,
  method = c("sann", "random", "sfta"),
  sfta_control = list()
)

Arguments

Details

This function constructs only the sequence portion of a quantitative- sequence initial design. Use magp_initial_design() to combine it with a quantitative Latin hypercube and improve the pairing of the two portions.

The SFTA option applies threshold accepting to the complete sequence-design criterion returned by magp_sequence_criterion(). It does not require responses and does not change the expected-improvement search used later in Bayesian optimization.

Phase I constructs nstarts complete designs by accepting candidate rows with probability equal to their minimum Hamming distance to already accepted rows divided by q. The best complete design, including the original starting design, enters Phase II. Sampling is bounded; when its budget or all ⁠q!⁠ permutations are exhausted, remaining rows are sampled uniformly. Repeated sequences are allowed (and unavoidable when ⁠n > q!⁠). The diagnostics record these fallback counts.

Phase II swaps two entries of one execution-order row at a time. Thresholds use type-7 empirical quantiles of absolute criterion differences at the fixed Phase-I winner, at probabilities 0.5 * (1 - r/nrounds). The last threshold is explicitly set to zero for a final greedy round; a move is accepted only when its criterion increase is strictly below the threshold. The best visited design is retained, so the result cannot be worse than initial. Adjacent-pair counts and Hamming-distance histograms are updated in C++ after each swap. Sequence columns always contain component positions; execution-order conversion is internal.

temp and tmax only affect "sann". Random generation avoids duplicate rows when ⁠n <= q!⁠. Neither optimized method guarantees unique rows or a global optimum.

Value

An object of class magp_sequence_design, containing the optimized sequence matrix, the starting matrix, criterion values, and search settings. The sequence matrix is ready to use as the sequence half of a magp input matrix. method_key records the method selector. For SFTA, sfta contains resolved controls, Phase-I criteria and sampling counts, thresholds, accepted moves, best-criterion history, and evaluation counts (including threshold calibration and the final independent check).

References

Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P. (1983). Optimization by Simulated Annealing. Science, 220, 671-680. doi:10.1126/science.220.4598.671.

Xiao, Q., Wang, Y., Mandal, A., and Deng, X. (2024). Modeling and Active Learning for Experiments with Quantitative-Sequence Factors. Journal of the American Statistical Association. doi:10.1080/01621459.2022.2123335.

Examples

design <- magp_sequence_design(
  n = 8,
  q = 4,
  maxit = 500,
  seed = 1
)
design$sequence
design$criterion

random <- magp_sequence_design(8, 4, method = "random", seed = 1)
sfta <- magp_sequence_design(
  8, 4, method = "sfta", maxit = 200, seed = 1,
  sfta_control = list(nstarts = 2, ncalibrate = 50)
)
c(random = random$criterion, sfta = sfta$criterion)

Fit a MaGP model with a full sequence map

Description

Fits the quantitative-sequence model with q - 1 latent mapping dimensions, giving the sequence positions a less constrained coordinate representation.

Usage

magpfull_fit(
  X,
  y = NULL,
  q = NULL,
  tau = 0.001,
  maxeval = 500,
  xtol_rel = 1e-05,
  lb_sigma = 10,
  ub_sigma = 1000,
  lb_theta = 0.5,
  ub_theta = 1000,
  lb_delta = -1,
  ub_delta = 1,
  seed = NULL,
  n_starts = 1,
  workers = 1
)

Arguments

Details

The data layout, quantitative scaling, and covariance construction are the same as in magp2d_fit(). The difference is the number of latent coordinates used to represent the sequence positions. With more than one start, the function keeps the converged result with the lowest objective and records the outcome of every start in fit$multistart$starts.

Value

An object of class magpfull.

Examples

train <- read.table(
  system.file("extdata", "example_train.txt", package = "magp"),
  header = TRUE
)
fit <- magpfull_fit(train, seed = 1, n_starts = 2)
fit

Predict outcomes from a fitted MaGP model

Description

Returns predictions for new quantitative-sequence inputs. Plug-in standard errors and variances can be returned with the predictions.

Usage

## S3 method for class 'magp2d'
predict(
  object,
  newdata,
  se.fit = FALSE,
  type = c("script", "response", "latent"),
  ...
)

## S3 method for class 'magpfull'
predict(
  object,
  newdata,
  se.fit = FALSE,
  type = c("script", "response", "latent"),
  ...
)

Arguments

Value

If se.fit = FALSE, a numeric vector of predictions. Otherwise, a list with components fit, se.fit, variance, and type.

Examples

train <- read.table(
  system.file("extdata", "example_train.txt", package = "magp"),
  header = TRUE
)
test <- read.table(
  system.file("extdata", "example_test.txt", package = "magp"),
  header = TRUE
)
fit <- magp2d_fit(train, seed = 1)
predict(fit, test[1:3, ], se.fit = TRUE, type = "response")

Summarize a fitted two-dimensional MaGP model

Description

Prints the model size, nugget variance, fitted mean, objective value, and optimizer status.

Usage

## S3 method for class 'magp2d'
print(x, ...)

Arguments

Value

x, invisibly.

Print a MaGP Bayesian optimization result

Description

Print a MaGP Bayesian optimization result

Usage

## S3 method for class 'magp_bayes_opt'
print(x, ...)

Arguments

Value

x, invisibly.

Print a quantitative-sequence initial design

Description

Print a quantitative-sequence initial design

Usage

## S3 method for class 'magp_initial_design'
print(x, ...)

Arguments

Value

x, invisibly.

Print a MaGP acquisition-search result

Description

Print a MaGP acquisition-search result

Usage

## S3 method for class 'magp_next_point'
print(x, ...)

Arguments

Value

x, invisibly.

Print a quantitative initial design

Description

Print a quantitative initial design

Usage

## S3 method for class 'magp_quantitative_design'
print(x, ...)

Arguments

Value

x, invisibly.

Print a sequence initial design

Description

Print a sequence initial design

Usage

## S3 method for class 'magp_sequence_design'
print(x, ...)

Arguments

Value

x, invisibly.

Summarize a fitted full-mapping MaGP model

Description

Prints the mapping dimension, model size, nugget variance, fitted mean, objective value, and optimizer status.

Usage

## S3 method for class 'magpfull'
print(x, ...)

Arguments

Value

x, invisibly.