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\HPC\57\day1\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day1\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day1\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day1\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day1\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day2\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day2\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day2\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day2\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day2\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day3\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day3\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day3\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day3\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day3\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day4\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day4\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day4\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day4\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day4\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day5\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day5\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day5\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day5\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day5\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day6\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day6\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day6\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day6\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day6\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day7\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day7\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day7\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day7\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day7\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day8\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day8\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day8\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day8\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day8\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day9\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\221\day9\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\851\day9\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\887\day9\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\987\day9\XYtotalspike1.csv';...
 'D:\ACC Data\HPC\57\day1\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day1\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day1\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day1\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day1\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day2\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day2\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day2\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day2\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day2\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day3\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day3\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day3\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day3\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day3\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day4\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day4\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day4\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day4\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day4\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day5\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day5\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day5\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day5\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day5\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day6\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day6\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day6\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day6\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day6\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day7\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day7\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day7\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day7\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day7\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day8\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day8\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day8\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day8\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day8\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\57\day9\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\221\day9\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\851\day9\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\887\day9\XYtotalspike3.csv';...
 'D:\ACC Data\HPC\987\day9\XYtotalspike3.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_basic_hpcacc;

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




%% 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),2,'filled','b');
scatter(tsne_results_place(find(all_cell_ids==2),1),tsne_results_place(find(all_cell_ids==2),2),5,'filled','g');



% 1. Setup figure and parameters
groups_to_plot = [1, 2];
group_names = {'ACC (1)', 'HPC (2)'};
figure('Color', 'w', 'Position', [100 100 1500 450]);

% Loop through each group to create subplots
for i = 1:2
    current_group = groups_to_plot(i);
    subplot(1, 2, 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');xlim([-65 65]);ylim([-65 65]);
    
    
    if i == 3
        cb = colorbar;
        cb.Label.String = 'Firing Position (Bins/cm)';
    end
end



rng(0,'twister');
%tsne_results = tsne(data_for_tsne, 'Perplexity', 40, 'Exaggeration',40,'LearnRate',2000,'Standardize', false);
tsne_results = tsne(data_for_tsne, 'Perplexity', 70, 'Exaggeration',40,'LearnRate',4000,'Standardize', false);
figure;hold on;

scatter(tsne_results(find(all_cell_ids==2),1),tsne_results(find(all_cell_ids==2),2),2,'filled','b');
scatter(tsne_results(find(all_cell_ids==1),1),tsne_results(find(all_cell_ids==1),2),2,'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 == 2, 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';



% 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');
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,:)))