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 = 45;
    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


% load cancer group data
% the R code Data_processing.R contains how dtcancer.mat and dtnocancer.mat
% are generated

load 'dtcancer.mat';

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

gamma_0 = 0.2;
Niter = 10000;

[sigmasq_est_cancer, CI_L_cancer, CI_m_cancer, CI_R_cancer, Conv_diag_cancer]...
    = Bayes_Neigh_Sel(dtcancer, gamma_0, Niter, sigmasq);


% load no-cancer group data
load 'dtnocancer.mat';

[sigmasq_est_no_cancer, CI_L_no_cancer, CI_m_no_cancer, CI_R_no_cancer, Conv_diag_no_cancer]...
    = Bayes_Neigh_Sel(dtnocancer, 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.


save res_real_data;

%%For post-processing of the results see Data_processing.m

delete(myPool);



