space-fillingness functions

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")

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: nlme

## Loading required package: rpart

## 
## spatstat 1.55-0       (nickname: 'Stunned Mullet') 
## For an introduction to spatstat, type 'beginner'

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)){
  if(dim == 2){w <- owin(c(0,1),c(0,1))}
    else if(dim == 3){w <- owin(c(0,1),c(0,1),c(0,1))}
    else if(dim == 4){w <- owin(c(0,1),c(0,1),c(0,1),c(0,1))}
  mse1 <- msev[nv==n & dimv == dim]
  rank1 <- c(1:length(mse1))[rank(mse1) <= range[2] & rank(mse1) > range[1]]
  for( i in rank1){
    design0 <- design[[i]]
    mypattern <- as.ppp(X = design0, W = w)
    k1 <- Kest(mypattern, correction = "isotropic", rmax = sqrt(2)/2)
    if( i == rank1[1]){plot(k1, main = paste0("n = ", n, ", dim = ", dim, ", rank: ", range[1], "~", range[2]))}
    else{ plot(k1, add = T)}
  }
}
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,3,c(0,5))
# f(8,3,c(9995,10000))
# f(16,3,c(0,5))
# f(16,3,c(9995,10000))
# f(8,4,c(0,5))
# f(8,4,c(9995,10000))
# f(16,4,c(0,5))
# f(16,4,c(9995,10000))

# G(r) = #(ri ≤ r)/n
# proportion of nearest neighbor distances that are less than r (over when r's small: cluster; below when r's small: regular process)
par(mfcol = c(1,2))
x <- matrix(runif(8*2), ncol = 2)
w <- owin(c(0,1),c(0,1)) # create a window
ppp <- as.ppp(X = x, W = w)
plot(ppp)
G <- Gest(ppp)
plot(G$r,G$han, type = "l")