% This script estimates the directional STTC values for a dataset and 
% computes the significant edges among them, using the null distribution 
% test

clear 
close all
% Load the dataset
load 'Biological_data\saad16_003'

% Set it to the  value  
dt = 2; 

% Find the number of neurons of the dataset
num_of_neurons = size(psm_avalanche,2);

% Size of pool for the null distribution
maxnull = 100;

cnt = 1;
% For all possible pairs of neurons we estimate the directiona STTC values
for neuronA=1:size(psm_avalanche,2);
    neuronA
    for neuronB=1:size(psm_avalanche,2);
        if neuronA ~= neuronB
            flag = 0;
            A=psm_avalanche(:,neuronA);
            B=psm_avalanche(:,neuronB);
            
            for i=1:maxnull
                I = randi(size(psm_avalanche,1),1);
                A_shift = circshift(A,I);
                [A_shift_values(i),~,~] = STTC_AB(A_shift,B,dt);
            end
            
            STTCcirc(cnt,1) = neuronA;
            STTCcirc(cnt,2) = neuronB;
            [STTCcirc(cnt,3),STTCcirc(cnt,4),STTCcirc(cnt,5)]=STTC_AB(A,B,dt);
            STTCcirc(cnt,6)=mean(A_shift_values);
            STTCcirc(cnt,7)=std(A_shift_values);
            if(STTCcirc(cnt,3) > (STTCcirc(cnt,6) + 3*STTCcirc(cnt,7)) )
                flag = 1;
            end
            STTCcirc(cnt,8)=flag;
            
            cnt = cnt + 1;
        end
    end
end
% save('directional_animal_1', 'STTCcirc')


% Estimate the ids and the percentage of significant edges
pos_edges = size(psm_avalanche,2) * (size(psm_avalanche,2)-1);
ids = find(STTCcirc(:,8) == 1 & (STTCcirc(:,4)) > 3 & (STTCcirc(:,5)) > 3 );
sign_ids = length(ids);

perc_sign = sign_ids / pos_edges

% Create a matrix that contains only the significant edges
STTC_sign = STTCcirc(ids,:);

% save('significant_directional_animal_1', 'STTC_sign')

% Set the ids for each type of neuron (interneuron, pyramidals & astrocytes)
inter = interneurons;
astr = astrocytes;
other = union(inter,astr);%redcellz;
all = 1:num_of_neurons;
pyr = setdiff(all, other);
final_indexes = setdiff(all, astr);

% Produce CDFs about the STTC values and the type of pairs

ind_x_pyr = [];
ind_y_pyr = [];
ind_x_inter = [];
ind_y_inter = [];


for i =1:length(pyr)
    ind_x_pyr = [ind_x_pyr; find(STTC_sign(:,1) == pyr(i))];
    ind_y_pyr = [ind_y_pyr; find(STTC_sign(:,2) == pyr(i))];
end

for i =1:length(inter)
    ind_x_inter = [ind_x_inter; find(STTC_sign(:,1) == inter(i))];
    ind_y_inter = [ind_y_inter; find(STTC_sign(:,2) == inter(i))];
end

only_pyr = intersect(ind_x_pyr, ind_y_pyr);
only_inter = intersect(ind_x_inter, ind_y_inter);
pyr_iner = setdiff((1:length(STTC_sign)), only_pyr);
pyr_iner = setdiff(pyr_iner, only_inter);

STTC_only_inter1 = STTC_sign(only_inter,3);
STTC_only_pyr1 = STTC_sign(only_pyr,3);
STTC_pyr_inter1 = STTC_sign(pyr_iner,3);


%CDF and statistics

figure;
[x, f] = ecdf(STTC_only_pyr1);
plot(f,x, 'k','LineWidth',2);

min(STTC_only_pyr1);
max(STTC_only_pyr1);
mean(STTC_only_pyr1);
median(STTC_only_pyr1);

hold on;
grid on;
[x1, f1] = ecdf(STTC_pyr_inter1);
plot(f1,x1, 'r', 'LineWidth',1.5);

min(STTC_pyr_inter1)
max(STTC_pyr_inter1)
mean(STTC_pyr_inter1)
median(STTC_pyr_inter1)

hold on;
% [x2, f2] = ecdf(STTC_only_inter1);
% plot(f2, x2, 'c', 'LineWidth',1.5);

min(STTC_only_inter1)
max(STTC_only_inter1)
mean(STTC_only_inter1)
median(STTC_only_inter1)

legend('Pyramidal pairs', 'Interneuron-pyramidal pairs');
set(gca,'FontSize',16);
title('BD')
ylabel('CDF')
xlabel('Directional STTC ')
