R/max_EQI.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

##' Maximizer of the Expected Quantile Improvement criterion function
##' 
##' Maximization, based on the package rgenoud of the Expected Quantile
##' Improvement (EQI) criterion.
##' 
##' 
##' @param model a Kriging model of "km" class
##' @param new.noise.var the (scalar) noise variance of an observation. Default
##' value is 0 (noise-free observation).
##' @param beta Quantile level (default value is 0.9)
##' @param q.min The current best kriging quantile. If not provided, this
##' quantity is evaluated inside the EQI function (may increase computational
##' time).
##' @param type Kriging type: "SK" or "UK"
##' @param lower vector containing the lower bounds of the variables to be
##' optimized over
##' @param upper optional vector containing the upper bounds of the variables
##' to be optimized over
##' @param parinit optional vector containing the initial values for the
##' variables to be optimized over
##' @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 EQI at par.}
##' @author Victor Picheny 
##' 
##' David Ginsbourger 
##' @examples
##' 
##' set.seed(10)
##' 
##' # Set test problem parameters
##' doe.size <- 10
##' dim <- 2
##' test.function <- get("branin2")
##' lower <- rep(0,1,dim)
##' upper <- rep(1,1,dim)
##' noise.var <- 0.2
##' 
##' # Generate DOE and response
##' doe <- as.data.frame(matrix(runif(doe.size*dim),doe.size))
##' y.tilde <- rep(0, 1, doe.size)
##' for (i in 1:doe.size)  {y.tilde[i] <- test.function(doe[i,]) 
##' + sqrt(noise.var)*rnorm(n=1)}
##' y.tilde <- as.numeric(y.tilde)
##' 
##' # Create kriging model
##' model <- km(y~1, design=doe, response=data.frame(y=y.tilde),
##'      covtype="gauss", noise.var=rep(noise.var,1,doe.size), 
##'      lower=rep(.1,dim), upper=rep(1,dim), control=list(trace=FALSE))
##' 
##' # Optimisation using max_EQI
##' res <- max_EQI(model, new.noise.var=noise.var, type = "UK", 
##' lower=c(0,0), upper=c(1,1)) 
##' X.genoud <- res$par
##' 
##' \dontrun{
##' # Compute actual function and criterion on a grid
##' n.grid <- 12 # Change to 21 for a nicer picture
##' x.grid <- y.grid <- seq(0,1,length=n.grid)
##' design.grid <- expand.grid(x.grid, y.grid)
##' names(design.grid) <- c("V1","V2")
##' nt <- nrow(design.grid)
##' crit.grid <- apply(design.grid, 1, EQI, model=model, new.noise.var=noise.var, beta=.9)
##' 
##' # # 2D plots
##' z.grid <- matrix(crit.grid, n.grid, n.grid)
##' tit <- "Green: best point found by optimizer"
##' filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = rainbow,
##' plot.axes = {title(tit);points(model@@X[,1],model@@X[,2],pch=17,col="blue"); 
##' points(X.genoud[1],X.genoud[2],pch=17,col="green");
##' axis(1); axis(2)})
##' }
##' 
##' @export
max_EQI <-function(model, new.noise.var=0, beta=0.9, q.min=NULL, type = "UK", lower, upper, parinit=NULL, control=NULL) {

  ######### Compute q.min if missing #########
  if (is.null(q.min))
  { pred <- predict.km(model, newdata=model@X, type=type, checkNames = FALSE)
    q.min <- min(pred$mean + qnorm(beta)*pred$sd)  }
  
	EQI.envir <- new.env()	
	environment(EQI) <- environment(EQI.grad) <- EQI.envir 
	gr = EQI.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(EQI, 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, new.noise.var=new.noise.var, beta=beta, q.min=q.min, type=type, envir=EQI.envir
		)
                            
    o$par <- t(as.matrix(o$par))
	colnames(o$par) <- colnames(model@X)
	o$value <- as.matrix(o$value)
	colnames(o$value) <- "EQI"  
	return(list(par=o$par, value=o$value)) 
}