close all; clear; clc

%-----------------------------------------------
% Specify the name of the COMPASS root directory
dataPath='/path/to/COMPASS/repository/';
%-----------------------------------------------
addpath(genpath([dataPath 'Scripts/']))

[~, txtData, allData]=xlsread([dataPath 'Data/COMPASS_SceneInfo.xlsx']); 
header=txtData(1,:); sceneInfo=cell2table(allData(2:end,:),'VariableNames',header);

images=sceneInfo.ImageFile;

%% Do any of the images change category when all trials are included (i.e., no RT cutoffs)?

[numData, ~, ~]=xlsread([dataPath 'Data/COMPASS_ValenceRatings.xlsx']); 
valRatings=numData(:,2:end);

[numData, ~, ~]=xlsread([dataPath 'Data/COMPASS_ArousalRatings.xlsx']); 
arousRatings=numData(:,2:end);

[numData, ~, ~]=xlsread([dataPath 'Data/COMPASS_ValenceRTs.xlsx']); 
valRTs=numData(:,2:end);

[numData, ~, ~]=xlsread([dataPath 'Data/COMPASS_ArousalRTs.xlsx']); 
arousRTs=numData(:,2:end);

exclValTrials=valRTs>4000 | valRTs<150;
exclArousTrials=arousRTs>4000 | arousRTs<150;

fprintf('Excluded %d (%.1f%%) valence ratings \n',sum(sum(exclValTrials)),...
    sum(sum(exclValTrials))/numel(exclValTrials)*100) 
fprintf('Excluded %d (%.1f%%) arousal ratings \n\n',sum(sum(exclArousTrials)),...
    sum(sum(exclArousTrials))/numel(exclArousTrials)*100) 

valRatingsExcl=valRatings;
valRatingsExcl(logical(exclValTrials))=NaN; 
arousRatingsExcl=arousRatings;
arousRatingsExcl(logical(exclArousTrials))=NaN; 

clearvars meanRatings
meanRatings.Image=images;meanRatings=struct2table(meanRatings);
meanRatings.ImageName=sceneInfo.ImageName;
meanRatings.Name=sceneInfo.Name;
meanRatings.ValenceMean=nanmean(valRatings)';
meanRatings.ValenceSD=nanstd(valRatings)';	
meanRatings.ArousalMean=nanmean(arousRatings)';	
meanRatings.ArousalSD=nanstd(arousRatings)';	
meanRatings.NumValRatings=sum(~isnan(valRatings))';	
meanRatings.NumArousRatings=sum(~isnan(arousRatings))';

clearvars meanRatingsExcl
meanRatingsExcl.Image=images;meanRatingsExcl=struct2table(meanRatingsExcl);
meanRatingsExcl.ImageName=sceneInfo.ImageName;
meanRatingsExcl.Name=sceneInfo.Name;
meanRatingsExcl.ValenceMean=nanmean(valRatingsExcl)';
meanRatingsExcl.ValenceSD=nanstd(valRatingsExcl)';	
meanRatingsExcl.ArousalMean=nanmean(arousRatingsExcl)';	
meanRatingsExcl.ArousalSD=nanstd(arousRatingsExcl)';	
meanRatingsExcl.NumValRatings=sum(~isnan(valRatingsExcl))';	
meanRatingsExcl.NumArousRatings=sum(~isnan(arousRatingsExcl))';

%% Plot changes in ratings from before to after trial exclusion (Figure S1).

lineWidth=1.5;fontSize=15;
screenSize=get(0,'ScreenSize');
socColor=[1 .6 .1];nonSocColor=[.1 .6 1];
socInd=strcmp(sceneInfo.SocialContent,'Social');

figure('color','w','units','normalized','Position',screenSize);hold on
l1=line([1 9],[5 5],'LineWidth',lineWidth,'Color',[.4 .4 .4]);
l2=line([5 5],[1 9],'LineWidth',lineWidth,'Color',[.4 .4 .4]);

for ii=1:numel(images)
    if socInd(ii)==1;pickColor=socColor; dispName='Social'; 
    elseif socInd(ii)==0;pickColor=nonSocColor; end
    plot([meanRatings.ValenceMean(ii) meanRatingsExcl.ValenceMean(ii)],...
        [meanRatings.ArousalMean(ii) meanRatingsExcl.ArousalMean(ii)],...
        'LineWidth',lineWidth*2,'Color',pickColor,'DisplayName',dispName); hold on
end

s1=scatter(meanRatings.ValenceMean,meanRatings.ArousalMean,...
    'Marker','.','DisplayName','All trials');
s2=scatter(meanRatingsExcl.ValenceMean,meanRatingsExcl.ArousalMean,...
    'Marker','.','MarkerEdgeColor','k','DisplayName','Trials excluded');

xlim([1 9]);ylim([1 9]); box on;grid on; axis square
set(gca,'XTick',[1:9],'YTick',[1:9],'LineWidth',lineWidth,'FontSize',fontSize)
xlabel('Valence');ylabel('Arousal'); 
title('Change in Mean Valence & Arousal Ratings After Trial Exclusion');l=legend([s1 s2]);

%% Calculate image complexity.
% Text S3.

inFolder=[dataPath 'Materials/Stimuli/'];
sceneInfo.Entropy=nan(size(sceneInfo,1),1);
sceneInfo.jpegSize=nan(size(sceneInfo,1),1);

% Calculate the entropy of each image.
for ii=1:height(sceneInfo)   
    scene=imread([inFolder sceneInfo.ImageName{ii} '.jpg']);
    grayScene=rgb2gray(scene);
    [counts,binLocations] = imhist(grayScene);
    p=counts./numel(grayScene);
    sceneInfo.Entropy(ii)=-nansum(p.*log2(p));
end

% Determine JPEG size of each image.
cd(inFolder)
jpegInfo=dir('*.jpg');
jpegInfo=struct2table(jpegInfo);
for ii=1:height(sceneInfo)   
    ind=strcmp([sceneInfo.ImageName{ii} '.jpg'],jpegInfo.name);
    sceneInfo.jpegSize(ii)=jpegInfo.bytes(ind);
end

%% Compare the complexity of social and nonsocial scenes using the TOST procedure.
% Text S3.

% Select the lower and upper equivalence bounds (if the effect size is
% significantly within these bounds, the two samples are considered equivalent)

dL=-.35; % Based on 80% power and alpha=0.05 given the sample size (Lakens, 2017)
dU=.35;

% Select (uncomment) a measure of complexity 

%%% JPEG size
sample1=sceneInfo.jpegSize(strcmp(sceneInfo.SocialContent,'Social'));
sample2=sceneInfo.jpegSize(strcmp(sceneInfo.SocialContent,'Nonsocial'));

%%% OR entropy
% sample1=sceneInfo.Entropy(strcmp(sceneInfo.SocialContent,'Social'));
% sample2=sceneInfo.Entropy(strcmp(sceneInfo.SocialContent,'Nonsocial'));

% Calculate NHST results & Cohen's d
[~,p,~,stats]=ttest2(sample1,sample2);
d=cohensD(sample1,sample2);
fprintf('NHST: t(%d)=%.3f, p=%.4f, d=%.3f\n',stats.df,stats.tstat,p,d)

% Calculate TOST results
results=tost(sample1,sample2,dL,dU);
if results.pL>results.pU 
    tstat=results.tL; pval=results.pL; 
else
    tstat=results.tU; pval=results.pU; 
end
fprintf('TOST: t(%d)=%.3f, p=%.4f\n',results.df,tstat,pval)

%% Test equivalence of low-level visual features between different affect categories.
% Text S3.

clc

% Select the equivalence bounds
dL=-.5;
dU=.5;

% Select (uncomment) the low-level visual feature of interest
meas='Contrast'
% meas='Luminance'

% Social Content
sample1=sceneInfo.(meas)(strcmp(sceneInfo.SocialContent,'Social'));
sample2=sceneInfo.(meas)(strcmp(sceneInfo.SocialContent,'Nonsocial'));

% Calculate NHST results & Cohen's d
[~,p,~,stats]=ttest2(sample1,sample2);d=cohensD(sample1,sample2);
fprintf('Social Content NHST: t(%d)=%.3f, p=%.4f, d=%.3f\n',stats.df,stats.tstat,p,d)

% Calculate TOST results
results=tost(sample1,sample2,dL,dU);
if results.pL>results.pU;tstat=results.tL; pval=results.pL; 
else tstat=results.tU; pval=results.pU; end
fprintf('Social Content TOST: t(%d)=%.3f, p=%.4f\n\n',results.df,tstat,pval)

% Arousal
ind11=strcmp(sceneInfo.Arousal,'High'); 
ind12=strcmp(sceneInfo.Arousal,'Mid') & strcmp(sceneInfo.Valence,'Neutral');
ind1=logical(ind11+ind12);sample1=sceneInfo.(meas)(ind1);
ind21=strcmp(sceneInfo.Arousal,'Low');
ind22=strcmp(sceneInfo.Arousal,'Mid') & ~strcmp(sceneInfo.Valence,'Neutral');
ind2=logical(ind21+ind22);sample2=sceneInfo.(meas)(ind2);

% Calculate NHST results & Cohen's d
[~,p,~,stats]=ttest2(sample1,sample2);d=cohensD(sample1,sample2);
fprintf('Arousal NHST: t(%d)=%.3f, p=%.4f, d=%.3f\n',stats.df,stats.tstat,p,d)

% Calculate TOST results
results=tost(sample1,sample2,dL,dU);
if results.pL>results.pU; tstat=results.tL; pval=results.pL; 
else tstat=results.tU; pval=results.pU; end
fprintf('Arousal TOST: t(%d)=%.3f, p=%.4f\n\n',results.df,tstat,pval)

% Valence
sample1=sceneInfo.(meas)(strcmp(sceneInfo.Valence,'Positive'));
sample2=sceneInfo.(meas)(strcmp(sceneInfo.Valence,'Neutral'));
sample3=sceneInfo.(meas)(strcmp(sceneInfo.Valence,'Negative'));

% Calculate NHST results & Cohen's d
[~,p,~,stats]=ttest2(sample1,sample2);d=cohensD(sample1,sample2);
fprintf('Pos vs Neut NHST: t(%d)=%.3f, p=%.4f, d=%.3f\n',stats.df,stats.tstat,p,d)

% Calculate TOST results
results=tost(sample1,sample2,dL,dU);
if results.pL>results.pU; tstat=results.tL; pval=results.pL; 
else tstat=results.tU; pval=results.pU; end
fprintf('Pos vs Neut TOST: t(%d)=%.3f, p=%.4f\n\n',results.df,tstat,pval)

% Calculate NHST results & Cohen's d
[~,p,~,stats]=ttest2(sample3,sample2);
d=cohensD(sample3,sample2);
fprintf('Neg vs Neut NHST: t(%d)=%.3f, p=%.4f, d=%.3f\n',stats.df,stats.tstat,p,d)

% Calculate TOST results
results=tost(sample3,sample2,dL,dU);
if results.pL>results.pU; tstat=results.tL; pval=results.pL; 
else tstat=results.tU; pval=results.pU; end
fprintf('Neg vs Neut TOST: t(%d)=%.3f, p=%.4f\n\n',results.df,tstat,pval)

% Calculate NHST results & Cohen's d
[~,p,~,stats]=ttest2(sample3,sample1);d=cohensD(sample3,sample1);
fprintf('Neg vs Pos NHST: t(%d)=%.3f, p=%.4f, d=%.3f\n',stats.df,stats.tstat,p,d)

% Calculate TOST results
results=tost(sample3,sample1,dL,dU);
if results.pL>results.pU; tstat=results.tL; pval=results.pL; 
else tstat=results.tU; pval=results.pU; end
fprintf('Neg vs Pos TOST: t(%d)=%.3f, p=%.4f\n',results.df,tstat,pval)
