# SPOTVignetteNutshell

- Package version SPOT should be larger than `2.6.0`.

    library("SPOT")
    packageVersion("SPOT")
    #> [1] '2.8.2'

## 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,
                -1,
                1)
    cbind(res$xbest, res$ybest)
    #>              [,1]         [,2]
    #> [1,] -0.001340476 1.796875e-06

### 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 `lower`vector
- 20 function evaluations
- surrogate model: `buildKriging`, which is the default
- acquisition function is based on expected improvement
- optimizer on the surrogate: `optimLBFGSB`

    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,] -6.191359e-05 3.833293e-09

### Plot results (basic plot only)

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

### Multi starts on the surrogate

    res <- spot(
      x = NULL,
      fun = funSphere,
      lower = -1,
      upper = 1,
      control = list(
        funEvals = 20,
        multiStart = 2,
        model = buildKriging,
        optimizer = optimDE,
        modelControl = list(target = "y"),
        designControl = list(size=10)
      )
    )
    cbind(res$xbest, res$ybest)
    #>              [,1]         [,2]
    #> [1,] 6.747016e-06 4.552222e-11

    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,] 6.747016e-06 4.552222e-11
    #>  [2,] 6.747016e-06 4.552222e-11
    #>  [3,] 6.747016e-06 4.552222e-11
    #>  [4,] 6.747016e-06 4.552222e-11
    #>  [5,] 6.747016e-06 4.552222e-11
    #>  [6,] 6.747016e-06 4.552222e-11
    #>  [7,] 6.747016e-06 4.552222e-11
    #>  [8,] 6.747016e-06 4.552222e-11
    #>  [9,] 6.747016e-06 4.552222e-11
    #> [10,] 6.747016e-06 4.552222e-11

## 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,
                designControl = list(size = 3)
                )
                )
    cbind(res$x, res$y)
    #>         [,1]       [,2]      [,3]
    #>   -0.3945281  0.3745242 0.2959208
    #>   -0.1988787 -0.8626836 0.7837758
    #>    0.9322598 -0.2156288 0.9156041
    #> 1 -0.3945281  0.3745242 0.2959208

### 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]
    #>    0.0500000  0.10000000 0.0125000
    #>   -0.1000000  0.01000000 0.0101000
    #>    0.0100000  0.02000000 0.0005000
    #>   -0.2557522  0.81641558 0.7319436
    #> 1  0.3967616 -0.01536565 0.1576559

- 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.01451384 -0.006339419 0.0002508398

- 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.00384341 -0.02166477 0.0004841339

### 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.09272048 0.07403003 0.01407753

### 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]
    #>  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]
    #>  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]
    #>  0.1531584 0.3294388

### Lasso and local optimizer NLOPTR

    res <- spot(,funSphere,c(-2,-3),c(1,2),
       control=list(model=buildLasso, optimizer = optimNLOPTR))
    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold

    #> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
    #> fold
    res$xbest
    #>       [,1]      [,2]
    #>  0.1531584 0.3294388

### 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.01573496 -0.008860252 0.000326093

### 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.01440329 0.0006858711 0.0002079252

### 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.03586733 -0.004407048

### 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("identity", length(lower))
    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]
    #>    -90.487117  -10.373697 5.763198e-28 7.537129e-04 5.680832e-07
    #>    -21.050860  -26.752346 4.603198e-07 8.845885e-09 2.119726e-13
    #>    -50.573494  -53.964230 5.968447e-16 5.690468e-17 3.594617e-31
    #>    -84.273153  -56.501840 4.278122e-26 9.800567e-18 9.605111e-35
    #>    -71.793098  -35.123667 2.444129e-22 2.671301e-11 7.135848e-22
    #>    -37.782215  -12.619332 4.230777e-12 1.589287e-04 2.525834e-08
    #>    -23.235905  -70.039402 1.012268e-07 8.242125e-22 1.024687e-14
    #>      1.553773    3.073928 2.935839e+00 8.420627e+00 7.952611e+01
    #>    -64.351072  -91.408807 4.250078e-20 3.042336e-28 1.806317e-39
    #>     -3.234803  -88.877355 1.062251e-01 1.758936e-27 1.128378e-02
    #> 1 -100.000000   10.000000 7.888609e-31 1.024000e+03 1.048576e+06
    #> 1   10.000000  -38.697957 1.024000e+03 2.242612e-12 1.048576e+06
    #> 1  -70.487451  -24.761405 6.041754e-22 3.516205e-08 1.236370e-15
    #> 1  -43.828184  -58.653331 6.403270e-14 2.205911e-18 4.100187e-27
    #> 1  -40.534528 -100.000000 6.279012e-13 7.888609e-31 3.942599e-25
    #> 1  -57.840985  -80.635742 3.873722e-18 5.323816e-25 1.500572e-35
    #> 1  -71.835137   10.000000 2.373936e-22 1.024000e+03 1.048576e+06
    #> 1   10.000000 -100.000000 1.024000e+03 7.888609e-31 1.048576e+06
    #> 1  -24.133535  -84.591626 5.433531e-08 3.430704e-26 2.952326e-15
    #> 1  -65.787072 -100.000000 1.570789e-20 7.888609e-31 2.467378e-40
    print(which.min(res$y[, 1, drop = FALSE]))
    #> [1] 20
    cbind(res$xbest, res$ybest)
    #>        [,1] [,2]         [,3]
    #> 1 -65.78707 -100 2.467378e-40

### With noise: (this takes some time)

    # 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))
    # noise with replicated evaluations
    res2 <-
      spot(x = NULL,
           fun = fNoise,
           lower = c(-2, -3),
           upper = c(1, 2),
        control = list(
          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(
        funEvals = 40,
        noise = TRUE,
        replicates = 2,
        OCBA = TRUE,
        OCBABudget = 1,
        designControl = list(replicates = 2)
      )
    )
    # Check results with non-noisy function:
    funSphere(res1$xbest)
    #>           [,1]
    #> [1,] 0.9811919
    funSphere(res2$xbest)
    #>         [,1]
    #> [1,] 1.74417
    funSphere(res3$xbest)
    #>          [,1]
    #> [1,] 1.084404

    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.9586859 0.4066022  0.1312009
    #>  [2,] 0.9586859 0.4066022  1.4909468
    #>  [3,] 0.9586859 0.4066022  3.3136663
    #>  [4,] 0.9586859 0.4066022 -0.4300930
    #>  [5,] 0.9586859 0.4066022  1.0226966
    #>  [6,] 0.9586859 0.4066022  0.9371333
    #>  [7,] 0.9586859 0.4066022  2.6259971
    #>  [8,] 0.9586859 0.4066022  0.1025484
    #>  [9,] 0.9586859 0.4066022  1.5809822
    #> [10,] 0.9586859 0.4066022  2.7813519

### 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: -1.296329 = -1.296329. Should be different:  -0.889934"

### 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 2.619386 2.725642    2
     res$ybest
    #>           [,1]
    #> [1,] 0.6777176

### 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: 0.988475, time: 8.426000"

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.988475, time: 0.800000"

### 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
    #>          X3         X4
    #>    8.648076 81.4784571
    #>   71.771945 66.5887761
    #>   44.933187 41.8506996
    #>   14.297134 39.5437814
    #>   25.642638 58.9784849
    #>   56.561623 79.4369705
    #>   69.785541 27.2369075
    #>   92.321611 93.7035707
    #>   32.408116  7.8101754
    #>   87.968361 10.1114951
    #> 1  3.152969  3.1352284
    #> 1  3.152969  3.1352284
    #> 1  3.152969  3.1352284
    #> 1  3.152969  3.1352284
    #> 1  3.377815  0.5194955
    #> 1  3.377815  0.5194955
    #> 1  3.377815  0.5194955
    #> 1  3.377815  0.5194955
    #> 1  3.377815  0.5194955
    #> 1  3.377815  0.5194955
    ## 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 12
    #>  $ xbest   : num [1, 1:2] 0 0
    #>  $ ybest   : num 0
    #>  $ x       : num [1:33, 1:2] -4.135 2.177 -0.507 -3.57 -2.436 ...
    #>   ..- attr(*, "dimnames")=List of 2
    #>   .. ..$ : chr [1:33] "" "" "" "" ...
    #>   .. ..$ : NULL
    #>  $ 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 ...
    #>   .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. ..$ : chr [1:19] "" "" "" "" ...
    #>   .. .. ..$ : NULL
    #>   ..$ 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 ...
    #>   .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. ..$ : chr [1:19] "" "" "" "" ...
    #>   .. .. ..$ : NULL
    #>   ..$ normalizexmin   : num [1:2] -4.14 -4.22
    #>   ..$ normalizexmax   : num [1:2] 4.23 4.37
    #>   ..$ dmodeltheta     : num [1:2] 3.54 1.71
    #>   ..$ Lambda          : num -5.96
    #>   ..$ dmodellambda    : num 1.09e-06
    #>   ..$ Theta           : num [1:2] 0.549 0.234
    #>   ..$ yonemu          : num [1:19, 1] -11.4 -30.9 -37.5 -24.5 -31.6 ...
    #>   ..$ ssq             : num 174
    #>   ..$ mu              : num 38.4
    #>   ..$ Psi             : num [1:19, 1:19] 1 0.127 0.357 0.654 0.768 ...
    #>   ..$ Psinv           : num [1:19, 1:19] 9.51 -2.16 203.14 3.54 -28.86 ...
    #>   ..$ nevals          : num 1200
    #>   ..$ like            : num [1, 1] -16.8
    #>   ..$ returnCrossCor  : logi FALSE
    #>   ..$ min             : num 0.0988
    #>   ..- 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 34
    #>   ..$ funEvals             : num 20
    #>   ..$ lower                : num [1:2] -5 -5
    #>   ..$ upper                : num [1:2] 5 5
    #>   ..$ types                : chr [1:2] "numeric" "numeric"
    #>   ..$ design               :function (x = NULL, lower, upper, control = list())
    #>   ..$ designControl        :List of 1
    #>   .. ..$ types: chr [1:2] "numeric" "numeric"
    #>   ..$ subsetSelect         :function (x, y, control)
    #>   ..$ subsetControl        : list()
    #>   ..$ direct               : logi FALSE
    #>   ..$ model                :function (x, y, control = list())
    #>   ..$ modelControl         :List of 2
    #>   .. ..$ modelInitialized: logi FALSE
    #>   .. ..$ types           : chr [1:2] "numeric" "numeric"
    #>   ..$ multiStart           : num 1
    #>   ..$ optimizer            :function (x = NULL, fun, lower, upper, control = list(), ...)
    #>   ..$ optimizerControl     :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
    #>   ..$ returnFullControlList: logi TRUE
    #>   ..$ plots                : logi FALSE
    #>   ..$ progress             : logi FALSE
    #>   ..$ OCBA                 : logi FALSE
    #>   ..$ OCBABudget           : num 1
    #>   ..$ replicates           : num 1
    #>   ..$ transformFun         : logi(0)
    #>   ..$ xNewActualSize       : num 0
    #>   ..$ noise                : logi FALSE
    #>   ..$ replicateResult      : logi FALSE
    #>   ..$ seedFun              : logi NA
    #>   ..$ seedSPOT             : num 1
    #>   ..$ infillCriterion      : NULL
    #>   ..$ verbosity            : num 0
    #>   ..$ time                 :List of 3
    #>   .. ..$ maxTime  : num Inf
    #>   .. ..$ startTime: POSIXct[1:1], format: "2022-03-17 23:26:16"
    #>   .. ..$ endTime  : POSIXct[1:1], format: "2022-03-17 23:26:17"
    #>   ..$ maxTime              : num Inf
    #>   ..$ tolerance            : num 1.49e-08
    #>   ..$ fixError             : logi FALSE

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