function [meanPupil, meanMissing] = study2_replayPassive


% function [meanPupil, meanMissing] = study2_replayPassive
%
% computes the curves of the "replay passive" condition used for
% correction of the replay active curves


% data to be used left and right of 0, in ms [note that the extra 2s
% to the right are to allow shifting the curve before subtraction]
maxDataLeft    = 2000;
maxDataRight   = 4000;
dataLength = maxDataLeft+maxDataRight+1;

% load the stimulus change data
load dataStudy2_stimPhases_replayPassive stimPhaseStart

% load the pupil data
load(sprintf('dataStudy2_pupil_replayPassive'),'pupilSize','xEye','yEye','tEye');

% number of subjects
nSubjects   = size(pupilSize,1);
nBlocks     = size(pupilSize,2);

% to collect data by subject
meanPupil   = nan(nSubjects,dataLength);
meanMissing = nan(nSubjects,dataLength);

for subjectCount = 1 : nSubjects
    
    % to collect data per subject
    allPupil    = nan(0,dataLength);
    allMissing  = nan(0,dataLength);
    
    for blockCount = 1 : nBlocks
        
        
        % short hand notations for data of the current block
        currentPupilSize = pupilSize{subjectCount,blockCount};
        currentXEye      = xEye{subjectCount,blockCount};
        currentYEye      = yEye{subjectCount,blockCount};
        
        % remove data outside a central 10deg square
        currentPupilSize((currentXEye>1285)|(currentXEye<635))=nan;
        currentPupilSize((currentYEye>865) |(currentYEye<215))=nan;
        clear currentXEye currentYEye
        
        
        % remove the first 3s from analysis
        currentPupilSize(1:3000)=nan;
        
        
        % with this data, do the z-normalization
        pupilZ = (currentPupilSize-nanmean(currentPupilSize))/nanstd(currentPupilSize);
        clear currentPupilSize;
        
        
        % unlike in the active conditions, we do not care about button
        % presses, but look for stimulus changes
        
        % get the transitions, the previous and the following one
        transitionTimes = stimPhaseStart{subjectCount,blockCount}(2:end); % the first onset is not a transition
        
        % time of preceding and next transition
        prevTimes       = [0, transitionTimes(1:end-1)]; % at onset, this is trial start
        nextTimes       = [transitionTimes(2:end), numel(pupilZ)]; % at offset, this is trial end
        
        
        
        
        % define the usable time window for each transition
        startIdx = nan(1,length(transitionTimes));
        endIdx   = nan(1,length(transitionTimes));
        idxLeft = nan(1,length(transitionTimes));
        idxRight = nan(1,length(transitionTimes));
        
        for k = 1 : length(transitionTimes)
            
            % start and end of analysed epoch (-2s:2s)
            startIdx(k) = max(round((prevTimes(k)+transitionTimes(k))/2),transitionTimes(k)-maxDataLeft);
            endIdx(k)   = min(round((nextTimes(k)+transitionTimes(k))/2),transitionTimes(k)+maxDataRight);
            
            % index in aggregated trace (1,4001 for full length, centred
            % around 2001 otherwise)
            idxLeft(k)  = 1+maxDataLeft+startIdx(k)-transitionTimes(k);
            idxRight(k) = endIdx(k)-transitionTimes(k)+maxDataLeft+1; % left is correct here, as it as about the shift relative to t=0, which corresponds to maxDataLeft+1;
            
        end
        
        % now generate the traces
        localPupil   = nan(length(transitionTimes),dataLength);
        localMissing = nan(length(transitionTimes),dataLength);
        
        
        for k = 1 : length(transitionTimes)
            
            localPupil(k,idxLeft(k):idxRight(k))    = pupilZ(startIdx(k):endIdx(k));
            localMissing(k,1:dataLength) = isnan(localPupil(k,:));
            
        end
        
        % aggregate data across blocks
        allPupil    = [allPupil;localPupil];
        allMissing  = [allMissing;localMissing];
        
    end
    
    
    % mean per subject
    meanPupil(subjectCount,:)   = nanmean(allPupil);
    
    % only use traces for the computation of missing data points over
    % time that include at least one data point
    idxGood = find(sum(allMissing')~=dataLength);
    meanMissing(subjectCount,:) = ...
        nanmean(allMissing(idxGood,:),1);
    
    % baseline correction
    baseline = mean(meanPupil(subjectCount,1:200));
    meanPupil(subjectCount,:) = meanPupil(subjectCount,:) - baseline;
    
end

