
clear all
close all
clc


%% Load dataset
disp('Loading MN64 dataset...')

load('EDS_MN64.mat')
npts = length(MN64.resp);

listenerLabels = unique(cellstr(MN64.listener(:,1:3)));
listenerLabels = listenerLabels(1:length(listenerLabels)-1); % get rid of 'nob'
nlisteners = length(listenerLabels);

respLabels = unique(MN64.resp(:,1));
respLabels = respLabels(2:length(respLabels)); % get rid of '?'
nresponses = length(respLabels);

talkerLabels = unique(cellstr(MN64.talker));
ntalkers = length(talkerLabels);

consLabels = unique(cellstr(MN64.stim(:,1)));
ncons = length(consLabels);

vowelLabels = unique(cellstr(MN64.stim(:,2)));
nvowels = length(vowelLabels);

utteranceLabels = unique(strcat(cellstr(MN64.talker),cellstr(MN64.stim)));
nutterances = length(utteranceLabels);

snrLabels = unique(cellstr(num2str(MN64.snr)));
nsnr = length(snrLabels);


%% Get response probabilities per consonant
disp('Getting responses...')

responses = zeros(nlisteners,nsnr,nutterances,nresponses);

for row = 1:npts

    rowdata = {MN64.listener(row,1:3) ... % listener
               strcat(MN64.talker(row,:),MN64.stim(row,:)) ... % utterance
               MN64.resp(row,1) ... % response
               MN64.snr(row) }; % snr

    listener = strcmp(rowdata(1),listenerLabels);
    utterance = strcmp(rowdata(2),utteranceLabels);
    response = strcmp(rowdata(3),respLabels);
    snr = cell2mat(rowdata(4))==str2num(char(snrLabels));

    responses(listener,snr,utterance,response) = responses(listener,snr,utterance,response) + 1;
                
end


% compute Pe and entropy
disp('Computing Pe and entropy...')

Pe = zeros(nsnr,nutterances);
Hs = zeros(nsnr,nutterances);
Hn = zeros(nsnr,nutterances);
Hgroup = zeros(nsnr,nutterances);

for snr = 1:nsnr
    for utterance = 1:nutterances
            
                n = sum(sum(responses(:,snr,utterance,:)));

                if n > 0
                
                    [v correctResp] = max(utteranceLabels{utterance}(5)==respLabels);
                    ncorrect = sum(sum(responses(:,snr,utterance,correctResp)));
                    nerror = n - ncorrect;
                
                    rowsum = sum(responses(:,snr,utterance,:));
                    rowprop = rowsum / n;

                    Pe(snr,utterance) = nerror / n;
                    Hs(snr,utterance) = -sum(rowprop(rowprop>0).*log2(rowprop(rowprop>0)));
                        
                    i = 1;
                    dummyResp = zeros(1,nresponses-1);
                    nerrorCounter = nerror;
                    while nerrorCounter>0
                        dummyResp(i) = dummyResp(i) + 1;
                        i = i + 1;
                        if i > nresponses -1
                            i = 1;
                        end
                        nerrorCounter = nerrorCounter - 1;
                    end
                    dummyResp(nresponses) = ncorrect;
                    dummyrowprop = dummyResp./n;
                    Hm = -sum(dummyrowprop(dummyrowprop>0).*log2(dummyrowprop(dummyrowprop>0)));
                    Hn(snr,utterance) = Hs(snr,utterance) ./ Hm;
                    
                    % 2 modes
                    dummyResp = zeros(1,2);
                    dummyResp(1) = ncorrect;
                    dummyResp(2) = nerror;
                    dummyrowprop = dummyResp./n;
                    H2 = -sum(dummyrowprop(dummyrowprop>0).*log2(dummyrowprop(dummyrowprop>0)));
                    
                    % 3 modes
                    dummyResp = zeros(1,3);
                    dummyResp(1) = ncorrect;
                    dummyResp(2) = nerror/2;
                    dummyResp(3) = nerror/2;
                    dummyrowprop = dummyResp./n;
                    H3 = -sum(dummyrowprop(dummyrowprop>0).*log2(dummyrowprop(dummyrowprop>0)));
                    
                    % 4 modes
                    %{
                    dummyResp = zeros(1,4);
                    dummyResp(1) = ncorrect;
                    dummyResp(2) = nerror/3;
                    dummyResp(3) = nerror/3;
                    dummyResp(4) = nerror/3;
                    dummyrowprop = dummyResp./n;
                    H4 = -sum(dummyrowprop(dummyrowprop>0).*log2(dummyrowprop(dummyrowprop>0)));
                    %}
                    
                    Hmodes = [0, H2, H3, Hm];
                    [v I] = min(abs(Hmodes-Hs(snr,utterance)));
                    Hgroup(snr,utterance) = I;
                    
                else
                    Pe(snr,utterance) = NaN;
                end
                
    end
end

% Output data for R
longresp = cell(nsnr*nutterances,8);
row = 1;
for snr = 1:nsnr
    for utterance = 1:nutterances
        longresp(row,1) = snrLabels(snr);
        longresp(row,2) = {utteranceLabels{utterance}(1:4)}; % talker
        longresp(row,3) = {utteranceLabels{utterance}(5)}; % cons
        longresp(row,4) = {utteranceLabels{utterance}(6)}; % vowel
        longresp(row,5) = {Pe(snr,utterance)};
        longresp(row,6) = {Hs(snr,utterance)};
        longresp(row,7) = {Hn(snr,utterance)};
        longresp(row,8) = {Hgroup(snr,utterance)};
        row = row + 1;
    end        
end

fid = fopen('mn64data.2.txt','w')
fprintf(fid,'snr\ttalker\tcons\tvowel\tpe\ths\thn\thgroup\n');
for i = 1:length(longresp)
    fprintf(fid,'%s\t%s\t%s\t%s\t%f\t%f\t%f\t%i\n',longresp{i,1},longresp{i,2},longresp{i,3},longresp{i,4},longresp{i,5},longresp{i,6},longresp{i,7},longresp{i,8});
end
fclose(fid);

