clear all

   mainfolder='';

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

%% Uncomment to load the RSA data and perform the cluster perm

% folder{1}=[mainfolder, 'results/mean_recentered_th100/euclideanmodels/single_trial_with_between_RSA_object_1_cued_cont_Pearson/']
% folder{2}=[mainfolder, 'results/mean_recentered_th100/euclideanmodels/single_trial_with_between_RSA_object_1_uncued_cont_Pearson/']
% folder{3}=[mainfolder, 'eyedata/results/mean_recentered_th100/euclideanmodels/single_trial_with_between_RSA_object_2_cued_cont_Pearson/']
% folder{4}=[mainfolder, 'eyedata/results/mean_recentered_th100/euclideanmodels/single_trial_with_between_RSA_object_2_uncued_cont_Pearson/']
% 
% % %2nd cue
% 
% folder2{1}=[mainfolder, 'results/mean_recentered_th100/euclideanmodels/2nd_cue/single_trial_with_between_RSA_object_1_cued_cont_Pearson/']
% folder2{2}=[mainfolder, 'results/mean_recentered_th100/euclideanmodels/2nd_cue/single_trial_with_between_RSA_object_1_uncued_cont_Pearson/']
% folder2{3}=[mainfolder, 'results/mean_recentered_th100/euclideanmodels/2nd_cue/single_trial_with_between_RSA_object_2_cued_cont_Pearson/']
% folder2{4}=[mainfolder, 'results/mean_recentered_th100/euclideanmodels/2nd_cue/single_trial_with_between_RSA_object_2_uncued_cont_Pearson/']
% 
% load([mainfolder, 'good_partis.mat']) %list of participant with a good performance
% % 
%  partis=good_partis'
% 
% %% Loading with the script
% [RHO_id_group_1, time_1]=loadwithbet(folder, partis)
% [RHO_id_group_2, time_2]=loadwithbet(folder2, partis)
% 
% %% Ordering variables
% 
% a1_1=[];
% a2_1=[];
% a3_1=[];
% 
% a1_1{1}=RHO_id_group_1{1};
% a1_1{2}=RHO_id_group_1{2};
% a1_1{3}=RHO_id_group_1{4};
% a1_1{4}=RHO_id_group_1{5};
% 
% a2_1{1}=RHO_id_group_1{7};
% a2_1{2}=RHO_id_group_1{8};
% a2_1{3}=RHO_id_group_1{10};
% a2_1{4}=RHO_id_group_1{11};
% 
% a3_1{1}=RHO_id_group_1{3};
% a3_1{2}=RHO_id_group_1{6};
% a3_1{3}=RHO_id_group_1{9};
% a3_1{4}=RHO_id_group_1{12};
% 
% a1_2=[];
% a2_2=[];
% a3_2=[];
% 
% a1_2{1}=RHO_id_group_2{1};
% a1_2{2}=RHO_id_group_2{2};
% a1_2{3}=RHO_id_group_2{4};
% a1_2{4}=RHO_id_group_2{5};
% 
% a2_2{1}=RHO_id_group_2{7};
% a2_2{2}=RHO_id_group_2{8};
% a2_2{3}=RHO_id_group_2{10};
% a2_2{4}=RHO_id_group_2{11};
% 
% a3_2{1}=RHO_id_group_2{3};
% a3_2{2}=RHO_id_group_2{6};
% a3_2{3}=RHO_id_group_2{9};
% a3_2{4}=RHO_id_group_2{12};
% 
% %% Cluster analysis of the difference within-between (delay 1)
% 
% %Setting times (from -100 stim onset until the end of the delay period)
% t1=find(time_1==-100)
% t2=find(time_1==5850)
% t1_2=find(time_1==900)
% [timesignclust cluspval] = cluster_dif(a1_1,a2_1,time_1,t1,t2,t1_2)
% 
% % mean([timesignclust(1,1) timesignclust(2,1)]) %average of "starting" points;
% 
% % Creating 2 variables to use it for the plots and respect the colours
% timesignclust_1 = nan(4,2);
% timesignclust_2 = nan(4,2);
% timesignclust_1([1 3],:)=timesignclust([1 3],:);
% timesignclust_2([1 3],:)=timesignclust([2 4],:);
% 
% %% Cluster analysis of the difference within-between (delay 2)
% 
% %Setting times (from -100 stim onset until the end of the delay period)
% t1=find(time_2==-100);
% t2=find(time_2==2850);
% t1_2=find(time_2==-100);
% [timesignclust2 cluspva2] = cluster_dif(a1_2,a2_2,time_2,t1,t2,t1_2)
% 
% % mean([timesignclust(1,1) timesignclust(2,1)]) %average of "starting" point;
% 
% % Creating 2 variables to use it for the plots and respect the colours
% timesignclust_1_2 = nan(4,2);
% timesignclust_2_2 = nan(4,2);
% timesignclust_1_2([1 3],:)=timesignclust2([1 3],:);
% timesignclust_2_2([1 3],:)=timesignclust2([2 4],:);

%% Optional -> Loading the data ready to be plotted

load('data_for_plots/Figure_3/data_n41_2022.mat');

%% And plot

%Preparing colors
colores1=(colormap((winter(4))));
colores2=(colormap((cool(4))));
colores=[colores2(4,:);colores2(3,:);colores1(4,:);colores1(3,:)];

colores3=(colormap((winter(8))));
colores4=(colormap((cool(8))));
coloresclus=[colores4(7,:);colores4(7,:);colores3(7,:);colores3(7,:)];

figure('Position',([200 800 750 700]))
subplot(2,1,1)

leng1='cued within'
leng2='cued between'
leng3='uncued within'
leng4='uncued between'
titu=''
xl=[-500 5850];
yl=[-0.05 0.17];
plot4lines_errorback(time_1,a1_1,colores,xl,yl,leng1,leng2,leng3,leng4,titu,1,0,1,timesignclust_1,coloresclus)
set(gca,'FontSize',22);

subplot(2,1,2)

titu=''
plot4lines_errorback(time_1,a2_1,colores,xl,yl,leng1,leng2,leng3,leng4,titu,1,0,1,timesignclust_2,coloresclus)
set(gca,'FontSize',22);
 
 
 %% Plotting the 2nd delay

figure('Position',([200 800 420 700]))

subplot(2,1,1)

leng1='cued 1st within'
leng2='cued 1st between'
leng3='cued 2nd within'
leng4='cued 2nd between'
titu=''
xl=[-300 2850];
yl=[-0.05 0.17];
plot4lines_2nd_errorback(time_2,a1_2,colores,xl,yl,leng1,leng2,leng3,leng4,titu,1,0,1,timesignclust_1_2,coloresclus)
set(gca,'FontSize',22);

subplot(2,1,2)

titu=''

plot4lines_2nd_errorback(time_2,a2_2,colores,xl,yl,leng1,leng2,leng3,leng4,titu,1,0,1,timesignclust_2_2,coloresclus)
set(gca,'FontSize',22);


%% STATS ANOVA
% Encoding

 toi1ori1 = [0 1000];%time of interest for ori 1
 toi1ori2 = [1000 2000];%time of interest for ori 1

% toi1 = [0 500];%time of interest for ori 1

t1_1=find(time_1==toi1ori1(1));
t2_1=find(time_1==toi1ori1(2));

t1_1ori2=find(time_1==toi1ori2(1));
t2_1ori2=find(time_1==toi1ori2(2));

%ori1 cued
av_a1cued1 =nanmean(a1_1{1}(t1_1:t2_1,:),1); %  with 1
av_a1cued2 =nanmean(a1_1{2}(t1_1:t2_1,:),1); %  bet 1

%ori2 cued
av_a2cued1 =nanmean(a2_1{1}(t1_1ori2:t2_1ori2,:),1); %  w 2
av_a2cued2 =nanmean(a2_1{2}(t1_1ori2:t2_1ori2,:),1); %  betwen 2

% addpath(genpath(['/rm_anova2']))

aa= cat(1,av_a1cued1,av_a1cued2,av_a2cued1,av_a2cued2)'
p=size(av_a1cued1,2);

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={'between-within', 'item order'}

stats = rm_anova2(aa,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]

%% STATS ANOVA
% 1st delay

  toi1 = [2350 5850];%time of interest for ori 1
% toi1 = [0 500];%time of interest for ori 1
%    toi1 = [3350 5850];%time of interest for delay 1 1s_after_cueonset


t1_1=find(time_1==toi1(1));
t2_1=find(time_1==toi1(2));

%ori1 cued
av_a1cued1 =nanmean(a1_1{1}(t1_1:t2_1,:),1); %  with 1
av_a1cued2 =nanmean(a1_1{2}(t1_1:t2_1,:),1); %  bet 1

%ori2 cued
av_a2cued1 =nanmean(a2_1{1}(t1_1:t2_1,:),1); %  w 2
av_a2cued2 =nanmean(a2_1{2}(t1_1:t2_1,:),1); %  betwen 2

addpath(genpath(['/rm_anova2']))

aa= cat(1,av_a1cued1,av_a1cued2,av_a2cued1,av_a2cued2)'
p=size(av_a1cued1,2);

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={'between-within', 'item order'}

stats = rm_anova2(aa,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]

%% T-test (delay1)

[h,p,ci,stats]= ttest(av_a2cued1,av_a2cued2) % stim 2 cued - withing vs between
[h,p,ci,stats]= ttest(av_a1cued1,av_a1cued2) % stim 2 cued - withing vs between


%% 2nd delay

% toi12 = [350 2850];%time of interest for ori 1
  toil2 = [1350 2850];%time of interest for delay 1 1s_after_cueonset


t1_2=find(time_2==toi12(1));
t2_2=find(time_2==toi12(2));

%ori1 cued
av_a1cued1 =nanmean(a1_2{3}(t1_2:t2_2,:),1); %  with 1
av_a1cued2 =nanmean(a1_2{4}(t1_2:t2_2,:),1); %  bet 1

%ori2 cued
av_a2cued1 =nanmean(a2_2{3}(t1_2:t2_2,:),1); %  w 2
av_a2cued2 =nanmean(a2_2{4}(t1_2:t2_2,:),1); %  betwen 2

addpath(genpath(['/rm_anova2']))

aa= cat(1,av_a1cued1,av_a1cued2,av_a2cued1,av_a2cued2)'
p=size(av_a1cued1,2);

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={'between-within', 'item order'}

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

%%
%% T-test (delay2)

[h,p,ci,stats]= ttest(av_a2cued1,av_a2cued2) % stim 2 cued - withing vs between
[h,p,ci,stats]= ttest(av_a1cued1,av_a1cued2) % stim 1 cued - withing vs between

%% On the DIF

toi11 = [2350 5850];%time of interest for delay 1 1s_after_cueonset
    
t1_1=find(time_1==toi11(1));
t2_1=find(time_1==toi11(2));

toil2 = [350 2850];%time of interest for delay 1 1s_after_cueonset
  
t1_2=find(time_2==toil2(1));
t2_2=find(time_2==toil2(2));

av_a1cued1=nanmean(a1_1{1}(t1_1:t2_1,:),1)-nanmean(a1_1{2}(t1_1:t2_1,:),1); %dif item 1, cue 1 b-w 1st delay
av_a2cued1=nanmean(a2_1{1}(t1_1:t2_1,:),1)-nanmean(a2_1{2}(t1_1:t2_1,:),1); %dif item 2, cue 1 b-w 1st delay

av_a1cued2=nanmean(a1_2{3}(t1_2:t2_2,:),1)-nanmean(a1_2{4}(t1_2:t2_2,:),1); %dif item 1, cue 1 b-w 1st delay
av_a2cued2=nanmean(a2_2{3}(t1_2:t2_2,:),1)-nanmean(a2_2{4}(t1_2:t2_2,:),1); %dif item 2, cue 1 b-w 1st delay

aa= cat(1,av_a1cued1,av_a1cued2,av_a2cued1,av_a2cued2)'
p=size(av_a1cued1,2);

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={'delay 1st or 2nd', 'item order'}

stats = rm_anova2(aa,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]

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


%% eye plot

labels={'stim 1/delay 1','stim 2/delay 1','stim 1/delay 2','stim 2/delay 2'};

figure
barplotbias(eye.all{1},[-.075,.2],' ','diff within - between',labels)

%% Linear trend regarding distance

%from closest to longest disntace

distances(:,1)=av_a2cued1 %stim 2 cued 1st
distances(:,2)=av_a1cued1 %stim 1 cued 1st
distances(:,3)=av_a2cued2 %stim 2 cued 2nd
distances(:,4)=av_a1cued2 %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)
[p,h,stats] = signrank(ball,0)


    %%
    all{1}= distances
    
    labels={'Stim 2','Stim 1','Stim 2','Stim 1'};

figure

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

plot(mean(yfit,1)','Color',[0 0 0],'LineWidth',4);hold on;
barplotbias(all{1},[-.07,.19],' ','\Delta (rho) within - between',labels)


%% subfuctions

function barplotbias(data,yl,tito,yla,labels)

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 = 10;
    ax(i).mu.LineWidth = 5;

    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)
ylabel(yla)
xlim([0.5 4.5])

sublabels{1}='Delay 1'
sublabels{2}='Delay 2'
xpos=0.15

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

set(gca,'FontSize',25);

title(tito)

end

function [timesignclust cluspv]=allclustperm2(RHO_all,time,az,pv_cluster)
timesignclust=nan(size(RHO_all,2),2); % output variable (rows -> lines; columns -> start and end (ms))
for r=1:size(RHO_all,2);      
    [clusters, p_values, t_sums, permutation_distribution] =permutest(RHO_all{r},az{r},1,0.05,20000)
    if p_values(1)<pv_cluster
        
        timesignclust(r,1)=time{r}(clusters{1}(1))
        timesignclust(r,2)=time{r}(clusters{1}(end))
        
    end
    
    cluspv(r)=p_values(1);
end
end

function [timesignclust cluspv] = cluster_dif(a1_1,a2_1,time_1,t1,t2,t1_2)

%ori1 cued
wi1 =a1_1{1}(t1:t2,:); %  with 1
bet1 =a1_1{2}(t1:t2,:); %  bet 1
%ori1 uncued
wi1u =a1_1{3}(t1:t2,:); %  with 1
bet1u =a1_1{4}(t1:t2,:); %  bet 1
%ori2 cued
wi2 =a2_1{1}(t1_2:t2,:); %  w 2
bet2 =a2_1{2}(t1_2:t2,:); %  betwen 2
%ori2 uncued
wi2u =a2_1{3}(t1_2:t2,:); %  w 2
bet2u =a2_1{4}(t1_2:t2,:); %  betwen 2

RHOD_wi{1}=wi1();
RHOD_wi{2}=wi2();
RHOD_wi{3}=wi1u();
RHOD_wi{4}=wi2u();

RHOD_be{1}=bet1();
RHOD_be{2}=bet2();
RHOD_be{3}=bet1u();
RHOD_be{4}=bet2u();

pv_cluster=0.0125; % p-vals for clust perm test
% az=zeros(size(RHOD_diff{1})); % matrix to be compare (in this case against 0);

timesall{1}=time_1(t1:t2);
timesall{2}=time_1(t1_2:t2);
timesall{3}=time_1(t1:t2);
timesall{4}=time_1(t1_2:t2);

[timesignclust cluspv]=allclustperm2(RHOD_wi,timesall,RHOD_be,pv_cluster) % output variable (rows -> lines; columns -> start and end (ms))

end