%% Pattern level transition with regularization experiment
% 12/15/2019
% The purpose is to check if regularization works by shrinking some
% impossible pathways to 0
% Compare to the previous version, this version contains the following
% updates
% (1) save ebic from group penalty per replication
% (2) add regular bic as a model selection criterion, two versions of bic
% will be considered: counting the number of non-zero elements vs. number
% of non-zero columns
% (3) change lambda from group penalty to -0.4:0.1:-2
% (4) After choosing the best model using one of the BIC criteria, instead
% of re-running the model again, we first check if the same model has
% appeared before or not.
% (5) Delete the results from individual penalty to same time for now.


clear all
clc
N=3200;
K1=5;
K=10;
truea1=zeros(2^K1,K1);
a=[0;1];
for k=1:K1
    true1=repmat(a,2^(k-1),2^(K1-k));
    truea1(:,k)=reshape(true1',2^K1,1);
end
prior1=ones(1,2^K1).*1/(2^K1);

Replications=25;
results=zeros(Replications,36);

bic_group=zeros(17,6,Replications);
%%
for rep=1:1
    [rep,1]
    % True alpha at time 1
    true_alpha_t1=zeros(N,K1);
    true_alpha_t1(:,1:K1)=repmat(truea1,100,1);  % 525=3360/2^K1
    % True alpha at time 2
    % Set up learning rate
    learning=zeros(K,1);
    forget=zeros(K1,1);
    for k=1:K
        learning(k)=0.6;
    end
    for k=1:K1
        forget(k)=0.2;
    end
    true_alpha_t2=zeros(N,K);
    for i=1:N
        for k=1:K1
            if(true_alpha_t1(i,k)==1)
                true_alpha_t2(i,k)=binornd(1,1-forget(k));
            end
            if(true_alpha_t1(i,k)==0)
                true_alpha_t2(i,k)=binornd(1,learning(k));
            end
        end
        if (true_alpha_t2(i,1)==1)
            true_alpha_t2(i,6)=binornd(1,learning(6));
        end
        if (true_alpha_t2(i,1)*true_alpha_t2(i,2)==1)
            true_alpha_t2(i,7)=binornd(1,learning(7));
        end
        if (true_alpha_t2(i,3)*true_alpha_t2(i,4)==1)
            true_alpha_t2(i,8)=binornd(1,learning(8));
        end
        if (true_alpha_t2(i,3)==1)
            true_alpha_t2(i,9)=binornd(1,learning(9));
        end
        if (true_alpha_t2(i,4)*true_alpha_t2(i,5)==1)
            true_alpha_t2(i,10)=binornd(1,learning(10));
        end
    end
    
    truea2=zeros(2^K,K);  % all possible patterns
    for k=1:K
        true1=repmat(a,2^(k-1),2^(K-k));
        truea2(:,k)=reshape(true1',2^K,1);
    end
    transitionm=zeros(2^K1,2^K);
    
    for i=1:2^K1
        for j=1:2^K
            part1=1.0;
            for k=1:K1
                temp1=((1-forget(k))^truea1(i,k)*learning(k)^(1-truea1(i,k)))^truea2(j,k);
                temp2=(forget(k)^truea1(i,k)*(1-learning(k))^(1-truea1(i,k)))^(1-truea2(j,k));
                part1=part1*temp1*temp2;
            end
            part21=((learning(6)^truea2(j,1)*0^(1-truea2(j,1)))^truea2(j,6))*(((1-learning(6))^truea2(j,1)*1^(1-truea2(j,1)))^(1-truea2(j,6)));%alpha6
            part23=((learning(9)^truea2(j,3)*0^(1-truea2(j,3)))^truea2(j,9))*(((1-learning(9))^truea2(j,3)*1^(1-truea2(j,3)))^(1-truea2(j,9)));%alpha9
            part22=(learning(7)^(truea2(j,1)*truea2(j,2))*0^(2-truea2(j,1)-truea2(j,2)))^truea2(j,7)*((1-learning(7))^(truea2(j,1)*truea2(j,2)))^(1-truea2(j,7)); %alpha7
            part24=(learning(8)^(truea2(j,3)*truea2(j,4))*0^(2-truea2(j,3)-truea2(j,4)))^truea2(j,8)*((1-learning(8))^(truea2(j,3)*truea2(j,4)))^(1-truea2(j,8)); %alpha8
            part25=(learning(10)^(truea2(j,4)*truea2(j,5))*0^(2-truea2(j,4)-truea2(j,5)))^truea2(j,10)*((1-learning(10))^(truea2(j,4)*truea2(j,5)))^(1-truea2(j,10)); %alpha10
            transitionm(i,j)=part1*part21*part22*part23*part24*part25;
        end
    end
    
    sum(sum(transitionm==0)) % 27392 out of 32768
    sum(sum(transitionm<=1/(2*N))) %
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Create item parameters and responses
    % Create a test length of 30 at time 1, and 54 at time 2
    % create q-matrix for the first stage: only measure a1~a3
    % find q-vectors that measure 2 out of 5 attributes
    combo1=combnk(1:K1,2);
    leng=size(combo1);
    temp=zeros(leng(1), K1);
    for l=1:leng(1)
        temp(l, combo1(l,:))=1;
    end
    
    % find q-vectors that measure 3 out of 5 attributes
    combo1=combnk(1:K1,3);
    leng=size(combo1);
    temp2=zeros(leng(1), K1);
    for l=1:leng(1)
        temp2(l, combo1(l,:))=1;
    end
    qtest=[eye(K1);eye(K1);temp; temp2]; % q-matrix for the first 30 items
    qk1=[qtest,zeros(size(qtest,1),size(qtest,2))];
    
    % create another 24 items measure a1~a10
    trueq3=zeros(24,K);
    for j=1:8
        temp1=randsample(1:5,1);
        temp2=randsample(6:10,1);
        trueq3(j,temp1)=1;
        trueq3(j,temp2)=1;
    end
    for j=9:16
        temp1=randsample(1:5,2);
        temp2=randsample(6:10,1);
        trueq3(j,temp1)=1;
        trueq3(j,temp2)=1;
    end
    for j=17:24
        temp1=randsample(1:5,1);
        temp2=randsample(6:10,2);
        trueq3(j,temp1)=1;
        trueq3(j,temp2)=1;
    end
    qk2=[zeros(size(qtest,1),size(qtest,2)), qtest];
    qtest2=[qk2;trueq3];
    
    trueq=[qk1;qtest2];
    L=size(trueq,1);
    itempara=zeros(L,K+1);
    itempara(:,1)=unifrnd(0.01,0.2,[L,1]); % intercept
    for j=1:L
        tq=sum(trueq(j,:));
        index=find(trueq(j,:)==1);
        itempara(j,index+1)=(0.8-0.05*3)/tq+unifrnd(0.01,0.05,[1,tq]); % Note: I used K1 here because otherwise, the item paraemter will be bad.
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Create response matrix at time 1
    L1=30;
    truealpha=true_alpha_t1;
    response1=zeros(N,L1); % set the response matrix
    for i=1:N
        for j=1:L1
            cprob=itempara(j,1:(K1+1))*[1,truealpha(i,1:K1)]';
            temp=rand;
            if(cprob>=temp)
                response1(i,j)=1;
            end
        end
    end
    % create response matrix at time 2
    L2=54;
    truealpha=true_alpha_t2;
    response2=zeros(N,L2); % set the response matrix
    for i=1:N
        for j=(L1+1):(L1+L2)
            cprob=itempara(j,:)*[1,truealpha(i,:)]';
            temp=rand;
            if(cprob>=temp)
                response2(i,j-L1)=1;
            end
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % EM without penalty
    % create a matrix of likelihood
    % because item parameters and q-matrix are all known in this case, we do
    % not need to update likelihood during estimation
    %%%%%%%%%%%%%%%%%%%%%Very time-consuming likelihood computation%%%%%%%%%%%
    likelihood=ones(2^K1, 2^K,N);
    for i=1:N
        for l1=1:2^K1
            for l2=1:2^K
                tempprob=itempara(1:30,1:(K1+1))*[1,truea1(l1,1:K1)]';
                temp1=response1(i,:)*log(tempprob)+(1-response1(i,:))*log(1-tempprob);
                tempprob=itempara(31:84,:)*[1,truea2(l2,1:K)]';
                temp2=response2(i,:)*log(tempprob)+(1-response2(i,:))*log(1-tempprob);
                likelihood(l1,l2,i)=exp(temp1+temp2);
            end
        end
    end
    
    %
    %pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05];
    tic
    pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
    tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
    Number_of_EM_Cycles=0;
    flag1=1;
    Accuracy_of_Iteration=0.0001;
    Maximum_EM_Cycles=150;
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        % compute posterior
        eta=ones(2^K1, 2^K,N);
        denominator_eta=zeros(N,1);
        for i=1:N
            for l1=1:2^K1
                for l2=1:2^K
                    eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                    denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                end
            end
            eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
        end
        n_l1l2=sum(eta,3);
        n_l1=sum(n_l1l2,2);
        n_l2=sum(n_l1l2);
        
        pi1_final=(n_l1./N)';
        tau12_final=n_l1l2./n_l1;
        % this stopping criterion does not need to be updated even if some
        % tau12 contains 0
        if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
        end
    end
    toc
    % It seems that without penalty, the zero-elements will shrink to be very
    % small.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    results(rep,1)=sum(abs(pi1_final-prior1))/(2^K1);  %0.002
    results(rep,2)=sum(sum(abs(tau12_final-transitionm)));
    %sum(sum(abs(tau12_final(transitionm>0)-transitionm(transitionm>0))))/27392 %0.0013
    tau12_0=(tau12_final<=1/2/N).*1.0;
    transition0=(transitionm==0);
    transition0=transition0.*1.0;
    tau12_0=(tau12_final<=1/2/N).*1.0;
    results(rep,3)=sum(sum(tau12_0))/2^15;
    results(rep,4)=sum(sum(abs(tau12_0-transition0)))/(2^15);
    true0=find(transitionm==0);
    estimate0=find(tau12_final<=1/2/N);
    results(rep,5)=length(intersect(true0,estimate0))/27392; % proportion of true 0's that are correctly recovered
    true1=find(transitionm>0);
    estimate1=find(tau12_final>1/2/N);
    results(rep,6)=length(intersect(true1,estimate1))/(32768-27392);
    %sum(sum(abs(tau12_final(transitionm>0.01)-transitionm(transitionm>0.01))))/length(find(transitionm>0.01)) %0.0157 seems large?
    

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Group penalty: Regularized EM
    pgamma=-0.8:-0.05:-1.8; % smaller pgamma induces more sparsity
    C=0.0001;
    rhoN=1/(50*N); % a threshold according to Gu & Xu
    Accuracy_of_Iteration=0.0001;
    ebic=zeros(length(pgamma),6);
    interimcheck=zeros(length(pgamma),2);
    for penalty=20:length(pgamma)
        Number_of_EM_Cycles=0;
        %%%%%%%Setting initial values
        pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
        tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
        flag1=1;
        tau12_final=zeros(2.^K1,2.^K);
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            % compute posterior
            % Add a group level check: if the column sum of all elements = 0
            nzerocolumn=ones(1,2^K);
            for l2=1:2^K
                if (sum(tau12_initial(:,l2))<=rhoN*2^K1)
                    tau12_initial(:,l2)=0;
                    nzerocolumn(l2)=0;  % record which columns are 0's
                end
            end
            n_l1p=sum(tau12_initial,2);
            tau12_initial=tau12_initial./n_l1p; % standardize it so that row sums are 1
            eta=ones(2^K1, 2^K,N);
            denominator_eta=zeros(N,1);
            for i=1:N
                for l1=1:2^K1
                    for l2=1:2^K
                        eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                        denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                    end
                end
                eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
            end
            
            n_l1l2=sum(eta,3);
            n_l1=sum(n_l1l2,2);% row sum
            n_l2=sum(n_l1l2); % column sum
            
            pi1_final=(n_l1./N)';
            
            temp1=find(nzerocolumn==1); % index non-zero columns
            weight=tau12_initial(:,temp1)./sum(tau12_initial(:,temp1));
            tau12_final(:,temp1)=max(C,(pgamma(penalty).*weight+n_l1l2(:,temp1)));
            
            temp2=find(nzerocolumn==0); % index zero columns
            tau12_final(:,temp2)=0;
            
            n_l1p=sum(tau12_final,2);
            tau12_final=tau12_final./n_l1p;
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
            end
            % Note that different from JASA code, there is no step of "finding active set"
        end
        %    tau12_final_p=max(rhoN.*ones(2^K1,2^K),tau12_final);
        %    tau12_final_p(tau12_final_p==rhoN)=0;
        interimcheck(penalty,1)=sum(sum(tau12_final==0));
        interimcheck(penalty,2)=Number_of_EM_Cycles;
        
        
        % After deciding on the sparsity structure of tau12, re-run EM to
        % obtain EBIC
        pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
        tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
        tau12_initial(tau12_final==0)=0;
        tempp=sum(tau12_initial,2);
        tau12_initial=tau12_initial./tempp; % standardization such as the row sum is still 0
        % Once the initials contain elements of 0, those 0's will be carried
        % forward in the following iterations?
        Number_of_EM_Cycles=0;
        flag1=1;
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            % compute posterior
            eta=ones(2^K1, 2^K,N);
            denominator_eta=zeros(N,1);
            for i=1:N
                for l1=1:2^K1
                    for l2=1:2^K
                        eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                        denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                    end
                end
                eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
            end
            n_l1l2=sum(eta,3);
            n_l1=sum(n_l1l2,2);
            n_l2=sum(n_l1l2);
            
            pi1_final=(n_l1./N)';
            tau12_final=n_l1l2./n_l1;
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
            end
        end
        % Compute EBIC
        marginal_log=0.0;
        for i=1:N
            marginal_log=marginal_log+log(sum(sum(pi1_final'.*tau12_final.*likelihood(:,:,i))));
        end
        count1=sum(sum(tau12_final>0));
        %temp=log(nchoosek(2^K, count1)); %maybe a very large value
        count_column=2^K-sum(sum(tau12_final)==0);
        temp1=log(nchoosek(2^K, count_column));
        % ebic(penalty,1)=-2*marginal_log+count_column*log(N)+2*temp1;
 %       ebic(penalty,1)=-2*marginal_log+count1*log(log(N))*log(N); %+2*temp1; %model1
        ebic(penalty,1)=-2*marginal_log+count_column*log(log(N))*log(N);
        ebic(penalty,2)=-2*marginal_log;
        ebic(penalty,3)=count1;
        ebic(penalty,4)=-2*marginal_log+count_column*log(N)+2*temp1; % model2
        ebic(penalty,5)=-2*marginal_log+count1*log(N) % model 3
        ebic(penalty,6)=-2*marginal_log+count_column*log(N) % model4
        [penalty, Number_of_EM_Cycles]
    end
    bic_group(:,:, rep)=ebic;
    
    best_model=zeros(4,1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Using ebic to find the best one, and re-run the regularized EM followed by EM
    penalty=max(find(ebic(:,1)==min(ebic(:,1))));
    best_model(1)=penalty;
    %%%%%%%%%%%%%%%%%%%%%%%%Below is the same as in previous chunk
    Number_of_EM_Cycles=0;
    %%%%%%%Setting initial values
    pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
    tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
    flag1=1;
    tau12_final=zeros(2.^K1,2.^K);
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        % compute posterior
        % Add a group level check: if the column sum of all elements = 0
        nzerocolumn=ones(1,2^K);
        for l2=1:2^K
            if (sum(tau12_initial(:,l2))<=rhoN*2^K1)
                tau12_initial(:,l2)=0;
                nzerocolumn(l2)=0;  % record which columns are 0's
            end
        end
        n_l1p=sum(tau12_initial,2);
        tau12_initial=tau12_initial./n_l1p; % standardize it so that row sums are 1
        eta=ones(2^K1, 2^K,N);
        denominator_eta=zeros(N,1);
        for i=1:N
            for l1=1:2^K1
                for l2=1:2^K
                    eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                    denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                end
            end
            eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
        end
        
        n_l1l2=sum(eta,3);
        n_l1=sum(n_l1l2,2);% row sum
        n_l2=sum(n_l1l2); % column sum
        
        pi1_final=(n_l1./N)';
        
        temp1=find(nzerocolumn==1); % index non-zero columns
        weight=tau12_initial(:,temp1)./sum(tau12_initial(:,temp1));
        tau12_final(:,temp1)=max(C,(pgamma(penalty).*weight+n_l1l2(:,temp1)));
        
        temp2=find(nzerocolumn==0); % index zero columns
        tau12_final(:,temp2)=0;
        
        n_l1p=sum(tau12_final,2);
        tau12_final=tau12_final./n_l1p;
        % this stopping criterion does not need to be updated even if some
        % tau12 contains 0
        if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % After deciding on the sparsity structure of tau12, re-run EM to
    % obtain final parameter estimates
    pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
    tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
    tau12_initial(tau12_final==0)=0;
    tempp=sum(tau12_initial,2);
    tau12_initial=tau12_initial./tempp; % standardization such as the row sum is still 0
    % Once the initials contain elements of 0, those 0's will be carried
    % forward in the following iterations?
    Number_of_EM_Cycles=0;
    flag1=1;
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        % compute posterior
        eta=ones(2^K1, 2^K,N);
        denominator_eta=zeros(N,1);
        for i=1:N
            for l1=1:2^K1
                for l2=1:2^K
                    eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                    denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                end
            end
            eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
        end
        n_l1l2=sum(eta,3);
        n_l1=sum(n_l1l2,2);
        n_l2=sum(n_l1l2);
        
        pi1_final=(n_l1./N)';
        tau12_final=n_l1l2./n_l1;
        % this stopping criterion does not need to be updated even if some
        % tau12 contains 0
        if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
        end
    end
    results(rep,13)=sum(abs(pi1_final-prior1))/(2^K1);  %0.002
    results(rep,14)=sum(sum(abs(tau12_final-transitionm)));
    %sum(sum(abs(tau12_final(transitionm>0)-transitionm(transitionm>0))))/27392 %0.0013
    results(rep,15)=sum(sum(tau12_final==0))/(2^15);
    tau12_0=(tau12_final==0).*1.0;
    transition0=(transitionm==0);
    transition0=transition0.*1.0;
    results(rep,16)=sum(sum(abs(tau12_0-transition0)))/(2^15);
    true0=find(transitionm==0);
    estimate0=find(tau12_final==0);
    results(rep,17)=length(intersect(true0,estimate0))/27392; % proportion of true 0's that are correctly recovered
    true1=find(transitionm>0);
    estimate1=find(tau12_final>0);
    results(rep,18)=length(intersect(true1,estimate1))/(32768-27392);
    
    
    

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Using ebic to find the best one, and re-run the regularized EM followed by EM
    penalty=max(find(ebic(:,6)==min(ebic(:,6))));
    best_model(4)=penalty;
    tempo=find(best_model(4)==best_model(1:3));
    if (length(tempo)>0)
        results(rep, 31:36)=results(rep, (tempo(1)*6+7):(tempo(1)*6+12));
    else
        Number_of_EM_Cycles=0;
        %%%%%%%Setting initial values
        pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
        tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
        flag1=1;
        tau12_final=zeros(2.^K1,2.^K);
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            % compute posterior
            % Add a group level check: if the column sum of all elements = 0
            nzerocolumn=ones(1,2^K);
            for l2=1:2^K
                if (sum(tau12_initial(:,l2))<=rhoN*2^K1)
                    tau12_initial(:,l2)=0;
                    nzerocolumn(l2)=0;  % record which columns are 0's
                end
            end
            n_l1p=sum(tau12_initial,2);
            tau12_initial=tau12_initial./n_l1p; % standardize it so that row sums are 1
            eta=ones(2^K1, 2^K,N);
            denominator_eta=zeros(N,1);
            for i=1:N
                for l1=1:2^K1
                    for l2=1:2^K
                        eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                        denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                    end
                end
                eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
            end
            
            n_l1l2=sum(eta,3);
            n_l1=sum(n_l1l2,2);% row sum
            n_l2=sum(n_l1l2); % column sum
            
            pi1_final=(n_l1./N)';
            
            temp1=find(nzerocolumn==1); % index non-zero columns
            weight=tau12_initial(:,temp1)./sum(tau12_initial(:,temp1));
            tau12_final(:,temp1)=max(C,(pgamma(penalty).*weight+n_l1l2(:,temp1)));
            
            temp2=find(nzerocolumn==0); % index zero columns
            tau12_final(:,temp2)=0;
            
            n_l1p=sum(tau12_final,2);
            tau12_final=tau12_final./n_l1p;
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % After deciding on the sparsity structure of tau12, re-run EM to
        % obtain final parameter estimates
        pi1_initial=[3/64.*ones(1,16), 1/64.*ones(1,16)];
        tau12_initial=repmat(0.00097656.*ones(1,2.^K),2.^K1,1);
        tau12_initial(tau12_final==0)=0;
        tempp=sum(tau12_initial,2);
        tau12_initial=tau12_initial./tempp; % standardization such as the row sum is still 0
        % Once the initials contain elements of 0, those 0's will be carried
        % forward in the following iterations?
        Number_of_EM_Cycles=0;
        flag1=1;
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            % compute posterior
            eta=ones(2^K1, 2^K,N);
            denominator_eta=zeros(N,1);
            for i=1:N
                for l1=1:2^K1
                    for l2=1:2^K
                        eta(l1,l2,i)=pi1_initial(l1)*tau12_initial(l1,l2)*likelihood(l1,l2,i);
                        denominator_eta(i)=denominator_eta(i)+eta(l1,l2,i);
                    end
                end
                eta(:,:,i)=eta(:,:,i)./denominator_eta(i);
            end
            n_l1l2=sum(eta,3);
            n_l1=sum(n_l1l2,2);
            n_l2=sum(n_l1l2);
            
            pi1_final=(n_l1./N)';
            tau12_final=n_l1l2./n_l1;
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            if ((max(max(abs(tau12_final-tau12_initial)),max(abs(pi1_final-pi1_initial))))<Accuracy_of_Iteration)
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
            end
        end
        results(rep,31)=sum(abs(pi1_final-prior1))/(2^K1);  %0.002
        results(rep,32)=sum(sum(abs(tau12_final-transitionm)));
        %sum(sum(abs(tau12_final(transitionm>0)-transitionm(transitionm>0))))/27392 %0.0013
        results(rep,33)=sum(sum(tau12_final==0))/(2^15);
        tau12_0=(tau12_final==0).*1.0;
        transition0=(transitionm==0);
        transition0=transition0.*1.0;
        results(rep,34)=sum(sum(abs(tau12_0-transition0)))/(2^15);
        true0=find(transitionm==0);
        estimate0=find(tau12_final==0);
        results(rep,35)=length(intersect(true0,estimate0))/27392; % proportion of true 0's that are correctly recovered
        true1=find(transitionm>0);
        estimate1=find(tau12_final>0);
        results(rep,36)=length(intersect(true1,estimate1))/(32768-27392);
        
    end

    %%%%%%%%%%%%%%%%%%%%%%%%Below is the same as in previous chunk
    filename = 'K=10_3200_low_learning_solution_path_1_6.mat';
    save(filename, 'results','bic_group','ebic','N','learning','forget')
    
    
end

