NHP_VisualSearch_Pop-in/CODE/figfunctions/figure_SRT.m

function figure_SRT(fld)

load(fullfile(fld.basedir, 'results','figure_behavior','RTcoll.mat'),'RTcoll');

%% RT distribution
green = [23, 105, 13]./255;
red = [201, 0, 34]./255;
cols = [green; 0.3 0.3 0.3; red];

minRT = 75;
minX2 = 120;

zoomY=[0.35 0.40];

maxRT = 750;
fastRT = 250;

frt = figure;
set(frt,'Position',[0,0,1900,800]);
frt2=figure;
set(frt2,'Position',[0,0,1800,700]);


for mi=1:2
    fastRT = ceil(prctile([RTcoll{mi,1};RTcoll{mi,3}],25));   
    fastRT = 10*ceil(fastRT/10);
    
    p25 = prctile([RTcoll{mi,1};RTcoll{mi,3}],25);
    p50 = median([RTcoll{mi,1};RTcoll{mi,3}]);
    p75 = prctile([RTcoll{mi,1};RTcoll{mi,3}],75);

    figure(frt);
    subplot(2,4,[(mi-1)*2+1 (mi-1)*2+2]); hold on;

    for c = [3 1]
        rts = RTcoll{mi,c};
        rts(rts>maxRT) = [];
        be=minRT:10:maxRT;
        [y,x]=histcounts(rts,be,'Normalization','cdf');
        x = x(1:end-1)+10/2;
        plot(x,y,'LineWidth',2,'Color',cols(c,:))
    end
    set(gca,'xlim',[minRT maxRT],'Box','off');
    ylimits=get(gca,'ylim');
    plot([p25 p25],[0 1.1],'k--');
    plot([p50 p50],[0 1.1],'k-');
    plot([p75 p75],[0 1.1],'k--');
    set(gca,'ylim',[0 1.1]);
    xlabel('RT (ms)');ylabel('CDF'); 
    title(['M' num2str(mi) ': all sacc > 75 ms'])

    %STATS
    % choices
    fprintf('\n==============================\n')
    fprintf(['Stats (RT all) for M' num2str(mi) '\n']);
    fprintf('==============================\n')
    rt1 = RTcoll{mi,1}; rt2=RTcoll{mi,3};

    [p,h,s]=ranksum(rt1,rt2);
    fprintf(['\nM' num2str(mi) ': Wilcoxon for all saccades -----\n'])
    fprintf(['z = ' num2str(s.zval) ', p = ' num2str(p) '\n']);
    fprintf(['TARGET median = ' num2str(median(rt1)) ', IQR = '...
        num2str(iqr(rt1)) '\n']);
    fprintf(['SALIENT DISTRACTOR median = ' num2str(median(rt2)) ', IQR = '...
        num2str(iqr(rt2)) '\n']);

    figure(frt);
    subplot(2,4,(mi-1)*2+5); hold off
    for c = [3 1]
        rts = RTcoll{mi,c};
        rts(rts>fastRT) = [];
        be=minRT:5:fastRT;
        histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
            'probability'); hold all
    end
    set(gca,'xlim',[minX2 fastRT],'ylim',[0 0.29],'Box','off');
    xlabel('RT (ms)');ylabel('Probability'); 
    title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])

    figure(frt2)
    subplot(2,4,4+mi);
    hold on;
    for c = [3 1]
        rts = RTcoll{mi,c};
        rts(rts>fastRT) = [];
        be=minRT:5:fastRT;
        histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
            'probability'); hold all
    end
    set(gca,'xlim',[130 fastRT],'ylim',[0 0.29],'Box','off');
    xlabel('RT (ms)');ylabel('Probability'); 
    title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])

    subplot(2,4,mi);
    hold on;
    nSRT = numel(RTcoll{mi,1})+numel(RTcoll{mi,2})+numel(RTcoll{mi,3});
    for c = [3 1]
        rts = RTcoll{mi,c};
        be=minRT:10:710;
        histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
            'probability'); hold all
    end
    %set(gca,'xlim',[130 715],'ylim',[0 0.29],'Box','off');
    xlabel('RT (ms)');ylabel('Probability'); 
    title(['M' num2str(mi) ': all sacc (norm by choice)'])

    subplot(2,4,2+mi);
    hold on;
    nSRT = numel(RTcoll{mi,1})+numel(RTcoll{mi,2})+numel(RTcoll{mi,3});
    for c = [3 1]
        rts = RTcoll{mi,c};
        be=minRT:10:710;
        H = histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
            'count'); hold all
        set(H,"BinCounts",H.BinCounts./nSRT);
    end
    %set(gca,'xlim',[130 715],'ylim',[0 0.29],'Box','off');
    xlabel('RT (ms)');ylabel('Probability'); 
    title(['M' num2str(mi) ': all sacc (norm by all)'])

    subplot(2,4,6+mi);
    hold on;
    for c = [3 1]
        rts = RTcoll{mi,c};
        rts(rts>fastRT) = [];
        be=minRT:5:fastRT;
        H2 = histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
            'count'); hold all
        set(H2,"BinCounts",H2.BinCounts./nSRT);
    end
    set(gca,'xlim',[130 fastRT],'ylim',[0 0.13],'Box','off');
    xlabel('RT (ms)');ylabel('Probability'); 
    title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])



    %STATS
    % choices
    fprintf('\n==============================\n')
    fprintf(['Stats (25prc fastest RTs) for M' num2str(mi) '\n']);
    fprintf('==============================\n')
    rt1 = RTcoll{mi,1}; rt1=rt1(rt1<fastRT & rt1>minRT);
    rt2 = RTcoll{mi,2}; rt2=rt2(rt2<fastRT & rt2>minRT);
    rt = [RTcoll{mi,1};RTcoll{mi,3}];
    sel = rt<fastRT & rt>minRT ; rt=rt(sel);
    s = categorical([ones(size(RTcoll{mi,1}));...
        2*ones(size(RTcoll{mi,3}))]); s=s(sel);
    
    [p,h,s]=ranksum(rt1,rt2);
    fprintf(['\nM' num2str(mi) ': Wilcoxon for fast saccades -----\n'])
    fprintf(['z = ' num2str(s.zval) ', p = ' num2str(p) '\n']);
    fprintf(['TARGET median = ' num2str(median(rt1)) ', IQR = '...
        num2str(iqr(rt1)) '\n']);
    fprintf(['SALIENT DISTRACTOR median = ' num2str(median(rt2)) ', IQR = '...
        num2str(iqr(rt2)) '\n']);
    fprintf(['Effect size: ' num2str(s.zval./sqrt(length(rt))) '\n']);

    figure(frt);
    subplot(2,4,(mi-1)*2+6);
    hold on;
    plot([p25 p25],[0 1],'k--');
    plot([p50 p50],[0 1],'k-');
    plot([p75 p75],[0 1],'k--');
    for c = [3 1]
        rts = RTcoll{mi,c};
        be=minRT:5:maxRT;
        hold all
        [y,x]=histcounts(rts,be,'Normalization','cdf');
        x = x(1:end-1)+5/2;
        plot(x,y,'LineWidth',2,'Color',cols(c,:))
    end
    set(gca,'xlim',[minX2 p25],'ylim',[0 zoomY(mi)],'Box','off');
    xlabel('RT (ms)');ylabel('CDF'); 
    title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])
end

%% RT_SlidingWindow ====
aRT = [...
    RTcoll{1,1} ones(size(RTcoll{1,1})) ones(size(RTcoll{1,1}));...
    RTcoll{2,1} 2.*ones(size(RTcoll{2,1})) ones(size(RTcoll{2,1}));...
    RTcoll{1,2} ones(size(RTcoll{1,2})) 2.*ones(size(RTcoll{1,2}));...
    RTcoll{2,2} 2.*ones(size(RTcoll{2,2})) 2.*ones(size(RTcoll{2,2}));...
    RTcoll{1,3} ones(size(RTcoll{1,3})) 3.*ones(size(RTcoll{1,3}));...
    RTcoll{2,3} 2.*ones(size(RTcoll{2,3})) 3.*ones(size(RTcoll{2,3}))];

DoPooled = false;
fsw=figure;
set(fsw,'Position',[0,0,900,350]);
if DoPooled
    set(fsw,'Position',[0,0,1350,350]);
end

ws = 20; st = 10; 
minw = 0; maxw = 500;

xmin=125; %minw+st/2; 
xmax=305; %455;

cw = [minw minw+ws];
dotsz = 90;

cRT_T = aRT(aRT(:,3)==1,1);
cRT_SD = aRT(aRT(:,3)==3,1);
cRT_ND = aRT(aRT(:,3)==2,1);

medRT = median(aRT(:,1));
p25 = prctile(aRT(:,1),25);
p75 = prctile(aRT(:,1),75);
mRT = mean(aRT(:,1));

smoothwin = 3;

pSD = []; pND = []; pSDR = [];
while cw(2) < maxw
    TSDratio = ...
        sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./...
        (sum(cRT_T>cw(1) & cRT_T<=cw(2)) + sum(cRT_SD>cw(1) & cRT_SD<=cw(2)));
    TNDratio = ...
        sum(cRT_ND>cw(1) & cRT_ND<=cw(2))./...
        (sum(cRT_T>cw(1) & cRT_T<=cw(2)) + sum(cRT_ND>cw(1) & cRT_ND<=cw(2)));
    SDR = sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./size(aRT,1);
    
    if ~isnan(TSDratio) 
        pSD=[pSD; mean(cw) TSDratio];
    end
    if ~isnan(TNDratio)
        pND=[pND; mean(cw) TNDratio];
    end
    if ~isnan(SDR)
        pSDR=[pSDR; mean(cw) SDR];
    end
    cw=cw+st;
end
if DoPooled
    subplot(1,3,1);hold on;
    plot([0 500],[0 0],'k-'); 
    plot([p25 p25],[-0.1 1.1],'k--');
    plot([medRT medRT],[-0.1 1.1],'k-');
    plot([p75 p75],[-0.1 1.1],'k--');
    scatter(pSDR(:,1),pSDR(:,2),'ko',...
        'MarkerFaceColor','k','MarkerFaceAlpha',0.9,'SizeData',dotsz);
    set(gca,'xlim',[xmin xmax],'ylim',[-0.02 1.02],'Box','off')
    title('Pooled');
    xlabel('RT (ms)'); ylabel('SD choice-ratio')
end

for mi=1:2
    cw = [minw minw+ws];

    cRT_T = aRT(aRT(:,2)==mi & aRT(:,3)==1,1); 
    cRT_SD = aRT(aRT(:,2)==mi & aRT(:,3)==3,1); 
    cRT_ND = aRT(aRT(:,2)==mi & aRT(:,3)==2,1); 
    cRT_A = aRT(aRT(:,2)==mi ,1); 

    medRT = median(aRT(aRT(:,2)==mi,1));
    p25 = prctile(aRT(aRT(:,2)==mi,1),25);
    p75 = prctile(aRT(aRT(:,2)==mi,1),75);
    mRT = mean(aRT(aRT(:,2)==mi,1));

    pSDT = [];pSDND = []; pSDR = [];
    while cw(2) < maxw
        SDTratio = ...
            sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./...
            (sum(cRT_T>cw(1) & cRT_T<=cw(2)) + sum(cRT_SD>cw(1) & cRT_SD<=cw(2)));
        SDNDratio = ...
            sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./...
            (sum(cRT_SD>cw(1) & cRT_SD<=cw(2)) + sum(cRT_ND>cw(1) & cRT_ND<=cw(2))/4);
        SDR = sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./(sum(cRT_A>cw(1) & cRT_A<=cw(2)));
        if ~isnan(SDTratio)
            pSDT=[pSDT; mean(cw) SDTratio];
        end
        if ~isnan(SDNDratio)
            pSDND=[pSDND; mean(cw) SDNDratio];
        end
        if ~isnan(SDR)
            pSDR=[pSDR; mean(cw) SDR];
        end
        cw=cw+st;
    end
    if DoPooled
        subplot(1,3,1+mi);
    else
        subplot(1,2,mi);
    end
    plot([0 500],[0 0],'k-'); hold on;
    plot([p25 p25],[-0.1 1.1],'k--');
    plot([medRT medRT],[-0.1 1.1],'k-');
    plot([p75 p75],[-0.1 1.1],'k--');
    scatter(pSDR(:,1),pSDR(:,2),'ko',...
        'MarkerFaceColor','k','MarkerFaceAlpha',0.8,'SizeData',dotsz);
    set(gca,'xlim',[xmin xmax],'ylim',[-0.02 1.02],'Box','off')
    title(['M' num2str(mi)])
    xlabel('RT (ms)'); ylabel('SD choice-ratio')
    legend({'','p25','median','p75','SD/ALL'});
end

%% stats on SD ratio
% proportions SD choices per octile
nbins=8; 
for m=1:2
    mo(m).cnt=[];
    for i=1:nbins
        STDR(m,i).m = m;
        STDR(m,i).prct = i;
        STDR(m,i).prctval =  prctile(aRT(aRT(:,2)==m,1),i*100/nbins);
        if i == 1
            STDR(m,i).SRT_SD = aRT( aRT(:,3)==3 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval,1);
            STDR(m,i).SRT_A = aRT( aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval,1);
            STDR(m,i).SRT_SDT = aRT( aRT(:,3)~=2 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval,1);
        else
            STDR(m,i).SRT_SD = aRT(aRT(:,3)==3 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval & aRT(:,1)> STDR(m,i-1).prctval,1);
            STDR(m,i).SRT_A = aRT(aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval & aRT(:,1)> STDR(m,i-1).prctval,1);
            STDR(m,i).SRT_SDT = aRT(aRT(:,3)~=2 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval & aRT(:,1)> STDR(m,i-1).prctval,1);
        end
        STDR(m,i).cnt = [length(STDR(m,i).SRT_SD) length(STDR(m,i).SRT_A) length(STDR(m,i).SRT_SDT)];
        mo(m).cnt=[mo(m).cnt; STDR(m,i).cnt];
    end

    % Observed data
    n1 = sum(mo(m).cnt(1,1)); % fastest 12.5%
    N1 = sum(mo(m).cnt(1,2));
    n2 = sum(mo(m).cnt(2:4,1)); % 12.5 - 25%
    N2 = sum(mo(m).cnt(2:4,2));
    % Pooled estimate of proportion
    p0 = (n1+n2) / (N1+N2);
    % Expected counts under H0 (null hypothesis)
    n10 = N1 * p0; n20 = N2 * p0;
    % Chi-square test, by hand
    observed = [n1 N1-n1 n2 N2-n2];
    expected = [n10 N1-n10 n20 N2-n20];
    chi2stat = sum((observed-expected).^2 ./ expected);
    p = 1 - chi2cdf(chi2stat,1);

    fprintf(['== Monkey ' num2str(m) ' pSD ===\n'])
    fprintf(['FASTEST 1/8 vs 1/8 to median\n'])
    fprintf(['CHISQ (1,' num2str(mo(m).cnt(1,2)) ') = ' ...
        num2str(chi2stat) ', p = ' num2str(p) '\n'])
end