clear all
addpath ./mat

CONS = {'C','D','G','J','S','T','Z','b','d','f','g','h','k','L','m','n','p',...
    'r','s','t','v','w','y','z'};
nCONS = length(CONS);
VOWS = {'@','A','E','I','O','R','U','W','Y','a','c','e','i','o','u'};
nVOWS = length(VOWS);

frictivesAffricates = { 'C', 'D', 'h', 'J', 'S', 'T', 'Z', 'f', 's', 'v', 'z'};
iFrictivesAffricates = [ismember(CONS,frictivesAffricates) zeros(1,nVOWS)];
stops = { 'b', 'd', 'g', 'k', 'p', 't' };
iStops = [ismember(CONS,stops)  zeros(1,nVOWS)];
nasalsLiquidsGlides = {'G', 'L', 'm', 'n', 'r', 'w', 'y'};
iNasalsLiquidsGlides = [ismember(CONS,nasalsLiquidsGlides)  zeros(1,nVOWS)];
vowels = VOWS;
iVowels = [zeros(1,nCONS)  ones(1,nVOWS)];
groups = {'frictivesAffricates','stops','nasalsLiquidsGlides','vowels'};

RD = {'Alina','Angela','Edward','Latisha','Laura','Norene',...
    'Shauna','Teddy','Tony'};
RC = {'Anton','Bob','Carly','Evan','Joanna','Miguel'};
SUBJECTS = {RD{:}, RC{:}};
nSUBJECTS = length(SUBJECTS);

table =0
if table
    output = cell(nCONS + nVOWS+2,nSUBJECTS*2+1); % add two rows and one column for titles
    % add target labels
    conIdx = 3:2+nCONS;
    output(conIdx,1) = CONS;
    vowIdx = 3+nCONS:size(output,1);
    output(vowIdx,1) = VOWS;
    output(1,2:2:end)=SUBJECTS;
    hIdx = 2:2:size(output,2);
    output(2,hIdx)=cellstr(repmat('H',nSUBJECTS,1));
    nIdx = 3:2:size(output,2);
    output(2,nIdx)=cellstr(repmat('N',nSUBJECTS,1));
    outFinal = output;
    outInitial = output;
    
    for i = 1 : nSUBJECTS
        subject = SUBJECTS(i);
        eval(['load ' subject{:} '.mat'])
        cm = CMs{1}; % CMs 1-2 Cons, 3-4 Vows (odd Initial, even Final),
        % Cons Initial
        N = sum(cm,2);
        p = cm ./ repmat(N,1,size(cm,2) );
        H = sum(log2(p.^(-p)),2);
        outInitial(conIdx,hIdx(i)) = num2cell(H);
        outInitial(conIdx,nIdx(i)) = num2cell(N);
        cm = CMs{2}; % CMs 1-2 Cons, 3-4 Vows (odd Initial, even Final),
        % Cons Final
        N = sum(cm,2);
        p = cm ./ repmat(N,1,size(cm,2) );
        H = sum(log2(p.^(-p)),2);
        outFinal(conIdx,hIdx(i)) = num2cell(H);
        outFinal(conIdx,nIdx(i)) = num2cell(N);
        
        cm = CMs{3}; % CMs 1-2 Cons, 3-4 Vows (odd Initial, even Final),
        % Vows Initial
        N = sum(cm,2);
        p = cm ./ repmat(N,1,size(cm,2) );
        H = sum(log2(p.^(-p)),2);
        outInitial(vowIdx,hIdx(i)) = num2cell(H);
        outInitial(vowIdx,nIdx(i)) = num2cell(N);
        cm = CMs{4}; % CMs 1-2 Cons, 3-4 Vows (odd Initial, even Final),
        % Cons Final
        N = sum(cm,2);
        p = cm ./ repmat(N,1,size(cm,2) );
        H = sum(log2(p.^(-p)),2);
        outFinal(vowIdx,hIdx(i)) = num2cell(H);
        outFinal(vowIdx,nIdx(i)) = num2cell(N);
    end
    
    %% write csv file (only works with database toolbox)
    % outFinal = cell2dataset(outFinal);
    % export(outFinal,'file','final.csv','delimiter',',')
    %
    % outInitial = cell2dataset(outInitial);
    % export(outInitial,'file','initial.csv','delimiter',',')
    
    %% write csv file (works in every Matlab)
    fid = fopen('final.csv','w')
    numColumns = size(outFinal,2);
    % use repmat to construct repeating formats
    % ( numColumns-1 because no comma on last one)
    headerFmt = repmat('%s,',1,numColumns-1);
    numFmt = repmat('%f,',1,numColumns-2);
    numFmt = ['%s,' numFmt]; % add targets
    fprintf(fid,[headerFmt,'%s\n'],outFinal{1,:})
    fprintf(fid,[headerFmt,'%s\n'],outFinal{2,:})
    for i = 3 : size(outFinal,1)
        fprintf(fid,[numFmt,'%f\n'],outFinal{i,:})
    end
    fclose(fid)
    
    fid = fopen('initial.csv','w')
    numColumns = size(outInitial,2);
    % use repmat to construct repeating formats
    % ( numColumns-1 because no comma on last one)
    headerFmt = repmat('%s,',1,numColumns-1);
    numFmt = repmat('%f,',1,numColumns-2);
    numFmt = ['%s,' numFmt]; % add targets
    fprintf(fid,[headerFmt,'%s\n'],outInitial{1,:})
    fprintf(fid,[headerFmt,'%s\n'],outInitial{2,:})
    for i = 3 : size(outInitial,1)
        fprintf(fid,[numFmt,'%f\n'],outInitial{i,:})
    end
    fclose(fid)
end

%% plot H(P_e)
set(0,'DefaultAxesFontSize',16)
ci = 1; % consonant initial
cf = 2; % consonant initial
vi = 3; % vowel initial
vf = 4; % vowel final
rd = 1; % reading disabeled
rc = 2; % reading control

% prepare gourp CMs
groupCM{rd,ci} = zeros(nCONS); groupCM{rc,ci} = zeros(nCONS);
groupCM{rd,cf} = zeros(nCONS); groupCM{rc,cf} = zeros(nCONS);
groupCM{rd,vi} = zeros(nVOWS); groupCM{rc,vi} = zeros(nVOWS);
groupCM{rd,vf} = zeros(nVOWS); groupCM{rc,vf} = zeros(nVOWS);

% create entropy and Pe matrix for RD group
for i = 1 : nSUBJECTS
    subject = SUBJECTS(i);
    eval(['load ' subject{:} '.mat'])
    if i > length(RD)
        idxRDRC = rc;
    else
        idxRDRC = rd;
    end
    % groupCM 1st dimension RD or RC, second dimension 1-2 Cons, 3-4 Vows (odd Initial, even Final)
    groupCM{idxRDRC,ci} = groupCM{idxRDRC,ci} + CMs{ci}; % CMs 1-2 Cons, 3-4 Vows (odd Initial, even Final),
    groupCM{idxRDRC,cf} = groupCM{idxRDRC,cf} + CMs{cf};
    groupCM{idxRDRC,vi} = groupCM{idxRDRC,vi} + CMs{vi};
    groupCM{idxRDRC,vf} = groupCM{idxRDRC,vf} + CMs{vf};
    
end

% Initial targets, first dimension RD or RC, second dimension target
entropy_initial = zeros(2, nCONS +nVOWS);
error_initial = zeros(2, nCONS +nVOWS);
entropy_final = zeros(2, nCONS +nVOWS);
error_final = zeros(2, nCONS +nVOWS);


% calculate H and Pe
for idxRDRC = 1:2
    % consonants initial
    cm = groupCM{idxRDRC,ci};
    N = sum(cm,2);
    p = cm ./ repmat(N,1,size(cm,2) ); % normalize CM by number of trials
    pc = diag(p);
    pe = ones(length(pc),1) - pc;
    H = sum(log2(p.^(-p)),2);
    entropy_initial(idxRDRC,1:nCONS) = H';
    error_initial(idxRDRC,1:nCONS) = pe';
    % vowels initial
    cm = groupCM{idxRDRC,vi};
    N = sum(cm,2);
    p = cm ./ repmat(N,1,size(cm,2) ); % normalize CM by number of trials
    pc = diag(p);
    pe = ones(length(pc),1) - pc;
    H = sum(log2(p.^(-p)),2);
    entropy_initial(idxRDRC,nCONS+1:end) = H';
    error_initial(idxRDRC,nCONS+1:end) = pe';
    % consonants final
    cm = groupCM{idxRDRC,cf};
    N = sum(cm,2);
    p = cm ./ repmat(N,1,size(cm,2) ); % normalize CM by number of trials
    pc = diag(p);
    pe = ones(length(pc),1) - pc;
    H = sum(log2(p.^(-p)),2);
    entropy_final(idxRDRC,1:nCONS) = H';
    error_final(idxRDRC,1:nCONS) = pe';
    % vowels final
    cm = groupCM{idxRDRC,vf};
    N = sum(cm,2);
    p = cm ./ repmat(N,1,size(cm,2) ); % normalize CM by number of trials
    pc = diag(p);
    pe = ones(length(pc),1) - pc;
    H = sum(log2(p.^(-p)),2);
    entropy_final(idxRDRC,nCONS+1:end) = H';
    error_final(idxRDRC,nCONS+1:end) = pe';
end

%% Plot

entropy_threshold = 1
error_threshold = 0.25
x=[1:100]/100;
H1=log2((x).^(-x))+log2(((1-x)).^(-(1-x)));
H2=log2((x).^(-x))+2*log2((((1-x)/2)).^(-(1-x)/2));
H3=log2((x).^(-x))+3*log2((((1-x)/3)).^(-(1-x)/3));
H4=log2((x).^(-x))+4*log2((((1-x)/4)).^(-(1-x)/4));
% initial
for group = 1:length(groups)
    iGroup = logical(eval(['i' ,regexprep(groups{group},'(\<[a-z])','${upper($1)}')]))
    figure
    plot(error_initial(rd,iGroup),entropy_initial(rd,iGroup),'xr');hold on;
    plot(error_initial(rc,iGroup),entropy_initial(rc,iGroup),'xb');
    l=legend('Reading Disabled','Reading Control')
    legend('boxoff')
    set(l,'Position',[0.5 0.15 0.035 0.05]);
    title(['H(p_e) for Targets (' groups{group} ') in Initial Position'])
    ylabel('H [bits]')
    xlabel('p_e')
    % plot reference lines
    plot(x,H1,'LineStyle','-','Color',[0.6 0.6 0.6]);hold on;
    plot(1-x,H2,'LineStyle','-','Color',[0.6 0.6 0.6]);
    plot(1-x,H3,'LineStyle','-','Color',[0.6 0.6 0.6]);
    plot(1-x,H4,'LineStyle','-','Color',[0.6 0.6 0.6]);
    % label data points
    labels = eval(groups{group});
    ixGroup= find(iGroup==1)
    for idx = find(entropy_initial(rd,iGroup) > entropy_threshold & error_initial(rd,iGroup) > error_threshold) 
        text(error_initial(rd,ixGroup(idx)) + 0.01, entropy_initial(rd,ixGroup(idx)) + 0.1 , labels{idx});
    end
    for idx = find(entropy_initial(rc,iGroup) > entropy_threshold & error_initial(rc,iGroup) > error_threshold)
        text(error_initial(rc,ixGroup(idx)) + 0.01, entropy_initial(rc,ixGroup(idx)) + 0.1, labels{idx});
    end
    
    fname = sprintf('HPeInitial%s.eps',regexprep(groups{group},'(\<[a-z])','${upper($1)}'));
    print(gcf, fullfile('./fig/', fname), '-depsc2');
end
% final
for group = 1:length(groups)
    iGroup = logical(eval(['i' regexprep(groups{group},'(\<[a-z])','${upper($1)}')]));
    
    figure
    plot(error_final(rd,iGroup),entropy_final(rd,iGroup),'xr');hold on;
    plot(error_final(rc,iGroup),entropy_final(rc,iGroup),'xb');
    l=legend('Reading Disabled','Reading Control')
    legend('boxoff')
    set(l,'Position',[0.5 0.15 0.035 0.05]);
    title(['H(p_e) for Targets (' groups{group} ') in Final Position'])
    ylabel('H [bits]')
    xlabel('p_e')
    % plot reference lines
    plot(x,H1,'LineStyle','-','Color',[0.6 0.6 0.6]);hold on;
    plot(1-x,H2,'LineStyle','-','Color',[0.6 0.6 0.6]);
    plot(1-x,H3,'LineStyle','-','Color',[0.6 0.6 0.6]);
    plot(1-x,H4,'LineStyle','-','Color',[0.6 0.6 0.6]);
    % label data points
    labels = eval(groups{group});
    ixGroup= find(iGroup==1)
    for idx = find(entropy_final(rd,iGroup) > entropy_threshold & error_final(rd,iGroup) > error_threshold)
        text(error_final(rd,ixGroup(idx)) + 0.01, entropy_final(rd,ixGroup(idx)) + 0.1 , labels{idx});
    end
    for idx = find(entropy_final(rc,iGroup) > entropy_threshold & error_final(rc,iGroup) > error_threshold)
        text(error_final(rc,ixGroup(idx)) + 0.01, entropy_final(rc,ixGroup(idx)) + 0.1, labels{idx});
    end
    fname = sprintf('HPeFinal%s.eps',regexprep(groups{group},'(\<[a-z])','${upper($1)}'));
    print(gcf, fullfile('./fig/', fname), '-depsc2');
end