ud.m

ud.m <- function(n, m, k, T = 10^4, dmax = 1, starter_num = 20){
  ksv <- rep(NA, starter_num)
  designs <- list()
  for(i in 1:starter_num){
   des <- ud(n,m,k,T,dmax)
   designs[[i]] <- des$X ## only output design and correponding k-s distance, could add more
   ksv[i] <- des$ksopt
  }
  min_loc <- which(ksv==min(ksv))
  des_min <- designs[[min_loc]]
  return(list(ksopt = min(ksv), X = des_min))
}

## compare mse of d on unifdist.rand and random design
library(plgp)
library(mvtnorm)
library(laGP)

library(foreach)
library(doParallel)
source("ud.R")
gpi.mse <- function(method, n, dim, drange = c(0,1), I = 1000){
  if(method ==1){
    x <- ud.m(n,dim, k = 1, T = 10^5, starter_num = 10)$X
  }
  else if(method == 2){
    x <- matrix(runif(n*dim), ncol = dim)
  }
  D <- distance(x)
  eps <- sqrt(.Machine$double.eps)
  dtrue <- runif(I, min = drange[1], max = drange[2])
  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))
}

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

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

  n dim method          
result.1  8   2      1 0.1587011
result.2  8   2      2 0.1933107
result.3  8   3      1 0.8462921
result.4  8   3      2 0.5986832
result.5  8   4      1 1.1894464
result.6  8   4      2 1.0139950
result.7  8   2      1 0.1713684
result.8  8   2      2 0.1495081
result.9  8   3      1 0.3153247
result.10 8   3      2 0.6629104
result.11 8   4      1 1.3678913
result.12 8   4      2 0.9040594