%This is the program used to calculate the location selectivity of units.
%By ploting the distribution of lr vs ud [preference indices and generates
%the quiver plot of figure 9
close all
clear all
filenames=dir('*12_output_v6.mat');
DirSTLM_pool=[];
dir_index_pool=[];
DirThetaLM_pool=[];
index_ud_pref_pool=[];
index_lr_pref_pool=[];
for i=1:length(filenames)
    filenames(i).name
    load(filenames(i).name)
    length(DirThetaLM(:,1))
    DirSTLM_pool=[DirSTLM_pool;DirThetaLM];
    index_ud_pref_pool=[index_ud_pref_pool,ud_pref]
    index_lr_pref_pool=[index_lr_pref_pool,lr_pref]
    dir_index_pool=[dir_index_pool,dir_index];
    
    figure(1)
    plot(ud_pref,lr_pref,'*')
    xlabel('ud')
    ylabel('lr')
    hold on
    figure(2)
    plot(dir_index,ud_pref,'*')
    hold on
    xlabel('dir')
    ylabel('ud')
%     figure(4)
%     plot(DirThetaLM(:,1),dir_index,'*')
%     hold on
    figure(3)
    plot(dir_index,lr_pref,'*')
    hold on
    xlabel('dir')
    ylabel('lr')
    [r,p]=corrcoef(dir_index,ud_pref)

    clear DirThetaLM
    clear dir_index
    clear ud_pref
    clear lr_pref
  length(index_ud_pref_pool)
end
signtest(dir_index_pool)
dfittool(dir_index_pool)

[p,a,stats]=kruskalwallis(DirSTLM_pool(:,2),DirSTLM_pool(:,1))
xlabel('Swim Direction')
ylabel('Peak direction stlm')
% ylim([-180,180])

% [p,a,stats]=kruskalwallis(DirSTLM_pool_flip,DirSTLM_pool(:,1))
% xlabel('Swim Direction')
% ylabel('Peak direction stlm')
% % ylim([-180,180])
% 
figure

index1=find(index_ud_pref_pool>0.03);
index2=find(index_ud_pref_pool<-0.03);
kruskalwallis([index_ud_pref_pool(index1),abs(index_ud_pref_pool(index2))],[ones(size(index1)),2*ones(size(index2))])
[p,h]=signtest(index_ud_pref_pool)
dfittool(index_ud_pref_pool)
mean(index_ud_pref_pool)

[p,h]=signtest(index_lr_pref_pool)
dfittool(index_lr_pref_pool)
[p,h]=signtest(index_lr_pref_pool,0.02)
mean(index_lr_pref_pool)


figure
plot(dir_index_pool,index_ud_pref_pool,'*')
hold on
mdlr = fitlm(dir_index_pool,index_ud_pref_pool,'RobustOpts','on');
slope=table2array(mdlr.Coefficients(2,1))
intercept=table2array(mdlr.Coefficients(1,1))
plot(dir_index_pool,dir_index_pool*slope+intercept)
figure
plot(index_lr_pref_pool,index_ud_pref_pool,'*')

figure
plot(dir_index_pool,index_lr_pref_pool,'*')


[r,p]=corrcoef(index_ud_pref_pool,index_lr_pref_pool)
xlabel('ud')
ylabel('lr')

[r,p]=corrcoef(index_ud_pref_pool([1:2,4:end]),dir_index_pool([1:2,4:end]))
[r,p]=corrcoef(index_lr_pref_pool,dir_index_pool)



%divide data into three categories based on the direction index
index1=find(index_ud_pref_pool<(-0.2))%group 1
index3=find(index_ud_pref_pool>(0.2))%group2
index2=find(index_ud_pref_pool<(0.2) & index_ud_pref_pool>(-0.2))%group3
group=[ones(size(index1)),2*ones(size(index2)),3*ones(size(index3))];

[p,a,stats]=kruskalwallis([dir_index_pool(index1),dir_index_pool(index2),dir_index_pool(index3)],group)



[p,a,stats]=kruskalwallis([index_lr_pref_pool(index1),index_lr_pref_pool(index2),index_lr_pref_pool(index3)],group)

[p,a,stats]=kruskalwallis([DirSTLM_pool(index1,2)',DirSTLM_pool(index2,2)',DirSTLM_pool(index3,2)'],group)



%%compare ud, lr bias
index1=find((dir_index_pool)<-0.1)%group 1
index2=find((dir_index_pool)>0.1)%group2
group=[ones(size(index1)),2*ones(size(index2))];
[p,a,stats]=kruskalwallis([index_ud_pref_pool(index1),index_ud_pref_pool(index2)],group)

[p,a,stats]=kruskalwallis([dir_index_pool(index1),dir_index_pool(index2)],group)

%%compare ud, lr bias
index1=find(DirSTLM_pool==180)%group 1
index2=find(DirSTLM_pool==0)%group2
group=[ones(size(index1')),2*ones(size(index2'))];
[p,a,stats]=kruskalwallis([index_ud_pref_pool(index1),index_ud_pref_pool(index2)],group)

[p,a,stats]=kruskalwallis([dir_index_pool(index1),dir_index_pool(index2)],group)

% 
figure
for i=1:length(index_ud_pref_pool)
    if i==17 % C4 (U3)
        quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','r')
    elseif i==3 %First of C3 (U2)
        quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','b')
    %elseif i==20 %U6 C4
        %quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','m')
    elseif i==13 %one to the last of A3 UA
        quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','c')
    elseif i==7 %A3 (UA)
        quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','y')
    elseif i==14 %last of A5 (UE)
        quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','g')
    else
        quiver(0,0,index_lr_pref_pool(i),index_ud_pref_pool(i),'Color','k')
    end
    hold on
    xlim([-1,1])
    ylim([-1,1])
    grid on
    
end
Vect_ave=mean([index_lr_pref_pool;index_ud_pref_pool],2)
Vect_std=std([index_lr_pref_pool;index_ud_pref_pool],[],2)

quiver(0,0,Vect_ave(1),Vect_ave(2),'y')
xlim([-0.4,0.4])
ylim([-0.4,0.4])
% quiver(0,0,0,mean(dir_index_pool),'m')


