clear all

mainfolder='';

addpath(genpath([mainfolder 'subfunctions/']))
addpath(genpath([mainfolder 'behavioural_modeling/']))
addpath(genpath('subfunctions/rm_anova2'))

%% Panel 1

clear all

%data for polarplots

[oriC, oriP, oriN]=orivectors(16,6.43);

theta1=oriC
theta=[0 0 pi/2 pi/2 pi/2 pi/2 pi pi pi pi pi*1.5 pi*1.5 pi*1.5 pi*1.5 0 0];
uno=ones(1,16)
thetaall{1}=[theta theta(1)];
uno=[uno uno(1)];

theta1=oriC
theta2=[0 0 0 0 pi/2 pi/2 pi/2 pi/2 pi pi pi pi pi*1.5 pi*1.5 pi*1.5 pi*1.5]-pi/4
theta3=[theta1;theta2]
theta=mean(theta3)
uno=ones(1,16)
thetaall{2}=[theta theta(1)];
uno=[uno uno(1)];

theta=oriC
uno=ones(1,16)
thetaall{3}=[theta theta(1)];
uno=[uno uno(1)];

theta1=oriC
theta2=[0 0 0 0 pi/2 pi/2 pi/2 pi/2 pi pi pi pi pi*1.5 pi*1.5 pi*1.5 pi*1.5]+pi/4
theta3=[theta1;theta2]
theta=mean(theta3)
uno=ones(1,16)
thetaall{4}=[theta theta(1)];

theta=[0 0 0 0 pi/2 pi/2 pi/2 pi/2 pi pi pi pi pi*1.5 pi*1.5 pi*1.5 pi*1.5]+pi/4
uno=ones(1,16)
thetaall{5}=[theta theta(1)];
uno=[uno uno(1)];

A=0.2
thetaall{4}=A.*thetaall{5}+(1-A)*thetaall{3};

A=1
thetaall{5}=A.*thetaall{5}+(1-A)*thetaall{3};

a{1}=creatcatcardeucl;
a{3}=creatcontmodel(1);
a{5}=creatcatcardeucl45;

A=0.5;

a{2}=A*a{1}+(1-A)*a{3}
a{4}=A*a{5}+(1-A)*a{3}

% Plotting

figure;

subplot(2,5,1);
polarmodelplot(thetaall{1},uno)

subplot(2,5,2);
polarmodelplot(thetaall{2},uno)

subplot(2,5,3);
polarmodelplot(thetaall{3},uno)

subplot(2,5,4);
polarmodelplot(thetaall{4},uno)

subplot(2,5,5);
polarmodelplot(thetaall{5},uno)

subplot(2,5,6);
moldelmatplot(a{5},theta1)

subplot(2,5,7);
moldelmatplot(a{4},theta1)

subplot(2,5,8);
moldelmatplot(a{3},theta1)

subplot(2,5,9);
moldelmatplot(a{2},theta1)

subplot(2,5,10);
moldelmatplot(a{1},theta1)



%% Panel 2

clear all

%Loading data

beh=load('data_for_plots/Figure_4/all_behaviour_n41_trialbytrial_C05.mat')
load('data_for_plots/Figure_4/avg_polarbias_n41_C1.mat')

%Loadin eyedata

load('data_for_plots/Figure_4/test1_meanrecent_100pth_n41_lastsec.mat')
eye.all{1}(:,1)=av_a1cued1(:)-av_a1cued2(:)
eye.all{1}(:,2)=av_a2cued1(:)-av_a2cued2(:)
load('data_for_plots/Figure_4/test2_meanrecent_100pth_n41_lastsec.mat')

eye.all{1}(:,3)=av_a1cued1(:)-av_a1cued2(:)
eye.all{1}(:,4)=av_a2cued1(:)-av_a2cued2(:)


fp=(~isnan(eye.all{1}(:,1))+~isnan(eye.all{1}(:,2))...
    +~isnan(eye.all{1}(:,3))+~isnan(eye.all{1}(:,4)))==4; %checking we have data from all participants

eye.all{1}=eye.all{1}(:,:)

figure('Position' ,[100 600 950 600]);

%behavioural plot

subplot(1,3,1)
colores=[0 0 0];

% barplotbias(beh.all{1},[-0.51,0.51],'behaviour','A (bias index)',1)
lineartrendana(beh.all,1,1,[-.4,0.7],'behaviour','B (bias index)',colores,-0.6)

% polarplot of behav bias

oi=avg_bCueor;
oiac=avg_aCueor;

yourBias=[];
yourAccu=[];

for ppp=1:size(oi,2)
    behB=oi(:,ppp)'; %BheaviouralBias
    
    s=avg_Paramsr(ppp,1);     % noise (in memory/decision-making)
    A=avg_Paramsr(ppp,2);      % Key Parameter (squircle) 0: all circle; 1: all square
    C=avg_Paramsr(ppp,3);
    makefig=0;
    [yourAccu(:,ppp) yourBias(:,ppp)]=squircleBehave2compB(s,A,C,behB,makefig);
    
end
counter=1;

subplot(1,3,2)

polarplot([0 pi],[1 1],'k-','LineWidth',1);hold on;
polarplot([pi/2 pi*1.5],[1 1],'k-','LineWidth',1);hold on;



Aave=round(mean(avg_Paramsr(:,2),1),3)
Astd=std(avg_Paramsr(:,2),0,1);

quickplotcompBerrors(bb',(yourBias),(oi),[''],1,Aave);
counter=counter+1,


set(gca,'GridAlpha',0.25);
set(gca,'FontSize',20);

%eye plot

subplot(1,3,3)

lineartrendana(eye.all,1,0,[-.1,.175],'gaze','\Delta (rho) repulsion - attraction',colores,-0.15)


%% Analysis

%T test againg 0
mmm=mean(beh.all{1}(:,:),2);
mean(beh.all{1}) %mean of B index in each condition

h = lillietest(mmm);
[h,p,ci,stats] = ttest(mmm,0)
[h,p,ci,stats] = ttest(beh.all{1},0)

[p,h,stats] = ranksum(mmm,0)


% Anova 

al=[]
al= cat(1,beh.all{1}(:,1),beh.all{1}(:,2),beh.all{1}(:,3),beh.all{1}(:,4))';

p=size(beh.all{1},1);

S=[1:p,1:p,1:p,1:p];
F1=[ones(p,1);ones(p,1)*2;ones(p,1);ones(p,1)*2]'
F2=[ones(p*2,1);ones(p*2,1)*2]'
FACTNAMES={'item order', 'test order'}

stats = rm_anova2(al,S,F1,F2,FACTNAMES)

totalSS=sum([stats{2,2},stats{3,2},stats{4,2},stats{5,2},stats{6,2},stats{7,2}]);
etaS=[stats{2,2}/totalSS,stats{3,2}/totalSS,stats{4,2}/totalSS]

%% Linear trend analysis beh

 b_beh=lineartrendana(beh.all,0);

%%

%T test againg 0
mmm=mean(eye.all{1}(:,:),2);

h = lillietest(mmm);
[h,p,ci,stats] = ttest(mmm,0)
[h,p,ci,stats] = ttest(eye.all{1},0)

[p,h,stats] = ranksum(mmm,0)
1
% Anova 

al=[]
al= cat(1,eye.all{1}(fp,1),eye.all{1}(fp,2),eye.all{1}(fp,3),eye.all{1}(fp,4))';

p=size(eye.all{1}(fp,1),1);

S=[1:p,1:p,1:p,1:p];
F1=[ones(p,1);ones(p,1)*2;ones(p,1);ones(p,1)*2]'
F2=[ones(p*2,1);ones(p*2,1)*2]'
FACTNAMES={'item order', 'test order'}

stats = rm_anova2(al,S,F1,F2,FACTNAMES)

totalSS=sum([stats{2,2},stats{3,2},stats{4,2},stats{5,2},stats{6,2},stats{7,2}]);
etaS=[stats{2,2}/totalSS,stats{3,2}/totalSS,stats{4,2}/totalSS]


%% Linear trend analysis eye
 b_eye=lineartrendana(eye.all,0,0,[-.075,.2]);

%% Correlation between slopes in beh and eyes

[rho pval]=corr(b_beh',b_eye','Type','Pearson')

%%

figure; % barplotbias(beh.all{1},[-0.51,0.51],'behaviour','A (bias index)',1)
lineartrendana(beh.all,1,1,[-.9,0.9],'behaviour','B (bias index)')

%% Subfunctions

function polarmodelplot(theta,uno)

    polarplot(theta,uno,'k-o','MarkerFaceColor', 'k', 'MarkerSize', 8);
    r=gca;
    r.FontSize=20;
    r.RTickLabel=[];
    r.RGrid='off';
    rlim([0 1.25])    
    c=creat_u_d_model;
    d=creat_l_r_model;
    aa=c+d;
    cont=creatcontmodel(1);
    A=0.5;
    a=A*aa+(1-A)*cont
    
end
        

function moldelmatplot(a,theta1)

    imagesc(a);hold on;
    xlabel('orientation (°)')
    ylabel('orientation (°)')

    
  xticks = 1:16;  %adjust as appropriate, positive integers only
  xlabels = round(rad2deg(theta1(1:4:16)));  %time labels
    set(gca, 'XTick', 1:4:16, 'XTickLabel', xlabels, ...
        'YTick', 1:4:16, 'YTickLabel', xlabels, 'YAxisLocation','left','FontSize',20);
    
end


function barplotbias(data,yl,tito,yla,labels,sublabels,colores,xpos)

ax=notBoxPlot(data,'style','sdline')

line([0,5], [0,0], 'Color', 'k','LineStyle',':','LineWidth',2);hold on;
line([2.5,2.5], [-10,10], 'Color', 'k','LineStyle','--','LineWidth',2);hold on;

for i=1:4
        
    ax(i).semPtch.FaceColor = [0.75 0.75 0.75];
    ax(i).semPtch.EdgeColor = [0.75 0.75 0.75];


    ax(i).semPtch.LineWidth = 3;
    ax(i).data.MarkerSize = 8;

    ax(i).mu.Color = [0 0 0];
    ax(i).sd.Color = [0 0 0];
   

end

xticks([1:4])
xticklabels(labels)

% xtickangle(45)
ylim(yl)
xlim([0.5 4.5])
ylabel(yla)
set(gca,'FontSize',20);
title(tito)

text(1.5,xpos,sublabels{1},'HorizontalAlignment','center','FontSize',20,'FontWeight','bold')
text(3.5,xpos,sublabels{2},'HorizontalAlignment','center','FontSize',20,'FontWeight','bold')



end

function ball=lineartrendana(all,plotyes,beh,ylimits,btit,ylab,colores,xpos)

distances(:,1)=all{1}(:,2) %stim 2 cued 1st
distances(:,2)=all{1}(:,1) %stim 1 cued 1st
distances(:,3)=all{1}(:,4)%stim 2 cued 2nd
distances(:,4)=all{1}(:,3) %stim 1 cued 2nd

%% Check linear trend

% ball=[];
% yfit = [];

for ppp = 1:size(distances,1)   
    
    tof=[];
    for c=1:size(distances,2);
    tof=[tof; distances(ppp,c)];    
    end
    
    [b,dev,stats] = glmfit(1:size(distances,2),tof,'normal');
    yfit(ppp,:) = polyval([b(2,1),b(1,1)],[1,2,3,4]);
    
    ball(ppp) = b(2);
end

[h,p,ci,stats] = ttest(ball,0)

if plotyes

if beh
    labels={'Stim 2 ','Stim 1','Stim 2','Stim 1'};
    sublabel={'test 1', 'test 2'};
    
    
   
else
        labels={'Stim 2 ','Stim 1','Stim 2','Stim 1'};

%     labels={'stim 2 - delay 1','stim 1 - delay 1','stim 2 - delay 2','stim 1 - delay 2'};
    sublabel={'delay 1', 'delay 2'};

end

plot(yfit','Color',[colores 0.25],'LineWidth',1);hold on;

plot(mean(yfit,1)','Color',[0 0 0],'LineWidth',3);hold on;


barplotbias(distances,ylimits,btit,ylab,labels,sublabel,colores,xpos)
    
end

end