plot_tuned_units.m

% plot tuned units

close all 
clear all
clc

%%
subject_id = 's3';
unit_region = 'SMG';

spike_sorting_type = '_unsorted_aligned_thr_-4.5';
flag_4S = true; % true = updated 4S action phase; false = original 2S action phase
flag_shuffled = false; % true = shuffled images task
flag_varied_sizes = false; % true for varied sizes task
flag_GB_images = false; % true for task using images of GB's own hands and real objects
flag_5050 = false; % true for 50% Go 50% NoGo trials
flag_combined = false; % true for combinations task

if ~flag_4S
    TaskCue = 'GraspObject';
    min_timebin_length = 134; % NOT VALID FOR 20230831    
elseif ~flag_shuffled
    TaskCue = 'GraspObject_4S_Action';
    min_timebin_length = 174;
else
    TaskCue = 'GraspObject_Shuffled';
    min_timebin_length = 174; 
end 

if flag_varied_sizes
    TaskCue = 'GraspObject_Varied_Size';
    min_timebin_length = 174;
end
if flag_GB_images
    TaskCue = 'GraspObject_GB_Images';
    min_timebin_length = 174;
end

if flag_5050
    TaskCue = 'GraspObject_5050';
    min_timebin_length = 174;
end

if flag_combined
    TaskCue = 'GraspObject_Combined';
    min_timebin_length = 174;
end

% Task Variables
% 4S data
Data = load(['C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\' subject_id '\Data\Table_' subject_id '_' TaskCue spike_sorting_type]);
Go_data = Data.Go_data;
color_info = {[.1176 .5333 .8980],[.8471 .1059 .3765],[1 .7569 .0275]};

% remove faulty sessions, if any
error_session = {};
if strcmp(subject_id, 's2')
    error_session = {'20231016'};
elseif strcmp(subject_id, 's3')
    error_session = {'20231207','20231212','20250212'};
elseif strcmp(subject_id, 's4')
    error_session = {'20240613'};
end 

if ~isempty(error_session)
    for session = 1:numel(error_session)
        err_session = error_session{session};
        condition = cellfun(@(x) strcmp(x, err_session), Go_data.session_date);
        Go_data = Go_data(~condition,:);
    end
end

% add Aperture Size column
if strcmp(TaskCue, 'GraspObject_Varied_Size')
    sizeKeywords = ['Small', 'Medium', 'Large'];
    Go_data.Aperture_Size = cell(height(Go_data),1);
    % Loop through each label and extract the size information
    for i = 1:height(Go_data)
        % Use regular expression to find the size keyword after the last underscore
        tokens = regexp(Go_data.LabelNames{i}, '_(Small|Medium|Large)$', 'tokens');
    
        if ~isempty(tokens)
            % tokens is a cell array; extract the size keyword from it
            Go_data.Aperture_Size{i} = tokens{1}{1};
        end
    end
end

if strcmp(TaskCue, 'GraspObject_Combined')
    Go_data.TrialType(strcmp(Go_data.TrialType, 'Unknown')) = {'Combined'}; % adds in Combined as Trial type
    % add in column with Object Type for Combined trials and original trial types (H, HO, O with Associated)
    % Loop through each label and extract the object information
    for i = 1:height(Go_data)
        % Use regular expression to find the size keyword after the last underscore
        tokens = regexp(Go_data.LabelNames{i}, '_(deck|block|rod|ball)$', 'tokens');
        
        if ~isempty(tokens)
            % tokens is a cell array; extract the size keyword from it
            Go_data.ObjectType{i} = tokens{1}{1};
        else 
            Go_data.ObjectType{i} = 'Associated';
        end
    end
    color_info = {[.3632 .2266 .6055],[.1176 .5333 .8980],[.8471 .1059 .3765],[1 .7569 .0275]}; % Combinations task (purple at beginning)
end

flagGoTrials = true; % false = No-Go

flagRegressionTuning = true;

if flagRegressionTuning
    analysis_type = 'LinearRegression';
else
    analysis_type = 'KruskalWallis';
end

flagBinPerBin = true;
multipleComparePhase = true;
flagTunedChannels = true;
flagSaveData = true;

%chose cue type:
taskCuesAll = {'Hand', 'Hand-Object', 'Object'};
if flag_combined
    taskCuesAll = {'Combined','Hand', 'Hand-Object', 'Object'};
end
sessions_all = unique(Go_data.session_date);
numSessions = numel(sessions_all);
phase_time_idx = Go_data.time_phase_labels{1,1};
numPhases = numel(unique(phase_time_idx));
phase_changes_idx = diff(phase_time_idx);
phase_changes(1) = 1;
phase_changes(2:numPhases) = find(phase_changes_idx) + 1;
phase_bin_ranges = {
    1:phase_changes(2)-1;                   % ITI
    phase_changes(2):phase_changes(3)-1;    % Cue
    phase_changes(3):phase_changes(4)-1;    % Delay
    phase_changes(4):numel(phase_time_idx)  % Action
};
phaseNames = {'ITI', 'Cue', 'Delay', 'Action'};

numUnitsPerSession = zeros(numSessions,1);

% Initialize cell arrays to store results
hand_ho_overlap_units_all = cell(numSessions,1);
object_ho_overlap_units_all = cell(numSessions,1);
object_hand_overlap_units_all = cell(numSessions,1);
object_hand_ho_overlap_units_all = cell(numSessions,1);

%uniqueApertureSize = unique(Go_data.Aperture_Size); % comment out when not varied sizes
%% load Data
goLabel = ["NoGo", "Go"];
goLabel = goLabel(flagGoTrials + 1);
directory = ['C:\Users\macthurston\Documents\GitHub\project_grasp_object_interaction\analyzedData\' subject_id];
analysis_type = 'LinearRegression'; % 'LinearRegression' or 'KW'
filename = "UPDATED_grasp_tuned_channels_per_condition_" + TaskCue + '_' + unit_region + "_" + analysis_type + "_" + goLabel + ".mat"; % just .mat for original data
%filename = "tuned_channels_" + TaskCue + '_' + unit_region + "_" + analysis_type + "_" + 'Pixelated'; % Pixelated/Regular
full_path = fullfile(directory, filename);
load(full_path);

%% Planning vs. Action
tuned_channels_per_phase; % total tuning per phase for each condition
tuned_channels_per_phase_vector; % how I can tell if they are tuned during both Planning and Action

% totals
for n_type = 1:numel(taskCuesAll)
    total_units_per_phase = sum(cell2mat(tuned_channels_per_phase(n_type,:)'));
    tuned_planning(n_type) = total_units_per_phase(2);
    tuned_action(n_type) = total_units_per_phase(4);
end

% overlap
tuned_overlap = cell(numSessions,1);
for n_session = 1:numSessions
    for n_type = 1:numel(taskCuesAll)
        tuning_per_session_per_condition = tuned_channels_per_phase_vector{n_type,n_session};
        tuned_overlap_per_session{n_type} = sum((tuning_per_session_per_condition(:,2) == 1) & (tuning_per_session_per_condition(:,4) == 1),1);
    end
    tuned_overlap{n_session} = cell2mat(tuned_overlap_per_session);
end
total_overlap = sum(cell2mat(tuned_overlap),1);

% plot venn
data_for_plot = cell(1,numel(taskCuesAll));
figure('units','normalized','outerposition',[0 0 0.27 .2])
sgtitle(['Tuned Overlap in ' unit_region])
% restructure data
for n_type = 1:numel(taskCuesAll)
    data_for_plot{:,n_type} = [tuned_planning(n_type) tuned_action(n_type) total_overlap(n_type)];
    subplot(1,3,n_type)
    venn(data_for_plot{n_type},'FaceColor',{'r','b'},'FaceAlpha',{1,0.6},'EdgeColor','black');
    xticks([])
    yticks([])
    title([taskCuesAll{n_type}])
    planning_only = tuned_planning(n_type) - total_overlap(n_type);
    action_only = tuned_action(n_type) - total_overlap(n_type);
    total_units = planning_only + action_only + total_overlap(n_type);
    planning_only_perc = round(planning_only/total_units * 100);
    action_only_perc = round(action_only/total_units * 100);
    total_overlap_perc = round(total_overlap(n_type)/total_units * 100);
    disp([unit_region ': ' taskCuesAll{n_type} ' - Planning: ' num2str(planning_only) ' (' num2str(planning_only_perc) '%). Action: ' num2str(action_only) ' (' num2str(action_only_perc) '%). Overlap: ' num2str(total_overlap(n_type)) ' (' num2str(total_overlap_perc) '%).' ])
end

%% Grasp Types
graspTypes = unique(Go_data.GraspType);
tuned_channels_per_grasp; % condition x sessions (each double is 4x4 phase x grasp type)
% adapt code below to color bars per grasp type rather than condition
% could also just 
%%
phase_tuned_mean_all = zeros(numel(taskCuesAll), numPhases); % Mean percentages
phase_yCI95_all = zeros(numel(taskCuesAll), numPhases); % 95% CI
percentage_tuned_all = cell(numel(taskCuesAll), 1);
%sessionToInclude = setdiff(1:numSessions,1);

% code for empty/missing session data
rowsToKeep = numUnitsPerSession ~= 0;
numUnitsPerSession = numUnitsPerSession(rowsToKeep);
sessionToInclude = 1:numel(numUnitsPerSession);
colsToKeep = true(1,numSessions);
for n_session = 1:numSessions
    if all(cellfun('isempty',tuned_channels_per_phase(:,n_session)))
        colsToKeep(n_session) = false;
    end
end
tuned_channels_per_phase = tuned_channels_per_phase(:,colsToKeep);

figure('units','normalized','outerposition',[0 0 .38 0.38]); %[0 0 .5 .5]); %[0 0 .38 0.38]);
for n_type = 1:numel(taskCuesAll)
    dataTmp = cell2mat(tuned_channels_per_phase(n_type,sessionToInclude)')*100;
    percentage_tuned = dataTmp./numUnitsPerSession(sessionToInclude);
    percentage_tuned_all{n_type} = percentage_tuned;
    yCI95tmp = utile.calculate_CI(percentage_tuned); % Calculate 95% Probability Intervals Of t-Distribution using SEM
    phase_yCI95_all(n_type,:) = yCI95tmp(2,:);
    phase_tuned_mean_all(n_type,:) = mean(percentage_tuned,1);

    subplot(1,numel(taskCuesAll),n_type)

    hold on
    bar(phase_tuned_mean_all(n_type,:),'FaceColor',color_info{n_type});

    hold on
    errorbar(phase_tuned_mean_all(n_type,:),phase_yCI95_all(n_type,:),'Color','k');

    % % Scatter plot of individual session data points
    % for n_phase = 1:4
    %     scatter(repmat(n_phase, numel(sessionToInclude), 1), percentage_tuned(:,n_phase), 50, 'k', 'filled', 'MarkerFaceAlpha', 0.6);
    % end 

    title(taskCuesAll(n_type));
    xticks(1:numel(phaseNames));
    xticklabels(phaseNames);
    xtickangle(45);
    ylabel('% of Total Units');
    ylim([0 100]);
    yticks([0 50 100]);
    sgtitle(['Tuned Units in ' unit_region])
    set(gca, 'FontSize', 12);
end
%% single vs multi-tuned
% need tuned_channels_per_grasp_vector