if (!require("BDgraph")) install.packages("BDgraph")
if (!require("mvtnorm")) install.packages("mvtnorm")

if (!require("dplyr")) install.packages("dplyr")


par_cor_CI <- function(X, alpha){

  # X: data frame
  if (!require("BDgraph")) install.packages("BDgraph")
  if (!require("qgraph")) install.packages("qgraph")
  if (!require("glasso")) install.packages("glasso")
  if (!require("mvtnorm")) install.packages("mvtnorm")


  n <- nrow(X)
  p <- ncol(X)
  I = diag(1,p)
  ## compute maximum likelihood estimator for covariance matrix
  mle_cov <- cov(X, use = "pairwise.complete")#crossprod(scale(X, scale = F)) / n

  ## compute maximum likelihood estimator of precision matrix
  ## (inverse covariance matrix)
  mle_inv <- solve(mle_cov)

  ## standardize and revese sign = partial correaltions
  pc  <-  as.matrix(wi2net(mle_inv))



  mle_parcors <- ci_par_cor(alpha = alpha, par_cors = pc, n = n, s = p - 1)

  mle_inv <- mle_parcors$sig_mat * mle_inv
  #if(det(mle_inv <= 0)) mle_inv <- mle_inv + (abs(min(eigen(mle_inv)$values)) + 10)*I
  list(mle_parcors = mle_parcors, mle_inv = mle_inv )

}


ci_par_cor <- function(alpha, par_cors, s, n) {

  # n: sample size
  # s: p - 1 (controlled for)
  # alpha: confidence level
  # par_cors: partial correlations
  mat <- matrix(0,nrow =  s + 1, ncol = s + 1)
  CI_ls <- list()
  par_cor <- par_cors[upper.tri(par_cors)]
  cov <- list()
  for(i in 1:length(par_cor)){
    # crtiical value
    z_crit <- qnorm(1 - alpha/2)
    # standard error
    se <- sqrt(1/((n - s - 3)))
    # z transformation
    z <- log((1 + par_cor[i])/(1 - par_cor[i]))/2
    # z lower bound
    Z_L <- z - z_crit * se
    # Z upper bound

    Z_U <- z + z_crit * se
    rho_L <- (exp(2*Z_L) - 1)/(exp(2*Z_L) + 1)
    rho_U <- (exp(2*Z_U) - 1)/(exp(2*Z_U) + 1)
    CI <- c(rho_L, rho_U)
    CI_ls[[i]] <- CI
    cov[[i]] <- ifelse(CI[1] < 0 & CI[2] > 0, 0, 1)
  }

  c_dat <- do.call(rbind.data.frame, CI_ls)
  colnames(c_dat) <- c("low", "up")
  c_dat$pcar <- unlist(par_cor)
  diag(mat) <- 1
  mat[upper.tri(mat)] <- unlist(cov)
  mat <- as.matrix(forceSymmetric(mat) )
  list(sig_mat = mat, par_cors = par_cors, par_sig = mat * par_cors,
       cis = c_dat, cov_prob = unlist(cov))
}

sim_coverage <- function(n, alpha){
  main <- bdgraph.sim(p = 20, n = n, prob = 0, b = 3)
  dat <-  rmvnorm(n, mean = rep(0, 20), sigma = main$sigma)
  test_run <-  par_cor_CI(dat,  alpha = alpha)
  spec <- compare(main, test_run$mle_parcors$sig_mat)[6,2]
  cov <- sum(ifelse(test_run$mle_parcors$cis$low < 0 & test_run$mle_parcors$cis$up > 0, 1, 0)) / 190
  cis <- test_run$mle_parcors$cis
  list(cis = cis, spec = spec, cov = cov)

}