GP_mse_rand_unifdist_beta_cp

## compare Imse of d on unifdist.rand, random design, maximin and beta(3,6)
library(plgp)
library(mvtnorm)
library(laGP)

library(foreach)
library(doParallel)
source("ud.R")
gpi.mse <- function(method, n, dim, I = 1000, seed = 0){
  set.seed(seed)
  if(method ==1){
    x <- ud.m(n, dim, k = 1, T = 10^5, starter_num = 5)$X
  }
  else if(method == 2){
    x <- matrix(runif(n*dim), ncol = dim)
  }
  else if(method == 3){
    x <- maximin(n, dim)
  }
  else if(method == 4){
    x <- bd(n, dim)$X
  }
  D <- distance(x)
  eps <- sqrt(.Machine$double.eps)
  dtrue <- runif(I, min = eps, max = sqrt(dim))
  dhat <- rep(NA, I)
  for(i in 1:I){
    sigma <- exp(-D/dtrue[i] + diag(eps, n))
    y <- rmvnorm(1, sigma = sigma)
    gpi <- newGP(x, y, d = 0.1, g = eps, dK = T)
    dhat[i] <- mleGP(gpi, param = "d", tmax = 10)$d
    deleteGP(gpi)
  }
  mse <- mean((dtrue - dhat)^2)
  return(c(n, dim, method, mse, seed))
}

A <- expand.grid(method = c(1,2,3,4), n = c(8,16), dim = 2:4)
set <- matrix(rep(t(A),100), ncol=ncol(A), byrow=TRUE)
set <- cbind(set, 1:nrow(set))
colnames(set) <- c("method", "n", "dim", "seed")
cl <- makeCluster(8)
registerDoParallel(cl)
out <- foreach(i = 1:nrow(set), .errorhandling = "remove", .combine='rbind', .packages=c('plgp', 'mvtnorm', 'laGP' )) %dopar% { 
              gpi.mse(method =set[i,1], n = set[i,2], dim = set[i,3], seed = set[i,4]) }

stopCluster(cl)
save.image('GP_mse_rand_unifdist_cp_2.RData') 

#-------------------------------------plot-----------------------------------------------------------------------------------                                                                                                                                                                                                                 
load('/home/boyazhang/repos/unifdist/code/GP_mse_rand_unifdist_cp_2.RData')
out <- as.data.frame(out)
colnames(out) <- c("n", "dim", "method", "ise")
out$method <- replace(out$method, out$method == 1, "5-start unifdist" )
out$method <- replace(out$method, out$method == 2, "random" )
out$method <- replace(out$method, out$method == 3, "maximin" )
out$method <- replace(out$method, out$method == 4, "beta(3,6)" )

par(mfcol = c(1,3))
boxplot(ise~method, data = out[out$dim == 2 & out$n == 16,], main = "dim = 2, n=16, ise comparison", ylab = "MSE")
boxplot(ise~method, data = out[out$dim == 3 & out$n == 16,], main = "dim = 3, n=16, ise comparison", ylab = "MSE")
boxplot(ise~method, data = out[out$dim == 4 & out$n == 16,], main = "dim = 4, n=16, ise comparison", ylab = "MSE")

boxplot(ise~method, data = out[out$dim == 2 & out$n == 8,], main = "dim = 2, n=8, ise comparison", ylab = "MSE")
boxplot(ise~method, data = out[out$dim == 3 & out$n == 8,], main = "dim = 3, n=8, ise comparison", ylab = "MSE")
boxplot(ise~method, data = out[out$dim == 4 & out$n == 8,], main = "dim = 4, n=8, ise comparison", ylab = "MSE")