% Andrea Trevino
% 11/19/13

% Plotting the confusions for RD vs RC children
clear

load DataNSCM_Alina.mat
load DataNSCM_Angela.mat
load DataNSCM_Anton.mat
load DataNSCM_Bob.mat
load DataNSCM_Carly.mat
load DataNSCM_Edward.mat
load DataNSCM_Evan.mat
load DataNSCM_Joanna.mat
load DataNSCM_Latisha.mat
load DataNSCM_Laura.mat
load DataNSCM_Miguel.mat
load DataNSCM_Norene.mat
load DataNSCM_Shauna.mat
load DataNSCM_Teddy.mat
load DataNSCM_Tony.mat

%% Compile the data

names={'Alina';'Angela';'Anton';'Bob';'Carly';'Edward';'Evan';'Joanna';'Latisha';'Laura';'Miguel';'Norene';'Shauna';'Teddy';'Tony'};
disabled={'Alina';'Angela';'Edward';'Latisha';'Laura';'Norene';'Shauna';'Teddy';'Tony'};
alldata={};

White=[210 210 210];
Gray=[110 110 110];
Colors=[White; Gray]/255;

% removed 'H','l','G'
cons_Dbet=['C';'S';'T';'f';'h';'k';'p';'s';'t';'D';'b';'d';'g';'J';'I';'m';'n';'r';'Z';'v';'w';'y';'z'];
vow_Dbet=['@';'A';'E';'I';'O';'R';'U';'W';'Y';'c';'a';'e';'i';'o';'u'];

total_CM=zeros(length(cons_Dbet)); % Add up all confusions to find MOST common confusions
total_CM_RC=zeros(length(cons_Dbet)); % reading control
total_CM_RD=zeros(length(cons_Dbet)); % reading disabled
% saved filled versions of above 3 matrices, loading result here:
load total_CM_norm % 'norm' means normalized
total_CM_norm_conf=total_CM_norm-(total_CM_norm.*eye(size(total_CM_norm)));
total_CM_RC_norm_conf=total_CM_RC_norm-(total_CM_RC_norm.*eye(size(total_CM_RC_norm)));
total_CM_RD_norm_conf=total_CM_RD_norm-(total_CM_RD_norm.*eye(size(total_CM_RD_norm)));
% the sorted vector of population confusions
population_CM_vec=sort(reshape(total_CM_norm_conf,length(cons_Dbet)^2,1),'descend');
RC_CM_vec=sort(reshape(total_CM_RC_norm_conf,length(cons_Dbet)^2,1),'descend');
RD_CM_vec=sort(reshape(total_CM_RD_norm_conf,length(cons_Dbet)^2,1),'descend');

figure(3)
plot(RC_CM_vec)
ylim([0 0.5])
xlim([-5 length(RC_CM_vec)])
set(gca,'FontSize',14)
ylabel('Probability of Error')
xlabel('Rank of Confusion')
title('Confusion Probabilities, Reading Control Group')
figure(4)
plot(RD_CM_vec)
ylim([0 0.5])
xlim([-5 length(RD_CM_vec)])
set(gca,'FontSize',14)
ylabel('Probability of Error')
xlabel('Rank of Confusion')
title('Confusion Probabilities, Reading Disabled Group')

% NO, find the most common reading disabled, and most common reading normal
% confusions, plot both in reading disabled order to show difference.
% Iterate through the subjects:
for n=1:length(names)
    filename=char(strcat('DataNSCM_',names(n)))
    eval(char(strcat('thisdata=',filename,';')))
    
    stims=unique(thisdata(2:end,4)); % get the list of stimuli
    % get the consonants tested for a CM
    cons_tested=[];
    for temp=1:length(stims) % pull consonant out of all stimuli
        test=stims{temp};
        %if any(test(1)==cons_Dbet) % if any C of a CV
            %cons_tested=[cons_tested; test(1)];
        cons_tested=[cons_tested; cons_Dbet(test(1)==cons_Dbet | test(2)==cons_Dbet)];   
        %end
    end
    cons_tested=unique(cons_tested);
    
    raw_stims=thisdata(2:end,4);
    %raw_resp=thisdata(2:end,5:7);
    % process data, pull out cons, vowel, get vowel Pc, cons Pc
    thisdata{1,12}='Stim Cons';
    thisdata{1,13}='Stim Vowel';
    thisdata{1,14}='Cons Pc';
    thisdata{1,15}='Vowel Pc';
    thisdata{1,16}='CV or VC';
    CV_consCM=zeros(length(cons_tested));
    CV_bigCM=zeros(length(cons_Dbet));
    VC_consCM=zeros(length(cons_tested));
    VC_bigCM=zeros(length(cons_Dbet));
    for temp=2:length(thisdata(:,1))
        thisstim=thisdata{temp,4};
        thisresp1=[];
        thisresp2=[];
        
        if any(thisstim(1)==cons_Dbet) % then a valid CV
            % populate thisdata
            thisdata{temp,12}=thisstim(1);
            thisdata{temp,13}=thisstim(2);
            thisdata{temp,16}='CV';
            
            flag=0; % remove and fix the case where the number of phones>2 later (never got implemented)
            for temp2=1:3
                for_parse=thisdata{temp,4+temp2};
                if any(~isnan(for_parse))&(length(for_parse)==1) % only one phone back
                    if any(for_parse==cons_Dbet)
                        thisresp1=[thisresp1;for_parse];
                        thisresp2=[thisresp2;'?']; % ? means phone is dropped
                    else
                        thisresp1=[thisresp1;'?']; % ? means phone is dropped
                        thisresp2=[thisresp2;for_parse];
                    end
                elseif any(~isnan(for_parse)) & length(for_parse)==2 % normal response CV
                    thisresp1=[thisresp1;for_parse(1)];
                    thisresp2=[thisresp2;for_parse(2)];  
                elseif any(~isnan(for_parse)) & length(for_parse)>2 % more than 2 phones
                    flag=1;
                end
            end
            
            if flag==0
            thisdata{temp,14}=sum(thisstim(1)==thisresp1)/length(thisresp1);
            thisdata{temp,15}=sum(thisstim(2)==thisresp2)/length(thisresp2);
            
            % populate CV CM for consonants
            for temp3=1:length(thisresp1)
                CV_consCM(find(thisstim(1)==cons_tested),find(thisresp1(temp3)==cons_tested))=CV_consCM(find(thisstim(1)==cons_tested),find(thisresp1(temp3)==cons_tested))+1;                
                total_CM(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))=total_CM(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))+1;
                CV_bigCM(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))=CV_bigCM(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))+1;
                if any(strcmp(names(n),disabled))
                    total_CM_RD(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))=total_CM_RD(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))+1; % reading control
                else
                    total_CM_RC(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))=total_CM_RC(find(thisstim(1)==cons_Dbet),find(thisresp1(temp3)==cons_Dbet))+1; % reading disabled
                end
            end
            end
        elseif any(thisstim(1)==vow_Dbet) % thisstim is a valid VC
            thisdata{temp,13}=thisstim(1);
            thisdata{temp,12}=thisstim(2);
            thisdata{temp,16}='VC';
            
            
            flag=0; % remove and fix>2 later
            for temp2=1:3
                for_parse=thisdata{temp,4+temp2};
                if any(~isnan(for_parse)) & length(for_parse)==1 % only one sound back
                    if any(for_parse==cons_Dbet)
                        thisresp1=[thisresp1;'?']; % ? means phone is dropped
                        thisresp2=[thisresp2;for_parse];
                    else
                        thisresp1=[thisresp1;for_parse];
                        thisresp2=[thisresp2;'?']; % ? means phone is dropped
                    end
                elseif any(~isnan(for_parse)) & length(for_parse)==2 % normal response VC
                    thisresp1=[thisresp1;for_parse(1)];
                    thisresp2=[thisresp2;for_parse(2)];     
                elseif any(~isnan(for_parse)) & length(for_parse)>2
                    flag=1;
                end
            end
            
            if flag==0
            thisdata{temp,15}=sum(thisstim(1)==thisresp1)/length(thisresp1);
            thisdata{temp,14}=sum(thisstim(2)==thisresp2)/length(thisresp2);         
            
            % populate VC CM for consonants
            for temp3=1:length(thisresp1)
                VC_consCM(find(thisstim(2)==cons_tested),find(thisresp2(temp3)==cons_tested))=VC_consCM(find(thisstim(2)==cons_tested),find(thisresp2(temp3)==cons_tested))+1;                
                total_CM(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))=total_CM(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))+1;
                VC_bigCM(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))=VC_bigCM(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))+1;
                if any(strcmp(names(n),disabled)) % if this child is learning disabled
                    total_CM_RD(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))=total_CM_RD(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))+1; % reading control
                else
                    total_CM_RC(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))=total_CM_RC(find(thisstim(2)==cons_Dbet),find(thisresp2(temp3)==cons_Dbet))+1; % reading disabled
                end
            end
            end
            %'check'
            %pause
        end % end if/then to determine CV vs VC
    end  % end filling in thisdata, iterating over all subjects
    'done thisdata'
    %pause
    
    % normalize the CMs, these CMs are for a single child
    CV_consCM=CV_consCM./repmat(sum(CV_consCM')',1,length(cons_tested));
    CV_consCM(isnan(CV_consCM))=0;
    CV_bigCM=CV_bigCM./repmat(sum(CV_bigCM')',1,length(cons_Dbet));
    CV_bigCM(isnan(CV_bigCM))=0;
    
    VC_consCM=VC_consCM./repmat(sum(VC_consCM')',1,length(cons_tested));
    VC_consCM(isnan(VC_consCM))=0;
    VC_bigCM=VC_bigCM./repmat(sum(VC_bigCM')',1,length(cons_Dbet));
    VC_bigCM(isnan(VC_bigCM))=0;
    
    total_consCM=CV_consCM+VC_consCM;
    Avg_consCM=total_consCM./repmat(sum(total_consCM')',1,length(cons_tested));
    
    
    
    % CM with diagonal zeroed out:
    conf=Avg_consCM-eye(length(cons_tested)).*Avg_consCM;
    %  CM with diagonal zeroed out:
    conf_cv=CV_consCM-eye(length(cons_tested)).*CV_consCM;
    conf_vc=VC_consCM-eye(length(cons_tested)).*VC_consCM;
    % get the child-specific order of confusions:
    most_conf=sort(reshape(conf,length(cons_tested)^2,1),'descend');
    most_conf=most_conf(1:15);
    
    matdim=length(cons_tested); % for easier use, the dimension of the matrix
    conf_mine=conf;
    labels={};
    cv_list=[];
    vc_list=[];
    % this loop just builds the label for this plot
    for plottemp=1:length(most_conf) % top confusions
        % find the next worst confusion
        loc=find(abs(conf_mine-most_conf(plottemp))<0.0001,1,'first');
        % locate x,y in conf matrix
        if mod(loc,matdim)==0
            yspot=matdim;
        else
            yspot=mod(loc,matdim);
        end
        xspot=ceil(loc/matdim);
        conf_mine(yspot,xspot)=0; % clear spot so that others can be found
        % store data for this confusion
        cv_list=[cv_list; conf_cv(yspot,xspot)];
        vc_list=[vc_list; conf_vc(yspot,xspot)];
        % set labels for bar plot
        labels=[labels; char(strcat(cons_tested(yspot),'>',cons_tested(xspot)))];
    end
    %
%     figure(1)
%     set(gca,'FontSize',14)
%     barh(cv_list)
%     xlim([0 1])
%     xlabel('Probability')
%     set(gca,'ydir','reverse')
%     set(gca,'YTickLabel',labels)
%     set(gca,'XTick',[0 .2 .4 .6 .8 1])
%     set(gca,'XGrid','on')
%     title(char(strcat('Frequent CV Confusions for ',{' '},names(n))))
    
    figure(5)
    set(gca,'FontSize',14)
    barh([cv_list vc_list],1)
    xlim([0 1])
    xlabel('Probability')
    set(gca,'ydir','reverse')
    set(gca,'YTickLabel',labels)
    set(gca,'XTick',[0 .2 .4 .6 .8 1])
    set(gca,'XGrid','on')
    title(char(strcat('Frequent Confusions for ',{' '},names(n))))
    legend({'CV','VC'},'Location','Best')
    colormap(Colors)
    
    if 1
        filename=char(strcat('figs/TopConf_CVandVC_15_IndividualTotalOrder_',names(n)));
        saveas(gcf,filename,'fig')
        set(gcf,'PaperPositionMode','auto')
        eval(char(strcat('print -depsc2 ',{' '},filename))) 
    end
    
    %pause
    
    % replacing above with population results
    %if any(strcmp(names(n),disabled)) % this is an RD child
        most_conf=RD_CM_vec(1:20); % THIS controls the order of the y axis
        group_conf=total_CM_RD_norm_conf;
    %else % child is control
    %    most_conf=RC_CM_vec(1:20);
    %    group_conf=total_CM_RC_norm_conf;
    %end
    matdim=length(cons_Dbet); % for easier use, this is the matrix dimension
    conf_mine=CV_bigCM;
    these_conf=[];
    labels={};
    cv_list=[];
    vc_list=[];
    nodata=[];
    for plottemp=1:length(most_conf) % top confusions        
        % find the next worst confusion
        loc=find(abs(group_conf-most_conf(plottemp))<0.0001,1,'first');
        % locate x,y in conf matrix
        if mod(loc,matdim)==0
            yspot=matdim;
        else
            yspot=mod(loc,matdim);
        end
        xspot=ceil(loc/matdim);
        % store this child's errors for plotting
        if ~isnan(conf_mine(yspot,xspot))
            these_conf=[these_conf; conf_mine(yspot,xspot)];
        else
            these_conf=[these_conf; 0];
            nodata=[nodata plottemp];
        end
        conf_mine(yspot,xspot)=0; % clear spot so that others can be found
        % store vc data for this confusion
        cv_list=[cv_list; CV_bigCM(yspot,xspot)];
        vc_list=[vc_list; VC_bigCM(yspot,xspot)];
        % set labels for bar plot
        labels=[labels; char(strcat(cons_Dbet(yspot),'>',cons_Dbet(xspot)))];
    end

%     figure(2)
%     set(gca,'FontSize',14)
%     barh(cv_list)
%     xlim([0 1])
%     ylim([0 length(cv_list)+1])
%     xlabel('Probability')
%     set(gca,'ydir','reverse')
%     set(gca,'YTick',[1:length(cv_list)])
%     set(gca,'YTickLabel',labels)
%     set(gca,'XTick',[0 .2 .4 .6 .8 1])
%     set(gca,'XGrid','on')
%     title(char(strcat('CV Confusions for ',{' '},names(n),{' '},'in RD Group Order')))
%     if ~isempty(nodata)
%         for temp=1:length(nodata)
%             text(nodata(temp),0.05,'no data','FontSize',14)
%         end
%     end
%     
%     if 0
%         %filename=char(strcat('figs/TopConf_CVonly_pop_',names(n)));
%         %filename=char(strcat('figs/TopConf_CVonly_grporder_',names(n)));
%         filename=char(strcat('figs/TopConf_CVonly_RDorder_',names(n)));
%         saveas(gcf,filename,'fig')
%         set(gcf,'PaperPositionMode','auto')
%         eval(char(strcat('print -depsc2 ',{' '},filename))) 
%     end
     
    figure(6)
    set(gca,'FontSize',14)
    barh([cv_list vc_list])
    xlim([0 1])
    ylim([0 length(cv_list)+1])
    xlabel('Probability')
    set(gca,'ydir','reverse')
    set(gca,'YTick',[1:length(cv_list)])
    set(gca,'YTickLabel',labels)
    set(gca,'XTick',[0 .2 .4 .6 .8 1])
    set(gca,'XGrid','on')
    title(char(strcat('Confusions for ',{' '},names(n),{' '},'in RD Group Order')))
    if ~isempty(nodata)
        for temp=1:length(nodata)
            text(nodata(temp),0.05,'no data','FontSize',14)
        end
    end
    legend({'CV','VC'},'Location','Best')
    colormap(Colors)
    
     if 1
        filename=char(strcat('figs/TopConf_CVandVC_20_RDtotalorder_',names(n)));
        saveas(gcf,filename,'fig')
        set(gcf,'PaperPositionMode','auto')
        eval(char(strcat('print -depsc2 ',{' '},filename))) 
    end
    
    'check'
    pause
        
end

% figure(3)
% total_CM_norm=total_CM./repmat(sum(total_CM')',1,length(cons_Dbet));
% imagesc(total_CM_norm)
% Colormap(1-gray)
% set(gca,'YTickLabel',cons_Dbet)
% set(gca,'XTickLabel',cons_Dbet)
% set(gca,'YTick',[1:length(cons_Dbet)])
% set(gca,'XTick',[1:length(cons_Dbet)])

'done with all subjects'
pause

figure(5)
%total_CM_RC_norm=total_CM_RC./repmat(sum(total_CM_RC')',1,length(cons_Dbet));
imagesc(total_CM_RC_norm)
set(gca,'FontSize',14)
Colormap(1-gray)
set(gca,'YTickLabel',cons_Dbet)
set(gca,'XTickLabel',cons_Dbet)
set(gca,'YTick',[1:length(cons_Dbet)])
set(gca,'XTick',[1:length(cons_Dbet)])
title('Population Confusion Matrix for Reading Controls')

figure(6)
%total_CM_RD_norm=total_CM_RD./repmat(sum(total_CM_RD')',1,length(cons_Dbet));
imagesc(total_CM_RD_norm)
set(gca,'FontSize',14)
Colormap(1-gray)
set(gca,'YTickLabel',cons_Dbet)
set(gca,'XTickLabel',cons_Dbet)
set(gca,'YTick',[1:length(cons_Dbet)])
set(gca,'XTick',[1:length(cons_Dbet)])
title('Population Confusion Matrix for Reading Disabled')



%% Plot the RD and RC group bar confusions (RD order)

% replacing above with population results
    most_conf=RD_CM_vec(1:20);
    group_conf=total_CM_RD_norm_conf; % to set ordering
    matdim=length(cons_Dbet); % for easier use
    conf_mine=total_CM_RC_norm_conf; % to plot
    these_conf=[];
    labels={};
    nodata=[];
    for plottemp=1:length(most_conf) % top confusions        
        % find the next  worst confusion
        loc=find(abs(group_conf-most_conf(plottemp))<0.0001,1,'first');
        % locate x,y in conf matrix
        if mod(loc,matdim)==0
            yspot=matdim;
        else
            yspot=mod(loc,matdim);
        end
        xspot=ceil(loc/matdim);
        % store this child's errors
        if ~isnan(conf_mine(yspot,xspot))
            these_conf=[these_conf; conf_mine(yspot,xspot)];
        else
            these_conf=[these_conf; 0];
            nodata=[nodata plottemp];
        end
        conf_mine(yspot,xspot)=0; % clear spot so that others can be found
        % set labels for bar plot
        labels=[labels; char(strcat(cons_Dbet(yspot),'>',cons_Dbet(xspot)))];
    end
    % Eventually have one CV and one VC bar
    figure(1)
    set(gca,'FontSize',14)
    barh(these_conf)
    xlim([0 1])
    ylim([0 length(these_conf)+1])
    xlabel('Probability')
    set(gca,'ydir','reverse')
    set(gca,'YTick',[1:length(these_conf)])
    set(gca,'YTickLabel',labels)
    set(gca,'XTick',[0 .2 .4 .6 .8 1])
    set(gca,'XGrid','on')
    title(char(strcat('Confusions for RC Group, RD Order')))




% replacing above with population results
    most_conf=RD_CM_vec(1:20);
    group_conf=total_CM_RD_norm_conf; % to set ordering
    matdim=length(cons_Dbet); % for easier use
    conf_mine=total_CM_RD_norm_conf; % to plot
    these_conf=[];
    labels={};
    nodata=[];
    for plottemp=1:length(most_conf) % top confusions        
        % find the next worst confusion
        loc=find(abs(group_conf-most_conf(plottemp))<0.0001,1,'first');
        % locate x,y in conf matrix
        if mod(loc,matdim)==0
            yspot=matdim;
        else
            yspot=mod(loc,matdim);
        end
        xspot=ceil(loc/matdim);
        % store this child's errors
        if ~isnan(conf_mine(yspot,xspot))
            these_conf=[these_conf; conf_mine(yspot,xspot)];
        else
            these_conf=[these_conf; 0];
            nodata=[nodata plottemp];
        end
        conf_mine(yspot,xspot)=0; % clear spot so that others can be found
        % set labels for bar plot
        labels=[labels; char(strcat(cons_Dbet(yspot),'>',cons_Dbet(xspot)))];
    end
    % Eventually have one CV and one VC bar
    figure(2)
    set(gca,'FontSize',14)
    barh(these_conf)
    xlim([0 1])
    ylim([0 length(these_conf)+1])
    xlabel('Probability')
    set(gca,'ydir','reverse')
    set(gca,'YTick',[1:length(these_conf)])
    set(gca,'YTickLabel',labels)
    set(gca,'XTick',[0 .2 .4 .6 .8 1])
    set(gca,'XGrid','on')
    title(char(strcat('Confusions for RD Group, RD Order')))

%% plot the major differences in confusions between RC and RD 

diff_CM_conf=total_CM_RD_norm_conf-total_CM_RC_norm_conf;
diff_CM_vec=sort(reshape(diff_CM_conf,length(cons_Dbet)^2,1),'descend');

% replacing above with population results
    most_conf=[diff_CM_vec(1:10); diff_CM_vec(end-10:end)];
    group_conf=diff_CM_conf; % to set ordering
    matdim=length(cons_Dbet); % for easier use
    conf_mine=diff_CM_conf; % to plot
    these_conf=[];
    labels={};
    nodata=[];
    for plottemp=1:length(most_conf) % top confusions        
        % find the next worst confusion
        loc=find(abs(group_conf-most_conf(plottemp))<0.0001,1,'first');
        % locate x,y in conf matrix
        if mod(loc,matdim)==0
            yspot=matdim;
        else
            yspot=mod(loc,matdim);
        end
        xspot=ceil(loc/matdim);
        % store this child's errors
        if ~isnan(conf_mine(yspot,xspot))
            these_conf=[these_conf; conf_mine(yspot,xspot)];
        else
            these_conf=[these_conf; 0];
            nodata=[nodata plottemp];
        end
        conf_mine(yspot,xspot)=0; % clear spot so that others can be found
        % set labels for bar plot
        labels=[labels; char(strcat(cons_Dbet(yspot),'>',cons_Dbet(xspot)))];
    end
    % Eventually have one CV and one VC bar
    figure(2)
    set(gca,'FontSize',14)
    barh(these_conf)
    xlim([-.5 .5])
    ylim([0 length(these_conf)+1])
    xlabel('Probability')
    set(gca,'ydir','reverse')
    set(gca,'YTick',[1:length(these_conf)])
    set(gca,'YTickLabel',labels)
    set(gca,'XTick',[-.5:.1:.5])
    set(gca,'XGrid','on')
    title(char(strcat('RD-RC group avg')))

    





