if (FALSE) {
#---------------------------------------------------------------------------
# EHI surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------
set.seed(25468)
library(DiceDesign)
d <- 2
n.obj <- 2
fname <- "P1"
n.grid <- 51
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
nappr <- 15
design.grid <- round(maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design, 1)
response.grid <- t(apply(design.grid, 1, fname))
paretoFront <- t(nondominated_points(t(response.grid)))
mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])
model <- list(mf1, mf2)
EHI_grid <- apply(test.grid, 1, crit_EHI, model = list(mf1, mf2),
critcontrol = list(refPoint = c(300, 0)))
lower <- rep(0, d)
upper <- rep(1, d)
omEGO <- crit_optimizer(crit = "EHI", model = model, lower = lower, upper = upper,
optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2),
critcontrol = list(refPoint = c(300, 0)))
print(omEGO)
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
matrix(EHI_grid, nrow = n.grid), main = "Expected Hypervolume Improvement",
xlab = expression(x[1]), ylab = expression(x[2]), color = terrain.colors,
plot.axes = {axis(1); axis(2);
points(design.grid[, 1], design.grid[, 2], pch = 21, bg = "white");
points(omEGO$par, col = "red", pch = 4)
}
)
#---------------------------------------------------------------------------
# SMS surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------
SMS_grid <- apply(test.grid, 1, crit_SMS, model = model,
critcontrol = list(refPoint = c(300, 0)))
lower <- rep(0, d)
upper <- rep(1, d)
omEGO2 <- crit_optimizer(crit = "SMS", model = model, lower = lower, upper = upper,
optimcontrol = list(method="genoud", pop.size = 200, BFGSburnin = 2),
critcontrol = list(refPoint = c(300, 0)))
print(omEGO2)
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
matrix(pmax(0,SMS_grid), nrow = n.grid), main = "SMS Criterion (>0)",
xlab = expression(x[1]), ylab = expression(x[2]), color = terrain.colors,
plot.axes = {axis(1); axis(2);
points(design.grid[, 1], design.grid[, 2], pch = 21, bg = "white");
points(omEGO2$par, col = "red", pch = 4)
}
)
#---------------------------------------------------------------------------
# Maximin Improvement surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------
EMI_grid <- apply(test.grid, 1, crit_EMI, model = model,
critcontrol = list(nb_samp = 20, type ="EMI"))
lower <- rep(0, d)
upper <- rep(1, d)
omEGO3 <- crit_optimizer(crit = "EMI", model = model, lower = lower, upper = upper,
optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2))
print(omEGO3)
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
matrix(EMI_grid, nrow = n.grid), main = "Expected Maximin Improvement",
xlab = expression(x[1]), ylab = expression(x[2]), color = terrain.colors,
plot.axes = {axis(1);axis(2);
points(design.grid[, 1], design.grid[, 2], pch = 21, bg = "white");
points(omEGO3$par, col = "red", pch = 4)
}
)
#---------------------------------------------------------------------------
# crit_SUR surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------
library(KrigInv)
integration.param <- integration_design_optim(lower = c(0, 0), upper = c(1, 1), model = model)
integration.points <- as.matrix(integration.param$integration.points)
integration.weights <- integration.param$integration.weights
precalc.data <- list()
mn.X <- sn.X <- matrix(0, n.obj, nrow(integration.points))
for (i in 1:n.obj){
p.tst.all <- predict(model[[i]], newdata = integration.points, type = "UK",
checkNames = FALSE)
mn.X[i,] <- p.tst.all$mean
sn.X[i,] <- p.tst.all$sd
precalc.data[[i]] <- precomputeUpdateData(model[[i]], integration.points)
}
critcontrol <- list(integration.points = integration.points,
integration.weights = integration.weights,
mn.X = mn.X, sn.X = sn.X, precalc.data = precalc.data)
EEV_grid <- apply(test.grid, 1, crit_SUR, model=model, paretoFront = paretoFront,
critcontrol = critcontrol)
lower <- rep(0, d)
upper <- rep(1, d)
omEGO4 <- crit_optimizer(crit = "SUR", model = model, lower = lower, upper = upper,
optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2))
print(omEGO4)
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid),
matrix(pmax(0,EEV_grid), n.grid), main = "EEV criterion", nlevels = 50,
xlab = expression(x[1]), ylab = expression(x[2]), color = terrain.colors,
plot.axes = {axis(1); axis(2);
points(design.grid[,1], design.grid[,2], pch = 21, bg = "white")
points(omEGO4$par, col = "red", pch = 4)
}
)
# example using user defined optimizer, here L-BFGS-B from base optim
userOptim <- function(par, fn, lower, upper, control, ...){
return(optim(par = par, fn = fn, method = "L-BFGS-B", lower = lower, upper = upper,
control = control, ...))
}
omEGO4bis <- crit_optimizer(crit = "SUR", model = model, lower = lower, upper = upper,
optimcontrol = list(method = "userOptim"))
print(omEGO4bis)
#---------------------------------------------------------------------------
# crit_SMS surface with problem "P1" with 15 design points, using cheapfn
#---------------------------------------------------------------------------
# Optimization with fastfun: SMS with discrete search
# Separation of the problem P1 in two objectives:
# the first one to be kriged, the second one with fastobj
# Definition of the fastfun
f2 <- function(x){
return(P1(x)[2])
}
SMS_grid_cheap <- apply(test.grid, 1, crit_SMS, model = list(mf1, fastfun(f2, design.grid)),
paretoFront = paretoFront, critcontrol = list(refPoint = c(300, 0)))
optimcontrol <- list(method = "pso")
model2 <- list(mf1)
omEGO5 <- crit_optimizer(crit = "SMS", model = model2, lower = lower, upper = upper,
cheapfn = f2, critcontrol = list(refPoint = c(300, 0)),
optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2))
print(omEGO5)
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid),
matrix(pmax(0, SMS_grid_cheap), nrow = n.grid), nlevels = 50,
main = "SMS criterion with cheap 2nd objective (>0)", xlab = expression(x[1]),
ylab = expression(x[2]), color = terrain.colors,
plot.axes = {axis(1); axis(2);
points(design.grid[,1], design.grid[,2], pch = 21, bg = "white")
points(omEGO5$par, col = "red", pch = 4)
}
)
}
Run the code above in your browser using DataLab