%% Pattern level transition with regularization
% Use DINA as measurement model
% Per reviewer's request, we now assume DINA model parameters are unknown
clear all
clc
N=800;
K=6;
K1=3;
truea1=zeros(2^K1,K1); % all possible patterns
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,30);

%%
for rep=1:Replications
    [rep,1]
    % True alpha at time 1
    true_alpha_t1=zeros(N,K1);
    true_alpha_t1(:,1:K1)=repmat(truea1,100,1);
    % True alpha at time 2
    % Set up learning rate
    learning=zeros(K,1);
    forget=zeros(K1,1);
    for k=1:K
        %learning(k)=unifrnd(0.6,0.7,1);
        learning(k)=0.8;
    end
    for k=1:K1
        %  forget(k)=unifrnd(0.1,0.15,1);
        forget(k)=0.15;
    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,4)=binornd(1,learning(4));
        end
        if (true_alpha_t2(i,1)*true_alpha_t2(i,2)==1)
            true_alpha_t2(i,5)=binornd(1,learning(5));
        end
        if (true_alpha_t2(i,3)==1)
            true_alpha_t2(i,6)=binornd(1,learning(6));
        end
    end
    
    
    % create a 8*64 transition matrix
    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(4)^truea2(j,1)*0^(1-truea2(j,1)))^truea2(j,4))*(((1-learning(4))^truea2(j,1)*1^(1-truea2(j,1)))^(1-truea2(j,4)));
            part23=((learning(6)^truea2(j,3)*0^(1-truea2(j,3)))^truea2(j,6))*(((1-learning(6))^truea2(j,3)*1^(1-truea2(j,3)))^(1-truea2(j,6)));
            part22=(learning(5)^(truea2(j,1)*truea2(j,2))*0^(2-truea2(j,1)-truea2(j,2)))^truea2(j,5)*((1-learning(5))^(truea2(j,1)*truea2(j,2)))^(1-truea2(j,5));
            transitionm(i,j)=part1*part21*part22*part23;
        end
    end
    transition0=(transitionm==0);
    transition0=transition0.*1.0;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % % 320 entries in transitionm=0;
    % sum(sum(transitionm==0)) %320
    % sum(sum(transitionm<=1/(2*N))) %320
    % % mininum non-zero transition is 0.001
    % The purpose is to recover the transition matrix (with proper 0's) from the two time point responses
    % Generate item parameters and responses
    % Create a test with length of 20 at time 1, and 40 at time 2.
    % Test length could be manipulated.
    % We want to have a long test length at first, so that measurement error is
    % reduced
    N=size(true_alpha_t1,1);
    % create item bank
    K1=3;  % 3 attributes measured in stage I
    % create q-matrix for the first stage: only measure a1~a3
    qtest=[eye(K1);eye(K1);1,1,0; 1,0,1;0,1,1; 1,1,1];
    qk1=[qtest,zeros(size(qtest,1),size(qtest,2))];
    trueq1=repmat(qk1,2,1); % the first 20 items only measure a1~a3
    qk2=[zeros(size(qtest,1),size(qtest,2)), qtest];
    trueq2=repmat(qk2,1,1); % the second 10 items only measure a4~a6
    % create another 30 items measure a1~a6
    trueq3=zeros(30,K);
    for j=1:10
        temp1=randsample(1:3,1);
        temp2=randsample(4:6,1);
        trueq3(j,temp1)=1;
        trueq3(j,temp2)=1;
    end
    for j=11:20
        temp1=randsample(1:3,2);
        temp2=randsample(4:6,1);
        trueq3(j,temp1)=1;
        trueq3(j,temp2)=1;
    end
    for j=21:30
        temp1=randsample(1:3,1);
        temp2=randsample(4:6,2);
        trueq3(j,temp1)=1;
        trueq3(j,temp2)=1;
    end
    trueq=[trueq1;trueq2;trueq3];
    L=size(trueq,1);
    itempara=zeros(L,2);
    itempara(:,1)=unifrnd(0.1,0.2,[L,1]); % guessing parameter
    itempara(:,2)=unifrnd(0.6,0.8,[L,1]); % slipping=1-itempara(j,1)-itempara(j,2)
    L1=20;
    truealpha=true_alpha_t1;
    response1=zeros(N,L1); % set the response matrix
    for i=1:N
        for j=1:L1
            multi=prod(truealpha(i,1:K1).^trueq(j,1:K1));
            cprob=itempara(j,:)*[1,multi]';
            temp=rand;
            if(cprob>=temp)
                response1(i,j)=1;
            end
        end
    end
    L2=40;
    truealpha=true_alpha_t2;
    response2=zeros(N,L2); % set the response matrix
    for i=1:N
        for j=(L1+1):(L1+L2)
            multi=prod(truealpha(i,1:K).^trueq(j,1:K));
            cprob=itempara(j,:)*[1,multi]';
            temp=rand;
            if(cprob>=temp)
                response2(i,j-L1)=1;
            end
        end
    end
    
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % EM without penalty (assume known learning trajectory, i.e.,
    % (we have to assume known impermissible latent classes at time 2)
    % Note: EM without penalty may not work because if there are
    % impermissible latent classes at time 2, rho's at time 2 cannot be
    % estimated because for impermissible latent classes, the denominator
    % in the rho equation is close to 0
    
    pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05];
    tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
    
   % pi1_initial=prior1;
   % tau12_initial=transitionm;
    
    % Use true item parameters as initial values%%%%%%%%%%
    for l1=1:2^K1
        for j=1:L1
            multi=prod(truea1(l1,1:K1).^trueq(j,1:K1));
            rho1_initial(j, l1)=itempara(j,:)*[1,multi]';
        end
    end
    rho1_true=rho1_initial;
    
    for l2=1:2^K
        for j=1:L2
            multi=prod(truea2(l2,1:K).^trueq(j+L1,1:K));
            rho2_initial(j, l2)=itempara(j+L1,:)*[1,multi]';
        end
    end
    rho2_true=rho2_initial;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % time 1 DINA model
    % rho1_initial=0.5.*ones(L1, 2.^K1)+unifrnd(0.05,0.2,[L1,2^K1]);
    rho1_final=zeros(L1,2^K1);
    % time 2 DINA model
    %rho2_initial=0.5.*ones(L2, 2.^K)+unifrnd(0.05,0.2,[L2,2^K]);
    rho2_final=zeros(L2,2^K);
    
    % Directly estimate slipping and guessing parameters at both time
    % points
    slip_initial=ones(L,1).*0.1;
    guess_initial=ones(L,1).*0.11;
    
    
    
    
    Number_of_EM_Cycles=0;
    flag1=1;
    Accuracy_of_Iteration=0.0001;
    Maximum_EM_Cycles=300;
    
    %%%%%% Assume we know impermissible patterns%%%%%%
    impermissible=find(sum(transitionm)==0);
    tau12_initial(impermissible)=0;
    permissible=find(sum(transitionm)>0);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    itempara_initial=[guess_initial,1-slip_initial-guess_initial];
    % create two matrices based on trueq and truea1/truea2
    prodm1=zeros(L1,K1);
    prodm2=zeros(L2,K);
    
    for j=1:L1
        for l1=1:2^K1
            prodm1(j,l1)=prod(truea1(l1,1:K1).^trueq(j,1:K1));
        end
    end
    
    for j=1:L2
        for l2=1:2^K
            prodm2(j,l2)=prod(truea2(l2,1:K).^trueq(j+L1,1:K));
        end
    end
    
    slip_est=zeros(L,1);
    guess_est=zeros(L,1);
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        likelihood=ones(2^K1, 2^K,N);
        %for i=1:N
        for l1=1:2^K1
            for l2=1:2^K
                %tempprob=itempara(1:20,1:(K1+1))*[1,truea1(l1,1:K1)]';
                %tempprob=rho1_initial(:,l1);
                tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                
                temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                %tempprob=itempara(21:60,:)*[1,truea2(l2,1:K)]';
                %tempprob=rho2_initial(:,l2);
                tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                %if(sum(tempprob)>0)
                temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                %else
                %    temp2=response2*log(rho2_true(:,l2))+(1-response2)*log(1-rho2_true(:,l2));
                % I think using temp2=0 should not change the result
                % though because of the zero elements in tau12
                %end
                likelihood(l1,l2,:)=exp(temp1+temp2);
            end
        end
        %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;
        
        %        pi1_final=prior1;
        %        tau12_final=transitionm;
        
        % update rho1 and rho2
        %         tempo=sum(eta,[2,3]);
        %         for l1=1:2^K1
        %              rho1_final(:,l1)=sum(squeeze(eta(l1,:,:))*response1,1)'./tempo(l1);
        %         end
        
        %         rho1_final=rho1_true;
        %
        %         tempo=sum(eta,[1,3]);
        %         for l2=1:2^K
        %             if(sum(transitionm(:,l2))>0)
        %             rho2_final(:,l2)=sum(squeeze(eta(:,l2,:))*response2,1)'./tempo(l2);
        %             else
        %                 rho2_final(:,l2)=rho2_true(:,l2);
        %             end
        %         end
        %  rho2_final=rho2_true;
        
        %Update slipping and guessing parameters from the DINA model
        %Time 1
        tempo=sum(eta,[2,3]); % expected number of persons in each latent class
        for j=1:L1
            master=find(prodm1(j,:)==1);
            non_master=find(prodm1(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,2));  %8*800
            tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j)=R_jl0/I_jl0;
        end
        
        %Time 2
        tempo=sum(eta,[1,3]);
        for j=1:L2
            master=find(prodm2(j,:)==1);
            non_master=find(prodm2(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,1));
            tempoj=tt*response2(:,j);
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j+L1)=R_jl0/I_jl0;
        end
        
        itempara_final=[guess_est, 1-slip_est-guess_est];
        
        
        
        % this stopping criterion does not need to be updated even if some
        % tau12 contains 0
        dif1=max(max(abs(tau12_final-tau12_initial)));
        %dif2=max(max(abs(rho1_final-rho1_initial)));
        %dif3=max(max(abs(rho2_final(:, permissible)-rho2_initial(:, permissible))));
        dif2=max(max(abs(itempara_final-itempara_initial)));
        dif4=max(abs(pi1_final-pi1_initial));
        if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
            % rho1_initial=rho1_final;
            % rho2_initial=rho2_final;
            itempara_initial=itempara_final;
        end
    end
    results(rep,1)=sum(abs(pi1_final-prior1))/(2^K1); %0.0035
    results(rep,2)=sum(sum(abs(tau12_final-transitionm))); %4.9218
    results(rep,3)=sum(sum(tau12_final==0))/512;
    results(rep,4)=sum(sum(tau12_final<=1/15/N))/512;
    tau12_0=(tau12_final<=1/15/N).*1.0;
    results(rep,5)=sum(sum(abs(tau12_0-transition0)))/512; % in general, the proportion of sparsity correctly identified
    true0=find(transitionm==0);
    estimate0=find(tau12_final<=1/15/N);
    results(rep,6)=length(intersect(true0,estimate0))/320; % proportion of true 0's that are correctly recovered
    true1=find(transitionm>0);
    estimate1=find(tau12_final>1/15/N);
    results(rep,7)=length(intersect(true1,estimate1))/(512-320);% proportion of true 1's that are correctly recovered (we prefer not over shrinkage)
    results(rep,8)=sum(sum(abs(itempara-itempara_final)))/(2*L);
    
    
    % Penalized EM (group level)
    % the targe parameters are pi1, and tau12, where pi1 is a 1-by-2^K1 vector,
    % and tau12 is a 2^K1-by-2^K matrix
    % true values are stored in "prior1" and "transitionm"
    %%%%%%%Starting EM algorithm
    [rep,2]
    pgamma=0:-0.2:-2.4; % smaller pgamma induces more sparsity
    C=0.0001;
    rhoN=1/(15*N); % a threshold according to Gu & Xu
    Accuracy_of_Iteration=0.0001;
    ebic=zeros(length(pgamma),7);
    interimcheck=zeros(length(pgamma),2);
    for penalty=1:length(pgamma)
        Number_of_EM_Cycles=0;
        %%%%%%%Setting initial values
        pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05]; % I intentionally picked values that are different from truth
        tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
        flag1=1;
        tau12_final=zeros(2.^K1,2.^K);
        slip_initial=ones(L,1).*0.1;
        guess_initial=ones(L,1).*0.11;
        
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            %compute likelihood
            likelihood=ones(2^K1, 2^K,N);
            for l1=1:2^K1
                for l2=1:2^K
                    tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                    temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                    tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                    temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                    likelihood(l1,l2,:)=exp(temp1+temp2);
                end
            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;
            
            
            %Update slipping and guessing parameters from the DINA model
            %Time 1
            tempo=sum(eta,[2,3]); % expected number of persons in each latent class
            for j=1:L1
                master=find(prodm1(j,:)==1);
                non_master=find(prodm1(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,2));  %8*800
                tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j)=R_jl0/I_jl0;
            end
            
            %Time 2
            tempo=sum(eta,[1,3]);
            for j=1:L2
                master=find(prodm2(j,:)==1);
                non_master=find(prodm2(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,1));
                tempoj=tt*response2(:,j);
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j+L1)=R_jl0/I_jl0;
            end
            
            itempara_final=[guess_est, 1-slip_est-guess_est];
            
            
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            
            dif1=max(max(abs(tau12_final-tau12_initial)));
            dif2=max(max(abs(itempara_final-itempara_initial)));
            dif4=max(abs(pi1_final-pi1_initial));
            if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
                
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
                itempara_initial=itempara_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=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05];
        tau12_initial=repmat(0.0156.*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;
        slip_initial=ones(L,1).*0.1;
        guess_initial=ones(L,1).*0.11;
        
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            
            likelihood=ones(2^K1, 2^K,N);
            %for i=1:N
            for l1=1:2^K1
                for l2=1:2^K
                    tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                    temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                    tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                    temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                    likelihood(l1,l2,:)=exp(temp1+temp2);
                end
            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;
            
            
            %Update slipping and guessing parameters from the DINA model
            %Time 1
            tempo=sum(eta,[2,3]); % expected number of persons in each latent class
            for j=1:L1
                master=find(prodm1(j,:)==1);
                non_master=find(prodm1(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,2));  %8*800
                tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j)=R_jl0/I_jl0;
            end
            
            %Time 2
            tempo=sum(eta,[1,3]);
            for j=1:L2
                master=find(prodm2(j,:)==1);
                non_master=find(prodm2(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,1));
                tempoj=tt*response2(:,j);
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j+L1)=R_jl0/I_jl0;
            end
            
            itempara_final=[guess_est, 1-slip_est-guess_est];
            
            
            
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            dif1=max(max(abs(tau12_final-tau12_initial)));
            dif2=max(max(abs(itempara_final-itempara_initial)));
            dif4=max(abs(pi1_final-pi1_initial));
            if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
                
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
                itempara_initial=itempara_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));
        count_column=2^K-sum(sum(tau12_final)==0);
        temp=log(nchoosek(2^K1*2^K, count1)); %maybe a very large value
        ebic(penalty,1)=-2*marginal_log+count1*log(N)+2*temp;
        ebic(penalty,2)=-2*marginal_log;
        ebic(penalty,3)=count1;
        ebic(penalty,4)=-2*marginal_log+count1*log(log(N))*log(N);%GIC2
        ebic(penalty,5)=-2*marginal_log+count_column*log(N);%BIC using columns
        ebic(penalty,6)=-2*marginal_log+count1*log(N);%BIC
        ebic(penalty,7)=-2*marginal_log+count_column*log(log(N))*log(N);%GIC2 using columns
        
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Using ebic to find the best one, and re-run the regularized EM followed by EM
    penalty=max(find(ebic(:,1)==min(ebic(:,1))));
    %%%%%%%%%%%%%%%%%%%%%%%%Below is the same as in previous chunk
    Number_of_EM_Cycles=0;
    %%%%%%%Setting initial values
    pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05]; % I intentionally picked values that are different from truth
    tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
    flag1=1;
    tau12_final=zeros(2.^K1,2.^K);
    slip_initial=ones(L,1).*0.1;
    guess_initial=ones(L,1).*0.11;
    
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        likelihood=ones(2^K1, 2^K,N);
        %for i=1:N
        for l1=1:2^K1
            for l2=1:2^K
                tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                likelihood(l1,l2,:)=exp(temp1+temp2);
            end
        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;
        
        
        %Update slipping and guessing parameters from the DINA model
        %Time 1
        tempo=sum(eta,[2,3]); % expected number of persons in each latent class
        for j=1:L1
            master=find(prodm1(j,:)==1);
            non_master=find(prodm1(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,2));  %8*800
            tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j)=R_jl0/I_jl0;
        end
        
        %Time 2
        tempo=sum(eta,[1,3]);
        for j=1:L2
            master=find(prodm2(j,:)==1);
            non_master=find(prodm2(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,1));
            tempoj=tt*response2(:,j);
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j+L1)=R_jl0/I_jl0;
        end
        
        itempara_final=[guess_est, 1-slip_est-guess_est];
        
        
        
        
        
        dif1=max(max(abs(tau12_final-tau12_initial)));
        dif2=max(max(abs(itempara_final-itempara_initial)));
        dif4=max(abs(pi1_final-pi1_initial));
        if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
            itempara_initial=itempara_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 final parameter estimates
    pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05];
    tau12_initial=repmat(0.0156.*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;
    slip_initial=ones(L,1).*0.1;
    guess_initial=ones(L,1).*0.11;
    
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        likelihood=ones(2^K1, 2^K,N);
        %for i=1:N
        for l1=1:2^K1
            for l2=1:2^K
                tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                likelihood(l1,l2,:)=exp(temp1+temp2);
            end
        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;
        
        
        %Update slipping and guessing parameters from the DINA model
        %Time 1
        tempo=sum(eta,[2,3]); % expected number of persons in each latent class
        for j=1:L1
            master=find(prodm1(j,:)==1);
            non_master=find(prodm1(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,2));  %8*800
            tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j)=R_jl0/I_jl0;
        end
        
        %Time 2
        tempo=sum(eta,[1,3]);
        for j=1:L2
            master=find(prodm2(j,:)==1);
            non_master=find(prodm2(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,1));
            tempoj=tt*response2(:,j);
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j+L1)=R_jl0/I_jl0;
        end
        
        itempara_final=[guess_est, 1-slip_est-guess_est];
        
        
        dif1=max(max(abs(tau12_final-tau12_initial)));
        dif2=max(max(abs(itempara_final-itempara_initial)));
        dif4=max(abs(pi1_final-pi1_initial));
        if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
            itempara_initial=itempara_final;
        end
    end
    true0=find(transitionm==0);
    estimate0=find(tau12_final==0);
    results(rep,9)=sum(abs(pi1_final-prior1))/(2^K1);
    results(rep,10)=sum(sum(abs(tau12_final-transitionm)));
    results(rep,11)=sum(sum(tau12_final==0))/512;
    tau12_0=(tau12_final==0).*1.0;
    results(rep,12)=sum(sum(abs(tau12_0-transition0)))/512; % in general, the proportion of sparsity correctly identified
    true0=find(transitionm==0);
    estimate0=find(tau12_final==0);
    results(rep,13)=length(intersect(true0,estimate0))/320; % proportion of true 0's that are correctly recovered
    true1=find(transitionm>0);
    estimate1=find(tau12_final>0);
    results(rep,14)=length(intersect(true1,estimate1))/(512-320);% proportion of true 0's that are correctly recovered (we prefer not over shrinkage)
    results(rep,15)=sum(sum(abs(itempara-itempara_final)))/(2*L);
    
    % Penalized EM (individual level)
    [rep,3]
    pgamma=0:-0.05:-1.0; % 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),3);
    interimcheck=zeros(length(pgamma),2);
    for penalty=1:length(pgamma)
        Number_of_EM_Cycles=0;
        %%%%%%%Setting initial values
        pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05]; % I intentionally picked values that are different from truth
        tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
        flag1=1;
        slip_initial=ones(L,1).*0.1;
        guess_initial=ones(L,1).*0.11;
        
        while (flag1==1)
            Number_of_EM_Cycles=Number_of_EM_Cycles+1;
            
            if (Number_of_EM_Cycles>Maximum_EM_Cycles)
                break;
            end
            
            likelihood=ones(2^K1, 2^K,N);
            %for i=1:N
            for l1=1:2^K1
                for l2=1:2^K
                    tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                    temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                    tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                    temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                    likelihood(l1,l2,:)=exp(temp1+temp2);
                end
            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=max(C,(pgamma(penalty)+n_l1l2)); %./(2.^K*pgamma(penalty)+n_l1);
            n_l1p=sum(tau12_final,2);
            tau12_final=tau12_final./n_l1p;
            
            
            %Update slipping and guessing parameters from the DINA model
            %Time 1
            tempo=sum(eta,[2,3]); % expected number of persons in each latent class
            for j=1:L1
                master=find(prodm1(j,:)==1);
                non_master=find(prodm1(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,2));  %8*800
                tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j)=R_jl0/I_jl0;
            end
            
            %Time 2
            tempo=sum(eta,[1,3]);
            for j=1:L2
                master=find(prodm2(j,:)==1);
                non_master=find(prodm2(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,1));
                tempoj=tt*response2(:,j);
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j+L1)=R_jl0/I_jl0;
            end
            
            itempara_final=[guess_est, 1-slip_est-guess_est];
            
            
            
            dif1=max(max(abs(tau12_final-tau12_initial)));
            dif2=max(max(abs(itempara_final-itempara_initial)));
            dif4=max(abs(pi1_final-pi1_initial));
            if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
                itempara_initial=itempara_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_p==0));
        interimcheck(penalty,2)=Number_of_EM_Cycles;
        
        % After deciding on the sparsity structure of tau12, re-run EM to
        % obtain EBIC
        pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05];
        tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
        tau12_initial(tau12_final_p==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
            
            
            likelihood=ones(2^K1, 2^K,N);
            %for i=1:N
            for l1=1:2^K1
                for l2=1:2^K
                    tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                    temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                    tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                    temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                    likelihood(l1,l2,:)=exp(temp1+temp2);
                end
            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;
            
            %Update slipping and guessing parameters from the DINA model
            %Time 1
            tempo=sum(eta,[2,3]); % expected number of persons in each latent class
            for j=1:L1
                master=find(prodm1(j,:)==1);
                non_master=find(prodm1(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,2));  %8*800
                tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j)=R_jl0/I_jl0;
            end
            
            %Time 2
            tempo=sum(eta,[1,3]);
            for j=1:L2
                master=find(prodm2(j,:)==1);
                non_master=find(prodm2(j,:)==0);
                I_jl0=sum(tempo(non_master));
                I_jl1=sum(tempo(master));
                tt=squeeze(sum(eta,1));
                tempoj=tt*response2(:,j);
                R_jl0=sum(tempoj(non_master));
                R_jl1=sum(tempoj(master));
                slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
                guess_est(j+L1)=R_jl0/I_jl0;
            end
            
            itempara_final=[guess_est, 1-slip_est-guess_est];
            
            
            
            % this stopping criterion does not need to be updated even if some
            % tau12 contains 0
            dif1=max(max(abs(tau12_final-tau12_initial)));
            dif2=max(max(abs(itempara_final-itempara_initial)));
            dif4=max(abs(pi1_final-pi1_initial));
            if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
                flag1=0;
            else
                tau12_initial=tau12_final;
                pi1_initial=pi1_final;
                itempara_initial=itempara_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^K1*2^K, count1)); %maybe a very large value
        ebic(penalty,1)=-2*marginal_log+count1*log(N);%+2*temp;
        ebic(penalty,2)=-2*marginal_log;
        ebic(penalty,3)=count1;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Using ebic to find the best one, and re-run the regularized EM followed by EM
    penalty=max(find(ebic(:,1)==min(ebic(:,1))));
    %%%%%%%%%%%%%%%%%%%%%%%%Below is the same as in previous chunk
    Number_of_EM_Cycles=0;
    %%%%%%%Setting initial values
    pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05]; % I intentionally picked values that are different from truth
    tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
    flag1=1;
    slip_initial=ones(L,1).*0.1;
    guess_initial=ones(L,1).*0.11;
    
    while (flag1==1)
        Number_of_EM_Cycles=Number_of_EM_Cycles+1;
        
        if (Number_of_EM_Cycles>Maximum_EM_Cycles)
            break;
        end
        
        likelihood=ones(2^K1, 2^K,N);
        %for i=1:N
        for l1=1:2^K1
            for l2=1:2^K
                tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                likelihood(l1,l2,:)=exp(temp1+temp2);
            end
        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=max(C,(pgamma(penalty)+n_l1l2)); %./(2.^K*pgamma(penalty)+n_l1);
        n_l1p=sum(tau12_final,2);
        tau12_final=tau12_final./n_l1p;
        
        %Update slipping and guessing parameters from the DINA model
        %Time 1
        tempo=sum(eta,[2,3]); % expected number of persons in each latent class
        for j=1:L1
            master=find(prodm1(j,:)==1);
            non_master=find(prodm1(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,2));  %8*800
            tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j)=R_jl0/I_jl0;
        end
        
        %Time 2
        tempo=sum(eta,[1,3]);
        for j=1:L2
            master=find(prodm2(j,:)==1);
            non_master=find(prodm2(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,1));
            tempoj=tt*response2(:,j);
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j+L1)=R_jl0/I_jl0;
        end
        
        itempara_final=[guess_est, 1-slip_est-guess_est];
        
        
        dif1=max(max(abs(tau12_final-tau12_initial)));
        dif2=max(max(abs(itempara_final-itempara_initial)));
        dif4=max(abs(pi1_final-pi1_initial));
        if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
            itempara_initial=itempara_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_p==0));
    interimcheck(penalty,2)=Number_of_EM_Cycles;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % After deciding on the sparsity structure of tau12, re-run EM to
    % obtain final parameter estimate
    pi1_initial=[0.15,0.15,0.15,0.15,0.05,0.05,0.05,0.05];
    tau12_initial=repmat(0.0156.*ones(1,2.^K),2.^K1,1);
    tau12_initial(tau12_final_p==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
        likelihood=ones(2^K1, 2^K,N);
        %for i=1:N
        for l1=1:2^K1
            for l2=1:2^K
                tempprob=sum(itempara_initial(1:L1,:).*[ones(L1,1),prodm1(:,l1)],2);
                temp1=response1*log(tempprob)+(1-response1)*log(1-tempprob);
                tempprob=sum(itempara_initial((L1+1):L,:).*[ones(L2,1),prodm2(:,l2)],2);
                temp2=response2*log(tempprob)+(1-response2)*log(1-tempprob);
                likelihood(l1,l2,:)=exp(temp1+temp2);
            end
        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;
        
        %Update slipping and guessing parameters from the DINA model
        %Time 1
        tempo=sum(eta,[2,3]); % expected number of persons in each latent class
        for j=1:L1
            master=find(prodm1(j,:)==1);
            non_master=find(prodm1(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,2));  %8*800
            tempoj=tt*response1(:,j); % expected number of persons in each latent class who answer the item right
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j)=R_jl0/I_jl0;
        end
        
        %Time 2
        tempo=sum(eta,[1,3]);
        for j=1:L2
            master=find(prodm2(j,:)==1);
            non_master=find(prodm2(j,:)==0);
            I_jl0=sum(tempo(non_master));
            I_jl1=sum(tempo(master));
            tt=squeeze(sum(eta,1));
            tempoj=tt*response2(:,j);
            R_jl0=sum(tempoj(non_master));
            R_jl1=sum(tempoj(master));
            slip_est(j+L1)=(I_jl1-R_jl1)/I_jl1;
            guess_est(j+L1)=R_jl0/I_jl0;
        end
        
        itempara_final=[guess_est, 1-slip_est-guess_est];
        
        dif1=max(max(abs(tau12_final-tau12_initial)));
        dif2=max(max(abs(itempara_final-itempara_initial)));
        dif4=max(abs(pi1_final-pi1_initial));
        if(max([dif1,dif2, dif4])<Accuracy_of_Iteration)
            flag1=0;
        else
            tau12_initial=tau12_final;
            pi1_initial=pi1_final;
            itempara_initial=itempara_final;
        end
    end
    true0=find(transitionm==0);
    tau12_final_p=max(rhoN.*ones(2^K1,2^K),tau12_final);
    tau12_final_p(tau12_final_p==rhoN)=0;
    estimate0=find(tau12_final==0);
    
    
    results(rep,20)=sum(abs(pi1_final-prior1))/(2^K1);
    results(rep,21)=sum(sum(abs(tau12_final_p-transitionm)));
    results(rep,22)=sum(sum(tau12_final_p==0))/512;
    tau12_0=(tau12_final_p==0).*1.0;
    results(rep,23)=sum(sum(abs(tau12_0-transition0)))/512; % in general, the proportion of sparsity correctly identified
    true0=find(transitionm==0);
    estimate0=find(tau12_final_p==0);
    results(rep,24)=length(intersect(true0,estimate0))/320; % proportion of true 0's that are correctly recovered
    true1=find(transitionm>0);
    estimate1=find(tau12_final_p>0);
    results(rep,25)=length(intersect(true1,estimate1))/(512-320);% proportion of true 0's that are correctly recovered (we prefer not over shrinkage)
    results(rep,26)=sum(sum(abs(itempara-itempara_final)))/(2*L);
    filename = 'K=6_800_DINA_high_learning_revised.mat';
    save(filename)
end