Nothing
R/ALMC.R
In RNAmf: Recursive Non-Additive Emulator for Multi-Fidelity Data
Defines functions
ALMC_RNAmf
Documented in
ALMC_RNAmf
#' @title find the next point by ALMC criterion
#'
#' @description The function acquires the new point by the hybrid approach,
#' referred to as Active learning MacKay-Cohn (ALMC) criterion.
#' It finds the optimal input location \eqn{\bm{x^*}}
#' by maximizing \eqn{\sigma^{*2}_L(\bm{x})},
#' the posterior predictive variance at the highest-fidelity level \eqn{L}.
#' After selecting \eqn{\bm{x^*}},
#' it finds the optimal fidelity level by maximizing ALC criterion at \eqn{\bm{x^*}},
#' \eqn{\text{argmax}_{l\in\{1,\ldots,L\}} \frac{\Delta \sigma_L^{2}(l,\bm{x^*})}{\sum^l_{j=1}C_j}},
#' where \eqn{C_j} is the simulation cost at level \eqn{j}.
#' See \code{\link{ALC_RNAmf}}.
#' For details, see Heo and Sung (2025, <\doi{https://doi.org/10.1080/00401706.2024.2376173}>).
#'
#' A new point is acquired on \code{Xcand}. If \code{Xcand=NULL} and \code{Xref=NULL}, a new point is acquired on unit hypercube \eqn{[0,1]^d}.
#'
#' @param Xref vector or matrix of reference location to approximate the integral of ALC. If \code{Xref=NULL}, \eqn{100 \times d} points from 0 to 1 are generated by Latin hypercube design. Default is \code{NULL}.
#' @param Xcand vector or matrix of candidate set which could be added into the current design only when \code{optim=FALSE}. \code{Xcand} is the set of the points where ALM criterion is evaluated. If \code{Xcand=NULL}, \code{Xref} is used. Default is \code{NULL}.
#' @param fit object of class \code{RNAmf}.
#' @param mc.sample a number of mc samples generated for the imputation through MC approximation. Default is \code{100}.
#' @param cost vector of the costs for each level of fidelity. If \code{cost=NULL}, total costs at all levels would be 1. \code{cost} is encouraged to have a ascending order of positive value. Default is \code{NULL}.
#' @param optim logical indicating whether to optimize AL criterion by \code{optim}'s gradient-based \code{L-BFGS-B} method. If \code{optim=TRUE}, \eqn{5 \times d} starting points are generated by Latin hypercube design for optimization. If \code{optim=FALSE}, AL criterion is optimized on the \code{Xcand}. Default is \code{TRUE}.
#' @param parallel logical indicating whether to compute the AL criterion in parallel or not. If \code{parallel=TRUE}, parallel computation is utilized. Default is \code{FALSE}.
#' @param ncore a number of core for parallel. It is only used if \code{parallel=TRUE}. Default is 1.
#' @param trace logical indicating whether to print the computational time for each step. If \code{trace=TRUE}, the computation time for each step is printed. Default is \code{TRUE}.
#' @return
#' \itemize{
#' \item \code{ALMC}: vector of ALMC criterion \eqn{ \frac{\Delta \sigma_L^{2}(l,\bm{x^*})}{\sum^l_{j=1}C_j}} for \eqn{1\leq l\leq L}.
#' \item \code{ALM}: vector of ALM criterion computed at each point of \code{Xcand} at the highest fidelity level if \code{optim=FALSE}. If \code{optim=TRUE}, \code{ALM} returns \code{NULL}.
#' \item \code{ALC}: list of ALC criterion integrated on \code{Xref} when each data point on \code{Xcand} is added at each level \eqn{l} if \code{optim=FALSE}. If \code{optim=TRUE}, \code{ALC} returns \code{NULL}.
#' \item \code{cost}: a copy of \code{cost}.
#' \item \code{Xcand}: a copy of \code{Xcand}.
#' \item \code{chosen}: list of chosen level and point.
#' \item \code{time}: a scalar of the time for the computation.
#' }
#' @importFrom plgp covar.sep
#' @importFrom stats rnorm
#' @importFrom lhs randomLHS
#' @importFrom lhs maximinLHS
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @importFrom doParallel registerDoParallel
#' @importFrom doParallel stopImplicitCluster
#' @usage ALMC_RNAmf(Xref = NULL, Xcand = NULL, fit, mc.sample = 100,
#' cost = NULL, optim = TRUE, parallel = FALSE, ncore = 1, trace=TRUE)
#' @export
#' @examples
#' \donttest{
#' library(lhs)
#' library(doParallel)
#' library(foreach)
#'
#' ### simulation costs ###
#' cost <- c(1, 3)
#'
#' ### 1-d Perdikaris function in Perdikaris, et al. (2017) ###
#' # low-fidelity function
#' f1 <- function(x) {
#' sin(8 * pi * x)
#' }
#'
#' # high-fidelity function
#' f2 <- function(x) {
#' (x - sqrt(2)) * (sin(8 * pi * x))^2
#' }
#'
#' ### training data ###
#' n1 <- 13
#' n2 <- 8
#'
#' ### fix seed to reproduce the result ###
#' set.seed(1)
#'
#' ### generate initial nested design ###
#' X <- NestedX(c(n1, n2), 1)
#' X1 <- X[[1]]
#' X2 <- X[[2]]
#'
#' ### n1 and n2 might be changed from NestedX ###
#' ### assign n1 and n2 again ###
#' n1 <- nrow(X1)
#' n2 <- nrow(X2)
#'
#' y1 <- f1(X1)
#' y2 <- f2(X2)
#'
#' ### n=100 uniform test data ###
#' x <- seq(0, 1, length.out = 100)
#'
#' ### fit an RNAmf ###
#' fit.RNAmf <- RNAmf(list(X1, X2), list(y1, y2), kernel = "sqex", constant=TRUE)
#'
#' ### predict ###
#' predy <- predict(fit.RNAmf, x)$mu
#' predsig2 <- predict(fit.RNAmf, x)$sig2
#'
#' ### active learning with optim=TRUE ###
#' almc.RNAmf.optim <- ALMC_RNAmf(
#' Xref = x, Xcand = x, fit.RNAmf, cost = cost,
#' optim = TRUE, parallel = TRUE, ncore = 2
#' )
#' print(almc.RNAmf.optim$time) # computation time of optim=TRUE
#'
#' ### active learning with optim=FALSE ###
#' almc.RNAmf <- ALMC_RNAmf(
#' Xref = x, Xcand = x, fit.RNAmf, cost = cost,
#' optim = FALSE, parallel = TRUE, ncore = 2
#' )
#' print(almc.RNAmf$time) # computation time of optim=FALSE
#'
#' ### visualize ALMC ###
#' oldpar <- par(mfrow = c(1, 2))
#' plot(x, almc.RNAmf$ALM,
#' type = "l", lty = 2,
#' xlab = "x", ylab = "ALM criterion at the high-fidelity level"
#' )
#' points(almc.RNAmf$chosen$Xnext,
#' almc.RNAmf$ALM[which(x == drop(almc.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' plot(x, almc.RNAmf$ALC$ALC1,
#' type = "l", lty = 2,
#' ylim = c(min(c(almc.RNAmf$ALC$ALC1, almc.RNAmf$ALC$ALC2)),
#' max(c(almc.RNAmf$ALC$ALC1, almc.RNAmf$ALC$ALC2))),
#' xlab = "x", ylab = "ALC criterion augmented at each level on the optimal input location"
#' )
#' lines(x, almc.RNAmf$ALC$ALC2, type = "l", lty = 2)
#' points(almc.RNAmf$chosen$Xnext,
#' almc.RNAmf$ALC$ALC1[which(x == drop(almc.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' points(almc.RNAmf$chosen$Xnext,
#' almc.RNAmf$ALC$ALC2[which(x == drop(almc.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' par(oldpar)}
#'
ALMC_RNAmf <- function(Xref = NULL, Xcand = NULL, fit, mc.sample = 100, cost = NULL, optim = TRUE, parallel = FALSE, ncore = 1, trace=TRUE) {
t0 <- proc.time()
### check the object ###
if (!inherits(fit, "RNAmf")) {
stop("The object is not of class \"RNAmf\" \n")
}
if (length(cost) != fit$level) stop("The length of cost should be the level of object")
if (is.null(dim(Xref)) & length(Xref) > 0) Xref <- matrix(Xref, ncol = 1)
### ALMC ###
if (fit$level == 2) { # level 2
if (!is.null(cost) & cost[1] >= cost[2]) {
warning("If the cost for high-fidelity is cheaper, acquire the high-fidelity")
} else if (is.null(cost)) {
cost <- c(1, 0)
}
if (is.null(Xref)) Xref <- randomLHS(100 * dim(fit$fits[[1]]$X)[2], dim(fit$fits[[1]]$X)[2])
if (parallel) registerDoParallel(ncore)
Icurrent <- mean(predict(fit, Xref)$sig2[[2]])
fit1 <- f1 <- fit$fits[[1]]
fit2 <- f2 <- fit$fits[[2]]
constant <- fit$constant
kernel <- fit$kernel
g <- fit1$g
# x.center1 <- attr(fit1$X, "scaled:center")
# x.scale1 <- attr(fit1$X, "scaled:scale")
# y.center1 <- attr(fit1$y, "scaled:center")
#
# x.center2 <- attr(fit2$X, "scaled:center")
# x.scale2 <- attr(fit2$X, "scaled:scale")
# y.center2 <- attr(fit2$y, "scaled:center")
### Generate the candidate set ###
if (optim){ # optim = TRUE
Xcand <- randomLHS(5*ncol(fit1$X), ncol(fit1$X))
}else{ # optim = FALSE
if (is.null(Xcand)){
Xcand <- Xref
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = ncol(fit1$X))
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Find the next point ###
t1 <- proc.time()
if (optim) {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_2, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
return(c(-optim.ALM$value, optim.ALM$par))
}
} else {
alm.mat <- matrix(c(rep(0, nrow(Xcand))), ncol=1)
x.mat <- matrix(0, ncol=ncol(Xcand), nrow=nrow(Xcand))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
alm.mat[i,] <- -optim.ALM$value
x.mat[i,] <- optim.ALM$par
}
optm.mat <- cbind(alm.mat, x.mat)
}
Xnext <- matrix(optm.mat[which.max(optm.mat[,1]), 2:(ncol(fit1$X)+1)], nrow = 1)
x.cand <- Xnext ### Anyway it will return NULL
# Xnext <- matrix(optm.mat[ncol(fit1$X) + 1, which.max(optm.mat[(1:ncol(fit1$X)), ])], nrow = 1)
# x.cand <- matrix(optm.mat[ncol(fit1$X) + 1, ], ncol = ncol(fit1$X))
ALM <- optm.mat[1,]/c(cost[1] + cost[2])
} else {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(-obj.ALM_level_2(fit = fit, newx)))
}
} else {
optm.mat <- c(rep(0, nrow(Xcand)))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optm.mat[i] <- -obj.ALM_level_2(fit = fit, newx)
}
}
Xnext <- matrix(Xcand[which.max(optm.mat), ], nrow = 1)
x.cand <- Xcand
ALM <- optm.mat/c(cost[1] + cost[2])
}
t2 <-proc.time()
if(trace) cat("Finding the next point:", (t2 - t1)[3], "seconds\n")
### Calculate the deduced variance at selected point###
ALMC.1 <- obj.ALC_level_1(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t3 <-proc.time()
if(trace) cat("Calculating deduced variance for level 1:", (t3 - t2)[3], "seconds\n")
ALMC.2 <- obj.ALC_level_2(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t4 <-proc.time()
if(trace) cat("Calculating deduced variance for level 2:", (t4 - t3)[3], "seconds\n")
### Compute the deduced variance at all candidate points ###
if (optim) {
intvar1 <- NULL
intvar2 <- NULL
} else {
if (parallel) {
pseudointvar <- foreach(i = 1:nrow(x.cand), .combine = cbind) %dopar% {
newx <- matrix(x.cand[i, ], nrow = 1)
return(c(
obj.ALC_level_1(newx, Xref, fit, mc.sample),
obj.ALC_level_2(newx, Xref, fit, mc.sample)
))
}
intvar1 <- as.matrix(pseudointvar)[1, ]
intvar2 <- as.matrix(pseudointvar)[2, ]
} else {
intvar1 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
intvar2 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
for (i in 1:nrow(x.cand)) {
print(paste(i, nrow(x.cand), sep = "/"))
newx <- matrix(x.cand[i, ], nrow = 1)
intvar1[i] <- obj.ALC_level_1(newx, Xref, fit, mc.sample)
intvar2[i] <- obj.ALC_level_2(newx, Xref, fit, mc.sample)
}
}
}
t5 <-proc.time()
if(trace) cat("Computing deduced variance for all points:", (t5 - t4)[3], "seconds\n")
ALMCvalue <- c(Icurrent - ALMC.1, Icurrent - ALMC.2) / c(cost[1], cost[1] + cost[2])
if (ALMCvalue[2] > ALMCvalue[1]) {
level <- 2
} else {
level <- 1
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALC <- list(ALC1 = intvar1 / cost[1], ALC2 = intvar2 / (cost[1] + cost[2]))
ALMC <- c(ALMC1 = ALMC.1, ALMC2 = ALMC.2)
} else if (fit$level == 3) { # level 3
if (!is.null(cost) & (cost[1] >= cost[2] | cost[2] >= cost[3])) {
warning("If the cost for high-fidelity is cheaper, acquire the high-fidelity")
} else if (is.null(cost)) {
cost <- c(1, 0, 0)
}
if (is.null(Xref)) Xref <- randomLHS(100 * dim(fit$fits[[1]]$X)[2], dim(fit$fits[[1]]$X)[2])
if (parallel) registerDoParallel(ncore)
Icurrent <- mean(predict(fit, Xref)$sig2[[3]])
fit1 <- fit$fits[[1]]
fit2 <- fit$fits[[2]]
fit3 <- fit$fits[[3]]
constant <- fit$constant
kernel <- fit$kernel
g <- fit1$g
# x.center1 <- attr(fit1$X, "scaled:center")
# x.scale1 <- attr(fit1$X, "scaled:scale")
# y.center1 <- attr(fit1$y, "scaled:center")
#
# x.center2 <- attr(fit2$X, "scaled:center")
# x.scale2 <- attr(fit2$X, "scaled:scale")
# y.center2 <- attr(fit2$y, "scaled:center")
#
# x.center3 <- attr(fit3$X, "scaled:center")
# x.scale3 <- attr(fit3$X, "scaled:scale")
# y.center3 <- attr(fit3$y, "scaled:center")
### Generate the candidate set ###
if (optim){ # optim = TRUE
Xcand <- randomLHS(5*ncol(fit1$X), ncol(fit1$X))
}else{ # optim = FALSE
if (is.null(Xcand)){
Xcand <- Xref
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = ncol(fit1$X))
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Find the next point ###
t1 <- proc.time()
if (optim) {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
return(c(-optim.ALM$value, optim.ALM$par))
}
} else {
alm.mat <- matrix(c(rep(0, nrow(Xcand))), ncol=1)
x.mat <- matrix(0, ncol=ncol(Xcand), nrow=nrow(Xcand))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
alm.mat[i,] <- -optim.ALM$value
x.mat[i,] <- optim.ALM$par
}
optm.mat <- cbind(alm.mat, x.mat)
}
Xnext <- matrix(optm.mat[which.max(optm.mat[,1]), 2:(ncol(fit1$X)+1)], nrow = 1)
x.cand <- Xnext ### Anyway it will return NULL
# Xnext <- matrix(optm.mat[ncol(fit1$X) + 1, which.max(optm.mat[(1:ncol(fit1$X)), ])], nrow = 1)
# x.cand <- matrix(optm.mat[ncol(fit1$X) + 1, ], ncol = ncol(fit1$X))
ALM <- optm.mat[1,]/c(cost[1] + cost[2] + cost[3])
} else {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(-obj.ALM_level_3(fit = fit, newx)))
}
} else {
optm.mat <- c(rep(0, nrow(Xcand)))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optm.mat[i] <- -obj.ALM_level_3(fit = fit, newx)
}
}
Xnext <- matrix(Xcand[which.max(optm.mat), ], nrow = 1)
x.cand <- Xcand
ALM <- optm.mat/c(cost[1] + cost[2] + cost[3])
}
t2 <-proc.time()
if(trace) cat("Finding the next point:", (t2 - t1)[3], "seconds\n")
### Calculate the deduced variance ###
ALMC.1 <- obj.ALC_level_3_1(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t3 <-proc.time()
if(trace) cat("Calculating deduced variance for level 1:", (t3 - t2)[3], "seconds\n")
ALMC.2 <- obj.ALC_level_3_2(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t4 <-proc.time()
if(trace) cat("Calculating deduced variance for level 2:", (t4 - t3)[3], "seconds\n")
ALMC.3 <- obj.ALC_level_3_3(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t5 <-proc.time()
if(trace) cat("Calculating deduced variance for level 3:", (t5 - t4)[3], "seconds\n")
### Compute the deduced variance at all candidate points ###
if (optim) {
intvar1 <- NULL
intvar2 <- NULL
intvar3 <- NULL
} else {
if (parallel) {
pseudointvar <- foreach(i = 1:nrow(x.cand), .combine = cbind) %dopar% {
newx <- matrix(x.cand[i, ], nrow = 1)
return(c(
obj.ALC_level_3_1(newx, Xref, fit, mc.sample),
obj.ALC_level_3_2(newx, Xref, fit, mc.sample),
obj.ALC_level_3_3(newx, Xref, fit, mc.sample)
))
}
intvar1 <- as.matrix(pseudointvar)[1, ]
intvar2 <- as.matrix(pseudointvar)[2, ]
intvar3 <- as.matrix(pseudointvar)[3, ]
} else {
intvar1 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
intvar2 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
intvar3 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
for (i in 1:nrow(x.cand)) {
print(paste(i, nrow(x.cand), sep = "/"))
newx <- matrix(x.cand[i, ], nrow = 1)
intvar1[i] <- obj.ALC_level_3_1(newx, Xref, fit, mc.sample)
intvar2[i] <- obj.ALC_level_3_2(newx, Xref, fit, mc.sample)
intvar3[i] <- obj.ALC_level_3_3(newx, Xref, fit, mc.sample)
}
}
}
t6 <-proc.time()
if(trace) cat("Computing deduced variance for all points:", (t6 - t5)[3], "seconds\n")
ALMCvalue <- c(Icurrent - ALMC.1, Icurrent - ALMC.2, Icurrent - ALMC.3) / c(cost[1], cost[1] + cost[2], cost[1] + cost[2] + cost[3])
if (ALMCvalue[3] > ALMCvalue[2]) {
level <- 3
} else if (ALMCvalue[2] > ALMCvalue[1]) {
level <- 2
} else {
level <- 1
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALC <- list(ALC1 = intvar1 / cost[1], ALC2 = intvar2 / (cost[1] + cost[2]), ALC3 = intvar3 / (cost[1] + cost[2] + cost[3]))
ALMC <- c(ALMC1 = ALMC.1, ALMC2 = ALMC.2, ALMC3 = ALMC.3)
} else {
stop("level is missing")
}
if (parallel) stopImplicitCluster()
if (optim) ALC <- ALM <- NULL
return(list(ALMC = ALMC, ALM = ALM, ALC=ALC,
cost = cost, Xcand = Xcand, chosen = chosen, time = (proc.time() - t0)[3]))
}
Try the RNAmf package in your browser
Any scripts or data that you put into this service are public.
RNAmf documentation built on June 8, 2025, 10:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ALMC_RNAmf
Documented in ALMC_RNAmf
#' @title find the next point by ALMC criterion
#'
#' @description The function acquires the new point by the hybrid approach,
#' referred to as Active learning MacKay-Cohn (ALMC) criterion.
#' It finds the optimal input location \eqn{\bm{x^*}}
#' by maximizing \eqn{\sigma^{*2}_L(\bm{x})},
#' the posterior predictive variance at the highest-fidelity level \eqn{L}.
#' After selecting \eqn{\bm{x^*}},
#' it finds the optimal fidelity level by maximizing ALC criterion at \eqn{\bm{x^*}},
#' \eqn{\text{argmax}_{l\in\{1,\ldots,L\}} \frac{\Delta \sigma_L^{2}(l,\bm{x^*})}{\sum^l_{j=1}C_j}},
#' where \eqn{C_j} is the simulation cost at level \eqn{j}.
#' See \code{\link{ALC_RNAmf}}.
#' For details, see Heo and Sung (2025, <\doi{https://doi.org/10.1080/00401706.2024.2376173}>).
#'
#' A new point is acquired on \code{Xcand}. If \code{Xcand=NULL} and \code{Xref=NULL}, a new point is acquired on unit hypercube \eqn{[0,1]^d}.
#'
#' @param Xref vector or matrix of reference location to approximate the integral of ALC. If \code{Xref=NULL}, \eqn{100 \times d} points from 0 to 1 are generated by Latin hypercube design. Default is \code{NULL}.
#' @param Xcand vector or matrix of candidate set which could be added into the current design only when \code{optim=FALSE}. \code{Xcand} is the set of the points where ALM criterion is evaluated. If \code{Xcand=NULL}, \code{Xref} is used. Default is \code{NULL}.
#' @param fit object of class \code{RNAmf}.
#' @param mc.sample a number of mc samples generated for the imputation through MC approximation. Default is \code{100}.
#' @param cost vector of the costs for each level of fidelity. If \code{cost=NULL}, total costs at all levels would be 1. \code{cost} is encouraged to have a ascending order of positive value. Default is \code{NULL}.
#' @param optim logical indicating whether to optimize AL criterion by \code{optim}'s gradient-based \code{L-BFGS-B} method. If \code{optim=TRUE}, \eqn{5 \times d} starting points are generated by Latin hypercube design for optimization. If \code{optim=FALSE}, AL criterion is optimized on the \code{Xcand}. Default is \code{TRUE}.
#' @param parallel logical indicating whether to compute the AL criterion in parallel or not. If \code{parallel=TRUE}, parallel computation is utilized. Default is \code{FALSE}.
#' @param ncore a number of core for parallel. It is only used if \code{parallel=TRUE}. Default is 1.
#' @param trace logical indicating whether to print the computational time for each step. If \code{trace=TRUE}, the computation time for each step is printed. Default is \code{TRUE}.
#' @return
#' \itemize{
#' \item \code{ALMC}: vector of ALMC criterion \eqn{ \frac{\Delta \sigma_L^{2}(l,\bm{x^*})}{\sum^l_{j=1}C_j}} for \eqn{1\leq l\leq L}.
#' \item \code{ALM}: vector of ALM criterion computed at each point of \code{Xcand} at the highest fidelity level if \code{optim=FALSE}. If \code{optim=TRUE}, \code{ALM} returns \code{NULL}.
#' \item \code{ALC}: list of ALC criterion integrated on \code{Xref} when each data point on \code{Xcand} is added at each level \eqn{l} if \code{optim=FALSE}. If \code{optim=TRUE}, \code{ALC} returns \code{NULL}.
#' \item \code{cost}: a copy of \code{cost}.
#' \item \code{Xcand}: a copy of \code{Xcand}.
#' \item \code{chosen}: list of chosen level and point.
#' \item \code{time}: a scalar of the time for the computation.
#' }
#' @importFrom plgp covar.sep
#' @importFrom stats rnorm
#' @importFrom lhs randomLHS
#' @importFrom lhs maximinLHS
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @importFrom doParallel registerDoParallel
#' @importFrom doParallel stopImplicitCluster
#' @usage ALMC_RNAmf(Xref = NULL, Xcand = NULL, fit, mc.sample = 100,
#' cost = NULL, optim = TRUE, parallel = FALSE, ncore = 1, trace=TRUE)
#' @export
#' @examples
#' \donttest{
#' library(lhs)
#' library(doParallel)
#' library(foreach)
#'
#' ### simulation costs ###
#' cost <- c(1, 3)
#'
#' ### 1-d Perdikaris function in Perdikaris, et al. (2017) ###
#' # low-fidelity function
#' f1 <- function(x) {
#' sin(8 * pi * x)
#' }
#'
#' # high-fidelity function
#' f2 <- function(x) {
#' (x - sqrt(2)) * (sin(8 * pi * x))^2
#' }
#'
#' ### training data ###
#' n1 <- 13
#' n2 <- 8
#'
#' ### fix seed to reproduce the result ###
#' set.seed(1)
#'
#' ### generate initial nested design ###
#' X <- NestedX(c(n1, n2), 1)
#' X1 <- X[[1]]
#' X2 <- X[[2]]
#'
#' ### n1 and n2 might be changed from NestedX ###
#' ### assign n1 and n2 again ###
#' n1 <- nrow(X1)
#' n2 <- nrow(X2)
#'
#' y1 <- f1(X1)
#' y2 <- f2(X2)
#'
#' ### n=100 uniform test data ###
#' x <- seq(0, 1, length.out = 100)
#'
#' ### fit an RNAmf ###
#' fit.RNAmf <- RNAmf(list(X1, X2), list(y1, y2), kernel = "sqex", constant=TRUE)
#'
#' ### predict ###
#' predy <- predict(fit.RNAmf, x)$mu
#' predsig2 <- predict(fit.RNAmf, x)$sig2
#'
#' ### active learning with optim=TRUE ###
#' almc.RNAmf.optim <- ALMC_RNAmf(
#' Xref = x, Xcand = x, fit.RNAmf, cost = cost,
#' optim = TRUE, parallel = TRUE, ncore = 2
#' )
#' print(almc.RNAmf.optim$time) # computation time of optim=TRUE
#'
#' ### active learning with optim=FALSE ###
#' almc.RNAmf <- ALMC_RNAmf(
#' Xref = x, Xcand = x, fit.RNAmf, cost = cost,
#' optim = FALSE, parallel = TRUE, ncore = 2
#' )
#' print(almc.RNAmf$time) # computation time of optim=FALSE
#'
#' ### visualize ALMC ###
#' oldpar <- par(mfrow = c(1, 2))
#' plot(x, almc.RNAmf$ALM,
#' type = "l", lty = 2,
#' xlab = "x", ylab = "ALM criterion at the high-fidelity level"
#' )
#' points(almc.RNAmf$chosen$Xnext,
#' almc.RNAmf$ALM[which(x == drop(almc.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' plot(x, almc.RNAmf$ALC$ALC1,
#' type = "l", lty = 2,
#' ylim = c(min(c(almc.RNAmf$ALC$ALC1, almc.RNAmf$ALC$ALC2)),
#' max(c(almc.RNAmf$ALC$ALC1, almc.RNAmf$ALC$ALC2))),
#' xlab = "x", ylab = "ALC criterion augmented at each level on the optimal input location"
#' )
#' lines(x, almc.RNAmf$ALC$ALC2, type = "l", lty = 2)
#' points(almc.RNAmf$chosen$Xnext,
#' almc.RNAmf$ALC$ALC1[which(x == drop(almc.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' points(almc.RNAmf$chosen$Xnext,
#' almc.RNAmf$ALC$ALC2[which(x == drop(almc.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' par(oldpar)}
#'
ALMC_RNAmf <- function(Xref = NULL, Xcand = NULL, fit, mc.sample = 100, cost = NULL, optim = TRUE, parallel = FALSE, ncore = 1, trace=TRUE) {
t0 <- proc.time()
### check the object ###
if (!inherits(fit, "RNAmf")) {
stop("The object is not of class \"RNAmf\" \n")
}
if (length(cost) != fit$level) stop("The length of cost should be the level of object")
if (is.null(dim(Xref)) & length(Xref) > 0) Xref <- matrix(Xref, ncol = 1)
### ALMC ###
if (fit$level == 2) { # level 2
if (!is.null(cost) & cost[1] >= cost[2]) {
warning("If the cost for high-fidelity is cheaper, acquire the high-fidelity")
} else if (is.null(cost)) {
cost <- c(1, 0)
}
if (is.null(Xref)) Xref <- randomLHS(100 * dim(fit$fits[[1]]$X)[2], dim(fit$fits[[1]]$X)[2])
if (parallel) registerDoParallel(ncore)
Icurrent <- mean(predict(fit, Xref)$sig2[[2]])
fit1 <- f1 <- fit$fits[[1]]
fit2 <- f2 <- fit$fits[[2]]
constant <- fit$constant
kernel <- fit$kernel
g <- fit1$g
# x.center1 <- attr(fit1$X, "scaled:center")
# x.scale1 <- attr(fit1$X, "scaled:scale")
# y.center1 <- attr(fit1$y, "scaled:center")
#
# x.center2 <- attr(fit2$X, "scaled:center")
# x.scale2 <- attr(fit2$X, "scaled:scale")
# y.center2 <- attr(fit2$y, "scaled:center")
### Generate the candidate set ###
if (optim){ # optim = TRUE
Xcand <- randomLHS(5*ncol(fit1$X), ncol(fit1$X))
}else{ # optim = FALSE
if (is.null(Xcand)){
Xcand <- Xref
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = ncol(fit1$X))
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Find the next point ###
t1 <- proc.time()
if (optim) {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_2, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
return(c(-optim.ALM$value, optim.ALM$par))
}
} else {
alm.mat <- matrix(c(rep(0, nrow(Xcand))), ncol=1)
x.mat <- matrix(0, ncol=ncol(Xcand), nrow=nrow(Xcand))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
alm.mat[i,] <- -optim.ALM$value
x.mat[i,] <- optim.ALM$par
}
optm.mat <- cbind(alm.mat, x.mat)
}
Xnext <- matrix(optm.mat[which.max(optm.mat[,1]), 2:(ncol(fit1$X)+1)], nrow = 1)
x.cand <- Xnext ### Anyway it will return NULL
# Xnext <- matrix(optm.mat[ncol(fit1$X) + 1, which.max(optm.mat[(1:ncol(fit1$X)), ])], nrow = 1)
# x.cand <- matrix(optm.mat[ncol(fit1$X) + 1, ], ncol = ncol(fit1$X))
ALM <- optm.mat[1,]/c(cost[1] + cost[2])
} else {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(-obj.ALM_level_2(fit = fit, newx)))
}
} else {
optm.mat <- c(rep(0, nrow(Xcand)))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optm.mat[i] <- -obj.ALM_level_2(fit = fit, newx)
}
}
Xnext <- matrix(Xcand[which.max(optm.mat), ], nrow = 1)
x.cand <- Xcand
ALM <- optm.mat/c(cost[1] + cost[2])
}
t2 <-proc.time()
if(trace) cat("Finding the next point:", (t2 - t1)[3], "seconds\n")
### Calculate the deduced variance at selected point###
ALMC.1 <- obj.ALC_level_1(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t3 <-proc.time()
if(trace) cat("Calculating deduced variance for level 1:", (t3 - t2)[3], "seconds\n")
ALMC.2 <- obj.ALC_level_2(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t4 <-proc.time()
if(trace) cat("Calculating deduced variance for level 2:", (t4 - t3)[3], "seconds\n")
### Compute the deduced variance at all candidate points ###
if (optim) {
intvar1 <- NULL
intvar2 <- NULL
} else {
if (parallel) {
pseudointvar <- foreach(i = 1:nrow(x.cand), .combine = cbind) %dopar% {
newx <- matrix(x.cand[i, ], nrow = 1)
return(c(
obj.ALC_level_1(newx, Xref, fit, mc.sample),
obj.ALC_level_2(newx, Xref, fit, mc.sample)
))
}
intvar1 <- as.matrix(pseudointvar)[1, ]
intvar2 <- as.matrix(pseudointvar)[2, ]
} else {
intvar1 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
intvar2 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
for (i in 1:nrow(x.cand)) {
print(paste(i, nrow(x.cand), sep = "/"))
newx <- matrix(x.cand[i, ], nrow = 1)
intvar1[i] <- obj.ALC_level_1(newx, Xref, fit, mc.sample)
intvar2[i] <- obj.ALC_level_2(newx, Xref, fit, mc.sample)
}
}
}
t5 <-proc.time()
if(trace) cat("Computing deduced variance for all points:", (t5 - t4)[3], "seconds\n")
ALMCvalue <- c(Icurrent - ALMC.1, Icurrent - ALMC.2) / c(cost[1], cost[1] + cost[2])
if (ALMCvalue[2] > ALMCvalue[1]) {
level <- 2
} else {
level <- 1
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALC <- list(ALC1 = intvar1 / cost[1], ALC2 = intvar2 / (cost[1] + cost[2]))
ALMC <- c(ALMC1 = ALMC.1, ALMC2 = ALMC.2)
} else if (fit$level == 3) { # level 3
if (!is.null(cost) & (cost[1] >= cost[2] | cost[2] >= cost[3])) {
warning("If the cost for high-fidelity is cheaper, acquire the high-fidelity")
} else if (is.null(cost)) {
cost <- c(1, 0, 0)
}
if (is.null(Xref)) Xref <- randomLHS(100 * dim(fit$fits[[1]]$X)[2], dim(fit$fits[[1]]$X)[2])
if (parallel) registerDoParallel(ncore)
Icurrent <- mean(predict(fit, Xref)$sig2[[3]])
fit1 <- fit$fits[[1]]
fit2 <- fit$fits[[2]]
fit3 <- fit$fits[[3]]
constant <- fit$constant
kernel <- fit$kernel
g <- fit1$g
# x.center1 <- attr(fit1$X, "scaled:center")
# x.scale1 <- attr(fit1$X, "scaled:scale")
# y.center1 <- attr(fit1$y, "scaled:center")
#
# x.center2 <- attr(fit2$X, "scaled:center")
# x.scale2 <- attr(fit2$X, "scaled:scale")
# y.center2 <- attr(fit2$y, "scaled:center")
#
# x.center3 <- attr(fit3$X, "scaled:center")
# x.scale3 <- attr(fit3$X, "scaled:scale")
# y.center3 <- attr(fit3$y, "scaled:center")
### Generate the candidate set ###
if (optim){ # optim = TRUE
Xcand <- randomLHS(5*ncol(fit1$X), ncol(fit1$X))
}else{ # optim = FALSE
if (is.null(Xcand)){
Xcand <- Xref
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = ncol(fit1$X))
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Find the next point ###
t1 <- proc.time()
if (optim) {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
return(c(-optim.ALM$value, optim.ALM$par))
}
} else {
alm.mat <- matrix(c(rep(0, nrow(Xcand))), ncol=1)
x.mat <- matrix(0, ncol=ncol(Xcand), nrow=nrow(Xcand))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optim.ALM <- optim(newx, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
alm.mat[i,] <- -optim.ALM$value
x.mat[i,] <- optim.ALM$par
}
optm.mat <- cbind(alm.mat, x.mat)
}
Xnext <- matrix(optm.mat[which.max(optm.mat[,1]), 2:(ncol(fit1$X)+1)], nrow = 1)
x.cand <- Xnext ### Anyway it will return NULL
# Xnext <- matrix(optm.mat[ncol(fit1$X) + 1, which.max(optm.mat[(1:ncol(fit1$X)), ])], nrow = 1)
# x.cand <- matrix(optm.mat[ncol(fit1$X) + 1, ], ncol = ncol(fit1$X))
ALM <- optm.mat[1,]/c(cost[1] + cost[2] + cost[3])
} else {
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(-obj.ALM_level_3(fit = fit, newx)))
}
} else {
optm.mat <- c(rep(0, nrow(Xcand)))
for (i in 1:nrow(Xcand)) {
print(paste(i, nrow(Xcand), sep = "/"))
newx <- matrix(Xcand[i, ], nrow = 1)
optm.mat[i] <- -obj.ALM_level_3(fit = fit, newx)
}
}
Xnext <- matrix(Xcand[which.max(optm.mat), ], nrow = 1)
x.cand <- Xcand
ALM <- optm.mat/c(cost[1] + cost[2] + cost[3])
}
t2 <-proc.time()
if(trace) cat("Finding the next point:", (t2 - t1)[3], "seconds\n")
### Calculate the deduced variance ###
ALMC.1 <- obj.ALC_level_3_1(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t3 <-proc.time()
if(trace) cat("Calculating deduced variance for level 1:", (t3 - t2)[3], "seconds\n")
ALMC.2 <- obj.ALC_level_3_2(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t4 <-proc.time()
if(trace) cat("Calculating deduced variance for level 2:", (t4 - t3)[3], "seconds\n")
ALMC.3 <- obj.ALC_level_3_3(Xnext, Xref = Xref, fit = fit, mc.sample = mc.sample, parallel = parallel, ncore = ncore)
t5 <-proc.time()
if(trace) cat("Calculating deduced variance for level 3:", (t5 - t4)[3], "seconds\n")
### Compute the deduced variance at all candidate points ###
if (optim) {
intvar1 <- NULL
intvar2 <- NULL
intvar3 <- NULL
} else {
if (parallel) {
pseudointvar <- foreach(i = 1:nrow(x.cand), .combine = cbind) %dopar% {
newx <- matrix(x.cand[i, ], nrow = 1)
return(c(
obj.ALC_level_3_1(newx, Xref, fit, mc.sample),
obj.ALC_level_3_2(newx, Xref, fit, mc.sample),
obj.ALC_level_3_3(newx, Xref, fit, mc.sample)
))
}
intvar1 <- as.matrix(pseudointvar)[1, ]
intvar2 <- as.matrix(pseudointvar)[2, ]
intvar3 <- as.matrix(pseudointvar)[3, ]
} else {
intvar1 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
intvar2 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
intvar3 <- c(rep(0, nrow(x.cand))) # IMSPE candidates
for (i in 1:nrow(x.cand)) {
print(paste(i, nrow(x.cand), sep = "/"))
newx <- matrix(x.cand[i, ], nrow = 1)
intvar1[i] <- obj.ALC_level_3_1(newx, Xref, fit, mc.sample)
intvar2[i] <- obj.ALC_level_3_2(newx, Xref, fit, mc.sample)
intvar3[i] <- obj.ALC_level_3_3(newx, Xref, fit, mc.sample)
}
}
}
t6 <-proc.time()
if(trace) cat("Computing deduced variance for all points:", (t6 - t5)[3], "seconds\n")
ALMCvalue <- c(Icurrent - ALMC.1, Icurrent - ALMC.2, Icurrent - ALMC.3) / c(cost[1], cost[1] + cost[2], cost[1] + cost[2] + cost[3])
if (ALMCvalue[3] > ALMCvalue[2]) {
level <- 3
} else if (ALMCvalue[2] > ALMCvalue[1]) {
level <- 2
} else {
level <- 1
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALC <- list(ALC1 = intvar1 / cost[1], ALC2 = intvar2 / (cost[1] + cost[2]), ALC3 = intvar3 / (cost[1] + cost[2] + cost[3]))
ALMC <- c(ALMC1 = ALMC.1, ALMC2 = ALMC.2, ALMC3 = ALMC.3)
} else {
stop("level is missing")
}
if (parallel) stopImplicitCluster()
if (optim) ALC <- ALM <- NULL
return(list(ALMC = ALMC, ALM = ALM, ALC=ALC,
cost = cost, Xcand = Xcand, chosen = chosen, time = (proc.time() - t0)[3]))
}
Try the RNAmf package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.