load('Fig1andSuppleFigs_Behaviour')

bhv = group_level_stat_input.bhv ;
 bhv(bhv(:,7)==0,7) = nan;
 %bhv(:,[7 13 14]) = (bhv(:,[7 13 14])/8); 
 includedtrial = find(bhv(:,12)==1 & bhv(:,6)~=0 & bhv(:,4)~=0 & bhv(:,5)==1  & (bhv(:,3)+1==bhv(:,4))); %   
 bhv = bhv(includedtrial,:);

   lme = fitlme(tbl,'logit(response) ~ 1 + val + aro + sim  + (1|subject) ' ,...
     'DummyVarCoding','full' )

 [B,Bnames] = randomEffects(lme);
[beta,betanames] = fixedEffects(lme);
% [feCI,reCI] = coefCI(lme);
[psi,mse,stats] = covarianceParameters(lme);
std_group = cellfun(@sqrt,psi);
figure, bar([1 2 3],beta([3 2 4]),'FaceColor',[0.7 0.7 0.7]) 
hold on, errorbar([1 2 3],beta([3 2 4]),std_group([3 2 4]),'LineStyle','none','Color','k','LineWidth',2)
hold on, scatter(ones(1,7)+0.25*rand(1,7), beta(3)+B(15:21),'filled') % aro
hold on, scatter(2*ones(1,7)+0.25*rand(1,7), beta(2)+B(8:14),'filled') % val
hold on, scatter(3*ones(1,7)+0.25*rand(1,7),( beta(4)+B(22:28)),'filled') % sim
,xlim([0.5 3.5]), ylim([-0.5 0.5])


figure, bar([1 2 3],beta([3 2 4]),'FaceColor',[0.7 0.7 0.7]) 
hold on, errorbar([1 2 3],beta([3 2 4]),std_group([3 2 4]),'LineStyle','none','Color','k','LineWidth',2)
hold on, scatter(ones(1,7)+0.25*rand(1,7), beta(3)+B(15:21),'filled') % aro
hold on, scatter(2*ones(1,7)+0.25*rand(1,7), beta(2)+B(8:14),'filled') % val
hold on, scatter(3*ones(1,7)+0.25*rand(1,7),( beta(4)+B(22:28)),'filled') % sim
,xlim([0.5 3.5]), ylim([-0.5 0.5])




figure
 for N=1:20
     new_v =(bhv(:,14))+0.5*rand(length(bhv),1) ; new_a=8*rand(length(bhv),1); new_s=8*rand(length(bhv),1);
tbl_new = table(new_v,new_a,new_s,bhv(:,6)-1,bhv(:,11),'VariableNames',{'val','aro','sim','response','subject'});
y_pred = predict(lme,tbl_new);
hold on, scatter3(new_v,new_a,new_s,20,y_pred,'filled')
xlabel('Val'),ylabel('Aro'),zlabel('Sim')
 end
caxis([0.3 1])
set(gca,'LineWidth',2)

% Low Aro vs High Aro 
 for N=1:20
     new_v =(bhv(:,14))+0.5*rand(length(bhv),1) ; 
     new_a_low = 4.28*rand(length(bhv),1); 
     new_a_high = (8-4.28)*rand(length(bhv),1) + 4.28; 
     new_s=8*rand(length(bhv),1);
     tbl_new_low = table(new_v,new_a_low,new_s,bhv(:,6)-1,bhv(:,11),'VariableNames',{'val','aro','sim','response','subject'});
     y_pred_low = predict(lme,tbl_new_low);
     tbl_new_high = table(new_v,new_a_high,new_s,bhv(:,6)-1,bhv(:,11),'VariableNames',{'val','aro','sim','response','subject'});
     y_pred_high = predict(lme,tbl_new_high);
     for sim = 1:8
         trial_sim = find(new_s>=(sim-1) & new_s<=(sim));
         median_discrimination_low(sim,N) = median(y_pred_low(trial_sim)) ;
         median_discrimination_high(sim,N) = median(y_pred_high(trial_sim)) ;
     end
 end
 figure,  plot((1:8)-0.5,median(median_discrimination_low,2),'b--','LineWidth',3)
hold on, plot((1:8)-0.5,median(median_discrimination_high,2),'r--','LineWidth',3), ylim([0.4 1.1])
caxis([0.3 1])
set(gca,'LineWidth',3)

for sim = 1:8
         trial_sim = find(bhv(:,7)>=(sim-1) & bhv(:,7)<=(sim) & bhv(:,13)<=4.28);
         if ~isempty(trial_sim)
         median_discrimination_low_true(sim) = sum(bhv(trial_sim,6)==2)/sum(bhv(trial_sim,6)==1 | bhv(trial_sim,6)==2) ;
         else
             median_discrimination_high_true(sim)=nan;
         end
         trial_sim = find(bhv(:,7)>=(sim-1) & bhv(:,7)<=(sim) & bhv(:,13)>=4.28);
         if ~isempty(trial_sim)
         median_discrimination_high_true(sim) = sum(bhv(trial_sim,6)==2)/sum(bhv(trial_sim,6)==1 | bhv(trial_sim,6)==2) ;
         else
             median_discrimination_high_true(sim)=nan;
         end
         trial_sim = find(bhv(:,7)>=(sim-1) & bhv(:,7)<=(sim));
         median_discrimination_all(sim) = sum(bhv(trial_sim,6)==2)/sum(bhv(trial_sim,6)==1 | bhv(trial_sim,6)==2) ;
     end
figure,  plot(-((1:8)-0.5),median_discrimination_low_true,'LineWidth',3),ylim([0.4 1.1])
hold on, plot(-((1:8)-0.5),median_discrimination_high_true,'LineWidth',3),ylim([0.4 1.1])
hold on, plot(-((1:8)-0.5),median_discrimination_all,'Color',[0.7 0.7 0.7],'LineStyle','--','LineWidth',3),ylim([0.4 1.1])
