% add function folder and clear
mfilePath = mfilename('fullpath');
if contains(mfilePath,'LiveEditorEvaluationHelper')
    mfilePath = matlab.desktop.editor.getActiveFilename;
end
disp(mfilePath);

mainfold   = mfilePath;


% go to data folder
warning off
tbl = readtable('Datatable_EXP2.csv');

% check and remove outliers
tbl           = tbl_subset(tbl,'ptrialout','<.2','errorsigm','<45');
actual_n      = numel(unique(tbl.obs));
n             = actual_n;

%% serial dependence (stimulus)
% split the conditions
stim_stable       = tbl_subset(tbl,'condition',1);
stim_volati       = tbl_subset(tbl,'condition',0);

% permutations
datasets                = {{stim_volati.delta,stim_volati.errorc} , ...
                           {stim_stable.delta,stim_stable.errorc}};
par                     = [2 .05; -20 .01; 20 .05];

rng(333)
[results, difference]   = sdp_permutation_stat('data',datasets,'nperm',10000,'start_par',par);

%% non param
lim               = [1 50]; % bin for model free bias

bias_stable_sti   = sdp_model_free_bias(stim_stable,'lim',lim,'variables',{'delta','errorc'});
bias_volati_sti   = sdp_model_free_bias(stim_volati,'lim',lim,'variables',{'delta','errorc'});

[~,p,~,stat_stable]      = ttest(bias_volati_sti, bias_stable_sti);
d                 = computeCohen_d(bias_volati_sti, bias_stable_sti,'paired');

%% serial dependence (decision)
% split the conditions
stim_stable       = tbl_subset(tbl,'condition',1);
stim_volati       = tbl_subset(tbl,'condition',0);

% permutations
datasets                = {{stim_volati.delta_btw,stim_volati.errorc} , ...
                           {stim_stable.delta_btw,stim_stable.errorc}};
par                     = [2 .05; -20 .01; 20 .05];

rng(333)
[results, difference]   = sdp_permutation_stat('data',datasets,'nperm',10000,'start_par',par);

%% non param
lim               = [1 50]; % bin for model free bias

bias_stable_dec   = sdp_model_free_bias(stim_stable,'lim',lim,'variables',{'delta_btw','errorc'});
bias_volati_dec   = sdp_model_free_bias(stim_volati,'lim',lim,'variables',{'delta_btw','errorc'});

[h,p,ci,stat_stable]     = ttest(bias_volati_dec, bias_stable_dec);
d                 = computeCohen_d(bias_volati_dec, bias_stable_dec,'paired');

%% ========================================================================
%% Sequence analysis
%% ========================================================================
nstim                   = unique(tbl.nstim);

%% effect of the last stimulus, depending on the sequence length and condition
stable_seq              = nan(n,numel(nstim));
volatile_seq            = nan(n,numel(nstim));
lim                     = [1 50]; % bin for model free bias
for k = 1:numel(nstim) % this is done with sliding averages (+/-1)
    stim_volati.sel     = ismember(stim_volati.nstim,nstim(k)-1:nstim(k)+1);
    volatile_seq(:,k)   = sdp_model_free_bias(tbl_subset(stim_volati,'sel',1),'lim',lim);
    stim_stable.sel     = ismember(stim_stable.nstim,nstim(k)-1:nstim(k)+1);
    stable_seq(:,k)     = sdp_model_free_bias(tbl_subset(stim_stable,'sel',1),'lim',lim);    
end

% rearrange for linear model
bias                    = cat(1,volatile_seq(:),stable_seq(:));
slength                 = cat(1,reshape(repmat(nstim,1,n)',numel(nstim)*n,1),reshape(repmat(nstim,1,n)',numel(nstim)*n,1));
stimstat                = categorical(reshape(repmat([0 1],numel(nstim)*n,1),numel(nstim)*n*2,1));
obs                     = reshape(repmat(1:n,numel(nstim)*2,1)',numel(nstim)*2*n,1);
tbl_sequence            = table(obs,slength,bias,stimstat);

%% effect of the last decision, depending on the sequence length and condition
stable_seq              = nan(n,numel(nstim));
volatile_seq            = nan(n,numel(nstim));
lim                     = [1 50]; % bin for model free bias
for k = 1:numel(nstim) % this is done with sliding averages (+/-1)
    stim_volati.sel     = ismember(stim_volati.nstim,nstim(k)-1:nstim(k)+1);
    volatile_seq(:,k)   = sdp_model_free_bias(tbl_subset(stim_volati,'sel',1),'variables',{'delta_btw','errorc'},'lim',lim);
    stim_stable.sel     = ismember(stim_stable.nstim,nstim(k)-1:nstim(k)+1);
    stable_seq(:,k)     = sdp_model_free_bias(tbl_subset(stim_stable,'sel',1),'variables',{'delta_btw','errorc'},'lim',lim);
end

% rearrange for linear model
bias                    = cat(1,volatile_seq(:),stable_seq(:));
slength                 = cat(1,reshape(repmat(nstim,1,n)',numel(nstim)*n,1),reshape(repmat(nstim,1,n)',numel(nstim)*n,1));
stimstat                = categorical(reshape(repmat([0 1],numel(nstim)*n,1),numel(nstim)*n*2,1));
obs                     = reshape(repmat(1:n,numel(nstim)*2,1)',numel(nstim)*2*n,1);
tbl_sequence            = table(obs,slength,bias,stimstat);

% run a linear model predicting the bias with slength (sequence length) and
% stimstat (condition)
model                   = fitlm(tbl_sequence,'bias ~ slength*stimstat');
disp(model)
