%clear all;

% This is an example of using a quasi-Bayesian neighborhood selection
% method for fitting Gaussian graphical models with an elastic net prior. 
% For more detail on the methods see the paper
% "A scalable quasi-Bayesian framework for Gaussian graphical models", arXiv: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This code requires Matlab Statistics and Parallel toolbox.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(getenv('PBS_NP'))
    NP = 2;
    profilename = 'local';
else
    NP = str2double(getenv('PBS_NP'));
    profilename = 'current';
end

myPool = parpool(profilename, NP);
if isempty(myPool)
    error('pctexample:backslashbench:poolClosed', ...
         ['This example requires a parallel pool. ' ...
          'Manually start a pool using the parpool command or set ' ...
          'your parallel preferences to automatically start a pool.']);
    exit
end


% Generate simulation data
%p=30; n=200;
%s = 2.0*p; %number of non-zeros
%[X,theta_star]=GenerateData_ggm(p,s,n);

sigmasq = []; %we do not know the diagonal. It is then estimated by empirical Bayes

gamma_0 = 0.2;
Niter = 2000;
[sigmasq_est, CI_L, CI_m, CI_R, Conv_diag] = Bayes_Neigh_Sel(X, gamma_0, Niter, sigmasq);

%Output:
% sigmasq_est: p by 1: inverse of the diagonal elements of the precision
% matrix
% CI_m: p by p, estimated (sparse) precision matrix (with zero diagonal).
% The diagonal is estimated separately in sigmasq_est
% CI_L and CI_R: p by p, boundaries of the 95% posterior intervals. 
% Conv_diag: an integer. The number of neigh selection regression MCMC
% where Geweke convergence fail.

delete(myPool);


