NHP_VisualSearch_Pop-in/CODE/figfunctions/figure_priming_color.m

function figure_priming_color(fld)

fprintf('\n======================================================\n');
fprintf('-- Creating Supp Figure priming --\n');
fprintf('-- COLOR SWAP & REWARD SIZE --\n');
fprintf('======================================================\n');

%% get sessions
datainfo = session_info();
datadir = fld.procdatadir;
savedir = fullfile(fld.basedir, 'results','figure_priming');

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

%% plots
% load traces
monkeys = {'M1','M2'};
f1 = figure; set(gcf,'Position',[100   100   600   1000]);
for mi = 1:2
    monkey = monkeys{mi};
    
    info = datainfo(mi);
    sessions = info.dates;
    
    clut = [];
    for sid = 1:length(sessions)
        session = sessions{sid};
        % session info
        q = strsplit(session,'_');
        cdate = q{2};
        block = q{3}(4);
        
        % load the data 
        session_id = [monkey '_' cdate '_' num2str(block)];
        sessiondir = fullfile(datadir, monkey, session_id);
        
        % load the data for this session for this channel
        load(fullfile(sessiondir, [session_id '_LUT.mat'])); % loads 'lut'

        clut = [clut; lut];
    end
    
    lut = clut;
    
    % calculate mean across all sessions
    [meanperf, condinfo] = calculate_perf(lut);
    
    % calculate mean across individual sessions
    incl_sessions = unique(lut.session);
    dayperf = [];
    for sid = 1:length(incl_sessions)
        
        clut = lut(strcmp(lut.session,incl_sessions{sid}),:);
        
        dayperf(end+1,:) = calculate_perf(clut);
    end
    
    % calculate sem based on individual sessions
    perfstd = std(dayperf)/sqrt(size(dayperf,1));
    
    yellow = [247 148 29]./255;
    blue = [27 117 188]./255;
    cols = [blue; yellow];
    
    subplot(2,2,mi)
    gap = 0.05; 
    for colswitch = [0 1]
        i = condinfo(:,1)==colswitch;
        mn = meanperf(i);
        st = perfstd(i); 
        plot([1 2]+gap*colswitch*-1,mn,'Color',cols(colswitch+1,:)); hold all;
        errorbar([1 2]+gap*colswitch*-1,mn,st,'o','Color',cols(colswitch+1,:),'MarkerFaceColor', cols(colswitch+1,:));
    end
    xlim([0.8 2.2]);
    set(gca,'XTick',[1 2],'XTickLabel',{'Low Reward','High Reward'},...
        'ylim',[0.4 1]);
    title(monkey);
    ylabel('Performance Ratio');
    legend({'no cswap','','cswap',''});

    
    %% stats 
    % let's start with a two-way anova. we have two conditions, each with
    % two levels (reward, color switch). in matlab, the 2d matrix supplying
    % the data needs to be a bit weird. columns contain the  first
    % condition (reward value in our case). rows contains the second
    % condition (color switch), and it concatenates the
    % repetitions for each level of that condition.
    nocolswitch = dayperf(:,1:2);
    colswitch = dayperf(:,3:4);
    Y = [nocolswitch; colswitch];
    [p,tbl,stats] = anova2(Y,size(colswitch,1),'off');
    tbl
    disp('-------performance-------');
    disp([monkey ', interaction colswap and rewval, p=' num2str(p(3))]);
    disp([monkey ', main effect for color swap, p=' num2str(p(2))]);
    disp([monkey ', main effect for reward value of previous trial, p=' num2str(p(1))]);
    
    %% reaction time
    % select trials that have: 
    % - current trial correct
    % - previous trial correct
    % - trial immediately followed previous trial
    % - color switch/ color same
    % - reward high/ reward low (on previous trial)
    rt = []; rtst = [];
    Y_rt = []; Y_colswitch = []; Y_rewval = [];
    for colswitch = [0 1]
        for rewval = [1 2]
            inclTrls = strcmp(lut.Status,'Target') & strcmp(lut.prevStatus,'Target') & lut.colSwitch==colswitch & lut.rewardVal==rewval & lut.immediateFollow;
            
            rt(end+1) = mean(lut.reactTime2(inclTrls));
            rtst(end+1) = std(lut.reactTime2(inclTrls))./sqrt(sum(inclTrls));
            
            % collect in a vector for statistics
            Y_rt = [Y_rt; lut.reactTime2(inclTrls)];
            Y_colswitch = [Y_colswitch; repmat(colswitch,sum(inclTrls),1)];
            Y_rewval = [Y_rewval; repmat(rewval,sum(inclTrls),1)]; 
            
        end
    end
    
    subplot(2,2,mi+2);
    for colswitch = [0 1]
        i = condinfo(:,1)==colswitch;
        mn = rt(i);
        st = rtst(i); 
        plot([1 2]+gap*colswitch*-1,mn,'Color',cols(colswitch+1,:)); hold all;
        errorbar([1 2]+gap*colswitch*-1,mn,st,'o','Color',cols(colswitch+1,:),'MarkerFaceColor', cols(colswitch+1,:));
    end
    xlim([0.8 2.2]);
    set(gca,'XTick',[1 2],'XTickLabel',{'Low Reward','High Reward'},...
        'ylim',[180 280]);
    title(monkey);
    ylabel('Reaction Time (ms)');
    legend({'no colorswap','','colorswap',''})
    
    %% stats
    % in this case, we have unbalanced design, because we take all the
    % trials as repetitions, and they are not equal amounts for each
    % condition, so we use anovan. here, Y is a vector, and the second
    % input gives the group number. 
    [p,tbl,stats] = anovan(Y_rt,{Y_colswitch, Y_rewval},'model','interaction','varnames',{'Color swap','Reward value of prevtrial'},'display','off');
    tbl
    disp('-------reaction time-------');
    disp([monkey ', interaction colswap and rewval, p=' num2str(p(3))]);
    disp([monkey ', main effect for reward value of previous trial, p=' num2str(p(2))]);
    disp([monkey ', main effect for color swap, p=' num2str(p(1))]);
    
end
suptitle('Behavioral effects of color swap & reward')
savefig(f1,fullfile(savedir, 'figure_priming_color'));
print(f1,fullfile(savedir, 'figure_priming_color'),'-dpng');
end


function [perf, cond] = calculate_perf(lut) 
    % select trials that have: 
    % - are not aborted
    % - current trial correct
    % - previous trial correct
    % - trial immediately followed previous trial
    % - color switch/ color same
    % - reward high/ reward low (on previous trial)
    perf = []; cond = [];
    for colswitch = [0 1]
        for rewval = [1 2]
            inclTrls = strcmp(lut.Status,'Aborted')==0 & strcmp(lut.prevStatus,'Target') & lut.colSwitch==colswitch & lut.rewardVal==rewval & lut.immediateFollow;
            curlut = lut(inclTrls,:);
            
            perf(end+1) = sum(strcmp(curlut.Status,'Target'))/size(curlut,1); 
            cond(end+1,:) = [colswitch rewval];
        end
    end
end