% routine to calculate the results for Figure3&4, requires the results from
% Ensemble_detection routine. it loads the raw data and results from the
% Ensemble_detection and calculates the cross/auto correlograms and the 
% synch index
%------------------------------------------------------------------------
% parameters: 
% binsize: binsize used for binning spike times.
% windows: how many steps in each direction to calculate the correlogram
% ---------------------------------------------------------
% important outputs: 
% xcorr_sig_big: cell array of size # of rats, inside of each binary matrices of
% size # of cross-correlograms and length(windows) + 1. each 1 means a
% significant excitatory interaction
% xcorr_sig_small: cell array of size # of rats, inside of each binary matrices of
% size # of cross-correlograms and length(windows) + 1. each 1 means a
% significant inhibitory interaction
% synch_index: synchrony index (refere to paper)
% acorr_sig_big: like the xcorr_sig_big for autocorrelation
% acorr_sig_small: like the xcorr_sig_small for autocorrelation


clear all 
clc
close all


filepath = fileparts(pwd);
addpath(strcat(filepath,filesep,'Data',filesep));
addpath(strcat(filepath,filesep,'Results',filesep));
addpath(strcat(filepath,filesep,'Aux_functions',filesep));
%% load the results
load('All_Coeff.mat');
load('Counted_Spikes.mat');
load('raw_data.mat');

Number_of_rats = length(AllSpikes);
rat_id = unique(recording_id,'stable');
binsize = 100; % binsize that was used in msec;

%% Determining the pairs correlograms 
windows = [-10 10]; % how many time steps to go 
merge = cell(1,Number_of_rats);
xcorr_sig_big = cell(1,Number_of_rats);
xcorr_sig_small = cell(1,Number_of_rats);
sig_count = 0;
total_number_of_xcorr = 0;
cross_corr_all = cell(1,Number_of_rats);
avg_rand_cross_all = cell(1,Number_of_rats);
pairwise_max_all = cell(1,Number_of_rats);
pairwise_min_all = cell(1,Number_of_rats);
for i=1:Number_of_rats
    Number_of_ensembles = size(Activation_coeff_th{i},2);
    Number_of_xcorr = (Number_of_ensembles*(Number_of_ensembles-1))/2;
    merge{i} = zeros(Number_of_ensembles);
    total_number_of_xcorr = total_number_of_xcorr + ((Number_of_ensembles)*(Number_of_ensembles-1))/2;
    count = 1;
    for j=1:Number_of_ensembles
        time_coeff_1 = Activation_coeff_th{i}(:,j);
        idx_act1 = find(time_coeff_1 == 1);
        ts1 = (idx_act1 - 1);
        for k=j+1:Number_of_ensembles
            time_coeff_2 = Activation_coeff_th{i}(:,k);
            idx_act2 = find(time_coeff_2 == 1);
            ts2 = (idx_act2 - 1);
            [tsOffsets] = crosscorrelogram(ts1, ts2, windows);
            tsOffsets_rand = [];
            h = [];
            bintimes = (-10.5:1:10.5);
            for l=1:1000
                ts_rand1 = ts1;
                ts_rand2 = sort(ts2 + randi([-10 10],size(ts2)));
                tsOffsets_rand = crosscorrelogram(ts_rand1,ts_rand2,windows);
                [h(l,:),edges] = histcounts(tsOffsets_rand,bintimes);
            end
            avg_rand_cross = mean(h);
            pairwise_max = prctile(h,99);
            pairwise_min = prctile(h,1);
            x = movsum(edges,2)/2;
            x(1) = [];
            cross_counts = histcounts(tsOffsets,edges);
            cross_corr_all{i}(count,:) = cross_counts;
            avg_rand_cross_all{i}(count,:) = avg_rand_cross;
            pairwise_max_all{i}(count,:) = pairwise_max;
            pairwise_min_all{i}(count,:) = pairwise_min;
            if (cross_counts(11) > pairwise_max(11)); merge{i}(j,k) = 1; end
            xcorr_sig_big{i}(count,:) = cross_counts > pairwise_max;
            xcorr_sig_small{i}(count,:) = cross_counts < pairwise_min;
            count = count + 1;
            if sum(cross_counts > pairwise_max) > 0 || sum(cross_counts < pairwise_min) > 0
               sig_count = sig_count + 1;
            end
        end
    end
end

%% reporting results of significant interactions

time_inhibit = zeros(1,21);
time_excite = zeros(1,21);
for i=1:Number_of_rats
    tmp_small = double(xcorr_sig_small{i});
    time_inhibit = time_inhibit + sum(tmp_small);
    tmp_big = double(xcorr_sig_big{i});
    time_excite = time_excite + sum(tmp_big);
end

subplot(2,1,1);bar(time_inhibit,'BarWidth',1);xlabel('Time (S)','FontSize',12);
ylabel('Significant interactions','FontSize',12);xlim([0,22]);
set(gca,'xtick',[1,8,11,14,21],'xticklabels',{'-1','-0.3','0','0.3','1'},'Fontsize',12);
title('Inhibitory Interactions');axis square 
subplot(2,1,2);bar(time_excite,'BarWidth',1);xlabel('Time (S)','FontSize',12);
ylabel('Significant interactions','FontSize',12);xlim([0,22]);
set(gca,'xtick',[1,4,11,18,21],'xticklabels',{'-1','-0.7','0','0.7','1'},'Fontsize',12);
title('Excitatory Interactions'); axis square

%% measuring relative number of coincidents interactions at zero to total number of activations
number_0_xcorr = 0;
total_number_of_xcorr = 0;
count = 1;
synch_index = [];
for i=1:Number_of_rats
    tmp_xcorr = merge{i};
    number_0_xcorr = number_0_xcorr + length(find(tmp_xcorr == 1));
    total_number_of_xcorr = total_number_of_xcorr + (length(tmp_xcorr)*(length(tmp_xcorr)-1))/2;
    tmp_0 = find(tmp_xcorr == 1);
    for j=1:length(tmp_0)
        [r,c] = ind2sub(size(tmp_xcorr),tmp_0(j));
        for k=1:length(r)
            time_coeff_1 = Activation_coeff_th{i}(:,r(k));
            for l=1:length(c)
                time_coeff_2 = Activation_coeff_th{i}(:,c(l));
                synch_index(count) = (2*(sum(time_coeff_1.*time_coeff_2))/(sum(time_coeff_1) + sum(time_coeff_2)))*100;
                count = count + 1;
            end
        end
    end
end
figure()
histogram(synch_index);xlabel('synch index (%)');xline(mean(synch_index),'--r','mean','LineWidth',1);axis square
ylabel('counts')
 
% save_dir = strcat(pwd,'/results/');
% save_file = strcat(save_dir,'cross_interactions.mat');
% save(save_file,'avg_rand_cross_all'',cross_corr_all','pairwise_max_all','pairwise_min_all','sig_count','total_number_of_xcorr','xcorr_sig_big','xcorr_sig_small');
%% Measuring signiicant autocorrelograms
windows = [-10 10]; % how many time steps to go
acorr_sig_big = cell(1,Number_of_rats);
acorr_sig_small = cell(1,Number_of_rats);
total_number_of_acorr = 0;
signif_count = 0;
pairwise_max_all = cell(1,Number_of_rats);
pairwise_min_all = cell(1,Number_of_rats);
auto_corr_all = cell(1,Number_of_rats);
avg_rand_auto_all = cell(1,Number_of_rats);
for i=1:Number_of_rats
    Number_of_ensembles = size(Activation_coeff_th{i},2);
    Number_of_acorr = Number_of_ensembles;
    count = 1;
    for j=1:Number_of_ensembles
        time_coeff_1 = Activation_coeff_th{i}(:,j);
        idx_act1 = find(time_coeff_1 == 1);
        ts1 = (idx_act1 - 1);
        [tsOffsets] = crosscorrelogram(ts1, ts1, windows);
        tsOffsets_rand = [];
        bintimes = (-10.5:1:10.5);
        for l=1:1000
            ts_rand1 = sort(ts1 + randi([-10 10],size(ts1)));
            tsOffsets_rand = crosscorrelogram(ts_rand1,ts_rand1,windows);
            [h(l,:),edges] = histcounts(tsOffsets_rand,bintimes);
        end
        avg_rand_auto = mean(h);
        pairwise_max = prctile(h,99);
        pairwise_min = prctile(h,1);
        x = movsum(edges,2)/2;
        x(1) = [];
        se_rand_auto = std(h);
        cross_counts = histcounts(tsOffsets,edges);
        avg_rand_auto_all{i}(count,:) = avg_rand_auto;
        pairwise_max_all{i}(count,:) = pairwise_max;
        pairwise_min_all{i}(count,:) = pairwise_min;
        auto_corr_all{i}(count,:) = cross_counts;
        cross_counts(11) = avg_rand_auto(11); %not to counting the zero interactions
        total_number_of_acorr = total_number_of_acorr + 1;
        acorr_sig_big{i}(count,:) = cross_counts > pairwise_max;
        acorr_sig_small{i}(count,:) = cross_counts < pairwise_min;
        count = count + 1;
        if sum(cross_counts > pairwise_max) > 0 || sum(cross_counts < pairwise_min) > 0
            signif_count = signif_count + 1;
        end
    end
end

time_inhibit = zeros(1,21);
time_excite = zeros(1,21);
for i=1:length(acorr_sig_big)
    tmp_small = double(acorr_sig_small{i});
    time_inhibit = time_inhibit + sum(tmp_small);
    tmp_big = double(acorr_sig_big{i});
    time_excite = time_excite + sum(tmp_big);
end

figure()
subplot(2,1,1);bar(time_inhibit(11:end),'BarWidth',1);xlabel('Time (S)','FontSize',12);
ylabel('Significant interactions','FontSize',12);xlim([0,11]);
set(gca,'xtick',[1,6,11],'xticklabels',{'0','0.5','1'},'Fontsize',12);
title('Inhibitory Interactions');axis square 
subplot(2,1,2);bar(time_excite(11:end),'BarWidth',1);xlabel('Time (S)','FontSize',12);
ylabel('Significant interactions','FontSize',12);xlim([0,11]);
set(gca,'xtick',[1,6,11],'xticklabels',{'0','0.5','1'},'Fontsize',12);
title('Excitatory Interactions'); axis square 

% save_dir = strcat(pwd,'/results/');
% save_file = strcat(save_dir,'auto_interactions.mat');
% save(save_file,'avg_rand_auto_all'',auto_corr_all','pairwise_max_all','pairwise_min_all','sig_count','total_number_of_acorr','acorr_sig_big','acorr_sig_small');