Sequential Parameter Optimization Examples: How to Call SPOT

Most simple example: Kriging + LHS + predicted mean optimization (not expected improvement)

The most simple spot() run requires the following parameters:
- fun: objective function, e.g., funSphere
- lower: vector of lower limits of the search space, one entry for each dimension. The length of the lower vector determines the problem dimension.
- upper: vector of upper limits of the search space, one entry for each dimension
A default number of 20 function evalutions is used (funEvals = 20).

res <- spot(,funSphere,
            c(-1,-1),
            c(1,1))
cbind(res$xbest, res$ybest)
#>             [,1]         [,2]         [,3]
#> [1,] -0.01032073 -0.001220482 0.0001080069

With y

Explicit starting point x=0.05, specified as a matrix
1-dim objective function funSphere: dimension is specified by the length of the lowervector
20 function evaluations
surrogate model: buildKriging, which is the default
acquisition function is based on expected improvement
optimizer on the surrogate: optimLBFGSB

set.seed(1)
res <- spot(
  x = matrix(c(0.05), 1, 1),
  fun = funSphere,
  lower = -1,
  upper = 1,
  control = list(
    funEvals = 20,
    model = buildKriging,
    optimizer = optimDE,
    modelControl = list(target = "y")
  )
)
cbind(res$xbest, res$ybest)
#>               [,1]         [,2]
#> [1,] -3.711794e-05 1.377741e-09

With expected improvement

Explicit starting point x=0.05, specified as a matrix
1-dim objective function funSphere: dimension is specified by the length of the lowervector
20 function evaluations
surrogate model: buildKriging, which is the default
acquisition function is based on expected improvement
optimizer on the surrogate: optimLBFGSB

set.seed(1)
res <- spot(
  x = matrix(c(0.05), 1, 1),
  fun = funSphere,
  lower = -1,
  upper = 1,
  control = list(
    funEvals = 20,
    model = buildKriging,
    optimizer = optimDE,
    modelControl = list(target = "ei")
  )
)
cbind(res$xbest, res$ybest)
#>               [,1]         [,2]
#> [1,] -8.664214e-06 7.506861e-11

Plot results (basic plot only)

plot(res$ybestVec, log = "y", type="b")

Multi starts on the surrogate

res <- spot(
   x = matrix(c(0.05), 1, 1),
  fun = funSphere,
  lower = -1,
  upper = 1,
  control = list(
    funEvals = 20,
    multiStart = 3,
    model = buildKriging,
    optimizer = optimDE,
    modelControl = list(target = "y"),
    designControl = list(size=10)
  )
)
cbind(res$xbest, res$ybest)
#>               [,1]         [,2]
#> [1,] -1.188579e-07 1.412721e-14

plot(res$ybestVec, log = "y", type="b")

Bayesian Optimization

res <- spot(
  x = NULL,
  fun = funSphere,
  lower = -1,
  upper = 1,
  control = list(
    funEvals = 20,
    multiStart = 2,
    model = buildBO,
    optimizer = optimDE,
    modelControl = list(target = "y"),
    designControl = list(size=10)
  )
)
cbind(res$xbest, res$ybest)

Left: best y val, right: search on the surrogate

par(mfrow=c(1,2))
plot(res$ybestVec, log = "y", type="b")
plot(res$ySurr[!is.na(res$ySurr)], type="b")

par(mfrow=c(1,1))

Replicate (reevaluate the best solution)

Results can be replicated (which is helpful for noisy objective functions)
Example below is based on the previous example, which uses a deterministic function, so all replicated values will be identical.
The result to be reevaluted is specified as the vector lower, upper has to be set to this value as well (which restricts the search space to one point only).
The number of replications is specified via funEvals.

lower <-res$xbest
resBestReplicated <- spot(x=NULL,
  fun = funSphere,
  lower = lower,
  upper = lower,
  control = list(
    replicateResult = TRUE,
    funEvals = 10
  )
)
cbind(resBestReplicated$x,
      resBestReplicated$y)
#>                [,1]         [,2]
#>  [1,] -1.188579e-07 1.412721e-14
#>  [2,] -1.188579e-07 1.412721e-14
#>  [3,] -1.188579e-07 1.412721e-14
#>  [4,] -1.188579e-07 1.412721e-14
#>  [5,] -1.188579e-07 1.412721e-14
#>  [6,] -1.188579e-07 1.412721e-14
#>  [7,] -1.188579e-07 1.412721e-14
#>  [8,] -1.188579e-07 1.412721e-14
#>  [9,] -1.188579e-07 1.412721e-14
#> [10,] -1.188579e-07 1.412721e-14

Modify the initial design size

Initial design size is size = 3
Altogether funEvals = 4 function evaluations are use: the surrogate model is build only once

res <- spot(,
            fun = funSphere,
            lower = c(-1,-1),
            upper = c(1,1),
            control = list( #funEvals=4,
                           #verbosity=0
                           #,designControl = list(size = 3)
            )
            )
cbind(res$x, res$y)
#>              [,1]         [,2]         [,3]
#>  [1,] -0.82703849  0.629569141 1.0803499676
#>  [2,]  0.43543891  0.331775522 0.2996820387
#>  [3,] -0.10133625 -0.162986007 0.0368334749
#>  [4,] -0.71405733 -0.209124373 0.5536108690
#>  [5,] -0.48714723  0.179569698 0.2695577016
#>  [6,]  0.13123246  0.588739411 0.3638360533
#>  [7,]  0.39571081 -0.455261850 0.3638503990
#>  [8,]  0.84643223  0.874071413 1.4804483548
#>  [9,] -0.35183768 -0.843796492 0.8357822734
#> [10,]  0.75936722 -0.797770098 1.2130756996
#> [11,]  0.09865883 -0.097796930 0.0192978047
#> [12,]  0.02070748 -0.103273036 0.0110941196
#> [13,] -0.01032073 -0.001220482 0.0001080069
#> [14,] -0.01867512  0.046198481 0.0024830599
#> [15,]  0.09245131  0.029345694 0.0094084142
#> [16,] -0.11676163  0.070734437 0.0186366381
#> [17,]  0.15592757 -0.015600590 0.0245567870
#> [18,] -0.05542781  0.010803192 0.0031889511
#> [19,] -0.02151719 -0.054124859 0.0033924901
#> [20,]  0.15854184 -0.069685662 0.0299916064

With additional start points:

Three starting points x= (0.05, 0.1), x=(-0.1, 0.01) and x=(0.01, 0.02), which are specified as a matrix. Note the byrow argument.
One additional starting point is used designControl = list(size = 1)
The surrogate model is build once, because the budget is set to funEvals=5

res <- spot(x = matrix(c(0.05,0.1,
                         -0.1, 0.01,
                         0.01, 0.02),3,2, byrow = TRUE),
            fun = funSphere,
            lower = c(-1,-1),
            upper = c(1,1),
             control = list(funEvals=5,
            designControl = list(size = 1)
            )
            )
cbind(res$x, res$y)
#>            [,1]       [,2]      [,3]
#> [1,]  0.0500000 0.10000000 0.0125000
#> [2,] -0.1000000 0.01000000 0.0101000
#> [3,]  0.0100000 0.02000000 0.0005000
#> [4,] -0.2557522 0.81641558 0.7319436
#> [5,]  0.8709325 0.06266123 0.7624498

comparison with optim (simulated annealing):

fs <- function(x){
  x[1]^2 + x[2]^2}
res <- optim(par=c(0.05, 0,1), 
             fn=fs, 
             method = "SANN",
             control = list(maxit = 20))
cbind(t(res$par),res$value)
#>      [,1] [,2] [,3]   [,4]
#> [1,] 0.05    0    1 0.0025

Surrogate-based optimization with optimLBFGSB instead of LHS

res <- spot(x = matrix(c(0.05,0.1),1,2),
            fun = funSphere,
            lower = c(-2,-3),
            upper = c(1,2),
            control=list(funEvals=20,
                         optimizer=optimLBFGSB,
                          modelControl=list(target="ei")
                         )
            )
cbind(res$xbest, res$ybest)
#>      [,1] [,2]   [,3]
#> [1,] 0.05  0.1 0.0125

Multi start might be useful for LBFGSB:

res <- spot(x = matrix(c(0.05,0.1),1,2),
            fun = funSphere,
            lower = c(-2,-3),
            upper = c(1,2),
            control=list(funEvals=20,
                         optimizer=optimLBFGSB,
                          modelControl=list(target="ei"),
                         multiStart = 2
                         )
            )
cbind(res$xbest, res$ybest)
#>             [,1]          [,2]         [,3]
#> [1,] -0.01265257 -5.745797e-05 0.0001600907

Random Forest instead of Kriging

res <- spot(,
            funSphere,
            c(-1,-1),
            c(1,1),
            control=list(model=buildRandomForest))
cbind(res$xbest, res$ybest)
#>              [,1]         [,2]         [,3]
#> [1,] 0.0004720909 -0.006394821 4.111661e-05

LM instead of Kriging

res <- spot(,
            funSphere,
            c(-2,-3),
            c(1,2),
     control=list(model=buildLM)) #lm as surrogate
cbind(res$xbest, res$ybest)
#>           [,1]      [,2]      [,3]
#> [1,] 0.1531584 0.3294388 0.1319874

Bayesian Optimization

res <- spot(
  x = matrix(c(0.05, 0.1), 1, 2),
  fun = funSphere,
  lower = c(-1, -1),
  upper = c(1, 1),
  control = list(
    funEvals = 20,
    model = buildBO,
    optimizer = optimLBFGSB,
    modelControl = list(target = "ei")
  )
)
cbind(res$xbest, res$ybest)
#>      [,1] [,2]   [,3]
#> [1,] 0.05  0.1 0.0125

LM and local optimizer (which for this simple example is perfect)

res <- spot(,funSphere,c(-2,-3),c(1,2),
   control=list(model=buildLM, optimizer=optimLBFGSB))
res$xbest
#>           [,1]      [,2]
#> [1,] 0.1531584 0.3294388

Lasso and local optimizer NLOPTR


res <- spot(,funSphere,
            lower = c(-2,-3),
            upper = c(1,2),
             control = list(funEvals=50,
                            model=buildLasso, 
                optimizer = optimNLOPTR,
            designControl = list(size = 20)
            ))
res$xbest
#>            [,1]       [,2]
#> [1,] 0.09877579 -0.1655962

Kriging and local optimizer LBFGSB

res <- spot(,funSphere,c(-2,-3),c(1,2), 
   control=list(model=buildKriging, optimizer = optimLBFGSB))
cbind(res$xbest, res$ybest)
#>             [,1]         [,2]         [,3]
#> [1,] -0.01773916 -0.001142165 0.0003159822

Kriging and local optimizer NLOPTR

res <- spot(,funSphere,c(-2,-3),c(1,2), 
     control=list(model=buildKriging, optimizer = optimNLOPTR))
cbind(res$xbest, res$ybest)
#>             [,1]        [,2]         [,3]
#> [1,] -0.02263374 0.002972108 0.0005211198

Or a different Kriging model:

res <- spot(,funSphere,c(-2,-3),c(1,2),
 control=list(model=buildKrigingDACE, optimizer=optimLBFGSB))
res$xbest
#>            [,1]        [,2]
#> [1,] 0.09403992 -0.04007123

Transform x values

x values are transformed via 2^x
output shows natural and transformed x values

# Use transformed input values
set.seed(1)
f2 <- function(x){2^x}
lower <- c(-100, -100)
upper <- c(10, 10)
transformFun <- rep("f2", length(lower))
res <- spot(x=NULL,
            fun=funSphere,
            lower=lower, 
            upper=upper,
             control=list(funEvals=20,
                          modelControl=list(target="ei"),
                          optimizer=optimLBFGSB,
                          transformFun=transformFun,
                          verbosity = 0,
                          progress = FALSE,
                          plots = FALSE))
print(cbind(res$x, res$xt, res$y))
#>              [,1]        [,2]         [,3]         [,4]         [,5]
#>  [1,]  -90.487117  -10.373697 5.763198e-28 7.537129e-04 5.680832e-07
#>  [2,]  -21.050860  -26.752346 4.603198e-07 8.845885e-09 2.119726e-13
#>  [3,]  -50.573494  -53.964230 5.968447e-16 5.690468e-17 3.594617e-31
#>  [4,]  -84.273153  -56.501840 4.278122e-26 9.800567e-18 9.605111e-35
#>  [5,]  -71.793098  -35.123667 2.444129e-22 2.671301e-11 7.135848e-22
#>  [6,]  -37.782215  -12.619332 4.230777e-12 1.589287e-04 2.525834e-08
#>  [7,]  -23.235905  -70.039402 1.012268e-07 8.242125e-22 1.024687e-14
#>  [8,]    1.553773    3.073928 2.935839e+00 8.420627e+00 7.952611e+01
#>  [9,]  -64.351072  -91.408807 4.250078e-20 3.042336e-28 1.806317e-39
#> [10,]   -3.234803  -88.877355 1.062251e-01 1.758936e-27 1.128378e-02
#> [11,] -100.000000  -18.335894 7.888609e-31 3.022359e-06 9.134652e-12
#> [12,] -100.000000  -45.470842 7.888609e-31 2.050749e-14 4.205573e-28
#> [13,] -100.000000  -78.918523 7.888609e-31 1.750481e-24 3.064184e-48
#> [14,] -100.000000  -14.180231 7.888609e-31 5.386729e-05 2.901685e-09
#> [15,]  -54.483350   10.000000 3.970793e-17 1.024000e+03 1.048576e+06
#> [16,]   -9.233011   10.000000 1.661830e-03 1.024000e+03 1.048576e+06
#> [17,]   10.000000  -40.224817 1.024000e+03 7.782579e-13 1.048576e+06
#> [18,]  -91.586630 -100.000000 2.689534e-28 7.888609e-31 7.233657e-56
#> [19,]   10.000000   -2.281238 1.024000e+03 2.057212e-01 1.048576e+06
#> [20,]  -31.756840 -100.000000 2.755743e-10 7.888609e-31 7.594120e-20
print(which.min(res$y[, 1, drop = FALSE]))
#> [1] 18
cbind(res$xbest, res$ybest)
#>           [,1] [,2]         [,3]
#> [1,] -91.58663 -100 7.233657e-56

With noise: (this takes some time)

Note: If control$OCBA then control$replicates and control$designControl$replicates should be larger than one. The parameter control$OCBABudget defines how many additional candidate solutions (from xnew) should be evaluated. So, the total number of new evaluations is the sum of control$designControl$replicates and control$OCBABudget.

# noisy objective
fNoise = function(x){funSphere(x) + rnorm(nrow(x))}
res1 <-
  spot(x = NULL, 
       fun = fNoise, 
       lower = c(-2, -3), 
       upper = c(1, 2),
    control = list(funEvals = 40, noise = TRUE, verbosity=0))
# noise with replicated evaluations
res2 <-
  spot(x = NULL, 
       fun = fNoise, 
       lower = c(-2, -3), 
       upper = c(1, 2),
    control = list(
      verbosity=0,
      funEvals = 40,
      noise = TRUE,
      replicates = 2,
      designControl = list(replicates = 2)
    ))
# and with OCBA
res3 <- spot(x = NULL, 
       fun = fNoise, 
       lower = c(-2, -3), 
       upper = c(1, 2),
  control = list(
    verbosity=0,
    funEvals = 40,
    noise = TRUE,
    replicates = 2,
    OCBA = TRUE,
    OCBABudget = 5,
    designControl = list(replicates = 2)
  )
)
# Check results with non-noisy function:
funSphere(res1$xbest)
#>           [,1]
#> [1,] 0.7806337
funSphere(res2$xbest)
#>           [,1]
#> [1,] 0.3128347
funSphere(res3$xbest)
#>           [,1]
#> [1,] 0.1319874

lower <-res3$xbest
resBestReplicated <- spot(
  ,
  fun = fNoise,
  lower = lower,
  upper = lower,
  control = list(
    funEvals = 10,
    replicateResults = TRUE
   )
)
cbind(resBestReplicated$x,
      resBestReplicated$y)
#>            [,1]      [,2]        [,3]
#>  [1,] 0.1531584 0.3294388 -0.82121574
#>  [2,] 0.1531584 0.3294388  0.53853014
#>  [3,] 0.1531584 0.3294388  2.36124961
#>  [4,] 0.1531584 0.3294388 -1.38250960
#>  [5,] 0.1531584 0.3294388  0.07027999
#>  [6,] 0.1531584 0.3294388 -0.01528338
#>  [7,] 0.1531584 0.3294388  1.67358048
#>  [8,] 0.1531584 0.3294388 -0.84986826
#>  [9,] 0.1531584 0.3294388  0.62856558
#> [10,] 0.1531584 0.3294388  1.82893529

Random number seed handling

The following is for demonstration only, to be used for random number seed handling in case of external noisy target functions.

res1a <- spot(,function(x,seed){set.seed(seed);funSphere(x)+rnorm(nrow(x))},
     c(-2,-3),c(1,2),control=list(funEvals=25,noise=TRUE,seedFun=1))

res1b <- spot(,function(x,seed){set.seed(seed);funSphere(x)+rnorm(nrow(x))},
     c(-2,-3),c(1,2),control=list(funEvals=25,noise=TRUE,seedFun=1))

res2 <- spot(,function(x,seed){set.seed(seed);funSphere(x)+rnorm(nrow(x))},
     c(-2,-3),c(1,2),control=list(funEvals=25,noise=TRUE,seedFun=2))

sprintf("Should be equal: %f = %f. Should be different:  %f", res1a$ybest, res1b$ybest, res2$ybest)
#> [1] "Should be equal: -0.574866 = -0.574866. Should be different:  -0.807379"

Handling factor variables

Note: factors should be coded as integer values, i.e., 1,2,…,n First, we create a test function with a factor variable:

braninFunctionFactor <- function (x) {
   y <- (x[2]  - 5.1/(4 * pi^2) * (x[1] ^2) + 5/pi * x[1]  - 6)^2 +
     10 * (1 - 1/(8 * pi)) * cos(x[1] ) + 10
   if(x[3]==1)
     y <- y +1
   else if(x[3]==2)
     y <- y -1
   return(y)
}

Vectorize the test function.

objFun <- function(x){apply(x,1,braninFunctionFactor)}

Run spot.

set.seed(1)
res <- spot(fun=objFun,lower=c(-5,0,1),upper=c(10,15,3),
            control=list(model=buildKriging,
                         types= c("numeric","numeric","factor"),
                         optimizer=optimLHD))
 res$xbest
#>          [,1]     [,2] [,3]
#> [1,] 3.116086 2.166734    3
 res$ybest
#>           [,1]
#> [1,] 0.4174566

High dimensional problem

n <- 10
a <- rep(0,n)
b <- rep(1,n)

First, we consider the default spot setting with buildKriging().

tic <- proc.time()[3]
res0 <- spot(x=NULL, funSphere, lower = a, upper = b, 
             control=list(funEvals=30))
toc <- proc.time()[3]
sprintf("value: %f, time: %f",  res0$ybest, toc-tic)
#> [1] "value: 1.026498, time: 3.622000"

Then, we use the buildGaussianProcess() model.

tic <- proc.time()[3]
res1 <-  spot(x=NULL, funSphere, lower = a, upper = b, 
             control=list(funEvals=30, 
                          model = buildGaussianProcess))
toc <- proc.time()[3]
sprintf("value: %f, time: %f",  res1$ybest, toc-tic)
#> [1] "value: 0.819676, time: 0.408000"

Run SPOT with logging

## run spot without log
res <- spot(fun = funSphere,
            lower=c(0,0),
            upper=c(100,100)
)
## run spot with log (3-dim "y" output: first is y, last are x val (here 2-dim))
funSphereLog <- function(x){
  cbind(funSphere(x),x)
}
res2 <- spot(fun = funSphereLog,
            lower=c(0,0),
            upper=c(100,100)
)
res$logInfo
#> [1] NA
res2$logInfo
#>             [,1]       [,2]
#>  [1,]  8.6480755 81.4784571
#>  [2,] 71.7719454 66.5887761
#>  [3,] 44.9331873 41.8506996
#>  [4,] 14.2971337 39.5437814
#>  [5,] 25.6426384 58.9784849
#>  [6,] 56.5616232 79.4369705
#>  [7,] 69.7855406 27.2369075
#>  [8,] 92.3216115 93.7035707
#>  [9,] 32.4081160  7.8101754
#> [10,] 87.9683608 10.1114951
#> [11,]  0.8695938  3.4789743
#> [12,]  6.7788538  7.4840977
#> [13,] 10.5474376  0.9580708
#> [14,]  0.4337100  9.4376420
#> [15,]  1.1522512 14.3301670
#> [16,]  2.2464471  6.1725581
#> [17,]  5.8038975  2.6715535
#> [18,]  4.9227178  8.2127192
#> [19,]  8.8485311 11.1649501
#> [20,]  3.2337531  2.8721535
## re-evaluation of the x-values and comparison with evaluated x-values:
funSphere(res2$logInfo) == res$y
#>       [,1]
#>  [1,] TRUE
#>  [2,] TRUE
#>  [3,] TRUE
#>  [4,] TRUE
#>  [5,] TRUE
#>  [6,] TRUE
#>  [7,] TRUE
#>  [8,] TRUE
#>  [9,] TRUE
#> [10,] TRUE
#> [11,] TRUE
#> [12,] TRUE
#> [13,] TRUE
#> [14,] TRUE
#> [15,] TRUE
#> [16,] TRUE
#> [17,] TRUE
#> [18,] TRUE
#> [19,] TRUE
#> [20,] TRUE

Hybrid optimization

res <- spot(fun = funSphere, lower = c(-5,-5),
                upper = c(5,5), 
                control = list(funEvals = 20,
                directOpt = optimNLOPTR,
                directOptControl = list(funEvals = 10)
                ))
str(res)
#> List of 14
#>  $ xbest    : num [1, 1:2] 0 0
#>  $ ybest    : num 0
#>  $ xBestOcba: logi NA
#>  $ yBestOcba: logi NA
#>  $ x        : num [1:33, 1:2] -4.135 2.177 -0.507 -3.57 -2.436 ...
#>  $ xt       : logi NA
#>  $ y        : num [1:33, 1] 27.009 7.492 0.921 13.84 6.739 ...
#>  $ logInfo  : logi NA
#>  $ count    : int 20
#>  $ msg      : chr "budget exhausted"
#>  $ modelFit :List of 33
#>   ..$ thetaLower      : num 1e-04
#>   ..$ thetaUpper      : num 100
#>   ..$ types           : chr [1:2] "numeric" "numeric"
#>   ..$ algTheta        :function (x = NULL, fun, lower, upper, control = list(), ...)  
#>   ..$ budgetAlgTheta  : num 200
#>   ..$ optimizeP       : logi FALSE
#>   ..$ useLambda       : logi TRUE
#>   ..$ lambdaLower     : num -6
#>   ..$ lambdaUpper     : num 0
#>   ..$ startTheta      : NULL
#>   ..$ reinterpolate   : logi TRUE
#>   ..$ target          : chr "y"
#>   ..$ modelInitialized: logi TRUE
#>   ..$ x               : num [1:19, 1:2] -4.135 2.177 -0.507 -3.57 -2.436 ...
#>   ..$ y               : num [1:19, 1] 27.009 7.492 0.921 13.84 6.739 ...
#>   ..$ normalizeymin   : num 0
#>   ..$ normalizeymax   : num 1
#>   ..$ scaledx         : num [1:19, 1:2] 0 0.7544 0.4337 0.0675 0.2031 ...
#>   ..$ normalizexmin   : num [1:2] -4.14 -4.22
#>   ..$ normalizexmax   : num [1:2] 4.23 4.37
#>   ..$ Theta           : num [1:2] -0.0511 -0.0104
#>   ..$ dmodeltheta     : num [1:2] 0.889 0.976
#>   ..$ Lambda          : num -6
#>   ..$ dmodellambda    : num 1e-06
#>   ..$ yonemu          : num [1:19, 1] -55 -74.5 -81.1 -68.1 -75.2 ...
#>   ..$ ssq             : num 749
#>   ..$ mu              : num 82
#>   ..$ Psi             : num [1:19, 1:19] 1 0.586 0.687 0.789 0.902 ...
#>   ..$ Psinv           : num [1:19, 1:19] 268 344 -376 340 -1781 ...
#>   ..$ nevals          : num 1200
#>   ..$ like            : num [1, 1] -14.9
#>   ..$ returnCrossCor  : logi FALSE
#>   ..$ min             : num 0.0027
#>   ..- attr(*, "class")= chr "kriging"
#>  $ ybestVec : num [1:20] 0.921 0.921 0.921 0.921 0.921 ...
#>  $ ySurr    : num [1:20] NA NA NA NA NA NA NA NA NA NA ...
#>  $ control  :List of 36
#>   ..$ funEvals             : num 20
#>   ..$ lower                : num [1:2] -5 -5
#>   ..$ upper                : num [1:2] 5 5
#>   ..$ design               :function (x = NULL, lower, upper, control = list())  
#>   ..$ designControl        :List of 1
#>   .. ..$ types: chr [1:2] "numeric" "numeric"
#>   ..$ directOpt            :function (x = NULL, fun, lower, upper, control = list(), ...)  
#>   ..$ directOptControl     :List of 1
#>   .. ..$ funEvals: num 10
#>   ..$ infillCriterion      : NULL
#>   ..$ maxTime              : num Inf
#>   ..$ model                :function (x, y, control = list())  
#>   ..$ modelControl         :List of 2
#>   .. ..$ modelInitialized: logi FALSE
#>   .. ..$ types           : chr [1:2] "numeric" "numeric"
#>   ..$ multiStart           : num 0
#>   ..$ noise                : logi FALSE
#>   ..$ OCBA                 : logi FALSE
#>   ..$ OCBABudget           : num 1
#>   ..$ optimizer            :function (x = NULL, fun, lower, upper, control = list(), ...)  
#>   ..$ optimizerControl     :List of 1
#>   .. ..$ types: chr [1:2] "numeric" "numeric"
#>   ..$ parNames             : chr [1:2] "x1" "x2"
#>   ..$ plots                : logi FALSE
#>   ..$ progress             : logi FALSE
#>   ..$ replicates           : num 1
#>   ..$ replicateResult      : logi FALSE
#>   ..$ returnFullControlList: logi TRUE
#>   ..$ seedFun              : logi NA
#>   ..$ seedSPOT             : num 1
#>   ..$ subsetSelect         :function (x, y, control)  
#>   ..$ subsetControl        : list()
#>   ..$ time                 :List of 3
#>   .. ..$ maxTime  : num Inf
#>   .. ..$ startTime: POSIXct[1:1], format: "2022-06-11 15:04:11"
#>   .. ..$ endTime  : POSIXct[1:1], format: "2022-06-11 15:04:12"
#>   ..$ tolerance            : num 1.49e-08
#>   ..$ transformFun         : logi(0) 
#>   ..$ types                : chr [1:2] "numeric" "numeric"
#>   ..$ verbosity            : num 0
#>   ..$ xNewActualSize       : num 0
#>   ..$ xBestOcba            : logi NA
#>   ..$ yBestOcba            : logi NA
#>   ..$ yImputation          :List of 3
#>   .. ..$ imputeCriteriaFuns:List of 3
#>   .. .. ..$ :function (x)  
#>   .. .. ..$ :function (x)  
#>   .. .. ..$ :function (x)  
#>   .. ..$ handleNAsMethod   : NULL
#>   .. ..$ penaltyImputation : num 3

Handling constraints

Definitions and constraints

library(babsim.hospital)
n <- 29 
reps <- 2
funEvals <- 3*n 
size <- 2*n
x0 <- matrix(as.numeric(babsim.hospital::getParaSet(5374)[1,-1]),1,)
bounds <- getBounds()
a <- bounds$lower
b <- bounds$upper
g <- function(x) {
      return(rbind(a[1] - x[1], x[1] - b[1], a[2] - x[2], x[2] - b[2], 
                   a[3] - x[3], x[3] - b[3], a[4] - x[4], x[4] - b[4], 
                   a[5] - x[5], x[5] - b[5], a[6] - x[6], x[6] - b[6], 
                   a[7] - x[7], x[7] - b[7], a[8] - x[8], x[8] - b[8], 
                   a[9] - x[9], x[9] - b[9], a[10] - x[10], x[10] - b[10],
                   a[11] - x[11], x[11] - b[11], a[12] - x[12],  x[12] - b[12],
                   a[13] - x[13], x[13] - b[13], a[14] - x[14],  x[14] - b[14],
                   a[15] - x[15], x[15] - b[15], a[16] - x[16],  x[16] - b[16],
                   a[17] - x[17], x[17] - b[17], a[18] - x[18],  x[18] - b[18],
                   a[19] - x[19], x[19] - b[19], a[20] - x[20],  x[20] - b[20],
                   a[21] - x[21], x[21] - b[21], a[22] - x[22],  x[22] - b[22],
                   a[23] - x[23], x[23] - b[23], a[24] - x[24],  x[24] - b[24],
                   a[25] - x[25], x[25] - b[25], a[26] - x[26],  x[26] - b[26],
                   a[27] - x[27], x[27] - b[27], x[15] + x[16] - 1, 
                   x[17] + x[18] + x[19] - 1, x[20] + x[21] - 1, x[23] + x[29] - 1)
      )
  }

res <- spot(
  x = x0,
  fun = funBaBSimHospital,
  lower = a,
  upper = b,
  verbosity = 0,
  control = list(
    funEvals = 2 * funEvals,
    noise = TRUE,
    designControl = list(# inequalityConstraint = g,
      size = size,
      retries = 1000),
    optimizer = optimNLOPTR,
    optimizerControl = list(
      opts = list(algorithm = "NLOPT_GN_ISRES"),
      eval_g_ineq = g
    ),
    model =  buildKriging,
    plots = FALSE,
    progress = TRUE,
    directOpt = optimNLOPTR,
    directOptControl = list(funEvals = 0),
    eval_g_ineq = g
  )
)
print(res)

GECCO Industrial Challenge 2021

A description of the challenge can be found here: GECCO Industrial Challenge 2021. In short the goal of the challenge is to find an optimal parameter configuration for the BabSim.Hospital simulator. This is a noisy and complex real-world problem.

Evaluation Using the Docker Container

In order to be able to execute the necessary code of the GECCO Industrial challenge 2021 you will need to have Docker installed in your machine. On your terminal console an evaluation of the BabSim.Hospital should looks like the command below. This command will automatically download the Docker image with the BabSim.Hospital code in it (may need sudo rights to download). Take care, the formatting of the symbols - and ’ can cause this command not to work on your terminal:

# docker run --rm mrebolle/r-geccoc:Track1 -c 'Rscript objfun.R "6,7,3,3,3,5,3,3,25,17,2,1,0.25,0.05,0.07,0.005,0.07,1e-04,0.08,0.25,0.08,0.5,1e-06,2,1e-06,1e-06,1,2,0.5"'

An optimization run with SPOT, using the Docker command call as objective function, can be directly implemented in R as follows:

library(SPOT)

evalFun <- function(candidateSolution){
    evalCommand <- paste0("docker run --rm mrebolle/r-geccoc:Track1 -c ", "'","Rscript objfun.R ")
    parsedCandidate <- paste(candidateSolution, sep=",", collapse = ",")
    return(as.numeric(system(paste0(evalCommand, '"', parsedCandidate, '"', "'"), intern = TRUE)))
}

#The BabSim.Hospital requires 29 parameters. Here we specify the upper and lower bounds
lower <- c(6,7,3,3,3,5,3,3,25,17,2,1,0.25,0.05,0.07,
           0.005,0.07,1e-04,0.08,0.25,0.08,0.5,1e-06,
           2,1e-06,1e-06,1,2,0.5)

upper<- c(14,13,7,9,7,9,5,7,35,25,5,7,2,0.15,0.11,0.02,
          0.13,0.002,0.12,0.35,0.12,0.9,0.01,4,1.1,0.0625,
          2,5,0.75)

wFun <- wrapFunction(evalFun)

n <- 29 
reps <- 2
funEvals <- 10*n 
size <- 2*n
x0<-matrix(lower,nrow = 1)

res <- spot(x = x0,
  fun = wFun,
  lower = lower,
  upper = upper,
  control = list(
    funEvals = 2 * funEvals,
    noise = TRUE,
    designControl = list(
      size = size,
      retries = 1000),
    optimizer = optimNLOPTR,
    optimizerControl = list(
      opts = list(algorithm = "NLOPT_GN_ISRES")
    ),
    model =  buildKriging,
    plots = TRUE,
    progress = TRUE,
    directOpt = optimNLOPTR,
    directOptControl = list(funEvals = 0)
  )
)

Evaluation Using the babsim.hospital R Package

The optimization of the BabSim.Hospital parameters can also be executed directly using the babsim.hospital package.

The babsim.hospital package can be installed by downloading the source from the Gitlab repository and building the package.

git clone http://owos.gm.fh-koeln.de:8055/bartz/babsim.hospital.git

library(SPOT)
library(babsim.hospital)

n <- 29 
reps <- 2
funEvals <- 3*n 
size <- 2*n
#Get suggested parameter values as initial point in the optimization run
x0 <- matrix(as.numeric(babsim.hospital::getParaSet(5374)[1,-1]),1,)
bounds <- getBounds()
a <- bounds$lower
b <- bounds$upper
g <- function(x) {
      return(rbind(a[1] - x[1], x[1] - b[1], a[2] - x[2], x[2] - b[2], 
                   a[3] - x[3], x[3] - b[3], a[4] - x[4], x[4] - b[4], 
                   a[5] - x[5], x[5] - b[5], a[6] - x[6], x[6] - b[6], 
                   a[7] - x[7], x[7] - b[7], a[8] - x[8], x[8] - b[8], 
                   a[9] - x[9], x[9] - b[9], a[10] - x[10], x[10] - b[10],
                   a[11] - x[11], x[11] - b[11], a[12] - x[12],  x[12] - b[12],
                   a[13] - x[13], x[13] - b[13], a[14] - x[14],  x[14] - b[14],
                   a[15] - x[15], x[15] - b[15], a[16] - x[16],  x[16] - b[16],
                   a[17] - x[17], x[17] - b[17], a[18] - x[18],  x[18] - b[18],
                   a[19] - x[19], x[19] - b[19], a[20] - x[20],  x[20] - b[20],
                   a[21] - x[21], x[21] - b[21], a[22] - x[22],  x[22] - b[22],
                   a[23] - x[23], x[23] - b[23], a[24] - x[24],  x[24] - b[24],
                   a[25] - x[25], x[25] - b[25], a[26] - x[26],  x[26] - b[26],
                   a[27] - x[27], x[27] - b[27], x[15] + x[16] - 1, 
                   x[17] + x[18] + x[19] - 1, x[20] + x[21] - 1, x[23] + x[29] - 1)
      )
  }

wrappedFunBab <- function(x){
  print(SPOT::funBaBSimHospital(x, region = 5374, nCores = 1))
}
res <- spot(
  x = x0,
  fun = wrappedFunBab,
  lower = a,
  upper = b,
  control = list(
    funEvals = 2 * funEvals,
    noise = TRUE,
    designControl = list(
      size = size,
      retries = 1000),
    optimizer = optimNLOPTR,
    optimizerControl = list(
      opts = list(algorithm = "NLOPT_GN_ISRES"),
      eval_g_ineq = g
    ),
    model =  buildKriging,
    plots = FALSE,
    progress = TRUE,
    directOpt = optimNLOPTR,
    directOptControl = list(funEvals = 0),
    eval_g_ineq = g
  )
)
print(res)

Benchmarking

res <- spot(
  x = NULL,
  fun = funGoldsteinPrice,
  lower = rep(0,2),
  upper = rep(1,2),
  control = list(
    funEvals = 20,
    multiStart = 5,
    model = buildKriging,
    optimizer = optimDE,
    modelControl = list(target = "ei"),
    designControl = list(size=10)
  )
)
cbind(res$xbest, res$ybest)

### Compare Performance
# rm(list=ls())
library(microbenchmark)
library(SPOT)
set.seed(1)
n <- 3
low = -2
up = 2
a = runif(n, low, 0)
b = runif(n, 0, up)
x0 = a + runif(n)*(b-a)
plot(a, type = "l", ylim=c(up,low))
lines(b)
lines(x0)
x0 = matrix( x0, nrow = 1)
  
set.seed(1)
  perf1 <- spot(x= x0, funSphere, a, b, control=list(time=list(maxTime = 0.25), funEvals=10*n, plots=FALSE,
                                                     model = buildKriging, optimizer=optimNLOPTR))
  set.seed(1)
  perf2 <- spot(x= x0, funSphere, a, b, control=list(time=list(maxTime = 0.25), funEvals=10*n, plots=FALSE,
                model = buildGaussianProcess, optimizer=optimNLOPTR, directOptControl = list(funEvals=0)))
  
  set.seed(1)
  perf3 <- spot(x= x0, funSphere, a, b, control=list(time=list(maxTime = 0.25), funEvals=10*n, plots=FALSE,
                                                     model = buildGaussianProcess, optimizer=optimNLOPTR,
                                                     directOptControl = list(funEvals=10)))

### Plot Repeats (Sphere Function)
#  rm(list=ls())
library(microbenchmark)
library(SPOT)
set.seed(1)
n <- 2
low = -2
up = 2
a = runif(n, low, 0)
b = runif(n, 0, up)
x0 = a + runif(n)*(b-a)
 #plot(a, type = "l", ylim=c(up,low))
#  #lines(b)
#  #lines(x0)
x0 = matrix( x0, nrow = 1)
#  
reps <- 10
end <- 10*n
ninit <- n
#  
  progSpot <- matrix(NA, nrow = reps, ncol = end)
  for(r in 1:reps){
    set.seed(r)
    x0 <- a + runif(n)*(b-a)
    x0 = matrix( x0, nrow = 1)
    sol <- spot(x= x0, funGoldsteinPrice, a, b, control=list(funEvals=end,
                                                     model = buildGaussianProcess,
                                                     optimizer=optimNLOPTR,
                                                     directOptControl = list(funEvals=0),
                                                     designControl = list(size = ninit)))
    progSpot[r, ] <- prepareBestObjectiveVal(sol$y, end)
  
}
#  
matplot(t(progSpot), type="l", col="gray", lty=1,
          xlab="n: blackbox evaluations", ylab="best objective value", log="y")
  abline(v=ninit, lty=2)
  legend("topright", "seed LHS", lty=2, bty="n")

Simple Benchmark

Setup

### Plot Repeats (Sphere Function)
#  rm(list=ls())
library(SPOT)
#
set.seed(1)
n <- 2 #dim
repeats <- 30 # repeats
#
end <- 50 * n
ninit <- 5* n
# objective fun
fun <- funGoldsteinPrice
# 
low = rep(-1, n)
up = rep(1,n)
x0 <- c()
for(i in 1:repeats){
x0 = rbind(x0, low + runif(n) * (up - low))
}

Optim L-BFGS-B

f <- fun
  fprime <- function(x) {
    x <- matrix( x, 1)
    ynew <- as.vector(f(x))
    y <<- c(y, ynew)
    return(ynew)
}
#  
progOptim <- matrix(NA, nrow=nrow(x0), ncol=end)
for (r in 1:nrow(x0)) {
    y <- c()
    os <- optim(x0[r, ,drop=FALSE], fprime, lower=low, upper=up, method="L-BFGS-B")
    progOptim[r,] <- prepareBestObjectiveVal(y, end)
  }
#

Bayesian Optimization with neg log EI

progSpotBOEI <- matrix(NA, nrow = nrow(x0), ncol = end)
for (r in 1:nrow(x0)) {
  set.seed(r)
 sol <- spot(
    x = x0[r, ,drop=FALSE],
    fun=fun,
    a,
    b,
    control = list(
      funEvals = end,
      multiStart = 2,
      model = buildBO,
      optimizer = optimDE,
      modelControl = list(target = "negLog10ei"),
      directOptControl = list(funEvals =
                                0),
      designControl = list(size = ninit)
    )
  )
  progSpotBOEI[r,] <- prepareBestObjectiveVal(sol$y, end)
}
#

Bayesian Optimization with Y

progSpotBOY <- matrix(NA, nrow = nrow(x0), ncol = end)
for (r in 1:nrow(x0)) {
  set.seed(r)
 sol <- spot(
    x = x0[r, ,drop=FALSE],
    fun=fun,
    lower = low,
    upper = up,
    control = list(
      funEvals = end,
      multiStart = 2,
      model = buildBO,
      optimizer = optimDE,
      modelControl = list(target = "y"),
      directOptControl = list(funEvals =
                                0),
      designControl = list(size = ninit)
    )
  )
  progSpotBOY[r,] <- prepareBestObjectiveVal(sol$y, end)
}
#

SPOT Kriging with EI

progSpotBuildKrigingEi <- matrix(NA, nrow = nrow(x0), ncol = end)
for (r in 1:nrow(x0)) {
  set.seed(r)
 sol <- spot(
    x = x0[r, ,drop=FALSE],
    fun=fun,
   lower = low,
    upper = up,
    control = list(
      funEvals = end,
      multiStart = 2,
      model = buildKriging,
      optimizer = optimDE,
      modelControl = list(target = "ei"),
      directOptControl = list(funEvals =
                                0),
      designControl = list(size = ninit)
    )
  )
  progSpotBuildKrigingEi[r,] <- prepareBestObjectiveVal(sol$y, end)
}
#

SPOT Kriging with Y

progSpotBuildKrigingY <- matrix(NA, nrow = nrow(x0), ncol = end)
for (r in 1:nrow(x0)) {
  set.seed(r)
 sol <- spot(
    x = x0[r, ,drop=FALSE],
    fun=fun,
   lower = low,
    upper = up,
    control = list(
      funEvals = end,
      multiStart = 2,
      model = buildKriging,
      optimizer = optimDE,
      modelControl = list(target = "y"),
      directOptControl = list(funEvals =
                                0),
      designControl = list(size = ninit)
    )
  )
  progSpotBuildKrigingY[r,] <- prepareBestObjectiveVal(sol$y, end)
}
#

Summary

print("optim:")
summary(progOptim[,end])
print("SPOT BO EI:")
summary(progSpotBOEI[,end])
print("SPOT BO Y:")
summary(progSpotBOY[,end])
print("SPOT Kriging EI:")
summary(progSpotBuildKrigingEi[,end])
print("SPOT Kriging Y:")
summary(progSpotBuildKrigingY[,end])

Plot Run Curves

matplot(
  colMeans(progOptim),
  type = "l",
  col = 1,
  lty = 1,
  xlab = "n: blackbox evaluations",
  ylab = "avg best objective value",
  ylim = c(-3,3)
)
abline(v = ninit, lty = 2)
lines(colMeans(progSpotBOEI, na.rm=TRUE), col = 2, lwd=2)
lines(colMeans(progSpotBOY, na.rm=TRUE), col = 3, lwd=2)
lines(colMeans(progSpotBuildKrigingEi, na.rm=TRUE), col = 4, lwd=2)
lines(colMeans(progSpotBuildKrigingY, na.rm=TRUE), col = 5, lwd=2)
legend("topright", c("optim", "BO EI", "BO Y", "Krig EI", "Krig Y", "initial design size"), col=c(1:5,1), lwd = c(rep(1,5), 1), lty = c(rep(1,5),2), bty="n")

Plot Boxplots

boxplot(progOptim[,end],
        progSpotBOEI[,end],
        progSpotBOY[,end],
        progSpotBuildKrigingEi[,end],
        progSpotBuildKrigingY[,end],
        names = c("optim", "BO EI", "BO Y", "Krig EI", "Krig Y"),
        xlab="algorithm",
        ylab = "best y value"
        )

Save Results

resBench01 <- list(progOptim=progOptim, progSpotBOEI=progSpotBOEI, progSpotBOY=progSpotBOY, progSpotBuildKrigingEi=progSpotBuildKrigingEi, progSpotBuildKrigingY=progSpotBuildKrigingY)
usethis::use_data(resBench01)

Handle NAs

dim = 2
lower = c(-2, -3)
upper = c(1, 2)

control <- spotControl(dimension = dim)
control$verbosity <- 0
control$designControl$size <- 10
control$funEvals <- 15
control$yImputation$handleNAsMethod <- handleNAsMean
res <- spot(x = NULL,
           fun = funError,
           lower = lower,
           upper = upper,
           control)
#> NAs were found and treated!
#> y before treatment:
#>           [,1]
#>  [1,] 4.182852
#>  [2,]      Inf
#>  [3,] 1.248602
#>  [4,] 3.514453
#>  [5,]       NA
#>  [6,] 1.036390
#>  [7,]       NA
#>  [8,] 3.432185
#>  [9,] 7.865728
#> [10,] 6.630543
#> y after treatment:
#>            [,1]
#>  [1,]  4.182852
#>  [2,] 11.616856
#>  [3,]  1.248602
#>  [4,]  3.514453
#>  [5,] 11.616855
#>  [6,]  1.036390
#>  [7,] 11.616856
#>  [8,]  3.432185
#>  [9,]  7.865728
#> [10,]  6.630543
#> NAs were found and treated!
#> y before treatment:
#>             [,1]
#>  [1,]  4.1828515
#>  [2,] 11.6168555
#>  [3,]  1.2486025
#>  [4,]  3.5144534
#>  [5,] 11.6168555
#>  [6,]  1.0363903
#>  [7,] 11.6168555
#>  [8,]  3.4321853
#>  [9,]  7.8657279
#> [10,]  6.6305433
#> [11,]  2.2613854
#> [12,]  0.2015209
#> [13,]  4.7535966
#> [14,]        Inf
#> y after treatment:
#>             [,1]
#>  [1,]  4.1828515
#>  [2,] 11.6168555
#>  [3,]  1.2486025
#>  [4,]  3.5144534
#>  [5,] 11.6168555
#>  [6,]  1.0363903
#>  [7,] 11.6168555
#>  [8,]  3.4321853
#>  [9,]  7.8657279
#> [10,]  6.6305433
#> [11,]  2.2613854
#> [12,]  0.2015209
#> [13,]  4.7535966
#> [14,] 17.8131733
#> NAs were found and treated!
#> y before treatment:
#>             [,1]
#>  [1,]  4.1828515
#>  [2,] 11.6168555
#>  [3,]  1.2486025
#>  [4,]  3.5144534
#>  [5,] 11.6168555
#>  [6,]  1.0363903
#>  [7,] 11.6168555
#>  [8,]  3.4321853
#>  [9,]  7.8657279
#> [10,]  6.6305433
#> [11,]  2.2613854
#> [12,]  0.2015209
#> [13,]  4.7535966
#> [14,] 17.8131733
#> [15,]         NA
#> y after treatment:
#>             [,1]
#>  [1,]  4.1828515
#>  [2,] 11.6168555
#>  [3,]  1.2486025
#>  [4,]  3.5144534
#>  [5,] 11.6168555
#>  [6,]  1.0363903
#>  [7,] 11.6168555
#>  [8,]  3.4321853
#>  [9,]  7.8657279
#> [10,]  6.6305433
#> [11,]  2.2613854
#> [12,]  0.2015209
#> [13,]  4.7535966
#> [14,] 17.8131733
#> [15,] 21.8256805

SPOTVignetteNutshell

Introduction

Sequential Parameter Optimization Examples: How to Call SPOT

Most simple example: Kriging + LHS + predicted mean optimization (not expected improvement)

With y

With expected improvement

Plot results (basic plot only)

Multi starts on the surrogate

Bayesian Optimization

Replicate (reevaluate the best solution)

Modify the initial design size

With additional start points:

Surrogate-based optimization with optimLBFGSB instead of LHS

Random Forest instead of Kriging

LM instead of Kriging

Bayesian Optimization

LM and local optimizer (which for this simple example is perfect)

Lasso and local optimizer NLOPTR

Kriging and local optimizer LBFGSB

Kriging and local optimizer NLOPTR

Or a different Kriging model:

Transform x values

With noise: (this takes some time)

Random number seed handling

Handling factor variables

High dimensional problem

Run SPOT with logging

Hybrid optimization

Handling constraints

GECCO Industrial Challenge 2021

Evaluation Using the Docker Container

Evaluation Using the babsim.hospital R Package

Benchmarking

Simple Benchmark

Setup

Optim L-BFGS-B

Bayesian Optimization with neg log EI

Bayesian Optimization with Y

SPOT Kriging with EI

SPOT Kriging with Y

Summary

Plot Run Curves

Plot Boxplots

Save Results

Handle NAs