# SPOTVignetteNutshell

• Package version SPOT should be larger than 2.6.0.
library("SPOT")
packageVersion("SPOT")
#> [1] '2.11.14'

## Introduction

• The performance of modern search heuristics relies crucially on their parameterizations—or, statistically speaking, on their factor settings.
• Finding good parameter settings for an optimization algorithm will be referred to as tuning.
• This vignette is an abbreviated version of a the technical report “In a Nutshell-The Sequential Parameter Optimization Toolbox”.
• It illustrates how spot can be called.
• The extended version of this vignette is available on arXiv and will be updated regularly.

## 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

### 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

• 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

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
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