tuning_analysis.m

clc
clear all
close all

subject_id = 's3';
unit_region = 'S1';

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_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 = true; % 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_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;

if flag_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
end

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

if ~isempty(error_session)
    condition = cellfun(@(x) strcmp(x, error_session), Go_data.session_date);
    Go_data = Go_data(~condition,:);
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;
phaseNames = {'ITI', 'Cue', 'Delay', 'Action'};
color_info = {[.1176 .5333 .8980],[.8471 .1059 .3765],[1 .7569 .0275]};
if flag_combined
    color_info = {[.3632 .2266 .6055],[.1176 .5333 .8980],[.8471 .1059 .3765],[1 .7569 .0275]};
end

numUnitsPerSession = zeros(numSessions,1);

% Initialize cell arrays to store results
hand_ho_overlap_units_all = cell(numSessions,1);
% hand_only_units_all = cell(numSessions,1);
% ho_only_units_h_all = cell(numSessions,1);
object_ho_overlap_units_all = cell(numSessions,1);
% object_only_units_all = cell(numSessions,1);
% ho_only_units_o_all = cell(numSessions,1);
object_hand_overlap_units_all = cell(numSessions,1);
% object_only_units_h_all = cell(numSessions,1);
% hand_only_units_o_all = cell(numSessions,1);
object_hand_ho_overlap_units_all = cell(numSessions,1);


%% Analysis
for n_session = 1:numSessions

    disp(['Classification session ' sessions_all{n_session} ]);  

    %find idx of current session day
    idxThisSession = ismember(Go_data.session_date, sessions_all(n_session));
        

    %extract data from selected brain area

     if strcmp('SMG', unit_region)
            SessionData = Go_data.SMG_Go(idxThisSession,:);
        elseif strcmp('PMV', unit_region) 
            SessionData = Go_data.PMV_Go(idxThisSession,:);
        elseif strcmp('S1', unit_region)
            SessionData = Go_data.S1X_Go(idxThisSession,:);
         elseif strcmp('M1', unit_region) 
            SessionData = Go_data.M1_Go(idxThisSession,:);
        elseif strcmp('AIP', unit_region)
            SessionData = Go_data.AIP_Go(idxThisSession,:);
        elseif strcmp('dlPFC', unit_region)
            SessionData = Go_data.dlPFC_Go(idxThisSession,:);
        else
            error([unit_region ' does not exist '])
     end

     % skip session days that are empty - relevant for S1 session 20230810
     if isempty(SessionData{1})
        continue
     end
    
    %labels 
    sessionLabels = Go_data.GoLabels(idxThisSession,:);
    
    %trialType
    trialTypeSession = Go_data.TrialType(idxThisSession,:);

    %get idx for Go or NoGo trials
    GoNoGoidx =  logical(cell2mat(Go_data.TrialCue(idxThisSession,:)));
    timePhaseLabels = Go_data.time_phase_labels(idxThisSession);
    

    if flagGoTrials
        SessionData = SessionData(GoNoGoidx);
        sessionLabels = sessionLabels(GoNoGoidx);
        timePhaseLabels = timePhaseLabels(GoNoGoidx);
        trialTypeSession = trialTypeSession(GoNoGoidx);
    else
        SessionData = SessionData(~GoNoGoidx);
        sessionLabels = sessionLabels(~GoNoGoidx);
        timePhaseLabels = timePhaseLabels(~GoNoGoidx);
        trialTypeSession = trialTypeSession(~GoNoGoidx);
    end
     
    %seperate data according to cue modality 

    unTrialType = unique(Go_data.TrialType);
    numUnitsPerSession(n_session) = size(SessionData{1},2);
    % loop through cue modalities 
    for n_type = 1:numel(unTrialType)

        % find idx of trial type 
        trialTypeIdx = ismember(trialTypeSession, unTrialType(n_type));

         if flagTunedChannels
            %Compute index of units that are tuned
            if flagRegressionTuning
            
                [tunedCombinedChannels, tunedChannelsPhase, tunedChannelsBin, sumPhase, sumBin,numTunedChannelsPerCategory,~,~,p_per_phase] ...
                     = classification.getRegressionTunedChannels_paper(SessionData(trialTypeIdx),sessionLabels(trialTypeIdx), ...
                 timePhaseLabels(trialTypeIdx), 'multcompare', multipleComparePhase, 'BinperBinTuning', flagBinPerBin);

                condToTest = arrayfun(@(x) preproc.image2class_simple(x),  unique(sessionLabels), 'UniformOutput', false);

                if nnz(sumBin) ~= 0
                    figure();
                    plot(sumBin);
                end

                tuned_channels_per_graps{n_type,n_session} = numTunedChannelsPerCategory;
            else

                tuned_channels_per_graps{n_type,n_session} = [];
                [tunedCombinedChannels, tunedChannelsPhase, tunedChannelsBin, sumPhase, sumBin]= classification.getTunedChannels(SessionData(trialTypeIdx),sessionLabels(trialTypeIdx), ...
                timePhaseLabels(trialTypeIdx), 'multcompare', multipleComparePhase,'removeITItuning', 'false', 'BinperBinTuning', flagBinPerBin);
                sumBin = sumBin';
            end

                if nnz(sumBin) > 0
                sum_bin_all{n_type, n_session } = sumBin;

            else
                sum_bin_all{n_type, n_session } = [];

                end

            tuned_channels_per_phase{n_type,n_session} = sumPhase;
            tuned_channels_per_phase_vector{n_type,n_session} = tunedChannelsPhase;


         
         end
    end 

    
    % calculating tuning overlap
    % H-HO overlap
    hand_ho_overlap_vector = (tuned_channels_per_phase_vector{1,n_session} == 1) & (tuned_channels_per_phase_vector{2,n_session} == 1); % this tells me the overlap between hand and hand-object units
    hand_ho_overlap_units = sum(hand_ho_overlap_vector, 1);
    hand_ho_overlap_units_all{n_session} = hand_ho_overlap_units;

    % hand_only_units = tuned_channels_per_phase{1,n_session} - hand_ho_overlap_units;
    % ho_only_units_h = tuned_channels_per_phase{2,n_session} - hand_ho_overlap_units;
    % hand_only_units_all{n_session} = hand_only_units;
    % ho_only_units_h_all{n_session} = ho_only_units_h;

    % O-HO overlap
    object_ho_overlap_vector = (tuned_channels_per_phase_vector{3,n_session} == 1) & (tuned_channels_per_phase_vector{2,n_session} == 1); % this tells me the overlap between object and hand-object units
    object_ho_overlap_units = sum(object_ho_overlap_vector, 1);
    object_ho_overlap_units_all{n_session} = object_ho_overlap_units;

    % object_only_units = tuned_channels_per_phase{3,n_session} - object_ho_overlap_units;
    % ho_only_units_o = tuned_channels_per_phase{2,n_session} - object_ho_overlap_units;
    % object_only_units_all{n_session} = object_only_units;
    % ho_only_units_o_all{n_session} = ho_only_units_o;

    % O-H overlap
    object_hand_overlap_vector = (tuned_channels_per_phase_vector{3,n_session} == 1) & (tuned_channels_per_phase_vector{1,n_session} == 1); % this tells me the overlap between object and hand units
    object_hand_overlap_units = sum(object_hand_overlap_vector, 1);
    object_hand_overlap_units_all{n_session} = object_hand_overlap_units;

    % object_only_units_h = tuned_channels_per_phase{3,n_session} - object_hand_overlap_units;
    % hand_only_units_o = tuned_channels_per_phase{1,n_session} - object_hand_overlap_units;
    % object_only_units_h_all{n_session} = object_only_units_h;
    % hand_only_units_o_all{n_session} = hand_only_units_o;

    % all 3 modalities overlap
    object_hand_ho_overlap_vector = (tuned_channels_per_phase_vector{3,n_session} == 1) & (tuned_channels_per_phase_vector{1,n_session} == 1) & (tuned_channels_per_phase_vector{2,n_session} == 1); % this tells me the overlap between object, HO, and hand units
    object_hand_ho_overlap_units = sum(object_hand_ho_overlap_vector, 1);
    object_hand_ho_overlap_units_all{n_session} = object_hand_ho_overlap_units;
end 

hand_ho_overlap_units_all_sessions = sum(cell2mat(hand_ho_overlap_units_all));
% hand_only_units_all = cell2mat(hand_only_units_all');
% ho_only_units_h_all = cell2mat(ho_only_units_h_all');
object_ho_overlap_units_all_sessions = sum(cell2mat(object_ho_overlap_units_all));
% object_only_units_all = cell2mat(object_only_units_all');
% ho_only_units_o_all = cell2mat(ho_only_units_o_all');
object_hand_overlap_units_all_sessions = sum(cell2mat(object_hand_overlap_units_all));
% object_only_units_h_all = cell2mat(object_only_units_h_all');
% hand_only_units_o_all = cell2mat(hand_only_units_o_all');
object_hand_ho_overlap_units_all_sessions = sum(cell2mat(object_hand_ho_overlap_units_all));

hand_total_units = sum(cell2mat(tuned_channels_per_phase(1,:)'));
ho_total_units = sum(cell2mat(tuned_channels_per_phase(2,:)'));
object_total_units = sum(cell2mat(tuned_channels_per_phase(3,:)'));

%% saving variables
goLabel = ["NoGo", "Go"];
goLabel = goLabel(flagGoTrials + 1);

% Create the filename using the brain region and analysis type
%filename = ['tuned_channels_' unit_region '_' TaskCue '_' analysis_type '_' goLabel '.mat']; % goLabel determines Go, NoGo label
filename = "tuned_channels_" + TaskCue + '_' + unit_region + "_" + analysis_type + "_" + goLabel + ".mat";

directory = ['C:\Users\macthurston\Documents\GitHub\project_grasp_object_interaction\analyzedData\' subject_id];
full_path = fullfile(directory, filename);

% Save the relevant variables with the dynamic filename
save(full_path, 'sum_bin_all', 'tuned_channels_per_phase', 'tuned_channels_per_phase_vector','numUnitsPerSession',...
    'hand_ho_overlap_units_all','object_ho_overlap_units_all','object_hand_overlap_units_all','object_hand_ho_overlap_units_all',...
    'hand_ho_overlap_units_all_sessions','object_ho_overlap_units_all_sessions','object_hand_overlap_units_all_sessions','object_hand_ho_overlap_units_all_sessions',...
    'hand_total_units','ho_total_units','object_total_units');

%% across sessions tuned units overlapping 

% tuned_channels_per_phase_vector; % 3 (modality) x 5 (sessions)
% % I can compare within sessions how much overlap there is and then average
% % the sessions together to get average overlap
% 
% % H-HO overlap
% hand_ho_overlap_vector = (tuned_channels_per_phase_vector{1,1} == 1) & (tuned_channels_per_phase_vector{2,1} == 1); % this tells me the overlap between hand and hand-object units
% % I can sum and then substract from the total to get the venn diagram
% hand_ho_overlap_units = sum(hand_ho_overlap_vector, 1);
% 
% tuned_channels_per_phase; % total units for each modality
% hand_only_units = tuned_channels_per_phase{1,1} - hand_ho_overlap_units;
% ho_only_units_h = tuned_channels_per_phase{2,1} - hand_ho_overlap_units;
% 
% % next find average by iterating through each session and then finding the
% % mean
% 
% % O-HO overlap
% object_ho_overlap_vector = (tuned_channels_per_phase_vector{3,1} == 1) & (tuned_channels_per_phase_vector{2,1} == 1); % this tells me the overlap between object and hand-object units
% % I can sum and then substract from the total to get the venn diagram
% object_ho_overlap_units = sum(object_ho_overlap_vector, 1);
% 
% tuned_channels_per_phase; % total units for each modality
% object_only_units = tuned_channels_per_phase{3,1} - object_ho_overlap_units;
% ho_only_units_o = tuned_channels_per_phase{2,1} - object_ho_overlap_units;
% 
% % next find average by iterating through each session and then finding the
% % mean
% 
% % O-H overlap
% object_hand_overlap_vector = (tuned_channels_per_phase_vector{3,1} == 1) & (tuned_channels_per_phase_vector{1,1} == 1); % this tells me the overlap between object and hand units
% % I can sum and then substract from the total to get the venn diagram
% object_hand_overlap_units = sum(object_hand_overlap_vector, 1);
% 
% tuned_channels_per_phase; % total units for each modality
% object_only_units_h = tuned_channels_per_phase{3,1} - object_hand_overlap_units;
% hand_only_units_o = tuned_channels_per_phase{1,1} - object_hand_overlap_units;
% 
% % next find average by iterating through each session and then finding the
% % mean
% 
% % all 3 modalities overlap
% object_hand_ho_overlap_vector = (tuned_channels_per_phase_vector{3,1} == 1) & (tuned_channels_per_phase_vector{1,1} == 1) & (tuned_channels_per_phase_vector{2,1} == 1); % this tells me the overlap between object and hand units
% % I can sum and then substract from the total to get the venn diagram
% object_hand_ho_overlap_units = sum(object_hand_ho_overlap_vector, 1);
% 
% tuned_channels_per_phase; % total units for each modality
% object_only_units_h = tuned_channels_per_phase{3,1} - object_hand_overlap_units;
% hand_only_units_o = tuned_channels_per_phase{1,1} - object_hand_overlap_units;

%% example (requires Statistics and Machine Learning Toolbox) => unsure if
% % can handle 3 inputs
% % Sample data (replace with your own data)
% set1 = randi([0, 1], 1, 100); % Binary data for set 1
% set2 = randi([0, 1], 1, 100); % Binary data for set 2
% 
% % Create a logical array for the Venn diagram
% venn_data = [sum(set1 & ~set2), sum(~set1 & set2), sum(set1 & set2)];
% 
% % Create a Venn diagram using vennplot
% figure;
% vennplot(venn_data, 'FaceColor', {'r', 'g', 'b'}, 'FaceAlpha', 0.5);
% 
% % Add labels
% vennlabel({'Set 1', 'Set 2'});
% 
% % Add a title
% title('Venn Diagram');
% 
% % Adjust the display
% axis equal;

%% 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 = "tuned_channels_" + TaskCue + '_' + unit_region + "_" + analysis_type + "_" + goLabel + ".mat";
full_path = fullfile(directory, filename);
load(full_path);

%% Overlap values
% Cue
all_overlap_cue = object_hand_ho_overlap_units_all_sessions(2);
h_o_overlap_cue = object_hand_overlap_units_all_sessions(2) - all_overlap_cue;
h_ho_overlap_cue = hand_ho_overlap_units_all_sessions(2) - all_overlap_cue;
h_only_units_cue = hand_total_units(2) - (all_overlap_cue + h_o_overlap_cue + h_ho_overlap_cue);

o_ho_overlap_cue = object_ho_overlap_units_all_sessions(2) - all_overlap_cue;
o_only_units_cue = object_total_units(2) - (all_overlap_cue + h_o_overlap_cue + o_ho_overlap_cue);

ho_only_units_cue = ho_total_units(2) - (all_overlap_cue + h_ho_overlap_cue + o_ho_overlap_cue);

cue_total_units = all_overlap_cue + h_only_units_cue + o_only_units_cue + ho_only_units_cue + h_ho_overlap_cue + h_o_overlap_cue + o_ho_overlap_cue;

% Action
all_overlap_action = object_hand_ho_overlap_units_all_sessions(4);
h_o_overlap_action = object_hand_overlap_units_all_sessions(4) - all_overlap_action;
h_ho_overlap_action = hand_ho_overlap_units_all_sessions(4) - all_overlap_action;
h_only_units_action = hand_total_units(4) - (all_overlap_action + h_o_overlap_action + h_ho_overlap_action);

o_ho_overlap_action = object_ho_overlap_units_all_sessions(4) - all_overlap_action;
o_only_units_action = object_total_units(4) - (all_overlap_action + h_o_overlap_action + o_ho_overlap_action);

ho_only_units_action = ho_total_units(4) - (all_overlap_action + h_ho_overlap_action + o_ho_overlap_action);

action_total_units = all_overlap_action + h_only_units_action + o_only_units_action + ho_only_units_action + h_ho_overlap_action + h_o_overlap_action + o_ho_overlap_action;

% Percentages
% Cue
all_overlap_cue_perc = (all_overlap_cue/cue_total_units)*100;
h_o_overlap_cue_perc = (h_o_overlap_cue/cue_total_units)*100;
h_ho_overlap_cue_perc = (h_ho_overlap_cue/cue_total_units)*100;
h_only_units_cue_perc = (h_only_units_cue/cue_total_units)*100;
o_ho_overlap_cue_perc = (o_ho_overlap_cue/cue_total_units)*100;
o_only_units_cue_perc = (o_only_units_cue/cue_total_units)*100;
ho_only_units_cue_perc = (ho_only_units_cue/cue_total_units)*100;

% Action
all_overlap_action_perc = (all_overlap_action/action_total_units)*100;
h_o_overlap_action_perc = (h_o_overlap_action/action_total_units)*100;
h_ho_overlap_action_perc = (h_ho_overlap_action/action_total_units)*100;
h_only_units_action_perc = (h_only_units_action/action_total_units)*100;
o_ho_overlap_action_perc = (o_ho_overlap_action/action_total_units)*100;
o_only_units_action_perc = (o_only_units_action/action_total_units)*100;
ho_only_units_action_perc = (ho_only_units_action/action_total_units)*100;

%% creating bar plots of overlap
% look at all conditions
conditionLabels = categorical({'H','H&HO','HO','H&HO&O','HO&O','O','O&H'});
conditionLabels = reordercats(conditionLabels,{'H','H&HO','HO','H&HO&O','HO&O','O','O&H'});
cueOverlap = [h_only_units_cue_perc h_ho_overlap_cue_perc ho_only_units_cue_perc all_overlap_cue_perc o_ho_overlap_cue_perc o_only_units_cue_perc h_o_overlap_cue_perc]; % Cue
actionOverlap = [h_only_units_action_perc h_ho_overlap_action_perc ho_only_units_action_perc all_overlap_action_perc o_ho_overlap_action_perc o_only_units_action_perc h_o_overlap_action_perc]; % Action

% Create the first bar chart
b1 = bar(conditionLabels, cueOverlap, 'FaceColor', 'b', 'FaceAlpha', 0.5); % Blue bars with transparency
hold on;

% Create the second bar chart
b2 = bar(conditionLabels, actionOverlap, 'FaceColor', 'r', 'FaceAlpha', 0.5); % Red bars with transparency

% Add labels and legend
xlabel('Condition');
ylabel('% of Tuned Units');
legend({'Cue', 'Action'});
title(['Overlapping Units Throughout Trial - ' unit_region]);

hold off;
%%
% look at overlap of all 3 conditions with SDs
% okay so hand_total_units gives the total number of tuned units during the
% H blocks across all sessions, same for the others (can also look at
% tuned_units_per_phase if needed). we also have the values of tuned units
% for the overlap of all 3 for each individual session stored in
% object_hand_ho_overlap_all. I want to look at the percentage of tuned
% units contained in the all overlapping condition, so I think I need to add the
% tuned units that were present in each condition and use that as the
% denomiator

allOverlap = cell2mat(object_hand_ho_overlap_units_all);
percentageOverlap = cell(numSessions,1);

for n_session = 1:numSessions
    overlap_cue = allOverlap(n_session,2);
    overlap_action = allOverlap(n_session,4);

    tmpData = cell2mat(tuned_channels_per_phase(:,n_session));
    total_units_cue = sum(tmpData(:,2));
    total_units_action = sum(tmpData(:,4));
    percentageOverlap{n_session} = [(overlap_cue/total_units_cue)*100 (overlap_action/total_units_action)*100];
end

% now find mean and SD


%% bar plot of tuned units w/ 95% CIs (work on getting them all on same plot)
%save('C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\Workspaces\LinearRegression\s2\s2_GraspObject_2S_unsorted_aligned_thr_-4.5_SMG_Example.mat','sum_bin_all')
%ExampleSMG = load('C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\Workspaces\LinearRegression\s2\s2_GraspObject_2S_unsorted_aligned_thr_-4.5_SMG_Example.mat');
% upload linear regression analysis
%load('C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\Workspaces\LinearRegression\s2\s2_GraspObject_2S_unsorted_aligned_thr_-4.5_PMV.mat');

% phase_yCI95 = [];
% phase_tuned_mean = [];
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]);
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);

    % figure()
    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');

    title(taskCuesAll(n_type));
    xticks(1:numel(phaseNames));
    xticklabels(phaseNames);
    xtickangle(45);
    ylabel('% of Total Units');
    ylim([0 70]);
    sgtitle(['Tuned Units in ' unit_region])
    set(gca, 'FontSize', 12);
end
%% ANOVA for percentage tuned
% formatting the data
numModalities = numel(taskCuesAll);

% Initialize the 3D array
percentage_tuned_3D = zeros(numSessions, numPhases, numModalities);

% Loop through each modality and fill the 3D array
for n_modality = 1:numModalities
    % Extract the data from the cell
    current_data = percentage_tuned_all{n_modality};
    
    % Check the size of current_data
    [numSessionsCheck, numPhasesCheck] = size(current_data);
    
    % Ensure the size matches the expected dimensions
    if numSessionsCheck ~= numSessions || numPhasesCheck ~= numPhases
        error('The dimensions of the data in cell %d do not match expected dimensions.', n_modality);
    end
    
    % Fill the 3D array
    percentage_tuned_3D(:, :, n_modality) = current_data;
end

[numSessions, numPhases, numModalities] = size(percentage_tuned_3D);

% Loop through each phase
for n_phase = 1:numPhases
    % Extract data for the current phase
    phaseData = squeeze(percentage_tuned_3D(:, n_phase, :)); % Shape: [numSessions, numModalities]

    % Reshape data for ANOVA
    data = [];
    modality = [];

    % Loop through each modality
    for n_modality = 1:numModalities
        % Append the classification percentage
        data = [data; phaseData(:, n_modality)];
        % Append modality identifier
        modality = [modality; repmat(n_modality, numSessions, 1)];
    end

    % Convert modality to categorical for ANOVA
    modality = categorical(modality);

    % Create a table with the reshaped data
    dataTable = table(data, modality);

    % Perform ANOVA for the current phase
    [p, tbl, stats] = anova1(dataTable.data, dataTable.modality, 'off');
    
    % Display ANOVA results
    disp(['Phase ' num2str(n_phase) ':']);
    disp(['p-value = ' num2str(p)]);
    
    % Perform post-hoc tests if ANOVA is significant
    if p < 0.05
        [c, m, h, gnames] = multcompare(stats, 'CType', 'tukey-kramer');
        disp('Post-hoc test results:');
        disp(c);
    end
end

%% bar plot w/o CIs

for n_type = 1:numel(unTrialType) 
    if numSessions ~= 1
        tunedUnitsPerType(n_type,:) = sum(cell2mat(tuned_channels_per_phase(n_type,:)'));

    else
        tunedUnitsPerType(n_type,:) = cell2mat(tuned_channels_per_phase(n_type,:)');

    end 
end

figure('units','normalized','outerposition',[0 0 0.2 0.35]);
b = bar((tunedUnitsPerType'./sum(numUnitsPerSession))*100);
%bar((((tunedUnitsPerType')*8)./sum(numUnitsPerSession))*100);
%bar(tunedUnitsPerType');
hold on;
for k = 1:numel(b)
    b(k).FaceColor = color_info{k}; % Assign colors from cell array
end
title(['Tuned Units in ' unit_region]);
xticks(1:numel(phaseNames));
xticklabels(phaseNames);
ylabel('% of Total Units');
%ylabel('# of Tuned Units');
ylim([0 70]);
yticks([0 35 70]);
%ylim([0 50]);
legend(taskCuesAll, 'Location', 'Best', 'Interpreter', 'none','FontSize',12);
set(gca, 'FontSize', 12);
hold off

% % for F30
% phaseNames = {'Action'};
% taskCuesAll = {'G', 'G+O'};
% %tuned_channels_per_phase = [23 36; 38 47]; % pulling out cue and action of H & H+O, specific session for proposal
% tuned_channels_per_phase = [50; 56]; % pulling out action of H & H+O, specific session for proposal
% %tuned_channels_per_phase = [20 39; 30 46]; % SMG - pulling out cue and
% %action of H & H+O, specific session for proposal, OVERLAPPING UNITS found
% %in hand_ho_overlap_units for each phase
% %tuned_channels_per_phase = [5 24; 16 26]; % AIP - cue and action for H & HO
% %tuned_channels_per_phase = [34 50; 38 56]; % M1 - cue and action for H & HO
% 
% tunedUnitsPerType = tuned_channels_per_phase;
% figure('Position',[500 500 200 300]);
% h = bar((tunedUnitsPerType'./sum(numUnitsPerSession))*100, 'FaceColor','flat');
% h.CData(1,:) = [0.1176, 0.5333, 0.8980];
% h.CData(2,:) = [0.8471, 0.1059, 0.3765];
% %bar((((tunedUnitsPerType')*8)./sum(numUnitsPerSession))*100);
% %bar(tunedUnitsPerType');
% hold on;
% %title([unit_region ' - ' phaseNames]);
% xticks(1:numel(taskCuesAll));
% xticklabels(taskCuesAll);
% xlim([0.5, 2.5]);
% %ylabel('% of Total Units');
% %ylabel('# of Tuned Units');
% ylim([0 100]);
% yticks([0 50 100]);
% %ylim([0 50]);
% %legend(taskCuesAll, 'Location', 'Best', 'Interpreter', 'none','FontSize',12);
% set(gca, 'FontSize', 12);
% hold off

%% line plot w/o CIs

for n_type = 1:numel(unTrialType)
    tunedUnitsPerTypeBin(n_type,:)  = sum(cell2mat(sum_bin_all(n_type,:)),2);

end 

figure('units','normalized','outerposition',[0 0 0.3 0.45]);
lp = plot((tunedUnitsPerTypeBin'./sum(numUnitsPerSession))*100,'LineWidth',2);
%plot((((tunedUnitsPerTypeBin')*8)./sum(numUnitsPerSession))*100,'LineWidth',2);
%plot(tunedUnitsPerTypeBin','LineWidth',2);
hold on
for n_phase = 1:numPhases
    xline(phase_changes(n_phase), 'k--', phaseNames{n_phase}, 'LineWidth', 1.5,'FontSize',12);
end
for k = 1:numel(lp)
    lp(k).Color = color_info{k}; % Assign colors manually
end
title(['Tuned Units Throughout Trial in ' unit_region]) % ' - ' sessions_all{n_session}]);
xlabel('Time Bins (50 ms)');
xlim([0 (min_timebin_length + 5)])
%xticks([0 50 100 150]);
ylabel('% of Total Units');
%ylabel('# of Tuned Units');
ylim([0 70]);
%yticks([0 20 40 60]);
%ylim([0 50]);
legend(taskCuesAll, 'Location', 'Best','FontSize',12);
set(gca, 'FontSize', 12);
hold off

%% for line plot w/ 95% CI
% Initialize variables
percentage_tuned_per_bin_all = cell(numel(taskCuesAll), 1);
per_bin_yCI95 = zeros(numel(taskCuesAll),min_timebin_length);
per_bin_tuned_mean = zeros(numel(taskCuesAll),min_timebin_length);
%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',sum_bin_all(:,n_session)))
        colsToKeep(n_session) = false;
    end
end
sum_bin_all = sum_bin_all(:,colsToKeep);

figure('units','normalized','outerposition',[0 0 0.5 0.4]);
err_bar = {};
for n_type = 1:numel(taskCuesAll)
    dataTmp = cell2mat(sum_bin_all(n_type,sessionToInclude))*100;
    percentage_tuned_per_bin = dataTmp./(numUnitsPerSession(sessionToInclude)');
    percentage_tuned_per_bin_all{n_type} = percentage_tuned_per_bin;
    yCI95tmp = utile.calculate_CI(percentage_tuned_per_bin');
    per_bin_yCI95(n_type,:) = yCI95tmp(2,:);
    per_bin_tuned_mean(n_type,:) = mean(percentage_tuned_per_bin,2);

    hold on
    err_bar{n_type} = plot(1:length(dataTmp),per_bin_tuned_mean(n_type,:),'Color', color_info{n_type},'LineWidth',2);

    ER = utile.shadedErrorBar(1:length(dataTmp),per_bin_tuned_mean(n_type,:),per_bin_yCI95(n_type,:));
    ER.mainLine.Color = color_info{n_type};
    ER.patch.FaceColor = color_info{n_type};
    ER.edge(1).LineStyle = 'none'; %color_info{n_type};
    ER.edge(2).LineStyle = 'none'; %color_info{n_type};
end

for n_phase = 1:numPhases
    xline(phase_changes(n_phase), 'k--', 'LineWidth', 1.5,'FontSize',12); %phaseNames{n_phase}
end

title(['Tuned Units Throughout Trial in ' unit_region]);
xlabel('Time Bins (50 ms)');
xlim([0 (min_timebin_length + 5)]);    %([30 134]); %shortened %(min_timebin_length + 5)]) % 5 chosen as a buffer
xticks(phase_changes);
xticklabels(phaseNames);
xtickangle(45);
ylabel('% of Total Units');
ylim([0 70]);
yticks([0 35 70]);
legend([err_bar{:}], taskCuesAll,'Location', 'Best','Interpreter', 'none','FontSize',12);
set(gca, 'FontSize', 12);

%% ANOVA for percentage tuned across first half of Cue
% I need to take the average percent of units across that 1st
% second for each session and then use those averages to run ANOVA and
% determine if there are any differences => Cue timebins = 42-82 (can find
% with var phase_changes)

% formatting the data
numModalities = numel(taskCuesAll);
numTimebins = min_timebin_length;
% Define the timebins of interest (1st half of Cue)
timebins_of_interest = 42:61;
num_timebins_of_interest = length(timebins_of_interest);

% Initialize storage for averaged data
halfCue_ave_percentage_tuned = zeros(numSessions, numModalities);

% Loop through each modality
for n_modality = 1:numModalities
    % Extract the data from the cell
    current_data = percentage_tuned_per_bin_all{n_modality}';
    
    % Average across the timebins of interest for each session
    halfCue_ave_percentage_tuned(:, n_modality) = mean(current_data(:, timebins_of_interest), 2);
end

% Reshape data for ANOVA
data = [];
modality = [];

% Loop through each modality
for n_modality = 1:numModalities
    % Append the averaged percentage of tuned units
    data = [data; halfCue_ave_percentage_tuned(:, n_modality)];
    % Append modality identifier
    modality = [modality; repmat(n_modality, numSessions, 1)];
end

% Convert modality to categorical for ANOVA
modality = categorical(modality);

% Create a table with the reshaped data
dataTable = table(data, modality);

% Perform ANOVA
[p, tbl, stats] = anova1(dataTable.data, dataTable.modality, 'off');

% Display ANOVA results
disp('ANOVA Results:');
disp(['p-value = ' num2str(p)]);

% Perform post-hoc tests if ANOVA is significant
if p < 0.05
    [c, m, h, gnames] = multcompare(stats, 'CType', 'tukey-kramer');
    disp('Post-hoc test results:');
    disp(c);
end

%%
% blub = percentage_tuned.*numUnitsPerSession';
% blub2= mean(blub');
% figure(); plot(zscore(blub2)); hold on; plot(zscore(per_bin_tuned_mean(n_type,:)))