%This is a program to further characterize EOD spike relationship pooling
%all trials together
close all
clearvars
exp_date=input('please type in the experiment date : ','s');
load(strcat(exp_date,'_Stage1_output.mat'));
load(strcat(exp_date,'_Stage2_output.mat'));
load(strcat(exp_date,'_Stage4_output.mat'));

[a,b]=size(uname);
w1=-4;% w1 second delay before each spike and w2 s delay after
w2=4;
res2=0.04;
interval=w1:res2:w2;

for tn=1:length(filename)

    for i=1:a
        if isempty(uname{i,tn})
            uname{i,tn}='NO';
        end
    end
end

k=1;
for tn=1:length(filename)
    for i=1:length(nuraw_count{tn})
        unitnames(k,:)=uname{i,tn};
        unitcounts(k)=nuraw_count{tn}(i);
        k=k+1;
    end
end

goodcells=unique(unitnames(:,2));
for w=1:length(goodcells)% go through all the cells
    cellname=strcat('U',goodcells(w));
    index=strmatch(cellname,uname);%find out which row/column that particular cell's data is located 
    [i,j]=ind2sub(size(uname),index);%i is trial number, j is the cell position in the matrix
    %for each cell, go through all the trials
    z=1;
    for k=1:length(index)
        EOD_rate_smooth=smoothdata(EOD_rate{j(k)},'movmean',15,'omitnan');
       % promin_thresh=0.2*nanstd((EOD_rate_smooth-nanmean(EOD_rate_smooth))/max(EOD_rate_smooth));%define prominent EOD peak threshold as 20% of the SD of the EOD over the whole trial
        %this version is for the non-normalized EOD rate plot 
        promin_thresh=0.2*nanstd(EOD_rate_smooth);%define prominent EOD peak threshold as 20% of the SD of the EOD over the whole trial
         
        speed_smooth=smoothdata(abs_swim_speed{j(k)},'movmean',15,'omitnan');
        promin_thresh_speed=0.2*nanstd(speed_smooth);
        
         for n=1:nuraw_count{j(k)}(i(k)) %for each single spike calculate EOD rate, swim direction, swim speed and sampling density in a 4 s window
            if nutimes{i(k),j(k)}(n)>=w2 && nutimes{i(k),j(k)}(n)+w2<=EOD_crct{j(k)}(end)
               lb=nutimes{i(k),j(k)}(n)+w1;%go w1 s before and w2 s after each spike
               rb=nutimes{i(k),j(k)}(n)+w2;
               %resample EOD_rate in the period
                temp_rate(z,:)=interp1(EOD_crct{j(k)}(find(EOD_crct{j(k)}<=rb & EOD_crct{j(k)}>lb)),EOD_rate{j(k)}(find(EOD_crct{j(k)}<=rb &EOD_crct{j(k)}>lb)),lb:res2:rb);
                %temp_rate(z,:)=temp_rate(z,:)/max(temp_rate(z,:));
                temp_angle(z,:)=interp1(xvalues{j(k)}(find(xvalues{j(k)}<=rb & xvalues{j(k)}>lb)),Theta_crct{j(k)}(find(xvalues{j(k)}<=rb &xvalues{j(k)}>lb)),lb:res2:rb);
                temp_speed(z,:)=interp1(xvalues{j(k)}(find(xvalues{j(k)}<=rb & xvalues{j(k)}>lb)),abs_swim_speed{j(k)}(find(xvalues{j(k)}<=rb &xvalues{j(k)}>lb)),lb:res2:rb);
                temp_speed_nozero(z,:)=temp_speed(z,:);
                temp_speed_nozero(z,temp_speed_nozero(z,:)<=1)=NaN;%take only speeds larger than 3 cm/s
                temp_sampdens(z,:)=temp_rate(z,:)./temp_speed_nozero(z,:);%calculate sampling density as rate divided by speed
                temp_rate_norm(z,:)=(temp_rate(z,:)-nanmean(temp_rate(z,:)))/max(temp_rate(z,:));%calculate mean subtracted max normalized rate
                z=z+1;
            else %if the spike time was such that it fell outside the overall bounds (This could happen for the first spikes or last spikes) set that window values to NaN
               temp_rate(z,:)=NaN(size(interval));               
               temp_angle(z,:)=NaN(size(interval));
               temp_speed(z,:)=NaN(size(interval));
               temp_sampdens(z,:)=NaN(size(interval));
               temp_rate_norm(z,:)=NaN(size(interval));
               z=z+1;
            end
            indexnoNaN=find(~isnan(temp_rate(z-1,:)));
            %do the subsequent analysis only if the current window contains
            %non - NaN values
            if ~isempty(indexnoNaN)
                %smooth the time series
                temp_rate(z-1,:)=smoothdata(temp_rate(z-1,:),'movmean',15,'omitnan');
                temp_rate_norm(z-1,:)=smoothdata(temp_rate_norm(z-1,:),'movmean',15,'omitnan');
                temp_angle(z-1,:)=smoothdata(temp_angle(z-1,:),'movmean',15,'omitnan');
                temp_speed(z-1,:)=smoothdata(temp_speed(z-1,:),'movmean',15,'omitnan');
                temp_sampdens(z-1,:)=smoothdata(temp_sampdens(z-1,:),'movmean',15,'omitnan');
%                 figure(1)
%                 plot(temp_rate(z-1,:),'r')
%                 hold on

                %find prominent peaks in the rate and speed trace
                [peakValues,indexOfPeak]=findpeaks((temp_rate(z-1,indexnoNaN(1):indexnoNaN(end))),'MinPeakProminence',promin_thresh);
                %find the location of the largest peak
                [aa,bb]=max(peakValues);
                if ~isempty(bb)
                    EODr_peak_index{w}(z-1)=indexOfPeak(bb);
                    EODr_peak_amp{w}(z-1)=peakValues(bb);
                else
                    EODr_peak_index{w}(z-1)=NaN;
                end

                % do the same for speed
                [peakValues,indexOfPeak]=findpeaks((temp_speed(z-1,indexnoNaN(1):indexnoNaN(end))),'MinPeakProminence',promin_thresh_speed);
                %find the location of the largest peak
                [aa,bb]=max(peakValues);
                if ~isempty(bb)
                    speed_peak_index{w}(z-1)=indexOfPeak(bb);
                    speed_peak_amp{w}(z-1)=peakValues(bb);
                else
                    speed_peak_index{w}(z-1)=NaN;
                end

            end
        end
    end
    %repeat the same procedure with shuffled spike times
    %uncomment for shuffle
    z=1;
    for k=1:length(index)
        for rep=1:100
            rep
            rng('shuffle') 
            [nurawshuffcnt,nurawshufftmp]=hist(nutimes{i(k),j(k)},time2{j(k)});
            %duration of the recording -30s
            shift_time=30:10:time2{j(k)}(end)-30;
            rnd_shift_time=shift_time(randi(length(shift_time)));
            nbins_shift=rnd_shift_time*fs;
            shuff=circshift(nurawshuffcnt,nbins_shift);
            shuffspktimes=nurawshufftmp(shuff==1);
            for n=1:length(shuffspktimes)
                if shuffspktimes(n)>=w2 && shuffspktimes(n)+w2<=EOD_crct{j(k)}(end)
                   lb=shuffspktimes(n)+w1;%go w1 s before and w2 s after each spike
                   rb=shuffspktimes(n)+w2;
                   %resample EOD_rate in the period
                   temp_rate_shuff(z,:)=interp1(EOD_crct{j(k)}(find(EOD_crct{j(k)}<=rb & EOD_crct{j(k)}>lb)),EOD_rate{j(k)}(find(EOD_crct{j(k)}<=rb &EOD_crct{j(k)}>lb)),lb:res2:rb);
                   temp_angle_shuff(z,:)=interp1(xvalues{j(k)}(find(xvalues{j(k)}<=rb & xvalues{j(k)}>lb)),Theta_crct{j(k)}(find(xvalues{j(k)}<=rb &xvalues{j(k)}>lb)),lb:res2:rb);
                   temp_speed_shuff(z,:)=interp1(xvalues{j(k)}(find(xvalues{j(k)}<=rb & xvalues{j(k)}>lb)),abs_swim_speed{j(k)}(find(xvalues{j(k)}<=rb &xvalues{j(k)}>lb)),lb:res2:rb);
                   temp_speed_nozero_shuff(z,:)=temp_speed_shuff(z,:);
                   temp_speed_nozero_shuff(z,temp_speed_nozero_shuff(z,:)<=1)=NaN;
                   temp_sampdens_shuff(z,:)=temp_rate_shuff(z,:)./temp_speed_nozero_shuff(z,:);
                   temp_rate_norm_shuff(z,:)=(temp_rate_shuff(z,:)-nanmean(temp_rate_shuff(z,:)))/max(temp_rate_shuff(z,:));
                   z=z+1;
                 else
                   temp_rate_shuff(z,:)=NaN(size(interval));
                   temp_rate_norm_shuff(z,:)=NaN(size(interval));
                   temp_angle_shuff(z,:)=NaN(size(interval));
                   temp_speed_shuff(z,:)=NaN(size(interval));
                   temp_sampdens_shuff(z,:)=NaN(size(interval));
                   z=z+1;
                end
               temp_rate_shuff(z-1,:)=smoothdata(temp_rate_shuff(z-1,:),'movmean',15,'omitnan');
               temp_rate_norm_shuff(z-1,:)=smoothdata(temp_rate_norm_shuff(z-1,:),'movmean',15,'omitnan');
               temp_angle_shuff(z-1,:)=smoothdata(temp_angle_shuff(z-1,:),'movmean',15,'omitnan');
               temp_speed_shuff(z-1,:)=smoothdata(temp_speed_shuff(z-1,:),'movmean',15,'omitnan');
               temp_sampdens_shuff(z-1,:)=smoothdata(temp_sampdens_shuff(z-1,:),'movmean',15,'omitnan');
%                figure(1)
%                plot(202:1:402,temp_rate_shuff(z-1,:),'b')

            end
        end
    end
    spktrig_EODr{w}=temp_rate;
    spktrig_EODr_norm{w}=temp_rate_norm;
    spktrig_Theta{w}=temp_angle;
    spktrig_speed{w}=temp_speed;
    spktrig_sampdens{w}=temp_sampdens;
    
    clear temp_rate;
    clear temp_rate_norm;
    clear temp_angle;
    clear temp_speed;
    clear temp_sampdens;
    %calculate averages for spike triggerred EOD rate, speed, angle and
    %sampling density
    spk_trig_ave_EOD_rate{w}=nanmean(spktrig_EODr{w});
    spk_trig_std_EOD_rate{w}=nanstd(spktrig_EODr{w});
    
    spk_trig_ave_EOD_rate_norm{w}=nanmean(spktrig_EODr_norm{w});
    spk_trig_std_EOD_rate_norm{w}=nanstd(spktrig_EODr_norm{w});
    
    spk_trig_ave_Theta{w}=nanmean(spktrig_Theta{w});
    spk_trig_std_Theta{w}=nanstd(spktrig_Theta{w});

    spk_trig_ave_speed{w}=nanmean(spktrig_speed{w});
    spk_trig_std_speed{w}=nanstd(spktrig_speed{w});
    
   spk_trig_ave_sampdens{w}=nanmean(spktrig_sampdens{w});
   spk_trig_std_sampdens{w}=nanstd(spktrig_sampdens{w});

   
    spk_trig_median_EOD_rate{w}=nanmedian(spktrig_EODr{w});
    spk_trig_median_EOD_rate_norm{w}=nanmedian(spktrig_EODr_norm{w});
    spk_trig_median_Theta{w}=nanmedian(spktrig_Theta{w});
    spk_trig_median_speed{w}=nanmedian(spktrig_speed{w});
    spk_trig_median_sampdens{w}=nanmedian(spktrig_sampdens{w});
    
        %do the same for shuffle
    spktrig_EODr_shuff{w}=temp_rate_shuff;
    spktrig_EODr_norm_shuff{w}=temp_rate_norm_shuff;
    spktrig_angle_shuff{w}=temp_angle_shuff;
    spktrig_speed_shuff{w}=temp_speed_shuff;
    spktrig_sampdens_shuff{w}=temp_sampdens_shuff;

    %calculate averages
    spk_trig_ave_EOD_rate_shuff{w}=nanmean(spktrig_EODr_shuff{w});
    spk_trig_std_EOD_rate_shuff{w}=nanstd(spktrig_EODr_shuff{w});
    spk_trig_median_EOD_rate_shuff{w}=nanmedian(spktrig_EODr_shuff{w});

    spk_trig_ave_EOD_rate_norm_shuff{w}=nanmean(spktrig_EODr_norm_shuff{w});
    spk_trig_std_EOD_rate_norm_shuff{w}=nanstd(spktrig_EODr_norm_shuff{w});
    spk_trig_median_EOD_rate_norm_shuff{w}=nanmedian(spktrig_EODr_norm_shuff{w});

    spk_trig_ave_angle_shuff{w}=nanmean(spktrig_angle_shuff{w});
    spk_trig_std_angle_shuff{w}=nanstd(spktrig_angle_shuff{w});
    spk_trig_median_angle_shuff{w}=nanmedian(spktrig_angle_shuff{w});
  
    spk_trig_ave_speed_shuff{w}=nanmean(spktrig_speed_shuff{w});
    spk_trig_std_speed_shuff{w}=nanstd(spktrig_speed_shuff{w});
    spk_trig_median_speed_shuff{w}=nanmedian(spktrig_speed_shuff{w});
    
    spk_trig_ave_sampdens_shuff{w}=nanmean(spktrig_sampdens_shuff{w});
    spk_trig_std_sampdens_shuff{w}=nanstd(spktrig_sampdens_shuff{w});
    spk_trig_median_sampdens_shuff{w}=nanmedian(spktrig_sampdens_shuff{w});
    
    
    clear temp_rate_shuff;
    clear temp_angle_shuff;
    clear temp_speed_shuff;
    clear temp_sampdens_shuff;
    clear temp_rate_norm_shuff;
    
    
    
    
%     figure(1)
%     plot(spk_trig_ave_EOD_rate{w},'k')
%     plot(202:1:402,spk_trig_ave_EOD_rate_shuff{w},'k')
  
    %calculate the significance of the peak
    [peakValues,indexOfPeak]=findpeaks(spk_trig_ave_EOD_rate{w});
    %find the location of the largest peak
    [aa,bb]=max(peakValues);
    %find amplitude at that index value across trials
    [r,p]=ttest(spktrig_EODr{w}(:,indexOfPeak(bb)));
    P_EOD_peak(w)=p;
    %calculate the average EODr peak time
    EODr_peak_time_ave(w)=mean(interval(EODr_peak_index{w}(~isnan(EODr_peak_index{w})&EODr_peak_index{w}>0)));
    EODr_peak_time_std(w)=std(interval(EODr_peak_index{w}(~isnan(EODr_peak_index{w})&EODr_peak_index{w}>0)));
    
    %do the same for speed
    %calculate the significance of the peak
    [peakValues,indexOfPeak]=findpeaks(spk_trig_ave_speed{w});
    %find the location of the largest peak
    [aa,bb]=max(peakValues);
    %find amplitude at that index value across trials
    [r,p]=ttest(spktrig_speed{w}(:,indexOfPeak(bb)),mean(spk_trig_ave_speed{w}));
    P_speed_peak(w)=p;
    
    speed_peak_time_ave(w)=mean(interval(speed_peak_index{w}(~isnan(speed_peak_index{w})&speed_peak_index{w}>0)));
    speed_peak_time_std(w)=std(interval(speed_peak_index{w}(~isnan(speed_peak_index{w})&speed_peak_index{w}>0)));

  



    %plot means and se
    figure
    title(strcat('Unit ',goodcells(w),'spike triggered EOD rate'))
    hold on
    plot(interval(3:end-2),spk_trig_ave_EOD_rate{w}(3:end-2),'k-');     
    hold on
    plot(interval(3:end-2),spk_trig_ave_EOD_rate{w}(3:end-2)+ spk_trig_std_EOD_rate{w}(3:end-2)/sqrt(length(spktrig_EODr{w})),'k:');
    plot(interval(3:end-2),spk_trig_ave_EOD_rate{w}(3:end-2)- spk_trig_std_EOD_rate{w}(3:end-2)/sqrt(length(spktrig_EODr{w})),'k:');
    
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_shuff{w}(3:end-2),'r-');
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_shuff{w}(3:end-2)+ spk_trig_std_EOD_rate_shuff{w}(3:end-2)/sqrt(length(spktrig_EODr_shuff{w})),'r:');
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_shuff{w}(3:end-2)- spk_trig_std_EOD_rate_shuff{w}(3:end-2)/sqrt(length(spktrig_EODr_shuff{w})),'r:');
   
    %plot the median
    plot(interval(3:end-2),spk_trig_median_EOD_rate{w}(3:end-2),'m-');
    hold on
    plot(interval(3:end-2),spk_trig_median_EOD_rate_shuff{w}(3:end-2),'b-');
   
    %saveas(gcf,strcat(exp_date,'spktrigEODrate mean U ',goodcells(w),'.pdf'))


    figure
    title(strcat('Unit ',goodcells(w),'spike triggered EOD rate normalized'))
    hold on
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_norm{w}(3:end-2),'k-');     
    hold on
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_norm{w}(3:end-2)+ spk_trig_std_EOD_rate_norm{w}(3:end-2)/sqrt(length(spktrig_EODr_norm{w})),'k:');
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_norm{w}(3:end-2)- spk_trig_std_EOD_rate_norm{w}(3:end-2)/sqrt(length(spktrig_EODr_norm{w})),'k:');
    
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_norm_shuff{w}(3:end-2),'r-');
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_norm_shuff{w}(3:end-2)+ spk_trig_std_EOD_rate_norm_shuff{w}(3:end-2)/sqrt(length(spktrig_EODr_norm_shuff{w})),'r:');
    plot(interval(3:end-2),spk_trig_ave_EOD_rate_norm_shuff{w}(3:end-2)- spk_trig_std_EOD_rate_norm_shuff{w}(3:end-2)/sqrt(length(spktrig_EODr_norm_shuff{w})),'r:');
    
    plot(interval(3:end-2),spk_trig_median_EOD_rate_norm{w}(3:end-2),'m-');
    hold on
    plot(interval(3:end-2),spk_trig_median_EOD_rate_norm_shuff{w}(3:end-2),'b-');

    %saveas(gcf,strcat(exp_date,'spktrigEODrate norm mean U ',goodcells(w),'.pdf'))

    figure
    plot(interval(3:end-2),spk_trig_ave_speed{w}(3:end-2),'k-');     
    hold on
    plot(interval(3:end-2),spk_trig_ave_speed{w}(3:end-2)+ spk_trig_std_speed{w}(3:end-2)/sqrt(length(spktrig_speed{w})),'k:');
    plot(interval(3:end-2),spk_trig_ave_speed{w}(3:end-2)- spk_trig_std_speed{w}(3:end-2)/sqrt(length(spktrig_speed{w})),'k:');
    
    plot(interval(3:end-2),spk_trig_median_speed{w}(3:end-2),'m-');     

    plot(interval(3:end-2),spk_trig_ave_speed_shuff{w}(3:end-2),'r-');     
    plot(interval(3:end-2),spk_trig_ave_speed_shuff{w}(3:end-2)+ spk_trig_std_speed_shuff{w}(3:end-2)/sqrt(length(spktrig_speed_shuff{w})),'r:');
    plot(interval(3:end-2),spk_trig_ave_speed_shuff{w}(3:end-2)- spk_trig_std_speed_shuff{w}(3:end-2)/sqrt(length(spktrig_speed_shuff{w})),'r:');

    plot(interval(3:end-2),spk_trig_median_speed_shuff{w}(3:end-2),'b-'); 
    
    %saveas(gcf,strcat(exp_date,'spktrigspeed median U ',goodcells(w),'.pdf'))

    figure
    plot(interval(3:end-2),spk_trig_ave_sampdens{w}(3:end-2),'k-');     
    hold on
    plot(interval(3:end-2),spk_trig_ave_sampdens{w}(3:end-2)+ spk_trig_std_sampdens{w}(3:end-2)/sqrt(length(spktrig_sampdens{w})),'k:');
    plot(interval(3:end-2),spk_trig_ave_sampdens{w}(3:end-2)- spk_trig_std_sampdens{w}(3:end-2)/sqrt(length(spktrig_sampdens{w})),'k:');
    
    plot(interval(3:end-2),spk_trig_median_sampdens{w}(3:end-2),'b-');     

    plot(interval(3:end-2),spk_trig_ave_sampdens_shuff{w}(3:end-2),'r-');     
    plot(interval(3:end-2),spk_trig_ave_sampdens_shuff{w}(3:end-2)+ spk_trig_std_sampdens_shuff{w}(3:end-2)/sqrt(length(spktrig_sampdens_shuff{w})),'r:');
    plot(interval(3:end-2),spk_trig_ave_sampdens_shuff{w}(3:end-2)- spk_trig_std_sampdens_shuff{w}(3:end-2)/sqrt(length(spktrig_sampdens_shuff{w})),'r:');
    plot(interval(3:end-2),spk_trig_median_sampdens_shuff{w}(3:end-2),'b-');     

    %saveas(gcf,strcat(exp_date,'spktrigSD median U ',goodcells(w),'.pdf'))
    
end


%save(strcat(exp_date,'_Stage8_V6_output.mat'),'spktrig_EODr','spktrig_EODr_norm','spktrig_Theta','spktrig_speed','spktrig_sampdens','spktrig_EODr_shuff','spktrig_EODr_norm_shuff','spktrig_angle_shuff','spktrig_speed_shuff','spktrig_sampdens_shuff','EODr_peak_amp','EODr_peak_index','speed_peak_amp','speed_peak_index','P_EOD_peak','P_speed_peak','speed_peak_time_ave','speed_peak_time_std','EODr_peak_time_ave','EODr_peak_time_ave','spktrig_EODr','spktrig_Theta','spktrig_speed');