#' @title Efficient Global Robust Inversion
#'
#' @description Sequential design of experiments for robust inversion.
#'
#' @param T Target threshold.
#' @param model The current kriging model. km object.
#' @param method Criterion used for choosing observations.
#' Currently, only the \code{"surnew"} criterion is available.
#' @param fun Objective function.
#' @param iter Number of iterations. At each iteration, a batch of batchsize point is evaluated in parallel.
#' @param batchsize The size of the batch of points evaluated at each iteration.
#' @param opt.index Array with integers corresponding to the indices of the nuisance parameters.
#' @param inv.index Array with integers corresponding to the indices of the controlled parameters.
#' @param lower Array of size d. Lower bound of the input domain.
#' @param upper Array of size d. Upper bound of the input domain.
#' @param new.noise.var Optional scalar value of the noise variance of the new observations.
#' @param integcontrol The integcontrol argument used in the \code{integration_design_robinv} function.
#' This applies for the following criterion: \code{"surnew"}. See the help of the
#' \code{integration_design_robinv} function.
#' @param optimcontrol A list used to set up the optimization procedure of the chosen sampling criterion.
#' For the \code{"surnew"} criterion, see the help of the \code{max_surnew_robinv} function.
#' @param kmcontrol Optional list representing the control variables for the re-estimation of the kriging model once new points are sampled.
#' The items are the same as in the \code{km} function of the DiceKriging package.
#' @param ... Other arguments of the objective function fun.
#' @return A list with the following fields.
#' (i) par: The added observations (iter*batchsize) x d matrix,
#' (ii) value: The value of the function fun at the added observations,
#' (iii) lastmodel: The current (last) kriging model of km class,
#' (iv) lastvalue: The value of the criterion at the last added point,
#' (v) allvalues: If an optimization on a discrete set of points is chosen, the value of the criterion at all these points, for the last iteration,
#' (vi) lastintegration.param: For debug. The last value returned by the integration_design_robinv function, if applicable.
#' @export
#' @author Clement Chevalier \email{clement.chevalier@@unine.ch}
#' @examples
#'
#' library(KrigInv)
#' myfun <- function(x) return(-1 * branin_robinv(x))
#' d <- 4
#'
#' set.seed(8)
#'
#' # an example with scaling
#' n0 <- 40
#' T <- -10
#' opt.index <- c(3,4)
#' inv.index <- c(1,2)
#' lower <- rep(0,times=d)
#' upper <- rep(1,times=d)
#'
#' design <- matrix(runif(d*n0),nrow=n0)
#' response <- myfun(design)
#'
#' knots.number <- c(2,2,2,2)
#' knots <- generate_knots(knots.number = knots.number , d = d)
#'
#' model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2",scaling = TRUE,knots = knots)
#'
#' integcontrol <- list(distrib = "surnew",n.points = 100,finaldistrib="surnew",n.candidates=300,
#' nsimu=100,n.optpoints = 50,choose_optpoints=TRUE,n.optpoints.candidates=200)
#'
#' optimcontrol <- list(method = "genoud", pop.size = 400, max.generations = 4, wait.generation = 1)
#' \dontrun{
#' obj <- EGRI(T = T,method="surnew",model = model,fun = myfun,iter = 2,
#' batchsize = 2,opt.index = opt.index,inv.index = inv.index,
#' integcontrol=integcontrol,optimcontrol=optimcontrol)
#' }
EGRI <- function(T, model, method=NULL, fun, iter=1, batchsize=1,
opt.index,inv.index,lower=NULL, upper=NULL,
new.noise.var=NULL,integcontrol=NULL,
optimcontrol=NULL, kmcontrol=NULL,...) {
n <- nrow(model@X)
d <- model@d
if (is.null(kmcontrol$penalty)) kmcontrol$penalty <- model@penalty
if (length(model@penalty==0)) kmcontrol$penalty <- NULL
if (is.null(kmcontrol$optim.method)) kmcontrol$optim.method <- model@optim.method
if (is.null(kmcontrol$parinit)) kmcontrol$parinit <- model@parinit
if (is.null(kmcontrol$control)) kmcontrol$control <- model@control
if (is.null(kmcontrol$cov.reestim)) kmcontrol$cov.reestim <- model@param.estim
if (is.null(kmcontrol$trend.reestim)) kmcontrol$trend.reestim <- model@param.estim
if (is.null(kmcontrol$nugget.reestim)) kmcontrol$nugget.reestim <- model@param.estim
if (is.null(optimcontrol$optim.option)) optimcontrol$optim.option <- 2
if (is.null(method)) method <- "surnew"
if (is.null(lower)) lower <- rep(0,times=d)
if (is.null(upper)) upper <- rep(1,times=d)
integration.param <- NULL
# so far, we force the criterion to be surnew :
method <- "surnew"
for (i in 1:iter) {
if(method== "surnew"){
if(i==1){
if(is.null(integcontrol)) integcontrol$distrib <- "surnew"
if(!is.null(integcontrol)){
integcontrol$finaldistrib <- integcontrol$distrib
integcontrol$distrib <- "surnew"
}
}
tmp_model <- model
result <- matrix(c(0),nrow=batchsize,ncol=d)
for(the_iter in 1:batchsize){
integration.param <- integration_design_robinv(integcontrol=integcontrol,d=d,lower=lower,upper=upper,
opt.index=opt.index,inv.index=inv.index,
model=tmp_model,T=T,seed=8)
oEGOI <- max_surnew_robinv(lower=lower, upper=upper, optimcontrol=optimcontrol, opt.index=opt.index,inv.index=inv.index,
integration.param=integration.param,T=T, model=tmp_model, new.noise.var=new.noise.var)
Xnext <- oEGOI$par
result[the_iter,] <- as.numeric(Xnext)
# need to set a lie
if(the_iter < batchsize){
apred <- predict(object = model,newdata = Xnext,type = "UK",checkNames = FALSE,light.return = TRUE,se.compute = TRUE)
the_lie <- apred$mean+3*apred$sd
tmp_model <- update(object = tmp_model,newX = Xnext,newy = the_lie,newX.alreadyExist = FALSE,cov.reestim = FALSE,trend.reestim = FALSE,nugget.reestim = FALSE,newnoise.var = new.noise.var,kmcontrol = kmcontrol)
}
}
}
print("New points"); print(result)
X.new <- result; y.new <- rep(0,times=nrow(X.new))
for (i in 1:nrow(X.new)) y.new[i] <- fun(X.new[i,],...)
model <- update(object=model,newX=X.new,newy=y.new,cov.reestim=kmcontrol$cov.reestim,trend.reestim=kmcontrol$trend.reestim,nugget.reestim=kmcontrol$nugget.reestim, newnoise.var=rep(new.noise.var,times=batchsize),kmcontrol=kmcontrol)
}
return(list(
par=model@X[(n+1):(n+iter*batchsize),, drop=FALSE],
value=model@y[(n+1):(n+iter*batchsize),, drop=FALSE],
lastmodel=model,
lastvalue=oEGOI$value,
allvalues=oEGOI$allvalues,
lastintegration.param=integration.param
)
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.