R/fastEGO.nsteps.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

##' Sequential EI maximization and model re-estimation, with a number of
##' iterations fixed in advance by the user
##' 
##' Executes \emph{nsteps} iterations of the EGO method to an object of class
##' \code{\link[DiceKriging]{km}}.  At each step, a kriging model is
##' re-estimated (including covariance parameters re-estimation) based on the
##' initial design points plus the points visited during all previous
##' iterations; then a new point is obtained by maximizing the Expected
##' Improvement criterion (\code{\link{EI}}).
##' 
##' 
##' @param model an object of class \code{\link[DiceKriging]{km}} ,
##' @param fun the objective function to be minimized,
##' @param nsteps an integer representing the desired number of iterations,
##' @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 control an optional list of control parameters for EGO. One
##' can control
##' 
##' \code{"warping"} whether or not a warping is applied to the outputs (default FALSE)
##' 
##' \code{"cov.reestim"} whether or not the covariance parameters are estimated at
##' each step (default TRUE)
##' \code{"gpmean.trick"} whether or not EI should be replaced periodically by the GP mean
##'  (default FALSE)
##'  
##' \code{"gpmean.freq"} frequency at which EI is replaced by the GP mean (default 1e4)
##' 
##' \code{"always.sample"} if TRUE, forces observation even if it creates poor conditioning
##' 
##' @param trace between -1 (no trace) and 3 (full messages)
##' @param n.cores number of cores used for EI maximisation
##' @param ... additional parameters to be given to \code{fun}
##' @return A list with components:
##' 
##' \item{par}{a data frame representing the additional points visited during
##' the algorithm,}
##' 
##' \item{value}{a data frame representing the response values at the points
##' given in \code{par},}
##' 
##' \item{npoints}{an integer representing the number of parallel computations
##' (=1 here),}
##' 
##' \item{nsteps}{an integer representing the desired number of iterations
##' (given in argument),}
##' 
##' \item{lastmodel}{an object of class \code{\link[DiceKriging]{km}}
##' corresponding to the last kriging model fitted. If warping is true, 
##' \code{y} values are normalized (warped) and will not match \code{value}.}
##' @author Victor Picheny
##' 
##' @seealso \code{\link{EI}}, \code{\link{max_crit}}, \code{\link{EI.grad}}
##' @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.
##' 
##' J. Mockus (1988), \emph{Bayesian Approach to Global Optimization}. Kluwer
##' academic publishers.
##' 
##' T.J. Santner, B.J. Williams, and W.J. Notz (2003), \emph{The design and
##' analysis of computer experiments}, Springer.
##' 
##' M. Schonlau (1997), \emph{Computer experiments and global optimization},
##' Ph.D. thesis, University of Waterloo.
##' @keywords optimize
##' @examples
##' 
##' set.seed(123)
##' ###############################################################
##' ### 	10 ITERATIONS OF EGO ON THE BRANIN FUNCTION, 	   ####
##' ###	 STARTING FROM 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 n steps
##' nsteps <- 5 
##' lower <- rep(0,d) 
##' upper <- rep(1,d)     
##' oEGO <- fastEGO.nsteps(model=fitted.model1, fun=branin, nsteps=nsteps, lower=lower, upper=upper)
##' print(oEGO$par)
##' print(oEGO$value)
##' 
##' # graphics
##' n.grid <- 15 # Was 20, reduced to 15 for speeding up compilation
##' 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, pch=19, col="red")
##' text(oEGO$par[,1], oEGO$par[,2], labels=1:nsteps, pos=3)
##' 
##' @export fastEGO.nsteps
fastEGO.nsteps <-function(model, fun, nsteps, lower, upper, control=NULL, trace=0, n.cores=1, ...) {
  
  n <- nrow(model@X)
  n.proxy <- n.unif <- 0
  d <- model@d
  added.obs <- c()
  added.X <- c()
  
  if (is.null(control$warping)) control$warping <- FALSE
  if (is.null(control$cov.reestim)) control$cov.reestim <- TRUE
  if (!is.null(control$gpmean.freq)) gpmean.freq <- control$gpmean.freq else gpmean.freq <- 1e4
  if (!is.null(control$gpmean.trick)) gpmean.trick <- control$gpmean.trick else gpmean.trick <- FALSE
  if (!is.null(control$always.sample)) always.sample <- control$always.sample else always.sample <- FALSE
  
  observations <- model@y
  
  if (trace > 1) message("Ite | max(EI) ||  ", colnames(model@X),"|  obj \n")
  
  for (i in 1:nsteps) {
    
    if (((i %% gpmean.freq)==0) && (gpmean.trick)) {
      # Bypass EI maximization to switch to gpmean
      sol <- list(value=0)
    } else {
      sol <- try(max_crit(model=model, lower=lower, upper=upper, control=control, n.cores=n.cores))  
    }
    
    ## Exit if optimization failed
    if (typeof(sol) == "character") {
      warning("Unable to maximize EI at iteration ", i, "- optimization stopped \n Last model returned \n")
      
      par <- values <- c()
      if (i > 1) {
        par <- model@X[(n+1):model@n,, drop = FALSE]
        values <- model@y[(n+1):model@n]
      }
      return(list(par=par, values=values, nsteps = i-1, lastmodel = model))
    }
    
    if (sol$value < max(1e-16, sqrt(model@covariance@sd2/1e16)) && gpmean.trick) {
      # Restart optimization with a proxy criterion
      if (trace>0)
        if ((i %% gpmean.freq)==0) {
          message("Forcing proxy acquisition criterion.\n")
        } else {
          message("No point found with non-null EI, restarted with proxy criterion.\n")
          n.proxy <- n.proxy+1
        }                
      
      sol <- try(max_crit(model=model, lower=lower, upper=upper, control=control, proxy=TRUE, n.cores=n.cores))

      ## Exit if optimization failed
      if (typeof(sol) == "character") {
        warning("Unable to maximize proxy criterion - optimization stopped \n Last model returned \n")
        
        par <- values <- c()
        if (i > 1) {
          par <- model@X[(n+1):model@n,, drop = FALSE]
          values <- model@y[(n+1):model@n]
        }
        return(list(par=par, values=values, nsteps = i-1, lastmodel = model))
      }
    } else {
      n.proxy <- 0
    }
    
    ## Check if optimization does not return already known point
    replace.point = FALSE
    
    if (!always.sample) {
      if (checkPredict(sol$par, model= list(model), threshold = 1e-6)){
        if (trace>0)
          message("Optimization failed (an existing point was selected), so a random point is taken instead.\n")
        sol$par <- matrix(runif(d) * (upper - lower) + lower, nrow = 1) 
        n.unif <- n.unif + 1
      } else {
        n.unif <- 0
      }
    } else {
      if (checkPredict(sol$par, model= list(model), threshold = 1e-6)){
        if (trace>0)
          message("Point too close to an existing one - replacement scheme activated.\n")
        replace.point = TRUE
      }
    }
    
    ## Update
    X.new <- matrix(as.numeric(sol$par), nrow=1, ncol=d)
    Y.new <- try(fun(as.numeric(sol$par),...))

    if (typeof(Y.new) == "character") {
      if (trace>0) {
        warning("Unable to compute objective function at iteration ", i, "- optimization stopped \n")
        warning("Problem occured for the design: ", X.new, "\n Last model returned \n")
      }
      
      par <- values <- c()
      if (i > 1) {
        par <- added.X 
        values <- added.obs
      }
      return(list(par=par, values=values, nsteps = i-1, lastmodel = model))
    }
    if (trace > 1) message( i, "|", signif(sol$val,3), "||", signif(X.new,3), "|", signif(Y.new,3), "\n")
    
    # Transform Ynew if needed
    added.obs <- c(added.obs, Y.new)
    added.X <- rbind(added.X, X.new)
    if (control$warping) {
      Y.new <- warp(Y.new, control$wparam)
    }
    
    # Update models
    observations <- rbind(observations, Y.new)
    newmodel <- model
    
    if (!replace.point) {
      if (n.proxy < 5 && n.unif <3) {
        newmodel <- try(update(object = model, newX = X.new, newy=Y.new, newX.alreadyExist=FALSE,
                               cov.reestim = control$cov.reestim), silent = TRUE)
  
        if (typeof(newmodel) == "character" && control$cov.reestim) {
          newmodel <- try(update(object = model, newX = X.new, newy=Y.new, newX.alreadyExist=FALSE, cov.reestim = FALSE), silent = TRUE)
          if (typeof(newmodel) != "character") {
            if (trace > 0) message("Error in hyperparameter estimation - old hyperparameter values used instead for the objective model \n")
          } else {
            if (length(model@covariance@nugget) == 0) {
              newmodel <- try(km(formula=model@trend.formula, design=rbind(model@X, X.new), response=c(model@y, Y.new), 
                                 covtype=model@covariance@name, coef.trend=model@trend.coef, nugget= model@covariance@sd2/1e12,
                                 multistart=model@control$multistart, 
                                 coef.cov=covparam2vect(model@covariance), coef.var=model@covariance@sd2, iso=is(model@covariance,"covIso")))
              if (typeof(newmodel) != "character") {
                if (trace > 0) message("Error in hyperparameter estimation - nugget added to the model \n")
              }
            }
          }
        }
      } else {
        # If too many proxy or random iterations => update model with smaller range
        newmodel <- try(km(formula=model@trend.formula, design=rbind(model@X, X.new), response=c(model@y, Y.new), 
                           covtype=model@covariance@name, coef.trend=model@trend.coef, nugget= model@covariance@nugget,
                           multistart=model@control$multistart, 
                           coef.cov=covparam2vect(model@covariance)/1.5, coef.var=model@covariance@sd2, iso=is(model@covariance,"covIso")))
  
        if (typeof(newmodel) == "character" && length(model@covariance@nugget) == 0) {
          newmodel <- try(km(formula=model@trend.formula, design=rbind(model@X, X.new), response=c(model@y, Y.new), 
                             covtype=model@covariance@name, coef.trend=model@trend.coef, nugget= model@covariance@sd2/1e12,
                             multistart=model@control$multistart, 
                             coef.cov=covparam2vect(model@covariance)/1.5, coef.var=model@covariance@sd2, iso=is(model@covariance,"covIso")))
          if (typeof(newmodel) != "character") {
            if (trace > 0) message("Too many proxy/random iterations - model updated with smaller covariance ranges and nugget \n")
            model@lower <- model@lower/1.5
            model@upper <- model@upper/1.5
          }
        } else {
          if (trace > 0) message("Too many proxy/random iterations - model updated with smaller covariance ranges \n")
        }
      }
    } else {
      # Replace current point in the model by a new point if better
      if (Y.new < min(model@y)) {
        ibest <- which.min(model@y)
        model@y[ibest] <- Y.new
        model@X[ibest,] <- X.new
        newmodel <- computeAuxVariables(model)
      }
    }
    
    if (typeof(newmodel) == "character") {
      warning("Unable to udpate kriging model at iteration", i-1, "- optimization stopped \n")
      warning("lastmodel is the model at iteration", i-1, "\n")
      warning("par and values contain the ",i, "th observation \n \n")
      if (i > 1) allX.new <- rbind(model@X[(n+1):model@n,, drop=FALSE], X.new)
      return(list(
        par    = added.X,
        values = added.obs,
        nsteps = i, 
        lastmodel = model))
    } else {
      newmodel@known.param <- model@known.param
      newmodel@case <- model@case
      newmodel@param.estim <- model@param.estim
      newmodel@method <- model@method
      newmodel@penalty <- model@penalty
      newmodel@optim.method <- model@optim.method
      newmodel@lower <- model@lower
      newmodel@upper <- model@upper
      newmodel@control <- model@control
      
      model <- newmodel
    }
    if (trace > 2) message("Covariance parameters: ", model@covariance@sd2, "|", covparam2vect(model@covariance), "\n")
  }
  
  if(trace>0) message("\n")
  #### End of main loop ################
  
  return(list(par=added.X, value=added.obs, npoints=1, nsteps=nsteps, lastmodel=model))
}