mykest plot of random designs

gpi.mse <- function(n, dim, I = 1000, seed = 0){
  set.seed(seed)
  x <- matrix(runif(n * dim), ncol = dim) # random design
  D <- distance(x)
  eps <- sqrt(.Machine$double.eps)
  dtrue <- runif(I, min = 0, 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(list(X = x, n=n, dim = dim, mse = mse, seed = seed))
}


library(plgp)
library(mvtnorm)
library(laGP)
source("ud.R")

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

stopCluster(cl)
save.image("space-fillingness_functions.RData")

source("/home/boyazhang/repos/unifdist/code/mykest.R")
load("/home/boyazhang/repos/unifdist/code/space-fillingness_functions.RData")
l <- length(out)
msev <- rep(NA, l)
nv <- rep(NA, l)
dimv <- rep(NA, l)
design <- list()
for(i in 1:l){
  design[[i]] <- out[[i]][[1]]
  nv[i] <- out[[i]][[2]]
  dimv[i] <- out[[i]][[3]]
  msev[i] <- out[[i]][[4]]
}

# Ripley’s K Function, the expected number of events within a distance r from an arbitrary event (over: cluster; below: regular process)
f <- function(n, dim,range = c(0, 5), overlay = F, color = 1){
  mse1 <- msev[nv==n & dimv == dim]
  subind <- c(1:l)[nv==n & dimv == dim]
  design_sub <- list()
  for(i in 1:length(subind)){
    design_sub[[i]] <- design[[subind[i]]]
  }
  rank1 <- c(1:length(mse1))[rank(mse1) <= range[2] & rank(mse1) > range[1]]
  for( i in rank1){
    design0 <- design_sub[[i]]
    if( i == rank1[1]){mykest(x = design0, add = overlay, col = color)}
    else{ mykest(design0, add = T, col = color)}
  }
}
par(mfcol = c(1,2))
# f(8,2,c(0,1))
# f(8,2,c(9999,10000))
# f(8,2,c(0,5))
# f(8,2,c(9995,10000))
# f(8,2,c(0,20))
# f(8,2,c(9980,10000))
# f(16,2,c(0,5))
# f(16,2,c(9995,10000))
# f(16,2,c(0,20))
# f(16,2,c(9980,10000))

f(8,2,c(0,20))
f(8,2,c(9980,10000), overlay = T, color = 2)
f(16,2,c(0,20))
f(16,2,c(9980,10000), overlay = T, color = 2)

f(8,3,c(0,20))
f(8,3,c(9980,10000), overlay = T, color = 2)
f(16,3,c(0,20))
f(16,3,c(9980,10000), overlay = T, color = 2)

f(8,4,c(0,20))
f(8,4,c(9980,10000), overlay = T, color = 2)
f(16,4,c(0,20))
f(16,4,c(9980,10000), overlay = T, color = 2)

# f(16,2,c(0,50))
# f(16,2,c(9950,10000))
# f(8,3,c(0,50))
# f(8,3,c(9950,10000))
# f(16,3,c(0,50))
# f(16,3,c(9950,10000))
# f(8,4,c(0,50))
# f(8,4,c(9950,10000))
# f(16,4,c(0,50))
# f(16,4,c(9950,10000))