macaque_MST_map1/code/technical_validation.m

% Note: This script is meant to be run as a whole;
%       assigning the directories only works if the script is run, not if
%       individual lines are executed.
%       If you want to run it line by line, set the directories manually 
%       (e.g. by setting codeDir = pwd if you are currently in the directory 
%       that contains this file).
%       On a late 2013 iMac with a 3.4 GHz Quad-Core Intel Core i5 this
%       takes around 1.5 minutes to run.
clc; 

% Specify directories
codeDir = fileparts(mfilename('fullpath'));
mainDir = codeDir(1:find(codeDir==filesep,1,'last'));
dataDir = [mainDir,'data/'];
% Add all folders with code and data to the path
addpath(genpath(mainDir));


% Specify and initiate some variables
experiments = {'MSTm', 'MSTt', 'MSTn'};
param2read  = {'/event_value', '/event_time'};

nCells = 172;
nExperiments = length(experiments);

%% Questions to be answered
% Are there negative spike times?
negative_spiketime_counter = 0;

% Are the number of timestamps for every event identical to the number of event values?
unequeal_value_times_counter = 0;

% Are trial end times are always later than trial start times, i.e., do all trials have a positive duration?
trial_durations = cell(nCells,nExperiments);
% What is the range and average of trial durations for hit trials (from trial start to trial outcome assignment?)
hit_trial_durations = cell(nCells,nExperiments);

% Does the average firing rate per trial differ between Mapping and Tuning?
% What is the average per trial firing rate? (for all trials and for trials >500ms)
avg_firing_rates = nan(nCells,nExperiments);
avg_firing_rates_500 = nan(nCells,nExperiments);

% What is the highest per trial firing rate? (for all trials >500ms)
highest_firing_rate = 0;
highest_firing_rate_trial = nan;
highest_firing_rate_datafile = '';


% How many trials lasting longer than 500ms have 0 spikes?
zero_spike_trials = 0;
trials_longer_500 = 0;

% What is the longest interspike interval?
ISIs = cell(nCells,nExperiments);



%% Read in and loop through data
% Read in meta data
meta_data = readtable([dataDir,'meta_data.txt']);


% Loop through all recordings
%  and all units within each recording for plausibility checks
cell_idx = 0;
for rec = 1:size(meta_data,1) % loop through all recordings
    for exp = 1:length(experiments) % loop through the three experiments
        % uncomment this line to get output of which file is currently
        % being loaded
        %fprintf("Loading experiment file %s\n", ['amm-',experiments{exp},'-',meta_data.monkey{rec},'-',strip(meta_data.session_number{rec},'both','"'),'-',meta_data.daily_count{rec},'-task.h5']);
        % get number of trials; will be 0 for nonexistant files
        numTrials = eval(['meta_data.Exp_',experiments{exp},'_tt(',num2str(rec),')']);
        % check if there is a file for this experiment in this recording
        if numTrials > 0
            % change into the directory of the respective experiment
            cd([dataDir,experiments{exp}]);
            % specify file name
            datafile_h5 = ['amm-',experiments{exp},'-',meta_data.monkey{rec},...
                '-',strip(meta_data.session_number{rec},'both','"'),'-',...
                meta_data.daily_count{rec},'-task.h5'];
            % get information about HDF5 file
            file_info = h5info(datafile_h5, '/event_value');
            % get names of events
            events2read = {file_info.Datasets.Name};
            % extract values and times of each event
            [event_value, event_time] = h5_extract(datafile_h5, param2read, events2read);
            
            % check that the number of trial starts, end, types, and
            % outcomes is the same
            if length(event_value.TRIAL_start) ~= length(event_value.TRIAL_end) || ...
                    length(event_value.TRIAL_start) ~= length(event_value.TRIAL_type) || ...
                    length(event_value.TRIAL_start) ~= length(event_value.TRIAL_outcome)
                warning("Different number of Trial starts, ends, types, and/or outomes in recordings %s",datafile_h5);
            end
            
            event_names = fieldnames(event_value);
            
            %% Validation of stimulus/behavioral data
            % check that number of event values == number of event times
            for fn = 1:length(event_names)
                if length(eval(['event_time.',event_names{fn}]))~=length(eval(['event_value.',event_names{fn}]))
                    warning('Number of event times and event values did not agree for %s in file %s\n',event_names{fn},datafile_h5);
                    unequeal_value_times_counter = unequeal_value_times_counter+1;
                end
            end
            
            % get trial durations for this recording and this experiment
            trial_durations{rec,exp} = (event_time.TRIAL_end - event_time.TRIAL_start)/1000; % in milliseconds
            
            
            % get time from trial start until reward for hit trials
            for tr = 1:numTrials
                if strcmp(event_value.TRIAL_outcome{tr},'hit')
                    tmp = (event_time.TRIAL_outcome(tr)-event_time.TRIAL_start(tr))/1000; % in ms
                    if tmp < 3000 || tmp > 4750 
                        warning("Short or long hit trial duration on trial %d of %s\n",tr,datafile_h5);
                    end
                    hit_trial_durations{rec,exp} = cat(2,hit_trial_durations{rec,exp},tmp);
                end
            end
            
            
            %% Validation of neural data            
            for u = 1:3 % loop through up to 3 offline sorted units per recording
                unit = eval(['meta_data.offline_units_',num2str(u),'(',num2str(rec),')']);
                % this will be NaN if no units have been offline sorted
                spike_trials = str2num(strip(eval(['meta_data.Exp_',experiments{exp},'_st_',num2str(u),'{',num2str(rec),'}']),"'")); %
                if ~isnan(unit) && length(spike_trials) > 1% if there is a unit
                    if exp==1
                        cell_idx = cell_idx+1; % increase cell counter only if it's the first experiment
                    end
                    % get the times of all spikes that were sorted to the current unit
                    spiketimes = event_time.SPIKE_channelUnit(event_value.SPIKE_channelUnit==unit);
                    
                    % Check for negative spike times
                    if sum(spiketimes<0)>0
                        fprintf("Negative spike times in datafile %s", datafile_h5);
                        negative_spiketime_counter = negative_spiketime_counter + sum(spiketimes<0);
                    end
                    
                    
                    % get per trial firing rates
                    tmp_firing_rates = [];
                    tmp_firing_rates_500 = [];
                    tmp_ISIs = [];
                    for tr = spike_trials
                        trialStart = event_time.TRIAL_start(tr);
                        trialEnd = event_time.TRIAL_end(tr);
                        trialDuration = (trialEnd-trialStart)/1000000; % in seconds
                        numSpikes = sum(spiketimes>trialStart & spiketimes < trialEnd);
                        tmp_ISIs = cat(1,tmp_ISIs,diff(spiketimes(spiketimes>trialStart & spiketimes < trialEnd))/1000);
                        if trialDuration > 0.5
                            trials_longer_500 = trials_longer_500+1;
                            if numSpikes == 0
                                zero_spike_trials = zero_spike_trials+1;
                            end
                            
                            if numSpikes/trialDuration > highest_firing_rate
                                highest_firing_rate = numSpikes/trialDuration;
                                highest_firing_rate_trial = tr;
                                highest_firing_rate_datafile = [datafile_h5,'-',num2str(unit)];
                            end
                        end
                        tmp_firing_rates = cat(2,tmp_firing_rates,numSpikes/trialDuration);
                        if trialDuration > 0.5
                            tmp_firing_rates_500 = cat(2,tmp_firing_rates_500,numSpikes/trialDuration);
                        end
                    end
                    
                    avg_firing_rates(rec,exp) = mean(tmp_firing_rates);
                    avg_firing_rates_500(rec,exp) = mean(tmp_firing_rates_500);
                    
                    
                    % interspike intervals
                    ISIs{cell_idx, exp} = tmp_ISIs; % in ms

                end % end of ~isnan(unit)
            end % end of u-loop
        end % end of if statement that checks whether file exists
    end % end of exp-loop
end % end of rec-loop


% all_firing_rates = [trial_firing_rates{:}];
all_trial_durations = vertcat(trial_durations{:});
all_hit_trial_durations = horzcat(hit_trial_durations{:});



%% Report results
fprintf('\n\nOverall Results\n');
fprintf('===============\n\n');


fprintf('General Points\n');
fprintf('---------------\n');


if negative_spiketime_counter == 0
    fprintf('There are no negative spike times.\n');
else
    fprintf('There are %d spikes with negative spiketimes.\n', negative_spiketime_counter);
end

if sum(all_trial_durations<0) == 0
    fprintf('Trial end times are always later than trial start times, i.e., all trials have a positive duration.\n')
else
    fprintf('%d trials have a negative duration.\n', sum(all_trial_durations<0));
end

if unequeal_value_times_counter == 0
    fprintf('The number of timestamps for every event is identical to the number of event values.\n');
else
    fprintf('For %d cases the number of event values did not agree with the number of timestamps.\n',unequeal_value_times_counter);
end

fprintf('\nTrial duration\n');
fprintf('---------------\n');

fprintf("For hit trials, the time from trial start to reward delivery ranged from %.2f to %.2f ms.\n", min(all_hit_trial_durations),max(all_hit_trial_durations));
fprintf("The average time from trial start to reward delivery was %.2f ms (SD: %.2f ms)\n",mean(all_hit_trial_durations), std(all_hit_trial_durations))




fprintf('\nFiring Rate\n');
fprintf('------------\n')
fprintf('The average firing rate in the Mapping experiment was %.2f, compared to %.2f for Tuning and %.2f for RC\n',mean(avg_firing_rates,'omitnan'));
[H,P,CI,STATS] = ttest(avg_firing_rates(:,2), avg_firing_rates(:,1));
if P < .05
    fprintf("The difference between Mapping and Tuning was significant, ");
else
    fprintf("The difference between Mapping and Tuning was not significant, ");
end
fprintf("t(%d) = %.2f, p = %.4f\n", STATS.df, STATS.tstat, P);


fprintf('Across all recordings and all experiments, %d trials (%.2f%%) that last at least 500ms have 0 spikes.\n', zero_spike_trials, 100*zero_spike_trials/trials_longer_500);



fprintf('Across all units, all experiments, and all trials that last at least 500ms,\n');
fprintf('the average firing rate (calculated per trial) is  %.2f spikes/s (SD: %.2f),\n', mean(avg_firing_rates_500(:),'omitnan'), std(avg_firing_rates_500(:),'omitnan'));
fprintf('and the highest firing rate calculated for one individual trial is  %.2f spikes/s on trial %d of cell %s.\n',highest_firing_rate,highest_firing_rate_trial,highest_firing_rate_datafile);

all_ISIs = vertcat(ISIs{:});
fprintf('\nInterspike Intervals\n');
fprintf('---------------------\n')
fprintf('The longest interspike interval is %.0f milliseconds\n',max(all_ISIs));
fprintf('Across all recordings and experiments, there were %d interspike intervals (%.4f%%) that were longer than 4s.\n',sum(all_ISIs>=4000), 100*sum(all_ISIs>=4000)/length(all_ISIs));
fprintf('Across all recordings and experiments, there were %d interspike intervals (%.2f%%) that were longer than 2s.\n',sum(all_ISIs>2000), 100*sum(all_ISIs>2000)/length(all_ISIs));
fprintf('%d interspike intervals (%.2f%%) were shorter than 1ms.\n',sum(all_ISIs<=1), 100*sum(all_ISIs<=1)/length(all_ISIs));