R/ALC.R

In deepICMGP: Deep Intrinsic Coregionalization Gaussian Processes for Multiple Outputs

# Function Contents -----------------------------------------------------------
# Internal:
#   alc.C: calculates ALC using C for deepICMGP
# External (see documentation below):
#   ALC_deepICMGP

# ALC C -----------------------------------------------------------------------
#' @useDynLib deepICMGP, .registration = TRUE

#' @title Internal R wrapper for the registered C active-learning routine.
#'
#' @description Validates latent coordinates and kernel quantities before calling
#' alcICMGP_R. User-facing active learning is provided by ALC_deepICMGP().
#' @noRd

alc.C <- function(X, Ki, theta, g, Xcand, Xref, tau2, quad, q, trace = 0) {
  X <- numeric_matrix(X, "X")
  Ki <- numeric_matrix(Ki, "Ki")
  Xcand <- numeric_matrix(Xcand, "Xcand")
  Xref <- numeric_matrix(Xref, "Xref")
  if (ncol(Xcand) != ncol(X) || ncol(Xref) != ncol(X) ||
      !identical(dim(Ki), c(nrow(X), nrow(X)))) stop("incompatible native input dimensions")
  scalar_number(theta, "theta", strict=TRUE)
  scalar_number(g, "g", strict=TRUE)
  scalar_number(tau2, "tau2", strict=TRUE)
  scalar_number(q, "q", 1, integer=TRUE)
  scalar_number(trace, "trace", integer=TRUE)
  if (!is.numeric(quad) || length(quad) != nrow(Xref) || any(!is.finite(quad))) stop("invalid quad")

  result <- .C("alcICMGP_R",
               X = as.double(t(X)),
               n = as.integer(nrow(X)),
               col = as.integer(ncol(X)),
               Ki = as.double(t(Ki)),
               d = as.double(theta),
               g = as.double(g),
               ncand = as.integer(nrow(Xcand)),
               Xcand = as.double(t(Xcand)),
               nref = as.integer(nrow(Xref)),
               Xref = as.double(t(Xref)),
               phi = as.double(tau2),
               quad = as.double(quad),
               q = as.integer(q),
               trace = as.integer(trace),
               alc = double(nrow(Xcand)),
               PACKAGE = "deepICMGP")
  return(result$alc)
}

# ALC deepICMGP Function ------------------------------------------------------
#' @title Active Learning Cohn for deepICMGP
#' @description Select the next design point by minimizing the Active Learning
#' Cohn (ALC) criterion implemented in deepICMGP over a set of candidate inputs.
#' The criterion is averaged over reference inputs and retained MCMC samples.
#'
#' @param object An object of class \code{deepICMGP} returned by
#' \code{fit_deepICMGP} or \code{trim_deepICMGP}, containing at least two
#' retained MCMC samples.
#' @param Xcand Numeric matrix of candidate inputs, with one row per candidate
#' and the same column order and original units as the training inputs.
#' A vector is accepted for one input dimension. Default is \code{NULL}, which
#' generates 100 candidates using \code{maximinLHS(100, d)}, where \code{d}
#' is the number of input dimensions, and maps them to the observed training
#' input ranges. 
#' @param Xref Numeric matrix of reference inputs over which the criterion
#' is averaged. Use the same column order and original units as the
#' training inputs. A vector is accepted for one input dimension.
#' Default is \code{NULL}, which uses \code{Xcand}. The number of reference
#' inputs can be different from the number of candidates.
#'
#' @details To select the next design point that provides information about
#' all outputs jointly, ALC uses a determinant-based measure of the predictive
#' uncertainty remaining after adding a candidate. This measure is averaged
#' over reference inputs and retained MCMC samples. The candidate with the
#' smallest value is selected.
#' 
#' The returned \code{x_next} is the row of \code{Xcand} corresponding to the 
#' lowest ALC value. Candidate inputs and the selected point are returned in 
#' their original units.
#'
#' @return A list containing:
#' \itemize{
#' \item \code{value}: Numeric vector of ALC values at the candidate points,
#' in the same row order as \code{Xcand}.
#' \item \code{Xcand}: Matrix of candidate inputs supplied by the user or
#' generated by the function, in original input units.
#' \item \code{x_next}: One-row matrix containing the candidate with the
#' smallest ALC value, in original input units.
#' \item \code{time}: Computation time in seconds.
#' }
#'
#' @importFrom lhs maximinLHS
#'
#' @section Code provenance:
#' This implementation adapts code from the \pkg{deepgp} package and extends
#' it to deep intrinsic coregionalization models. See Sauer, Gramacy, and
#' Higdon (2023) for the underlying deep GP methodology, and
#' \url{https://CRAN.R-project.org/package=deepgp} for the source package.
#'
#' @references Chang, C.-Y. and Sung, C.-L. (2026). Deep Intrinsic Coregionalization Multi-Output Gaussian Process Surrogate with Active Learning. International Journal for Uncertainty Quantification, 16(4), 29--52. \doi{10.1615/Int.J.UncertaintyQuantification.2026062570}.
#'
#' Sauer, A., Gramacy, R. B., and Higdon, D. (2023). Active Learning for Deep Gaussian Process Surrogates. Technometrics, 65(1), 4--18. \doi{10.1080/00401706.2021.2008505}.
#'
#' @seealso \code{\link{pred_deepICMGP}}, \code{\link{trim_deepICMGP}}
#'
#' @examples
#' library(lhs)
#'
#' # Branin Function: 2d input with 3 outputs
#' f1 <- function(x) {
#'   x1 <- x[, 1] * 15 - 5
#'   x2 <- x[, 2] * 15
#'   10 * sqrt(
#'     (-1.275 * (1.2 * x1 + 0.4)^2 / pi^2 +
#'        5 * (1.2 * x1 + 0.4) / pi + (1.2 * x2 + 0.4) - 6)^2 +
#'       (10 - 5 / (4 * pi)) * cos(1.2 * x1 + 0.4) + 10
#'   ) + 2 * (1.2 * x1 + 1.9) - 3 * (3 * (1.2 * x2 + 2.4) - 1) -
#'     1 - 3 * x2 + 1
#' }
#' f2 <- function(x) {
#'   x1 <- x[, 1] * 15 - 5
#'   x2 <- x[, 2] * 15
#'   10 * sqrt(
#'     (-1.275 * (x1 + 2)^2 / pi^2 + 5 * (x1 + 2) / pi + (x2 + 2) - 6)^2 +
#'       (10 - 5 / (4 * pi)) * cos(x1 + 2) + 10
#'   ) + 2 * (x1 - 0.5) - 3 * (3 * x2 - 1) - 1
#' }
#' f3 <- function(x) {
#'   x1 <- x[, 1] * 15 - 5
#'   x2 <- x[, 2] * 15
#'   (-1.275 * x1^2 / pi^2 + 5 * x1 / pi + x2 - 6)^2 +
#'     (10 - 5 / (4 * pi)) * cos(x1) + 10
#' }
#' response <- function(x) cbind(f1(x), f2(x), f3(x))
#' set.seed(2025)
#' X <- maximinLHS(30, 2)
#' Y <- response(X)
#'
#' Xcand <- maximinLHS(50, 2)
#'
#' # Fit a deepICMGP model with 500 MCMC samples and set nugget to 1e-8 for deterministic case.
#' # For actual analysis, use a longer MCMC chain and assess convergence.
#' fit <- fit_deepICMGP(X, Y, nmcmc = 500, true_g = 1e-8)
#' 
#' # Compute ALC and select the next design point from Xcand
#' alc <- ALC_deepICMGP(fit, Xcand = Xcand)
#' alc$x_next
#'
#' \donttest{
#' # Optionally, supply a separate reference grid.
#' Xref <- maximinLHS(200, 2)
#' alc <- ALC_deepICMGP(fit, Xcand = Xcand, Xref = Xref)
#' alc$x_next
#'
#' # If candidates are not given, 100 candidate inputs will be generated automatically.
#' alc <- ALC_deepICMGP(fit)
#' nrow(alc$Xcand)
#' alc$x_next
#' }
#' @export


ALC_deepICMGP <- function(object, Xcand = NULL, Xref = NULL) {
  tic <- proc.time()[3]
  check_object(object, 2L)
  if (any(c("x_new", "mean", "s2", "C2", "w_new", "quad", "s2_all") %in% names(object)))
    stop("object must be a fitted model, not a prediction object. Use the fit returned by fit_deepICMGP() or trim_deepICMGP().")
  d <- ncol(object$x)
  if (is.null(Xcand)) {
    # Draw in the unit cube, then restore the observed training input ranges.
    # The observed span is zero for constant columns, whose scaling divisor is one.
    input_span <- (apply(object$x, 2L, max) - apply(object$x, 2L, min)) * object$scaling$x_range
    Xcand <- sweep(sweep(maximinLHS(100, d), 2L, input_span, "*"),
                   2L, object$scaling$x_min, "+")
    colnames(Xcand) <- colnames(object$x)
  }
  Xcand <- numeric_matrix(Xcand, "Xcand", vector_ok=TRUE)
  if (ncol(Xcand) != d) stop("Xcand must have the same number of columns as x")
  if (is.null(Xref)) Xref <- Xcand
  Xref <- numeric_matrix(Xref, "Xref", vector_ok=TRUE)
  if (ncol(Xref) != d) stop("Xref must have the same number of columns as x")

  candidate <- pred_deepICMGP(object, Xcand, store_latent=TRUE)
  if (nrow(Xref) == nrow(Xcand) && all(Xref == Xcand)) {
    reference <- candidate
  } else {
    reference <- pred_deepICMGP(object, Xref, store_latent=TRUE)
  }
  nt <- object$nmcmc
  value <- numeric(nrow(Xcand))
  for (t in seq_len(nt)) {
    w <- object$w[[t]]
    w_cand <- candidate$w_new[[t]]
    w_ref <- reference$w_new[[t]]
    K <- kernel_icm(sq_dist(w), 1, object$theta_y[t], object$g[t])
    value <- value + alc.C(w, invdet_icm(K)$Mi, object$theta_y[t], object$g[t],
                           w_cand, w_ref, 1, reference$quad[,t], ncol(object$y))
  }
  value <- value / nt
  list(value=value, Xcand=Xcand, x_next=Xcand[which.min(value), , drop=FALSE],
       time=unname(proc.time()[3]-tic))
}