NHP_VisualSearch_Pop-in/CODE/figfunctions/figure_neuropriming_prepare.m

function figure_neuropriming_prepare(fld,monkey)

fprintf('\n======================================================\n');
fprintf(['-- Running neuropriming for ' monkey ' --\n']);
fprintf('======================================================\n');

% settings 
snrthres = 2.5;
daythres = 3;
cols = [1 0 0;...
    0 0 0;...
    0 1 0];
do_cleantraces = true; % if true we throw out artefact traces

%% find which data to include 
datainfo = session_info();

if strcmp(monkey,'M1')
    info = datainfo(1);
else 
    info = datainfo(2);
end
sessions = info.dates;

MP = get_mp(fld,monkey,sessions);

%% set SNR threshold and find channels and sessions
% limit MP to snr threshold
incl = MP.snr>snrthres;
MP_ = MP; % copy for backup

% only include channels that have at least 5 sessions
MP = MP(incl,:);
[chid, ch_list] = findgroups(MP.chan_id);
session_per_chan = splitapply(@length,MP.session_id,chid);

% make a list of channels that we will include
chans = ch_list(session_per_chan>daythres);

%% load and concatenate the data we need
datadir = fld.procdatadir;

% first loop all channels 
numtrls = [];
targmod = []; targmag = [];
distmod = []; distmag = [];
tracesLUT = [];
for ch = chans'

    chidx = find(chans==ch);
    disp(['analysing channel ' num2str(ch) ', ' ...
        num2str(chidx) '/' num2str(size(chans,1))]);
    
    incl_sessions = MP.chan_id==ch & MP.snr>snrthres;
    cur_mp = MP(incl_sessions,:);
    numtrls(chidx) = sum(cellfun(@sum,cur_mp.goodtrials));
    
    % look at each session for each channel individually 
    env = []; clut = [];
    for sid = 1:size(cur_mp,1)

        % session info
        cdate = cur_mp.date{sid};
        block = cur_mp.block(sid);
        
        % load the data 
        session_id = [monkey '_' cdate '_' num2str(block)];
        sessiondir = fullfile(datadir, monkey, session_id);
        channeldir = fullfile(sessiondir, 'individual_channels');
        chanfile = fullfile(channeldir, ['Env_preprocessed_ch' ...
            num2str(ch) '_' session_id '.mat']);
        if ~exist(chanfile)
            continue
        end
        
        % load the data for this session for this channel
        load(chanfile); % loads 'data'
        load(fullfile(sessiondir, [session_id '_LUT.mat'])); % loads 'lut'
        load(fullfile(sessiondir, [session_id '.mat']),'env_t'); 
        % this is overwritten every time but it's always the same so who cares
        
        if do_cleantraces % throw out artefact traces
            prewin = env_t > -0.1 & env_t < 0;
            postwin = env_t > 0 & env_t < 0.1;
            trace = mean(data');

            sd_pre = std(trace(prewin));
            m_post = mean(trace(postwin));

            if m_post < 3*sd_pre
                data = nan(size(data));
            end
        end

        env = [env; data'];
        clut = [clut; lut];
    end
    
    lut = clut;
    
    % add a shapeswitch column, because apparently no one ever looked at
    % this before
    shapeSwitch = zeros(size(lut,1),1);
    shapeSwitch(1) = nan;
    shapeSwitch(abs(lut.TarShape - [nan; lut.TarShape(1:end-1)]) > 0) = 1;
    lut = addvars(lut, shapeSwitch,'after','colSwitch');  
    
    targetInRF = lut.TarPos==MP.rfpos(ch);
    nontargInRF = lut.TarPos~=MP.rfpos(ch) & lut.DistPos~=MP.rfpos(ch);
    distInRF = lut.DistPos==MP.rfpos(ch);
    
    goodtrials = [cur_mp.sample_wise_goodtrials{:}]';
    
    for shapeswitch = [0 1]
        incl = lut.TarChoice==lut.TarPos & ...
            targetInRF & goodtrials & lut.shapeSwitch==shapeswitch;
        mn1 = nanmean(env(incl,:));

        % same, when a non-target is in the RF
        incl = lut.TarChoice==lut.TarPos & ...
            nontargInRF & goodtrials & lut.shapeSwitch==shapeswitch;
        mn2 = nanmean(env(incl,:));

        % same, for distractor
        incl = lut.TarChoice==lut.TarPos & ...
            distInRF & goodtrials & lut.shapeSwitch==shapeswitch;
        mn3 = nanmean(env(incl,:));

        % store the difference in traces
        targmod(end+1,:) = mn1-mn2;
        distmod(end+1,:) = mn3-mn2;
        targmag(end+1,:) = mn1;
        distmag(end+1,:) = mn3;
        tracesLUT(end+1,:) = [ch,shapeswitch,nan,sum(incl)];
    end
    
    for colorswitch = [0 1]
        incl = lut.TarChoice==lut.TarPos & ...
            targetInRF & goodtrials & lut.colSwitch==colorswitch;
        mn1 = nanmean(env(incl,:));

        % same, when a non-target is in the RF
        incl = lut.TarChoice==lut.TarPos & ...
            nontargInRF & goodtrials & lut.colSwitch==colorswitch;
        mn2 = nanmean(env(incl,:));

        % same, for distractor
        incl = lut.TarChoice==lut.TarPos & ...
            distInRF & goodtrials & lut.colSwitch==colorswitch;
        mn3 = nanmean(env(incl,:));

        % store the difference in traces
        targmod(end+1,:) = mn1-mn2;
        distmod(end+1,:) = mn3-mn2;
        targmag(end+1,:) = mn1;
        distmag(end+1,:) = mn3;
        tracesLUT(end+1,:) = [ch,nan,colorswitch,sum(incl)];
    end
end

tracesLUT = array2table(tracesLUT,'VariableNames',{'channel',...
    'shapeswitch','colorswitch','numtrials'});

% save output to make a combined plot of all channels
savedir = fullfile(fld.basedir,'results','figure_neuropriming');
save(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) ...
    '_mindays' num2str(daythres) '.mat']));