R/tol.mc.fix.R

In sparr: Spatial and Spatiotemporal Relative Risk

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