clear all

load('behavdata/datperitemandtrial.mat')
addpath(genpath([''])); %Add here the subfolder of subfunction (root folder)

% This data contains participant responses (1 =CW and 0 = CCW)

aCue{1,1}=squeeze(cw_resp_cu_uncue{1}(:,1,:,:)); % 1st obj and 1st Test
aCue{1,2}=squeeze(cw_resp_cu_uncue{1}(:,2,:,:)); % 1st obj 2nd Test
aCue{2,1}=squeeze(cw_resp_cu_uncue{2}(:,1,:,:)); % 2st obj and 1st Test
aCue{2,2}=squeeze(cw_resp_cu_uncue{2}(:,2,:,:)); % 2st obj 2nd Test

% This data contains the correct responses per trial (1 =CW and 0 = CCW)

tCue{1,1}=squeeze(cw_test_cu_uncue{1}(:,1,:,:)); % 1st obj and 1st Test
tCue{1,2}=squeeze(cw_test_cu_uncue{1}(:,2,:,:)); % 1st obj 2nd Test
tCue{2,1}=squeeze(cw_test_cu_uncue{2}(:,1,:,:)); % 2st obj and 1st Test
tCue{2,2}=squeeze(cw_test_cu_uncue{2}(:,2,:,:)); % 2st obj 2nd Test

% groupave=0 % If we want to calculate A using the group average
% 
% if groupave 
%     for item=1:2
%         for test=1:2
%             aCue{item,test}=mean(aCue{item,test}(:,:),2)
%         end
%     end   
% end

%% order the data in standard presentation

ao=deg2rad(cw_resp_cu_uncue{3}(:));
ao2=abs(2*pi-ao); %move anticlockwise
ao3=(angdiff(-ao2,-(ones(16,1).*pi*1.5)))+pi;

[bb ind]=sort(ao3);

for item=1:2
    for test=1:2
        aCueo{item,test}=aCue{item,test}(ind,:,:);
        tCueo{item,test}=tCue{item,test}(ind,:,:);

    end
end


%% Obtaining the best parameters

BestParams=nan(size(aCueo{1,1},2),4,2,2);
CvetIn=-0.5:0.1:0.5;%vector grid search of C

parfor ppp =1:size(aCueo{1,1},2)
    
    for item = 1:2
        for test=1:2
            
            aCueI=squeeze(aCueo{item,test}(:,ppp,:));
            tCueI=squeeze(tCueo{item,test}(:,ppp,:));

                        
            BestParams(ppp,:,item,test)=gridsearch_jld_trial(aCueI,tCueI,0,0,CvetIn,0.01); % This fuction is the the same as original, but we can add some input (ssqmat only based on accu,granularity for the grid search)
            
        end
    end
    ppp
end

%% Changing format

for item = 1:2
    for test=1:2
        Params{item,test}= squeeze(BestParams(:,:,item,test));% info about test 1
        Params{item,test}= squeeze(BestParams(:,:,item,test));% info about test 1
    end
end

%% Plotting all parameters

outputlabels={'s (noise)',' A (squircle factor)', 'C (cardinal precision)',' min SSQ'};

figure;

for lp = 1:4
all{1}= (cat(2,Params{1,1}(:,lp),Params{2,1}(:,lp),Params{1,2}(:,lp),Params{2,2}(:,lp)))

subplot(2,2,lp)
barplotbias(all{1},[min(all{1}(:)-0.0001),max(all{1}(:))+0.0001],' ',outputlabels{lp},1)

end


for lp = 2
all{1}= (cat(2,Params{1,1}(:,lp),Params{2,1}(:,lp),Params{1,2}(:,lp),Params{2,2}(:,lp)))

end


%% C against zero

allc{1}= (cat(2,Params{1,1}(:,3),Params{2,1}(:,3),Params{1,2}(:,3),Params{2,2}(:,3)))

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

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

%% C - Linear trend

cc=allc{1}

for ppp = 1:size(cc,1)   
    
    tof=[];
    for c=1:size(cc,2);
    tof=[tof; cc(ppp,c)];    
    end
    
    [b,dev,stats] = glmfit(1:size(cc,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)

%% Anova 

addpath '\rm_anova2'

removingout=0

if removingout ==1 %if removing the outlier
[m I]=max(Params{1,2}(:,4))
ind = 1:length((Params{1,2}(:,4)))~=I;
else
ind = 1:length((Params{1,2}(:,4)))~=0;
end

a1= (cat(1,Params{1,1}(ind,2),Params{2,1}(ind,2),Params{1,2}(ind,2),Params{2,2}(ind,2))')


p=(sum(ind));

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(a1,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]

%% Testing A index vs. 0

all{1}= (cat(2,Params{1,1}(:,2),Params{2,1}(:,2),Params{1,2}(:,2),Params{2,2}(:,2)))

for rrr=1:4
[h,p,ci,stats] = ttest(all{1}(:,rrr),0)
ttesre(rrr,1)=stats.tstat;
ttesre(rrr,2)=p;
end
%% Plotting A index (bars)

all{1}= (cat(2,Params{1,1}(:,2),Params{2,1}(:,2),Params{1,2}(:,2),Params{2,2}(:,2)))

names{1}='all'

figure;

for t=1
    
    avg=nanmean(all{t},1);
    sem=std(all{t})/sqrt(size(all{t},1));
    
%     subplot(1,3,t)
    % bar(rot_dif_list,avg);hold on,
    if groupave
     
        plot(avg,'o')
        
    else
    er=errorbar([],avg,sem,sem,'o');
    er.Color = [0 0 0];
    er.LineStyle = 'none';
    end
%     ylim([0 0.34])
    xlim([0 5])
%     line([0,5],[0.5 0.5], 'Color', 'k','LineStyle','--');hold on;
    
    xticks([1:4])
    xticklabels({'ori 1 1st test','ori 2 1st test','ori 1 2nd test','ori 2 2nd test'})
    
    % xlabel('distance between 1stcued minus 2nd cued ori (degrees)')
    ylabel('cardinal repulsion index (bias index)')
      
    title(['bias index - ' names{t}])
    
end


%% Plotting A index (bars) -version 2

figure('Position' ,[100 600 950 600]);
ax=notBoxPlot(all{1},'style','sdline')

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).mu.Color = [0 0 0];
ax(i).sd.Color = [0 0 0];
end

 xticks([1:4])
    xticklabels({'object 1 - test 1','object 2 - test 1','object 1 - test 2','object 2- test 2'})
        ylim([-0.41 .615])
        
        ylabel('A (bias index)')
       set(gca,'FontSize',13);


%% and groupal polar plot per item and test

for item = 1:2

for test=1:2
    oi=bCueo{item,test};
    oiac=aCueo{item,test};
    
    yourBias=[];
    yourAccu=[];
    for ppp=1:size(oi,2)
        behB=oi(:,ppp)'; %BheaviouralBias
        
        s=Params{item,test}(ppp,1);     % noise (in memory/decision-making)
        A=Params{item,test}(ppp,2);      % Key Parameter (squircle) 0: all circle; 1: all square
        C=Params{item,test}(ppp,3);      % Reduce noise near cardinal by this factor (1: off)
        makefig=0;
        [yourAccu(:,ppp) yourBias(:,ppp)]=squircleBehave2compB(s,A,C,behB,makefig);
        
    end
    counter=1;
figure;
    subplot(1,2,counter)
    
    quickplotcompBerrors(bb',(yourBias),(oi),['bias - object - '  num2str(item) ' test- ' num2str(test)],1);
    counter=counter+1,

       set(gca,'FontSize',12);

    
    subplot(1,2,counter)
    
    quickplotcompBerrors(bb',(yourAccu),(oiac),['accuracy - item- ' num2str(item) ' test- ' num2str(test)],0);
    counter=counter+1,

       set(gca,'FontSize',12);
    
end
end

%%
bCueor=cat(3,bCueo{1,1},bCueo{1,2},bCueo{2,1},bCueo{2,2});
avg_bCueor=mean(bCueor,3);

aCueor=cat(3,aCueo{1,1},aCueo{1,2},aCueo{2,1},aCueo{2,2});
avg_aCueor=mean(aCueor,3);

Paramsr=cat(3,Params{1,1},Params{1,2},Params{2,1},Params{2,2});
avg_Paramsr=mean(Paramsr,3);

%% and general polar plot of bias (averaged across items and tests)

    oi=avg_bCueor;
    oiac=avg_aCueor;
    
    yourBias=[];
    yourAccu=[];
    for ppp=1:size(oi,2)
        behB=oi(:,ppp)'; %Bevah Bias
        
        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);      % Reduce noise near cardinal by this factor (1: off)
        makefig=0;
        [yourAccu(:,ppp) yourBias(:,ppp)]=squircleBehave2compB(s,A,C,behB,makefig);
        
    end
    counter=1;
figure;
    
    quickplotcompBerrors(bb',(yourBias),(oi),['bias'],1);
    counter=counter+1,

       set(gca,'FontSize',12);


    %%

%% subfucntions

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

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).mu.Color = [0 0 0];
    ax(i).sd.Color = [0 0 0];
end

xticks([1:4])
if beh
xticklabels({'stim 1/test 1','stim 2/test 1','stim 1/test 2','stim 2/test 2'})
else
xticklabels({'stim 1/delay 1','stim 2/delay 1','stim 1/delay 2','stim 2/delay 2'})
end
xtickangle(45)
ylim(yl)
xlim([0.5 4.5])
ylabel(yla)
set(gca,'FontSize',20);
title(tito)

end