library(XRJulia)
library(MASS)

bb_fit <- juliaEval("
function BB_est(x, sims)
X = center_data(x)
n = size(X, 1)
p = size(X, 2)

offdiag = (p * (p - 1)) / 2

loop_max = convert(Int64, offdiag)

temp    = Exponential(1)
exp_drw = rand(temp, n * sims)

rexp_mat = reshape(exp_drw, sims, n)
dir_weights = rexp_mat ./sum(rexp_mat, 2)

A = zeros(sims, offdiag)
B = zeros((1, offdiag))
C = zeros(sims, p)
mat_point = zeros(p,p)
mat_sig = zeros(p, p)
idx = tril(trues(size(mat_point)), -1);

for i = 1:sims
t = dir_weights[i, 1:n]
d = X .* sqrt.(t)
bb_cov = (d' * d)
    inv_cov =  inv(bb_cov)
    d = diag(inv_cov)
    inv_cov[diagind(inv_cov)] = 0

    test = inv_cov[idx]
    rmv_0 = filter!(x->x≠ 0,test)
    A[i,:] = rmv_0
    C[i,:] = d
    end

   for i in 1:loop_max
    temp1 =  quantile(A[:,i], 0.025)
    temp2 =  quantile(A[:,i], 0.975)
    B[1,i] = ifelse(temp1 < 0 && temp2 > 0, 0, 1)
 end

    mat_sig[idx] = B
    mat_point[idx] = mean(A, 1)
    mat_point[diagind(mat_point)] = mean(C, 1)
    force_sym(mat_sig), force_sym(mat_point)
     mat_sig, mat_point
end")

n<- 200
p <- 20
X <- mvrnorm(n, rep(0, p), Sigma = diag(p))


bb_fit_function=JuliaFunction(bb_fit)


test <-bb_fit_function(main$data, 5000)


JLreg=juliaGet(test)