NHP_VisualSearch_Pop-in/CODE/functions/analyse_individual_sessions.m

function analyse_individual_sessions(fld, monkey, sessions)

datadir = fullfile(fld.basedir,'Datasets', monkey ,'MAT_eye');
logdir = fullfile(fld.basedir,'Datasets', monkey ,'log');

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

% loop sessions
for sid = 1:length(sessions)
    disp(['analysing session ' num2str(sid) ' out of ' num2str(length(sessions))]);
    
    % session info
    session = sessions{sid};
    q = strsplit(session,'_');
    cdate = q{2};
    
    % find files
    FN = dir(fullfile(datadir, [monkey cdate '-block*']));
    % we will have 2 envs for each block, and at least 1 block for each
    % date. so, if length(FN)>2, we have more than 1 block. 
    available_blocks = [];
    for f = 1:length(FN)
        available_blocks(end+1) = str2num(FN(f).name(end-9));
    end
    blocks = unique(available_blocks);
    
    % let's loop blocks, and concatenate Env within each block
    for b = blocks
        savedir = fullfile(fld.procdatadir, ...
            monkey, [monkey '_' cdate '_' num2str(b)]);
        mkdir(savedir); % create a directory to save processed data
        mkdir(fullfile(savedir, 'individual_channels'));
        
        disp(['block ' num2str(b) '...']);
        % make sure we have 2 files, Env1 and Env2
        fn = dir(fullfile(datadir, [monkey cdate '-block-' num2str(b) '*']));
        available_envs = [];
        for ff = 1:length(fn)
            available_envs(end+1) = str2num(fn(ff).name(end-4));
        end
        envs = unique(available_envs);
        
        % group the Env files so we have all channels together
        env = cell(0);
        for f = envs
            % each file contains two structures:
            % - Data (contains ENV)
            % - Info (unused)
            load(fullfile(datadir,[monkey cdate '-block-' num2str(b) '-ENV' num2str(f)]),'Data');

            % convert ENV to better format
            env_t = Data.Time;
            for ch = 1:size(Data.Signal,2)
                env{ch + (f-1)*96} = squeeze(Data.Signal(:,ch,:));
            end
        end
        
        % TrackerFile is the name of the directory, as well as the filename
        % the mat file will give us LOG, Par, and STIM.
        TrackerFile = [monkey '_' cdate '_B00' num2str(b)];
        load(fullfile(logdir,TrackerFile,[TrackerFile '.mat']),'LOG','Par','STIM');
        
        %% convert LOG to LUT
        lut = struct2table(LOG.Trial);
        
        monkey = repmat({monkey},size(lut,1),1);
        session = repmat({cdate},size(lut,1),1);
        block = repmat(b,size(lut,1),1);
        lut = addvars(lut,monkey,session,block,'Before','TrialNr');
        monkey = monkey{1};
        
        % we need to compute some extra variables
        statusMapping = {'Aborted','Target','Nontarget','Distractor','NoHit'};
        trialPerf = cellfun(@(x) find(strcmp(statusMapping,x))-1,lut.Status);
        prevTrialPerf = [nan; trialPerf(1:end-1)];
        prevStatus = [{''}; lut.Status(1:end-1)];
        
        % colswitch says whether the ColScheme of current and previous
        % trials were the same or not
        colswitch = zeros(size(lut,1),1);
        colswitch(1) = nan;
        colswitch(abs(lut.ColScheme - [nan; lut.ColScheme(1:end-1)]) > 0) = 1;
        
        % make variables containing previous reward values and current
        % reward values (1 = low, 2 = high), based on the reward value
        % settings for this session
        rewardvals = sort(STIM.RewardValues);
        rewardval = ones(size(lut,1),1);
        rewardval(lut.RewValue==rewardvals(2)) = 2;
        prevrewardval = [nan; rewardval(1:end-1)];
        
        % check whether the current trial immediately followed the previous
        % one (so diff between trial numbers should be 1, not more)
        immediateFollow = zeros(size(lut,1),1);
        immediateFollow(find(diff(lut.TrialNr)==1)+1) = 1; % we have to add 1 because of the offset of the diff function
        
        % now add the previously computed vectors to the LUT
        lut.trialPerf = trialPerf; % range: 0 - 4, denotes the chosen target (see statusmapping)
        lut.prevTrialPerf = prevTrialPerf; % same as trialPerf, but for previous trial
        lut.colSwitch = colswitch; % boolean: 0 = no switch
        lut.rewardVal = rewardval; % nan = no reward, 1 = low reward, 2 = high reward
        lut.prevRewardVal = prevrewardval;
        lut.immediateFollow = immediateFollow; % boolean, 0 = there were excluded trials in between
        lut = addvars(lut,prevStatus,'After','Status');
        
        %% select trials based on human analogs
        % 1.	RT_colswitch_HR and Perf_colswitch_HR
        % 2.	RT_colswitch_LR and Perf_colswitch_LR
        % 3.	RT_colsame_HR and Perf_colsame_HR
        % 4.	RT_colsame_LR and Perf_colsame_LR
        lut.RT_ColSwitch_HR = strcmp(lut.Status,'Target') & strcmp(lut.prevStatus,'Target') & lut.colSwitch==1 & lut.prevRewardVal==2 & lut.immediateFollow==1;
        lut.RT_ColSwitch_LR = strcmp(lut.Status,'Target') & strcmp(lut.prevStatus,'Target') & lut.colSwitch==1 & lut.prevRewardVal==1 & lut.immediateFollow==1;
        lut.RT_ColSame_HR = strcmp(lut.Status,'Target') & strcmp(lut.prevStatus,'Target') & lut.colSwitch==0 & lut.prevRewardVal==2 & lut.immediateFollow==1;
        lut.RT_ColSame_LR = strcmp(lut.Status,'Target') & strcmp(lut.prevStatus,'Target') & lut.colSwitch==0 & lut.prevRewardVal==1 & lut.immediateFollow==1;
        
        lut.Perf_ColSwitch_HR = strcmp(lut.prevStatus,'Target') & lut.colSwitch==1 & lut.prevRewardVal==2 & lut.immediateFollow==1;
        lut.Perf_ColSwitch_LR = strcmp(lut.prevStatus,'Target') & lut.colSwitch==1 & lut.prevRewardVal==1 & lut.immediateFollow==1;
        lut.Perf_ColSame_HR = strcmp(lut.prevStatus,'Target') & lut.colSwitch==0 & lut.prevRewardVal==2 & lut.immediateFollow==1;
        lut.Perf_ColSame_LR = strcmp(lut.prevStatus,'Target') & lut.colSwitch==0 & lut.prevRewardVal==1 & lut.immediateFollow==1;
        
        %% correct the reaction time (make reactiontime2)
        % there are three vectors in the LUT that have timing information
        % that we need to calculate reactiontime2. however, these are cell
        % arrays, that have some empty cells, so we can't use cell2mat
        % immediately, because that will skip the empty ones. let's fix,
        % and calculate reactTime2.
        rt = nan(length(lut.ReactTime),1);
        rt(cellfun('isempty',lut.ReactTime)==0,1) = cell2mat(lut.ReactTime);
        lut.ReactTime = rt;
        sacc = nan(size(lut,1),1);
        sacc(cellfun('isempty',lut.SaccTime)==0,1) = cell2mat(lut.SaccTime);
        lut.SaccTime = sacc;
        
        % apparently, StimTarDur is usually a cell array, but not always
        if iscell(lut.StimTarDur)
            stimtar = nan(size(lut,1),1);
            stimtar(cellfun('isempty',lut.StimTarDur)==0,1) = cell2mat(lut.StimTarDur);
            lut.StimTarDur = stimtar;
        end 
        
        if STIM.KeepFixT < 175
            lut.reactTime2 = lut.ReactTime + lut.SaccTime + lut.StimTarDur*1000;
        end
        
        %% check which stimulus was in RF when
        % we need to define 6 variables
        % Tar - list of trials in which Target is in RF
        % iTar - same as Tar, but with proper trial number, corresponding to LUT
        % Dist - etc...
        % iDist
        % NonTar
        % iNonTar
        
        % we first load some RF information. this is a structure with 3
        % maps, one for M1, and 2 for M2. Apparently, there are two
        % configurations for M2, one up to 20141101, and another from
        % that date onwards. we need to make sure we use the right ones.
        % since only M2 has two configurations, and the date at which
        % they switch is hardcoded in Hickey_RFMAPS, I figured we can just
        % hardcode it here instead, to keep it simple.
        if strcmp(monkey,'M1')
            mp = RFMAPS_function(1);
        elseif strcmp(monkey,'M2')
            dt = datetime(cdate,'InputFormat','yyyyMMdd');
            if dt > datetime('20141101','InputFormat','yyyyMMdd') % the original script uses >, not >=. However, as far as I know, there is no recording from that actual date, so it doesn't really matter
                mp = RFMAPS_function(3);
            else
                mp = RFMAPS_function(2);
            end
        else
            error('Wrong monkey?');
        end
        
        % now we have a table with 1 row per channel. let's loop them, and
        % find, per channel, which trials we should include
        qn = size(mp,1);
        mp.iTar = cell(qn,1);
        mp.iDist = cell(qn,1);
        mp.iNonTar = cell(qn,1);
        mp.Tar = cell(qn,1);
        mp.Dist = cell(qn,1);
        mp.NonTar = cell(qn,1);
        for ch = 1:size(mp,1)
            if ~isempty(mp.rfpos(ch))
                mp.iTar{ch} = find(lut.TarPos==mp.rfpos(ch));
                mp.iDist{ch} = find(lut.DistPos==mp.rfpos(ch));
                mp.iNonTar{ch} = find(lut.TarPos~=mp.rfpos(ch) & lut.DistPos~=mp.rfpos(ch));
                mp.Tar{ch} = lut.TrialNr(mp.iTar{ch});
                mp.Dist{ch} = lut.TrialNr(mp.iDist{ch});
                mp.NonTar{ch} = lut.TrialNr(mp.iNonTar{ch});
            end
        end
        
        %% make some figures
%         plot_rt(lut);
%         plot_perf(lut);
%         plot_perf_fraction(lut);
        
        %% process MUA data
        
        % find nonaborted trials
        nonaborted = ~strcmp(lut.Status,'Aborted') & ~strcmp(lut.Status,'NoHit');
        
        RTMask=[];
        latencydelay = 50;
        for t=1:size(lut,1)
            RTMask = [RTMask (env_t >= (lut.reactTime2(t)+latencydelay)/1000)'];
        end
        
        % loop channels 
        snr = nan(size(mp,1),1);
        snr_cons = nan(size(mp,1),1);
        allgoodtrials = nan(size(lut,1),size(mp,1));
        allsamplewisegoodtrials = nan(size(lut,1),size(mp,1));
        baseline_value = nan(size(mp,1),1);
        peakresp = nan(size(mp,1),1);
        for ch = 1:size(mp,1)
            if ch/10 == ceil(ch/10) || ch == 1
                disp(['analysing channel ' num2str(ch) '...']);
            end
            
            % check if we have data for this channel
            if length(env) < ch
                break
            end
            if isempty(env{ch})
                continue
            end
            
            % make sure that the lut trials line up with the env trials
            % (or at least that we have the same amount)
            if (size(lut,1) ~= size(env{ch},2))
                error('env and lut don''t align');
            end
            
            % preprocess 
            data = env{ch};
            raw_data = data;
            [goodtrials, sample_wise_goodtrials] = remove_outliers(data,3);
            [data, blval] = correct_baseline(data,sample_wise_goodtrials'&nonaborted,env_t,[-0.1 0]);
            baseline_value(ch) = blval;
            
            % HickeyAnalysis_CK uses masked data to calculate the peak
            % response for normalization. i tried both, and it doesn't seem
            % to matter (i.e. the peak occurs before the RT)
            data_masked = data;
            data_masked(RTMask==1) = nan;

            % use the masked_data to get the peak response and normalize
            % the masked data 
            [data_masked, cpeakresp] = normalize(data_masked,sample_wise_goodtrials'&nonaborted,env_t,[0 0.25]);
            peakresp(ch) = cpeakresp;
            
            % now also normalize the unmasked dataset
            data = data./cpeakresp;
            
            % get SNR based on the non-masked data 
            [snr(ch), snr_cons(ch)] = get_snr(data,sample_wise_goodtrials'&nonaborted,env_t,[-0.1 0]);
            
            allgoodtrials(:,ch) = goodtrials;
            allsamplewisegoodtrials(:,ch) = sample_wise_goodtrials;
%             Env{ch} = data;
%             Env_masked{ch} = data_masked;
%             Env_raw{ch} = raw_data;
            
            % save individual channels
            save(fullfile(savedir,'individual_channels',...
                ['Env_preprocessed_ch' num2str(ch) '_' monkey '_' cdate '_' num2str(b)]),...
                'data','data_masked','raw_data','cpeakresp','blval','goodtrials','sample_wise_goodtrials');
        end
        
        mp.snr = snr;
        mp.snr_conservative = snr_cons;
        mp.goodtrials = allgoodtrials';
        mp.sample_wise_goodtrials = allsamplewisegoodtrials';
        mp.baseline = baseline_value;
        mp.peak = peakresp;
        
        % save lut and mp
        save(fullfile(savedir,[monkey '_' cdate '_' num2str(b) '_LUT']),'lut');
        save(fullfile(savedir,[monkey '_' cdate '_' num2str(b) '_channels']),'mp');
        
        % clear some big variables and then save everything
        clear Env env data;
        save(fullfile(savedir, [monkey '_' cdate '_' num2str(b) '.mat']));
    end
end
end


function fig = plot_rt(lut)
    fig = figure; set(gcf,'Position',[ 616   234   480   577]);
        
    % colswitch
    lr = lut.reactTime2(lut.RT_ColSwitch_LR);
    hr = lut.reactTime2(lut.RT_ColSwitch_HR);
    mn_lr = mean(lr);
    mn_hr = mean(hr);
    std_lr = std(lr)./sqrt(length(lr));
    std_hr = std(hr)./sqrt(length(hr));
    errorbar([.99 1.99],[mn_lr mn_hr],[std_lr,std_hr], 'bo-'); hold all;

    % colsame
    lr = lut.reactTime2(lut.RT_ColSame_LR);
    hr = lut.reactTime2(lut.RT_ColSame_HR);
    mn_lr = mean(lr);
    mn_hr = mean(hr);
    std_lr = std(lr)./sqrt(length(lr));
    std_hr = std(hr)./sqrt(length(hr));
    errorbar([1.01 2.01],[mn_lr mn_hr],[std_lr,std_hr], 'ro-');

    hl=legend({'Color swap','Color same'},'Location','SouthWest');
    set(gca,'xtick',1:2,'xticklabel',{'Low','High'},'FontSize',14);
    xlabel('Previous Trial Reward','FontSize',16);
    ylabel('Reaction time (ms)','FontSize',16);
    title(sprintf(['RT ' lut.monkey{1} ' ' lut.session{1} ' B' num2str(lut.block(1))]),'FontSize',18,'FontWeight','Demi');
    xlim([0.9 2.1]);
end

function fig = plot_perf(lut)
    figure; set(gcf,'Position',[574   212   560   420]);
        
    % colswitch
    lr = lut.Status(lut.Perf_ColSwitch_LR);
    lrp = sum(strcmp(lr,'Target'))./length(lr);
    hr = lut.Status(lut.Perf_ColSwitch_HR);
    hrp = sum(strcmp(hr,'Target'))./length(hr);
    bar([1:2],[lrp hrp],'FaceColor','b'); hold all;

    % colsame
    lr = lut.Status(lut.Perf_ColSame_LR);
    lrp = sum(strcmp(lr,'Target'))./length(lr);
    hr = lut.Status(lut.Perf_ColSame_HR);
    hrp = sum(strcmp(hr,'Target'))./length(hr);
    bar([3:4],[lrp hrp],'FaceColor','r');

    hl=legend({'Color swap','Color same'},'Location','NorthEastOutside', 'FontSize',14);
    set(gca,'xtick',1:4,'xticklabel',{sprintf('Low'),sprintf('High'),sprintf('Low'),sprintf('High')},'FontSize',14);
    xlabel('Previous Reward','FontSize',16);
    ylabel('Fraction correct','FontSize',16);
    title(['RT ' lut.monkey{1} ' ' lut.session{1} ' B' num2str(lut.block(1))],'FontSize',18,'FontWeight','Demi');
end

function fig = plot_perf_fraction(lut)
    fig = figure; set(gcf,'Position',[435   267   996   477]);
    usedata = {lut.Perf_ColSwitch_LR, lut.Perf_ColSwitch_HR, lut.Perf_ColSame_LR, lut.Perf_ColSame_HR};
    for i = 1:4
        u = usedata{i};
        subplot(2,2,i)

        b = []; status_list = {'Target','Nontarget','Distractor','NoHit'};
        for bi = 1:4
            b(bi) = sum(strcmp(lut.Status(u),status_list{bi}))/sum(u);
        end

        bar(1:4, b, 'FaceColor','k');
        set(gca,'xtick',1:4,'xticklabel',{'Targ','nonTarg','Distr','NoHit'},'FontSize',11,'ylim',[0 1]);
        ylabel('Fraction of trials','FontSize',12);
        title('Color Swap / High Reward','FontSize',14)
        
    end
    st=mtit(['Choices ' lut.monkey{1} ' ' lut.session{1} ' B' num2str(lut.block(1))],'xoff',0,'yoff',0.04);
end