# Example code for talk at GECCO 2017: Simulation-based Test Functions for Optimization Algorithms

You can find the corresponding article in the publications section of this website, and the slides in the presentation section.

The R-code below can be used to reproduce the simple example illustrated in the slides and in the article.

```
#
# To run this script, you need to install R, the language for statistical computing.
# https://cran.r-project.org/
# For ease of use, you may consider RStudio IDE, but this is not required.
# https://www.rstudio.com/
# For a tutorial/introduction to R, see e.g.:
# https://cran.r-project.org/doc/manuals/r-release/R-intro.html

# Note: the following line loads the package CEGO
# If it is not installed, please do so first, using:
# install.packages("CEGO")
# of the package CEGO.
require(CEGO)

nsim <- 10 # number of simulations / realizations
seed <- 12121212 # RNG seed
set.seed(seed)
m <- 100 # number of samples to be simulated at

# objective  function:
fun <- function(x){
exp(-20* x) + sin(6*x^2) + x
}
# "vectorize" target
f <- function(x){sapply(x,fun)}

# distance function (for model/kernel)
dF <- function(x,y)(sum((x-y)^2))

# plot parameters
par(mfrow=c(3,1),mar=c(2.3,2.5,0.2,0.2),mgp=c(1.4,0.5,0))

# create test samples for plots
xtest <- as.list(seq(from=-0,by=0.005,to=1))
plot(xtest,f(xtest),type="l",lty=2,xlab="x",
ylab="Estimation",ylim=c(-0.5,1.7))

# evaluation samples (training data)
xb <- c(0.1,0.4,0.54,0.6,0.8,0.99)
yb <- f(xb)

# create support samples for simulation
x <- as.list(sort(c(runif(m-length(xb)),unlist(xb))))

# fit the model
fit <- modelKriging(xb,yb,
control=list(distanceFunction=dF,
algThetaControl=
list(method="NLOPT_GN_DIRECT_L",
funEvals=100),useLambda=F))
fit

# predicted obj. function values
ypred <- predict(fit,as.list(xtest))\$y
lines(unlist(xtest),ypred)
points(unlist(xb),yb,pch=19)

##############################
# create test functions with non conditional sim.
##############################

fun <- createSimulatedTestFunction(x,fit,nsim,F,seed=1)

ynew <- NULL
for(i in 1:nsim)
ynew <- cbind(ynew,fun[[i]](xtest))

rangeY <- range(ynew)

plot(unlist(xtest),ynew[,1],type="l",
xlab="x",ylab="Non-conditional",ylim=c(-0.5,1.7))
for(i in 2:nsim){
lines(unlist(xtest),ynew[,i],col=i,type="l")
}

##############################
# create test functions with conditional sim.
##############################

fun <- createSimulatedTestFunction(x,fit,nsim,T,seed=1)

ynew <- NULL
for(i in 1:nsim)
ynew <- cbind(ynew,fun[[i]](xtest))

rangeY <- range(ynew)

plot(unlist(xtest),ynew[,1],type="l",
xlab="x",ylab="Conditional",ylim=c(-0.5,1.7))
for(i in 2:nsim){
lines(unlist(xtest),ynew[,i],col=i,type="l")
}
points(unlist(xb),yb,pch=19)

```

# Updated – on CRAN: CEGO Version 2.2.0

A new version (2.2.0) of my R-package CEGO was just uploaded to CRAN.

This update contained some various fixes to code and documentation. Also, the interfacing of C code was reworked, following recent changes to CRAN checks (registering of entry points to compiled code). This may have yielded a speed up for some of the distance calculations. The update also contains code that has been employed in the article Simulation-based Test Functions for Optimization Algorithms (see Publications for the pre-print PDF). This mostly includes functions for Kriging simulation as well as for the generation of corresponding test functions.

# PPSN Tutorial 2016: Surrogate Model Optimization

>>>Tutorial PPSN 2016<<<

In the above attached PDF you can find our slides for the tutorial “Meta-Model Assisted (Evolutionary) Optimization”, held at PPSN 2016.

In addition, please find below the code that is used in the tutorial’s example with R:

```
##
## To run this script, you need to install R, the language for statistical computing.
## https://cran.r-project.org/
## For ease of use, you may consider RStudio IDE, but this is not required.
## https://www.rstudio.com/
## For a tutorial/introduction to R, see e.g.:
## https://cran.r-project.org/doc/manuals/r-release/R-intro.html

## Preparation:
##
## You may need to install some of the packages that this example
## for the PPSN 2016 tutorial is using. If they are available,
## the library("...") command used in the following will
## To install the required packages, uncomment the following lines:
# install.packages("SPOT")

## Initialize random number generator seed. Reproducibility.
set.seed(1)

## Main part
## This example should demonstrate the use of surrogate-modeling
## in optimization. To that end, we first need to define an
## optimization problem to be solved by such methods.
## Clearly, we lack the time to consider a real-world, expensive
## optimization problem. Hence, use the following simple,
## one-dimensional test function, from the book
## A. I. J. Forrester, A. Sobester, A. J. Keane;
## "Engineering Design via Surrogate Modeling";
## Wiley (2008)

objectFun <- function(x){
(6*x-2)^2 * sin(12*x-4)
}

## Plot the function:
par(mar=c(4,4,0.5,4),mgp=c(2,1,0))
curve(objectFun(x),0,1)

## Now, let us assume objectFun is expensive.
## design of experiment, which in this case
## is simply a regular grid:
x <- seq(from=0, by=0.3,to=1)

## Evaluate with objective:
y <- sapply(x,objectFun)

points(x,y)

## Build a model (here: Kriging, with the SPOT package.
## But plenty of alternatives available)
fit <- forrBuilder(as.matrix(x),as.matrix(y),
control=list(uselambda=FALSE #do not use nugget effect (regularization)
))

## Evaluate prediction based on model fit
xtest <- seq(from=0, by=0.001,to=1)
pred <- predict(fit,as.matrix(xtest),predictAll=T)
ypred <- pred\$f
spred <- pred\$s

## Plot the prediction of the model:
lines(xtest,ypred,lty=2)

## Plot suggested candidate solution
points(xtest[which.min(ypred)],ypred[which.min(ypred)],col="black",pch=20)

## Calculate expected improvement (EI)
ei <- 10^(-spotInfillExpImp(ypred,spred,min(y)))
par(new = T)
plot(xtest,ei,lty=3, type="l", axes=F, xlab=NA, ylab=NA,
ylim=rev(range(ei)))
axis(side = 4); mtext(side = 4, line = 2, 'EI')
## but note: EI is on a different scale

## Plot suggested candidate solution, based on EI
points(xtest[which.max(ei)],ei[which.max(ei)],col="red",pch=20)
newx <- xtest[which.max(ei)]

x <- c(x,newx)
y <- c(y,objectFun(newx))

## Now repeat the same as often as necessary:
repeatThis <- expression({
curve(objectFun(x),0,1)
points(x,y)
fit <- forrBuilder(as.matrix(x),as.matrix(y),
control=list(uselambda=FALSE
))
xtest <- seq(from=0, by=0.001,to=1)
pred <- predict(fit,as.matrix(xtest),predictAll=T)
ypred <- pred\$f
spred <- pred\$s
lines(xtest,ypred,lty=2)
points(xtest[which.min(ypred)],ypred[which.min(ypred)],col="black",pch=20)
ei <- 10^(-spotInfillExpImp(ypred,spred,min(y)))
par(new = T)
plot(xtest,ei,lty=3, type="l", axes=F, xlab=NA, ylab=NA,
ylim=rev(range(ei)))
axis(side = 4); mtext(side = 4, line = 2, 'EI')
points(xtest[which.max(ei)],ei[which.max(ei)],col="red",pch=20)
newx <- xtest[which.max(ei)]
x <- c(x,newx)
y <- c(y,objectFun(newx)) })
eval(repeatThis)
eval(repeatThis)
eval(repeatThis)
eval(repeatThis)
eval(repeatThis)
eval(repeatThis)
eval(repeatThis)
eval(repeatThis)
## Observation:
## EI looks noisy, strange.
## Predicted mean has low accuracy.
## Why?
## If repeated to often -> Numerical issues:
## Due to close spacing of candidates -> Problem for Kriging model
## Potential remedy: use regularization with nugget and reinterpolation.
## Note: Other interpretation of such an issue may be convergence of
## the optimization process. But this is not necessarily correct.

## repeat as often as necessary (but now with regularization):
repeatThis <- expression({
curve(objectFun(x),0,1)
points(x,y)
fit <- forrBuilder(as.matrix(x),as.matrix(y),
control=list(
uselambda=TRUE, # Use nugget (parameter lambda)
reinterpolate=T # Reinterpolation, to fix uncertainty estimates, etc.
))
xtest <- seq(from=0, by=0.001,to=1)
pred <- predict(fit,as.matrix(xtest),predictAll=T)
ypred <- pred\$f
spred <- pred\$s
lines(xtest,ypred,lty=2)
points(xtest[which.min(ypred)],ypred[which.min(ypred)],col="black",pch=20)
ei <- 10^(-spotInfillExpImp(ypred,spred,min(y)))
par(new = T)
plot(xtest,ei,lty=3, type="l", axes=F, xlab=NA, ylab=NA,
ylim=rev(range(ei)))
axis(side = 4); mtext(side = 4, line = 2, 'EI')
points(xtest[which.max(ei)],ei[which.max(ei)],col="red",pch=20)
newx <- xtest[which.max(ei)]
x <- c(x,newx)
y <- c(y,objectFun(newx))
})
eval(repeatThis)
eval(repeatThis)

```

# Updated – on CRAN: CEGO Version 2.1.0

The new version 2.1.0 of my R-package CEGO was just uploaded to CRAN.

This update mostly contained methods that were used in the upcoming paper Efficient Global Optimization with Indefinite Kernels (see Publications for the pre-print PDF). That includes methods to check whether matrices are (conditionally) definite, methods to repair definiteness and other conditions and their integration into a Kriging model to enable modeling and optimization with indefinite kernels.