% Code to analyse test results from MUSHRA/HEMVIP
% (c) Gustav Eje Henter 2014, 2015, 2016, 2018, 2019, 2020, 2022, 2023

function analyse_hemvip_data(infile,outpath,typestr,...
    pval,nexp,labelstrs,neworder,includesubj,aspect,ignoreinitial)

converttoranks = false;
normaliserange = false;
normalisefixed = [];
prunetype = 0; % 0 => no pruning; 1 => prune ref=0; 2 => prune ref<100

showsubjbox = false; % Set to true to plot each subject's ratings

whiskerq = 0.025; % Fraction of datapoints outside each whisker
showboxmeans = true; % Show squares for means in box plot
showboxmconfs = true; % Show confidence intervals in box plot
%outliersymbol = '+'; % Set to '' to disable plotting outliers
outliersymbol = '';
jitterfrac = 0.5;
bpfontsize = 14;
xangle = 0; % Rotation angle for system labels on x-axis

% Permute variables to use new order

ncond = numel(neworder);

labelstrs = labelstrs(neworder);

invorder(neworder) = (1:ncond);

% Extract input file name without path or suffix
filestr = split(infile,filesep);
filestr = filestr{end}; % Only keep filename, not path
filestr = filestr(1:max(strfind(filestr,'.'))-1); % Remove (last) suffix

jsn = jsondecode(fileread(infile));

nsubj = numel(jsn);

subjresults = zeros(nexp,ncond,nsubj);
stimnums = zeros(nexp,ncond,nsubj);
allresults = [];
allstimnums = [];

for n = 1:nsubj
    [subjresults(:,:,n),stimnums(:,:,n)]...
        = json2result(jsn(n).trials,labelstrs,ignoreinitial);
    
    if showsubjbox
        enhanced_boxplot(subjresults(:,:,n),labelstrs,...
            whiskerq,showboxmeans,outliersymbol,jitterfrac);
        
        axis([0.5 (ncond+0.5) 0 100]);
        title(['Subject number ' int2str(n)]);
        
        subjmultimin = sum(sum(subjresults(:,:,n) == 0,2) > 1);
        subjmultimax = sum(sum(subjresults(:,:,n) == 100,2) > 1);
        if (subjmultimin > 0) || (subjmultimax > 0)
            disp(['Subject number: ' int2str(n)]);
            disp(['Number of screens with more than one 100 rating: '...
                int2str(subjmultimax)]);
            disp(['Number of screens with more than one   0 rating: '...
                int2str(subjmultimin)]);
        end
        
        disp('Paused. Press any key to continue...');
        pause;
    end
    
    if includesubj(n)
        allresults = [allresults;subjresults(:,:,n)];
        allstimnums = [allstimnums;stimnums(:,:,n)];
    end
end

% Prune and process data

unnat = (allresults(:,1) == 0);
incorr = (allresults(:,1) < 100); % No system rated as perfectly natural

if (prunetype == 2) % Only use screens where natural == 100
    corrresults = allresults(~incorr,:);
    corrstimnums = allstimnums(~incorr,:);
elseif (prunetype == 1) % Only use screens where natural > 0
    corrresults = allresults(~unnat,:);
    corrstimnums = allstimnums(~unnat,:);
else % Use all screens/pages
    corrresults = allresults;
    corrstimnums = allstimnums;
end

nkept = sum(includesubj);
nretained = size(corrresults,1);

% Round to integer and limit range
nanmask = isnan(corrresults);
corrresults = round(corrresults);
corrresults = max(0,corrresults);
corrresults = min(100,corrresults);

corrresults(nanmask) = NaN;
corrstimnums(nanmask) = NaN;

scalerange = [0 100];

% Convert to ranks
if converttoranks
    for n = 1:nretained
        corrresults(n,:) = tiedrank(corrresults(n,:));
    end
    scalerange = [1 ncond];
end

% Scale normalisation
if ~isempty(normalisefixed)
    % "Fixed" normalisation
    mincr = corrresults(:,normalisefixed(1));
    maxcr = corrresults(:,normalisefixed(2));
    
    corrresults = 100*(corrresults - repmat(mincr,1,ncond))...
        ./repmat(maxcr - mincr,1,ncond);
elseif normaliserange
    
    % Adaptive normalisation
    mincr = min(corrresults,[],2);
    maxcr = max(corrresults,[],2);
    
    corrresults = 100*(corrresults - repmat(mincr,1,ncond))...
        ./repmat(maxcr - mincr,1,ncond);
end

% Create output directory
if ~exist(outpath,'dir')
    mkdir('.',outpath);
end

% Compute result statistics for each stimulus
sns = unique(stimnums);
sns = sns(~isnan(sns))';

stimmeans = zeros(numel(sns),ncond);
stimstds = zeros(numel(sns),ncond);
stimmedians = zeros(numel(sns),ncond);
stimnumresults = zeros(numel(sns),ncond);
for sn = sns
    for c = 1:ncond
        mask = (corrstimnums(:,c) == sn);
        stimresults = corrresults(mask,c);
        
        % Add 1 to sn here since sequence numbers are zero-indexed
        stimnumresults(sn+1,c) = numel(stimresults);
        stimmeans(sn+1,c) = mean(stimresults);
        stimstds(sn+1,c) = std(stimresults,1);
        stimmedians(sn+1,c) = median(stimresults);
    end
end

% Create a boxplot of ratings and confidence intervals
clf;
absstats = enhanced_boxplot(corrresults(:,1:end),labelstrs(1:end),...
    whiskerq,showboxmeans,outliersymbol,jitterfrac,pval*showboxmconfs);

axis([0.5 (ncond + 0.5) scalerange]);

xtickangle(xangle);

set(gca,'FontSize',bpfontsize);
txts = findobj(gca,'Type','text');
set(txts,'FontSize',bpfontsize);
set(txts,'VerticalAlignment','Middle');
clear txts;

if ~converttoranks
    if normaliserange
        ylabel(['Normalised ' typestr ' rating'],'FontSize',bpfontsize);
    else
        ylabel([upper(typestr(1)) typestr(2:end) ' rating'],...
            'FontSize',bpfontsize);
    end
else
    ylabel([upper(typestr(1)) typestr(2:end) ' rank (highest is best)']...
        ,'FontSize',bpfontsize);
end

saveas(gcf,[outpath filesep filestr '_boxplot.fig'],'fig');
%saveas(gcf,[outpath filesep filestr '_boxplot.eps'],'epsc');
saveas(gcf,[outpath filesep filestr '_boxplot.png'],'png');

fig2fm([outpath filesep filestr '_boxplot'],'article',...
    'a4paper,british,12pt',17,aspect,9,[],[],5,[],[],true);

% Print sort order to use in order to sort by median (largest to smallest)
if ~all(sort(absstats(:,4)','descend') == absstats(:,4))
    [~,sortorder] = sort(absstats(:,4)','descend');
    fprintf(1,'neworder = [%s];\n', join(string(sortorder), ' '));
end

% Print basic LaTeX code for results table
disp('Code for creating LaTeX code for the tables in the paper:');
disp('');

disp('\begin{tabular}{@{}l|cc@{}}');
disp('\toprule');
disp(['& \multicolumn{2}{c@{}}{' upper(typestr(1)) typestr(2:end) '} \\']);
disp('ID & Median & Mean\\');
disp('\midrule');
for c = 1:ncond
    fprintf('%s & $%d\\in{}[%d,\\,%d]$ & $%.1f{\\pm}%.1f$\\\\\n',...
        labelstrs{c},[absstats(c,(4:6)),round(absstats(c,1),1),...
        ceil(10*(absstats(c,3)-absstats(c,1)))/10]);
end
disp('\bottomrule');
disp('\end{tabular}');
disp(' ');

% Perform statistical analysis of system pairs
[numcmps,dmean,dmedian,nxgty] = deal(zeros(ncond));
[pt,pwilc,pmwupaired,...
    ptunpaired,pmwuunpaired,psgn,prefnoties] = deal(eye(ncond)/2);

ncondratings = sum(~isnan(corrresults),1);

for c1 = 1:ncond
    for c2 = (c1+1):ncond
        ds = corrresults(:,c1) - corrresults(:,c2);
        numcmps(c1,c2) = sum(~isnan(ds));
        ds = ds(~isnan(ds));
        
        dmean(c1,c2) = mean(ds);
        dmedian(c1,c2) = median(ds);
        
        validpairs = [corrresults(:,c1),corrresults(:,c2)];
        validpairs = validpairs(~any(isnan(validpairs),2),:);
        
        [~,pt(c1,c2)] = ttest(validpairs(:,1),validpairs(:,2));
        [pwilc(c1,c2)] = signrank(validpairs(:,1),validpairs(:,2));        
        [pmwupaired(c1,c2)] = ranksum(validpairs(:,1),validpairs(:,2));
        
        [~,ptunpaired(c1,c2)]...
            = ttest2(corrresults(~isnan(corrresults(:,c1)),c1),...
            corrresults(~isnan(corrresults(:,c2)),c2),...
            'Vartype','unequal');
        [pmwuunpaired(c1,c2)]...
            = ranksum(corrresults(~isnan(corrresults(:,c1)),c1),...
            corrresults(~isnan(corrresults(:,c2)),c2));
        
        nxgty(c1,c2) = sum(validpairs(:,1) > validpairs(:,2));
        nxgty(c2,c1) = sum(validpairs(:,1) < validpairs(:,2));
        
        nxy = nxgty(c1,c2) + nxgty(c2,c1);
        minxy = min(nxgty(c1,c2),nxgty(c2,c1));
        [psgn(c1,c2)] = binocdf(minxy,nxy,0.5)...
            + binocdf(nxy - minxy,nxy,0.5,'upper');
        
        prefnoties(c1,c2) = nxgty(c1,c2)/nxy;
        prefnoties(c2,c1) = nxgty(c2,c1)/nxy;
    end
end

mincmps = min(numcmps(numcmps > 0)),
maxcmps = max(numcmps(numcmps > 0)),
numcmps = numcmps + numcmps';

dmean = dmean - dmean';
dmedian = dmedian - dmedian';

prefties = nxgty./nretained + eye(ncond)./2;
effsizesgn = prefties - prefties';

pt = pt + pt';
pwilc = pwilc + pwilc';
pmwupaired = pmwupaired + pmwupaired';
ptunpaired = ptunpaired + ptunpaired';
pmwuunpaired = pmwuunpaired + pmwuunpaired';
psgn = psgn + psgn';

ncmp = ncond*(ncond-1)/2;
[rejt,pcorrt] = holmbonferroni(pt,pval);
[rejwilc,pcorrwilc] = holmbonferroni(pwilc,pval);
[rejmwupaired,pcorrmwupaired] = holmbonferroni(pmwupaired,pval);
[rejtunpaired,pcorrtunpaired] = holmbonferroni(ptunpaired,pval);
[rejmwuunpaired,pcorrmwuunpaired] = holmbonferroni(pmwuunpaired,pval);
[rejsgn,pcorrsgn] = holmbonferroni(psgn,pval);

triumask = triu(true(ncond),1);

% Plot the balancing (number of comparisons for distinct condition pairs)
plot(sort(numcmps(triumask(:))));

xlabel('Pair comparison number (sorted)');
ylabel('Number of contrasts gathered');

saveas(gcf,[outpath filesep filestr '_balancing.fig'],'fig');
%saveas(gcf,[outpath filesep filestr '_balancing.eps'],'epsc');
saveas(gcf,[outpath filesep filestr '_balancing.png'],'png');

% Plot the distribution of p-values for different tests
semilogy(1:ncmp,sort(pcorrt(triumask(:)),'descend'),'-',...
    1:ncmp,sort(pcorrwilc(triumask(:)),'descend'),'-',...
    1:ncmp,sort(pcorrmwupaired(triumask(:)),'descend'),'-',...
    1:ncmp,sort(pcorrsgn(triumask(:)),'descend'),'-',...
    1:ncmp,sort(pcorrtunpaired(triumask(:)),'descend'),'.-',...
    1:ncmp,sort(pcorrmwuunpaired(triumask(:)),'descend'),'.-');
hold on;
semilogy([1 ncmp],pval*[1 1],'k:');
hold off;

curraxis = axis;
axis([1 ncmp 10e-8 1]);

xlabel('Pair comparison number (sorted)');
ylabel('Corrected p-value');
legend('t-test (paired)','Wilcoxon signed-rank (paired)',...
    'Mann-Whitney U ("paired")','Sign test (paired)',...
    't-test (unpaired)','Mann-Whitney U (unpaired)',...
    'Limit of significance','Location','NorthEast');

saveas(gcf,[outpath filesep filestr '_corr_pvals.fig'],'fig');
%saveas(gcf,[outpath filesep filestr '_corr_pvals.eps'],'epsc');
saveas(gcf,[outpath filesep filestr '_corr_pvals.png'],'png');

%pcorrwilc, % Holm-Bonferroni corrected Wilcoxon signed-rank test p-values

% Print the results of the statistical tests and their differences

%rejwilc, % Show a matrix of rejected Wilcoxon signed rank null hypotheses
rejmwuunpaired, % Show a matrix of rejected Mann-Whitney U null hypotheses

%disp('Non-significant pairs (Wilcoxon signed-rank test):');
%flaggedpairs(~rejwilc,labelstrs,true);

disp('Non-significant pairs (unpaired Mann-Whitney U test):');
flaggedpairs(~rejmwuunpaired,labelstrs,true);

%if converttoranks || normaliserange
%    mwupairedwilcdisagree = rejmwupaired - rejwilc,
%    disp(['Pairs where Wilcoxon signed-rank test '...
%        'and paired Mann-Whitney U test disagree:']);
%    flaggedpairs(mwupairedwilcdisagree,labelstrs,true);
%end

%if converttoranks || normaliserange
    tmwuunpaireddisagree = rejtunpaired - rejmwuunpaired,
    disp(['Pairs where unpaired Mann-Whitney U test '...
        'and unpaired t-test disagree:']);
    flaggedpairs(tmwuunpaireddisagree,labelstrs,true);
%end

%if converttoranks || normaliserange
    wilcmwuunpaireddisagree = rejwilc - rejmwuunpaired,
    disp(['Pairs where unpaired Mann-Whitney U test '...
        'and Wilcoxon signed-rank test disagree:']);
    flaggedpairs(wilcmwuunpaireddisagree,labelstrs,true);
%end

if ~converttoranks
    twilcdisagree = rejt - rejwilc,
    disp(['Pairs where paired Wilcoxon signed-rank test '...
        'and paired t-test disagree:']);
    flaggedpairs(twilcdisagree,labelstrs,true);
end

sgnmwuunpaireddisagree = rejsgn - rejmwuunpaired,
disp('Pairs where unpaired Mann-Whitney U test and sign test disagree:');
flaggedpairs(sgnmwuunpaireddisagree,labelstrs,true);

% Display system differences graphically

clf;
imagesc(sign(rejmwuunpaired.*dmedian));

colormap(gca,'gray');
%colorbar;

set(gca,'XTick',1:ncond);
set(gca,'XTickLabel',labelstrs);
set(gca,'YTick',1:ncond);
set(gca,'YTickLabel',labelstrs);

ylabel('Significant preference for condition $y$...','FontSize',bpfontsize);
xlabel(['...over condition $x$, in terms of ' typestr],'FontSize',bpfontsize);

xtickangle(xangle);

set(gca,'FontSize',bpfontsize);
txts = findobj(gca,'Type','text');
set(txts,'FontSize',bpfontsize);
set(txts,'VerticalAlignment','Middle');
clear txts;

curraxis = axis;
hold on;
plot(curraxis(1:2),curraxis(3:4),'k-');
hold off;

saveas(gcf,[outpath filesep filestr '_significance.fig'],'fig');
%saveas(gcf,[outpath filesep filestr '_significance.eps'],'epsc');
saveas(gcf,[outpath filesep filestr '_significance.png'],'png');

fig2fm([outpath filesep filestr '_significance'],'article',...
    'a4paper,british,12pt',17,aspect,9,[],[],5,[],[],true,false);

%prefnoties, % Print empirical preferences after removing ties
%effsizesgn, % Print signed effect size

% The last output is the signed effect size. See this article for details: 
% D. S. Kerby, "The simple difference formula: an approach to teaching 
% nonparametric correlation," Comprehensive Psychology, vol. 3, 2014.

save([outpath filesep filestr '_analysis.mat']); % Save results