%plot histogram of peak times
close all
clearvars
w1=-4;% w1 second delay before each spike and w2 s delay after
w2=4;
res2=0.04;
interval=w1:res2:w2;
filenames=dir('*_V6_output.mat');
warning('off','all')
%cells with less than 10 spikes
badcells{1}=[];
badcells{2}=['6';'7'];
badcells{3}=[];
badcells{4}=['8';'9'];
badcells{5}=[];
for i=1:length(filenames)
    temp_EOD_peak=[];
    temp_EOD_peak_shuff=[];
    temp_EOD_peak_amp=[];
    temp_EOD_peak_amp_shuff=[];
    
    temp_EOD_dip=[];
    temp_EOD_dip_shuff=[];
    temp_EOD_dip_amp=[];
    temp_EOD_dip_amp_shuff=[];

    temp_speed_peak=[];
    temp_speed_peak_shuff=[];
    temp_speed_peak_amp=[];
    temp_speed_peak_amp_shuff=[];
    
    
    temp_speed_dip=[];
    temp_speed_dip_shuff=[];
    temp_speed_dip_amp=[];
    temp_speed_dip_amp_shuff=[];
    
    temp_SD_dip=[];
    temp_SD_dip_shuff=[];
    temp_SD_dip_amp=[];
    temp_SD_dip_amp_shuff=[];
    
    temp_SD_peak=[];
    temp_SD_peak_shuff=[];
    temp_SD_peak_amp=[];
    temp_SD_peak_amp_shuff=[];

 
    load(filenames(i).name);
    %put all the spktrig eod r peak times for all units of this fish in
    %arrays
    filenames(i).name
    %save peak times for units with significant peak in EODr
    %     Sig_mean_time{i}=EOD_Shuff_P(find(EOD_Shuff_P<0.05));
    vgoodcells=ones(size(goodcells));
    [a,b]=ismember(badcells{i},goodcells);
    vgoodcells(b)=0;
    Sig_peak_time{i}=EOD_P(EOD_P(:,2)<0.05 & vgoodcells);
    Sig_speed_peak_time{i}=speed_P(speed_P(:,2)<0.05 & vgoodcells);
        %calculate the average dip time
    [a,b]=size((spktrig_sampdens))
    for k=1:b
        %calculate the average
        Ave_spktrigSD=nanmean(spktrig_sampdens{k});
        %find the smallest dip of the average
        TF=islocalmin(Ave_spktrigSD);
        all_TF=find(TF==1);
        [y,ind]=min(Ave_spktrigSD(TF));
        temp(k)=interval(all_TF(ind));
        %find the largest peak in the average
        [PKS,LOCS]=findpeaks(Ave_spktrigSD);
        peak_index=LOCS(PKS==max(PKS));
        temp2(k)=interval(peak_index);
    end
    SD_ave_diptime{i}=temp(speed_P(:,2)<0.05 & EOD_P(:,2)<0.05 & vgoodcells);
    SD_ave_peaktime{i}=temp2(speed_P(:,2)<0.05 & EOD_P(:,2)<0.05 & vgoodcells);
    
    for k=1:length(spktrig_EODr)
        bb=find((badcells{i})==goodcells(k))
        if isempty(bb)
            if EOD_P(k,2)<0.05
                [a,b]=size(spktrig_EODr{k});
                for j=1:a
                    [PKS,LOCS]=findpeaks(spktrig_EODr{k}(j,:));
                    peak_index=LOCS(PKS==max(PKS));
                    temp_EOD_peak=[temp_EOD_peak,interval(peak_index)];
                    temp_EOD_peak_amp=[temp_EOD_peak_amp,max(PKS)];
                    %do the same for dip
                    
                    TF=islocalmin(spktrig_EODr{k}(j,:));
                    all_TF=find(TF==1);
                    [y,ind]=min(spktrig_EODr{k}(j,TF));
                    temp_EOD_dip=[temp_EOD_dip,interval(all_TF(ind))];
                    temp_EOD_dip_amp=[temp_EOD_dip_amp,y];

                    %calculate 
                end
                %do the same for the shuffle
                [a,b]=size(spktrig_EODr_shuff{k});
                for j=1:a
                    [PKS,LOCS]=findpeaks(spktrig_EODr_shuff{k}(j,:));
                    peak_index=LOCS(PKS==max(PKS));
                    temp_EOD_peak_shuff=[temp_EOD_peak_shuff,interval(peak_index)];
                    temp_EOD_peak_amp_shuff=[temp_EOD_peak_amp_shuff,max(PKS)];
                    
                    TF=islocalmin(spktrig_EODr_shuff{k}(j,:));
                    all_TF=find(TF==1);
                    [y,ind]=min(spktrig_EODr_shuff{k}(j,TF));
                    temp_EOD_dip_shuff=[temp_EOD_dip_shuff,interval(all_TF(ind))];
                    temp_EOD_dip_amp_shuff=[temp_EOD_dip_amp_shuff,y];
       
                end
                
            end
            if speed_P(k,2)<0.05
                [a,b]=size(spktrig_EODr{k});
                for j=1:a
                    [PKS,LOCS]=findpeaks(spktrig_speed{k}(j,:));
                    peak_index=LOCS(PKS==max(PKS));
                    temp_speed_peak=[temp_speed_peak,interval(peak_index)];
                    temp_speed_peak_amp=[temp_speed_peak_amp,max(PKS)];
                    
                    %find the dips
                                        
                    TF=islocalmin(spktrig_speed{k}(j,:));
                    all_TF=find(TF==1);
                    [y,ind]=min(spktrig_speed{k}(j,TF));
                    temp_speed_dip=[temp_speed_dip,interval(all_TF(ind))];
                    temp_speed_dip_amp=[temp_speed_dip_amp,y];


                end
                %do the same for the shuffle
                [a,b]=size(spktrig_EODr_shuff{k});
                for j=1:a
                    [PKS,LOCS]=findpeaks(spktrig_speed_shuff{k}(j,:));
                    peak_index=LOCS(PKS==max(PKS));
                    temp_speed_peak_shuff=[temp_speed_peak_shuff,interval(peak_index)];
                    temp_speed_peak_amp_shuff=[temp_speed_peak_amp_shuff,max(PKS)];
                    
                    TF=islocalmin(spktrig_speed_shuff{k}(j,:));
                    all_TF=find(TF==1);
                    [y,ind]=min(spktrig_speed_shuff{k}(j,TF));
                    temp_speed_dip_shuff=[temp_speed_dip_shuff,interval(all_TF(ind))];
                    temp_speed_dip_amp_shuff=[temp_speed_dip_amp_shuff,y];

                end
                
            end
            %for dip do it only for those with significant EOD and speed peak
            if EOD_P(k,2)<0.05 && speed_P(k,2)<0.05
                [a,b]=size(spktrig_EODr{k});
                for j=1:a
                    
                    TF=islocalmin(spktrig_sampdens{k}(j,:));
                    all_TF=find(TF==1);
                    [y,ind]=min(spktrig_sampdens{k}(j,TF));
                    temp_SD_dip=[temp_SD_dip,interval(all_TF(ind))];
                    temp_SD_dip_amp=[temp_SD_dip_amp,y];
                    %calculate peak times as well
                    [PKS,LOCS]=findpeaks(spktrig_sampdens{k}(j,:));
                    peak_index=LOCS(PKS==max(PKS));
                    temp_SD_peak=[temp_SD_peak,interval(peak_index)];
                    temp_SD_peak_amp=[temp_SD_peak_amp,max(PKS)];

                end
                %do the same for the shuffle
                [a,b]=size(spktrig_EODr_shuff{k});
                for j=1:a
                    
                    TF=islocalmin(spktrig_sampdens_shuff{k}(j,:));
                    all_TF=find(TF==1);
                    [y,ind]=min(spktrig_sampdens_shuff{k}(j,TF));
                    temp_SD_dip_shuff=[temp_SD_dip_shuff,interval(all_TF(ind))];
                    temp_SD_dip_amp_shuff=[temp_SD_dip_amp_shuff,y];
                    
                    [PKS,LOCS]=findpeaks(spktrig_sampdens_shuff{k}(j,:));
                    peak_index=LOCS(PKS==max(PKS));
                    temp_SD_peak_shuff=[temp_SD_peak_shuff,interval(peak_index)];
                    temp_SD_peak_amp_shuff=[temp_SD_peak_amp_shuff,max(PKS)];

                end
                
            end
            SD_dip{i}=temp_SD_dip;
            SD_dip_amp{i}=temp_SD_dip_amp;
            SD_dip_shuff{i}=temp_SD_dip_shuff;
            SD_dip_amp_shuff{i}=temp_SD_dip_amp_shuff;
            
            SD_peak{i}=temp_SD_peak;
            SD_peak_amp{i}=temp_SD_peak_amp;
            SD_peak_shuff{i}=temp_SD_peak_shuff;
            SD_peak_amp_shuff{i}=temp_SD_peak_amp_shuff;

            
            EOD_peak{i}=temp_EOD_peak;
            EOD_peak_shuff{i}=temp_EOD_peak_shuff;
            EOD_peak_amp{i}=temp_EOD_peak_amp;
            EOD_peak_amp_shuff{i}=temp_EOD_peak_amp_shuff;
            
            EOD_dip{i}=temp_EOD_dip;
            EOD_dip_shuff{i}=temp_EOD_dip_shuff;
            EOD_dip_amp{i}=temp_EOD_dip_amp;
            EOD_dip_amp_shuff{i}=temp_EOD_dip_amp_shuff;

            speed_peak{i}=temp_speed_peak;
            speed_peak_shuff{i}=temp_speed_peak_shuff;
            speed_peak_amp{i}=temp_speed_peak_amp;
            speed_peak_amp_shuff{i}=temp_speed_peak_amp_shuff;
            
            speed_dip{i}=temp_speed_dip;
            speed_dip_shuff{i}=temp_speed_dip_shuff;
            speed_dip_amp{i}=temp_speed_dip_amp;
            speed_dip_amp_shuff{i}=temp_speed_dip_amp_shuff;

            
            i
            %         Sig_mean_time{i}=EOD_Shuff_P(k,1);
        end
        
        clear temp
    end

    
end
% for i=1:5
%
%
% dfittool((EOD_peak{i}))
% dfittool(EOD_peak_shuff{i})
%
% % dfittool(EOD_peak_amp{i})
% % dfittool(EOD_peak_amp_shuff{i})
% end
dfittool(cell2mat(EOD_peak))
dfittool(cell2mat(EOD_peak_shuff))

dfittool(cell2mat(speed_peak))
dfittool(cell2mat(speed_peak_shuff))

dfittool(cell2mat(SD_dip))
dfittool(cell2mat(SD_dip_shuff))
dfittool(cell2mat(SD_peak))
dfittool(cell2mat(SD_peak_shuff))

%test whether the distributions of timing are different from uniform
% for i=1:5
%     [h,p]=ttest(EOD_peak{i})
%     [h,p]=ttest(EOD_peak_shuff{i})
%
% end
%compare the timing of speed peak SD dip and EOD peak
X0=cell2mat(EOD_peak);
X=X0(~isnan(X0));
Y0=cell2mat(speed_peak);
Y=Y0(~isnan(Y0));
Z0=cell2mat(SD_dip);
Z=Z0(~isnan(Z0));

W0=cell2mat(SD_peak);
W=W0(~isnan(W0));

A=[X,W];
AG=[ones(size(X)),2*ones(size(W))];
[p,a,s]=kruskalwallis(A,AG);
multcompare(s)
save('spktrigstats.mat')