Nothing
R/ALM.R
In RNAmf: Recursive Non-Additive Emulator for Multi-Fidelity Data
Defines functions
ALM_RNAmf
obj.ALM_level_3
obj.ALM_level_2
obj.ALM_level_1
Documented in
ALM_RNAmf
#' object to optimize the point by ALM criterion updating at level 1 with two levels of fidelity
#'
#' @param Xcand candidate data point to be optimized.
#' @param fit an object of class RNAmf.
#' @return A negative predictive posterior variance at Xcand.
#' @noRd
#'
obj.ALM_level_1 <- function(fit, Xcand) {
newx <- matrix(Xcand, nrow = 1)
fit1 <- fit$fits[[1]]
kernel <- fit$kernel
### calculate the posterior predictive variance ###
if (kernel == "sqex") {
predsig2 <- pred.GP(fit1, newx)$sig2
} else if (kernel == "matern1.5") {
predsig2 <- pred.matGP(fit1, newx)$sig2
} else if (kernel == "matern2.5") {
predsig2 <- pred.matGP(fit1, newx)$sig2
}
-predsig2 # to maximize the current variance.
}
#' object to optimize the point by ALM criterion updating at level 2 with two levels of fidelity
#'
#' @param Xcand candidate data point to be optimized.
#' @param fit an object of class RNAmf.
#' @return A negative predictive posterior variance at Xcand.
#' @noRd
#'
obj.ALM_level_2 <- function(fit, Xcand) {
newx <- matrix(Xcand, nrow = 1)
-predict(fit, newx)$sig2[[2]] # to maximize the current variance.
}
#' object to optimize the point by ALM criterion updating at level 3 with three levels of fidelity
#'
#' @param Xcand candidate data point to be optimized.
#' @param fit an object of class RNAmf.
#' @return A predictive posterior variance at Xcand.
#' @noRd
#'
obj.ALM_level_3 <- function(fit, Xcand) {
newx <- matrix(Xcand, nrow = 1)
-predict(fit, newx)$sig2[[3]] # to maximize the current variance.
}
#' @title find the next point by ALM criterion
#'
#' @description The function acquires the new point by the Active learning MacKay (ALM) criterion.
#' It calculates the ALM criterion \eqn{\frac{\sigma^{*2}_l(\bm{x})}{\sum^l_{j=1}C_j}},
#' where \eqn{\sigma^{*2}_l(\bm{x})} is the posterior predictive variance
#' at each fidelity level \eqn{l} and \eqn{C_j} is the simulation cost at level \eqn{j}.
#' 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}, a new point is acquired on unit hypercube \eqn{[0,1]^d}.
#'
#' @param Xcand vector or matrix of candidate set which could be added into the current design only used when \code{optim=FALSE}. \code{Xcand} is the set of the points where ALM criterion is evaluated. If \code{Xcand=NULL}, \eqn{100 \times d} number of points from 0 to 1 are generated by Latin hypercube design. Default is \code{NULL}.
#' @param fit object of class \code{RNAmf}.
#' @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{FALSE}.
#' @return
#' \itemize{
#' \item \code{ALM}: list of ALM criterion computed at each point of \code{Xcand} at each level if \code{optim=FALSE}. If \code{optim=TRUE}, \code{ALM} 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 lhs maximinLHS
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @importFrom doParallel registerDoParallel
#' @importFrom doParallel stopImplicitCluster
#' @usage ALM_RNAmf(Xcand = NULL, fit, 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 ###
#' alm.RNAmf.optim <- ALM_RNAmf(
#' Xcand = x, fit.RNAmf, cost = cost,
#' optim = TRUE, parallel = TRUE, ncore = 2
#' )
#' print(alm.RNAmf.optim$time) # computation time of optim=TRUE
#'
#' ### active learning with optim=FALSE ###
#' alm.RNAmf <- ALM_RNAmf(
#' Xcand = x, fit.RNAmf, cost = cost,
#' optim = FALSE, parallel = TRUE, ncore = 2
#' )
#' print(alm.RNAmf$time) # computation time of optim=FALSE
#'
#' ### visualize ALM ###
#' oldpar <- par(mfrow = c(1, 2))
#' plot(x, alm.RNAmf$ALM$ALM1,
#' type = "l", lty = 2,
#' xlab = "x", ylab = "ALM criterion at the low-fidelity level",
#' ylim = c(min(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)),
#' max(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)))
#' )
#' points(alm.RNAmf$chosen$Xnext,
#' alm.RNAmf$ALM$ALM1[which(x == drop(alm.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' plot(x, alm.RNAmf$ALM$ALM2,
#' type = "l", lty = 2,
#' xlab = "x", ylab = "ALM criterion at the high-fidelity level",
#' ylim = c(min(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)),
#' max(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)))
#' )
#' par(oldpar)}
#'
ALM_RNAmf <- function(Xcand = NULL, fit, 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")
### ALM ###
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 (parallel) registerDoParallel(ncore)
fit1 <- fit$fits[[1]]
fit2 <- 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 <- randomLHS(100*ncol(fit1$X), ncol(fit1$X))
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = 1)
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Calculate the current variance at each level ###
t1 <- proc.time()
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(
-obj.ALM_level_1(fit = fit, newx),
-obj.ALM_level_2(fit = fit, newx)
))
}
} else {
optm.mat <- rbind(c(rep(0, nrow(Xcand))), 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[1, i] <- -obj.ALM_level_1(fit = fit, newx)
optm.mat[2, i] <- -obj.ALM_level_2(fit = fit, newx)
}
}
t2 <- proc.time()
if(trace) cat("Calculating the current variance at each level:", (t2 - t1)[3], "seconds\n")
### Find the next point ###
if (optim) {
X.start <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_1, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.1 <- optim.out$par
ALM.1 <- -optim.out$value
t3 <- proc.time()
if(trace) cat("Running optim for level 1:", (t3 - t2)[3], "seconds\n")
X.start <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_2, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.2 <- optim.out$par
ALM.2 <- -optim.out$value
t4 <- proc.time()
if(trace) cat("Running optim for level 2:", (t4 - t3)[3], "seconds\n")
ALMvalue <- c(ALM.1, ALM.2) / c(cost[1], cost[1] + cost[2])
if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- Xnext.2
} else {
level <- 1
Xnext <- Xnext.1
}
} else {
ALMvalue <- c(max(optm.mat[1, ]), max(optm.mat[2, ])) / c(cost[1], cost[1] + cost[2])
if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
} else {
level <- 1
Xnext <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
}
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALM <- list(ALM1 = optm.mat[1, ] / cost[1], ALM2 = optm.mat[2, ] / (cost[1] + cost[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 (parallel) registerDoParallel(ncore)
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 <- randomLHS(100*ncol(fit1$X), ncol(fit1$X))
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = 1)
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Calculate the current variance at each level ###
t1 <- proc.time()
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(
-obj.ALM_level_1(fit = fit, newx),
-obj.ALM_level_2(fit = fit, newx),
-obj.ALM_level_3(fit = fit, newx)
))
}
} else {
optm.mat <- rbind(c(rep(0, nrow(Xcand))), c(rep(0, nrow(Xcand))), 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[1, i] <- -obj.ALM_level_1(fit = fit, newx)
optm.mat[2, i] <- -obj.ALM_level_2(fit = fit, newx)
optm.mat[3, i] <- -obj.ALM_level_3(fit = fit, newx)
}
}
t2 <- proc.time()
if(trace) cat("Calculating the current variance at each level:", (t2 - t1)[3], "seconds\n")
### Find the next point ###
if (optim) {
X.start <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_1, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.1 <- optim.out$par
ALM.1 <- -optim.out$value
t3 <- proc.time()
if(trace) cat("Running optim for level 1:", (t3 - t2)[3], "seconds\n")
X.start <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_2, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.2 <- optim.out$par
ALM.2 <- -optim.out$value
t4 <- proc.time()
if(trace) cat("Running optim for level 2:", (t4 - t3)[3], "seconds\n")
X.start <- matrix(Xcand[which.max(optm.mat[3, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.3 <- optim.out$par
ALM.3 <- -optim.out$value
t5 <- proc.time()
if(trace) cat("Running optim for level 3:", (t5 - t4)[3], "seconds\n")
ALMvalue <- c(ALM.1, ALM.2, ALM.3) / c(cost[1], cost[1] + cost[2], cost[1] + cost[2] + cost[3])
if (ALMvalue[3] > ALMvalue[2]) {
level <- 3
Xnext <- Xnext.3
} else if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- Xnext.2
} else {
level <- 1
Xnext <- Xnext.1
}
} else {
ALMvalue <- c(max(optm.mat[1, ]), max(optm.mat[2, ]), max(optm.mat[3, ])) / c(cost[1], cost[1] + cost[2], cost[1] + cost[2] + cost[3])
if (ALMvalue[3] > ALMvalue[2]) {
level <- 3
Xnext <- matrix(Xcand[which.max(optm.mat[3, ]), ], nrow = 1)
} else if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
} else {
level <- 1
Xnext <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
}
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALM <- list(ALM1 = optm.mat[1, ] / cost[1], ALM2 = optm.mat[2, ] / (cost[1] + cost[2]), ALM3 = optm.mat[3, ] / (cost[1] + cost[2] + cost[3]))
} else {
stop("level is missing")
}
if (parallel) stopImplicitCluster()
if (optim) ALM <- NULL
return(list(ALM = ALM, 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 ALM_RNAmf obj.ALM_level_3 obj.ALM_level_2 obj.ALM_level_1
Documented in ALM_RNAmf
#' object to optimize the point by ALM criterion updating at level 1 with two levels of fidelity
#'
#' @param Xcand candidate data point to be optimized.
#' @param fit an object of class RNAmf.
#' @return A negative predictive posterior variance at Xcand.
#' @noRd
#'
obj.ALM_level_1 <- function(fit, Xcand) {
newx <- matrix(Xcand, nrow = 1)
fit1 <- fit$fits[[1]]
kernel <- fit$kernel
### calculate the posterior predictive variance ###
if (kernel == "sqex") {
predsig2 <- pred.GP(fit1, newx)$sig2
} else if (kernel == "matern1.5") {
predsig2 <- pred.matGP(fit1, newx)$sig2
} else if (kernel == "matern2.5") {
predsig2 <- pred.matGP(fit1, newx)$sig2
}
-predsig2 # to maximize the current variance.
}
#' object to optimize the point by ALM criterion updating at level 2 with two levels of fidelity
#'
#' @param Xcand candidate data point to be optimized.
#' @param fit an object of class RNAmf.
#' @return A negative predictive posterior variance at Xcand.
#' @noRd
#'
obj.ALM_level_2 <- function(fit, Xcand) {
newx <- matrix(Xcand, nrow = 1)
-predict(fit, newx)$sig2[[2]] # to maximize the current variance.
}
#' object to optimize the point by ALM criterion updating at level 3 with three levels of fidelity
#'
#' @param Xcand candidate data point to be optimized.
#' @param fit an object of class RNAmf.
#' @return A predictive posterior variance at Xcand.
#' @noRd
#'
obj.ALM_level_3 <- function(fit, Xcand) {
newx <- matrix(Xcand, nrow = 1)
-predict(fit, newx)$sig2[[3]] # to maximize the current variance.
}
#' @title find the next point by ALM criterion
#'
#' @description The function acquires the new point by the Active learning MacKay (ALM) criterion.
#' It calculates the ALM criterion \eqn{\frac{\sigma^{*2}_l(\bm{x})}{\sum^l_{j=1}C_j}},
#' where \eqn{\sigma^{*2}_l(\bm{x})} is the posterior predictive variance
#' at each fidelity level \eqn{l} and \eqn{C_j} is the simulation cost at level \eqn{j}.
#' 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}, a new point is acquired on unit hypercube \eqn{[0,1]^d}.
#'
#' @param Xcand vector or matrix of candidate set which could be added into the current design only used when \code{optim=FALSE}. \code{Xcand} is the set of the points where ALM criterion is evaluated. If \code{Xcand=NULL}, \eqn{100 \times d} number of points from 0 to 1 are generated by Latin hypercube design. Default is \code{NULL}.
#' @param fit object of class \code{RNAmf}.
#' @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{FALSE}.
#' @return
#' \itemize{
#' \item \code{ALM}: list of ALM criterion computed at each point of \code{Xcand} at each level if \code{optim=FALSE}. If \code{optim=TRUE}, \code{ALM} 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 lhs maximinLHS
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @importFrom doParallel registerDoParallel
#' @importFrom doParallel stopImplicitCluster
#' @usage ALM_RNAmf(Xcand = NULL, fit, 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 ###
#' alm.RNAmf.optim <- ALM_RNAmf(
#' Xcand = x, fit.RNAmf, cost = cost,
#' optim = TRUE, parallel = TRUE, ncore = 2
#' )
#' print(alm.RNAmf.optim$time) # computation time of optim=TRUE
#'
#' ### active learning with optim=FALSE ###
#' alm.RNAmf <- ALM_RNAmf(
#' Xcand = x, fit.RNAmf, cost = cost,
#' optim = FALSE, parallel = TRUE, ncore = 2
#' )
#' print(alm.RNAmf$time) # computation time of optim=FALSE
#'
#' ### visualize ALM ###
#' oldpar <- par(mfrow = c(1, 2))
#' plot(x, alm.RNAmf$ALM$ALM1,
#' type = "l", lty = 2,
#' xlab = "x", ylab = "ALM criterion at the low-fidelity level",
#' ylim = c(min(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)),
#' max(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)))
#' )
#' points(alm.RNAmf$chosen$Xnext,
#' alm.RNAmf$ALM$ALM1[which(x == drop(alm.RNAmf$chosen$Xnext))],
#' pch = 16, cex = 1, col = "red"
#' )
#' plot(x, alm.RNAmf$ALM$ALM2,
#' type = "l", lty = 2,
#' xlab = "x", ylab = "ALM criterion at the high-fidelity level",
#' ylim = c(min(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)),
#' max(c(alm.RNAmf$ALM$ALM1, alm.RNAmf$ALM$ALM2)))
#' )
#' par(oldpar)}
#'
ALM_RNAmf <- function(Xcand = NULL, fit, 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")
### ALM ###
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 (parallel) registerDoParallel(ncore)
fit1 <- fit$fits[[1]]
fit2 <- 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 <- randomLHS(100*ncol(fit1$X), ncol(fit1$X))
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = 1)
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Calculate the current variance at each level ###
t1 <- proc.time()
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(
-obj.ALM_level_1(fit = fit, newx),
-obj.ALM_level_2(fit = fit, newx)
))
}
} else {
optm.mat <- rbind(c(rep(0, nrow(Xcand))), 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[1, i] <- -obj.ALM_level_1(fit = fit, newx)
optm.mat[2, i] <- -obj.ALM_level_2(fit = fit, newx)
}
}
t2 <- proc.time()
if(trace) cat("Calculating the current variance at each level:", (t2 - t1)[3], "seconds\n")
### Find the next point ###
if (optim) {
X.start <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_1, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.1 <- optim.out$par
ALM.1 <- -optim.out$value
t3 <- proc.time()
if(trace) cat("Running optim for level 1:", (t3 - t2)[3], "seconds\n")
X.start <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_2, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.2 <- optim.out$par
ALM.2 <- -optim.out$value
t4 <- proc.time()
if(trace) cat("Running optim for level 2:", (t4 - t3)[3], "seconds\n")
ALMvalue <- c(ALM.1, ALM.2) / c(cost[1], cost[1] + cost[2])
if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- Xnext.2
} else {
level <- 1
Xnext <- Xnext.1
}
} else {
ALMvalue <- c(max(optm.mat[1, ]), max(optm.mat[2, ])) / c(cost[1], cost[1] + cost[2])
if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
} else {
level <- 1
Xnext <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
}
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALM <- list(ALM1 = optm.mat[1, ] / cost[1], ALM2 = optm.mat[2, ] / (cost[1] + cost[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 (parallel) registerDoParallel(ncore)
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 <- randomLHS(100*ncol(fit1$X), ncol(fit1$X))
}else if(is.null(dim(Xcand))){
Xcand <- matrix(Xcand, ncol = 1)
}
}
# if (ncol(Xcand) != dim(fit$fit1$X)[2]) stop("The dimension of candidate set should be equal to the dimension of the design")
### Calculate the current variance at each level ###
t1 <- proc.time()
if (parallel) {
optm.mat <- foreach(i = 1:nrow(Xcand), .combine = cbind) %dopar% {
newx <- matrix(Xcand[i, ], nrow = 1)
return(c(
-obj.ALM_level_1(fit = fit, newx),
-obj.ALM_level_2(fit = fit, newx),
-obj.ALM_level_3(fit = fit, newx)
))
}
} else {
optm.mat <- rbind(c(rep(0, nrow(Xcand))), c(rep(0, nrow(Xcand))), 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[1, i] <- -obj.ALM_level_1(fit = fit, newx)
optm.mat[2, i] <- -obj.ALM_level_2(fit = fit, newx)
optm.mat[3, i] <- -obj.ALM_level_3(fit = fit, newx)
}
}
t2 <- proc.time()
if(trace) cat("Calculating the current variance at each level:", (t2 - t1)[3], "seconds\n")
### Find the next point ###
if (optim) {
X.start <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_1, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.1 <- optim.out$par
ALM.1 <- -optim.out$value
t3 <- proc.time()
if(trace) cat("Running optim for level 1:", (t3 - t2)[3], "seconds\n")
X.start <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_2, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.2 <- optim.out$par
ALM.2 <- -optim.out$value
t4 <- proc.time()
if(trace) cat("Running optim for level 2:", (t4 - t3)[3], "seconds\n")
X.start <- matrix(Xcand[which.max(optm.mat[3, ]), ], nrow = 1)
optim.out <- optim(X.start, obj.ALM_level_3, method = "L-BFGS-B", lower = 0, upper = 1, fit = fit)
Xnext.3 <- optim.out$par
ALM.3 <- -optim.out$value
t5 <- proc.time()
if(trace) cat("Running optim for level 3:", (t5 - t4)[3], "seconds\n")
ALMvalue <- c(ALM.1, ALM.2, ALM.3) / c(cost[1], cost[1] + cost[2], cost[1] + cost[2] + cost[3])
if (ALMvalue[3] > ALMvalue[2]) {
level <- 3
Xnext <- Xnext.3
} else if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- Xnext.2
} else {
level <- 1
Xnext <- Xnext.1
}
} else {
ALMvalue <- c(max(optm.mat[1, ]), max(optm.mat[2, ]), max(optm.mat[3, ])) / c(cost[1], cost[1] + cost[2], cost[1] + cost[2] + cost[3])
if (ALMvalue[3] > ALMvalue[2]) {
level <- 3
Xnext <- matrix(Xcand[which.max(optm.mat[3, ]), ], nrow = 1)
} else if (ALMvalue[2] > ALMvalue[1]) {
level <- 2
Xnext <- matrix(Xcand[which.max(optm.mat[2, ]), ], nrow = 1)
} else {
level <- 1
Xnext <- matrix(Xcand[which.max(optm.mat[1, ]), ], nrow = 1)
}
}
chosen <- list(
"level" = level, # next level
"Xnext" = Xnext
) # next point
ALM <- list(ALM1 = optm.mat[1, ] / cost[1], ALM2 = optm.mat[2, ] / (cost[1] + cost[2]), ALM3 = optm.mat[3, ] / (cost[1] + cost[2] + cost[3]))
} else {
stop("level is missing")
}
if (parallel) stopImplicitCluster()
if (optim) ALM <- NULL
return(list(ALM = ALM, 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.