
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% A simulation example of the random design where
%%% the predictor matrix X ~ N(0, Omega/n) and Omega is a tri-diagonal matrix.
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear

%% Generate X
p = 5000;
theta =0.95;
n = round(p^theta); 
%%% Construct a tri-diagonal matrix
rho = 0.45;
II = [1:p,1:(p-1),2:p];
JJ = [1:p,2:p,1:(p-1)];
VV = [ones(1,p),rho*repmat([-1,1,1],1,(p-2)/3),-rho,rho*repmat([-1,1,1],1,(p-2)/3),-rho];%,rho2*ones(1,2*(p-2))];
Omega = sparse(II,JJ,VV);
shuffler = randperm(p);
Omega = Omega(shuffler,shuffler);
OmegaRoot =  chol(Omega);
%%% generate X from a gaussian random design------
 X = randn(n,p);
 X =  X* OmegaRoot/sqrt(n);

%% generate the signal beta 
v = 0.4;
r = 2;
tau = sqrt(2*r*log(p))
sp = p^(1-v)

%%% generate the support of the signals
sigseed = rand(p,1);
signal = (abs(sigseed-0.5)<p^(-v)/2);
sumsig = sum(signal)

%%% generate the signs of the signals
signbeta = 2*(randi(2,[p,1])-1.5).*signal;

%%% generate signals with different strength.
tt = 1;
uun = randn(p,tt).*(rand(p,1)<0.2);
uun2 = (sum(uun.^2,2)/tt/6);
beta = tau*(1+uun2).*signbeta;

%% generate y
z = randn(n,1);
y = X*beta + z;

%% Optional block: Lasso procedure for this problem. 
% %%% you can compare our result with the lasso's.
% qV = v+(r-v).*(1:50)./50;
% hamLassoV = zeros(1,50);
% for jjj=1:50
% qlasso=qV(jjj);
% lambdaLasso =  sqrt(2*qlasso*log(p));
% lambdaLasso = lambdaLasso/sqrt(n); % rescaling needed by glmnet
% options.lambda = lambdaLasso;
% options.type = 'naive';
% options = glmnetSet(options);
% % run lasso
% fitLasso = glmnet(X,y,'gaussian', options);
% estimateLasso = fitLasso.beta;
% hamLassoV(jjj)=sum(signbeta~= sign(estimateLasso));
% end
% %%% We pick the best lasso can do in the above 50 runs.
% hamLasso=min(hamLassoV)


%% univariate hard thresholding
yt=X'*y;
betaHard = tau.*sign(yt).*(abs(yt)>tau);
hamHard = sum(signbeta ~= sign(betaHard)) 

%% prepare for the Iterative graphlet screening
%yt=X'*y;
%%% Estiamte the gram matrix, and find all the small connected subgraphs
 gram = X'*X;
[gram, gramBias] = ThresholdGram(gram);
neighbor = sparse((gram ~= 0));
neighbor = neighbor-sparse(1:p, 1:p, 1);
%%hist(diag(neighbor))
nm = 3;
allCG = FindAllCG(neighbor, nm);


%% perform the iterative graphlet screening
betaRgs = IterGS(yt, gram, gramBias, allCG, sp, tau, nm);
hamRgs = sum(sign(betaRgs)~= signbeta)
sum(betaRgs ~= 0)


