%This is to plot summary data pooled across all fish on the properties of
%firing rate, and also EOD rate and firing rate relationship.
clearvars
close all
trialcnt=zeros(1,5);
trialduration=zeros(1,5);
dataF1_S1=load('20160713_Stage1_output.mat');
dataF1_S2=load('20160713_Stage2_output.mat');
dataF1_S4=load('20160713_Stage4_output.mat');
dataF1_S8_2=load('20160713_Stage8_output2.mat');

%calculate the total number of units
uname=dataF1_S2.uname;
[a,b]=size(dataF1_S2.uname);
m=1;
for i=1:a
    for j=1:b
        if isempty(uname{i,j})
            uname{i,j}='NO';
        else
            unitname(1,m)=uname{i,j}(2);
            m=m+1;
        end
        
    end
end
unitcnt(1)=length(unique(unitname(1,:)))
trialcnt(1)=length(dataF1_S2.avefr);
temp=0;
for i=1:length(dataF1_S2.time)
    temp=temp+dataF1_S2.time{i}(end)/3600;
end
trialduration(1)=temp;


dataF2_S1=load('20160822_Stage1_output.mat');
dataF2_S2=load('20160822_Stage2_output.mat');
dataF2_S4=load('20160822_Stage4_output.mat');
dataF2_S8_2=load('20160822_Stage8_output2.mat');

uname=dataF2_S2.uname;
[a,b]=size(dataF2_S2.uname);
m=1;
for i=1:a
    for j=1:b
        if isempty(uname{i,j})
            uname{i,j}='NO';
        else
            unitname(2,m)=uname{i,j}(2);
            m=m+1;
        end
        
    end
end
unitcnt(2)=length(unique(unitname(2,:)))
trialcnt=trialcnt+length(dataF2_S2.avefr);
temp=0;
for i=1:length(dataF2_S2.time)
    temp=temp+dataF2_S2.time{i}(end)/3600;
end
trialduration(2)=temp;

dataF3_S1=load('20160825_Stage1_output.mat');
dataF3_S2=load('20160825_Stage2_output.mat');
dataF3_S4=load('20160825_Stage4_output.mat');
dataF3_S8_2=load('20160825_Stage8_output2.mat');

uname=dataF3_S2.uname;
[a,b]=size(dataF3_S2.uname);
m=1;
for i=1:a
    for j=1:b
        if isempty(uname{i,j})
            uname{i,j}='NO';
        else
            unitname(3,m)=uname{i,j}(2);
            m=m+1;
        end
        
    end
end
unitcnt(3)=length(unique(unitname(3,:)))
trialcnt=trialcnt+length(dataF3_S2.avefr);
temp=0;
for i=1:length(dataF3_S2.time)
    temp=temp+dataF3_S2.time{i}(end)/3600;
end
trialduration(3)=temp;

dataF4_S1=load('20160607_Stage1_output.mat');
dataF4_S2=load('20160607_Stage2_output.mat');
dataF4_S4=load('20160607_Stage4_output.mat');
dataF4_S8_2=load('20160607_Stage8_output2.mat');

uname=dataF4_S2.uname;
[a,b]=size(dataF4_S2.uname);
m=1;
for i=1:a
    for j=1:b
        if isempty(uname{i,j})
            uname{i,j}='NO';
        else
            unitname(4,m)=uname{i,j}(2);
            m=m+1;
        end
        
    end
end
unitcnt(4)=length(unique(unitname(4,:)))
trialcnt=trialcnt+length(dataF4_S2.avefr);
temp=0;
for i=1:length(dataF4_S2.time)
    if ~isempty(dataF4_S2.time{i})
        temp=temp+dataF4_S2.time{i}(end)/3600;
    end
end
trialduration(4)=temp;

dataF5_S1=load('20160212_Stage1_output.mat');
dataF5_S2=load('20160212_Stage2_output.mat');
dataF5_S4=load('20160212_Stage4_output.mat');
dataF5_S8_2=load('20160212_Stage8_output2.mat');

uname=dataF5_S2.uname;
[a,b]=size(dataF5_S2.uname);
m=1;
for i=1:a
    for j=1:b
        if isempty(uname{i,j})
            uname{i,j}='NO';
        else
            unitname(5,m)=uname{i,j}(2);
            m=m+1;
        end
        
    end
end
unitcnt(5)=length(unique(unitname(5,:)))
trialcnt=trialcnt+length(dataF5_S2.avefr);

temp=0;
for i=1:length(dataF5_S2.time)
    temp=temp+dataF5_S2.time{i}(end)/3600;
end
trialduration(5)=temp;


avefrall=[cell2mat(dataF1_S2.avefr),cell2mat(dataF2_S2.avefr),cell2mat(dataF3_S2.avefr),cell2mat(dataF4_S2.avefr),cell2mat(dataF5_S2.avefr)];
histogram(avefrall,20,'Normalization','probability')
xlabel('Average firing rate Hz');
ylabel('Probability');
title(strcat('Ave FR, ',num2str(length(avefrall)),' trials'))
%%saveas(gcf,strcat('Hist_AveFRAllunits','.pdf'))
%Information content histogram

Infcontall=[cell2mat(dataF1_S4.Infcont),cell2mat(dataF2_S4.Infcont),cell2mat(dataF3_S4.Infcont),cell2mat(dataF4_S4.Infcont),cell2mat(dataF5_S4.Infcont)];
mean(Infcontall)
std(Infcontall)
histogram(Infcontall,20,'Normalization','probability')
xlabel('Place information bits/spike');
ylabel('Probability');
title(strcat('Inf cont, ',num2str(length(Infcontall)),' trials'))
%saveas(gcf,strcat('Hist_InfcontAllunits','.pdf'))

%caculate the directionality index -1: backwards, +1 forward
F1_index=(1-dataF1_S8_2.dir_ratio)./(1+dataF1_S8_2.dir_ratio);
F2_index=(1-dataF2_S8_2.dir_ratio)./(1+dataF2_S8_2.dir_ratio);
F3_index=(1-dataF3_S8_2.dir_ratio)./(1+dataF3_S8_2.dir_ratio);
F4_index=(1-dataF4_S8_2.dir_ratio)./(1+dataF4_S8_2.dir_ratio);
F5_index=(1-dataF5_S8_2.dir_ratio)./(1+dataF5_S8_2.dir_ratio);

Dirindexall=[F1_index,F2_index,F3_index,F4_index,F5_index];
mean(Dirindexall)
std(Dirindexall)
signtest(Dirindexall)
dfittool(Dirindexall)
% %histogram(DirRatioall,30,'Normalization','probability')
% %xlabel('Dir ratio all');
% %ylabel('Probability');
% title(strcat('Dir ratio, ',num2str(length(Infcontall)),' trials'))
% %saveas(gcf,strcat('Hist_DirratioAllunits','.pdf'))

EOD_rate_all=[cell2mat(dataF1_S1.EOD_rate');cell2mat(dataF2_S1.EOD_rate');cell2mat(dataF3_S1.EOD_rate');cell2mat(dataF4_S1.EOD_rate');cell2mat(dataF5_S1.EOD_rate')];

%pool all the EOD rates at the time of the spikes
figure

EOD_rate_interp_all=[cell2mat(reshape(dataF1_S1.EOD_rate_interp,[1,numel(dataF1_S1.EOD_rate_interp)])');cell2mat(reshape(dataF2_S1.EOD_rate_interp,[1,numel(dataF2_S1.EOD_rate_interp)])');cell2mat(reshape(dataF3_S1.EOD_rate_interp,[1,numel(dataF3_S1.EOD_rate_interp)])');cell2mat(reshape(dataF4_S1.EOD_rate_interp,[1,numel(dataF4_S1.EOD_rate_interp)])');cell2mat(reshape(dataF5_S1.EOD_rate_interp,[1,numel(dataF5_S1.EOD_rate_interp)])')];
histogram(EOD_rate_all,'Normalization', 'Probability','BinWidth',0.2)
hold on
histogram(EOD_rate_interp_all,'Normalization', 'Probability','BinWidth',0.2);
%xlim([45,75])
xlabel('EOD rate (Hz)')
ylabel('Probability')
%saveas(gcf,strcat('Hist_EODRateAllunits','.pdf'))

dfittool(EOD_rate_interp_all)

%pool all the swim speeds
allspeeds=[cell2mat(dataF1_S1.abs_swim_speed),cell2mat(dataF2_S1.abs_swim_speed),cell2mat(dataF3_S1.abs_swim_speed),cell2mat(dataF4_S1.abs_swim_speed),cell2mat(dataF5_S1.abs_swim_speed)];
allspikespeeds=[abs(cell2mat(dataF1_S1.speed_spike_POP'));abs(cell2mat(dataF2_S1.speed_spike_POP'));abs(cell2mat(dataF3_S1.speed_spike_POP'));abs(cell2mat(dataF4_S1.speed_spike_POP'));abs(cell2mat(dataF5_S1.speed_spike_POP'))];
%remove outliers, these are due to tracking errors
allspeeds(allspeeds>50)=[];
%remove negative speeds, this is a result of smoothing of distance
%travelled curve, otherwise speed is calculated as unsigned amplitude of
%displacement
allspeeds(allspeeds<0)=[];


figure
histogram(allspeeds,'Normalization', 'Probability','BinWidth',1)
hold on
histogram(allspikespeeds,'Normalization', 'Probability','BinWidth',1);
xlabel('swim speed cm/s')
ylabel('Probability')
%xlim([-0.5,30])
%saveas(gcf,strcat('Hist_speedAllunits','.pdf'))



dfittool(log(allspeeds))
