addpath('C:\Users\Steve\Desktop\Neuro_Utilities\Matlab Scripts\Hayashi Lab Code');

clear all; clc; close all;

filelist={'D:\ACC Data\ACC\AllDays\635\day1\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day1\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day1\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day1\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day2\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day2\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day2\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day2\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day3\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day3\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day3\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day3\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day4\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day4\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day4\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day4\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day5\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day5\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day5\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day5\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day6\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day6\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day6\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day6\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day7\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day7\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day7\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day7\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day8\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day8\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day8\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day8\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day9\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\681\day9\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\691\day9\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\747\day9\XYtotalspike3.csv';...
 'D:\ACC Data\ACC\AllDays\635\day1\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day1\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day1\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day1\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day2\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day2\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day2\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day2\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day3\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day3\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day3\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day3\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day4\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day4\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day4\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day4\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day5\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day5\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day5\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day5\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day6\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day6\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day6\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day6\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day7\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day7\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day7\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day7\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day8\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day8\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day8\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day8\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\635\day9\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\681\day9\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\691\day9\XYtotalspike1.csv';...
 'D:\ACC Data\ACC\AllDays\747\day9\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day2\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day2\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day2\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day3\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day3\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day3\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day4\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day4\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day4\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day5\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day5\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day5\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day6\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day6\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day6\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day7\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day7\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day7\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day8\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day8\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day8\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day9\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day9\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day9\XYtotalspike3.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day2\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day2\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day2\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day3\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day3\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day3\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day4\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day4\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day4\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day5\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day5\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day5\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day6\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day6\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day6\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day7\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day7\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day7\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day8\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day8\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day8\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\431\supression\day9\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\546\supression\day9\XYtotalspike1.csv';...
 'D:\ACC Data\Ripple Suppression\551\supression\day9\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day2\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day2\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day2\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day3\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day3\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day3\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day4\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day4\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day4\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day5\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day5\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day5\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day6\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day6\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day6\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day7\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day7\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day7\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day8\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day8\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day8\XYtotalspike3.csv';... 
 'D:\ACC Data\DelayedRippleStim\173\day9\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day9\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day9\XYtotalspike3.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day2\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day2\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day2\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day3\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day3\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day3\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day4\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day4\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day4\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day5\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day5\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day5\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day6\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day6\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day6\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day7\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day7\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day7\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\173\day8\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day8\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day8\XYtotalspike1.csv';... 
 'D:\ACC Data\DelayedRippleStim\173\day9\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\184\day9\XYtotalspike1.csv';...
 'D:\ACC Data\DelayedRippleStim\235\day9\XYtotalspike1.csv'};
 



Alldata4=[];All_firing_position=[];All_firing_rate=[];All_stability=[];All_SI=[];

for session=1:size(filelist,1);
[norm_ratemap, firingrate, firing_position, occupancy_per_cell, stability, spatial_info, best_direction] = ca_imag_to_ratemap_v3_directional(filelist{session});
Alldata4=[Alldata4,norm_ratemap];
All_firing_rate=[All_firing_rate;firingrate];
All_firing_position=[All_firing_position;firing_position];
All_stability=[All_stability;stability];
All_SI=[All_SI;spatial_info];
clear firing_position firingrate norm_ratemap;
end;

%[all_cell_ids] = create_cell_labels_new_version_basic2;
[all_cell_ids SCcellIDs] = create_cell_labels_new_version_basic3;


data=Alldata4;
% Clear out junk
bad=find(All_stability<0.3);
bad=[bad;find(sum(data,2)==0)];
bad=[bad;find(sum(data)<0.5)'];
bad=unique(bad);

data(:,bad)=[];
all_cell_ids(bad)=[];
All_SI(bad)=[];
All_firing_position(bad)=[];
All_firing_rate(bad)=[];
SCcellIDs(bad)=[];



%% 1. Compute Autocorrelograms
[num_bins, num_neurons] = size(data);
ac_matrix = zeros(num_neurons, 143); % Observations x Features

fprintf('Computing autocorrelograms...\n');
for nm = 1:num_neurons
    rm = data(:, nm);
    if sum(rm) > 0 && ~any(isnan(rm))
        [ac, ~] = xcorr(rm, 'coeff');
        ac_matrix(nm, :) = ac;
    else
        ac_matrix(nm, :) = NaN;
    end
end

% Remove neurons with no activity (NaNs)
valid_idx = find(~any(isnan(ac_matrix), 2));
ac_clean = ac_matrix(valid_idx, :);
%ac_clean=ac_clean(:,32:112);


%% 2. Dimensionality Reduction (PCA)
% We reduce 143 bins to 30 Principal Components to denoise
fprintf('Running PCA...\n');
[coeff, score, latent, tsquared, explained] = pca(ac_clean);
num_pcs = 30; 
data_for_tsne = score(:, 1:num_pcs);

rng(0,'twister');
tsne_results_place = tsne(data', 'Perplexity', 40, 'Exaggeration',40,'LearnRate',2000,'Standardize', false);

figure;hold on;

scatter(tsne_results_place(find(all_cell_ids==1),1),tsne_results_place(find(all_cell_ids==1),2),5,'filled','b');
scatter(tsne_results_place(find(all_cell_ids==6),1),tsne_results_place(find(all_cell_ids==6),2),5,'filled','g');
scatter(tsne_results_place(find(all_cell_ids==7),1),tsne_results_place(find(all_cell_ids==7),2),5,'filled','r');


% 1. Setup figure and parameters
groups_to_plot = [1, 6, 7];
group_names = {'Control (1)', 'Stimulated (6)', 'Delayed (7)'};
figure('Color', 'w', 'Position', [100 100 1500 450]);

% Loop through each group to create subplots
for i = 1:3
    current_group = groups_to_plot(i);
    subplot(1, 3, i);
    hold on;
    
    % 2. Plot "Background" (All neurons in light grey)
    % This provides context for the overall manifold shape
    scatter(tsne_results_place(:,1), tsne_results_place(:,2), 2, ...
        [0.9 0.9 0.9], 'filled', 'MarkerFaceAlpha', 0.1);
    
    % 3. Extract indices for the current group
    group_idx = (all_cell_ids == current_group);
    
    % 4. Plot Group-Specific neurons colored by firing position
    % We use All_firing_position to define the 'C' (color) input
    scatter(tsne_results_place(group_idx, 1), tsne_results_place(group_idx, 2), ...
        6, All_firing_position(group_idx), 'filled');
    
    % 5. Formatting
    colormap(jet); % Jet is standard for circular/spatial position
    caxis([min(All_firing_position) max(All_firing_position)]); % Ensure scale is consistent
    title(group_names{i});
    xlabel('t-SNE 1'); ylabel('t-SNE 2');
    axis tight;
    
    if i == 3
        cb = colorbar;
        cb.Label.String = 'Firing Position (Bins/cm)';
    end
end



rng(0,'twister');
tsne_results = tsne(data_for_tsne, 'Perplexity', 38, 'Exaggeration',40,'LearnRate',2000,'Standardize', false);%tsne_results = tsne(data_for_tsne, 'Perplexity', 40, 'Exaggeration',40,'LearnRate',2000,'Standardize', false);

figure;hold on;
scatter(tsne_results(find(all_cell_ids==1),1),tsne_results(find(all_cell_ids==1),2),5,'filled','b');
scatter(tsne_results(find(all_cell_ids==6),1),tsne_results(find(all_cell_ids==6),2),5,'filled','g');
scatter(tsne_results(find(all_cell_ids==7),1),tsne_results(find(all_cell_ids==7),2),5,'filled','r');


figure;hold on;
scatter(tsne_results(find(SCcellIDs==0),1),tsne_results(find(SCcellIDs==0),2),2,'filled','b');
scatter(tsne_results(find(SCcellIDs==1),1),tsne_results(find(SCcellIDs==1),2),5,'filled','r');



% 1. Setup Data
% tsne_out: 9066x2 matrix
num_total = size(tsne_results, 1);
num_pop1 = find(all_cell_ids == 1, 1);
num_pop2 = find(all_cell_ids == 6, 1);

% Create a label vector: 1 for Population 1, 0 for Population 2
labels = zeros(num_total, 1);
labels(num_pop1:num_pop2) = 1;

% 2. Find 50 Nearest Neighbors
% 'K', 51 because the point itself is usually included as the first neighbor
k = 150;
[indices, ~] = knnsearch(tsne_results, tsne_results, 'K', k + 1);

% Remove the first column (the point itself) to get the true neighbors
neighbor_indices = indices(:, 2:end);

% 3. Calculate Enrichment Score
% Extract labels of the neighbors for every point
neighbor_labels = labels(neighbor_indices);

% Calculate the percentage of neighbors belonging to Population 1
% enrichment_score will range from 0 to 1
enrichment_score = mean(neighbor_labels, 2);

% 4. Plotting
figure('Color', 'w');
scatter(tsne_results(:,1), tsne_results(:,2), 3, enrichment_score, 'filled');

% Apply Jet colormap as requested
colormap(jet);
colorbar;
title('Population 1 Neighborhood Enrichment');
xlabel('t-SNE 1');
ylabel('t-SNE 2');

% Label the colorbar for clarity
c = colorbar;
c.Label.String = 'Proportion of Population 1 Neighbors';











%%

% 1. Setup Data
num_total = size(tsne_results, 1);

% Create a binary label vector: 1 for Population 6, 0 for everything else
% This ensures we only count the stimulated group
labels_pop6 = (all_cell_ids == 6);

% 2. Find Nearest Neighbors
% Using k = 150 as in your previous version for consistency
k = 150;
fprintf('Calculating nearest neighbors for %d neurons...\n', num_total);
[indices, ~] = knnsearch(tsne_results, tsne_results, 'K', k + 1);

% Remove the first column (the point itself)
neighbor_indices = indices(:, 2:end);

% 3. Calculate Enrichment Score for Population 6
% Extract labels of the neighbors for every point
neighbor_labels = labels_pop6(neighbor_indices);

% Calculate the proportion of neighbors belonging to Population 6
% enrichment_score_pop6 will range from 0 (no Pop 6) to 1 (only Pop 6)
enrichment_score_pop6 = mean(neighbor_labels, 2);

% 4. Plotting
figure('Color', 'w', 'Position', [100 100 850 700]);
% Plot all neurons, colored by the local density of Population 6
scatter(tsne_results(:,1), tsne_results(:,2), 4, enrichment_score_pop6, 'filled');

% Apply Jet colormap
colormap(jet);
c = colorbar;
c.Label.String = 'Proportion of Population 6 Neighbors';
clim([0 0.3]); % Force scale from 0 to 100%

title('Population 6 (Stimulated) Neighborhood Enrichment');
subtitle('Red: High Concentration of Pop 6 | Blue: Depletion of Pop 6');
xlabel('t-SNE 1');
ylabel('t-SNE 2');
axis tight;
grid off;

%%












% Circle neurons

[isInB, indexInB] = ismember(brushedData, tsne_results, 'rows');
figure;imagesc(data(:,indexInB)); colormap(jet)




% 1. Extract the coordinates of the points in your cluster
cluster_coords = tsne_results(indexInB, :);

% 2. Calculate the Convex Hull indices
% 'k' will be the indices of points that form the "corners" of the polygon
k = convhull(cluster_coords(:,1), cluster_coords(:,2));

% 3. Plot the boundary on your t-SNE
hold on; % Ensure you are plotting on top of your existing scatter
hPoly = plot(cluster_coords(k,1), cluster_coords(k,2), 'k-', 'LineWidth', 0.5);

% Optional: Fill the polygon with a light color to make it stand out
fill(cluster_coords(k,1), cluster_coords(k,2), 'k', 'FaceAlpha', 0.05, 'EdgeColor', 'none');



tempmaps=data(:,indexInB);
% 1. Find the index of the maximum value for each column (neuron)
[~, peak_indices] = max(tempmaps, [], 1);

% 2. Sort those indices in ascending order
[~, sort_order] = sort(peak_indices);

% 3. Reorder the columns of your matrix
sorted_maps = tempmaps(:, sort_order);

%% Optional: Visualization
figure;
imagesc(sorted_maps'); % Transpose so neurons are rows and bins are columns
colormap(jet);
xlabel('Spatial Bins (1-72)');
ylabel('Neurons (Sorted)');
title('Sorted Population Ratemap');
colorbar;



[isInB2, indexInB2] = ismember(brushedData1, tsne_results, 'rows');
figure;imagesc(data(:,indexInB2)); colormap(jet)


tempmaps=data(:,indexInB2);
% 1. Find the index of the maximum value for each column (neuron)
[~, peak_indices] = max(tempmaps, [], 1);

% 2. Sort those indices in ascending order
[~, sort_order] = sort(peak_indices);

% 3. Reorder the columns of your matrix
sorted_maps = tempmaps(:, sort_order);

%% Optional: Visualization
figure;
imagesc(sorted_maps'); % Transpose so neurons are rows and bins are columns
colormap(jet);
xlabel('Spatial Bins (1-72)');
ylabel('Neurons (Sorted)');
title('Sorted Population Ratemap');
colorbar;




















[isInB2, indexInB2] = ismember(brushedData1, tsne_results, 'rows');
targetclu=all_cell_ids(indexInB);
remainderclu=all_cell_ids(indexInB2);


% 1. Extract Labels for the two groups
labels_target = all_cell_ids(indexInB);
labels_rem = all_cell_ids(indexInB2);

% 2. Calculate Counts for 1, 6, and 7
% [Pop 1, Pop 6, Pop 7]
counts_target = [sum(labels_target == 1), sum(labels_target == 6), sum(labels_target == 7)];
counts_rem = [sum(labels_rem == 1), sum(labels_rem == 6), sum(labels_rem == 7)];

% 3. Normalize to percentages for plotting
perc_target = (counts_target / sum(counts_target)) * 100;
perc_rem = (counts_rem / sum(counts_rem)) * 100;
data_to_plot = [perc_target; perc_rem];

% 4. Grouped Bar Chart
figure('Color', 'w');
b = bar(data_to_plot, 'grouped');
set(gca, 'XTickLabel', {'Target Cluster (indexInB)', 'Remainder (indexInB2)'});
ylabel('Percentage of Population (%)');
legend({'Control (1)', 'Stimulated (6)', 'Delayed (7)'}, 'Location', 'northeast');
title('Distribution of Populations across Clusters');
grid on;

% 5. Stacked Bar Chart (Composition)
figure('Color', 'w');
bar(data_to_plot, 'stacked');
set(gca, 'XTickLabel', {'Target Cluster', 'Remainder'});
ylabel('Cumulative Percentage (%)');
legend({'Control (1)', 'Stimulated (6)', 'Delayed (7)'});
title('Cluster Composition Comparison');

% 6. Statistics: Chi-Square Test for 3x2 Contingency Table
% Rows: Pop 1, Pop 6, Pop 7 | Columns: Target, Remainder
contingency_table = [counts_target', counts_rem'];

% Calculate Chi-Square
observed = contingency_table;
row_totals = sum(observed, 2);
col_totals = sum(observed, 1);
grand_total = sum(observed(:));
expected = (row_totals * col_totals) / grand_total;

chi2_stat = sum((observed(:) - expected(:)).^2 ./ expected(:));
df = (size(observed, 1) - 1) * (size(observed, 2) - 1);
p_val = 1 - chi2cdf(chi2_stat, df);

% Print Results
fprintf('\n--- Cluster Comparison Statistics ---\n');
fprintf('Target Cluster Counts: [1: %d, 6: %d, 7: %d]\n', counts_target);
fprintf('Remainder Cluster Counts: [1: %d, 6: %d, 7: %d]\n', counts_rem);
fprintf('Chi-Square Statistic: %.4f\n', chi2_stat);
fprintf('Degrees of Freedom: %d\n', df);
fprintf('p-value: %e\n', p_val);

if p_val < 0.05
    fprintf('Result: Significant difference in population distribution between clusters.\n');
else
    fprintf('Result: No significant difference found.\n');
end


figure;plot(mean(ac_clean(indexInB,:)))
hold on;plot(mean(ac_clean(indexInB2,:)))



save 'D:\ACC Data\ACC\WorkingACCTSNEOUTPUT.mat'

% Fig 1j stats
% 1. Determine Center and Lags
% If ac_clean is 73 bins, the center (lag 0) is index 37.
num_bins = size(ac_clean, 2);
center_idx = ceil(num_bins/2); 
lags_vector = (0:(num_bins - center_idx)) * 2; % 2cm per bin

% 2. Extract Widths for both groups
% We'll use a helper to find the half-width
get_hwhm = @(profile) find(profile(center_idx:end) <= 0.5, 1, 'first');

widths_c1 = [];
for i = 1:length(indexInB)
    idx = get_hwhm(ac_clean(indexInB(i), :));
    if ~isempty(idx)
        widths_c1(i,1) = lags_vector(idx);
    else
        widths_c1(i,1) = max(lags_vector); % Cap at max
    end
end

widths_rem = [];
for i = 1:length(indexInB2)
    idx = get_hwhm(ac_clean(indexInB2(i), :));
    if ~isempty(idx)
        widths_rem(i,1) = lags_vector(idx);
    else
        widths_rem(i,1) = max(lags_vector);
    end
end

% 3. Simple Statistical Test (Mann-Whitney U)
[p_val, ~, stats] = ranksum(widths_c1, widths_rem);
fprintf('Median Widths - Cluster 1: %.1f cm, Remainder: %.1f cm\n', median(widths_c1), median(widths_rem));
fprintf('Mann-Whitney p-value: %e\n', p_val);

% 4. Nature-Style Violin Plot
figure('Color', 'w', 'Name', 'Spatial Tuning Widths');
data_combined = [widths_rem; widths_c1];
group_idx = [zeros(length(widths_rem),1); ones(length(widths_c1),1)];

% If you don't have violinplot.m, a boxplot with swarm is a great alternative:
boxplot(data_combined, group_idx, 'Labels', {'Remainder', 'Cluster 1'}, 'Symbol', '');
hold on;
% Add a bit of jittered raw data points to be transparent (Nature loves raw data)
scatter(group_idx+1 + (rand(size(group_idx))-0.5)*0.2, data_combined, 5, [0.5 0.5 0.5], 'filled', 'MarkerFaceAlpha', 0.1);
ylabel('Field Half-Width (cm)');
title('Spatial Scale Divergence');
grid on;





%extra
% 1. Identify the center of the AC matrix
[num_neurons, num_bins] = size(ac_clean);
center_bin = ceil(num_bins / 2);

% 2. Initialize the width vector
field_widths = zeros(num_neurons, 1);

% 3. Loop through each neuron to find the HWHM
for i = 1:num_neurons
    % Get the right half of the AC starting from the peak
    right_half = ac_clean(i, center_bin:end);
    
    % Find the first bin where the value drops below 0.5
    % We assume the peak is normalized to 1.0
    half_max_idx = find(right_half < 0.5, 1, 'first');
    
    if ~isempty(half_max_idx)
        % The index represents the number of bins from the center
        field_widths(i) = half_max_idx - 1; 
    else
        % If it never drops below 0.5, assign the max possible width
        field_widths(i) = length(right_half);
    end
end

% 4. Convert bins to spatial units (cm) if necessary
% Example: if each bin is 2cm
% field_widths_cm = field_widths * 2;

% 4. Plotting
figure('Color', 'w');
scatter(tsne_results(:,1),tsne_results(:,2), 5, field_widths, 'filled');

% Apply Jet colormap as requested
colormap(jet);
colorbar;


figure;hold on;

scatter(tsne_results_place(find(all_cell_ids==1),1),tsne_results_place(find(all_cell_ids==1),2),2,'filled','b');
scatter(tsne_results_place(find(all_cell_ids==6),1),tsne_results_place(find(all_cell_ids==6),2),5,'filled','g');
scatter(tsne_results_place(find(all_cell_ids==7),1),tsne_results_place(find(all_cell_ids==7),2),2,'filled','r');


% 4. Plotting
figure('Color', 'w');
scatter(tsne_results_place(find(all_cell_ids==1),1),tsne_results_place(find(all_cell_ids==1),2), 3, field_widths(find(all_cell_ids==1)), 'filled');

% Apply Jet colormap as requested
colormap(jet);
colorbar;
title('Population 1 Neighborhood Enrichment');
xlabel('t-SNE 1');
ylabel('t-SNE 2');

% Label the colorbar for clarity
c = colorbar;
c.Label.String = 'Proportion of Population 1 Neighbors';


% 4. Plotting
figure('Color', 'w');
scatter(tsne_results_place(find(all_cell_ids==6),1),tsne_results_place(find(all_cell_ids==6),2), 5, field_widths(find(all_cell_ids==6)), 'filled');

% Apply Jet colormap as requested
colormap(jet);
colorbar;
title('Population 1 Neighborhood Enrichment');
xlabel('t-SNE 1');
ylabel('t-SNE 2');

% Label the colorbar for clarity
c = colorbar;
c.Label.String = 'Proportion of Population 1 Neighbors';


% 4. Plotting
figure('Color', 'w');
scatter(tsne_results_place(find(all_cell_ids==7),1),tsne_results_place(find(all_cell_ids==7),2), 5, field_widths(find(all_cell_ids==7)), 'filled');

% Apply Jet colormap as requested
colormap(jet);
colorbar;
title('Population 1 Neighborhood Enrichment');
xlabel('t-SNE 1');
ylabel('t-SNE 2');

% Label the colorbar for clarity
c = colorbar;
c.Label.String = 'Proportion of Population 1 Neighbors';





% 4. Plotting
figure('Color', 'w');
scatter(tsne_results_place(find(all_cell_ids==1),1),tsne_results_place(find(all_cell_ids==1),2), 5, All_firing_position(find(all_cell_ids==1)), 'filled');
colormap(jet);
colorbar;xlim([-60 60]);ylim([-60 60]);

figure('Color', 'w');
scatter(tsne_results_place(find(all_cell_ids==6),1),tsne_results_place(find(all_cell_ids==6),2), 5, All_firing_position(find(all_cell_ids==6)), 'filled');
colormap(jet);
colorbar;xlim([-60 60]);ylim([-60 60]);

figure('Color', 'w');
scatter(tsne_results_place(find(all_cell_ids==7),1),tsne_results_place(find(all_cell_ids==7),2), 5, All_firing_position(find(all_cell_ids==7)), 'filled');
colormap(jet);
colorbar;xlim([-60 60]);ylim([-60 60]);


hold on;
scatter(tsne_results_place(find(all_cell_ids==6),1),tsne_results_place(find(all_cell_ids==6),2), 3, All_firing_position(find(all_cell_ids==6)), 'filled');
scatter(tsne_results_place(find(all_cell_ids==7),1),tsne_results_place(find(all_cell_ids==7),2), 5, All_firing_position(find(all_cell_ids==7)), 'filled');

% Apply Jet colormap as requested
colormap(jet);
colorbar;
title('Population 1 Neighborhood Enrichment');
xlabel('t-SNE 1');
ylabel('t-SNE 2');

% Label the colorbar for clarity
c = colorbar;
c.Label.String = 'Proportion of Population 1 Neighbors';