optim_on_diff_initial_designs, ninit = 8, dim = 2, theta_max = 10*sqrt(m)

library(lhs)
library(laGP)
source("/home/boyazhang/repos/unifdist/code/ud_new.R")

## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.

## 

## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.

## 
## Attaching package: 'gdata'

## The following object is masked from 'package:stats':
## 
##     nobs

## The following object is masked from 'package:utils':
## 
##     object.size

## The following object is masked from 'package:base':
## 
##     startsWith

load("/home/boyazhang/repos/unifdist/code/optim_2d_8.RData")

f <- function(X)
{
  if(is.null(nrow(X))) X <- matrix(X, nrow=1)
  m <- 8.6928
  s <- 2.4269
  x1 <- 4 * X[,1] - 2
  x2 <- 4 * X[,2] - 2
  a <- 1 + (x1 + x2 + 1)^2 * (19 - 14 * x1 + 3 * x1^2 - 14 *
                                x2 + 6 * x1 * x2 + 3 * x2^2)
  b <- 30 + (2 * x1 - 3 * x2)^2 * (18 - 32 * x1 + 12 * x1^2 +
                                     48 * x2 - 36 * x1 * x2 + 27 * x2^2)
  f <- log(a * b)
  f <- (f - m)/s
  return(f)
}


obj.mean <- function(x, gpi) predGP(gpi, matrix(x, nrow=1), lite=TRUE)$mean


optim.surr <- function(f, ninit, m, stop, tol=1e-4, theta_max, ab = NULL, design)
{
  if(design == "maximin"){
    X <- maximin(ninit, m)
  } else if(design == "dist_beta"){
    X <- bd(ninit, m, T = 10^5, shape1 = 2.22, shape2 = 5.19)$X
  } else if(design == "LHS"){
    X <- randomLHS(ninit, m)
  } else if(design == "random"){
    X <- matrix(runif(ninit*m), ncol = m)
  }
  da <- darg(list(mle=TRUE, max=theta_max,  min = sqrt(.Machine$double.eps)), X) ## max = theta_max = sqrt(m)  min = sqrt(.Machine$double.eps)) 
  y <- f(X)
  gpi <- newGP(X, y, d=0.1, g=1e-7, dK=TRUE)
  if(!is.null(ab)) { da$ab <- ab } ## ab = c(0,0)
  mleGP(gpi, param="d", tmin=da$min, tmax=da$max, ab=da$ab)$msg
  for(i in (ninit+1):stop) { ## loop start from ninit+1
    m <- which.min(y)
    opt <- optim(X[m,], obj.mean, lower=0, upper=1, 
                 method="L-BFGS-B", gpi=gpi)
    ynew <- f(opt$par)
    if(abs(ynew - y[length(y)]) < tol) break;
    updateGP(gpi, matrix(opt$par, nrow=1), ynew)
    mle <- mleGP(gpi, param="d", tmin=da$min, tmax=da$max, ab=da$ab)
    X <- rbind(X, opt$par)
    y <- c(y, ynew)
  }
  deleteGP(gpi)
  
  return(list(X=X, y=y))## output y: first four numbers are f(X_init), the rest is the ymin for each iteration
}

Sys.time()
n <- 50
m <- 2
ninit <- 8
reps <- 100
prog <- matrix(NA, nrow=reps, ncol=n)
for(r in 1:reps) {
  os <- optim.surr(f, ninit, m, n, theta_max = 10*sqrt(m), ab = c(0,0), design = "dist_beta")
  prog[r,1:length(os$y)] <- os$y
  for(i in 2:n) { 
    if(is.na(prog[r,i]) || prog[r,i] > prog[r,i-1]) 
      prog[r,i] <- prog[r,i-1] 
  }
}
prog1 <- prog

prog <- matrix(NA, nrow=reps, ncol=n)
for(r in 1:reps) {
  os <- optim.surr(f,  ninit, m, n, theta_max = 10*sqrt(m), ab = c(0,0), design = "maximin")
  prog[r,1:length(os$y)] <- os$y
  for(i in 2:n) { 
    if(is.na(prog[r,i]) || prog[r,i] > prog[r,i-1]) 
      prog[r,i] <- prog[r,i-1] 
  }
}
prog2 <- prog

prog <- matrix(NA, nrow=reps, ncol=n)
for(r in 1:reps) {
  os <- optim.surr(f, ninit, m, n, theta_max = 10*sqrt(m), ab = c(0,0), design = "LHS")
  prog[r,1:length(os$y)] <- os$y
  for(i in 2:n) { 
    if(is.na(prog[r,i]) || prog[r,i] > prog[r,i-1]) 
      prog[r,i] <- prog[r,i-1] 
  }
}
prog3 <- prog

prog <- matrix(NA, nrow=reps, ncol=n)
for(r in 1:reps) {
  os <- optim.surr(f, ninit, m, n, theta_max = 10*sqrt(m), ab = c(0,0), design = "random")
  prog[r,1:length(os$y)] <- os$y
  for(i in 2:n) { 
    if(is.na(prog[r,i]) || prog[r,i] > prog[r,i-1]) 
      prog[r,i] <- prog[r,i-1] 
  }
}
prog4 <- prog
save.image('/home/boyazhang/repos/unifdist/code/optim_2d_8.RData') 

par(mfcol = c(1,4))
matplot(t(prog1), type="l", col="gray", lty=1, ylab="progress", xlab="its", ylim = c(-3,1.5), main = "dist_beta")
matplot(t(prog2), type="l", col="gray", lty=1, ylab="progress", xlab="its", ylim = c(-3,1.5), main = "maximin")
matplot(t(prog3), type="l", col="gray", lty=1, ylab="progress", xlab="its", ylim = c(-3,1.5), main = "LHS")
matplot(t(prog4), type="l", col="gray", lty=1, ylab="progress", xlab="its", ylim = c(-3,1.5), main = "random")

par(mfcol = c(1,1))
n <- 25
boxplot(prog1[,n], prog2[,n], prog3[,n], prog4[,n], names = c("dist_beta", "maximin", "LHS", "random"), main = paste("boxplot of fmin after", n, "searches"))

n <- 50
boxplot(prog1[,n], prog2[,n], prog3[,n], prog4[,n], names = c("dist_beta", "maximin", "LHS", "random"), main = paste("boxplot of fmin after", n, "searches"))

Sys.time()

## [1] "2018-05-01 13:26:45 EDT"