diff --git a/DESCRIPTION b/DESCRIPTION index 0de5635..497dd53 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,14 +1,14 @@ Package: EPBO Title: Bayesian Optimization via Exact Penalty -Version: 0.0.11 -Date: 2023-11-20 +Version: 0.1.0 +Date: 2023-12-14 Authors@R: person(given = "Jiangyan", family = "Zhao", email = "zhaojy2017@126.com", role = c("aut", "cre")) -Description: Wrapper routines for blackbox optimization under mixed equality and inequality - constraints via an exact penalty. +Description: Wrapper routines for blackbox optimization under mixed (equality and/or inequality) + constraints via exact penalty. License: AGPL (>= 3) URL: https://github.com/Jiangyan-Zhao/EPBO BugReports: https://github.com/Jiangyan-Zhao/EPBO/issues @@ -16,7 +16,7 @@ Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) RoxygenNote: 7.2.3 -Imports: laGP, tgp, DiceKriging, matrixStats, mvtnorm +Imports: laGP, tgp Depends: R (>= 2.14) Suggests: stats, diff --git a/NAMESPACE b/NAMESPACE index 2870745..690fbbc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,38 +1,18 @@ # Generated by roxygen2: do not edit by hand -export(AF_AE) -export(AF_ConVar) -export(AF_EY) -export(AF_EY_Kernel) -export(AF_LCB) -export(AF_OOSS) -export(AF_ScaledEI) -export(AF_ScaledEI_Kernel) -export(AF_ScaledEI_MC) -export(AF_TS) -export(AF_UEI) -export(AF_UEI_Kernel) -export(AF_UEI_MC) -export(AF_VEI) -export(bilog) -export(normalize) export(optim.AE) export(optim.BM) export(optim.EP) -export(optim.EP.MC) -export(optim.EP.kernel) -export(optim.EP4HD) -export(unbilog) -export(unnormalize) -import(DiceKriging) import(laGP) -import(matrixStats) import(tgp) +importFrom(graphics,image) +importFrom(graphics,par) +importFrom(graphics,points) importFrom(stats,dnorm) +importFrom(stats,median) importFrom(stats,optim) importFrom(stats,pnorm) importFrom(stats,quantile) -importFrom(stats,rnorm) importFrom(stats,runif) importFrom(stats,sd) importFrom(utils,tail) diff --git a/R/AE.R b/R/AE.R index be7e795..4207695 100644 --- a/R/AE.R +++ b/R/AE.R @@ -46,7 +46,7 @@ #' \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} #' \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } #' -#' @seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.BM}}, \code{\link[EPBO]{optim.EP}} +#' @seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.BM}}, \code{\link[EPBO]{optim.EP}} #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} #' @@ -68,10 +68,9 @@ #' #' @examples #' # search space -#' B = rbind(c(-2.25, 2.5), c(-2.25, 1.75)) +#' B = rbind(c(-2.25, 2.5), c(-2.5, 1.75)) #' # Objective and constraint function for use -#' MTP = function(x, B=rbind(c(-2.25, 2.5), c(-2.5, 1.75))){ -#' x = pmin(pmax(B[,1], x), B[,2]) +#' MTP = function(x){ #' f = -(cos((x[1]-0.1)*x[2]))^2 - x[1]*sin(3*x[1]+x[2]) #' t = atan2(x[1], x[2]) #' c = x[1]^2 + x[2]^2 -((2*cos(t)-1/2*cos(2*t)-1/4*cos(3*t)-1/8*cos(4*t))^2) - ((2*sin(t))^2) @@ -84,14 +83,13 @@ #' AE$xbest -optim.AE = function(blackbox, B, - start=10, end=100, - Xstart=NULL, urate=10, - ncandf=function(k) { k }, - alpha1=1, alpha2=5, omega=2/3, - dg.start=c(0.1,1e-6), ab=c(3/2,8), - dlim=sqrt(ncol(B))*c(1/100,10), - verb=2, ...){ +optim.AE = function( + blackbox, B, start=10, end=100, + Xstart=NULL, urate=10, ncandf=function(k) { k }, + alpha1=1, alpha2=5, omega=2/3, + dg.start=c(0.1,1e-6), ab=c(3/2,8), dlim=sqrt(ncol(B))*c(1/100,10), + verb=2, ...) +{ ## check start if(start >= end) stop("must have start < end") ## get initial design diff --git a/R/AF_AE.R b/R/AF_AE.R index 97b76aa..e5f15f2 100644 --- a/R/AF_AE.R +++ b/R/AF_AE.R @@ -1,26 +1,21 @@ -#' @title Asymmetric Entropy +#' @title Asymmetric entropy acquisition function #' -#' @description Asymmetric Entropy +#' @description The asymmetric entropy (AE) acquisition function of the AE method #' -#' @param x description +#' @param x a vector containing a single candidate point; or a \code{matrix} with +#' multiple candidate points +#' @param fgpi the GP surrogate model of the objective function +#' @param fnorm the maximum of the objective +#' @param Cgpi the GP surrogate models of the constraints +#' @param Cnorm the maxima of the constraints +#' @param fmin the best objective value obtained so far +#' @param alpha1 a specified weight for the EI part +#' @param alpha2 a specified weight for the AE part +#' @param omega a mode location parameter #' -#' @param fgpi description +#' @returns The AE value(s) at \code{x}. #' -#' @param fnorm description -#' -#' @param Cgpi description -#' -#' @param Cnorm description -#' -#' @param fmin description -#' -#' @param alpha1 description -#' -#' @param alpha2 description -#' -#' @param omega description -#' -#' @returns AF +#' @seealso \code{\link[EPBO]{AF_ScaledEI}}, \code{\link[EPBO]{AF_EY}}, \code{\link[EPBO]{AF_OOSS}} #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} #' @@ -31,15 +26,7 @@ #' @importFrom stats dnorm #' @importFrom stats pnorm #' @importFrom stats quantile -#' -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' + AF_AE = function(x, fgpi, fnorm, Cgpi, Cnorm, fmin, alpha1=1, alpha2=5, omega=2/3) diff --git a/R/AF_ConVar.R b/R/AF_ConVar.R deleted file mode 100644 index 7c748f4..0000000 --- a/R/AF_ConVar.R +++ /dev/null @@ -1,66 +0,0 @@ -#' @title Constained variance acquisition function -#' -#' @description Predictive mean -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_ConVar = function(x, fgpi, Cgpi, equal) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x, lite=TRUE) - sigma_f = sqrt(pred_f$s2) - - ## constraints - nc = length(Cgpi) # number of the constraint - LCB = matrix(NA, ncand, nc) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x, lite=TRUE) - mu_C = pred_C$mean - sigma_C = sqrt(pred_C$s2) - if(equal[j]){ - LCB[,j] = abs(mu_C) - 6*sigma_C - }else{ - LCB[,j] = mu_C - 6*sigma_C - } - } - - ## Acquaisition function - infeasible = apply(LCB > 0, 1, any) - if(all(infeasible)){ - AF = sigma_f - }else{ - sigma_f[infeasible] = 0 - AF = sigma_f - } - - return(AF) -} diff --git a/R/AF_EY.R b/R/AF_EY.R index a99e978..d819c0f 100644 --- a/R/AF_EY.R +++ b/R/AF_EY.R @@ -1,40 +1,28 @@ #' @title Predictive mean acquisition function #' -#' @description Predictive mean -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param fmean description -#' -#' @param fsd description -#' -#' @param Cgpi description -#' -#' @param rho description -#' +#' @description The predictive mean acquisition function of the EPBO method +#' +#' @param x a vector containing a single candidate point; or a \code{matrix} with +#' multiple candidate points +#' @param fgpi the GP surrogate model of the objective function +#' @param fmean the mean of the objective value +#' @param fsd the standard deviation of the objective value +#' @param Cgpi the GP surrogate models of the constraints +#' @param rho the penalty parameters #' @param equal an optional vector containing zeros and ones, whose length equals the number of #' constraints, specifying which should be treated as equality constraints (\code{1}) and #' which as inequality (\code{0}) #' +#' @returns The Predictive Mean at \code{x}. #' -#' @returns AF +#' @seealso \code{\link[EPBO]{AF_ScaledEI}}, \code{\link[EPBO]{AF_OOSS}}, \code{\link[EPBO]{AF_AE}} #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} #' -#' #' @import laGP #' @importFrom stats dnorm #' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' + AF_EY = function(x, fgpi, fmean, fsd, Cgpi, rho, equal) { diff --git a/R/AF_EY_Kernel.R b/R/AF_EY_Kernel.R deleted file mode 100644 index daadd2f..0000000 --- a/R/AF_EY_Kernel.R +++ /dev/null @@ -1,59 +0,0 @@ -#' @title Predictive mean acquisition function -#' -#' @description Predictive mean -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_EY_Kernel = function(x, fgpi, fmean, fsd, Cgpi, rho, equal, type="UK") -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predict(object=fgpi, newdata=x, type=type, checkNames = FALSE, light.return = TRUE) - mu_f = pred_f$mean * fsd + fmean - - ## constraints - nc = length(Cgpi) # number of the constraint - EV = matrix(NA, nc, nrow(x)) - for (j in 1:nc) { - pred_C = predict(object=Cgpi[[j]], newdata=x, type=type, checkNames = FALSE, light.return = TRUE) - mu_C = pred_C$mean - sigma_C = pred_C$sd - dC = mu_C/sigma_C - EV[j,] = mu_C*((equal[j]+1)*pnorm(dC) - equal[j]) + (equal[j]+1)*sigma_C*dnorm(dC) - } - - ## Acquaisition function - AF = mu_f + rho%*%EV - - return(AF) -} diff --git a/R/AF_LCB.R b/R/AF_LCB.R deleted file mode 100644 index 4a96ee3..0000000 --- a/R/AF_LCB.R +++ /dev/null @@ -1,64 +0,0 @@ - -#' @title ScaledEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_LCB = function(x, fgpi, fmean, fsd, Cgpi, rho, equal) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x, lite=TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = sqrt(pred_f$s2) * fsd - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = omega = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x, lite=TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = sqrt(pred_C$s2) - omega[j,] = (equal[j]+1)*pnorm(mu_C[j,]/sigma_C[j,]) - equal[j] - } - - ## Acquaisition function - mu_ep = mu_f + rho%*%(omega*mu_C) - sigma_ep = sqrt(sigma_f^2 + (rho^2)%*%((omega*sigma_C)^2)) - LCB = mu_ep - 3*sigma_ep - return(LCB) -} - diff --git a/R/AF_OOSS.R b/R/AF_OOSS.R index 36ebf05..b6ef99e 100644 --- a/R/AF_OOSS.R +++ b/R/AF_OOSS.R @@ -1,19 +1,17 @@ -#' @title Asymmetric Entropy +#' @title One over sigma squared (OOSS) acquisition function #' -#' @description Asymmetric Entropy +#' @description The one over sigma squared (OOSS) acquisition function of the barrier method #' -#' @param x description +#' @param x a vector containing a single candidate point; or a \code{matrix} with +#' multiple candidate points +#' @param fgpi the GP surrogate model of the objective function +#' @param fnorm the maximum of the objective +#' @param Cgpi the GP surrogate models of the constraints +#' @param Cnorm the maxima of the constraints #' -#' @param fgpi description +#' @returns The OOSS value(s) at \code{x}. #' -#' @param fnorm description -#' -#' @param Cgpi description -#' -#' @param Cnorm description -#' -#' -#' @returns AF +#' @seealso \code{\link[EPBO]{AF_ScaledEI}}, \code{\link[EPBO]{AF_EY}}, \code{\link[EPBO]{AF_AE}} #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} #' @@ -21,15 +19,7 @@ #' \emph{Journal of Computational and Graphical Statistics} 31(1), 74-83. #' #' @import laGP -#' -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' + AF_OOSS = function(x, fgpi, fnorm, Cgpi, Cnorm) { diff --git a/R/AF_ScaledEI.R b/R/AF_ScaledEI.R index e99dc70..a74d51d 100644 --- a/R/AF_ScaledEI.R +++ b/R/AF_ScaledEI.R @@ -1,27 +1,22 @@ #' @title ScaledEI acquisition function #' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param fmean description -#' -#' @param fsd description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' +#' @description The Scaled expected improvement acquisition function of the EPBO method +#' +#' @param x a vector containing a single candidate point; or a \code{matrix} with +#' multiple candidate points +#' @param fgpi the GP surrogate model of the objective function +#' @param fmean the mean of the objective value +#' @param fsd the standard deviation of the objective value +#' @param Cgpi the GP surrogate models of the constraints +#' @param epbest the best exact penalty value obtained so far +#' @param rho the penalty parameters #' @param equal an optional vector containing zeros and ones, whose length equals the number of #' constraints, specifying which should be treated as equality constraints (\code{1}) and #' which as inequality (\code{0}) #' +#' @returns The ScaledEI value(s) at \code{x}. #' -#' @returns AF +#' @seealso \code{\link[EPBO]{AF_OOSS}}, \code{\link[EPBO]{AF_EY}}, \code{\link[EPBO]{AF_AE}} #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} #' @@ -31,14 +26,6 @@ #' @import laGP #' @importFrom stats dnorm #' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' AF_ScaledEI = function(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) { diff --git a/R/AF_ScaledEI_Kernel.R b/R/AF_ScaledEI_Kernel.R deleted file mode 100644 index d2ac978..0000000 --- a/R/AF_ScaledEI_Kernel.R +++ /dev/null @@ -1,73 +0,0 @@ -#' @title ScaledEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param fmean description -#' -#' @param fsd description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_ScaledEI_Kernel = function(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, type="UK") -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predict(object=fgpi, newdata=x, type=type, checkNames = FALSE, light.return = TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = pred_f$sd * fsd - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = omega = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predict(object=Cgpi[[j]], newdata=x, type=type, checkNames = FALSE, light.return = TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = pred_C$sd - omega[j,] = (equal[j]+1)*pnorm(mu_C[j,]/sigma_C[j,]) - equal[j] - } - - ## Acquaisition function - mu_ep = mu_f + rho%*%(omega*mu_C) - sigma_ep = sqrt(sigma_f^2 + (rho^2)%*%((omega*sigma_C)^2)) - d = (epbest - mu_ep)/sigma_ep - EI = d*pnorm(d) + dnorm(d) # expected improvement (remove sigma_ep) - VI = pnorm(d) + d*EI - EI^2 # variance of the improvement (remove sigma_ep) - VI = pmax(.Machine$double.xmin, VI) - ScaledEI = EI/sqrt(VI) # Scaled expected improvement - ScaledEI[is.nan(ScaledEI)] = 0 - return(ScaledEI) -} diff --git a/R/AF_ScaledEI_MC.R b/R/AF_ScaledEI_MC.R deleted file mode 100644 index 363b7a5..0000000 --- a/R/AF_ScaledEI_MC.R +++ /dev/null @@ -1,81 +0,0 @@ -#' @title ScaledEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param fmean description -#' -#' @param fsd description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' @param N positive scalar integer indicating the number of Monte Carlo samples to be used -#' for calculating ScaledEI -#' -#' @returns ScaledEI -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats rnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_ScaledEI_MC = function(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N=1000) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x, lite=TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = sqrt(pred_f$s2) * fsd - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x, lite=TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = sqrt(pred_C$s2) - } - - ## expected improvement - EP = matrix(NA, nrow = N, ncol = ncand) - for (n in 1:N) { - EP[n,] = rnorm(ncand, mu_f, sigma_f) - for (j in 1:nc) { - if(equal[j]){ - EP[n,] = EP[n,] + rho[j]*abs(rnorm(ncand, mu_C[j,], sigma_C[j,])) - }else{ - EP[n,] = EP[n,] + rho[j]*pmax(0, rnorm(ncand, mu_C[j,], sigma_C[j,])) - } - } - } - improvement = matrix(pmax(0, epbest-EP), nrow = N) - EI = colMeans(improvement) # expected improvement - SI = colSds(improvement) # Standard deviation of the improvement - ScaledEI = EI/SI # Scaled expected improvement - ScaledEI[SI == 0] = 0 - return(ScaledEI) -} diff --git a/R/AF_TS.R b/R/AF_TS.R deleted file mode 100644 index c2bc368..0000000 --- a/R/AF_TS.R +++ /dev/null @@ -1,73 +0,0 @@ -#' @title ScaledEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' @param batch_size aaa -#' -#' @returns ScaledEI -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats rnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_TS = function(x, fgpi, Cgpi, rho, equal, batch_size=1) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x) - mu_f = pred_f$mean - sigma_f = pred_f$Sigma - df_f = pred_f$df - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = matrix(NA, nc, ncand) - sigma_C = vector("list", length = nc) - df_C = rep(NA, nc) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x) - mu_C[j,] = pred_C$mean - sigma_C[[j]] = pred_C$Sigma - df_C[j] = pred_C$df - } - - ## expected improvement - AF = matrix(NA, batch_size, ncand) - for (b in 1:batch_size) { - AF[b,] = rmvt(1, sigma_f, df_f) + mu_f - for (j in 1:nc) { - if(equal[j]){ - AF[b,] = AF[b,] + rho[j]*abs(rmvt(1, sigma_C[[j]], df_C[j])+mu_C[j,]) - }else{ - AF[b,] = AF[b,] + rho[j]*pmax(0, rmvt(1, sigma_C[[j]], df_C[j])+mu_C[j,]) - } - } - } - return(AF) -} diff --git a/R/AF_UEI.R b/R/AF_UEI.R deleted file mode 100644 index e000d68..0000000 --- a/R/AF_UEI.R +++ /dev/null @@ -1,69 +0,0 @@ -#' @title UEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_UEI = function(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, beta=3) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x, lite=TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = sqrt(pred_f$s2) * fsd - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = omega = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x, lite=TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = sqrt(pred_C$s2) - omega[j,] = (equal[j]+1)*pnorm(mu_C[j,]/sigma_C[j,]) - equal[j] - } - - ## Acquaisition function - mu_ep = mu_f + rho%*%(omega*mu_C) - sigma_ep = sqrt(sigma_f^2 + (rho^2)%*%((omega*sigma_C)^2)) - d = (epbest - mu_ep)/sigma_ep - EI = sigma_ep * (d*pnorm(d) + dnorm(d)) # expected improvement - VI = sigma_ep^2 * ((d^2+1)*pnorm(d) + d*dnorm(d)) - EI^2 # variance of the improvement (remove sigma_ep) - VI = pmax(.Machine$double.xmin, VI) - UEI = EI + beta*sqrt(VI) - EI[is.nan(UEI)] = 0 - return(UEI) -} diff --git a/R/AF_UEI_Kernel.R b/R/AF_UEI_Kernel.R deleted file mode 100644 index aa2fdcb..0000000 --- a/R/AF_UEI_Kernel.R +++ /dev/null @@ -1,69 +0,0 @@ -#' @title UEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_UEI_Kernel = function(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, beta=3, type="UK") -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predict(object=fgpi, newdata=x, type=type, checkNames = FALSE, light.return = TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = pred_f$sd * fsd - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = omega = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predict(object=Cgpi[[j]], newdata=x, type=type, checkNames = FALSE, light.return = TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = pred_C$sd - omega[j,] = (equal[j]+1)*pnorm(mu_C[j,]/sigma_C[j,]) - equal[j] - } - - ## Acquaisition function - mu_ep = mu_f + rho%*%(omega*mu_C) - sigma_ep = sqrt(sigma_f^2 + (rho^2)%*%((omega*sigma_C)^2)) - d = (epbest - mu_ep)/sigma_ep - EI = sigma_ep * (d*pnorm(d) + dnorm(d)) # expected improvement - VI = sigma_ep^2 * ((d^2+1)*pnorm(d) + d*dnorm(d)) - EI^2 # variance of the improvement (remove sigma_ep) - VI = pmax(.Machine$double.xmin, VI) - UEI = EI + beta*sqrt(VI) - UEI[is.nan(UEI)] = 0 - return(UEI) -} diff --git a/R/AF_UEI_MC.R b/R/AF_UEI_MC.R deleted file mode 100644 index 34bb200..0000000 --- a/R/AF_UEI_MC.R +++ /dev/null @@ -1,83 +0,0 @@ -#' @title UEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param fmean description -#' -#' @param fsd description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' @param N description -#' -#' @param beta description -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_UEI_MC = function(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N=1000, beta=3) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x, lite=TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = sqrt(pred_f$s2) * fsd - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x, lite=TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = sqrt(pred_C$s2) - } - - ## expected improvement - EP = matrix(NA, nrow = N, ncol = ncand) - for (n in 1:N) { - EP[n,] = rnorm(ncand, mu_f, sigma_f) - for (j in 1:nc) { - if(equal[j]){ - EP[n,] = EP[n,] + rho[j]*abs(rnorm(ncand, mu_C[j,], sigma_C[j,])) - }else{ - EP[n,] = EP[n,] + rho[j]*pmax(0, rnorm(ncand, mu_C[j,], sigma_C[j,])) - } - } - } - improvement = matrix(pmax(0, epbest-EP), nrow = N) - EI = colMeans(improvement) # expected improvement - SI = colSds(improvement) # Standard deviation of the improvement - UEI = EI + beta * SI - EI[is.nan(UEI)] = 0 - return(UEI) -} diff --git a/R/AF_VEI.R b/R/AF_VEI.R deleted file mode 100644 index 6f46bbb..0000000 --- a/R/AF_VEI.R +++ /dev/null @@ -1,68 +0,0 @@ -#' @title VEI acquisition function -#' -#' @description Scaled expected improvement -#' -#' @param x description -#' -#' @param fgpi description -#' -#' @param Cgpi description -#' -#' @param epbest description -#' -#' @param rho description -#' -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' -#' -#' @returns AF -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -#' for Bayesian optimization. \emph{arXiv:1808.06918}. -#' -#' @import laGP -#' @importFrom stats dnorm -#' @importFrom stats pnorm -#' -#' -#' @export -#' -#' @examples -#' B = rbind(c(0, 1), c(0, 1)) -#' -#' - -AF_VEI = function(x, fgpi, Cgpi, epbest, rho, equal) -{ - if(is.null(nrow(x))) x = matrix(x, nrow=1) - ncand = nrow(x) # number of the candidate points - - ## objective - pred_f = predGPsep(fgpi, x, lite=TRUE) - mu_f = pred_f$mean - sigma_f = sqrt(pred_f$s2) - - ## constraints - nc = length(Cgpi) # number of the constraint - mu_C = sigma_C = omega = matrix(NA, nc, ncand) - for (j in 1:nc) { - pred_C = predGPsep(Cgpi[j], x, lite=TRUE) - mu_C[j,] = pred_C$mean - sigma_C[j,] = sqrt(pred_C$s2) - omega[j,] = (equal[j]+1)*pnorm(mu_C[j,]/sigma_C[j,]) - equal[j] - } - - ## Acquaisition function - mu_ep = mu_f + rho%*%(omega*mu_C) - sigma_ep = sqrt(sigma_f^2 + (rho^2)%*%((omega*sigma_C)^2)) - d = (epbest - mu_ep)/sigma_ep - EI = sigma_ep * (d*pnorm(d) + dnorm(d)) # expected improvement - VI = sigma_ep^2 * ((d^2+1)*pnorm(d) + d*dnorm(d)) - EI^2 # variance of the improvement (remove sigma_ep) - VI = pmax(0, VI) - VEI = EI - 0.5*VI # Scaled expected improvement - return(VEI) -} diff --git a/R/BM.R b/R/BM.R index 5c4473a..b25d948 100644 --- a/R/BM.R +++ b/R/BM.R @@ -65,10 +65,9 @@ #' #' @examples #' # search space -#' B = rbind(c(-2.25, 2.5), c(-2.25, 1.75)) +#' B = rbind(c(-2.25, 2.5), c(-2.5, 1.75)) #' # Objective and constraint function for use -#' MTP = function(x, B=rbind(c(-2.25, 2.5), c(-2.5, 1.75))){ -#' x = pmin(pmax(B[,1], x), B[,2]) +#' MTP = function(x){ #' f = -(cos((x[1]-0.1)*x[2]))^2 - x[1]*sin(3*x[1]+x[2]) #' t = atan2(x[1], x[2]) #' c = x[1]^2 + x[2]^2 -((2*cos(t)-1/2*cos(2*t)-1/4*cos(3*t)-1/8*cos(4*t))^2) - ((2*sin(t))^2) @@ -81,13 +80,12 @@ #' BM$xbest -optim.BM = function(blackbox, B, - start=10, end=100, - Xstart=NULL, urate=10, - ncandf=function(k) { k }, - dg.start=c(0.1,1e-6), ab=c(3/2,8), - dlim=sqrt(ncol(B))*c(1/100,10), - verb=2, ...){ +optim.BM = function( + blackbox, B, start=10, end=100, + Xstart=NULL, urate=10, ncandf=function(k) { k }, + dg.start=c(0.1,1e-6), ab=c(3/2,8), dlim=sqrt(ncol(B))*c(1/100,10), + verb=2, ...) +{ ## check start if(start >= end) stop("must have start < end") ## get initial design diff --git a/R/EPBO.R b/R/EPBO.R index 6d9c7c6..e2ec4cb 100644 --- a/R/EPBO.R +++ b/R/EPBO.R @@ -3,7 +3,7 @@ #' @description Black-box optimization under mixed equality and inequality constraints via an exact penalty. #' #' @param blackbox blackbox of an input (\code{x}), facilitating vectorization on a -#' \code{matrix} \code{X} thereof, returning a \code{list} +#' \code{matrix} \code{X} thereof, returning a \code{list} #' with elements \code{$obj} containing the (scalar) objective value and \code{$c} #' containing a vector of evaluations of the (multiple) constraint function at \code{x}. #' @param B 2-column \code{matrix} describing the bounding box. The number of rows @@ -26,21 +26,28 @@ #' @param urate positive integer indicating how many optimization trials should pass before #' each MLE/MAP update is performed for GP correlation lengthscale #' parameter(s) -#' @param rho positive vector initial exact penalty parameter in the exact penalty function; -#' the default setting of \code{rho = NULL} causes an automatic starting value to be chosen +#' @param rho positive vector initial exact penalty parameters in the exact penalty function; +#' the default setting of \code{rho = NULL} causes the automatic starting values to be chosen #' @param ncandf function taking a single integer indicating the optimization trial number \code{t}, where #' \code{start < t <= end}, and returning the number of search candidates (e.g., for #' expected improvement calculations) at round \code{t}; the default setting #' allows the number of candidates to grow linearly with \code{t} #' @param dg_start 2-vector giving starting values for the lengthscale and nugget parameters -#' of the GP surrogate model(s) for constraints +#' of the GP surrogate model(s) for the objective and constraints #' @param dlim 2-vector giving bounds for the lengthscale parameter(s) under MLE inference +#' @param plotprog \code{logical} indicating if progress plots should be made after each inner iteration; +#' the plots show four panels tracking the best valid objective, the ScaledEI or EY surface, +#' the predictive mean and standard deviation of the objective function, over the first two input variables +#' @param ey_tol a scalar proportion indicating how many of the ScaledEIs +#' at \code{ncandf(t)} candidate locations must be non-zero to \dQuote{trust} +#' that metric to guide search, reducing to an EY-based search instead +#' (choosing that proportion to be one forces EY-based search) #' @param verb a non-negative integer indicating the verbosity level; the larger the value the #' more that is printed to the screen #' @param ... additional arguments passed to \code{blackbox} #' #' @details -#' Additional details... +#' To be added. #' #' #' @returns The output is a \code{list} summarizing the progress of the evaluations of the @@ -50,8 +57,8 @@ #' \item{obj }{ vector giving the value of the objective for the input under consideration at each trial } #' \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} #' \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -#' \item{idx_v}{index of the valid solutions} -#' \item{AF_time}{the time required to compute the acquisition function} +#' \item{feasibility}{vector giving the feasibility for the input under consideration at each trial } +#' \item{rho}{vector giving the penalty parameters } #' #' @seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} #' @@ -62,7 +69,11 @@ #' #' @import laGP #' @import tgp -#' @importFrom stats dnorm +#' @importFrom graphics image +#' @importFrom graphics par +#' @importFrom graphics points +#' @importFrom stats dnorm +#' @importFrom stats median #' @importFrom stats optim #' @importFrom stats pnorm #' @importFrom stats sd @@ -77,8 +88,6 @@ #' B = rbind(c(0, 1), c(0, 1)) #' # Objective and constraint function for use #' HSQ = function(x){ -#' B = rbind(c(0, 1), c(0, 1)) -#' x = pmin(pmax(B[,1], x), B[,2]) #' t = function(z) #' return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) #' herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) @@ -157,8 +166,9 @@ optim.EP = function( blackbox, B, equal=FALSE, ethresh=1e-2, Xstart=NULL, start=10, end=100, urate=10, rho=NULL, ncandf=function(k) { k }, - dg_start=c(0.1, 1e-6), dlim=c(1e-4, 10), - plotprog=FALSE, ey_tol=1e-2, AF_tol=sqrt(.Machine$double.eps), + dg_start=c(1e-2*sqrt(nrow(B)), 1e-6), + dlim=c(1e-4, 1)*sqrt(nrow(B)), + plotprog=FALSE, ey_tol=1e-2, verb=2, ...) { ## check start @@ -251,13 +261,11 @@ optim.EP = function( fmean = mean(obj); fsd = sd(obj) # for standard normalization on objective values fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d = dg_start[1], g = dg_start[2], dK = TRUE) df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = verb-1)$d - if(urate > 1){ - deleteGPsep(fgpi) - df[dfdlim[2]] = dlim[2]/10 - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d=df, g=dg_start[2], dK=TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - } + deleteGPsep(fgpi) + df[dfdlim[2]] = dlim[2]/10 + fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d=df, g=dg_start[2], dK=TRUE) + df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d ## initializing constraint surrogates Cgpi = rep(NA, nc) @@ -265,13 +273,11 @@ optim.EP = function( for (j in 1:nc) { Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dg_start[1], g=dg_start[2], dK=TRUE) dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - if(urate > 1){ - deleteGPsep(Cgpi[j]) - dc[j, dc[j,]dlim[2]] = dlim[2]/10 - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dc[j,], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - } + deleteGPsep(Cgpi[j]) + dc[j, dc[j,]dlim[2]] = dlim[2]/10 + Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dc[j,], g=dg_start[2], dK=TRUE) + dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d } ## printing initial design @@ -322,66 +328,35 @@ optim.EP = function( ## calculate composite surrogate, and evaluate SEI and/or EY ncand = ncandf(k) cands = lhs(ncand, Hypercube) # random candidate grid - tic = proc.time()[3] # Start time - if(is.finite(m2) && since > 5 && since %% 2 == 0){ # maximize UEI approach - # by = "ConVar" - # TRlower = pmax(xbest_unit - 0.2, 0) - # TRupper = pmin(xbest_unit + 0.2, 1) - # TRspace = cbind(TRlower, TRupper) - # cands = lhs(ncand, TRspace) - # AF = AF_ConVar(cands, fgpi, Cgpi, equal) - # m = which.max(AF) - # out_AF = optim(par=cands[m, ], fn=AF_ConVar, method="L-BFGS-B", - # lower=TRlower, upper=TRupper, - # control = list(fnscale = -1), # maximization problem - # fgpi=fgpi, Cgpi=Cgpi, equal=equal) - - # by = "LCB" - # AF = AF_LCB(cands, fgpi, fmean, fsd, Cgpi, rho, equal) - # m = which.min(AF) - # out_AF = optim(par=cands[m, ], fn=AF_LCB, method="L-BFGS-B", - # lower=0, upper=1, - # fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - # rho=rho, equal=equal) - - by = "UEI" - AF = AF_UEI(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) + # tic = proc.time()[3] # Start time + AF = AF_ScaledEI(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) + nzsei = sum(AF > sqrt(.Machine$double.eps)) + # Augment the candidate points + if((ey_tol*ncand < nzsei && nzsei <= 0.1*ncand) || (since>1 && since %% 5 == 0)){ + cands = rbind(cands, lhs(10*ncand, Hypercube)) + AF = c(AF, AF_ScaledEI(cands[-(1:ncand),], fgpi, fmean, fsd, Cgpi, epbest, rho, equal)) + nzsei = sum(AF > sqrt(.Machine$double.eps)) + ncand = 11*ncand + } + if(nzsei <= ey_tol*ncand){ # minimize predictive mean approach + by = "EY" + AF = AF_EY(cands, fgpi, fmean, fsd, Cgpi, rho, equal) + m = which.min(AF) + out_AF = optim(par=cands[m, ], fn=AF_EY, method="L-BFGS-B", + lower=0, upper=1, + fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, + rho=rho, equal=equal) + }else{# maximize scaled expected improvement approach + by = "ScaledEI" m = which.max(AF) - out_AF = optim(par=cands[m, ], fn=AF_UEI, method="L-BFGS-B", + out_AF = optim(par=cands[m, ], fn=AF_ScaledEI, method="L-BFGS-B", lower=0, upper=1, fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, control = list(fnscale = -1), # maximization problem epbest=epbest, rho=rho, equal=equal) - }else{ - AF = AF_ScaledEI(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) - nzsei = sum(AF > AF_tol) - # Augment the candidate points - if((ey_tol*ncand < nzsei && nzsei <= 0.1*ncand) || (since>1 && since %% 5 == 0)){ - cands = rbind(cands, lhs(10*ncand, Hypercube)) - AF = c(AF, AF_ScaledEI(cands[-(1:ncand),], fgpi, fmean, fsd, Cgpi, epbest, rho, equal)) - nzsei = sum(AF > AF_tol) - ncand = 11*ncand - } - if(nzsei <= ey_tol*ncand){ # minimize predictive mean approach - by = "EY" - AF = AF_EY(cands, fgpi, fmean, fsd, Cgpi, rho, equal) - m = which.min(AF) - out_AF = optim(par=cands[m, ], fn=AF_EY, method="L-BFGS-B", - lower=0, upper=1, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - rho=rho, equal=equal) - }else{# maximize scaled expected improvement approach - by = "ScaledEI" - m = which.max(AF) - out_AF = optim(par=cands[m, ], fn=AF_ScaledEI, method="L-BFGS-B", - lower=0, upper=1, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - control = list(fnscale = -1), # maximization problem - epbest=epbest, rho=rho, equal=equal) - } } - toc = proc.time()[3] # End time - AF_time = AF_time + toc - tic # AF running time + # toc = proc.time()[3] # End time + # AF_time = AF_time + toc - tic # AF running time ## calculate next point xnext_unit = matrix(out_AF$par, nrow = 1) @@ -439,12 +414,26 @@ optim.EP = function( } ## update GP fits - updateGPsep(fgpi, xnext_unit, (obj[k]-fmean)/fsd, verb = verb-2) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d + updateGPsep(fgpi, xnext_unit, (obj[k]-fmean)/fsd, verb = 0) + df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d + sigma_f = sqrt(predGPsep(fgpi, X_unit, lite=TRUE)$s2) #* fsd + if(median(sigma_f) > 0.1){# rebuild GP as sigma is large at observed points + if(verb > 0) cat(" rebuild objective surrogate as sigma is large at observed points \n") + deleteGPsep(fgpi) + fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d = dlim[1], g = dg_start[2], dK=TRUE) + df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d + } for(j in 1:nc) { - updateGPsep(Cgpi[j], xnext_unit, C_bilog[k,j], verb = verb-2) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d + updateGPsep(Cgpi[j], xnext_unit, C_bilog[k,j], verb = 0) + dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d + sigma_c = sqrt(predGPsep(Cgpi[j], X_unit, lite=TRUE)$s2) + if(median(sigma_c) > 0.1){ # rebuild GP as sigma is large at observed points + if(verb > 0) cat(" rebuild constrained", j, "surrogate as sigma is large at observed points \n") + deleteGPsep(Cgpi[j]) + Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d = dlim[1], g = dg_start[2], dK=TRUE) + dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d + } } ## plot progress @@ -464,7 +453,7 @@ optim.EP = function( image(graphic, xlim=range(X[,1]), ylim=range(X[,2]), main=by) points(X[1:start,1:2], col=feasibility[1:start]+3) points(X[-(1:start),1:2, drop=FALSE], col=feasibility[-(1:start)]+3, pch=19) - points(X[k,,drop=FALSE], col="red", pch=18, cex=1.5) + points(X[k,,drop=FALSE], col="black", pch=18, cex=1.6) ## mean of objective function pred_f = predGPsep(fgpi, cands, lite=TRUE) mu_f = pred_f$mean * fsd + fmean @@ -474,14 +463,14 @@ optim.EP = function( image(graphic, xlim=range(X[,1]), ylim=range(X[,2]), main="mu_f") points(X[1:start,1:2], col=feasibility[1:start]+3) points(X[-(1:start),1:2, drop=FALSE], col=feasibility[-(1:start)]+3, pch=19) - points(X[k,,drop=FALSE], col="red", pch=18, cex=1.5) + points(X[k,,drop=FALSE], col="black", pch=18, cex=1.6) ## standard deviation of objective function graphic = interp.loess(cands_unnormalize[,1], cands_unnormalize[,2], sigma_f, span=span) image(graphic, xlim=range(X[,1]), ylim=range(X[,2]), main="sd_f") points(X[1:start,1:2], col=feasibility[1:start]+3) points(X[-(1:start),1:2, drop=FALSE], col=feasibility[-(1:start)]+3, pch=19) - points(X[k,,drop=FALSE], col="red", pch=18, cex=1.5) + points(X[k,,drop=FALSE], col="black", pch=18, cex=1.6) } } @@ -491,5 +480,5 @@ optim.EP = function( return(list(prog = prog, xbest = xbest, obj = obj, C=C, X = X, - feasibility=feasibility, rho=rho, AF_time=AF_time)) + feasibility=feasibility, rho=rho)) #AF_time = AF_time } diff --git a/R/EPBO4HD.R b/R/EPBO4HD.R deleted file mode 100644 index e51f9b2..0000000 --- a/R/EPBO4HD.R +++ /dev/null @@ -1,601 +0,0 @@ -#' @title Exact Penalty Bayesian Optimization -#' -#' @description Black-box optimization under mixed equality and inequality constraints via an exact penalty. -#' -#' @param blackbox blackbox of an input (\code{x}), facilitating vectorization on a -#' \code{matrix} \code{X} thereof, returning a \code{list} -#' with elements \code{$obj} containing the (scalar) objective value and \code{$c} -#' containing a vector of evaluations of the (multiple) constraint function at \code{x}. -#' @param B 2-column \code{matrix} describing the bounding box. The number of rows -#' of the \code{matrix} determines the input dimension (\code{length(x)} in \code{blackbox(x)}); -#' the first column gives lower bounds and the second gives upper bounds -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' @param ethresh a threshold used for equality constraints to determine validity for -#' progress measures; ignored if there are no equality constraints -#' @param Xstart optional matrix of starting design locations in lieu of, or in addition to, -#' \code{start} random ones; we recommend \code{nrow(Xstart) + start >= 6}; also must -#' have \code{ncol(Xstart) = nrow(B)} -#' @param start positive integer giving the number of random starting locations before -#' sequential design (for optimization) is performed; \code{start >= 6} is -#' recommended unless \code{Xstart} is non-\code{NULL}; in the current version -#' the starting locations come from a space-filling design via \code{\link[tgp]{dopt.gp}} -#' @param end positive integer giving the total number of evaluations/trials in the -#' optimization; must have \code{end > start} -#' @param urate positive integer indicating how many optimization trials should pass before -#' each MLE/MAP update is performed for GP correlation lengthscale -#' parameter(s) -#' @param rho positive vector initial exact penalty parameter in the exact penalty function; -#' the default setting of \code{rho = NULL} causes an automatic starting value to be chosen -#' @param ncandf function taking a single integer indicating the optimization trial number \code{t}, where -#' \code{start < t <= end}, and returning the number of search candidates (e.g., for -#' expected improvement calculations) at round \code{t}; the default setting -#' allows the number of candidates to grow linearly with \code{t} -#' @param dg_start 2-vector giving starting values for the lengthscale and nugget parameters -#' of the GP surrogate model(s) for constraints -#' @param dlim 2-vector giving bounds for the lengthscale parameter(s) under MLE inference -#' @param verb a non-negative integer indicating the verbosity level; the larger the value the -#' more that is printed to the screen -#' @param ... additional arguments passed to \code{blackbox} -#' -#' @details -#' Additional details... -#' -#' -#' @returns The output is a \code{list} summarizing the progress of the evaluations of the -#' blackbox under optimization: -#' \item{prog }{ vector giving the best feasible (\code{g(x) <= 0 && |h(x)| <= ethresh}) value of the objective over the trials } -#' \item{xbest }{ vector giving the recommended solution} -#' \item{obj }{ vector giving the value of the objective for the input under consideration at each trial } -#' \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} -#' \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -#' \item{idx_v}{index of the valid solutions} -#' \item{AF_time}{the time required to compute the acquisition function} -#' -#' @seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @keywords optimize -#' @keywords design -#' -#' @import laGP -#' @import tgp -#' @importFrom stats dnorm -#' @importFrom stats optim -#' @importFrom stats pnorm -#' @importFrom stats sd -#' @importFrom stats runif -#' @importFrom utils tail -#' -#' @export -#' -#' @examples -#' ## Inequality constrained problem -#' # search space -#' B = rbind(c(0, 1), c(0, 1)) -#' # Objective and constraint function for use -#' HSQ = function(x){ -#' B = rbind(c(0, 1), c(0, 1)) -#' x = pmin(pmax(B[,1], x), B[,2]) -#' t = function(z) -#' return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) -#' herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) -#' c1 = 1.5 - x[1] - 2*x[2] - 0.5*sin(2*pi*(x[1]^2 - 2*x[2])) -#' c2 = x[1]^2 + x[2]^2 - 1.5 -#' return(list(obj=herbtooth, c=cbind(c1,c2))) -#' } -#' EPBO = optim.EP(HSQ, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, verb = 0) -#' # progress, best feasible value of the objective over the trials -#' EPBO$prog -#' # the recommended solution -#' EPBO$xbest -#' -#' ## Mixed constrained problem -#' # search space -#' B = rbind(c(0, 1), c(0, 1)) -#' # Objective and constraint function for use -#' goldstein.price = function(X) -#' { -#' if(is.null(nrow(X))) X <- matrix(X, nrow=1) -#' m <- 8.6928 -#' s <- 2.4269 -#' x1 <- 4 * X[,1] - 2 -#' x2 <- 4 * X[,2] - 2 -#' a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2) -#' b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 + 27 * x2^2) -#' f <- log(a * b) -#' f <- (f - m)/s -#' return(f) -#' } -#' -#' toy.c1 <- function(X) { -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' c1 = 3/2 - X[,1] - 2*X[,2] - .5*sin(2*pi*(X[,1]^2 - 2*X[,2])) -#' } -#' -#' branin.c = function(X){ -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' x1 = X[,1] * 15 - 5 -#' x2 = X[,2] * 15 -#' f = (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10 -#' return(-f+25) -#' } -#' -#' parr <- function(X){ -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' x1 <- (2 * X[,1] - 1) -#' x2 <- (2 * X[,2] - 1) -#' g <- (4-2.1*x1^2+1/3*x1^4)*x1^2 + x1*x2 + (-4+4*x2^2)*x2^2 + -#' 3*sin(6*(1-x1)) + 3*sin(6*(1-x2)) -#' return(-g+4) -#' } -#' -#' gsbp.constraints <- function(x){ -#' return(cbind(toy.c1(x), branin.c(x), parr(x))) -#' } -#' -#' # problem definition -#' gsbpprob <- function(X, known.only=FALSE) -#' { -#' if(is.null(nrow(X))) X <- matrix(X, nrow=1) -#' if(known.only) stop("known.only not supported for this example") -#' f = goldstein.price(X) -#' C = gsbp.constraints(X) -#' return(list(obj=f, c=cbind(C[,1], C[,2]/100, C[,3]/10))) -#' } -#' EPBO = optim.EP(gsbpprob, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, -#' ethresh = 1e-2, equal = c(0,1,1), verb = 0) -#' -#' # progress, best feasible value of the objective over the trials -#' EPBO$prog -#' # the recommended solution -#' EPBO$xbest - -optim.EP4HD = function(blackbox, B, - equal=FALSE, ethresh=1e-2, - Xstart=NULL, start=10, end=100, - urate=10, rho=NULL, - dg_start=c(0.1, 1e-6), - dlim=c(1e-4, 10), - ncandf=function(k) { k }, - # ncand=min(5000, max(2000, 200*nrow(B))), - batch_size=1, - trcontrol=list( - length=0.8, length_min=0.5^7, length_max=1.6, - failure_counter=0, failure_tolerance=max(4, nrow(B)), - success_counter=0, success_tolerance=3), - plotprog=FALSE, - opt=TRUE, - ey.tol=1e-2, AF.tol=sqrt(.Machine$double.eps), - verb=1, ...) -{ - ## check start - if(start >= end) stop("must have start < end") - if(start == 0 & is.null(Xstart)) stop("must have start>0 or given Xstart") - - dim = nrow(B) # dimension - # Check trcontrol - if (is.null(trcontrol$length)) trcontrol$length = 0.8 - trcontrol$length0 = trcontrol$length - if (is.null(trcontrol$length_min)) trcontrol$length_min = 0.5^7 - if (is.null(trcontrol$length_max)) trcontrol$length_max = 1.6 - if (is.null(trcontrol$failure_counter)) trcontrol$failure_counter = 0 - if (is.null(trcontrol$failure_tolerance)) trcontrol$failure_tolerance = max(4, dim) - if (is.null(trcontrol$success_counter)) trcontrol$success_counter = 0 - if (is.null(trcontrol$success_tolerance)) trcontrol$success_tolerance = 3 - - # Hypercube [0, 1]^d - Hypercube = matrix(c(rep(0,dim), rep(1,dim)), ncol=2) - - ## get initial design - X = dopt.gp(start, Xcand=lhs(10*start, B))$XX - ## X = lhs(start, B) - X = rbind(Xstart, X) - X_unit = normalize(X, B) - start = nrow(X) - - - ## first run to determine the number of constraints - out = blackbox(X[1,], ...) - nc = length(out$c) - - ## check equal argument - if(length(equal) == 1) { - if(!(equal %in% c(0,1))) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - if(equal) equal = rep(1, nc) - else equal = rep(0, nc) - } else { - if(length(equal) != nc) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - if(! any(equal != 0 | equal != 1)) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - } - equal = as.logical(equal) - - ## allocate progress objects, and initialize - prog = rep(NA, start) # progress - obj = rep(NA, start) # objective values - feasibility = rep(NA, start) # Whether the point is feasible? - C = matrix(NA, nrow=start, ncol=nc) # constraint values - C_bilog = matrix(NA, nrow=start, ncol=nc) # bi-log transformation of constraint values - CV = matrix(NA, nrow=start, ncol=nc) # constraint violations in terms of bi-log transformation - - obj[1] = out$obj; C[1, ] = out$c; C_bilog[1, ] = bilog(C[1, ]) - CV[1, !equal] = pmax(0, C_bilog[1, !equal]) # CV for inequality - CV[1, equal] = abs(C_bilog[1, equal]) # CV for equality - feasibility[1] = all(C[1,!equal] <= 0) && all(abs(C[1,equal]) <= ethresh) - prog[1] = ifelse(feasibility[1], obj[1], Inf) - - ## remainder of starting run - for(k in 2:start) { ## now that problem is vectorized we can probably remove for - out = blackbox(X[k,], ...) - obj[k] = out$obj; C[k,] = out$c; C_bilog[k, ] = bilog(C[k, ]) - CV[k, !equal] = pmax(0, C_bilog[k, !equal]) - CV[k, equal] = abs(C_bilog[k, equal]) - feasibility[k] = all(C[k,!equal] <= 0) && all(abs(C[k,equal]) <= ethresh) - prog[k] = ifelse(feasibility[k] && obj[k] < prog[k-1], obj[k], prog[k-1]) - } - - ## handle initial rho values - if(is.null(rho)){ - if(all(feasibility)){ # - rho = rep(0, nc) - }else { - ECV = colMeans(CV) # averaged CV - rho = mean(abs(obj)) * ECV/sum(ECV^2) - } - if(any(equal)) rho[equal] = pmax(1/ethresh/sum(equal), rho[equal]) - }else{ - if(length(rho) != nc || any(rho < 0)) stop("rho should be a non-negative vector whose length is equal to the number of constraints.") - } - - ## calculate EP for data seen so far - scv = CV%*%rho # weighted sum of the constraint violations - ep = obj + scv # the EP values - epbest = min(ep) # best EP seen so far - m2 = prog[start] # BOFV - since = 0 - ## best solution so far - if(is.finite(m2)){ # if at least one feasible solution was found - xbest = X[which.min(prog),] - }else{ # if none of the solutions are feasible - xbest = X[which.min(scv),] - } - xbest_unit = as.vector(normalize(xbest, B)) - - ## initialize objective surrogate - ab = darg(NULL, X_unit)$ab - fmean = mean(obj); fsd = sd(obj) # for standard normalization on objective values - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d = dg_start[1], g = dg_start[2], dK = TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - if(urate > 1){ - deleteGPsep(fgpi) - df[dfdlim[2]] = dlim[2]/10 - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d=df, g=dg_start[2], dK=TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - } - - ## initializing constraint surrogates - Cgpi = rep(NA, nc) - dc = matrix(NA, nrow=nc, ncol=dim) - for (j in 1:nc) { - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dg_start[1], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - if(urate > 1){ - deleteGPsep(Cgpi[j]) - dc[j, dc[j,]dlim[2]] = dlim[2]/10 - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dc[j,], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - } - } - - ## printing initial design - if(verb > 0) { - cat("The initial design: ") - cat("ab=[", paste(signif(ab,3), collapse=", "), sep="") - cat("]; rho=[", paste(signif(rho,3), collapse=", "), sep="") - cat("]; xbest=[", paste(signif(xbest,3), collapse=" "), sep="") - cat("]; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - - AF_time = 0 # AF running time - - ## iterating over the black box evaluations - for (k in (start+1):end) { - ## rebuild surrogates periodically under new normalized responses - if(k > (start+1) && (k %% urate == 0)) { - ## objective surrogate - deleteGPsep(fgpi) - df[dfdlim[2]] = dlim[2]/10 - fmean = mean(obj); fsd = sd(obj) - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d=df, g=dg_start[2], dK=TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - - ## constraint surrogates - for(j in 1:nc) { - deleteGPsep(Cgpi[j]) - dc[j, dc[j,]dlim[2]] = dlim[2]/10 - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dc[j,], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - } - } - - - ## Update the rho when the smallest EP is not feasible - ck = C[which.min(ep),] - invalid = rep(NA, nc) - invalid[!equal] = (ck[!equal] > 0); invalid[equal] = (abs(ck[equal]) > ethresh) - if(is.finite(m2) && any(invalid) && since > 2){ - rho[invalid] = rho[invalid]*2 - scv = CV%*%rho; ep = obj + scv; epbest = min(ep) - if(verb > 0) cat(" the smallest EP is not feasible, updating rho=(", - paste(signif(rho,3), collapse=", "), ")\n", sep="") - } - - ## random candidate grid - weights = colMeans(rbind(df, dc)) - weights = weights / mean(weights) - weights = weights / prod(weights)^(1/dim) - TRlower = pmax(xbest_unit - weights*trcontrol$length/2, 0) - TRupper = pmin(xbest_unit + weights*trcontrol$length/2, 1) - # TRlower = pmax(xbest_unit - trcontrol$length/2, 0) - # TRupper = pmin(xbest_unit + trcontrol$length/2, 1) - TRspace = cbind(TRlower, TRupper) - ncand = ncandf(k) - if(dim <= 20){ - cands = lhs(ncand, TRspace) - }else{ - pert = lhs(ncand, TRspace) - # create a perturbation mask - prob_perturb = min(20/dim, 1) - mask = matrix(runif(ncand*dim), ncol = dim) <= prob_perturb - ind_mask = which(rowSums(mask) == 0) - mask[ind_mask, sample.int(dim, size = length(ind_mask), replace = TRUE)] = 1 - cands = matrix(rep(xbest_unit, ncand), nrow = ncand, byrow = TRUE) - cands[which(mask==1)] = pert[which(mask==1)] - } - - if(verb == 1) { - cat("k=", k, sep="") - cat(" length=", trcontrol$length, sep="") - } - if(verb == 2){ - cat("k=", k, sep="") - cat(" trcontrol(length=", trcontrol$length, sep="") - cat(", weights=[", paste(signif(weights, 3), collapse=" "), sep="") - cat("], lower=[", paste(signif(unnormalize(TRlower, B),3), collapse=" "), sep="") - cat("], upper=[", paste(signif(unnormalize(TRupper, B),3), collapse=" "), "])\n", sep="") - } - - ## calculate composite surrogate, and evaluate SEI and/or EY - tic = proc.time()[3] # Start time - if(is.finite(m2) && since > 5 && since %% 2 == 0){ # maximize UEI approach - by = "UEI" - AF = AF_UEI(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) - m = which.max(AF) - if(opt && max(AF) > AF.tol){ - out_AF = optim(par=cands[m, ], fn=AF_UEI, method="L-BFGS-B", - lower=TRlower, upper=TRupper, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - control = list(fnscale = -1), # maximization problem - epbest=epbest, rho=rho, equal=equal) - }else{ out_AF = list(par=cands[m, ], value=max(AF)) } - }else{ - AF = AF_ScaledEI(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) - nzsei = sum(AF > AF.tol) - # Augment the candidate points - if((0.01*ncand < nzsei && nzsei <= 0.1*ncand) || (since>1 && since %% 5 == 0)){ - cands = rbind(cands, lhs(10*ncand, TRspace)) - AF = c(AF, AF_ScaledEI(cands[-(1:ncand),], fgpi, fmean, fsd, Cgpi, epbest, rho, equal)) - nzsei = sum(AF > AF.tol) - ncand = 11*ncand - } - if(nzsei <= ey.tol*ncand){ # minimize predictive mean approach - by = "EY" - AF = AF_EY(cands, fgpi, fmean, fsd, Cgpi, rho, equal) - m = which.min(AF) - if(opt){ - out_AF = optim(par=cands[m, ], fn=AF_EY, method="L-BFGS-B", - lower=TRlower, upper=TRupper, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - rho=rho, equal=equal) - }else{ out_AF = list(par=cands[m, ], value=min(AF)) } - }else{ # maximize scaled expected improvement approach - by = "ScaledEI" - m = which.max(AF) - if(opt){ - out_AF = optim(par=cands[m, ], fn=AF_ScaledEI, method="L-BFGS-B", - lower=TRlower, upper=TRupper, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - control = list(fnscale = -1), # maximization problem - epbest=epbest, rho=rho, equal=equal) - }else{ out_AF = list(par=cands[m, ], value=max(AF)) } - # if(verb>0){ - # cat("out_AF(value=", out_AF$value, "par=[", - # paste(signif(out_AF$par,3), collapse=" "), "]\n") - # cat("max_AF(value=", max(AF), "par=[", - # paste(signif(cands[m, ],3), collapse=" "), "]\n") - # } - } - } - toc = proc.time()[3] # End time - AF_time = AF_time + toc - tic # AF running time - - - # AF = AF_TS(cands, fgpi, Cgpi, rho, equal) - # by = "TS" - # m = which.min(AF) - # out = list(par=cands[m,]) - - ## calculate next point - xnext_unit = matrix(out_AF$par, nrow = 1) - X_unit = rbind(X_unit, xnext_unit) - xnext = unnormalize(xnext_unit, B) - X = rbind(X, xnext) - - ## new run - out = blackbox(xnext, ...) - fnext = out$obj; obj = c(obj, fnext); C = rbind(C, out$c) - C_bilog = rbind(C_bilog, bilog(out$c)) - CV = rbind(CV, rep(NA, nc)) - CV[k, !equal] = pmax(0, C_bilog[k, !equal]) # constraint violations for inequality - CV[k, equal] = abs(C_bilog[k, equal]) # constraint violations for equality - - ## check if best valid has changed - feasibility = c(feasibility, all(out$c[!equal] <= 0) && all(abs(out$c[equal]) <= ethresh)) - since = since + 1 - if(feasibility[k] && fnext < prog[k-1]) { - m2 = fnext; since = 0 - } # otherwise m2 unchanged; should be the same as prog[k-1] - prog = c(prog, m2) - - ## rho update - if(all(feasibility)){ # - rho_new = rep(0, nc) - }else { - ECV = colMeans(CV) - rho_new = mean(abs(obj)) * ECV/sum(ECV^2) - } - if(verb > 0 && any(rho_new > rho)){ # printing progress - cat(" updating rho=[", paste(signif(pmax(rho_new, rho),3), collapse=", "), "]\n", sep="") - } - rho = pmax(rho_new, rho) - # if(is.finite(m2) && since > 1 && since %% 5 == 0){ - # rho = rho*2 - # if(verb > 0){cat(" double rho if since%%5==0")} - # } - - ## calculate EP for data seen so far - scv = CV%*%rho; ep = obj + scv; epbest = min(ep) - if(is.finite(m2)){ # best solution so far - xbest = X[which.min(prog),] - }else{ - xbest = X[which.min(scv),] - } - xbest_unit = as.vector(normalize(xbest, B)) - - ## progress meter - if(verb == 1) { - cat(" ", by, "=", out_AF$value, ", ncand=", ncand, sep="") - cat("; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - if(verb == 2) { - cat(" ", by, "=", out_AF$value, ", ncand=", ncand, sep="") - cat("; xnext ([", paste(signif(xnext,3), collapse=" "), - "], feasibility=", feasibility[k], ")\n", sep="") - cat(" xbest=[", paste(signif(xbest,3), collapse=" "), sep="") - cat("]; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - - ## update the control parameter of trust region - if(is.finite(m2)){ - # if at least one feasible solution was found - if(feasibility[k]){ - # The new candidate is feasible - if(is.infinite(prog[k-1]) || fnext < prog[k-1] - 1e-3 * abs(prog[k-1])){ # sufficient decrease condition - trcontrol$success_counter = trcontrol$success_counter + 1 - trcontrol$failure_counter = 0 - }else{ - trcontrol$failure_counter = trcontrol$failure_counter + 1 - trcontrol$success_counter = 0 - } - }else{ - # The new candidate is infeasible - if(is.infinite(prog[k-1]) && which.min(scv) == k){ - trcontrol$success_counter = trcontrol$success_counter + 1 - trcontrol$failure_counter = 0 - }else{ - trcontrol$failure_counter = trcontrol$failure_counter + 1 - trcontrol$success_counter = 0 - } - } - }else{ - # if none of the solutions is feasible - if(which.min(scv) == k){ - trcontrol$success_counter = trcontrol$success_counter + 1 - trcontrol$failure_counter = 0 - }else{ - trcontrol$failure_counter = trcontrol$failure_counter + 1 - trcontrol$success_counter = 0 - } - } - - ## update the trust region - if(trcontrol$success_counter == trcontrol$success_tolerance){ - # Expand trust region - trcontrol$length = min(trcontrol$length * 2, trcontrol$length_max) - trcontrol$success_counter = 0 - }else if(trcontrol$failure_counter == trcontrol$failure_tolerance){ - # Shrink trust region - trcontrol$length = trcontrol$length / 2 - trcontrol$failure_counter = 0 - } - if(trcontrol$length < trcontrol$length_min){ - # Restart when trust region becomes too small - trcontrol$length = trcontrol$length0 - trcontrol$success_counter = 0 - trcontrol$failure_counter = 0 - } - - ## update GP fits - updateGPsep(fgpi, xnext_unit, (obj[k]-fmean)/fsd, verb = 0) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - - for(j in 1:nc) { - updateGPsep(Cgpi[j], xnext_unit, C_bilog[k,j], verb = 0) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = 0)$d - } - - ## plot progress - if(plotprog) { - par(ps=16, mfrow=c(2,2)) - ## progress - if(is.finite(m2)){ - plot(prog, type="l", lwd=1.6, main="progress") - }else{ - plot(prog, type="l", ylim=range(obj), lwd=1.6, main="progress") - } - ## acquisition function - cands_unnormalize = unnormalize(cands, B) - span = ifelse(length(AF) < 30, 0.5, 0.1) - graphic = interp.loess(cands_unnormalize[,1], cands_unnormalize[,2], - as.vector(AF), span=span) - image(graphic, xlim=range(X[,1]), ylim=range(X[,2]), main=by) - points(X[1:start,1:2], col=feasibility[1:start]+3) - points(X[-(1:start),1:2, drop=FALSE], col=feasibility[-(1:start)]+3, pch=19) - points(X[k,,drop=FALSE], col="red", pch=18, cex=1.5) - ## mean of objective function - pred_f = predGPsep(fgpi, cands, lite=TRUE) - mu_f = pred_f$mean * fsd + fmean - sigma_f = sqrt(pred_f$s2) * fsd - graphic = interp.loess(cands_unnormalize[,1], cands_unnormalize[,2], - mu_f, span=span) - image(graphic, xlim=range(X[,1]), ylim=range(X[,2]), main="mu_f") - points(X[1:start,1:2], col=feasibility[1:start]+3) - points(X[-(1:start),1:2, drop=FALSE], col=feasibility[-(1:start)]+3, pch=19) - points(X[k,,drop=FALSE], col="red", pch=18, cex=1.5) - ## standard deviation of objective function - graphic = interp.loess(cands_unnormalize[,1], cands_unnormalize[,2], - sigma_f, span=span) - image(graphic, xlim=range(X[,1]), ylim=range(X[,2]), main="sd_f") - points(X[1:start,1:2], col=feasibility[1:start]+3) - points(X[-(1:start),1:2, drop=FALSE], col=feasibility[-(1:start)]+3, pch=19) - points(X[k,,drop=FALSE], col="red", pch=18, cex=1.5) - } - } - - ## delete GP surrogates - deleteGPsep(fgpi) - for(j in 1:nc) deleteGPsep(Cgpi[j]) - - return(list(prog = prog, xbest = xbest, - obj = obj, C=C, X = X, - feasibility=feasibility, rho=rho, AF_time=AF_time)) -} diff --git a/R/EPBO_Kernel.R b/R/EPBO_Kernel.R deleted file mode 100644 index f07da6c..0000000 --- a/R/EPBO_Kernel.R +++ /dev/null @@ -1,456 +0,0 @@ -#' @title Exact Penalty Bayesian Optimization -#' -#' @description Black-box optimization under mixed equality and inequality constraints via an exact penalty. -#' -#' @param blackbox blackbox of an input (\code{x}), facilitating vectorization on a -#' \code{matrix} \code{X} thereof, returning a \code{list} -#' with elements \code{$obj} containing the (scalar) objective value and \code{$c} -#' containing a vector of evaluations of the (multiple) constraint function at \code{x}. -#' @param B 2-column \code{matrix} describing the bounding box. The number of rows -#' of the \code{matrix} determines the input dimension (\code{length(x)} in \code{blackbox(x)}); -#' the first column gives lower bounds and the second gives upper bounds -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' @param ethresh a threshold used for equality constraints to determine validity for -#' progress measures; ignored if there are no equality constraints -#' @param Xstart optional matrix of starting design locations in lieu of, or in addition to, -#' \code{start} random ones; we recommend \code{nrow(Xstart) + start >= 6}; also must -#' have \code{ncol(Xstart) = nrow(B)} -#' @param start positive integer giving the number of random starting locations before -#' sequential design (for optimization) is performed; \code{start >= 6} is -#' recommended unless \code{Xstart} is non-\code{NULL}; in the current version -#' the starting locations come from a space-filling design via \code{\link[tgp]{dopt.gp}} -#' @param end positive integer giving the total number of evaluations/trials in the -#' optimization; must have \code{end > start} -#' @param urate positive integer indicating how many optimization trials should pass before -#' each MLE/MAP update is performed for GP correlation lengthscale -#' parameter(s) -#' @param rho positive vector initial exact penalty parameter in the exact penalty function; -#' the default setting of \code{rho = NULL} causes an automatic starting value to be chosen -#' @param ncandf function taking a single integer indicating the optimization trial number \code{t}, where -#' \code{start < t <= end}, and returning the number of search candidates (e.g., for -#' expected improvement calculations) at round \code{t}; the default setting -#' allows the number of candidates to grow linearly with \code{t} -#' @param dg_start 2-vector giving starting values for the lengthscale and nugget parameters -#' of the GP surrogate model(s) for constraints -#' @param dlim 2-vector giving bounds for the lengthscale parameter(s) under MLE inference -#' @param verb a non-negative integer indicating the verbosity level; the larger the value the -#' more that is printed to the screen -#' @param ... additional arguments passed to \code{blackbox} -#' -#' @details -#' Additional details... -#' -#' -#' @returns The output is a \code{list} summarizing the progress of the evaluations of the -#' blackbox under optimization: -#' \item{prog }{ vector giving the best feasible (\code{g(x) <= 0 && |h(x)| <= ethresh}) value of the objective over the trials } -#' \item{xbest }{ vector giving the recommended solution} -#' \item{obj }{ vector giving the value of the objective for the input under consideration at each trial } -#' \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} -#' \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -#' \item{idx_v}{index of the valid solutions} -#' \item{AF_time}{the time required to compute the acquisition function} -#' -#' @seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @keywords optimize -#' @keywords design -#' -#' @import DiceKriging -#' @import tgp -#' @importFrom stats dnorm -#' @importFrom stats optim -#' @importFrom stats pnorm -#' @importFrom stats sd -#' @importFrom stats runif -#' @importFrom utils tail -#' -#' @export -#' -#' @examples -#' ## Inequality constrained problem -#' # search space -#' B = rbind(c(0, 1), c(0, 1)) -#' # Objective and constraint function for use -#' HSQ = function(x){ -#' B = rbind(c(0, 1), c(0, 1)) -#' x = pmin(pmax(B[,1], x), B[,2]) -#' t = function(z) -#' return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) -#' herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) -#' c1 = 1.5 - x[1] - 2*x[2] - 0.5*sin(2*pi*(x[1]^2 - 2*x[2])) -#' c2 = x[1]^2 + x[2]^2 - 1.5 -#' return(list(obj=herbtooth, c=cbind(c1,c2))) -#' } -#' EPBO = optim.EP(HSQ, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, verb = 0) -#' # progress, best feasible value of the objective over the trials -#' EPBO$prog -#' # the recommended solution -#' EPBO$xbest -#' -#' ## Mixed constrained problem -#' # search space -#' B = rbind(c(0, 1), c(0, 1)) -#' # Objective and constraint function for use -#' goldstein.price = function(X) -#' { -#' if(is.null(nrow(X))) X <- matrix(X, nrow=1) -#' m <- 8.6928 -#' s <- 2.4269 -#' x1 <- 4 * X[,1] - 2 -#' x2 <- 4 * X[,2] - 2 -#' a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2) -#' b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 + 27 * x2^2) -#' f <- log(a * b) -#' f <- (f - m)/s -#' return(f) -#' } -#' -#' toy.c1 <- function(X) { -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' c1 = 3/2 - X[,1] - 2*X[,2] - .5*sin(2*pi*(X[,1]^2 - 2*X[,2])) -#' } -#' -#' branin.c = function(X){ -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' x1 = X[,1] * 15 - 5 -#' x2 = X[,2] * 15 -#' f = (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10 -#' return(-f+25) -#' } -#' -#' parr <- function(X){ -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' x1 <- (2 * X[,1] - 1) -#' x2 <- (2 * X[,2] - 1) -#' g <- (4-2.1*x1^2+1/3*x1^4)*x1^2 + x1*x2 + (-4+4*x2^2)*x2^2 + -#' 3*sin(6*(1-x1)) + 3*sin(6*(1-x2)) -#' return(-g+4) -#' } -#' -#' gsbp.constraints <- function(x){ -#' return(cbind(toy.c1(x), branin.c(x), parr(x))) -#' } -#' -#' # problem definition -#' gsbpprob <- function(X, known.only=FALSE) -#' { -#' if(is.null(nrow(X))) X <- matrix(X, nrow=1) -#' if(known.only) stop("known.only not supported for this example") -#' f = goldstein.price(X) -#' C = gsbp.constraints(X) -#' return(list(obj=f, c=cbind(C[,1], C[,2]/100, C[,3]/10))) -#' } -#' EPBO = optim.EP(gsbpprob, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, -#' ethresh = 1e-2, equal = c(0,1,1), verb = 0) -#' -#' # progress, best feasible value of the objective over the trials -#' EPBO$prog -#' # the recommended solution -#' EPBO$xbest - -optim.EP.kernel = function( - blackbox, B, equal=FALSE, ethresh=1e-2, Xstart=NULL, start=10, end=100, - urate=10, rho=NULL, ncandf=function(k) { k }, - kmcontrol=list(formula=~1, #coef.trend=0, - covtype="matern5_2", nugget=1e-6, - parinit=rep(0.1, nrow(B)), - lower = rep(1e-4, nrow(B)), upper = rep(10, nrow(B))), - verb=2, ...) -{ - ## check start - if(start >= end) stop("must have start < end") - if(start == 0 & is.null(Xstart)) stop("must have start>0 or given Xstart") - - dim = nrow(B) # dimension - - # Check kmcontrol - if (is.null(kmcontrol$formula)) kmcontrol$formula = ~1 - # if (is.null(kmcontrol$coef.trend)) kmcontrol$coef.trend = 0 - if (is.null(kmcontrol$covtype)) kmcontrol$covtype = "matern5_2" - if (is.null(kmcontrol$nugget)) kmcontrol$nugget = 1e-6 - if (is.null(kmcontrol$parinit)) kmcontrol$parinit = rep(0.1, dim) - if (is.null(kmcontrol$lower)) kmcontrol$lower = rep(1e-4, dim) - if (is.null(kmcontrol$upper)) kmcontrol$upper = rep(10, dim) - - - # Hypercube [0, 1]^d - Hypercube = matrix(c(rep(0,dim), rep(1,dim)), ncol=2) - - ## get initial design - X = dopt.gp(start, Xcand=lhs(10*start, B))$XX - ## X = lhs(start, B) - X = rbind(Xstart, X) - X_unit = normalize(X, B) - start = nrow(X) - - - ## first run to determine the number of constraints - out = blackbox(X[1,], ...) - nc = length(out$c) - - ## check equal argument - if(length(equal) == 1) { - if(!(equal %in% c(0,1))) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - if(equal) equal = rep(1, nc) - else equal = rep(0, nc) - } else { - if(length(equal) != nc) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - if(! any(equal != 0 | equal != 1)) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - } - equal = as.logical(equal) - - ## allocate progress objects, and initialize - prog = rep(NA, start) # progress - obj = rep(NA, start) # objective values - feasibility = rep(NA, start) # Whether the point is feasible? - C = matrix(NA, nrow=start, ncol=nc) # constraint values - C_bilog = matrix(NA, nrow=start, ncol=nc) # bi-log transformation of constraint values - CV = matrix(NA, nrow=start, ncol=nc) # constraint violations in terms of bi-log transformation - - obj[1] = out$obj; C[1, ] = out$c; C_bilog[1, ] = bilog(C[1, ]) - CV[1, !equal] = pmax(0, C_bilog[1, !equal]) # CV for inequality - CV[1, equal] = abs(C_bilog[1, equal]) # CV for equality - feasibility[1] = all(C[1,!equal] <= 0) && all(abs(C[1,equal]) <= ethresh) - prog[1] = ifelse(feasibility[1], obj[1], Inf) - - ## remainder of starting run - for(k in 2:start) { ## now that problem is vectorized we can probably remove for - out = blackbox(X[k,], ...) - obj[k] = out$obj; C[k,] = out$c; C_bilog[k, ] = bilog(C[k, ]) - CV[k, !equal] = pmax(0, C_bilog[k, !equal]) - CV[k, equal] = abs(C_bilog[k, equal]) - feasibility[k] = all(C[k,!equal] <= 0) && all(abs(C[k,equal]) <= ethresh) - prog[k] = ifelse(feasibility[k] && obj[k] < prog[k-1], obj[k], prog[k-1]) - } - - ## handle initial rho values - if(is.null(rho)){ - if(all(feasibility)){ # - rho = rep(0, nc) - }else { - ECV = colMeans(CV) # averaged CV - rho = mean(abs(obj)) * ECV/sum(ECV^2) - } - if(any(equal)) rho[equal] = pmax(1/ethresh/sum(equal), rho[equal]) - }else{ - if(length(rho) != nc || any(rho < 0)) stop("rho should be a non-negative vector whose length is equal to the number of constraints.") - } - - ## calculate EP for data seen so far - scv = CV%*%rho # weighted sum of the constraint violations - ep = obj + scv # the EP values - epbest = min(ep) # best EP seen so far - m2 = prog[start] # BOFV - since = 0 - ## best solution so far - if(is.finite(m2)){ # if at least one feasible solution was found - xbest = X[which.min(prog),] - }else{ # if none of the solutions are feasible - xbest = X[which.min(scv),] - } - - ## initialize objective surrogate - fmean = mean(obj); fsd = sd(obj) # for standard normalization on objective values - fgpi = km(formula = kmcontrol$formula, design = X_unit, response = (obj-fmean)/fsd, - covtype = kmcontrol$covtype, nugget = kmcontrol$nugget, - # coef.trend = kmcontrol$coef.trend, - parinit = kmcontrol$parinit, lower = kmcontrol$lower, upper = kmcontrol$upper, - control=list(trace=0)) - df = fgpi@covariance@range.val - if(urate > 1){ - fgpi = km(formula = kmcontrol$formula, design = X_unit, response = (obj-fmean)/fsd, - covtype = kmcontrol$covtype, nugget = kmcontrol$nugget, - # coef.trend = kmcontrol$coef.trend, - parinit = df, lower = kmcontrol$lower, upper = kmcontrol$upper, - control=list(trace=0)) - df = fgpi@covariance@range.val - } - - ## initializing constraint surrogates - Cgpi = vector("list", nc) - dc = matrix(NA, nrow=nc, ncol=dim) - for (j in 1:nc) { - Cgpi[[j]] = km(formula = kmcontrol$formula, design = X_unit, response = C_bilog[,j], - covtype = kmcontrol$covtype, nugget = kmcontrol$nugget, - # coef.trend = kmcontrol$coef.trend, - parinit = kmcontrol$parinit, lower = kmcontrol$lower, upper = kmcontrol$upper, - control=list(trace=0)) - dc[j,] = Cgpi[[j]]@covariance@range.val - if(urate > 1){ - Cgpi[[j]] = km(formula = kmcontrol$formula, design = X_unit, response = C_bilog[,j], - covtype = kmcontrol$covtype, nugget = kmcontrol$nugget, - # coef.trend = kmcontrol$coef.trend, - parinit = dc[j,], lower = kmcontrol$lower, upper = kmcontrol$upper, - control=list(trace=0)) - dc[j,] = Cgpi[[j]]@covariance@range.val - } - } - - ## printing initial design - if(verb > 0) { - cat("The initial design: ") - cat("rho=[", paste(signif(rho,3), collapse=", "), sep="") - cat("]; xbest=[", paste(signif(xbest,3), collapse=" "), sep="") - cat("]; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - - AF_time = 0 # AF running time - - ## iterating over the black box evaluations - for (k in (start+1):end) { - ## rebuild surrogates periodically under new normalized responses - if(k > (start+1) && (k %% urate == 0)) { - ## objective surrogate - fmean = mean(obj); fsd = sd(obj) - fgpi = km(formula = kmcontrol$formula, design = X_unit, response = (obj-fmean)/fsd, - covtype = kmcontrol$covtype, nugget = kmcontrol$nugget, - # coef.trend = kmcontrol$coef.trend, - parinit = df, lower = kmcontrol$lower, upper = kmcontrol$upper, - control=list(trace=0)) - df = fgpi@covariance@range.val - - ## constraint surrogates - for(j in 1:nc) { - Cgpi[[j]] = km(formula = kmcontrol$formula, design = X_unit, response = C_bilog[,j], - covtype = kmcontrol$covtype, nugget = kmcontrol$nugget, - # coef.trend = kmcontrol$coef.trend, - parinit = dc[j,], lower = kmcontrol$lower, upper = kmcontrol$upper, - control=list(trace=0)) - dc[j,] = Cgpi[[j]]@covariance@range.val - } - } - - - ## Update the rho when the smallest EP is not feasible - ck = C[which.min(ep),] - invalid = rep(NA, nc) - invalid[!equal] = (ck[!equal] > 0); invalid[equal] = (abs(ck[equal]) > ethresh) - if(is.finite(m2) && any(invalid) && since > 2){ - rho[invalid] = rho[invalid]*2 - scv = CV%*%rho; ep = obj + scv; epbest = min(ep) - if(verb > 0) cat(" the smallest EP is not feasible, updating rho=(", - paste(signif(rho,3), collapse=", "), ")\n", sep="") - } - - ## calculate composite surrogate, and evaluate SEI and/or EY - ncand = ncandf(k) - cands = lhs(ncand, Hypercube) # random candidate grid - tic = proc.time()[3] # Start time - if(is.finite(m2) && since > 5 && since %% 2 == 0){ # maximize UEI approach - by = "UEI" - AF = AF_UEI_Kernel(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) - m = which.max(AF) - out_AF = optim(par=cands[m, ], fn=AF_UEI_Kernel, method="L-BFGS-B", - lower=0, upper=1, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - control = list(fnscale = -1), # maximization problem - epbest=epbest, rho=rho, equal=equal) - }else{ - AF = AF_ScaledEI_Kernel(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) - nzsei = sum(AF > sqrt(.Machine$double.eps)) - # Augment the candidate points - if((0.01*ncand < nzsei && nzsei <= 0.1*ncand) || (since>1 && since %% 5 == 0)){ - cands = rbind(cands, lhs(10*ncand, Hypercube)) - AF = c(AF, AF_ScaledEI_Kernel(cands[-(1:ncand),], fgpi, fmean, fsd, Cgpi, epbest, rho, equal)) - nzsei = sum(AF > sqrt(.Machine$double.eps)) - ncand = 11*ncand - } - if(nzsei <= 0.01*ncand){ # minimize predictive mean approach - by = "ey" - AF = AF_EY_Kernel(cands, fgpi, fmean, fsd, Cgpi, rho, equal) - m = which.min(AF) - out_AF = optim(par=cands[m, ], fn=AF_EY_Kernel, method="L-BFGS-B", - lower=0, upper=1, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - rho=rho, equal=equal) - }else{# maximize scaled expected improvement approach - by = "sei" - m = which.max(AF) - out_AF = optim(par=cands[m, ], fn=AF_ScaledEI_Kernel, method="L-BFGS-B", - lower=0, upper=1, - fgpi=fgpi, Cgpi=Cgpi, fmean=fmean, fsd=fsd, - control = list(fnscale = -1), # maximization problem - epbest=epbest, rho=rho, equal=equal) - } - } - toc = proc.time()[3] # End time - AF_time = AF_time + toc - tic # AF running time - - ## calculate next point - xnext_unit = matrix(out_AF$par, nrow = 1) - X_unit = rbind(X_unit, xnext_unit) - xnext = unnormalize(xnext_unit, B) - X = rbind(X, xnext) - - ## new run - out = blackbox(xnext, ...) - fnext = out$obj; obj = c(obj, fnext); C = rbind(C, out$c) - C_bilog = rbind(C_bilog, bilog(out$c)) - CV = rbind(CV, rep(NA, nc)) - CV[k, !equal] = pmax(0, C_bilog[k, !equal]) # constraint violations for inequality - CV[k, equal] = abs(C_bilog[k, equal]) # constraint violations for equality - - ## check if best valid has changed - feasibility = c(feasibility, all(out$c[!equal] <= 0) && all(abs(out$c[equal]) <= ethresh)) - since = since + 1 - if(feasibility[k] && fnext < prog[k-1]) { - m2 = fnext; since = 0 - } # otherwise m2 unchanged; should be the same as prog[k-1] - prog = c(prog, m2) - - ## rho update - if(all(feasibility)){ # - rho_new = rep(0, nc) - }else { - ECV = colMeans(CV) - rho_new = mean(abs(obj)) * ECV/sum(ECV^2) - } - if(verb > 0 && any(rho_new > rho)){ # printing progress - cat(" updating rho=(", paste(signif(pmax(rho_new, rho),3), collapse=", "), ")\n", sep="") - } - rho = pmax(rho_new, rho) - - ## calculate EP for data seen so far - scv = CV%*%rho; ep = obj + scv; epbest = min(ep) - if(is.finite(m2)){ # best solution so far - xbest = X[which.min(prog),] - }else{ - xbest = X[which.min(scv),] - } - - ## progress meter - if(verb > 0) { - cat("k=", k, " ", sep="") - cat(by, "=", out_AF$value, sep="") - cat("; xnext ([", paste(signif(xnext,3), collapse=" "), - "], feasibility=", feasibility[k], ")\n", sep="") - cat(" xbest=[", paste(signif(xbest,3), collapse=" "), sep="") - cat("]; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - - ## update GP fits - fgpi = update(object = fgpi, newX = xnext_unit, newy = (obj[k]-fmean)/fsd, - cov.reestim = TRUE, #trend.reestim = FALSE, - newX.alreadyExist = FALSE) - df = fgpi@covariance@range.val - - for(j in 1:nc) { - Cgpi[[j]] = update(object = Cgpi[[j]], newX = xnext_unit, newy = C_bilog[k,j], - cov.reestim = TRUE, #trend.reestim = FALSE, - newX.alreadyExist = FALSE) - dc[j,] = Cgpi[[j]]@covariance@range.val - } - } - - return(list(prog = prog, xbest = xbest, - obj = obj, C=C, X = X, fgpi=fgpi, Cgpi=Cgpi, - feasibility=feasibility, rho=rho, AF_time=AF_time)) -} diff --git a/R/EPBO_MC.R b/R/EPBO_MC.R deleted file mode 100644 index 5830dfb..0000000 --- a/R/EPBO_MC.R +++ /dev/null @@ -1,430 +0,0 @@ -#' @title Exact Penalty Bayesian Optimization -#' -#' @description Black-box optimization under mixed equality and inequality constraints via an exact penalty. -#' -#' @param blackbox blackbox of an input (\code{x}), facilitating vectorization on a -#' \code{matrix} \code{X} thereof, returning a \code{list} -#' with elements \code{$obj} containing the (scalar) objective value and \code{$c} -#' containing a vector of evaluations of the (multiple) constraint function at \code{x}. -#' @param B 2-column \code{matrix} describing the bounding box. The number of rows -#' of the \code{matrix} determines the input dimension (\code{length(x)} in \code{blackbox(x)}); -#' the first column gives lower bounds and the second gives upper bounds -#' @param equal an optional vector containing zeros and ones, whose length equals the number of -#' constraints, specifying which should be treated as equality constraints (\code{1}) and -#' which as inequality (\code{0}) -#' @param ethresh a threshold used for equality constraints to determine validity for -#' progress measures; ignored if there are no equality constraints -#' @param Xstart optional matrix of starting design locations in lieu of, or in addition to, -#' \code{start} random ones; we recommend \code{nrow(Xstart) + start >= 6}; also must -#' have \code{ncol(Xstart) = nrow(B)} -#' @param start positive integer giving the number of random starting locations before -#' sequential design (for optimization) is performed; \code{start >= 6} is -#' recommended unless \code{Xstart} is non-\code{NULL}; in the current version -#' the starting locations come from a space-filling design via \code{\link[tgp]{dopt.gp}} -#' @param end positive integer giving the total number of evaluations/trials in the -#' optimization; must have \code{end > start} -#' @param urate positive integer indicating how many optimization trials should pass before -#' each MLE/MAP update is performed for GP correlation lengthscale -#' parameter(s) -#' @param rho positive vector initial exact penalty parameter in the exact penalty function; -#' the default setting of \code{rho = NULL} causes an automatic starting value to be chosen -#' @param ncandf function taking a single integer indicating the optimization trial number \code{t}, where -#' \code{start < t <= end}, and returning the number of search candidates (e.g., for -#' expected improvement calculations) at round \code{t}; the default setting -#' allows the number of candidates to grow linearly with \code{t} -#' @param dg_start 2-vector giving starting values for the lengthscale and nugget parameters -#' of the GP surrogate model(s) for constraints -#' @param dlim 2-vector giving bounds for the lengthscale parameter(s) under MLE inference -#' @param N positive scalar integer indicating the number of Monte Carlo samples to be used -#' for calculating ScaledEI -#' @param verb a non-negative integer indicating the verbosity level; the larger the value the -#' more that is printed to the screen -#' @param ... additional arguments passed to \code{blackbox} -#' -#' @details -#' Additional details... -#' -#' -#' @returns The output is a \code{list} summarizing the progress of the evaluations of the -#' blackbox under optimization: -#' \item{prog }{ vector giving the best feasible (\code{g(x) <= 0 && |h(x)| <= ethresh}) value of the objective over the trials } -#' \item{xbest }{ vector giving the recommended solution} -#' \item{obj }{ vector giving the value of the objective for the input under consideration at each trial } -#' \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} -#' \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -#' \item{idx_v}{index of the valid solutions} -#' \item{AF_time}{the time required to compute the acquisition function} -#' -#' @seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @keywords optimize -#' @keywords design -#' -#' @import laGP -#' @import tgp -#' @import matrixStats -#' @importFrom stats dnorm -#' @importFrom stats optim -#' @importFrom stats pnorm -#' @importFrom stats sd -#' @importFrom stats runif -#' @importFrom utils tail -#' -#' @export -#' -#' @examples -#' ## Inequality constrained problem -#' # search space -#' B = rbind(c(0, 1), c(0, 1)) -#' # Objective and constraint function for use -#' HSQ = function(x){ -#' B = rbind(c(0, 1), c(0, 1)) -#' x = pmin(pmax(B[,1], x), B[,2]) -#' t = function(z) -#' return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) -#' herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) -#' c1 = 1.5 - x[1] - 2*x[2] - 0.5*sin(2*pi*(x[1]^2 - 2*x[2])) -#' c2 = x[1]^2 + x[2]^2 - 1.5 -#' return(list(obj=herbtooth, c=cbind(c1,c2))) -#' } -#' EPBO = optim.EP(HSQ, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, verb = 0) -#' # progress, best feasible value of the objective over the trials -#' EPBO$prog -#' # the recommended solution -#' EPBO$xbest -#' -#' ## Mixed constrained problem -#' # search space -#' B = rbind(c(0, 1), c(0, 1)) -#' # Objective and constraint function for use -#' goldstein.price = function(X) -#' { -#' if(is.null(nrow(X))) X <- matrix(X, nrow=1) -#' m <- 8.6928 -#' s <- 2.4269 -#' x1 <- 4 * X[,1] - 2 -#' x2 <- 4 * X[,2] - 2 -#' a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2) -#' b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 + 27 * x2^2) -#' f <- log(a * b) -#' f <- (f - m)/s -#' return(f) -#' } -#' -#' toy.c1 <- function(X) { -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' c1 = 3/2 - X[,1] - 2*X[,2] - .5*sin(2*pi*(X[,1]^2 - 2*X[,2])) -#' } -#' -#' branin.c = function(X){ -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' x1 = X[,1] * 15 - 5 -#' x2 = X[,2] * 15 -#' f = (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10 -#' return(-f+25) -#' } -#' -#' parr <- function(X){ -#' if (is.null(dim(X))) X <- matrix(X, nrow=1) -#' x1 <- (2 * X[,1] - 1) -#' x2 <- (2 * X[,2] - 1) -#' g <- (4-2.1*x1^2+1/3*x1^4)*x1^2 + x1*x2 + (-4+4*x2^2)*x2^2 + -#' 3*sin(6*(1-x1)) + 3*sin(6*(1-x2)) -#' return(-g+4) -#' } -#' -#' gsbp.constraints <- function(x){ -#' return(cbind(toy.c1(x), branin.c(x), parr(x))) -#' } -#' -#' # problem definition -#' gsbpprob <- function(X, known.only=FALSE) -#' { -#' if(is.null(nrow(X))) X <- matrix(X, nrow=1) -#' if(known.only) stop("known.only not supported for this example") -#' f = goldstein.price(X) -#' C = gsbp.constraints(X) -#' return(list(obj=f, c=cbind(C[,1], C[,2]/100, C[,3]/10))) -#' } -#' EPBO = optim.EP(gsbpprob, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, -#' ethresh = 1e-2, equal = c(0,1,1), verb = 0) -#' -#' # progress, best feasible value of the objective over the trials -#' EPBO$prog -#' # the recommended solution -#' EPBO$xbest - -optim.EP.MC = function(blackbox, B, - equal=FALSE, ethresh=1e-2, - Xstart=NULL, start=10, end=100, - urate=10, rho=NULL, - ncandf=function(k) { k }, - dg_start=c(0.1, 1e-6), - dlim=c(1e-4, 10), - N = 1000, - plotprog=FALSE, - verb=2, ...) -{ - ## check start - if(start >= end) stop("must have start < end") - if(start == 0 & is.null(Xstart)) stop("must have start>0 or given Xstart") - - dim = nrow(B) # dimension - # Hypercube [0, 1]^d - Hypercube = matrix(c(rep(0,dim), rep(1,dim)), ncol=2) - - ## get initial design - X = dopt.gp(start, Xcand=lhs(10*start, B))$XX - ## X = lhs(start, B) - X = rbind(Xstart, X) - X_unit = normalize(X, B) - start = nrow(X) - - - ## first run to determine the number of constraints - out = blackbox(X[1,], ...) - nc = length(out$c) - - ## check equal argument - if(length(equal) == 1) { - if(!(equal %in% c(0,1))) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - if(equal) equal = rep(1, nc) - else equal = rep(0, nc) - } else { - if(length(equal) != nc) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - if(! any(equal != 0 | equal != 1)) - stop("equal should be a logical scalar or vector of lenghth(blackbox(x)$c)") - } - equal = as.logical(equal) - - ## allocate progress objects, and initialize - prog = rep(NA, start) # progress - obj = rep(NA, start) # objective values - feasibility = rep(NA, start) # Whether the point is feasible? - C = matrix(NA, nrow=start, ncol=nc) # constraint values - C_bilog = matrix(NA, nrow=start, ncol=nc) # bi-log transformation of constraint values - CV = matrix(NA, nrow=start, ncol=nc) # constraint violations in terms of bi-log transformation - - obj[1] = out$obj; C[1, ] = out$c; C_bilog[1, ] = bilog(C[1, ]) - CV[1, !equal] = pmax(0, C_bilog[1, !equal]) # CV for inequality - CV[1, equal] = abs(C_bilog[1, equal]) # CV for equality - feasibility[1] = all(C[1,!equal] <= 0) && all(abs(C[1,equal]) <= ethresh) - prog[1] = ifelse(feasibility[1], obj[1], Inf) - - ## remainder of starting run - for(k in 2:start) { ## now that problem is vectorized we can probably remove for - out = blackbox(X[k,], ...) - obj[k] = out$obj; C[k,] = out$c; C_bilog[k, ] = bilog(C[k, ]) - CV[k, !equal] = pmax(0, C_bilog[k, !equal]) - CV[k, equal] = abs(C_bilog[k, equal]) - feasibility[k] = all(C[k,!equal] <= 0) && all(abs(C[k,equal]) <= ethresh) - prog[k] = ifelse(feasibility[k] && obj[k] < prog[k-1], obj[k], prog[k-1]) - } - - ## handle initial rho values - if(is.null(rho)){ - if(all(feasibility)){ # - rho = rep(0, nc) - }else { - ECV = colMeans(CV) # averaged CV - rho = mean(abs(obj)) * ECV/sum(ECV^2) - } - if(any(equal)) rho[equal] = pmax(1/ethresh/sum(equal), rho[equal]) - }else{ - if(length(rho) != nc || any(rho < 0)) stop("rho should be a non-negative vector whose length is equal to the number of constraints.") - } - - ## calculate EP for data seen so far - scv = CV%*%rho # weighted sum of the constraint violations - ep = obj + scv # the EP values - epbest = min(ep) # best EP seen so far - m2 = prog[start] # BOFV - since = 0 - ## best solution so far - if(is.finite(m2)){ # if at least one feasible solution was found - xbest = X[which.min(prog),] - }else{ # if none of the solutions are feasible - xbest = X[which.min(scv),] - } - - ## initialize objective surrogate - ab = darg(NULL, X_unit)$ab - fmean = mean(obj); fsd = sd(obj) # for standard normalization on objective values - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d = dg_start[1], g = dg_start[2], dK = TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb = verb-1)$d - if(urate > 1){ - deleteGPsep(fgpi) - df[dfdlim[2]] = dlim[2]/10 - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d=df, g=dg_start[2], dK=TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - } - - ## initializing constraint surrogates - Cgpi = rep(NA, nc) - dc = matrix(NA, nrow=nc, ncol=dim) - for (j in 1:nc) { - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dg_start[1], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - if(urate > 1){ - deleteGPsep(Cgpi[j]) - dc[j, dc[j,]dlim[2]] = dlim[2]/10 - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dc[j,], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - } - } - - ## printing initial design - if(verb > 0) { - cat("The initial design: ") - cat("ab=[", paste(signif(ab,3), collapse=", "), sep="") - cat("]; rho=[", paste(signif(rho,3), collapse=", "), sep="") - cat("]; xbest=[", paste(signif(xbest,3), collapse=" "), sep="") - cat("]; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - - AF_time = 0 # AF running time - - ## iterating over the black box evaluations - for (k in (start+1):end) { - ## rebuild surrogates periodically under new normalized responses - if(k > (start+1) && (k %% urate == 0)) { - ## objective surrogate - deleteGPsep(fgpi) - df[dfdlim[2]] = dlim[2]/10 - fmean = mean(obj); fsd = sd(obj) - fgpi = newGPsep(X_unit, (obj-fmean)/fsd, d=df, g=dg_start[2], dK=TRUE) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - - ## constraint surrogates - for(j in 1:nc) { - deleteGPsep(Cgpi[j]) - dc[j, dc[j,]dlim[2]] = dlim[2]/10 - Cgpi[j] = newGPsep(X_unit, C_bilog[,j], d=dc[j,], g=dg_start[2], dK=TRUE) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - } - } - - - ## Update the rho when the smallest EP is not feasible - ck = C[which.min(ep),] - invalid = rep(NA, nc) - invalid[!equal] = (ck[!equal] > 0); invalid[equal] = (abs(ck[equal]) > ethresh) - if(is.finite(m2) && any(invalid) && since > 2){ - rho[invalid] = rho[invalid]*2 - scv = CV%*%rho; ep = obj + scv; epbest = min(ep) - if(verb > 0) cat(" the smallest EP is not feasible, updating rho=(", - paste(signif(rho,3), collapse=", "), ")\n", sep="") - } - - ## calculate composite surrogate, and evaluate SEI and/or EY - ncand = ncandf(k) - cands = lhs(ncand, Hypercube) # random candidate grid - tic = proc.time()[3] # Start time - if(is.finite(m2) && since > 5 && since %% 2 == 0){ # maximize UEI approach - by = "UEI" - AF = AF_UEI_MC(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N) - m = which.max(AF) - out_AF = list(par=cands[m, ], value=AF[m]) - }else{ - AF = AF_ScaledEI_MC(cands, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N) - nzsei = sum(AF > sqrt(.Machine$double.eps)) - # Augment the candidate points - if((0.01*ncand < nzsei && nzsei <= 0.1*ncand) || (since>1 && since %% 5 == 0)){ - cands = rbind(cands, lhs(10*ncand, Hypercube)) - AF = c(AF, AF_ScaledEI_MC( - cands[-(1:ncand),], fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N)) - nzsei = sum(AF > sqrt(.Machine$double.eps)) - ncand = 11*ncand - } - if(nzsei <= 0.01*ncand){ # minimize predictive mean approach - by = "ey" - AF = AF_EY(cands, fgpi, fmean, fsd, Cgpi, rho, equal) - m = which.min(AF) - out_AF = list(par=cands[m, ], value=AF[m]) - }else{# maximize scaled expected improvement approach - by = "sei" - m = which.max(AF) - out_AF = list(par=cands[m, ], value=AF[m]) - } - } - toc = proc.time()[3] # End time - AF_time = AF_time + toc - tic # AF running time - - ## calculate next point - xnext_unit = matrix(out_AF$par, nrow = 1) - X_unit = rbind(X_unit, xnext_unit) - xnext = unnormalize(xnext_unit, B) - X = rbind(X, xnext) - - ## new run - out = blackbox(xnext, ...) - fnext = out$obj; obj = c(obj, fnext); C = rbind(C, out$c) - C_bilog = rbind(C_bilog, bilog(out$c)) - CV = rbind(CV, rep(NA, nc)) - CV[k, !equal] = pmax(0, C_bilog[k, !equal]) # constraint violations for inequality - CV[k, equal] = abs(C_bilog[k, equal]) # constraint violations for equality - - ## check if best valid has changed - feasibility = c(feasibility, all(out$c[!equal] <= 0) && all(abs(out$c[equal]) <= ethresh)) - since = since + 1 - if(feasibility[k] && fnext < prog[k-1]) { - m2 = fnext; since = 0 - } # otherwise m2 unchanged; should be the same as prog[k-1] - prog = c(prog, m2) - - ## rho update - if(all(feasibility)){ # - rho_new = rep(0, nc) - }else { - ECV = colMeans(CV) - rho_new = mean(abs(obj)) * ECV/sum(ECV^2) - } - if(verb > 0 && any(rho_new > rho)){ # printing progress - cat(" updating rho=(", paste(signif(pmax(rho_new, rho),3), collapse=", "), ")\n", sep="") - } - rho = pmax(rho_new, rho) - - ## calculate EP for data seen so far - scv = CV%*%rho; ep = obj + scv; epbest = min(ep) - if(is.finite(m2)){ # best solution so far - xbest = X[which.min(prog),] - }else{ - xbest = X[which.min(scv),] - } - - ## progress meter - if(verb > 0) { - cat("k=", k, " ", sep="") - cat(by, "=", out_AF$value, sep="") - cat("; xnext ([", paste(signif(xnext,3), collapse=" "), - "], feasibility=", feasibility[k], ")\n", sep="") - cat(" xbest=[", paste(signif(xbest,3), collapse=" "), sep="") - cat("]; ybest (prog=", m2, ", ep=", epbest, ", since=", since, ")\n", sep="") - } - - ## update GP fits - updateGPsep(fgpi, xnext_unit, (obj[k]-fmean)/fsd, verb = verb-2) - df = mleGPsep(fgpi, param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - - for(j in 1:nc) { - updateGPsep(Cgpi[j], xnext_unit, C_bilog[k,j], verb = verb-2) - dc[j,] = mleGPsep(Cgpi[j], param = "d", tmin = dlim[1], tmax = dlim[2], ab = ab, verb=verb-1)$d - } - } - - ## delete GP surrogates - deleteGPsep(fgpi) - for(j in 1:nc) deleteGPsep(Cgpi[j]) - - return(list(prog = prog, xbest = xbest, - obj = obj, C=C, X = X, - feasibility=feasibility, rho=rho, AF_time=AF_time)) -} diff --git a/R/bilog.R b/R/bilog.R index 4045470..391bd34 100644 --- a/R/bilog.R +++ b/R/bilog.R @@ -1,4 +1,4 @@ -#' @title Bilog-transform outcomes +#' @title (un-)bilog-transform outcomes #' #' @description #' The Bilog transform is useful for modeling outcome constraints @@ -6,20 +6,22 @@ #' #' @param Y A training targets #' -#' @returns The transformed outcome observations. -#' -#' @seealso \code{\link[EPBO]{bilog}} +#' @returns The (un-)transformed outcome observations. #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} #' #' @references Eriksson, D., & Poloczek, M. (2021, March). Scalable constrained Bayesian optimization. #' In \emph{International Conference on Artificial Intelligence and Statistics} (pp. 730-738). PMLR. #' -#' @export -#' bilog = function(Y) { Y_tf = sign(Y) * log(abs(Y) + 1) return(Y_tf) +} + +unbilog = function(Y) +{ + Y_utf = sign(Y) * (exp(abs(Y)) - 1) + return(Y_utf) } \ No newline at end of file diff --git a/R/normalize.R b/R/normalize.R index 89231f8..27e83e5 100644 --- a/R/normalize.R +++ b/R/normalize.R @@ -1,28 +1,20 @@ -#' @title Min-max normalize function +#' @title Min-max normalize and un-normalize function #' #' @description Min-max normalize X w.r.t. the provided bounds. #' -#' @param X \code{matrix} +#' @param X \code{matrix} the candidate points #' #' @param bounds 2-column \code{matrix} describing the bounding box. #' the first column gives lower bounds and the second gives upper bounds #' -#' @returns Min-max normalize X w.r.t. the provided bounds. -#' If all elements of \code{X} are contained within \code{bounds}, -#' the normalized values will be contained within \code{[0, 1]^d}. -#' -#' @seealso \code{\link[EPBO]{unnormalize}} +#' @returns Min-max normalize or un-normalized X w.r.t. the provided bounds. +#' The normalized values will be contained within \code{[0, 1]^d}, +#' if all elements of \code{X} are contained within \code{bounds}. +#' the un-normalized values will be contained within \code{bounds}, +#' if all elements of \code{X} are contained within \code{[0, 1]^d}, #' #' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' -#' @export -#' -#' -#' @examples -#' X = matrix(rnorm(6), ncol=2) -#' bounds = cbind(c(-1,-1), c(1,1)) -#' X_normalized = normalize(X, bounds) + normalize = function(X, bounds) { @@ -31,3 +23,11 @@ normalize = function(X, bounds) Xnorm = t(Xnorm) return(Xnorm) } + +unnormalize = function(X, bounds) +{ + if(is.null(nrow(X))) X = matrix(X, nrow=1) + Xnorm = t(X) * (bounds[,2] - bounds[,1]) + bounds[,1] + Xnorm = t(Xnorm) + return(Xnorm) +} diff --git a/R/unbilog.R b/R/unbilog.R deleted file mode 100644 index 1ee4cfa..0000000 --- a/R/unbilog.R +++ /dev/null @@ -1,34 +0,0 @@ -#' @title Un-transform bilog-transformed outcomes -#' -#' @description -#' The Bilog transform is useful for modeling outcome constraints -#' as it magnifies values near zero and flattens extreme values. -#' -#' @param Y A bilog-transfomred targets -#' -#' @returns The un-transformed outcome observations. -#' -#' @seealso \code{\link[EPBO]{bilog}} -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' @references Eriksson, D., & Poloczek, M. (2021, March). Scalable constrained Bayesian optimization. -#' In \emph{International Conference on Artificial Intelligence and Statistics} (pp. 730-738). PMLR. -#' -#' @export -#' - -unbilog = function(Y) -{ - Y_utf = sign(Y) * (exp(abs(Y)) - 1) - return(Y_utf) -} - - - - - - - - - diff --git a/R/unnormalize.R b/R/unnormalize.R deleted file mode 100644 index fd8719a..0000000 --- a/R/unnormalize.R +++ /dev/null @@ -1,33 +0,0 @@ -#' @title Un-normalize min-max function. -#' -#' @description Un-normalizes X w.r.t. the provided bounds. -#' -#' @param X \code{matrix} -#' -#' @param bounds 2-column \code{matrix} describing the bounding box. -#' the first column gives lower bounds and the second gives upper bounds -#' -#' @returns Un-normalized X w.r.t. the provided bounds. -#' If all elements of \code{X} are contained within \code{[0, 1]^d}, -#' the un-normalized values will be contained within \code{bounds}. -#' -#' @seealso \code{\link[EPBO]{normalize}} -#' -#' @author Jiangyan Zhao \email{zhaojy2017@126.com} -#' -#' -#' @export -#' -#' -#' @examples -#' X_normalized = matrix(runif(6), ncol=2) -#' bounds = cbind(c(-1,-1), c(1,1)) -#' X = unnormalize(X_normalized, bounds) - -unnormalize = function(X, bounds) -{ - if(is.null(nrow(X))) X = matrix(X, nrow=1) - Xnorm = t(X) * (bounds[,2] - bounds[,1]) + bounds[,1] - Xnorm = t(Xnorm) - return(Xnorm) -} diff --git a/man/AF_AE.Rd b/man/AF_AE.Rd index ea2927f..1dfd5be 100644 --- a/man/AF_AE.Rd +++ b/man/AF_AE.Rd @@ -2,44 +2,43 @@ % Please edit documentation in R/AF_AE.R \name{AF_AE} \alias{AF_AE} -\title{Asymmetric Entropy} +\title{Asymmetric entropy acquisition function} \usage{ AF_AE(x, fgpi, fnorm, Cgpi, Cnorm, fmin, alpha1 = 1, alpha2 = 5, omega = 2/3) } \arguments{ -\item{x}{description} +\item{x}{a vector containing a single candidate point; or a \code{matrix} with +multiple candidate points} -\item{fgpi}{description} +\item{fgpi}{the GP surrogate model of the objective function} -\item{fnorm}{description} +\item{fnorm}{the maximum of the objective} -\item{Cgpi}{description} +\item{Cgpi}{the GP surrogate models of the constraints} -\item{Cnorm}{description} +\item{Cnorm}{the maxima of the constraints} -\item{fmin}{description} +\item{fmin}{the best objective value obtained so far} -\item{alpha1}{description} +\item{alpha1}{a specified weight for the EI part} -\item{alpha2}{description} +\item{alpha2}{a specified weight for the AE part} -\item{omega}{description} +\item{omega}{a mode location parameter} } \value{ -AF +The AE value(s) at \code{x}. } \description{ -Asymmetric Entropy -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - +The asymmetric entropy (AE) acquisition function of the AE method } \references{ Lindberg, D. V. and H. K. Lee (2015). Optimization under constraints by applying an asymmetric entropy measure. \emph{Journal of Computational and Graphical Statistics} 24(2), 379-393. } +\seealso{ +\code{\link[EPBO]{AF_ScaledEI}}, \code{\link[EPBO]{AF_EY}}, \code{\link[EPBO]{AF_OOSS}} +} \author{ Jiangyan Zhao \email{zhaojy2017@126.com} } diff --git a/man/AF_ConVar.Rd b/man/AF_ConVar.Rd deleted file mode 100644 index 5f6cd0b..0000000 --- a/man/AF_ConVar.Rd +++ /dev/null @@ -1,33 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_ConVar.R -\name{AF_ConVar} -\alias{AF_ConVar} -\title{Constained variance acquisition function} -\usage{ -AF_ConVar(x, fgpi, Cgpi, equal) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Predictive mean -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_EY.Rd b/man/AF_EY.Rd index f2645ac..fc0dc4b 100644 --- a/man/AF_EY.Rd +++ b/man/AF_EY.Rd @@ -7,32 +7,31 @@ AF_EY(x, fgpi, fmean, fsd, Cgpi, rho, equal) } \arguments{ -\item{x}{description} +\item{x}{a vector containing a single candidate point; or a \code{matrix} with +multiple candidate points} -\item{fgpi}{description} +\item{fgpi}{the GP surrogate model of the objective function} -\item{fmean}{description} +\item{fmean}{the mean of the objective value} -\item{fsd}{description} +\item{fsd}{the standard deviation of the objective value} -\item{Cgpi}{description} +\item{Cgpi}{the GP surrogate models of the constraints} -\item{rho}{description} +\item{rho}{the penalty parameters} \item{equal}{an optional vector containing zeros and ones, whose length equals the number of constraints, specifying which should be treated as equality constraints (\code{1}) and which as inequality (\code{0})} } \value{ -AF +The Predictive Mean at \code{x}. } \description{ -Predictive mean +The predictive mean acquisition function of the EPBO method } -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - +\seealso{ +\code{\link[EPBO]{AF_ScaledEI}}, \code{\link[EPBO]{AF_OOSS}}, \code{\link[EPBO]{AF_AE}} } \author{ Jiangyan Zhao \email{zhaojy2017@126.com} diff --git a/man/AF_EY_Kernel.Rd b/man/AF_EY_Kernel.Rd deleted file mode 100644 index 20a00a4..0000000 --- a/man/AF_EY_Kernel.Rd +++ /dev/null @@ -1,35 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_EY_Kernel.R -\name{AF_EY_Kernel} -\alias{AF_EY_Kernel} -\title{Predictive mean acquisition function} -\usage{ -AF_EY_Kernel(x, fgpi, fmean, fsd, Cgpi, rho, equal, type = "UK") -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Predictive mean -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_LCB.Rd b/man/AF_LCB.Rd deleted file mode 100644 index 93634a0..0000000 --- a/man/AF_LCB.Rd +++ /dev/null @@ -1,39 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_LCB.R -\name{AF_LCB} -\alias{AF_LCB} -\title{ScaledEI acquisition function} -\usage{ -AF_LCB(x, fgpi, fmean, fsd, Cgpi, rho, equal) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_OOSS.Rd b/man/AF_OOSS.Rd index ce223c5..77d5c48 100644 --- a/man/AF_OOSS.Rd +++ b/man/AF_OOSS.Rd @@ -2,36 +2,35 @@ % Please edit documentation in R/AF_OOSS.R \name{AF_OOSS} \alias{AF_OOSS} -\title{Asymmetric Entropy} +\title{One over sigma squared (OOSS) acquisition function} \usage{ AF_OOSS(x, fgpi, fnorm, Cgpi, Cnorm) } \arguments{ -\item{x}{description} +\item{x}{a vector containing a single candidate point; or a \code{matrix} with +multiple candidate points} -\item{fgpi}{description} +\item{fgpi}{the GP surrogate model of the objective function} -\item{fnorm}{description} +\item{fnorm}{the maximum of the objective} -\item{Cgpi}{description} +\item{Cgpi}{the GP surrogate models of the constraints} -\item{Cnorm}{description} +\item{Cnorm}{the maxima of the constraints} } \value{ -AF +The OOSS value(s) at \code{x}. } \description{ -Asymmetric Entropy -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - +The one over sigma squared (OOSS) acquisition function of the barrier method } \references{ Pourmohamad, T. and H. K. H. Lee (2022). Bayesian optimization via barrier functions. \emph{Journal of Computational and Graphical Statistics} 31(1), 74-83. } +\seealso{ +\code{\link[EPBO]{AF_ScaledEI}}, \code{\link[EPBO]{AF_EY}}, \code{\link[EPBO]{AF_AE}} +} \author{ Jiangyan Zhao \email{zhaojy2017@126.com} } diff --git a/man/AF_ScaledEI.Rd b/man/AF_ScaledEI.Rd index 2381b4e..e1fdf55 100644 --- a/man/AF_ScaledEI.Rd +++ b/man/AF_ScaledEI.Rd @@ -7,39 +7,38 @@ AF_ScaledEI(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal) } \arguments{ -\item{x}{description} +\item{x}{a vector containing a single candidate point; or a \code{matrix} with +multiple candidate points} -\item{fgpi}{description} +\item{fgpi}{the GP surrogate model of the objective function} -\item{fmean}{description} +\item{fmean}{the mean of the objective value} -\item{fsd}{description} +\item{fsd}{the standard deviation of the objective value} -\item{Cgpi}{description} +\item{Cgpi}{the GP surrogate models of the constraints} -\item{epbest}{description} +\item{epbest}{the best exact penalty value obtained so far} -\item{rho}{description} +\item{rho}{the penalty parameters} \item{equal}{an optional vector containing zeros and ones, whose length equals the number of constraints, specifying which should be treated as equality constraints (\code{1}) and which as inequality (\code{0})} } \value{ -AF +The ScaledEI value(s) at \code{x}. } \description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - +The Scaled expected improvement acquisition function of the EPBO method } \references{ Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function for Bayesian optimization. \emph{arXiv:1808.06918}. } +\seealso{ +\code{\link[EPBO]{AF_OOSS}}, \code{\link[EPBO]{AF_EY}}, \code{\link[EPBO]{AF_AE}} +} \author{ Jiangyan Zhao \email{zhaojy2017@126.com} } diff --git a/man/AF_ScaledEI_Kernel.Rd b/man/AF_ScaledEI_Kernel.Rd deleted file mode 100644 index 60754a7..0000000 --- a/man/AF_ScaledEI_Kernel.Rd +++ /dev/null @@ -1,45 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_ScaledEI_Kernel.R -\name{AF_ScaledEI_Kernel} -\alias{AF_ScaledEI_Kernel} -\title{ScaledEI acquisition function} -\usage{ -AF_ScaledEI_Kernel(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, type = "UK") -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{fmean}{description} - -\item{fsd}{description} - -\item{Cgpi}{description} - -\item{epbest}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_ScaledEI_MC.Rd b/man/AF_ScaledEI_MC.Rd deleted file mode 100644 index 4c3b7fe..0000000 --- a/man/AF_ScaledEI_MC.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_ScaledEI_MC.R -\name{AF_ScaledEI_MC} -\alias{AF_ScaledEI_MC} -\title{ScaledEI acquisition function} -\usage{ -AF_ScaledEI_MC(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N = 1000) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{fmean}{description} - -\item{fsd}{description} - -\item{Cgpi}{description} - -\item{epbest}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} - -\item{N}{positive scalar integer indicating the number of Monte Carlo samples to be used -for calculating ScaledEI} -} -\value{ -ScaledEI -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_TS.Rd b/man/AF_TS.Rd deleted file mode 100644 index 0c2ed3a..0000000 --- a/man/AF_TS.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_TS.R -\name{AF_TS} -\alias{AF_TS} -\title{ScaledEI acquisition function} -\usage{ -AF_TS(x, fgpi, Cgpi, rho, equal, batch_size = 1) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} - -\item{batch_size}{aaa} -} -\value{ -ScaledEI -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_UEI.Rd b/man/AF_UEI.Rd deleted file mode 100644 index bb73593..0000000 --- a/man/AF_UEI.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_UEI.R -\name{AF_UEI} -\alias{AF_UEI} -\title{UEI acquisition function} -\usage{ -AF_UEI(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, beta = 3) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{epbest}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_UEI_Kernel.Rd b/man/AF_UEI_Kernel.Rd deleted file mode 100644 index 50e5eb9..0000000 --- a/man/AF_UEI_Kernel.Rd +++ /dev/null @@ -1,52 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_UEI_Kernel.R -\name{AF_UEI_Kernel} -\alias{AF_UEI_Kernel} -\title{UEI acquisition function} -\usage{ -AF_UEI_Kernel( - x, - fgpi, - fmean, - fsd, - Cgpi, - epbest, - rho, - equal, - beta = 3, - type = "UK" -) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{epbest}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_UEI_MC.Rd b/man/AF_UEI_MC.Rd deleted file mode 100644 index d8506c0..0000000 --- a/man/AF_UEI_MC.Rd +++ /dev/null @@ -1,49 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_UEI_MC.R -\name{AF_UEI_MC} -\alias{AF_UEI_MC} -\title{UEI acquisition function} -\usage{ -AF_UEI_MC(x, fgpi, fmean, fsd, Cgpi, epbest, rho, equal, N = 1000, beta = 3) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{fmean}{description} - -\item{fsd}{description} - -\item{Cgpi}{description} - -\item{epbest}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} - -\item{N}{description} - -\item{beta}{description} -} -\value{ -AF -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/AF_VEI.Rd b/man/AF_VEI.Rd deleted file mode 100644 index 2152399..0000000 --- a/man/AF_VEI.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AF_VEI.R -\name{AF_VEI} -\alias{AF_VEI} -\title{VEI acquisition function} -\usage{ -AF_VEI(x, fgpi, Cgpi, epbest, rho, equal) -} -\arguments{ -\item{x}{description} - -\item{fgpi}{description} - -\item{Cgpi}{description} - -\item{epbest}{description} - -\item{rho}{description} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} -} -\value{ -AF -} -\description{ -Scaled expected improvement -} -\examples{ -B = rbind(c(0, 1), c(0, 1)) - - -} -\references{ -Noe, U. and D. Husmeier (2018). On a new improvement-based acquisition function -for Bayesian optimization. \emph{arXiv:1808.06918}. -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/bilog.Rd b/man/bilog.Rd index e835137..9b91e99 100644 --- a/man/bilog.Rd +++ b/man/bilog.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/bilog.R \name{bilog} \alias{bilog} -\title{Bilog-transform outcomes} +\title{(un-)bilog-transform outcomes} \usage{ bilog(Y) } @@ -10,7 +10,7 @@ bilog(Y) \item{Y}{A training targets} } \value{ -The transformed outcome observations. +The (un-)transformed outcome observations. } \description{ The Bilog transform is useful for modeling outcome constraints @@ -20,9 +20,6 @@ as it magnifies values near zero and flattens extreme values. Eriksson, D., & Poloczek, M. (2021, March). Scalable constrained Bayesian optimization. In \emph{International Conference on Artificial Intelligence and Statistics} (pp. 730-738). PMLR. } -\seealso{ -\code{\link[EPBO]{bilog}} -} \author{ Jiangyan Zhao \email{zhaojy2017@126.com} } diff --git a/man/normalize.Rd b/man/normalize.Rd index 2e57d7e..a90a743 100644 --- a/man/normalize.Rd +++ b/man/normalize.Rd @@ -2,32 +2,26 @@ % Please edit documentation in R/normalize.R \name{normalize} \alias{normalize} -\title{Min-max normalize function} +\title{Min-max normalize and un-normalize function} \usage{ normalize(X, bounds) } \arguments{ -\item{X}{\code{matrix}} +\item{X}{\code{matrix} the candidate points} \item{bounds}{2-column \code{matrix} describing the bounding box. the first column gives lower bounds and the second gives upper bounds} } \value{ -Min-max normalize X w.r.t. the provided bounds. -If all elements of \code{X} are contained within \code{bounds}, -the normalized values will be contained within \code{[0, 1]^d}. +Min-max normalize or un-normalized X w.r.t. the provided bounds. +The normalized values will be contained within \code{[0, 1]^d}, +if all elements of \code{X} are contained within \code{bounds}. +the un-normalized values will be contained within \code{bounds}, +if all elements of \code{X} are contained within \code{[0, 1]^d}, } \description{ Min-max normalize X w.r.t. the provided bounds. } -\examples{ -X = matrix(rnorm(6), ncol=2) -bounds = cbind(c(-1,-1), c(1,1)) -X_normalized = normalize(X, bounds) -} -\seealso{ -\code{\link[EPBO]{unnormalize}} -} \author{ Jiangyan Zhao \email{zhaojy2017@126.com} } diff --git a/man/optim.AE.Rd b/man/optim.AE.Rd index 9cb221f..ef05166 100644 --- a/man/optim.AE.Rd +++ b/man/optim.AE.Rd @@ -83,18 +83,15 @@ blackbox under optimization: \item{obj }{ vector giving the value of the objective for the input under consideration at each trial } \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } - -@seealso \code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.BM}}, \code{\link[EPBO]{optim.EP}} } \description{ Black-box optimization under inequality constraints via an asymmetric entropy. } \examples{ # search space -B = rbind(c(-2.25, 2.5), c(-2.25, 1.75)) +B = rbind(c(-2.25, 2.5), c(-2.5, 1.75)) # Objective and constraint function for use -MTP = function(x, B=rbind(c(-2.25, 2.5), c(-2.5, 1.75))){ - x = pmin(pmax(B[,1], x), B[,2]) +MTP = function(x){ f = -(cos((x[1]-0.1)*x[2]))^2 - x[1]*sin(3*x[1]+x[2]) t = atan2(x[1], x[2]) c = x[1]^2 + x[2]^2 -((2*cos(t)-1/2*cos(2*t)-1/4*cos(3*t)-1/8*cos(4*t))^2) - ((2*sin(t))^2) @@ -110,6 +107,9 @@ AE$xbest Lindberg, D. V. and H. K. Lee (2015). Optimization under constraints by applying an asymmetric entropy measure. \emph{Journal of Computational and Graphical Statistics} 24(2), 379-393. } +\seealso{ +\code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.BM}}, \code{\link[EPBO]{optim.EP}} +} \author{ Jiangyan Zhao \email{zhaojy2017@126.com} } diff --git a/man/optim.BM.Rd b/man/optim.BM.Rd index 5265b54..8c185e5 100644 --- a/man/optim.BM.Rd +++ b/man/optim.BM.Rd @@ -80,10 +80,9 @@ Black-box optimization under inequality constraints via an barrier functions. } \examples{ # search space -B = rbind(c(-2.25, 2.5), c(-2.25, 1.75)) +B = rbind(c(-2.25, 2.5), c(-2.5, 1.75)) # Objective and constraint function for use -MTP = function(x, B=rbind(c(-2.25, 2.5), c(-2.5, 1.75))){ - x = pmin(pmax(B[,1], x), B[,2]) +MTP = function(x){ f = -(cos((x[1]-0.1)*x[2]))^2 - x[1]*sin(3*x[1]+x[2]) t = atan2(x[1], x[2]) c = x[1]^2 + x[2]^2 -((2*cos(t)-1/2*cos(2*t)-1/4*cos(3*t)-1/8*cos(4*t))^2) - ((2*sin(t))^2) diff --git a/man/optim.EP.MC.Rd b/man/optim.EP.MC.Rd deleted file mode 100644 index 3a8fffd..0000000 --- a/man/optim.EP.MC.Rd +++ /dev/null @@ -1,188 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/EPBO_MC.R -\name{optim.EP.MC} -\alias{optim.EP.MC} -\title{Exact Penalty Bayesian Optimization} -\usage{ -optim.EP.MC( - blackbox, - B, - equal = FALSE, - ethresh = 0.01, - Xstart = NULL, - start = 10, - end = 100, - urate = 10, - rho = NULL, - ncandf = function(k) { - k - }, - dg_start = c(0.1, 1e-06), - dlim = c(1e-04, 10), - N = 1000, - plotprog = FALSE, - verb = 2, - ... -) -} -\arguments{ -\item{blackbox}{blackbox of an input (\code{x}), facilitating vectorization on a -\code{matrix} \code{X} thereof, returning a \code{list} -with elements \code{$obj} containing the (scalar) objective value and \code{$c} -containing a vector of evaluations of the (multiple) constraint function at \code{x}.} - -\item{B}{2-column \code{matrix} describing the bounding box. The number of rows -of the \code{matrix} determines the input dimension (\code{length(x)} in \code{blackbox(x)}); -the first column gives lower bounds and the second gives upper bounds} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} - -\item{ethresh}{a threshold used for equality constraints to determine validity for -progress measures; ignored if there are no equality constraints} - -\item{Xstart}{optional matrix of starting design locations in lieu of, or in addition to, -\code{start} random ones; we recommend \code{nrow(Xstart) + start >= 6}; also must -have \code{ncol(Xstart) = nrow(B)}} - -\item{start}{positive integer giving the number of random starting locations before -sequential design (for optimization) is performed; \code{start >= 6} is -recommended unless \code{Xstart} is non-\code{NULL}; in the current version -the starting locations come from a space-filling design via \code{\link[tgp]{dopt.gp}}} - -\item{end}{positive integer giving the total number of evaluations/trials in the -optimization; must have \code{end > start}} - -\item{urate}{positive integer indicating how many optimization trials should pass before -each MLE/MAP update is performed for GP correlation lengthscale -parameter(s)} - -\item{rho}{positive vector initial exact penalty parameter in the exact penalty function; -the default setting of \code{rho = NULL} causes an automatic starting value to be chosen} - -\item{ncandf}{function taking a single integer indicating the optimization trial number \code{t}, where -\code{start < t <= end}, and returning the number of search candidates (e.g., for -expected improvement calculations) at round \code{t}; the default setting -allows the number of candidates to grow linearly with \code{t}} - -\item{dg_start}{2-vector giving starting values for the lengthscale and nugget parameters -of the GP surrogate model(s) for constraints} - -\item{dlim}{2-vector giving bounds for the lengthscale parameter(s) under MLE inference} - -\item{N}{positive scalar integer indicating the number of Monte Carlo samples to be used -for calculating ScaledEI} - -\item{verb}{a non-negative integer indicating the verbosity level; the larger the value the -more that is printed to the screen} - -\item{...}{additional arguments passed to \code{blackbox}} -} -\value{ -The output is a \code{list} summarizing the progress of the evaluations of the -blackbox under optimization: -\item{prog }{ vector giving the best feasible (\code{g(x) <= 0 && |h(x)| <= ethresh}) value of the objective over the trials } -\item{xbest }{ vector giving the recommended solution} -\item{obj }{ vector giving the value of the objective for the input under consideration at each trial } -\item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} -\item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -\item{idx_v}{index of the valid solutions} -\item{AF_time}{the time required to compute the acquisition function} -} -\description{ -Black-box optimization under mixed equality and inequality constraints via an exact penalty. -} -\details{ -Additional details... -} -\examples{ -## Inequality constrained problem -# search space -B = rbind(c(0, 1), c(0, 1)) -# Objective and constraint function for use -HSQ = function(x){ - B = rbind(c(0, 1), c(0, 1)) - x = pmin(pmax(B[,1], x), B[,2]) - t = function(z) - return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) - herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) - c1 = 1.5 - x[1] - 2*x[2] - 0.5*sin(2*pi*(x[1]^2 - 2*x[2])) - c2 = x[1]^2 + x[2]^2 - 1.5 - return(list(obj=herbtooth, c=cbind(c1,c2))) -} -EPBO = optim.EP(HSQ, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, verb = 0) -# progress, best feasible value of the objective over the trials -EPBO$prog -# the recommended solution -EPBO$xbest - -## Mixed constrained problem -# search space -B = rbind(c(0, 1), c(0, 1)) -# Objective and constraint function for use -goldstein.price = function(X) -{ - if(is.null(nrow(X))) X <- matrix(X, nrow=1) - m <- 8.6928 - s <- 2.4269 - x1 <- 4 * X[,1] - 2 - x2 <- 4 * X[,2] - 2 - a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2) - b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 + 27 * x2^2) - f <- log(a * b) - f <- (f - m)/s - return(f) -} - -toy.c1 <- function(X) { - if (is.null(dim(X))) X <- matrix(X, nrow=1) - c1 = 3/2 - X[,1] - 2*X[,2] - .5*sin(2*pi*(X[,1]^2 - 2*X[,2])) -} - -branin.c = function(X){ - if (is.null(dim(X))) X <- matrix(X, nrow=1) - x1 = X[,1] * 15 - 5 - x2 = X[,2] * 15 - f = (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10 - return(-f+25) -} - -parr <- function(X){ - if (is.null(dim(X))) X <- matrix(X, nrow=1) - x1 <- (2 * X[,1] - 1) - x2 <- (2 * X[,2] - 1) - g <- (4-2.1*x1^2+1/3*x1^4)*x1^2 + x1*x2 + (-4+4*x2^2)*x2^2 + - 3*sin(6*(1-x1)) + 3*sin(6*(1-x2)) - return(-g+4) -} - -gsbp.constraints <- function(x){ - return(cbind(toy.c1(x), branin.c(x), parr(x))) -} - -# problem definition -gsbpprob <- function(X, known.only=FALSE) -{ - if(is.null(nrow(X))) X <- matrix(X, nrow=1) - if(known.only) stop("known.only not supported for this example") - f = goldstein.price(X) - C = gsbp.constraints(X) - return(list(obj=f, c=cbind(C[,1], C[,2]/100, C[,3]/10))) -} -EPBO = optim.EP(gsbpprob, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, - ethresh = 1e-2, equal = c(0,1,1), verb = 0) - -# progress, best feasible value of the objective over the trials -EPBO$prog -# the recommended solution -EPBO$xbest -} -\seealso{ -\code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} -\keyword{design} -\keyword{optimize} diff --git a/man/optim.EP.Rd b/man/optim.EP.Rd index 38d45fa..4cb6663 100644 --- a/man/optim.EP.Rd +++ b/man/optim.EP.Rd @@ -17,18 +17,17 @@ optim.EP( ncandf = function(k) { k }, - dg_start = c(0.1, 1e-06), - dlim = c(1e-04, 10), + dg_start = c(0.01 * sqrt(nrow(B)), 1e-06), + dlim = c(1e-04, 1) * sqrt(nrow(B)), plotprog = FALSE, ey_tol = 0.01, - AF_tol = sqrt(.Machine$double.eps), verb = 2, ... ) } \arguments{ \item{blackbox}{blackbox of an input (\code{x}), facilitating vectorization on a -\code{matrix} \code{X} thereof, returning a \code{list} +\code{matrix} \code{X} thereof, returning a \code{list} with elements \code{$obj} containing the (scalar) objective value and \code{$c} containing a vector of evaluations of the (multiple) constraint function at \code{x}.} @@ -59,8 +58,8 @@ optimization; must have \code{end > start}} each MLE/MAP update is performed for GP correlation lengthscale parameter(s)} -\item{rho}{positive vector initial exact penalty parameter in the exact penalty function; -the default setting of \code{rho = NULL} causes an automatic starting value to be chosen} +\item{rho}{positive vector initial exact penalty parameters in the exact penalty function; +the default setting of \code{rho = NULL} causes the automatic starting values to be chosen} \item{ncandf}{function taking a single integer indicating the optimization trial number \code{t}, where \code{start < t <= end}, and returning the number of search candidates (e.g., for @@ -68,10 +67,19 @@ expected improvement calculations) at round \code{t}; the default setting allows the number of candidates to grow linearly with \code{t}} \item{dg_start}{2-vector giving starting values for the lengthscale and nugget parameters -of the GP surrogate model(s) for constraints} +of the GP surrogate model(s) for the objective and constraints} \item{dlim}{2-vector giving bounds for the lengthscale parameter(s) under MLE inference} +\item{plotprog}{\code{logical} indicating if progress plots should be made after each inner iteration; +the plots show four panels tracking the best valid objective, the ScaledEI or EY surface, +the predictive mean and standard deviation of the objective function, over the first two input variables} + +\item{ey_tol}{a scalar proportion indicating how many of the ScaledEIs +at \code{ncandf(t)} candidate locations must be non-zero to \dQuote{trust} +that metric to guide search, reducing to an EY-based search instead +(choosing that proportion to be one forces EY-based search)} + \item{verb}{a non-negative integer indicating the verbosity level; the larger the value the more that is printed to the screen} @@ -85,14 +93,14 @@ blackbox under optimization: \item{obj }{ vector giving the value of the objective for the input under consideration at each trial } \item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} \item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -\item{idx_v}{index of the valid solutions} -\item{AF_time}{the time required to compute the acquisition function} +\item{feasibility}{vector giving the feasibility for the input under consideration at each trial } +\item{rho}{vector giving the penalty parameters } } \description{ Black-box optimization under mixed equality and inequality constraints via an exact penalty. } \details{ -Additional details... +To be added. } \examples{ ## Inequality constrained problem @@ -100,8 +108,6 @@ Additional details... B = rbind(c(0, 1), c(0, 1)) # Objective and constraint function for use HSQ = function(x){ - B = rbind(c(0, 1), c(0, 1)) - x = pmin(pmax(B[,1], x), B[,2]) t = function(z) return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) diff --git a/man/optim.EP.kernel.Rd b/man/optim.EP.kernel.Rd deleted file mode 100644 index 7b482a6..0000000 --- a/man/optim.EP.kernel.Rd +++ /dev/null @@ -1,183 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/EPBO_Kernel.R -\name{optim.EP.kernel} -\alias{optim.EP.kernel} -\title{Exact Penalty Bayesian Optimization} -\usage{ -optim.EP.kernel( - blackbox, - B, - equal = FALSE, - ethresh = 0.01, - Xstart = NULL, - start = 10, - end = 100, - urate = 10, - rho = NULL, - ncandf = function(k) { - k - }, - kmcontrol = list(formula = ~1, covtype = "matern5_2", nugget = 1e-06, parinit = - rep(0.1, nrow(B)), lower = rep(1e-04, nrow(B)), upper = rep(10, nrow(B))), - verb = 2, - ... -) -} -\arguments{ -\item{blackbox}{blackbox of an input (\code{x}), facilitating vectorization on a -\code{matrix} \code{X} thereof, returning a \code{list} -with elements \code{$obj} containing the (scalar) objective value and \code{$c} -containing a vector of evaluations of the (multiple) constraint function at \code{x}.} - -\item{B}{2-column \code{matrix} describing the bounding box. The number of rows -of the \code{matrix} determines the input dimension (\code{length(x)} in \code{blackbox(x)}); -the first column gives lower bounds and the second gives upper bounds} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} - -\item{ethresh}{a threshold used for equality constraints to determine validity for -progress measures; ignored if there are no equality constraints} - -\item{Xstart}{optional matrix of starting design locations in lieu of, or in addition to, -\code{start} random ones; we recommend \code{nrow(Xstart) + start >= 6}; also must -have \code{ncol(Xstart) = nrow(B)}} - -\item{start}{positive integer giving the number of random starting locations before -sequential design (for optimization) is performed; \code{start >= 6} is -recommended unless \code{Xstart} is non-\code{NULL}; in the current version -the starting locations come from a space-filling design via \code{\link[tgp]{dopt.gp}}} - -\item{end}{positive integer giving the total number of evaluations/trials in the -optimization; must have \code{end > start}} - -\item{urate}{positive integer indicating how many optimization trials should pass before -each MLE/MAP update is performed for GP correlation lengthscale -parameter(s)} - -\item{rho}{positive vector initial exact penalty parameter in the exact penalty function; -the default setting of \code{rho = NULL} causes an automatic starting value to be chosen} - -\item{ncandf}{function taking a single integer indicating the optimization trial number \code{t}, where -\code{start < t <= end}, and returning the number of search candidates (e.g., for -expected improvement calculations) at round \code{t}; the default setting -allows the number of candidates to grow linearly with \code{t}} - -\item{verb}{a non-negative integer indicating the verbosity level; the larger the value the -more that is printed to the screen} - -\item{...}{additional arguments passed to \code{blackbox}} - -\item{dg_start}{2-vector giving starting values for the lengthscale and nugget parameters -of the GP surrogate model(s) for constraints} - -\item{dlim}{2-vector giving bounds for the lengthscale parameter(s) under MLE inference} -} -\value{ -The output is a \code{list} summarizing the progress of the evaluations of the -blackbox under optimization: -\item{prog }{ vector giving the best feasible (\code{g(x) <= 0 && |h(x)| <= ethresh}) value of the objective over the trials } -\item{xbest }{ vector giving the recommended solution} -\item{obj }{ vector giving the value of the objective for the input under consideration at each trial } -\item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} -\item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -\item{idx_v}{index of the valid solutions} -\item{AF_time}{the time required to compute the acquisition function} -} -\description{ -Black-box optimization under mixed equality and inequality constraints via an exact penalty. -} -\details{ -Additional details... -} -\examples{ -## Inequality constrained problem -# search space -B = rbind(c(0, 1), c(0, 1)) -# Objective and constraint function for use -HSQ = function(x){ - B = rbind(c(0, 1), c(0, 1)) - x = pmin(pmax(B[,1], x), B[,2]) - t = function(z) - return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) - herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) - c1 = 1.5 - x[1] - 2*x[2] - 0.5*sin(2*pi*(x[1]^2 - 2*x[2])) - c2 = x[1]^2 + x[2]^2 - 1.5 - return(list(obj=herbtooth, c=cbind(c1,c2))) -} -EPBO = optim.EP(HSQ, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, verb = 0) -# progress, best feasible value of the objective over the trials -EPBO$prog -# the recommended solution -EPBO$xbest - -## Mixed constrained problem -# search space -B = rbind(c(0, 1), c(0, 1)) -# Objective and constraint function for use -goldstein.price = function(X) -{ - if(is.null(nrow(X))) X <- matrix(X, nrow=1) - m <- 8.6928 - s <- 2.4269 - x1 <- 4 * X[,1] - 2 - x2 <- 4 * X[,2] - 2 - a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2) - b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 + 27 * x2^2) - f <- log(a * b) - f <- (f - m)/s - return(f) -} - -toy.c1 <- function(X) { - if (is.null(dim(X))) X <- matrix(X, nrow=1) - c1 = 3/2 - X[,1] - 2*X[,2] - .5*sin(2*pi*(X[,1]^2 - 2*X[,2])) -} - -branin.c = function(X){ - if (is.null(dim(X))) X <- matrix(X, nrow=1) - x1 = X[,1] * 15 - 5 - x2 = X[,2] * 15 - f = (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10 - return(-f+25) -} - -parr <- function(X){ - if (is.null(dim(X))) X <- matrix(X, nrow=1) - x1 <- (2 * X[,1] - 1) - x2 <- (2 * X[,2] - 1) - g <- (4-2.1*x1^2+1/3*x1^4)*x1^2 + x1*x2 + (-4+4*x2^2)*x2^2 + - 3*sin(6*(1-x1)) + 3*sin(6*(1-x2)) - return(-g+4) -} - -gsbp.constraints <- function(x){ - return(cbind(toy.c1(x), branin.c(x), parr(x))) -} - -# problem definition -gsbpprob <- function(X, known.only=FALSE) -{ - if(is.null(nrow(X))) X <- matrix(X, nrow=1) - if(known.only) stop("known.only not supported for this example") - f = goldstein.price(X) - C = gsbp.constraints(X) - return(list(obj=f, c=cbind(C[,1], C[,2]/100, C[,3]/10))) -} -EPBO = optim.EP(gsbpprob, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, - ethresh = 1e-2, equal = c(0,1,1), verb = 0) - -# progress, best feasible value of the objective over the trials -EPBO$prog -# the recommended solution -EPBO$xbest -} -\seealso{ -\code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} -\keyword{design} -\keyword{optimize} diff --git a/man/optim.EP4HD.Rd b/man/optim.EP4HD.Rd deleted file mode 100644 index 8c5aaaf..0000000 --- a/man/optim.EP4HD.Rd +++ /dev/null @@ -1,190 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/EPBO4HD.R -\name{optim.EP4HD} -\alias{optim.EP4HD} -\title{Exact Penalty Bayesian Optimization} -\usage{ -optim.EP4HD( - blackbox, - B, - equal = FALSE, - ethresh = 0.01, - Xstart = NULL, - start = 10, - end = 100, - urate = 10, - rho = NULL, - dg_start = c(0.1, 1e-06), - dlim = c(1e-04, 10), - ncandf = function(k) { - k - }, - batch_size = 1, - trcontrol = list(length = 0.8, length_min = 0.5^7, length_max = 1.6, failure_counter = - 0, failure_tolerance = max(4, nrow(B)), success_counter = 0, success_tolerance = 3), - plotprog = FALSE, - opt = TRUE, - ey.tol = 0.01, - AF.tol = sqrt(.Machine$double.eps), - verb = 2, - ... -) -} -\arguments{ -\item{blackbox}{blackbox of an input (\code{x}), facilitating vectorization on a -\code{matrix} \code{X} thereof, returning a \code{list} -with elements \code{$obj} containing the (scalar) objective value and \code{$c} -containing a vector of evaluations of the (multiple) constraint function at \code{x}.} - -\item{B}{2-column \code{matrix} describing the bounding box. The number of rows -of the \code{matrix} determines the input dimension (\code{length(x)} in \code{blackbox(x)}); -the first column gives lower bounds and the second gives upper bounds} - -\item{equal}{an optional vector containing zeros and ones, whose length equals the number of -constraints, specifying which should be treated as equality constraints (\code{1}) and -which as inequality (\code{0})} - -\item{ethresh}{a threshold used for equality constraints to determine validity for -progress measures; ignored if there are no equality constraints} - -\item{Xstart}{optional matrix of starting design locations in lieu of, or in addition to, -\code{start} random ones; we recommend \code{nrow(Xstart) + start >= 6}; also must -have \code{ncol(Xstart) = nrow(B)}} - -\item{start}{positive integer giving the number of random starting locations before -sequential design (for optimization) is performed; \code{start >= 6} is -recommended unless \code{Xstart} is non-\code{NULL}; in the current version -the starting locations come from a space-filling design via \code{\link[tgp]{dopt.gp}}} - -\item{end}{positive integer giving the total number of evaluations/trials in the -optimization; must have \code{end > start}} - -\item{urate}{positive integer indicating how many optimization trials should pass before -each MLE/MAP update is performed for GP correlation lengthscale -parameter(s)} - -\item{rho}{positive vector initial exact penalty parameter in the exact penalty function; -the default setting of \code{rho = NULL} causes an automatic starting value to be chosen} - -\item{dg_start}{2-vector giving starting values for the lengthscale and nugget parameters -of the GP surrogate model(s) for constraints} - -\item{dlim}{2-vector giving bounds for the lengthscale parameter(s) under MLE inference} - -\item{ncandf}{function taking a single integer indicating the optimization trial number \code{t}, where -\code{start < t <= end}, and returning the number of search candidates (e.g., for -expected improvement calculations) at round \code{t}; the default setting -allows the number of candidates to grow linearly with \code{t}} - -\item{verb}{a non-negative integer indicating the verbosity level; the larger the value the -more that is printed to the screen} - -\item{...}{additional arguments passed to \code{blackbox}} -} -\value{ -The output is a \code{list} summarizing the progress of the evaluations of the -blackbox under optimization: -\item{prog }{ vector giving the best feasible (\code{g(x) <= 0 && |h(x)| <= ethresh}) value of the objective over the trials } -\item{xbest }{ vector giving the recommended solution} -\item{obj }{ vector giving the value of the objective for the input under consideration at each trial } -\item{C }{ \code{matrix} giving the value of the constraint function for the input under consideration at each trial} -\item{X }{ \code{matrix} giving the input values at which the blackbox function was evaluated } -\item{idx_v}{index of the valid solutions} -\item{AF_time}{the time required to compute the acquisition function} -} -\description{ -Black-box optimization under mixed equality and inequality constraints via an exact penalty. -} -\details{ -Additional details... -} -\examples{ -## Inequality constrained problem -# search space -B = rbind(c(0, 1), c(0, 1)) -# Objective and constraint function for use -HSQ = function(x){ - B = rbind(c(0, 1), c(0, 1)) - x = pmin(pmax(B[,1], x), B[,2]) - t = function(z) - return(exp(-(z-1)^2) + exp(-0.8*(z+1)^2) - 0.05*sin(8*(z+0.1))) - herbtooth = -t(4*x[1]-2)*t(4*x[2]-2) - c1 = 1.5 - x[1] - 2*x[2] - 0.5*sin(2*pi*(x[1]^2 - 2*x[2])) - c2 = x[1]^2 + x[2]^2 - 1.5 - return(list(obj=herbtooth, c=cbind(c1,c2))) -} -EPBO = optim.EP(HSQ, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, verb = 0) -# progress, best feasible value of the objective over the trials -EPBO$prog -# the recommended solution -EPBO$xbest - -## Mixed constrained problem -# search space -B = rbind(c(0, 1), c(0, 1)) -# Objective and constraint function for use -goldstein.price = function(X) -{ - if(is.null(nrow(X))) X <- matrix(X, nrow=1) - m <- 8.6928 - s <- 2.4269 - x1 <- 4 * X[,1] - 2 - x2 <- 4 * X[,2] - 2 - a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2) - b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 + 27 * x2^2) - f <- log(a * b) - f <- (f - m)/s - return(f) -} - -toy.c1 <- function(X) { - if (is.null(dim(X))) X <- matrix(X, nrow=1) - c1 = 3/2 - X[,1] - 2*X[,2] - .5*sin(2*pi*(X[,1]^2 - 2*X[,2])) -} - -branin.c = function(X){ - if (is.null(dim(X))) X <- matrix(X, nrow=1) - x1 = X[,1] * 15 - 5 - x2 = X[,2] * 15 - f = (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10 - return(-f+25) -} - -parr <- function(X){ - if (is.null(dim(X))) X <- matrix(X, nrow=1) - x1 <- (2 * X[,1] - 1) - x2 <- (2 * X[,2] - 1) - g <- (4-2.1*x1^2+1/3*x1^4)*x1^2 + x1*x2 + (-4+4*x2^2)*x2^2 + - 3*sin(6*(1-x1)) + 3*sin(6*(1-x2)) - return(-g+4) -} - -gsbp.constraints <- function(x){ - return(cbind(toy.c1(x), branin.c(x), parr(x))) -} - -# problem definition -gsbpprob <- function(X, known.only=FALSE) -{ - if(is.null(nrow(X))) X <- matrix(X, nrow=1) - if(known.only) stop("known.only not supported for this example") - f = goldstein.price(X) - C = gsbp.constraints(X) - return(list(obj=f, c=cbind(C[,1], C[,2]/100, C[,3]/10))) -} -EPBO = optim.EP(gsbpprob, B, ncandf = function(k){ 1e3 }, start = 20, end = 120, - ethresh = 1e-2, equal = c(0,1,1), verb = 0) - -# progress, best feasible value of the objective over the trials -EPBO$prog -# the recommended solution -EPBO$xbest -} -\seealso{ -\code{\link[laGP]{optim.auglag}}, \code{\link[laGP]{optim.efi}}, \code{\link[EPBO]{optim.AE}}, \code{\link[EPBO]{optim.BM}} -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} -\keyword{design} -\keyword{optimize} diff --git a/man/unbilog.Rd b/man/unbilog.Rd deleted file mode 100644 index 2613a76..0000000 --- a/man/unbilog.Rd +++ /dev/null @@ -1,28 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/unbilog.R -\name{unbilog} -\alias{unbilog} -\title{Un-transform bilog-transformed outcomes} -\usage{ -unbilog(Y) -} -\arguments{ -\item{Y}{A bilog-transfomred targets} -} -\value{ -The un-transformed outcome observations. -} -\description{ -The Bilog transform is useful for modeling outcome constraints -as it magnifies values near zero and flattens extreme values. -} -\references{ -Eriksson, D., & Poloczek, M. (2021, March). Scalable constrained Bayesian optimization. -In \emph{International Conference on Artificial Intelligence and Statistics} (pp. 730-738). PMLR. -} -\seealso{ -\code{\link[EPBO]{bilog}} -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -} diff --git a/man/unnormalize.Rd b/man/unnormalize.Rd deleted file mode 100644 index 66315c5..0000000 --- a/man/unnormalize.Rd +++ /dev/null @@ -1,33 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/unnormalize.R -\name{unnormalize} -\alias{unnormalize} -\title{Un-normalize min-max function.} -\usage{ -unnormalize(X, bounds) -} -\arguments{ -\item{X}{\code{matrix}} - -\item{bounds}{2-column \code{matrix} describing the bounding box. -the first column gives lower bounds and the second gives upper bounds} -} -\value{ -Un-normalized X w.r.t. the provided bounds. -If all elements of \code{X} are contained within \code{[0, 1]^d}, -the un-normalized values will be contained within \code{bounds}. -} -\description{ -Un-normalizes X w.r.t. the provided bounds. -} -\examples{ -X_normalized = matrix(runif(6), ncol=2) -bounds = cbind(c(-1,-1), c(1,1)) -X = unnormalize(X_normalized, bounds) -} -\seealso{ -\code{\link[EPBO]{normalize}} -} -\author{ -Jiangyan Zhao \email{zhaojy2017@126.com} -}