R/mmtabu.boot.R

In pchc: Bayesian Network Learning with the PCHC and Related Algorithms

mmtabu.boot <- function(x, method = "pearson", max_k = 3, alpha = 0.05, ini.stat = NULL,
                        R = NULL, tabu = 10, score = "bic-g", blacklist = NULL,
                        whitelist = NULL, B = 200, ncores = 1) {

  mod <- pchc::mmtabu(x, method = method, max_k = max_k, alpha = alpha, ini.stat = ini.stat, R = R,  tabu = tabu,
                      score = score, blacklist = blacklist, whitelist = whitelist)
  dm <- dim(x)
  n <- dm[1]   ;     p <- dm[2]

  runtime <- proc.time()
  if ( ncores == 1 ) {
    Gboot <- matrix(0, p, p)
    for (i in 1:B) {
      id <- Rfast2::Sample.int(n, n, replace = TRUE)
      gb <- pchc::mmtabu(x[id, ], method = method, max_k = max_k, alpha = alpha, tabu = tabu,
                         score = score, blacklist = blacklist, whitelist = whitelist)
      Gboot <- Gboot + pchc::bnmat(gb$dag)
    }  ## end for (i in 1:B)
  } else {
    cl <- parallel::makePSOCKcluster(ncores)
    doParallel::registerDoParallel(cl)
    Gboot <- foreach::foreach(i = 1:B, .combine = rbind,
                     .export = c("mmtabu", "Sample.int"),
                     .packages = c("pchc", "Rfast2") ) %dopar% {
      id <- Rfast2::Sample.int(n, n, replace = TRUE)
      gb <- pchc::mmtabu(x[id, ], method = method, alpha = alpha, tabu = tabu,
                          score = score, blacklist = blacklist, whitelist = whitelist)
      return( as.vector( pchc::bnmat(gb$dag) ) )
    }
    parallel::stopCluster(cl)
    Gboot <- Rfast::colsums(Gboot)
    Gboot <- matrix(Gboot, nrow = p, ncol = p)
  }
  runtime <- proc.time() - runtime

  list(mod = mod, Gboot = Gboot/B, runtime = runtime)
}