unreg_fit <- function(x, n_cl){
# x: data
# n_cl: number of clusters

# clusters
cl <- makeCluster(n_cl)
registerDoSNOW(cl)

# matrices
mat_p <- mat_aic <- mat_bic <- matrix(0, ncol = ncol(x), nrow = ncol(x))

# scale data
x <- as.data.frame(scale(x))
# number of columns
p <- ncol(x)
# column names
colnames(x) <- letters[1:p]

# name matrix columns
colnames(mat_p) <- colnames(mat_aic) <- colnames(mat_bic)  <- letters[1:p]
temp <-   noquote(paste(letters[1:p], "~ .", sep = " "))
# parallel loop
fit <- foreach(i = 1:p, .packages =  c("SignifReg", "dplyr"), rbind) %dopar%{

                 scope <- as.formula(temp[i])

                 # forward selection
                 #  1) p-value with no correction
                 coef_p <- coefficients(SignifReg(scope = scope, data = x, alpha = 0.05,
                                        direction = "forward", correction = "None"))[-1]
                 names_p <- names(coef_p)

                 #  2) aic
                 coef_aic <- coefficients(SignifReg(scope = scope, data = x, criterion = "AIC",
                                          direction = "forward", alpha = 0.99, correction = "None"))[-1]
                 names_aic <- names(coef_aic)

                 #  3) bic
                 coef_bic <- coefficients(SignifReg(scope = scope, data = x, criterion = "BIC",
                                          alpha = 0.99, direction = "forward", correction = "None"))[-1]
                 names_bic <- names(coef_bic)

                 # returned values
                 p <-   data.frame(coef_p, names_p)
                 aic <- data.frame(coef_aic, names_aic)
                 bic <- data.frame(coef_bic, names_bic)
                 list(p = p, aic = aic, bic = bic)
                 }

# storage lists
p_l <- aic_l <- bic_l <- list()

for(i in 1:length(fit)){
  p_l[[i]] <- fit[[i]]$p
  aic_l[[i]] <- fit[[i]]$aic
  bic_l[[i]] <- fit[[i]]$bic
}

# betas into data frame
p_beta <- suppressMessages(reshape2::melt(p_l))
aic_beta <- suppressMessages(reshape2::melt(aic_l))
bic_beta <- suppressMessages(reshape2::melt(bic_l))

# place betas in matrix
for(i in 1:ncol(x)){
  temp_p <-   p_beta %>% filter(L1 == i)
  temp_aic <- aic_beta %>% filter(L1 == i)
  temp_bic <- bic_beta %>% filter(L1 == i)

  mat_p[i, as.character(temp_p$names_p)] <- temp_p$value
  mat_aic[i, as.character(temp_aic$names_aic)] <- temp_aic$value
  mat_bic[i, as.character(temp_bic$names_bic)] <- temp_bic$value
}
stopCluster(cl)

# compute partial correlations
p_pcor <- parcor::Beta2parcor(mat_p)
aic_pcor <- parcor::Beta2parcor(mat_aic)
bic_pcor <- parcor::Beta2parcor(mat_bic)

# return these
list(p_pcor = p_pcor, aic_pcor = aic_pcor, bic_pcor = bic_pcor)

}