
%% Analysis of Experiment 3 and 4
% to use, put this file in the common root directory of the downloaded
% folder, set the variables expname and comment in the correct participant
% selection (var ptcs), or set your own selections for participants,
% inclusions etc... 
% Maëlle Lerebourg (2022)

%%
addpath(genpath(fullfile(pwd, 'support_funcs')))


%'3', '4' or '34_comb' (combined analysis of all participants)
expname = '34_comb';

if contains(expname, 'comb')
    resdir = {fullfile(pwd, 'data', 'Exp3'),fullfile(pwd, 'data', 'Exp4')};
else
    resdir = {fullfile(pwd, 'data', sprintf('Exp%s', expname)), {}};
end



for ne = 1:numel(resdir)
    if exist(resdir{ne}, 'dir')
        ptcfiles{ne,:}         = listfiles({'csv'}, resdir{ne});
    else
        ptcfiles{ne,:} = {};
    end
end
exps = [ones(size(ptcfiles{1})), 2 * ones(size(ptcfiles{2}))];
ptcfiles =  horzcat(ptcfiles{:});


for ne = 1:numel(resdir)
    if exist(resdir{ne}, 'dir')
        load(fullfile(resdir{ne}, 'targetdistr_ims.mat'));
        ims_targdistr{ne,:}         = ims_targdistr_os;
    else
        ims_targdistr{ne,:} = {};
    end
end



%
% %exclusion criteria
min_acc = 55;
rt_min   = 150;
rt_thr   = 3;   %std deviation within participants
%
%plot figures
save_fig = 1;
col_indivs = 0;




%ptcs = [1:12,14:17, 19:37, 39:46, 48:69, 71:85]; % Exp3
%ptcs = [1:6, 8:11,14:19, 21, 23:30, 33:52, 55:68, 70:73, 75:78, 80:85, 87:90, 93:95] %Exp4
ptcs = [1:12,14:17, 19:37, 39:46, 48:69, 71:85, [1:6, 8:11,14:19, 21, 23:30, 33:52, 55:68, 70:73, 75:78, 80:85, 87:90, 93:95] + 85]; %Exp3 and 4
%ptcs = [8];


acc_tr_m = nan(numel(ptcs), 40, 8);
acc_tr_mm = nan(numel(ptcs), 40, 8);
rt_tr_m  = nan(numel(ptcs), 40, 8);
rt_tr_mm  = nan(numel(ptcs), 40, 8);

pf = figure;

sigtxt = ['*', '*', '*'];
sigthr = [0.05, 0.01, 0.001];

acc.lim = [0, 100];
acc.dlim = [-30, 30];

rt.lim  = [200, 1600];
rt.dlim = [-200, 200];

L.lim = [200, 1800];
L.dlim = [-200, 400];

for np = 1:numel(ptcs)
    
    data = readtable(fullfile(resdir{exps(ptcs(np))}, ptcfiles{ptcs(np)}), 'ReadVariableNames', true);
    %shifted, to represent the preceding trial
    
    %select only trials of main experiment (not practice)
    main = data.main == 1;
    data = data(main,:);
    
    data_s = data([end, 1:end-1], :);
    rep_dim = contains(data.im, 'ori') == contains(data_s.im, 'ori');
    
    for nb =  1:max(data.blockn)
        blocki = find(data.blockn == nb);
        rep_dim(blocki(1)) = 0; 
    end 
    
    
    resp = ~ cellfun( @isempty, data.key_resp_keys );
    crt  = data.key_resp_corr == 1;
    crt_s = data_s.key_resp_corr == 1; 
    
    %distractor contexts
    r_c =  ismember(data.blockdistr, {'curv'});  % rectilinearity context
    r_o =  ismember(data.blockdistr, {'ori'});   % orientation context
    
    match = data.tvalid == 1 ;
    mismatch = data.tvalid == 0;
    
    cb = data.blockn <=4;
    ci = data.blockn >4;
    
    %code trialtypes,
    rm  = r_c & match;    %unique rectilinear trials, rect-blocks, matching
    rmm = r_o & mismatch; %unique rectilinearity trials, ori-blocks, mismatching
    om  = r_o & match;    %unique orientation trials, ori-blocks, matching
    omm = r_c & mismatch; %unique orientation trials, rect-blocks, mismatching
    
    
    %compute some general measures
    resprate(np)    = (sum(resp)/sum(main)) * 100;
    acc.gen(np)     = (sum(crt & resp )/ sum(resp)) * 100;
    
    if acc.gen(np) < min_acc
        warning( 'Participant %g was < %u %% correct, accuracy: %2.1f%% ', ptcs(np), min_acc, acc.gen(np));
    end
    
    all_rt     = mean(data.key_resp_rt(crt & resp)) * 1000 ;
    all_rt_std = std(data.key_resp_rt(crt & resp)) * 1000;
    
    rt_min = max(150, all_rt - (rt_thr * all_rt_std));
    rt_max = all_rt  + (rt_thr * all_rt_std);
    
    val_rt = crt & resp  & (data.key_resp_rt * 1000 >= rt_min) & (data.key_resp_rt * 1000 <=  rt_max);
    rts    = data.key_resp_rt * 1000;
    
    rt.gen(np)     = mean(rts(val_rt));
    rt.gen_std(np,3) = std(rts(val_rt));
    rt.gen_std(np,1) = std(rts(val_rt & cb));
    rt.gen_std(np,2) = std(rts(val_rt & ci));
    
    rt.rej(np)     = (sum(crt & resp & ((data.key_resp_rt * 1000 <= rt_min) | (data.key_resp_rt * 1000 >= rt_max )))./ sum(crt & resp)) * 100;
    
    L.gen(np)      = rt.gen(np) +((1 - acc.gen(np)./100) * (rt.gen_std(np)./std(crt)));
    
    targdistr = ismember(data.im,ims_targdistr{exps(ptcs(np))});
    acc.targdistr(np) =   (sum(crt & resp & targdistr) / sum(resp & targdistr)) * 100;
    
    totnum = sum(resp & targdistr);
    cornum = sum(crt  & resp & targdistr);
    uplim  = binoinv(1-0.05, totnum, 0.5);
    lowlim = binoinv(0.05  , totnum, 0.5);
    
    if cornum < uplim
        warning( 'Participant %g was %2.1f%% correct with familiar shape distractor, which is not better than chance (%g/%g correct, limit:%g)', ptcs(np),  acc.targdistr(np), cornum, totnum, uplim);
    end
    
    for cnd = 1:2
        
        if cnd == 1
            idxc = cb;
        elseif cnd == 2
            idxc = ci;
        end
        
        acc.rm(np,cnd) = (sum(crt & resp & rm & idxc)/sum(rm & resp & idxc)) * 100;
        acc.rmm(np,cnd)= (sum(crt & resp & rmm & idxc)/sum(rmm & resp & idxc)) * 100;
        acc.om(np,cnd) = (sum(crt & resp & om & idxc)/sum(om & resp & idxc)) * 100;
        acc.omm(np,cnd)= (sum(crt & resp & omm & idxc)/sum(omm & resp & idxc)) * 100;
        acc.dr(np,cnd) = (sum(crt & crt_s & resp & rep_dim & idxc)/sum(crt_s & resp & rep_dim & idxc)) * 100;
        acc.ds(np,cnd) = (sum(crt & crt_s & resp & ~rep_dim & idxc)/sum(crt_s & resp & ~rep_dim & idxc)) * 100; 
        acc.r(np, cnd)  = (sum(crt & resp & (rm|rmm) & idxc)/ sum((resp & (rm|rmm) & idxc))) * 100; 
        acc.o(np, cnd)  = (sum(crt & resp & (om|omm) & idxc)/ sum((resp & (om|omm) & idxc))) * 100;
        
        rt.rm(np,cnd)     = mean(rts( rm & val_rt & idxc));
        rt.rmm(np,cnd)    = mean(rts(rmm & val_rt & idxc));
        rt.om(np,cnd)     = mean(rts( om & val_rt & idxc));
        rt.omm(np,cnd)    = mean(rts(omm & val_rt & idxc));
        rt.dr(np, cnd)    = mean(rts(rep_dim & val_rt & idxc));
        rt.ds(np, cnd)    = mean(rts(~rep_dim & val_rt & idxc));
        rt.r(np, cnd)     = mean(rts((rm| rmm) & val_rt & idxc));
        rt.o(np, cnd)     = mean(rts((om| omm) & val_rt & idxc));
        
        L.rm(np,cnd)      = rt.rm(np,cnd) +((1 - acc.rm(np,cnd)./100) * (rt.gen_std(np,cnd)./std(crt(idxc))));
        L.rmm(np,cnd)     = rt.rmm(np,cnd)+((1 - acc.rmm(np,cnd)./100) * (rt.gen_std(np,cnd)./std(crt(idxc))));
        L.om(np,cnd)      = rt.om(np,cnd) +((1 - acc.om(np,cnd)./100) * (rt.gen_std(np,cnd)./std(crt(idxc))));
        L.omm(np,cnd)     = rt.omm(np,cnd)+((1 - acc.omm(np,cnd)./100) * (rt.gen_std(np,cnd)./std(crt(idxc))));
        L.dr(np, cnd)     = rt.dr(np, cnd)+((1 - acc.dr(np,cnd)./100) * (rt.gen_std(np,cnd)./std(crt(idxc))));
        L.ds(np, cnd)     = rt.ds(np, cnd)+((1 - acc.ds(np,cnd)./100) * (rt.gen_std(np,cnd)./std(crt(idxc))));
        L.r(np, cnd)      = rt.r(np, cnd) +((1 - acc.r(np,cnd)./100) * (rt.gen_std(np, cnd)./std(crt(idxc))));
        L.o(np, cnd)      = rt.o(np, cnd) +((1 - acc.o(np,cnd)./100) * (rt.gen_std(np, cnd)./std(crt(idxc))));
        
    end
    
    %measures per block
    
    %context effects are measured for every 2 blocks combined (1-2, 3-4,...), to compare
    %performance on the same trialtype when they are matching or
    %mismatching blockcontext
    
    for nb =  1:max(data.blockn)/2
        
        blocks = [(2 * nb) - 1, 2 * nb];
        acc.rm_bl(np,nb)  = (sum(crt & resp & rm & ismember(data.blockn, blocks))/sum(rm & resp &  ismember(data.blockn, blocks))) * 100;
        acc.rmm_bl(np,nb) = (sum(crt & resp & rmm & ismember(data.blockn, blocks))/sum(rmm & resp & ismember(data.blockn, blocks))) * 100;
        acc.om_bl(np,nb)  = (sum(crt & resp & om & ismember(data.blockn, blocks))/sum(om & resp & ismember(data.blockn, blocks))) * 100;
        acc.omm_bl(np,nb) = (sum(crt & resp & omm & ismember(data.blockn, blocks))/sum(omm & resp & ismember(data.blockn, blocks))) * 100;
        
        rt.rm_bl(np,nb)     = mean(rts( rm & val_rt &  ismember(data.blockn, blocks)));
        rt.rmm_bl(np,nb)    = mean(rts(rmm & val_rt &  ismember(data.blockn, blocks)));
        rt.om_bl(np,nb)     = mean(rts( om & val_rt &  ismember(data.blockn, blocks)));
        rt.omm_bl(np,nb)    = mean(rts(omm & val_rt &  ismember(data.blockn, blocks)));
        
        if nb <=2
            idxc = cb;
        elseif nb >= 3
            idxc = ci;
        end
        
        L.rm_bl(np,nb)      = rt.rm_bl(np,nb) +((1 - acc.rm_bl(np,nb)./100) * (rt.gen_std(np,ceil(nb/2))./std(crt(idxc))));
        L.rmm_bl(np,nb)     = rt.rmm_bl(np,nb)+((1 - acc.rmm_bl(np,nb)./100) * (rt.gen_std(np,ceil(nb/2))./std(crt(idxc))));
        L.om_bl(np,nb)      = rt.om_bl(np,nb) +((1 - acc.om_bl(np,nb)./100) * (rt.gen_std(np,ceil(nb/2))./std(crt(idxc))));
        L.omm_bl(np,nb)     = rt.omm_bl(np,nb)+((1 - acc.omm_bl(np,nb)./100) * (rt.gen_std(np,ceil(nb/2))./std(crt(idxc))));
        
    end
    
  
    
    for nb = 1:max(data.blockn)
        acc.block(np, nb) = (sum(crt & resp & data.blockn == nb)/sum(resp & data.blockn == nb))* 100;
        acc_tr_m(np,resp(data.blockn == nb) & match(data.blockn == nb),nb) = crt(data.blockn == nb & match & resp);
        acc_tr_mm(np,resp(data.blockn == nb)& mismatch(data.blockn == nb),nb) = crt(data.blockn == nb & mismatch & resp);
        
        rt.block(np, nb) =  mean(rts(  val_rt &  data.blockn == nb));
        rt_tr_m(np, val_rt(data.blockn == nb) & match(data.blockn == nb), nb) = rts(val_rt & data.blockn == nb & match);
        rt_tr_mm(np, val_rt(data.blockn == nb) & mismatch(data.blockn == nb), nb) = rts(val_rt & data.blockn == nb & mismatch);
        %
    end
    
    
    
end


perf = struct('Accuracy', acc,'RT',  rt, 'LISAS', L);
perf.Accuracy.unit = '(%)';
perf.RT.unit = '(ms)';
perf.LISAS.unit  = ' ';

if col_indivs
    cols = distinguishable_colors(numel(ptcs), {'w', 'k' });
else
    cols = repmat([0.75, 0.75, 0.75], numel(ptcs),1);
end

res = table();
f = fieldnames(perf);
conditions = {'blocked', 'interleaved'};

for nms = 1:numel(f)
    ms = f{nms};
    
    
    fprintf('\n %s\n', ms)
    fprintf( 'overall: %.2f %s (sd %.2f)\n', mean(perf.(ms).gen), replace(perf.(ms).unit,{'(', ')'}, {''}) , std(perf.(ms).gen));
    
    figure(pf)
    subplot(1,numel(f),nms)
    pl = plotSpread(perf.(ms).gen' , 'CategoryIdx', (1:numel(ptcs))' , 'categoryColors',cols )';
    errorbar(1, mean(perf.(ms).gen), std(perf.(ms).gen)/sqrt(numel(ptcs)), 'ok', 'MarkerSize', 12, 'LineWidth', 2, 'CapSize', 20);
    xticks([]);
    pbaspect([1,2.5,1 ])
    ylabel(sprintf('%s %s', ms, perf.(ms).unit))
    set(gcf,'color','w');
    set(gca,'linewidth',0.75)
    set(gca, 'FontSize', 13);
    
    
    mbl(:,:,1) = perf.(ms).omm_bl - perf.(ms).om_bl;
    mbl(:,:,2) = perf.(ms).rmm_bl - perf.(ms).rm_bl;
    d_bl = nanmean(mbl,3);
    
    if strcmp(ms, 'Accuracy')
        d_bl = - d_bl;
    end
    
    perf.(ms).d_bl= d_bl;
    
    
    m = nanmean(d_bl);
    sem = nanstd(d_bl)./sqrt(size(d_bl,1));
    CI  = m + 1.96 .* [sem; -sem];
    
    figure,
    h = plot(m, 'k', 'LineWidth', 1 );
    hold on
    xplot = [1:size(m,2), fliplr(1:size(m,2))];
    hold on
    fill(xplot,[CI(1,:), fliplr(CI(2,:))] ,1, 'facecolor', 'k', 'edgecolor', 'none', 'facealpha', 0.1)
    
    
    plot( [0.75, 4.25] , [0 0], 'k.:', 'LineWidth', 1)
    xlim([0.75,4.25] )
    xticks([1: 4 ])
    xticklabels({'1-2','3-4', '5-6', '7-8'})
    xlabel('Blocks')
    set(gcf,'color','w');
    ylabel(sprintf('Context facilitation %s %s', ms, perf.(ms).unit))
    box off
    set(gca,'linewidth',0.75)
    set(gca, 'FontSize', 13);
    if save_fig
        set( gcf, 'Renderer', 'painters');
        figname = sprintf('d_%s_blocks_%s', ms, expname);
        print( gcf, fullfile( pwd, figname), '-dpng');
        pause(0.5);
        print( gcf, fullfile( pwd, figname), '-dpdf');
        pause(0.5);
    end
    
    
    for cnd = 1:numel(conditions)
        fprintf('\n%s\n', conditions{cnd});
        
        
        X = [perf.(ms).om(:,cnd)', perf.(ms).omm(:,cnd)', perf.(ms).rm(:,cnd)', perf.(ms).rmm(:,cnd)']';
        %first factor is match, second trialtype
        group = [ones(size(perf.(ms).om(:,cnd)')), zeros(size(perf.(ms).omm(:,cnd)))', ones(size(perf.(ms).rm(:,cnd)))', zeros(size(perf.(ms).rmm(:,cnd)))'; ones(size(perf.(ms).om(:,cnd)))', ones(size(perf.(ms).omm(:,cnd)))', zeros(size(perf.(ms).rmm(:,cnd)))', zeros(size(perf.(ms).rm(:,cnd)))']';
        stats = mes2way(X,group, 'partialeta2', 'isDep',  [1 1 ] );
        d = mean([perf.(ms).omm(:,cnd)' - perf.(ms).om(:,cnd)'; perf.(ms).rmm(:,cnd)'- perf.(ms).rm(:,cnd)'],1);
        d_r =  mean([ perf.(ms).rmm(:,cnd)'- perf.(ms).rm(:,cnd)'],1);
        d_o =  mean([perf.(ms).omm(:,cnd)' - perf.(ms).om(:,cnd)'],1);
        d_int = d_r - d_o;
        d_tr = mean([perf.(ms).omm(:,cnd)'; perf.(ms).om(:,cnd)'],1) - mean([perf.(ms).rmm(:,cnd); perf.(ms).rm(:,cnd)],1);
    
    
    
        
        if strcmp(ms, 'Accuracy')
            d = -d;
            d_r = -d_r;
            d_o = -d_o;
        end
        
        perf.(ms).d(:,cnd) = d;
        perf.(ms).d_r(:,cnd) = d_r;
        perf.(ms).d_o(:,cnd) = d_o;
        
        [hd, pd, CId, st_d] = ttest(d);
        [hr, pr, CIr, str] = ttest(d_r);
        [ho, po, CIo, sto] = ttest(d_o);

        
        fprintf('overall context facilitation: %.2f %s (sd %.2f)\n', mean(d), replace(perf.(ms).unit,{'(', ')'}, {''}), std(d));
        fprintf('- rect trials: %.2f %s(sd %.2f), ([%.2f, %.2f], t(%u) = %.2f, p = %.3f)\n', mean(d_r), replace(perf.(ms).unit,{'(', ')'}, {''}), std(d_r), CIr, str.df, str.tstat, pr);
        fprintf('- ori  trials: %.2f %s(sd %.2f), ([%.2f, %.2f], t(%u) = %.2f, p = %.3f)\n', mean(d_o), replace(perf.(ms).unit,{'(', ')'}, {''}), std(d_o), CIo, sto.df, sto.tstat, po);
        
        
        res = [table(perf.(ms).om(:,cnd), perf.(ms).omm(:,cnd), perf.(ms).rm(:,cnd), perf.(ms).rmm(:,cnd), d', d_r', d_o', perf.(ms).dr(:,cnd) , perf.(ms).ds(:,cnd),  d_tr',d_int' ,'VariableNames',strcat( {ms}, '_', {'om', 'omm', 'rm', 'rmm', 'd', 'd_r', 'd_o', 'd_repdim', 'd_switchdim', 'd_tr', 'd_int'}, '_', conditions{cnd})) , res];
        
        %make figure
        figure
        hold on
        subplot(1,2,1)
        hold on
        
        rv = violinplot(perf.(ms).rm(:,cnd)  ,{1}, 'xpos', [0.75], 'ViolinColor', [0 0.4470 0.7410], 'Width', 0.2);
        oi = violinplot(perf.(ms).rmm(:,cnd),{2}, 'xpos', [1.25], 'ViolinColor',[0.6350 0.0780 0.1840], 'Width', 0.2);
        ri = violinplot(perf.(ms).om(:,cnd),{3}, 'xpos', [2.25], 'ViolinColor', [0 0.4470 0.7410] , 'Width', 0.2);
        ov = violinplot(perf.(ms).omm(:,cnd), {4}, 'xpos', [2.75], 'ViolinColor', [0.6350 0.0780 0.1840], 'Width', 0.2);
        xticks([1, 2.5])
        ylim(perf.(ms).lim)
        xlim([0.5, 3])
        set(gca,'linewidth',0.75)
        set(gca, 'FontSize', 13);
        xticklabels({'rectilinearity', 'orientation'})
        xlabel('trial type')
        ylabel(sprintf('%s %s', ms, perf.(ms).unit))
        title(conditions{cnd})
        
        
        p2 = subplot(1,2,2);
        hold on
        violinplot(d',  {1} , 'ViolinColor', [0, 0, 0]);
        plot( [0.5,1.5] , [0 0], 'k.:', 'LineWidth', 1);
        xticks([])
        ylim(perf.(ms).dlim)
        set(gca,'linewidth',0.75)
        set(gca, 'FontSize', 13);
        legend(p2, [rv.MedianPlot, ov.MedianPlot], {'matching context', 'mismatching context'}, 'Location', 'northoutside', 'FontSize', 14  )
        legend boxoff
        pbaspect([1,2.5,1])
        xlim([0.5, 1.5])
        ylabel('Context facilitation')
        
        lims = ylim;
        text(1,min([lims(2), max(d) + 0.7 * max(d)]), sigtxt(pd < sigthr), 'HorizontalAlignment', 'center', 'Fontsize', 20)
        
        
        if save_fig
            set( gcf, 'Renderer', 'painters');
            figname = sprintf('%s_trialsxblocks_c_%s_%s', ms, conditions{cnd},expname);
            print( gcf, fullfile( pwd, figname), '-dpng');
            pause(0.5);
            print( gcf, fullfile( pwd, figname), '-dpdf');
            pause(0.5);
        end
        
        

    end
    
        figure, 
        h = scatter( perf.(ms).d(:,1), perf.(ms).d(:,2)  , 'filled', 'jitter', 'on','MarkerFaceColor' , 'k');
        h.MarkerFaceAlpha = 0.5;
        %ylim([0,6])
        xlabel('Context facilitation blocked')
        ylabel('Context facilitation interleaved')
       % xticks([1:5])
        set(gca,'linewidth',0.75)
        set(gca, 'FontSize', 13)
        set( gcf, 'Renderer', 'painters');
        
    
    
    d1 = mean([ nanmean(perf.(ms).rmm_bl(:,[1,3]),2) - nanmean(perf.(ms).rm_bl(:,[1,3]),2) , nanmean(perf.(ms).omm_bl(:,[1,3]),2) - nanmean(perf.(ms).om_bl(:,[1,3]),2) ],2);
    d2 = mean([ nanmean(perf.(ms).rmm_bl(:,[2,4]),2) - nanmean(perf.(ms).rm_bl(:,[2,4]),2) , nanmean(perf.(ms).omm_bl(:,[2,4]),2) - nanmean(perf.(ms).om_bl(:,[2,4]),2) ],2);
    
    
    if strcmp(ms, 'Accuracy')
        d1 = -d1;
        d2 = -d2;
    end
    
    

    [cr, pcr] = corr(perf.(ms).d(:,1),perf.(ms).d(:,2), 'Type', 'Kendall' )
    fprintf('\ncorrelations blocked and interleaved: %.2f, p = %.2f\n', cr, pcr);

    
    fprintf('\n Intertrial effects (interleaved condition)\n')
    fprintf(' repeated dim: %.2f (sd %.2f) \n ', mean(perf.(ms).dr(:,2)), std(perf.(ms).dr(:,2)));
    fprintf('switched dim: %.2f (sd %.2f) \n ', mean(perf.(ms).ds(:,2)), std(perf.(ms).ds(:,2)));
    [h,p,ci,stats] = ttest(perf.(ms).dr(:,2) , perf.(ms).ds(:,2)); 
    fprintf('diff: %.2f (sd %.2f), [%.2f, %.2f],  t(%u) = %.2f, p = %.3f\n', mean(perf.(ms).dr(:,2) - perf.(ms).ds(:,2)),  std((perf.(ms).dr(:,2) - perf.(ms).ds(:,2))), ci, stats.df, stats.tstat, p)
end


save(sprintf('Exp_%s.mat', expname ), 'perf')

save( sprintf('Exp_%s_match_trials.mat', expname), 'acc_tr_m', 'acc_tr_mm', 'rt_tr_m', 'rt_tr_mm' );

writetable(res, sprintf('Exp_%s.csv', expname))
