library(GaSP) # install from CRAN
library(EGO)  # install manually for the .tar.gz

# log of Goldstein-Price function
log_goldpr  <- function(xx){  #min is log(3) at all x = (0,-1)
  x1 <- xx[1]
  x2 <- xx[2]
  fact1a <- (x1 + x2 + 1)^2
  fact1b <- 19 - 14*x1 + 3*x1^2 - 14*x2 + 6*x1*x2 + 3*x2^2
  fact1 <- 1 + fact1a*fact1b
  fact2a <- (2*x1 - 3*x2)^2
  fact2b <- 18 - 32*x1 + 12*x1^2 + 48*x2 - 36*x1*x2 + 27*x2^2
  fact2 <- 30 + fact2a*fact2b
  y <- fact1*fact2
  return(log(y))
}

# Min and max bounds for Goldstein-Price
x_min <- c(-2, -2); x_max <- c(2, 2);

# Set up optimization region
d <- 2 # Dimension of x
x_names <- paste("x",1:d, sep = "")
x_support <- rep("Continuous", d)
x_levels <- rep(0, d)
x_describe<- DescribeX(x_names, x_min, x_max, x_support, x_levels)
x_describe
# NumberLevels not used for continuous variables - don't worry

# Get ready for EGO 
# EGO will create a starting design of size n_design 
# if x_design and y_design are NULL
# n_design = 20 for 20 iterations
set.seed(1234)
ego_init <- Initialize(x_design = NULL, y_design = NULL, n_design =  20, 
                       x_describe = x_describe, fun = log_goldpr)

control <- EGO.control(alg = "genoud", rel_tol = 0, 
                      genoud_control = list(pop_size = 3*2^d, 
                                            max_generations = 5))

# Run EGO
log_goldpr_ego <- EGO(log_goldpr, reg_model = ~1, ego_init, x_describe, 
                      nsteps = 20, control = control)
# x for the min y found
log_goldpr_ego$opt_x

# Locations of initial evaluations
plot(log_goldpr_ego$x[1:20, ], pch = 20)
# Locations chosen by EGO
points(log_goldpr_ego$x[21:40, ], pch = "+", cex = 1.5, col = "red")