using MATLAB
using Distributions


function hs_ggm(X, nmc, interval)
    S = X' * X
   p = size(S, 2)
   sims = convert(Int64, nmc)
   @mput X; @mput S; @mput nmc
   mat"n = size(X, 1)"
   mat"[GHS_omega_save,GHS_lambda_sq_save,GHS_tau_sq_save] = GHS(S, n, 100, nmc)"
   mat_point = mat"mean(GHS_omega_save,3)"
   samps =   @mget GHS_omega_save
   #n = @mget n
   offdiag = (p * (p - 1)) / 2

   mat_select =  zeros(p, p)
   idx = tril(trues(size(mat_point)), -1)
   A = zeros(nmc, offdiag)
   for i = 1:sims
       samp = samps[:,:,i]
       test = samp[idx]
       A[i,:] = test
   end
  interval_select = ggm_select(A, interval)
  mat_select[idx] = interval_select
  force_sym(mat_select), force_sym(mat_point)
end


function ggm_select(post_mat, interval)
   low = (1 - interval) / 2
   up =  1 - low
   select = zeros(size(post_mat, 2))
   for i in 1:size(post_mat, 2)
   bounds = quantile(post_mat[:,i], [low, up]);
   select[i] =  ifelse(bounds[1] < 0 && bounds[2] > 0, 0, 1)
end
   select
end

function force_sym(x)
   x + tril(x,-1).'
end


############################
###### call from R #########
############################
fit <- julia_call("hs_ggm", X, 1000, .50)