


s=5; a=3; sigStrength=1; rho=0.9; n=100; 
noise=1;
R=500;
var=1;
L2=0.001;
alpha=0.01;
isInSet2=zeros(R,1);
isInSet21=zeros(R,1);
isInSet22=zeros(R,1);
proportionFound2=zeros(R,1);
proportionFound21=zeros(R,1);
proportionFound22=zeros(R,1);
sizeSet2=zeros(R,1);
sizeSet21=zeros(R,1);
sizeSet22=zeros(R,1);
fNoReject21=zeros(R,1);
%fNoReject22=zeros(R,1);
numberFalseModelsGoodFitF21=zeros(R,1);
%numberFalseModelsGoodFitF22=zeros(R,1);
d=1000; % total number of variables
cov1=rho.*ones(a+s)+(1-rho).*eye(a+s);
covMatrixInit=[cov1, zeros(s+a,d-(s+a));zeros(d-(s+a),s+a), eye(d-(s+a))];
covMatrix=diag(sqrt(var).*ones(d,1))*covMatrixInit*diag(sqrt(var).*ones(d,1));
intercept=5;
trueBetaInit=[sigStrength.*ones(s,1);zeros(d-s,1)];
seed=11;
rng(seed)
for r=1:R % replication index
    r
    permuteVec=randperm(d);
    trueBeta=trueBetaInit(permuteVec); % permute the rows of the initial beta vector
    I=eye(d);
    permMatrix=I(permuteVec,:);
    covPerm=permMatrix*covMatrix*(inv(permMatrix)); % permute rows and columns of the covariance matrix accordingly.
    XAll=mvnrnd(zeros(d,1),covPerm,n);
    activeSet=find(trueBeta);
    epsilon=sqrt(noise).*randn(n,1);
    YAll=intercept.*ones(n,1)+XAll*trueBeta+epsilon;
    nHalf=7*n/10;
    X=XAll(1:nHalf,:);
    Y=YAll(1:nHalf,1);
    setSelected2in3= phase1OnlySelect2(Y,X);
    isInSet2(r,1)=(all(ismember(activeSet,setSelected2in3))==1);
    proportionFound2(r,1)=sum(ismember(activeSet,setSelected2in3))/s;
    sizeSet2(r,1)=size(setSelected2in3,1);
    [setSelected22,setSelected21]=phase2Only(Y,X,setSelected2in3,L2);
    isInSet22(r,1)=(all(ismember(activeSet,setSelected22))==1);
    proportionFound22(r,1)=sum(ismember(activeSet,setSelected22))/s;
    sizeSet22(r,1)=size(setSelected22,1);
    %[fNoRejectSim,numberFalseModelsGoodFitFSim]=fTestModelsv2(YAll(nHalf+1:end,1),XAll(nHalf+1:end,:),setSelected22,alpha,activeSet);
    %fNoReject22(r,1)=fNoRejectSim;
    %numberFalseModelsGoodFitF22(r,1)=numberFalseModelsGoodFitFSim;
    isInSet21(r,1)=(all(ismember(activeSet,setSelected21))==1);
    proportionFound21(r,1)=sum(ismember(activeSet,setSelected21))/s;
    sizeSet21(r,1)=size(setSelected21,1);
    [fNoRejectSim,numberFalseModelsGoodFitFSim]=fTestModelsv2(YAll(nHalf+1:end,1),XAll(nHalf+1:end,:),setSelected21,alpha,activeSet);
    fNoReject21(r,1)=fNoRejectSim;
    numberFalseModelsGoodFitF21(r,1)=numberFalseModelsGoodFitFSim;
end% r


%mClock=ceil(clock);
%tStamp=1e6*mClock(2)+1e4*mClock(3)+100*mClock(4)+mClock(5);
fnme=['simuAug24Splitv3_n' num2str(n) '_s' num2str(s) '_sigToNoise' num2str(sigStrength/noise) '_a' num2str(a) '_rho' num2str(rho) '.mat'];
save(fnme); % and save the workspace





mean(isInSet21)
mean(fNoReject21)
mean(numberFalseModelsGoodFitF21)




mean(isInSet2)
mean(isInSet21)
mean(isInSetLasso)
mean(proportionFound2)
mean(proportionFound21)
mean(proportionFoundLasso)
mean(sizeSet2)
mean(sizeSet21)
mean(sizeSetLasso)

std(isInSet2)
std(isInSet21)
std(isInSetLasso)
std(proportionFound2)
std(proportionFound21)
std(proportionFoundLasso)
std(sizeSet2)
std(sizeSet21)
std(sizeSetLasso)



