R/corfbed.network.R

In MXM: Feature Selection (Including Multiple Solutions) and Bayesian Networks

corfbed.network <- function(x, threshold = 0.05, symmetry = TRUE, nc = 1) {
  
  dm <- dim(x)
  n <- dm[1]
  D <- dm[2]
  G <- matrix(0, D, D)
  counter <- 0
  x <- Rfast::standardise(x)
  if ( nc <= 1  || is.null(nc) ) {
    pa <- proc.time()
    for (i in 1:D) {
      id <- 1:D
      id <- id[-i] 
      a <- Rfast::cor.fbed(x[, i], x[, -i], ystand = FALSE, xstand = FALSE, alpha = threshold)$res[, 1]
      sel <- id[a]        
      G[i, sel] <- 1 
      counter <- counter + sum(D - 0:a)
    } 
    runtime <- proc.time() - pa
    
  } else {
    pa <- proc.time() 
    cl <- makePSOCKcluster(nc)
    registerDoParallel(cl)
    sel <- numeric(n)
    mod <- foreach(i = 1:D, .combine = rbind, .export = "cor.fbed", .packages = "Rfast" ) %dopar% {
      id <- 1:D
      id <- id[-i] 
      sela <- numeric(D)  
      a <- Rfast::cor.fbed(x[, i], x[, -i], ystand = FALSE, xstand = FALSE, alpha = threshold)$res[, 1]
      sel <- id[a]  
      sela[sel] <- 1
      return( sum(D - 0:a), sela)
    }    
    stopCluster(cl)
    G <- as.matrix(mod[, -1])
    counter <- sum(mod[, 1])
    runtime <- proc.time() - pa
  }	   
  
  diag(G) <- 0	   
  if ( symmetry ) {
    a <- which( G == 1  &  t(G) == 1 ) 
    G[ -a ] <- 0
  } else {
    G <- G + t(G)
    G[ G > 0 ] <- 1
  }
  
  info <- summary( Rfast::rowsums(G) )
  density <- sum(G) / D / ( D - 1 )  
  if ( is.null( colnames(x) ) ) {
    colnames(G) <- rownames(G) <- paste("X", 1:D, sep = "")
  } else  colnames(G) <- rownames(G) <- colnames(x) 
  
  list(runtime = runtime, density = density, info = info, ntests = counter, G = G)
}