clear all; clc;close all

files = dir('*WAKE.mat')
psth_spikes_PyN = [];
psth_spikes_nInN = [];

for ii = 1:length(files)
    load(files(ii).name)
   
    figure;
    colormap('jet')
    imagesc(zscore(squeeze(sum(PETHSpikes.PyN.counts,1))',[],2))
    ylabel('Cell Number')
    xlabel('Seconds from Triggers')
    title('Average repsonse per cell across all triggers')
    tickplaces = get(gca,'XTick');
    secondticks = tickplaces*PETHSpikes.PyN.parameters.binWidth-PETHSpikes.PyN.parameters.secondsBefore;
    ticktext = cellstr(num2str(secondticks'));
    set(gca,'XTickLabel',ticktext)
    hold on
    yl = ylim;
    zx = find(secondticks==0);%where zero is
    zx = tickplaces(zx);
    if ~isempty(zx)
        plot([zx zx],yl,'k')
    end
    data = squeeze(sum(PETHSpikes.PyN.counts,1));
    psth_spikes_PyN =[ psth_spikes_PyN, data];
    data = squeeze(sum(PETHSpikes.nInN.counts,1))
    psth_spikes_nInN =[ psth_spikes_nInN, data];
end


%%
firing_rates = psth_spikes_PyN';
Ncell = size(firing_rates, 1);
N2 = size(firing_rates, 2)/100;
windowed_data = reshape(firing_rates, [Ncell, N2, 100]);
% Initialize an empty matrix to store the correlation coefficients
pair_corr = zeros(N2,1);

% Loop over each window
for i = 1:size(windowed_data, 2)
    % Extract the data for this window
    window_data = squeeze(windowed_data(:, i, :));
    
    % Calculate the pairwise correlation for each pair of neurons
    cell_corr = corr(window_data');
    tmp= cell_corr;  % Create a copy of the matrix
    tmp(eye(size(cell_corr)) == 1) = NaN;  % Set the diagonal elements to NaN
    pair_corr(i) = nanmean(tmp(:));  % Calculate the mean of the non-diagonal elements
end




