function [] = figureS5_CE

% function [] = figureS5_CE
%
% plots panels c and e of supplemental figure S5 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

cutStart = 3000;
winSize = 300;
% to obtain the uncorrected replay active figure, set this to 1

conditionName = {'multistable','replay','random'};

% number of subjects
nSubs = 19;

% interval to analyse
lowerBound = -1000;
upperBound =  0;

uniqueCount = zeros(nSubs,length(conditionName));
eventCount = zeros(nSubs,length(conditionName));
eventSize  = cell(nSubs,length(conditionName));
totalTime  = zeros(nSubs,length(conditionName));
totalPDEvents =  zeros(nSubs,length(conditionName));
pressCount = zeros(nSubs,length(conditionName));
eventTime  = cell(nSubs,length(conditionName));

for conditionCount = 1:length(conditionName)
  
  currentCond = conditionName{conditionCount};
  
  % load button press data
  load(sprintf('./dataStudy2_buttons_%s',conditionName{conditionCount}),'buttonStates','buttonTimes');

  if strcmp(currentCond,'replay')
	load ./dataStudy2_stimPhases_replay stimPhaseStart 
  elseif strcmp(currentCond,'replayPassive')
	load ./dataStudy2_stimPhases_replayPassive stimPhaseStart
  end
	
	
  load (sprintf('dataStudy2_events_%3.3d_%s',winSize,currentCond),'dilationEvent','slope');
  
  % shorthands
  nSubjects = size(dilationEvent,1);
  nBlocks   = size(dilationEvent,2);
  
  % sanity check
  if nSubjects~=nSubs
	error('number of subjects in data file does not match expectation');
  end
  
  for subjectCount = 1 : nSubjects
	
	
	for blockCount = 1 : nBlocks
	
	  currentEvents = dilationEvent{subjectCount,blockCount};
	  currentSlope  = slope{subjectCount,blockCount};
	  
	  % process the button presses
	  currentButtonStates = buttonStates{subjectCount,blockCount};
	  currentButtonTimes  = buttonTimes{subjectCount,blockCount};
	  
	  % find all state changes that correspond to a button press (not a
	  % release), irrespective of the subsequent state
	  idx = find(diff(currentButtonStates)>0)+1;

	  % these are now defined as the transition times (i.e., switch times)
	  transitionTimes = currentButtonTimes(idx);	  
	  
	  
	  totalTime(subjectCount,conditionCount) = ...
		  totalTime(subjectCount,conditionCount) + ...
		  length(currentSlope) - ...
		  cutStart;
	  
	  totalPDEvents(subjectCount,conditionCount) = ...
		  totalPDEvents(subjectCount,conditionCount) + ...
		  sum(currentEvents>cutStart);
	  
	  
	  % for the replay condition, correction with the replay active
	  % condition needs to be applied, and some button presses need to
	  % be removed

	  if strcmp(currentCond,'replay')|strcmp(currentCond,'replayPassive')
	  
		% for the replay condition, some further conditions apply
		currentStimPhaseStart = stimPhaseStart{subjectCount,blockCount};
		removeTransitionIdx = [];

		% loop over all transition times
		for k = 1 : length(transitionTimes)
		
		  % preceding stimulus change
		  ix = max(find(currentStimPhaseStart<transitionTimes(k)));
		  
		  if isempty(ix)
			% there was no preceding stimulus change, remove transition
			removeTransitionIdx(end+1) = k;
		  else
		  
			alreadyRemoved = 0;
			
			if (k~=1)
			  
			  % this is not the first transition time
			  
			  if transitionTimes(k-1)>currentStimPhaseStart(ix)
				
				% the preceding transition followed the same stimulus
				% change (only use the first button press after a
				% stimulus change)
				
				removeTransitionIdx(end+1) = k;
				alreadyRemoved=1;
				
			  end
			  
			end
			
			
			if ((transitionTimes(k)-currentStimPhaseStart(ix))>2000)&(~alreadyRemoved)
			  
			  removeTransitionIdx(end+1) = k;
			  
			end
			
		  end
		  
		end

		% now remove the transitions that violate inclusion criteria
		transitionTimes(removeTransitionIdx)=[];
	  end
		
	  
	  for k = 1 : length(transitionTimes)
		
		if transitionTimes(k)>=cutStart-lowerBound
		  
		  idx1 = transitionTimes(k)+lowerBound;
		  idx2 = transitionTimes(k)+upperBound;
		  
		  eventsInInterval = ...
			  currentEvents((currentEvents>=idx1)&(currentEvents<=idx2));
		  
		  nEvents = numel(eventsInInterval);
		  
		  eventCount(subjectCount,conditionCount) = ...
			  eventCount(subjectCount,conditionCount)+nEvents;
		  
		  uniqueCount(subjectCount,conditionCount) = ...
			  uniqueCount(subjectCount,conditionCount)+(nEvents>0);
		  
		  eventSize{subjectCount,conditionCount}(end+1:end+nEvents) = ...
			  currentSlope(eventsInInterval);
		    eventTime{subjectCount,conditionCount}(end+1:end+nEvents) = ...
			  eventsInInterval-transitionTimes(k);
		  pressCount(subjectCount,conditionCount) = ...
			  pressCount(subjectCount,conditionCount)+1;
		  
		end
		
	  end
	  
	end % blockCount
	
  end % subjectCount

end % conditionCount


medEvent = nan(size(eventSize));

for k1 = 1 : size(medEvent,1)
  for k2 = 1 : size(medEvent,2)

	medEvent(k1,k2) = median(eventSize{k1,k2});
	
  end
end

% compute baseline (in percent)
intervalLength = upperBound-lowerBound;
baseline = 100*intervalLength*totalPDEvents./totalTime;

percWithEvent = 100*uniqueCount./pressCount;


figure(55);
subplot(4,1,2);
cla;
bar([1:3],mean(percWithEvent));
hold on;
errorbar([1:3],mean(percWithEvent),std(percWithEvent)/sqrt(size(percWithEvent,1)),'ko');
ylim([0 100]);
xlabel('condition');
set(gca,'XTick',1:3);
set(gca,'XTickLabel',conditionName);
ylabel('% intervals with event');
title('study 2');


subplot(4,1,4);
cla;
bar([1:3],mean(medEvent));
hold on;
errorbar([1:3],mean(medEvent),std(medEvent)/sqrt(size(medEvent,1)),'ko');
ylim([0 1.4]);
xlabel('condition');
set(gca,'XTick',1:3);
set(gca,'XTickLabel',conditionName);
ylabel('event magnitude');
title('study 2');
