plot_execution_latency.m

% Plot execution latency
clc
clear all
close all

% parameters
spike_sorting_type = '_unsorted_aligned_thr_-4.5';
taskName = 'GraspObject_4S_Action';
%taskName = 'GraspObject_5050';
subject_id = 's4';
classifierType = 'diaglinear';

% Data = load('C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\s3\Data\IndividualFiles\GraspObject\unsorted_aligned_thr_-4.5\s3_20230803_unsorted_aligned_thr_-4.5_GraspObject');
%Data = load('C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\s3\Data\IndividualFiles\GraspObject\unsorted_aligned_thr_-4.5\s3_20230724_unsorted_aligned_thr_-4.5_GraspObject');
% Data = load('C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\s3\Data\IndividualFiles\GraspObject\unsorted_aligned_thr_-4.5\s3_20230721_unsorted_aligned_thr_-4.5_GraspObject');

Data = load(['C:\Users\macthurston\OneDrive - Kaiser Permanente\CaltechData\GraspObject_project\' subject_id '\Data\Table_' subject_id '_' taskName spike_sorting_type]);

% keyboard

Go_data = Data.Go_data;

% 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

brainAreas = Go_data.frPerChannel{7}; % 7 for Ripple, 6 for Blackrock
phase_time_idx = Go_data.time_phase_labels{1,1};
numPhases = numel(unique(phase_time_idx));
phaseTimeTmp = diff(phase_time_idx);
phase_changes(1) = 1;
phase_changes(2:numPhases) = find(phaseTimeTmp) + 1;
phaseNames = {'ITI', 'Planning', 'Delay', 'Action'};
sessions_all = unique(Go_data.session_date);
numSessions = numel(sessions_all);
uniqueCueTypes = {'Hand','Hand-Object','Object'}; 
n_regions = length(brainAreas);
num_timebins = numel(phase_time_idx);

%% load analyzed data
directory = ['C:\Users\macthurston\Documents\GitHub\project_grasp_object_interaction\analyzedData\' subject_id];
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored.mat"; 
full_path = fullfile(directory, filename);
load(full_path);

% Go/NoGo Classification: every region x every session %
G_all_regions = all_session_acc;

% SHUFFLED %
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored_shuffled.mat"; 
full_path = fullfile(directory, filename);
load(full_path);

% Go/NoGo Classification: every region x every session %
S_all_regions = all_session_acc;

%% plot
ci_bounds = [2.5 97.5];
for r = 1:n_regions
    G_regionData = G_all_regions(r,:);
    S_regionData = S_all_regions(r,:);
    figure('units','normalized','outerposition',[0 0 0.25 1])
    sgtitle(['Go/NoGo Classification Over Time - ' brainAreas{r}]);
    for s = 1:numSessions
        % Go-NoGo %
        sessionData = cell2mat(G_regionData(s));
        mean_acc = mean(sessionData, 1);
        ci_low  = prctile(sessionData, ci_bounds(1), 1);
        ci_high = prctile(sessionData, ci_bounds(2), 1);
        subplot(numSessions, 1, s); 
        hold on;
        
        x = 1:num_timebins;
        
        % shaded CI
        fill([x fliplr(x)], ...
             [ci_low fliplr(ci_high)], ...
             [0.8 0.8 1], 'EdgeColor', 'none');
        
        % mean line
        plot(x, mean_acc, 'b', 'LineWidth', 2);

        % SHUFFLED %
        S_sessionData = cell2mat(S_regionData(s));
        S_mean_acc = mean(S_sessionData, 1);
        S_ci_low  = prctile(S_sessionData, ci_bounds(1), 1);
        S_ci_high = prctile(S_sessionData, ci_bounds(2), 1);
        subplot(numSessions, 1, s); 
        hold on;

        x = 1:num_timebins;

        % shaded CI
        fill([x fliplr(x)], ...
             [S_ci_low fliplr(S_ci_high)], ...
             [0.8 0.8 0.8], 'EdgeColor', 'none');

        % mean line
        plot(x, S_mean_acc, 'Color', [0.5,0.5,0.5], 'LineWidth', 2);

        % Add phase lines
        for n_phase = 1:numPhases
            xline(phase_changes(n_phase), 'k--', phaseNames{n_phase}, 'LineWidth', 1.5,'FontSize',9);
        end
        
        xlim([0 num_timebins + 1]);
        ylim([0 100]);
        %yline(50,'--r','LineWidth',1.5);
        yticks([0 25 50 75 100]);
        xlabel('Timebin');
        ylabel('Classification Accuracy [%]');
        title([sessions_all{s}]);
        set(gca, 'FontSize', 12);
        hold off;
    end
end

%% onset time

first_sig_bin_all = NaN(n_regions,numSessions);
effect_size_all = NaN(n_regions,numSessions);

for r = 1:n_regions

    G_regionData = G_all_regions(r,:);
    S_regionData = S_all_regions(r,:);

    for s = 1:numSessions
    
        boot_real = G_regionData{s};
        boot_null = S_regionData{s};
    
        if isempty(boot_real); continue; end
    
        pvals = NaN(1, num_timebins);
    
        for t = 1:num_timebins
            pvals(t) = ranksum(boot_real(:,t), boot_null(:,t), 'tail', 'right');
        end
    
        sig = mean(boot_real > boot_null, 1) > 0.975;
        first_sig_idx = find(sig, 1, 'first');

        if isempty(first_sig_idx)
            first_sig_bin = NaN;
        else
            first_sig_bin = first_sig_idx;
        end
    
        effect = mean(boot_real,1) - mean(boot_null,1);
        first_effect_idx = effect(first_sig_idx);

        if isempty(first_effect_idx)
            first_effect_size = NaN;
        else
            first_effect_size = first_effect_idx;
        end
    
        % % consecutive bins criterion
        % k = 3;
        % sig_conv = conv(double(sig), ones(1,k), 'same') >= k;
        % 
        first_sig_bin_all(r,s) = first_sig_bin;
        effect_size_all(r,s) = first_effect_size;
    end
end

%% box plots
figure('units','normalized','outerposition',[0 0 0.27 0.25])
%color_info = [0.3359, 0.7031, 0.9101;0.8984 0.6211 0;0.8320 0.3672 0;0.7969, 0.4726, 0.6523;0, 0.6171, 0.4492;.9961 .6875 0]; % SMG, PMV, S1, AIP, M1, dlPFC

% % AN % %
% groupOrder = ["S1","AIP","PMV","SMG","M1"];
% colorOrder = [0.8320 0.3672 0;0.7969, 0.4726, 0.6523;0.8984 0.6211 0;0.3359, 0.7031, 0.9101;0, 0.6171, 0.4492]; % SMG,PMV,S1,AIP,M1
% % GB % %
% groupOrder = ["dlPFC","AIP","SMG","S1","M1","PMV"];
% colorOrder = [.9961 .6875 0;0.7969, 0.4726, 0.6523;0.3359, 0.7031, 0.9101;0.8320 0.3672 0; 0, 0.6171, 0.4492; 0.8984 0.6211 0]; % AIP,SMG,M1,S1,dlPFC
% % FG % %
groupOrder = ["PMV","SMG","S1"];
colorOrder = [0.8984 0.6211 0;0.3359, 0.7031, 0.9101;0.8320 0.3672 0]; % SMG,PMV,S1
%sgtitle('Planning Stats');
[~, reorderedIdx] = ismember(brainAreas, groupOrder);

% Planning onset timing
%subplot(2,2,1);
h = boxplot(first_sig_bin_all', brainAreas(1:3),"Orientation", "horizontal", "BoxStyle","outline","Colors",colorOrder, "GroupOrder",groupOrder); %, 'Notch', 'on', 'Labels', {'SMG','PMv','S1','AIP','M1'})
hold on
% Scatter plot of individual session data points
for n_region = 1:3%numel(brainAreas)
    region_bins = first_sig_bin_all';
    x_vals = region_bins(:, n_region);   
    y_vals = repmat(reorderedIdx(n_region), length(x_vals), 1); 
    
    scatter(x_vals, y_vals, 10, 'k', 'o', ...
            'MarkerFaceAlpha', 0.2);
end
title('Execution Onset Response');
set(h,{'linew'},{1.5})
xlim([94 174]);
xticks([94 104 114 124 134 144 154 164 174]);
xticklabels([0 0.5 1 1.5 2 2.5 3 3.5 4]);
xlabel('Time (s)');
ylabel('Brain Region');

%% sig testing between regions
p_mat = NaN(n_regions);

for r1 = 1:n_regions
    for r2 = r1+1:n_regions

        x = first_sig_bin_all(r1,:);
        y = first_sig_bin_all(r2,:);

        % paired removal of NaNs
        valid = ~isnan(x) & ~isnan(y);
        x = x(valid);
        y = y(valid);

        if isempty(x)
            continue
        end

        p = signrank(x, y);

        p_mat(r1,r2) = p;
        p_mat(r2,r1) = p;
    end
end

% --- FDR correction ---
mask = triu(true(n_regions),1);
p_vec = p_mat(mask);

% remove NaNs before FDR
valid_idx = ~isnan(p_vec);
p_vec_valid = p_vec(valid_idx);

p_fdr_valid = mafdr(p_vec_valid, 'BHFDR', true);

% reconstruct full matrix
p_mat_fdr = NaN(size(p_mat));
p_vec_fdr = NaN(size(p_vec));
p_vec_fdr(valid_idx) = p_fdr_valid;

p_mat_fdr(mask) = p_vec_fdr;
p_mat_fdr = p_mat_fdr + p_mat_fdr'; % symmetric

% --- Effect sizes (median differences) ---
effect_mat = NaN(n_regions);

for r1 = 1:n_regions
    for r2 = r1+1:n_regions

        x = first_sig_bin_all(r1,:);
        y = first_sig_bin_all(r2,:);

        valid = ~isnan(x) & ~isnan(y);
        x = x(valid);
        y = y(valid);

        if isempty(x)
            continue
        end

        diff_xy = x - y;
        med_diff = median(diff_xy);

        effect_mat(r1,r2) = med_diff;
        effect_mat(r2,r1) = -med_diff;
    end
end

%% plot all subjects together 
% load S2
% load DATA 
subject_id = 's2';
directory = ['C:\Users\macthurston\Documents\GitHub\project_grasp_object_interaction\analyzedData\' subject_id];
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored.mat"; 
full_path = fullfile(directory, filename);
load(full_path);
n_regions = length(brainAreas(1:3));
numSessions = size(all_session_acc, 2);

% Go/NoGo Classification: every region x every session %
G_all_regions = all_session_acc;

% SHUFFLED %
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored_shuffled.mat"; 
full_path = fullfile(directory, filename);
load(full_path);

% Go/NoGo Classification: every region x every session %
S_all_regions = all_session_acc;

first_sig_bin_all = NaN(n_regions,numSessions);
effect_size_all = NaN(n_regions,numSessions);

for r = 1:n_regions

    G_regionData = G_all_regions(r,:);
    S_regionData = S_all_regions(r,:);

    for s = 1:numSessions
    
        boot_real = G_regionData{s};
        boot_null = S_regionData{s};
    
        if isempty(boot_real); continue; end
    
        pvals = NaN(1, num_timebins);
    
        for t = 1:num_timebins
            pvals(t) = ranksum(boot_real(:,t), boot_null(:,t), 'tail', 'right');
        end
    
        sig = mean(boot_real > boot_null, 1) > 0.975;
        first_sig_idx = find(sig, 1, 'first');

        if isempty(first_sig_idx)
            first_sig_bin = NaN;
        else
            first_sig_bin = first_sig_idx;
        end
    
        effect = mean(boot_real,1) - mean(boot_null,1);
        first_effect_idx = effect(first_sig_idx);

        if isempty(first_effect_idx)
            first_effect_size = NaN;
        else
            first_effect_size = first_effect_idx;
        end
    
        % % consecutive bins criterion
        % k = 3;
        % sig_conv = conv(double(sig), ones(1,k), 'same') >= k;
        % 
        first_sig_bin_all(r,s) = first_sig_bin;
        effect_size_all(r,s) = first_effect_size;
    end
end

% keep desired brain regions and reorder
S2_cue_onset = first_sig_bin_all;
S2_effect_size = effect_size_all;
%brainAreas(:,[2,4,5]) = []; % add/remove 6 for Ripple
S2_brainAreas = brainAreas(1:3);
S2_regions = length(S2_brainAreas);

%% load S3
% load DATA 
subject_id = 's3';
directory = ['C:\Users\macthurston\Documents\GitHub\project_grasp_object_interaction\analyzedData\' subject_id];
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored.mat"; 
full_path = fullfile(directory, filename);
load(full_path);
n_regions = length(brainAreas(1:5));
numSessions = size(all_session_acc, 2);

% Go/NoGo Classification: every region x every session %
G_all_regions = all_session_acc;

% SHUFFLED %
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored_shuffled.mat"; 
full_path = fullfile(directory, filename);
load(full_path);

% Go/NoGo Classification: every region x every session %
S_all_regions = all_session_acc;

first_sig_bin_all = NaN(n_regions,numSessions);
effect_size_all = NaN(n_regions,numSessions);

for r = 1:n_regions

    G_regionData = G_all_regions(r,:);
    S_regionData = S_all_regions(r,:);

    for s = 1:numSessions
    
        boot_real = G_regionData{s};
        boot_null = S_regionData{s};
    
        if isempty(boot_real); continue; end
    
        pvals = NaN(1, num_timebins);
    
        for t = 1:num_timebins
            pvals(t) = ranksum(boot_real(:,t), boot_null(:,t), 'tail', 'right');
        end
    
        sig = mean(boot_real > boot_null, 1) > 0.975;
        first_sig_idx = find(sig, 1, 'first');

        if isempty(first_sig_idx)
            first_sig_bin = NaN;
        else
            first_sig_bin = first_sig_idx;
        end
    
        effect = mean(boot_real,1) - mean(boot_null,1);
        first_effect_idx = effect(first_sig_idx);

        if isempty(first_effect_idx)
            first_effect_size = NaN;
        else
            first_effect_size = first_effect_idx;
        end
    
        % % consecutive bins criterion
        % k = 3;
        % sig_conv = conv(double(sig), ones(1,k), 'same') >= k;
        % 
        first_sig_bin_all(r,s) = first_sig_bin;
        effect_size_all(r,s) = first_effect_size;
    end
end

% keep desired brain regions and reorder
S3_cue_onset_tmp = first_sig_bin_all;
%brainAreas(:,2) = []; % remove PMV, add/remove 6 for Ripple
S3_cue_onset = S3_cue_onset_tmp([4,1,5,3,2],:);
S3_effect_size_tmp = effect_size_all;
S3_effect_size = S3_effect_size_tmp([4,1,5,3,2],:);
S3_brainAreas = brainAreas(:,[4,1,5,3,2]);
S3_regions = length(S3_brainAreas);

%% load S4
% load Data GRASP
subject_id = 's4';
directory = ['C:\Users\macthurston\Documents\GitHub\project_grasp_object_interaction\analyzedData\' subject_id];
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored.mat"; 
full_path = fullfile(directory, filename);
load(full_path);
n_regions = length(brainAreas);
numSessions = size(all_session_acc, 2);

% Go/NoGo Classification: every region x every session %
G_all_regions = all_session_acc;

% SHUFFLED %
filename = "GoNoGo_Classification_" + taskName + '_ALL_REGIONS' + "_LDAfit" + "_z_scored_shuffled.mat"; 
full_path = fullfile(directory, filename);
load(full_path);

% Go/NoGo Classification: every region x every session %
S_all_regions = all_session_acc;

first_sig_bin_all = NaN(n_regions,numSessions);
effect_size_all = NaN(n_regions,numSessions);

for r = 1:n_regions

    G_regionData = G_all_regions(r,:);
    S_regionData = S_all_regions(r,:);

    for s = 1:numSessions
    
        boot_real = G_regionData{s};
        boot_null = S_regionData{s};
    
        if isempty(boot_real); continue; end
    
        pvals = NaN(1, num_timebins);
    
        for t = 1:num_timebins
            pvals(t) = ranksum(boot_real(:,t), boot_null(:,t), 'tail', 'right');
        end
    
        sig = mean(boot_real > boot_null, 1) > 0.975;
        first_sig_idx = find(sig, 1, 'first');

        if isempty(first_sig_idx)
            first_sig_bin = NaN;
        else
            first_sig_bin = first_sig_idx;
        end
    
        effect = mean(boot_real,1) - mean(boot_null,1);
        first_effect_idx = effect(first_sig_idx);

        if isempty(first_effect_idx)
            first_effect_size = NaN;
        else
            first_effect_size = first_effect_idx;
        end
    
        % % consecutive bins criterion
        % k = 3;
        % sig_conv = conv(double(sig), ones(1,k), 'same') >= k;
        % 
        first_sig_bin_all(r,s) = first_sig_bin;
        effect_size_all(r,s) = first_effect_size;
    end
end

% keep desired brain regions and reorder
S4_cue_onset_tmp = first_sig_bin_all;
%brainAreas(:,2) = []; % remove PMV, add/remove 6 for Ripple
S4_cue_onset = S4_cue_onset_tmp([4,1,5,3,6],:);
S4_effect_size_tmp = effect_size_all;
S4_effect_size = S4_effect_size_tmp([4,1,5,3,6],:);
S4_brainAreas = brainAreas(:,[4,1,5,3,6]);
S4_regions = length(S4_brainAreas);

%% plot all together by region
participants = {S2_cue_onset, S3_cue_onset, S4_cue_onset};
participant_names = {'S2','S3','S4'};
participant_regions = {S2_brainAreas, S3_brainAreas, S4_brainAreas};
% all_regions = {'AIP','SMG','M1','PMV','S1','dlPFC'};
% color_info = {[0.7969, 0.4726, 0.6523],[0.3359, 0.7031, 0.9101],[0, 0.6171, 0.4492],[0.8984 0.6211 0],[0.8320 0.3672 0],[0 0 0]}; % SMG, AIP, M1, PMV, S1, dlPFC
all_regions = {'M1','S1','PMV','SMG','AIP','dlPFC'};
color_info = {[0, 0.6171, 0.4492],[0.8320 0.3672 0],[0.8984 0.6211 0],[0.3359, 0.7031, 0.9101],[0.7969, 0.4726, 0.6523],[0 0 0]}; % 'M1','PMV','S1','SMG','AIP','dlPFC'
brainAreas = Go_data.frPerChannel{7}; % 7 for Ripple, 6 for Blackrock
n_regions = length(all_regions);
values = [];
region_grp = [];
participant_grp = [];

% execution means (94 = timebins before Action phase; *50 bc 50 ms per
% timebin)
S2_mean = (nanmean(S2_cue_onset,2) - 94)*50; % order: SMG,PMv,S1
S3_mean = (nanmean(S3_cue_onset,2) - 94)*50; % AIP,SMG,M1,S1,PMV
S4_mean = (nanmean(S4_cue_onset,2) - 94)*50; % AIP,SMG,M1,S1,dlPFC

AIP_mean = (S3_mean(1) + S4_mean(1)) / 2;
SMG_mean = (S3_mean(2) + S2_mean(1)) / 2;
M1_mean= (S3_mean(3) + S4_mean(3)) / 2;
S1_mean = (S3_mean(4) + S4_mean(4)) / 2;
PMV_mean = S3_mean(5);
dlPFC_mean = S4_mean(5);

for p_i = 1:numel(participants)
    
    data = participants{p_i};           % sessions × regions
    regions_p = participant_regions{p_i};
    
    for r = 1:numel(regions_p)
        
        region_idx = find(strcmp(all_regions, regions_p{r}));
        
        y = data(:,r);   % all sessions for this participant + region
        
        n = numel(y);
        
        values = [values; y];
        region_grp = [region_grp; repmat(region_idx,n,1)];
        participant_grp = [participant_grp; repmat(p_i,n,1)];
        
    end
end

figure('units','normalized','outerposition',[0 0 0.362 0.36]); % 0 0 0.3 0.32
hold on

offset = 0.25;   % spacing between participants
width = 0.2;

for r = 1:numel(all_regions)
    
    for p_i = 1:numel(participants)
        
        regions_p = participant_regions{p_i};
        
        if any(strcmp(all_regions{r}, regions_p))
            
            % get data
            data = participants{p_i};
            row_idx = find(strcmp(regions_p, all_regions{r}));
            y = data(row_idx,:);   % time values
            
            % y-position (region index with participant offset)
            y_pos = r + (p_i - (numel(participants)+1)/2)*offset;
            
            b = boxchart(ones(size(y))*y_pos, y, ...
                'BoxWidth', width, ...
                'BoxFaceColor', color_info{r}, ...
                'Orientation','horizontal');

            % % outliers as red crosses
            % b.MarkerStyle = 'x';
            % b.MarkerColor = [1 0 0];
            % b.MarkerSize = 6;
            
        end
    end
end

title('Execution Onset Response');
%set(h,{'linew'},{1.5})
xlim([94 174]);
xticks([94 104 114 124 134 144 154 164 174]);
xticklabels([0 0.5 1 1.5 2 2.5 3 3.5 4]);
xlabel('Time from Go Cue [s]');
ylabel('Brain Region');
yticks([1 2 3 4 5 6]);
yticklabels([all_regions]);
ylabel('Brain regions')
set(gca,'YDir','reverse', 'FontSize', 20);