function [] = figure3and6

% function [] = figure3and6
%
% plots figures 3 and 6 of
% Grenzebach, J., Wegner, T. G. G., Einhäuser, W., & Bendixen, A. (2021).
% Pupillometry in auditory multistability.
% PLoS ONE, 16(6), e0252370. https://doi.org/10.1371/journal.pone.0252370

conditionName = {'multistable','replay','random'};
figureNum     = [3,6];

% loop over the two studies
for studyNum = 1:2
    
    % variables to colect mean and sem for plotting
    meanPerceptPercentage = nan(length(conditionName),4);
    meanMedianDuration    = nan(1,length(conditionName));
    meanCV                = nan(1,length(conditionName));
    semPerceptPercentage  = nan(length(conditionName),4);
    semMedianDuration     = nan(1,length(conditionName));
    semCV                 = nan(1,length(conditionName));
    
    for conditionCount = 1:length(conditionName)
        
        % load button press data
        load(sprintf('dataStudy%d_buttons_%s',studyNum,conditionName{conditionCount}),'buttonStates','buttonTimes');
        
        % shorthands for the number of subjects and blocks
        nSubjects = size(buttonStates,1);
        nBlocks   = size(buttonStates,2);
        
        % to store the total time a percept was experienced
        totalTime = zeros(nSubjects,4); % 2nd dimension percept type
        
        % to store the relative duration
        perceptPercentage = nan(nSubjects,4); % 2nd dimension percept type
        
        % to store all Durations of exclusive percepts (irrespective of which)
        allDurations = cell(nSubjects,1);
        
        % median
        medianDuration = nan(nSubjects,1);
        
        % CV
        CV = nan(nSubjects,1);
        
        for subjectCount = 1 : nSubjects
            
            for blockCount = 1 : nBlocks
                
                % the button data for the current subject and block
                currentStates = buttonStates{subjectCount,blockCount};
                currentTimes  = buttonTimes{subjectCount,blockCount};
                
                % compute the durations of all phases [in ms]
                durations    = diff(currentTimes);
                
                % exclude the last percept from analysis by setting it a unused
                % state (note that buttonTimes(end) and buttonStates(end)
                currentStates(end-1) = nan;
                
                % restrict the analysis to periods in which a button was
                % exclusively pressed
                
                % the total time each percept was experienced
                for k = 1 : 4
                    totalTime(subjectCount,k) = ....
                        totalTime(subjectCount,k) + ...
                        sum(durations(find(currentStates==2^(k-1))));
                    % (2^(k-1) is the code for the respective exclusive button)
                end
                
                
                % all dominance durations of exclusive percepts
                allDurations{subjectCount} = [allDurations{subjectCount},...
                    durations(find(currentStates==1)),...
                    durations(find(currentStates==2)),...
                    durations(find(currentStates==4)),...
                    durations(find(currentStates==8))];
                
            end
            
            
            % now compute the relevant stats for this subject
            perceptPercentage(subjectCount,:) = 100*totalTime(subjectCount,:)/sum(totalTime(subjectCount,:),2);
            medianDuration(subjectCount)     = ...
                median(allDurations{subjectCount})/1000; % convert to seconds
            CV(subjectCount)                 = std(allDurations{subjectCount})/mean(allDurations{subjectCount});
            
            
        end
        
        % compute mean and sem over subjects for later plotting
        meanPerceptPercentage(conditionCount,:) = mean(perceptPercentage);
        meanMedianDuration(conditionCount)     = mean(medianDuration);
        meanCV(conditionCount)                 = mean(CV);
        
        semPerceptPercentage(conditionCount,:)  = std(perceptPercentage)/sqrt(size(perceptPercentage,1));
        semMedianDuration(conditionCount)      = std(medianDuration)/sqrt(size(medianDuration,1));
        semCV(conditionCount)                   = std(CV)/sqrt(size(CV,1));
        
        
    end
    
    
    % do the plotting for respective study
    figure(figureNum(studyNum));
    clf;
    
    subplot(1,3,1);
    bar([1:3],meanPerceptPercentage,'stacked');
    ylabel('% exclusive percept');
    xlabel('condition');
    set(gca,'XTick',[1:3]);
    set(gca,'XTickLabel',conditionName);
    xtickangle(90);
    
    subplot(1,3,2);
    bar([1:3],meanMedianDuration);
    hold on;
    errorbar([1:3],meanMedianDuration,semMedianDuration,'ko');
    ylabel('median dominance duration [s]');
    xlabel('condition');
    set(gca,'XTick',[1:3]);
    set(gca,'XTickLabel',conditionName);
    xtickangle(90);
    
    subplot(1,3,3);
    bar([1:3],meanCV);
    hold on;
    errorbar([1:3],meanCV,semCV,'ko');
    ylabel('CV');
    xlabel('condition');
    set(gca,'XTick',[1:3]);
    set(gca,'XTickLabel',conditionName);
    xtickangle(90);
    
end
