R/sampleFromEI.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

#' Samples \code{n} points from a distribution proportional to the expected
#' improvement (EI) computed from a \code{km} object.
#' 
#' @title Sampling points according to the expected improvement criterion
#' @param model an object of class \code{\link[DiceKriging]{km}},
#' @param minimization logical specifying if EI is used in minimiziation or in
#' maximization,
#' @param n number of points to be sampled,
#' @param initdistrib matrix of candidate points. 
#' @param lower vector of lower bounds,
#' @param upper vector of upper bounds,
#' @param T optional scalar : if provided, it replaces the current minimum (or
#' maximum) of observations.
#' @return A \code{n*d} matrix containing the sampled points. If \code{NULL}, \code{1000*d} points 
#' are obtained by latin hypercube sampling,
#' @author Sebastien Marmin 
#' 
#' Clement Chevalier 
#' 
#' David Ginsbourger 
#' @seealso \code{\link{EI}}, \code{\link[DiceKriging]{km}}, \code{\link{qEI}}
#' @references
#' 
#' 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.
#' @keywords optimization
#' @examples
#' 
#' 
#' \donttest{
#' 
#' set.seed(004)
#' 
#' # a 9-points factorial design, and the corresponding responses
#' 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) 
#' lower <- c(0,0)
#' upper <- c(1,1)
#' names(response.branin) <- "y" 
#' 
#' 
#' # model identification
#' fitted.model <- km(~1, design=design.fact, response=response.branin, 
#'                    covtype="gauss", control=list(pop.size=50,trace=FALSE), parinit=c(0.5, 0.5))
#' 
#' # sample a 30 point batch
#' batchSize <- 30
#' x <- sampleFromEI(model = fitted.model, n = batchSize, lower = lower, upper = upper)
#' 
#' # graphics 
#' # displays the EI criterion, the design points in black and the EI-sampled points in red.
#' nGrid <- 15
#' gridAxe1 <- seq(lower[1],upper[1],length=nGrid)
#' gridAxe2 <- seq(lower[2],upper[2],length=nGrid)
#' grid <- expand.grid(gridAxe1,gridAxe2)
#' aa <- apply(grid,1,EI,model=fitted.model)
#' myMat <- matrix(aa,nrow=nGrid)
#' image(x = gridAxe1, y = gridAxe2, z = myMat, 
#'       col = colorRampPalette(c("darkgray","white"))(5*10), 
#'       ylab = names(design.fact)[1], xlab=names(design.fact)[2], 
#'       main = "Sampling from the expected improvement criterion", 
#'       axes = TRUE, zlim = c(min(myMat), max(myMat)))
#' contour(x = gridAxe1, y = gridAxe2, z = myMat, 
#'         add = TRUE, nlevels = 10)
#' points(x[,1],x[,2],pch=19,col='red')
#' points(fitted.model@@X[,1],fitted.model@@X[,2],pch=19)
#' }
#' 
#' @export sampleFromEI
sampleFromEI <- function(model,minimization=TRUE,n= 1,initdistrib=NULL,lower=rep(0,model@d),upper=rep(1,model@d),T=NULL){
	
  #Generic function to sample from a density proportional to EI
  
  #### Arguments:
  # model : a km object
  # minimization : do we perform a global minimization or maximization ?
  # n : the total size of the sample distributed from a density proportional to EI
  # initdistrib : how these points are generated (default: "LHS" sequence)
  # d : dimension of the input set
  # lower, upper : arrays of size d with the lower (resp. upper) bounds in each dimension
  # T: the threshold used to compute the EI
  
  if(is.null(model)){
    message("Error in samplefromEI. A km object needs to be given")
    return(NULL)
  }
  
  d <- model@d
  
  if (length(lower) != length(upper) ){
    message("Error in samplefromEI : lower and upper must have the same length")
    return(NULL)
  }
  
  if(is.null(initdistrib)) {
	initial.points <- t(lower+t(randLHS(n=d*1000,k=d))*(upper-lower))
	} else {
	initial.points <- initdistrib
	}
	if(d==1) initial.points <- matrix(initial.points,ncol=1)
	#prediction on these initial candidate points
	predictions <- predict(object=model,newdata=initial.points,type="UK",checkNames=FALSE,cov.compute = FALSE,se.compute = TRUE,light.return = TRUE)
	  
  mn <- predictions$mean
  sn <- predictions$sd
	
  if(is.null(T)){
    if(minimization) T <- min(model@y)
    if(!minimization) T <- max(model@y)
  }
  
  uu <- (mn-T)/sn
  if(minimization) uu <- -1*uu
  ei <- sn*( uu*pnorm(uu) + dnorm(uu) )
  
  if ( minimization)  ei[is.nan(ei)] <- (T-mn)
  if (!minimization)  ei[is.nan(ei)] <- (mn-T)  
  ei <- ei*(ei>0)
  if(sum(ei!=0)<n) {
    maxEI <- max(ei)
    if (maxEI > 0) {
      ei <- ei/maxEI+1/(1000*d)
    } else {
      ei <- rep(1,(1000*d))
    }
  }
  my.indices <- sample(1:(1000*d), n, replace = FALSE, prob = ei)
  my.points <- initial.points[my.indices,]
		
	if(d==1) my.points <- matrix(my.points,ncol=1)
	if(n==1) my.points <- matrix(my.points,ncol=d)
    
	return(my.points)
}

ranperm <- function(X, N) order(runif(N))
randLHS <- function(n,k) {
  P <- matrix(nrow = n, ncol = k)
  P <- apply(P, 2, ranperm, N = n)
  P <- P - 1 + matrix(runif(n * k), nrow = n, ncol = k)
  return(P/n)
}