%find whether or not swim speed is significantly lower during backward
%swims compared to forward swims. Also plot the correlation between swim
%Speed and EODr
close all
clear all
%% Calculate the average firing rate of cells with more than 10 spikes in each fish
filenames1=dir('*_Stage1_output.mat');
speed_s=[];
Theta_s=[];
speed_all=[];
Theta_all=[];
for i=1:length(filenames1)
        load(filenames1(i).name)
        %find if speed is significantly lower for backward swims
        %for each fish pool speed and theta points together 
        speed_s=cell2mat(abs_swim_speed);
        Theta_s=(cell2mat(Theta_crct'));

        Theta_s(speed_s>50 | speed_s<0)=[];
        speed_s(speed_s>50 | speed_s<0)=[];%these are tracking errors or interpolation artifacts
        speed_all=[speed_all,speed_s];
        Theta_all=[Theta_all;Theta_s];
        %plot(Theta_s,speed_s)
        speeds=[speed_s(Theta_s<20),speed_s(Theta_s>=150)];
        groups=[ones(size(speed_s(Theta_s<20))), 2*ones(size(speed_s(Theta_s>=150)))];
        speed_dir_comp(i,1)=kruskalwallis(speeds, groups','off')%p value
        speed_dir_comp(i,2)=mean(speed_s(Theta_s>=150))-mean(speed_s(Theta_s<30))
        
        speed_forward(i,1)=mean(speed_s(Theta_s<20));
        speed_forward(i,2)=std(speed_s(Theta_s<20));

        speed_backward(i,1)=mean(speed_s(Theta_s>=150));
        speed_backward(i,2)=std(speed_s(Theta_s>=150));
        
        speed_forward_med(i)=median(speed_s(Theta_s<20))
        speed_backward_med(i)=mean(speed_s(Theta_s>=150));
        
        %calculate the correlation between EOD rate and swim speed during
        %up state EODr>50 Hz
        speed_interp_pool=[];
        EOD_rate_pool=[];
        speed_pool=[];
        EOD_rate_interp_pool=[];
        for k=1:length(abs_swim_speed)
            if ~isempty(xvalues{k})
                speed_interp_pool=[speed_interp_pool;interp1(xvalues{k},abs_swim_speed{k},EOD_crct{k})];
                EOD_rate_pool=[EOD_rate_pool; EOD_rate{k}];
                
                speed_pool=[speed_pool,abs_swim_speed{k}];
                EOD_rate_interp_pool=[EOD_rate_interp_pool, interp1(EOD_crct{k},EOD_rate{k},xvalues{k})];

            end
        end
        [r,p]=corrcoef(speed_interp_pool(speed_interp_pool>0 & speed_interp_pool<=50 & EOD_rate_pool>40), EOD_rate_pool(speed_interp_pool>0 & speed_interp_pool<=50 & EOD_rate_pool>40));

        figure(1)
        subplot(5,1,i)
        [f,g]=fit(speed_interp_pool(speed_interp_pool>0 & speed_interp_pool<=50& EOD_rate_pool>40), EOD_rate_pool(speed_interp_pool>0 & speed_interp_pool<=50& EOD_rate_pool>40),'poly1');
        EOD_speed_r_p(i,:)=[r(1,2),p(1,2),g.rsquare];

        plot(speed_interp_pool(speed_interp_pool>0 & speed_interp_pool<=50),EOD_rate_pool(speed_interp_pool>0 &speed_interp_pool<=50),'.')
        L(i)=length(speed_interp_pool(speed_interp_pool>0 & speed_interp_pool<=50))
        hold on

        %plot(f)

        plot(f,speed_interp_pool(speed_interp_pool>0 & speed_interp_pool<=50 & EOD_rate_pool>40), EOD_rate_pool(speed_interp_pool>0 & speed_interp_pool<=50 & EOD_rate_pool>40));
        legend 'off'
        axis square
        %ylim([45,100]);
       % xlim([-1,55]);
       figure(2)
       plot(Theta_s, speed_s,'*')
       hold on

end
mean(speed_dir_comp(:,2))
std(speed_dir_comp(:,2))