1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
tol.mc.fix <- function(rs,ITER,parallel,verbose,...){
fd <- rs$f
gd <- rs$g
pool <- suppressWarnings(superimpose(fd$pp,gd$pp))
nf <- npoints(fd$pp)
ng <- npoints(gd$pp)
res <- dim(fd$z)[1]
fh <- fd$h0
gh <- gd$h0
indx <- 1:npoints(pool)
if(is.null(fd$q)){
edg <- "none"
} else if(is.im(fd$q)){
edg <- "uniform"
} else {
edg <- "diggle"
}
logt <- !all(rs$rr>=0)
rmat <- as.matrix(rs$rr)
mcmat <- matrix(1,res,res)
if(is.null(parallel)){
if(verbose) pb <- txtProgressBar(0,ITER-1,style=3)
for(i in 1:(ITER-1)){
shuff <- sample(indx)
ftemp <- bivariate.density(pool[shuff[1:nf]],h0=fh,resolution=res,edge=edg)
gtemp <- bivariate.density(pool[shuff[(nf+1):(nf+ng)]],h0=gh,resolution=res,edge=edg)
rtemp <- as.matrix(risk(ftemp,gtemp,log=logt,verbose=FALSE,...)$rr)
mcmat <- mcmat+(rtemp>=rmat)
if(verbose) setTxtProgressBar(pb,i)
}
if(verbose) close(pb)
} else {
ncores <- detectCores()
if(verbose) cat(paste("Running MC iterations on",parallel,"/",ncores,"cores..."))
if(parallel>ncores) stop("Parallel cores requested exceeds available count")
registerDoParallel(cores=parallel)
mclist <- foreach(i=1:(ITER-1),.packages=c("spatstat","sparr")) %dopar% {
shuff <- sample(indx)
ftemp <- bivariate.density(pool[shuff[1:nf]],h0=fh,resolution=res,edge=edg)
gtemp <- bivariate.density(pool[shuff[(nf+1):(nf+ng)]],h0=gh,resolution=res,edge=edg)
rtemp <- as.matrix(risk(ftemp,gtemp,log=logt,verbose=FALSE,...)$rr)
return(rtemp>=rmat)
}
mcmat <- Reduce("+",mclist) + mcmat
if(verbose) cat("Done.\n")
}
return(mcmat/ITER)
}
|
