R/max_EI.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

#' Maximization of the Expected Improvement criterion
#' 
#' Given an object of class \code{\link[DiceKriging]{km}} and a set of tuning
#' parameters (\code{lower},\code{upper},\code{parinit}, and \code{control}),
#' \code{max_EI} performs the maximization of the Expected Improvement
#' criterion and delivers the next point to be visited in an EGO-like
#' procedure.
#' 
#' The latter maximization relies on a genetic algorithm using derivatives,
#' \code{\link[rgenoud]{genoud}}. This function plays a central role in the
#' package since it is in constant use in the proposed algorithms. It is
#' important to remark that the information needed about the objective
#' function reduces here to the vector of response values embedded in
#' \code{model} (no call to the objective function or simulator).
#' 
#' The current minimum of the observations can be replaced by an arbitrary
#' value (plugin), which is usefull in particular in noisy frameworks.
#' 
#' 
#' @param model an object of class \code{\link[DiceKriging]{km}} ,
#' @param plugin optional scalar: if provided, it replaces the minimum of the
#' current observations,
#' @param type Kriging type: "SK" or "UK"
#' @param lower vector of lower bounds for the variables to be optimized over,
#' @param upper vector of upper bounds for the variables to be optimized over,
#' @param parinit optional vector of initial values for the variables to be
#' optimized over,
#' @param minimization logical specifying if EI is used in minimiziation or in
#' maximization,
#' @param control optional list of control parameters for optimization.  One
#' can control \code{"pop.size"} (default : [N=3*2^dim for dim<6 and N=32*dim
#' otherwise]), \code{"max.generations"} (12), \code{"wait.generations"} (2)
#' and \code{"BFGSburnin"} (2) of function \code{"genoud"} (see
#' \code{\link[rgenoud]{genoud}}).  Numbers into brackets are the default
#' values
#' @return A list with components:
#' 
#' \item{par}{The best set of parameters found.}
#' 
#' \item{value}{The value of expected improvement at par.}
#' @author David Ginsbourger 
#' 
#' Olivier Roustant
#' 
#' Victor Picheny 
#' @references
#' 
#' D. Ginsbourger (2009), \emph{Multiples metamodeles pour l'approximation et
#' l'optimisation de fonctions numeriques multivariables}, Ph.D. thesis, Ecole
#' Nationale Superieure des Mines de Saint-Etienne, 2009.
#' 
#' D.R. Jones, M. Schonlau, and W.J. Welch (1998), Efficient global
#' optimization of expensive black-box functions, \emph{Journal of Global
#' Optimization}, 13, 455-492.
#' 
#' W.R. Jr. Mebane and J.S. Sekhon (2009), in press, Genetic optimization
#' using derivatives: The rgenoud package for R, \emph{Journal of Statistical
#' Software}.
#' @keywords optimize
#' @examples
#' 
#' set.seed(123)
#' ##########################################################
#' ### "ONE-SHOT" EI-MAXIMIZATION OF THE BRANIN FUNCTION ####
#' ### 	KNOWN AT A 9-POINTS FACTORIAL DESIGN          ####
#' ##########################################################
#' 
#' # a 9-points factorial design, and the corresponding response
#' d <- 2 
#' n <- 9
#' design.fact <- expand.grid(seq(0,1,length=3), seq(0,1,length=3)) 
#' names(design.fact) <- c("x1", "x2")
#' design.fact <- data.frame(design.fact) 
#' names(design.fact) <- c("x1", "x2")
#' response.branin <- apply(design.fact, 1, branin)
#' response.branin <- data.frame(response.branin) 
#' names(response.branin) <- "y" 
#' 
#' # model identification
#' fitted.model1 <- km(~1, design=design.fact, response=response.branin, 
#' covtype="gauss", control=list(pop.size=50,trace=FALSE), parinit=c(0.5, 0.5))
#' 
#' # EGO one step
#' library(rgenoud)
#' lower <- rep(0,d) 
#' upper <- rep(1,d)     # domain for Branin function
#' oEGO <- max_EI(fitted.model1, lower=lower, upper=upper, 
#' control=list(pop.size=20, BFGSburnin=2))
#' print(oEGO)
#' 
#' # graphics
#' n.grid <- 20
#' x.grid <- y.grid <- seq(0,1,length=n.grid)
#' design.grid <- expand.grid(x.grid, y.grid)
#' response.grid <- apply(design.grid, 1, branin)
#' z.grid <- matrix(response.grid, n.grid, n.grid)
#' contour(x.grid,y.grid,z.grid,40)
#' title("Branin Function")
#' points(design.fact[,1], design.fact[,2], pch=17, col="blue")
#' points(oEGO$par[1], oEGO$par[2], pch=19, col="red")
#' 
#' 
#' #############################################################
#' ### "ONE-SHOT" EI-MAXIMIZATION OF THE CAMELBACK FUNCTION ####
#' ###	KNOWN AT A 16-POINTS FACTORIAL DESIGN            ####
#' #############################################################
#' \dontrun{
#' # a 16-points factorial design, and the corresponding response
#' d <- 2 
#' n <- 16
#' design.fact <- expand.grid(seq(0,1,length=4), seq(0,1,length=4)) 
#' names(design.fact)<-c("x1", "x2")
#' design.fact <- data.frame(design.fact) 
#' names(design.fact) <- c("x1", "x2")
#' response.camelback <- apply(design.fact, 1, camelback)
#' response.camelback <- data.frame(response.camelback) 
#' names(response.camelback) <- "y" 
#' 
#' # model identification
#' fitted.model1 <- km(~1, design=design.fact, response=response.camelback, 
#' covtype="gauss", control=list(pop.size=50,trace=FALSE), parinit=c(0.5, 0.5))
#' 
#' # EI maximization
#' library(rgenoud)
#' lower <- rep(0,d) 
#' upper <- rep(1,d)   
#' oEGO <- max_EI(fitted.model1, lower=lower, upper=upper, 
#' control=list(pop.size=20, BFGSburnin=2))
#' print(oEGO)
#' 
#' # graphics
#' n.grid <- 20
#' x.grid <- y.grid <- seq(0,1,length=n.grid)
#' design.grid <- expand.grid(x.grid, y.grid)
#' response.grid <- apply(design.grid, 1, camelback)
#' z.grid <- matrix(response.grid, n.grid, n.grid)
#' contour(x.grid,y.grid,z.grid,40)
#' title("Camelback Function")
#' points(design.fact[,1], design.fact[,2], pch=17, col="blue")
#' points(oEGO$par[1], oEGO$par[2], pch=19, col="red")
#' }
#' 
#' ####################################################################
#' ### "ONE-SHOT" EI-MAXIMIZATION OF THE GOLDSTEIN-PRICE FUNCTION #####
#' ### 	     KNOWN AT A 9-POINTS FACTORIAL DESIGN              #####
#' ####################################################################
#' 
#' \dontrun{
#' # a 9-points factorial design, and the corresponding response
#' d <- 2 
#' n <- 9
#' design.fact <- expand.grid(seq(0,1,length=3), seq(0,1,length=3)) 
#' names(design.fact)<-c("x1", "x2")
#' design.fact <- data.frame(design.fact) 
#' names(design.fact)<-c("x1", "x2")
#' response.goldsteinPrice <- apply(design.fact, 1, goldsteinPrice)
#' response.goldsteinPrice <- data.frame(response.goldsteinPrice) 
#' names(response.goldsteinPrice) <- "y" 
#' 
#' # model identification
#' fitted.model1 <- km(~1, design=design.fact, response=response.goldsteinPrice, 
#' covtype="gauss", control=list(pop.size=50, max.generations=50, 
#' wait.generations=5, BFGSburnin=10, trace=FALSE), parinit=c(0.5, 0.5), optim.method="gen")
#' 
#' # EI maximization
#' library(rgenoud)
#' lower <- rep(0,d); upper <- rep(1,d);     # domain for Branin function
#' oEGO <- max_EI(fitted.model1, lower=lower, upper=upper, control
#' =list(pop.size=50, max.generations=50, wait.generations=5, BFGSburnin=10))
#' print(oEGO)
#' 
#' # graphics
#' n.grid <- 20
#' x.grid <- y.grid <- seq(0,1,length=n.grid)
#' design.grid <- expand.grid(x.grid, y.grid)
#' response.grid <- apply(design.grid, 1, goldsteinPrice)
#' z.grid <- matrix(response.grid, n.grid, n.grid)
#' contour(x.grid,y.grid,z.grid,40)
#' title("Goldstein-Price Function")
#' points(design.fact[,1], design.fact[,2], pch=17, col="blue")
#' points(oEGO$par[1], oEGO$par[2], pch=19, col="red")
#' }
#' 
#' @export max_EI
max_EI <-function(model, plugin=NULL, type = "UK", lower, upper, parinit=NULL, minimization = TRUE, control=NULL) {

  if (is.null(plugin)){ plugin <- min(model@y) }
  
	EI.envir <- new.env()	
	environment(EI) <- environment(EI.grad) <- EI.envir 
	gr = EI.grad
  
	d <- ncol(model@X)

	if (is.null(control$print.level)) control$print.level <- 1
	if(d<=6) N <- 3*2^d else N <- 32*d 
	if (is.null(control$BFGSmaxit)) control$BFGSmaxit <- N
	if (is.null(control$pop.size))  control$pop.size <- N
	if (is.null(control$solution.tolerance))  control$solution.tolerance <- 1e-21
	if (is.null(control$max.generations))  control$max.generations <- 12
	if (is.null(control$wait.generations))  control$wait.generations <- 2
	if (is.null(control$BFGSburnin)) control$BFGSburnin <- 2
	if (is.null(parinit))  parinit <- lower + runif(d)*(upper-lower)
     
	domaine <- cbind(lower, upper)

	o <- genoud(EI, nvars=d, max=TRUE,
	            pop.size=control$pop.size, max.generations=control$max.generations, wait.generations=control$wait.generations,
	            hard.generation.limit=TRUE, starting.values=parinit, MemoryMatrix=TRUE, 
	            Domains=domaine, default.domains=10, solution.tolerance=control$solution.tolerance,
	            gr=gr, boundary.enforcement=2, lexical=FALSE, gradient.check=FALSE, BFGS=TRUE,
	            data.type.int=FALSE, hessian=FALSE, unif.seed=floor(runif(1,max=10000)), int.seed=floor(runif(1,max=10000)), 
	            print.level=control$print.level,  
	            share.type=0, instance.number=0, output.path="stdout", output.append=FALSE, project.path=NULL,
	            P1=50, P2=50, P3=50, P4=50, P5=50, P6=50, P7=50, P8=50, P9=0, P9mix=NULL, 
	            BFGSburnin=control$BFGSburnin, BFGSfn=NULL, BFGShelp=NULL, control=list("maxit"=control$BFGSmaxit), 
	            cluster=FALSE, balance=FALSE, debug=FALSE,
	            model=model, plugin=plugin, type=type,minimization = minimization, envir=EI.envir
	)
  
# 	o <- genoud(EI, nvars=d, max=TRUE, 
# 		pop.size=control$pop.size,
# 		max.generations=control$max.generations, 
# 		wait.generations=control$wait.generations,
#     hard.generation.limit=TRUE, starting.values=parinit, MemoryMatrix=TRUE, 
#     Domains=domaine, default.domains=10, solution.tolerance=0.00001, 
#     gr=gr, boundary.enforcement=2, lexical=FALSE, gradient.check=TRUE, BFGS=TRUE,
#     data.type.int=FALSE, hessian=FALSE, unif.seed=floor(runif(1,max=10000)), int.seed=floor(runif(1,max=10000)),
# 	  print.level=control$print.level, 
#     share.type=0, instance.number=0, output.path="stdout", output.append=FALSE, project.path=NULL,
#     P1=50, P2=50, P3=50, P4=50, P5=50, P6=50, P7=50, P8=50, P9=0,	P9mix=NULL, 
#     BFGSburnin=control$BFGSburnin,BFGSfn=NULL, BFGShelp=NULL,
#     cluster=FALSE, balance=FALSE, debug=TRUE, 
#     model=model, plugin=plugin, type=type, envir=EI.envir 
# 		)
                            
  o$par <- t(as.matrix(o$par))
	colnames(o$par) <- colnames(model@X)
	o$value <- as.matrix(o$value)
	colnames(o$value) <- "EI"  
	return(list(par=o$par, value=o$value)) 
}