clear all

mainfolder='';

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

%% Loading the data 
folder1=[mainfolder,'results/mean_recentered_th100/euclideanmodels/single_trial_RSA_first_check_object_1_Pearson/']
folder2=[mainfolder,'results/mean_recentered_th100/euclideanmodels/single_trial_RSA_first_check_object_2_Pearson/']

load([mainfolder 'good_partis.mat']) %list of participant with a good performance
partis=good_partis'

counter=1;
for ppp=1:length(partis)
      
  toload=[folder1, partis{ppp,:},'.mat']
  load(toload)
  RHO1{counter}(:,:)=(RSA_obj_rot.rho);
  
  toload=[folder2, partis{ppp,:},'.mat']
  load(toload)
  RHO2{counter}(:,:)=(RSA_obj_rot.rho);
  
  counter=counter+1;

end

time=RSA_obj_rot.time;
%% Select, cut and average time window of interest

p1=1250; %start time
p2=2000; %end time

for ppp=1:size(RHO2,2)

[rho1cut{ppp},timecut1]=selecttimewindow(RHO1{ppp},time,p1,p2);
[rho2cut{ppp},timecut2]=selecttimewindow(RHO2{ppp},time,p1,p2);

rho1cutavg{ppp}=nanmean(rho1cut{ppp},2);
rho2cutavg{ppp}=nanmean(rho2cut{ppp},2);

end

%% Exploring the data
colores=(colormap((parula(size(RHO2,2)))));
gooddecoding=[]

figure;
for ppp=1:size(RHO2,2)
    
    subplot(2,1,1)
    title(['single trial ori decoding per parti - stim 1'])
    a=plot(timecut1,nanmean(rho1cut{ppp},1),'Color',colores(ppp,:),'LineWidth',2);hold on
    subplot(2,1,2)
    title(['single trial ori decoding per parti - stim 2'])
    b=plot(timecut1,nanmean(rho2cut{ppp},1),'Color',colores(ppp,:),'LineWidth',2);hold on
        
    if nanmean(nanmean(rho2cut{ppp},1))>= 0.02
      
      gooddecoding=[gooddecoding ppp];
      
    end
end

xlabel('time (ms)')
ylabel('correlation with model')

%% Calculating  the simult

nbins=9; %number of bins
edges=[];
rho1aveperbin=[];
counter=1;

for ppp= 1:size(RHO2,2)%gooddecoding %
    
    [N,edges(counter,:)] = histcounts(rho2cutavg{ppp},nbins);    
    Y = discretize(rho2cutavg{ppp},edges(counter,:));  
      
    for nb = 1:nbins
        
        rho2cutavg{counter}(Y==nb);
        rho1aveperbin(counter,nb)=nanmean(rho1cutavg{ppp}(Y==nb));
        
    end
    counter=counter+1;
end

%% Analysis and plot

[h,p,ci,stats]=ttest(rho1aveperbin,0);
corrected_p=double(p*(length(p))<0.05); % Bonfe correction and p<0.05
foundp=find(corrected_p==1);

mean_da=mean(rho1aveperbin,1)
sem_da=nanstd(rho1aveperbin,0,1)/sqrt(sum(~isnan(rho1aveperbin(1,:))))
meanedges=mean(edges,1)
figure('Position',([200 800 380 360]))
errorbar(meanedges(2:end),mean_da,sem_da,'-k','LineWidth', 2.5);hold on;

line([1,-1], [0,0], 'Color', 'k','LineStyle',':','LineWidth', 3);hold on;

scatter(meanedges(foundp+1),corrected_p(foundp)*0.15,65,'k*');hold on;

xlim([-.7 0.8])
ylim([-.05 .2])

xlabel('Stimulus 2 (rho)')
ylabel('Stimulus 1 (rho)')

 set(gca,'FontSize',22);

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


%% Check linear trend

for ppp = 1:size(rho1aveperbin,1)   
    
    tof=[];
    for c=1:size(rho1aveperbin,2);
    tof=[tof; rho1aveperbin(ppp,c)];    
    end
    
    [b,dev,stats] = glmfit(1:size(rho1aveperbin,2),tof,'normal');
    
    ball(ppp) = b(2);
end

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