function [dataOut,dataOutI] = findConS(dataIn,seq,nTrials,nConS,exactNConSI)
% IDENTIFY consecutive standard responses
% This function identifies consecutive standard responses in the oddball
% paradigm. The number of responses is defined by nConS

% define some variables
pSmpl = size(dataIn,1);
nBlcks = size(dataIn,3);
nFiles = size(dataIn,4);
nFilt = size(dataIn,5);

dataOut = cell(nBlcks,nFiles); % preallocate
dataOutI = cell(nBlcks,nFiles); % preallocate
for f = 1:nFiles
    for b = 1:nBlcks
        % create std sequences
        seqStd = ~seq(1:nTrials(b,f),b,f)';

        % create vector that contains all clusters of consecutive standards
        % code from forum: replace ones by number of consecutive repetitions
        switch exactNConSI
            case 0
                vecClI = double(seqStd);
            case 1
                d = [1,diff(seqStd)~=0,1]; % TRUE if values change
                n = diff(find(d)); % Number of repetitions
                X = repelem(n,n);
                vecClI = X.*seqStd; % accept only clusters of consecutive standards and not deviants
                vecClI(vecClI<nConS) = 0; % replace numbers that are smaller than criterion (less than X consecutive repetitions) by 0
                vecClI(vecClI>=nConS) = 1; % replace numbers that are equal or larger than criterion (X consecutive repetitions or more) by 1
        end

        % count appearances of consecutive standards in each accepted
        % cluster (code from forum)
        B = vecClI; 
        d2 = [true,diff(B)~=0];
        r = find(d2);
        B(~d2) = 1;
        B(d2) = [1,1-diff(r)];
        B(1) = 1;
        B = cumsum(B);

        % remove standards that appear in a cluster at a higher position
        % than nConS
        B(B<=nConS) = 1; % replace numbers that are equal or smaller than criterion by 1
        B(B>nConS) = 0; % replace numbers that are larger than criterion by 0
        vecClI = logical(vecClI.*B); % accept only the first X standard responses of clusters of consecutive standards that are equal or larger than X

        % identify number and size of accepted clusters
        bCl = diff([0,vecClI,0]); % vector that marks borders of clusters
        nCl = sum(bCl==1); % number of accepted clusters (1 marks the starting point of a cluster)
        sCl = find(bCl==1); % starting points of clusters
        eCl = find(bCl==-1)-1; % ending points of clusters
        lCl = eCl-sCl+1; % size of clusters

        % produce matrices that contain all accepted vectors
        dataOutI{b,f} = zeros(nConS,nCl); % preallocate
        dataOut{b,f} = zeros(pSmpl,nConS,nCl,nFilt); % preallocate
        for n = 1:nCl % once for each cluster
            dataOutI{b,f}(1:lCl(n),n) = vecClI(sCl(n):eCl(n)); % replace 0s by 1s for each std in the given cluster
            dataOut{b,f}(:,1:lCl(n),n,:) = dataIn(:,sCl(n):eCl(n),b,f,:); % create new variable that contains all std responses of the given cluster
        end
        dataOutI{b,f} = logical(dataOutI{b,f}); % convert to logical
    end
end