Looking-to-cue-sides-and-modulation-of-value-in-OFC/code/behavioralAnalyses.m

clear all; close all; clc;

addpath('./support_functions');

figures_folder='../figures';

plotopts.Fig1B=0; % includes SuppFigS3A
plotopts.Fig1C=0; % handle memory usage
plotopts.Fig1D=0; % includes SuppFigS2
plotopts.Fig2A=0;
plotopts.Fig2B=0;
plotopts.Fig2C=0; % includes SuppFigS4
plotopts.Fig2D=0;

twinNames={'preoffer1','offer1','delay1','offer2','delay2','re-fixate','choice-go','choice-hold'};

subjsData(1).sbdata=getSubjectData(1);
subjsData(2).sbdata=getSubjectData(2);

s1gbi=getBehavioralInfos(subjsData(1).sbdata);
s2gbi=getBehavioralInfos(subjsData(2).sbdata);

s1gem1=getEyeTracks(subjsData(1).sbdata,'all','preoffer1','offer1');
s1gem2=getEyeTracks(subjsData(1).sbdata,'all','offer1','delay1');
s1gem3=getEyeTracks(subjsData(1).sbdata,'all','delay1','offer2');
s1gem4=getEyeTracks(subjsData(1).sbdata,'all','offer2','delay2');
s1gem5=getEyeTracks(subjsData(1).sbdata,'all','delay2','startfix');
s1gem6=getEyeTracks(subjsData(1).sbdata,'all','startfix','choicego');
s1gem7=getEyeTracks(subjsData(1).sbdata,'all','choicego','choicesacc');
s1gem8=getEyeTracks(subjsData(1).sbdata,'all','choicesacc','choicemade');

s2gem1=getEyeTracks(subjsData(2).sbdata,'all','preoffer1','offer1');
s2gem2=getEyeTracks(subjsData(2).sbdata,'all','offer1','delay1');
s2gem3=getEyeTracks(subjsData(2).sbdata,'all','delay1','offer2');
s2gem4=getEyeTracks(subjsData(2).sbdata,'all','offer2','delay2');
s2gem5=getEyeTracks(subjsData(2).sbdata,'all','delay2','startfix');
s2gem6=getEyeTracks(subjsData(2).sbdata,'all','startfix','choicego');
s2gem7=getEyeTracks(subjsData(2).sbdata,'all','choicego','choicesacc');
s2gem8=getEyeTracks(subjsData(2).sbdata,'all','choicesacc','choicemade');


if plotopts.Fig1B
    %% Fig. 1B  and Supp. Fig. S3A - Choice vs difference in EV
    nluqEV=-3:.05:3; % binning EV difference in 0.05 nominal units
    s1o1evsd1=s1gbi.offer1ev(s1gbi.ifLsR); s1o1evsd2=s1gbi.offer1ev(s1gbi.ifRsL);
    s1o2evsd1=s1gbi.offer2ev(s1gbi.ifLsR); s1o2evsd2=s1gbi.offer2ev(s1gbi.ifRsL);
    s2o1evsd1=s2gbi.offer1ev(s2gbi.ifLsR); s2o1evsd2=s2gbi.offer1ev(s2gbi.ifRsL);
    s2o2evsd1=s2gbi.offer2ev(s2gbi.ifLsR); s2o2evsd2=s2gbi.offer2ev(s2gbi.ifRsL);
    s1chsd1=s1gbi.choose12(s1gbi.ifLsR)==1; s1chsd2=s1gbi.choose12(s1gbi.ifRsL)==1;
    s2chsd1=s2gbi.choose12(s2gbi.ifLsR)==1; s2chsd2=s2gbi.choose12(s2gbi.ifRsL)==1;
    
    x1=[s1o2evsd1-s1o1evsd1; -s1o2evsd2+s1o1evsd2];
    y1=[1-s1chsd1; s1chsd2];
    g1=fitglm(x1,y1,'Distribution','binomial','link','logit');
    g1cci=g1.coefCI;
    bs1.p=g1.Coefficients.pValue;
    b1=g1.Coefficients.Estimate;
    bx1 = b1(1) + (x1*b1(2));
    hx1 = 1 ./ (1 + exp(-sort(bx1)));
    
    x2=[s2o2evsd1-s2o1evsd1; -s2o2evsd2+s2o1evsd2];
    y2=[1-s2chsd1; s2chsd2];
    g2=fitglm(x2,y2,'Distribution','binomial','link','logit');
    g2cci=g2.coefCI;
    bs2.p=g2.Coefficients.pValue;
    b2=g2.Coefficients.Estimate;
    bx2 = b2(1) + (x2*b2(2));
    hx2 = 1 ./ (1 + exp(-sort(bx2)));
    
    x3=[ x1; x2];
    y3=[ y1; y2];
    g3=fitglm(x3,y3,'Distribution','binomial','link','logit');
    g3cci=g3.coefCI;
    bs3.p=g3.Coefficients.pValue;
    b3=g3.Coefficients.Estimate;
    bx3 = b3(1) + (x3*b3(2));
    hx3 = 1 ./ (1 + exp(-sort(bx3)));
    
    [zq1,xq1]=nanmeanquantized(nluqEV,x1,y1);
    [zq2,xq2]=nanmeanquantized(nluqEV,x2,y2);
    [zq3,xq3]=nanmeanquantized(nluqEV,x3,y3);
    
    hfig=figure('units','normalized','outerposition',[0 0 .3 .3]);
    subplot(1,2,1);
    plot(xq3,zq3,'.','color',[0 150 0]./255); hold on;
    plot(sort(x3),hx3,'-','color',[0 150 0]./255);  plot([0 0],[0 1],'k-');
    plot([-3 0],[1 1]./(1+exp(-b3(1))),'--','color',[0 150 0]./255)
    fill([xq3 -xq3],[glmval(g3cci(:,1),xq3,'logit')' flip(glmval(g3cci(:,2),xq3,'logit')')],[0 150 0]./255,'edgecolor','none','facealpha',.15);
    xlabel('EV_R-EV_L'); ylabel('chR');
    title('Subject 1,2: all sessions'); xlim([-3 3]); daspect([6 1 1]);
    text(4,.9,sprintf(['\\beta_0 %1.4f' getpstars(bs3.p(1)) ';'],b3(1)));
    text(4,.8,sprintf(['\\beta_1 %1.4f' getpstars(bs3.p(2)) ';'],b3(2)));
    text(4,.7,sprintf('n= %d;',length(x3)));
    
    saveas(hfig,[figures_folder '/Fig1B_ratio_chosen_right_vs_offer_diff_ev.fig']);
    saveas(hfig,[figures_folder '/Fig1B_ratio_chosen_right_vs_offer_diff_ev.png']);
    saveas(hfig,[figures_folder '/Fig1B_ratio_chosen_right_vs_offer_diff_ev.svg']);
    
    hfig=figure('units','normalized','outerposition',[0 0 1 .3]);
    subplot(1,2,1);
    plot(xq1,zq1,'.','color',[250 0 0]./255); hold on;
    plot(sort(x1),hx1,'-','color',[250 0 0]./255);  plot([0 0],[0 1],'k-');
    plot([-3 0],[1 1]./(1+exp(-b1(1))),'--','color',[250 0 0]./255)
    fill([xq1 -xq1],[glmval(g1cci(:,1),xq1,'logit')' flip(glmval(g1cci(:,2),xq1,'logit')')],[250 0 0]./255,'edgecolor','none','facealpha',.15);
    xlabel('EV_R-EV_L'); ylabel('chR'); title('Subject 1: all sessions'); xlim([-3 3]); daspect([6 1 1]);
    text(4,.9,sprintf(['\\beta_0 %1.4f' getpstars(bs1.p(1)) ';'],b1(1)));
    text(4,.8,sprintf(['\\beta_1 %1.4f' getpstars(bs1.p(2)) ';'],b1(2)));
    text(4,.7,sprintf('n= %d;',length(x1)));
    
    subplot(1,2,2);
    plot(xq2,zq2,'.','color',[0 0 250]./255); hold on;
    plot(sort(x2),hx2,'-','color',[0 0 250]./255);  plot([0 0],[0 1],'k-');
    plot([-3 0],[1 1]./(1+exp(-b2(1))),'--','color',[0 0 250]./255)
    fill([xq2 -xq2],[glmval(g2cci(:,1),xq2,'logit')' flip(glmval(g2cci(:,2),xq1,'logit')')],[0 0 250]./255,'edgecolor','none','facealpha',.15);
    xlabel('EV_R-EV_L'); ylabel('chR'); title('Subject 2: all sessions'); xlim([-3 3]); daspect([6 1 1]);
    text(4,.9,sprintf(['\\beta_0 %1.4f' getpstars(bs2.p(1)) ';'],b2(1)));
    text(4,.8,sprintf(['\\beta_1 %1.4f' getpstars(bs2.p(2)) ';'],b2(2)));
    text(4,.7,sprintf('n= %d;',length(x2)));
    
    saveas(hfig,[figures_folder '/SuppFigS3A_ratio_chosen_right_vs_offer_diff_ev.fig']);
    saveas(hfig,[figures_folder '/SuppFigS3A_ratio_chosen_right_vs_offer_diff_ev.png']);
    saveas(hfig,[figures_folder '/SuppFigS3A_ratio_chosen_right_vs_offer_diff_ev.svg']);
end


if plotopts.Fig1C
    %% Fig. 1C - Histogram of eye position during task execution and count of trials LookL/LookR
    
    poolstackvars=getPoolStackVars(subjsData);
    poolstackepos=getPoolStackEyePos(subjsData);
    ncells=arrayfun(@(sn) poolstackvars(sn).ncells,1:8);
    poolstacktweyepos=getTimewinsPoolEyePos();%poolstackvars,poolstackepos);
    poolstackcctweyeposx=squeeze(cat(3,poolstacktweyepos(1).posX(:,cumsum(ncells),:),poolstacktweyepos(2).posX(:,cumsum(ncells),:),...
                                       poolstacktweyepos(3).posX(:,cumsum(ncells),:),poolstacktweyepos(4).posX(:,cumsum(ncells),:),...
                                       poolstacktweyepos(5).posX(:,cumsum(ncells),:),poolstacktweyepos(6).posX(:,cumsum(ncells),:),...
                                       poolstacktweyepos(7).posX(:,cumsum(ncells),:),poolstacktweyepos(8).posX(:,cumsum(ncells),:)));
    %[poolstackcctweyeposx, epxcovhist]=getCellsTimewinsPoolEyePos(poolstackvars,poolstacktweyepos);
    epxcovhist=squeeze(sum(sum(~isnan(poolstackcctweyeposx))))./sum(arrayfun(@(sn) size(poolstackvars(sn).offer1sd,1),1:8));
    ntptws=arrayfun(@(sn) size(poolstacktweyepos(sn).posX,3),1:8);
                  
    horiz_edges=min(poolstackcctweyeposx(:)):1e2:max(poolstackcctweyeposx(:));
    hcount=nan(length(horiz_edges)-1,sum(ntptws));
    hL=nan(1,sum(ntptws)); hR=hL;
    for tt=1:sum(ntptws)
        currepx=poolstackcctweyeposx(:,:,tt);
        hcount(:,tt)=histcounts(currepx(:),horiz_edges);
        hL(tt)=mean(currepx(~isnan(currepx))<0);
        hR(tt)=mean(currepx(~isnan(currepx))>0);
    end
    
    hfig= figure('units','normalized','position',[0 0 1 1]);
    subplot(2,1,1);
    surf(1:sum(ntptws),linspace(-1,1,length(horiz_edges)-1),hcount./max(hcount(:)),'edgecolor','none'); hold on
    for jj=1:length(twinNames); plot([1 1].*sum(ntptws(1:jj)),[-1 1],'k:'); text(sum(ntptws(1:jj-1)),1.05,twinNames{jj}); end
    view(2); colormap(1-gray(255)); caxis([0 .2]); xlim([0 sum(ntptws)]); ylim([-1 1]);
    xlabel('time aligned to epochs start (ms)'); ylabel('horizontal eye ');
    ntptws=arrayfun(@(tw) size(poolstacktweyepos(tw).posX,3), 1:8);
    
    subplot(2,1,2);
    coveragemask=nan(1,sum(ntptws)); coveragemask(epxcovhist>.999)=1;
    plot(1:sum(ntptws),hL.*coveragemask,'k','linewidth',2); hold on
    plot(1:sum(ntptws),hR.*coveragemask,'color',[.5 .5 .5],'linewidth',2);
    for jj=1:length(twinNames); plot([1 1].*sum(ntptws(1:jj)),[-1 1],'k:'); text(sum(ntptws(1:jj-1)),1.05,twinNames{jj}); end
    xlabel('time aligned to epochs start (ms)');  ylabel('% trials'); xlim([0 sum(ntptws)]); ylim([0 1]);
    legend('LookL', 'LookR');
    
    saveas(hfig,[figures_folder '/Fig1C_coverage_of_horizontal_coordinates.fig']);
    saveas(hfig,[figures_folder '/Fig1C_coverage_of_horizontal_coordinates.png']);
    saveas(hfig,[figures_folder '/Fig1C_coverage_of_horizontal_coordinates.svg']);
    
end



if plotopts.Fig1D
    %% Fig. 1D and Supp. Fig. S2 - Heatmaps
    
    rmabs=0; rmradius=0;
    dxy=1.5e2; xaxis=-(4e4+dxy/2):dxy:(4e4+dxy/2);
    xbins=-(4e4+dxy):dxy:(4e4+dxy); yaxis=xaxis; ybins=xbins; cxs=[0 .5];
    cmt=hot(255);
    
    srr={'all','bestL','bestR','chsnL','chsnR'};
    % all, chsnL, chsnR, firstL, firstR, bestL, bestR, subj1, subj2, bestLR, chsnLR
    
    hfig=figure('units','normalized','outerposition',[0 0 1 1]);
    for ssr=1:length(srr)
        strcond=srr{ssr};
        for tw=1:8
            subplot(length(srr),8,8*(ssr-1)+tw); hold on;
            if strcmpi(strcond,'all')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                s1cpx2=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                s2cpx2=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
            elseif strcmpi(strcond,'subj1')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                s1cpx2=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
                s2cpx1=[];%eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                s2cpx2=[];%eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                s2cpy1=[];%eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
                s2cpy2=[];%eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
            elseif strcmpi(strcond,'subj2')
                s1cpx1=[];%eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                s1cpx2=[];%eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                s1cpy1=[];%eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
                s1cpy2=[];%eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                s2cpx2=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
            elseif strcmpi(strcond,'bestL')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s1cpx2=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
                s2cpx2=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
            elseif strcmpi(strcond,'bestR')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s1cpx2=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s2cpx2=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
            elseif strcmpi(strcond,'bestLR')
                s1cpx1bL=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s1cpx2bL=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s1cpy1bL=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s1cpy2bL=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s2cpx1bL=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
                s2cpx2bL=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s2cpy1bL=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
                s2cpy2bL=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s1cpx1bR=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s1cpx2bR=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s1cpy1bR=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.offerRev>s1gbi.offerLev)),:))']);
                s1cpy2bR=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.offerLev>s1gbi.offerRev)),:))']);
                s2cpx1bR=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s2cpx2bR=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
                s2cpy1bR=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.offerRev>s2gbi.offerLev)),:))']);
                s2cpy2bR=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.offerLev>s2gbi.offerRev)),:))']);
            elseif strcmpi(strcond,'chsnLR')
                s1cpx1bL=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.choose12==1)),:))']);
                s1cpx2bL=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.choose12==1)),:))']);
                s1cpy1bL=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.choose12==1)),:))']);
                s1cpy2bL=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.choose12==1)),:))']);
                s2cpx1bL=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.choose12==1)),:))']);
                s2cpx2bL=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.choose12==1)),:))']);
                s2cpy1bL=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.choose12==1)),:))']);
                s2cpy2bL=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.choose12==1)),:))']);
                s1cpx1bR=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.choose12==2)),:))']);
                s1cpx2bR=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.choose12==2)),:))']);
                s1cpy1bR=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.choose12==2)),:))']);
                s1cpy2bR=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.choose12==2)),:))']);
                s2cpx1bR=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.choose12==2)),:))']);
                s2cpx2bR=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.choose12==2)),:))']);
                s2cpy1bR=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.choose12==2)),:))']);
                s2cpy2bR=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.choose12==2)),:))']);
            elseif strcmpi(strcond,'firstL')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                s1cpx2=[];%eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
                s1cpy2=[];%eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                s2cpx2=[];%eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
                s2cpy2=[];%eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
            elseif strcmpi(strcond,'firstR')
                s1cpx1=[];%eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                s1cpx2=-eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                s1cpy1=[];%eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
                s2cpx1=[];%eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                s2cpx2=-eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                s2cpy1=[];%eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
            elseif strcmpi(strcond,'chsnL')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.choose12==1)),:))']);
                s1cpx2=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.choose12==1)),:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.choose12==1)),:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.choose12==1)),:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.choose12==1)),:))']);
                s2cpx2=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.choose12==1)),:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.choose12==1)),:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.choose12==1)),:))']);
            elseif strcmpi(strcond,'chsnR')
                s1cpx1=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifLsR,find(s1gbi.choose12==2)),:))']);
                s1cpx2=eval(['(s1gem' num2str(tw) '.posX(intersect(s1gbi.ifRsL,find(s1gbi.choose12==2)),:))']);
                s1cpy1=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifLsR,find(s1gbi.choose12==2)),:))']);
                s1cpy2=eval(['(s1gem' num2str(tw) '.posY(intersect(s1gbi.ifRsL,find(s1gbi.choose12==2)),:))']);
                s2cpx1=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifLsR,find(s2gbi.choose12==2)),:))']);
                s2cpx2=eval(['(s2gem' num2str(tw) '.posX(intersect(s2gbi.ifRsL,find(s2gbi.choose12==2)),:))']);
                s2cpy1=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifLsR,find(s2gbi.choose12==2)),:))']);
                s2cpy2=eval(['(s2gem' num2str(tw) '.posY(intersect(s2gbi.ifRsL,find(s2gbi.choose12==2)),:))']);
            end
            % comment here to see subjects separately
            %s1cpx1=[]; s1cpx2=[]; s1cpy1=[]; s1cpy2=[];
            %s2cpx1=[]; s2cpx2=[]; s2cpy1=[]; s2cpy2=[];
            
            if strcmpi(strcond,'bestLR') || strcmpi(strcond,'chsnLR')
                nx11=size(s1cpx1bL,2); nx12=size(s1cpx2bL,2); nx21=size(s2cpx1bL,2); nx22=size(s2cpx2bL,2);
                nx=max([nx11 nx12 nx21 nx22]);
                ny11=size(s1cpy1bL,2); ny12=size(s1cpy2bL,2); ny21=size(s2cpy1bL,2); ny22=size(s2cpy2bL,2);
                ny=max([ny11 ny12 ny21 ny22]);
                cpxbL=[s1cpx1bL nan(size(s1cpx1bL,1),nx-nx11); -s1cpx2bL nan(size(s1cpx2bL,1),nx-nx12); ...
                    s2cpx1bL nan(size(s2cpx1bL,1),nx-nx21); -s2cpx2bL nan(size(s2cpx2bL,1),nx-nx22)];
                cpybL=[s1cpy1bL nan(size(s1cpy1bL,1),ny-ny11);  s1cpy2bL nan(size(s1cpy2bL,1),ny-ny12); ...
                    s2cpy1bL nan(size(s2cpy1bL,1),ny-ny21);  s2cpy2bL nan(size(s2cpy2bL,1),ny-ny22)];
                
                nx11=size(s1cpx1bR,2); nx12=size(s1cpx2bR,2); nx21=size(s2cpx1bR,2); nx22=size(s2cpx2bR,2);
                nx=max([nx11 nx12 nx21 nx22]);
                ny11=size(s1cpy1bR,2); ny12=size(s1cpy2bR,2); ny21=size(s2cpy1bR,2); ny22=size(s2cpy2bR,2);
                ny=max([ny11 ny12 ny21 ny22]);
                cpxbR=[s1cpx1bR nan(size(s1cpx1bR,1),nx-nx11); -s1cpx2bR nan(size(s1cpx2bR,1),nx-nx12); ...
                    s2cpx1bR nan(size(s2cpx1bR,1),nx-nx21); -s2cpx2bR nan(size(s2cpx2bR,1),nx-nx22)];
                cpybR=[s1cpy1bR nan(size(s1cpy1bR,1),ny-ny11);  s1cpy2bR nan(size(s1cpy2bR,1),ny-ny12); ...
                    s2cpy1bR nan(size(s2cpy1bR,1),ny-ny21);  s2cpy2bR nan(size(s2cpy2bR,1),ny-ny22)];
                
                if tw==1
                    cpxbL=cpxbL(:,end-400+1:end); cpybL=cpybL(:,end-400+1:end);
                    cpxbR=cpxbR(:,end-400+1:end); cpybR=cpybR(:,end-400+1:end);
                else
                    cpxbL=cpxbL(:,mean(~isnan(cpxbL),1)>=.999);
                    cpybL=cpybL(:,mean(~isnan(cpxbL),1)>=.999);
                    cpxbR=cpxbR(:,mean(~isnan(cpxbR),1)>=.999);
                    cpybR=cpybR(:,mean(~isnan(cpxbR),1)>=.999);
                end
                %s3=scatter(nanmean(cpxbL,2),nanmean(cpybL,2),'filled','SizeData',1); s3.MarkerFaceColor='r';
                %hold on; alpha(s3,.3);%p3.MarkerFaceColor=[1 1 1 ];
                %s3=scatter(nanmean(cpxbR,2),nanmean(cpybR,2),'filled','SizeData',1); s3.MarkerFaceColor='b';
                %hold off; alpha(s3,.3);%p3.MarkerFaceColor=[1 1 1 ];
                hxybL=imgaussfilt(histcounts2(cpxbL(:),cpybL(:),xbins,ybins,'Normalization','probability'),5);
                hxybR=imgaussfilt(histcounts2(cpxbR(:),cpybR(:),xbins,ybins,'Normalization','probability'),5);
                %hxybLth=normalize(hxybL,'range',[0 1]);%(hxybL./max(hxybL(:)));
                %hxybRth=normalize(hxybR,'range',[0 1]);%(hxybR./max(hxybR(:)));
                hxybLth=(hxybL./max(hxybL(:)));
                hxybRth=(hxybR./max(hxybR(:)));
                hxyz=normalize(hxybLth'-hxybRth','center',mode(hxybLth(:)'-hxybRth(:)'));
                %min(hxybLth(:))
                %max(hxybLth(:))
                sf=surf(xaxis,yaxis,hxyz,'edgecolor','none'); view(2); %daspect([1 1 1]);
                %colormap(getdivergentcolmap(1-[.2 .2 0; 1 1 0; 1 1 1; 0 1 1; 0 .2 .2],[256 64 64 256]));
                caxis([-1 1].*max(abs([min(hxyz(:)) max(hxyz(:))])));
                
            else
                nx11=size(s1cpx1,2); nx12=size(s1cpx2,2); nx21=size(s2cpx1,2); nx22=size(s2cpx2,2);
                nx=max([nx11 nx12 nx21 nx22]);
                ny11=size(s1cpy1,2); ny12=size(s1cpy2,2); ny21=size(s2cpy1,2); ny22=size(s2cpy2,2);
                ny=max([ny11 ny12 ny21 ny22]);
                cpx=[s1cpx1 nan(size(s1cpx1,1),nx-nx11); -s1cpx2 nan(size(s1cpx2,1),nx-nx12); ...
                    s2cpx1 nan(size(s2cpx1,1),nx-nx21); -s2cpx2 nan(size(s2cpx2,1),nx-nx22)];
                cpy=[s1cpy1 nan(size(s1cpy1,1),ny-ny11); s1cpy2 nan(size(s1cpy2,1),ny-ny12); ...
                    s2cpy1 nan(size(s2cpy1,1),ny-ny21); s2cpy2 nan(size(s2cpy2,1),ny-ny22)];
                %cpn=(isnan(cpx) | isnan(cpy)); cp0=(cpx==0 & cpy==0); cpx(cpn | cp0)=nan; cpy(cpn | cp0)=nan;
                
                if rmabs %cpx=rmnansandthreshold(cpx,5000,abs(cpx));
                    cpy=rmnansandthresholdreversed(cpy,6000,abs(cpy));
                    cpx=rmnansandthresholdreversed(cpx,12000,abs(cpx));
                elseif rmradius
                    cpx=rmnansandthreshold(cpx,5000,sqrt(cpx.^2+cpy.^2));
                    cpy=rmnansandthreshold(cpy,5000,sqrt(cpx.^2+cpy.^2));
                end
                
                if tw==1
                    cpx=cpx(:,end-400+1:end); cpy=cpy(:,end-400+1:end);
                else
                    cpx=cpx(:,mean(~isnan(cpx),1)>=.999);
                    cpy=cpy(:,mean(~isnan(cpx),1)>=.999);
                end
                hxy=imgaussfilt(histcounts2(cpx(:),cpy(:),xbins,ybins,'Normalization','probability'),5);
                hxyth=(hxy./max(hxy(:)));
                surf(xaxis,yaxis,hxyth','edgecolor','none'); view(2); %daspect([1 1 1]);
                %caxis([0 1].*max(abs(max(hxy(:)))));
                colormap(cmt); caxis(cxs);
            end
            
            %hold on; xlim([-4e4 4e4]); ylim([-4e4 4e4]); daspect([1 1 1]);
            title([twinNames{tw} ', ' strcond]); hold on;
            %plot([-4e4 4e4],[0 0],'w:'); plot([0 0],[-4e4 4e4],'w:');
            xlim([-4e4 4e4]/2); ylim([-4e4 4e4]/2); daspect([1 1 1]);
            %plot3([-.82e4 .82e4], [0 0], [1 1], 'kx');
            %   pl1=plot3(-.81395e4+[-1 -1 1 1 -1]*.2e4,[-1 1 1 -1 -1]*.8e4,-[1 1 1 1 1],'w:'); % stimulus frame borders
            %   pl2=plot3(+.81395e4+[-1 -1 1 1 -1]*.2e4,[-1 1 1 -1 -1]*.8e4,-[1 1 1 1 1],'w:'); % stimulus frame borders
            %plot3([-.85e4-(1.7e3) .85e4+1.7e3], [0 0], [1 1], 'kx');
            %plot3(-.81395e4,0,1,'kx');
            %plot3(+.81395e4,0,1,'kx');
            %plot(-0.6e4,-.8e4,'rx')
            %plot(-1.0e4,-.8e4,'rx')
            %plot(+0.6e4+.2e4,-.8e4,'bx')
            %plot(+1.0e4-.2e4,-.8e4,'bx')
            % x:max(cpx(:))=10.995:.81395e4
            % max(cpx(:)) : 43.98 = 10.995 : .81395e4
            %pl1b=plot3(-.85e4+[-1 -1 1 1 -1]*4.08*8500/(27.3),[-1 1 1 -1 -1]*11.35*8500/(27.3),-[1 1 1 1 1],'w:');
            %pl2b=plot3(+.85e4+[-1 -1 1 1 -1]*4.08*8500/(27.3),[-1 1 1 -1 -1]*11.35*8500/(27.3),-[1 1 1 1 1],'w:');
        end
        cb=colorbar(); cb.Position(1)=.96; daspect([1 1 1]);
    end
    saveas(hfig,[figures_folder '/Fig1D-Eyetrack-heatmaps.png']);
    saveas(hfig,[figures_folder '/Fig1D-Eyetrack-heatmaps.fig']);
    saveas(hfig,[figures_folder '/Fig1D-Eyetrack-heatmaps.svg']);
    
 end
 
 if plotopts.Fig2A
    %% Fig. 2A - EVdiff vs chR for fR>0.5 and fR<0.5
    
    hfig=figure('units','normalized','position',[0 0 1 .2]);
    for tw=1:8
        subplot(1,8,tw); hold on;
        s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
        s1cpx2=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
        s1cpy1=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
        s1cpy2=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
        s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
        s2cpx2=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
        s2cpy1=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
        s2cpy2=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
        
        nx11=size(s1cpx1,2); nx12=size(s1cpx2,2); nx21=size(s2cpx1,2); nx22=size(s2cpx2,2);
        nx=max([nx11 nx12 nx21 nx22]);
        ny11=size(s1cpy1,2); ny12=size(s1cpy2,2); ny21=size(s2cpy1,2); ny22=size(s2cpy2,2);
        ny=max([ny11 ny12 ny21 ny22]);
        cpx=[s1cpx1 nan(size(s1cpx1,1),nx-nx11); -s1cpx2 nan(size(s1cpx2,1),nx-nx12); ...
            s2cpx1 nan(size(s2cpx1,1),nx-nx21); -s2cpx2 nan(size(s2cpx2,1),nx-nx22)];
        cpy=[s1cpy1 nan(size(s1cpy1,1),ny-ny11); s1cpy2 nan(size(s1cpy2,1),ny-ny12); ...
            s2cpy1 nan(size(s2cpy1,1),ny-ny21); s2cpy2 nan(size(s2cpy2,1),ny-ny22)];
        cpn=(isnan(cpx) | isnan(cpy)); cp0=(cpx==0 & cpy==0); cpx(cpn | cp0)=nan; cpy(cpn | cp0)=nan;
        
        if rmabs; cpx=rmnansandthreshold(cpx,5000,abs(cpx));
        elseif rmradius; cpx=rmnansandthreshold(cpx,5000,sqrt(cpx.^2+cpy.^2)); else; end
        
        if tw==1
            cpx=cpx(:,end-400+1:end); cpy=cpy(:,end-400+1:end);
        else
            cpx=cpx(:,mean(~isnan(cpx),1)>=.999);
            cpy=cpy(:,mean(~isnan(cpx),1)>=.999);
        end
        
        sum_L=sum(~isnan(rmnansandthreshold(-cpx,0)),2);
        sum_R=sum(~isnan(rmnansandthreshold(+cpx,0)),2);
        frac_L=sum_L./(sum_L+sum_R); frac_R=1-frac_L;
        x1=[s1gbi.offer2ev(s1gbi.ifLsR,:)-s1gbi.offer1ev(s1gbi.ifLsR,:);...
            s1gbi.offer2ev(s1gbi.ifRsL,:)-s1gbi.offer1ev(s1gbi.ifRsL,:);...
            s2gbi.offer2ev(s2gbi.ifLsR,:)-s2gbi.offer1ev(s2gbi.ifLsR,:);...
            s2gbi.offer2ev(s2gbi.ifRsL,:)-s2gbi.offer1ev(s2gbi.ifRsL,:)];
        %[s1o2evsd1-s1o1evsd1; s1os2evsd2-s1o1evsd2];
        y1=[s1gbi.choose12(s1gbi.ifLsR,:)==2; s1gbi.choose12(s1gbi.ifRsL,:)==2;
            s2gbi.choose12(s2gbi.ifLsR,:)==2; s2gbi.choose12(s2gbi.ifRsL,:)==2];
        
        x1_gR=x1(frac_R>0.5); y1_gR=y1(frac_R>0.5);
        x1_lR=x1(frac_R<0.5); y1_lR=y1(frac_R<0.5);
        
        %y1=double(frac_R>0.5); %nluqEV=-3:.2:3;%linspace(-max(abs(x1)),max(abs(x1)),150);
        
        %[zq1,xq1]=nanmeanquantized(nluqEV,x1,y1); plot(xq1,zq1,'.','color',[0 0 .75]);
        %hold on; [sq1,st1,fsigm]=sigm_fit(xq1,zq1,[0 1 NaN NaN],[],0);
        %yq1=st1.ypred; [~,im1]=min(abs(xq1));     plot([-3 xq1(im1)],yq1(im1).*[1 1],'k:','linewidth',1.5);
        %plot(linspace(-3,3,length(yq1)),st1.ypred,'b-','linewidth',1.5);
        %fill([linspace(-3,3,length(yq1))'; linspace(3,-3,length(yq1))'],[st1.ypredupperCI; flipud(st1.ypredlowerCI)],'b','edgecolor','none','facealpha',.3);
        %NuqEV=31;
        %nluqEV=linspace(-max(abs(x1)),max(abs(x1)),NuqEV);
        nluqEV=-3:.2:3;
        [tested,]=nanhistquantized(nluqEV,x1,ones(1,length(y1)));
        [frac_R_g05,uq]=nanhistquantized(nluqEV,x1,y1);
        stst=statset('fitglm');
        stst.tolX=1e-20;
        logitMdl= fitglm(x1,y1,'Distribution','binomial','Link','logit');
        logitCci = coefCI(logitMdl);
        logitFit = glmval(logitMdl.Coefficients.Estimate,uq,'logit'); % logitFit=1./(1+exp(-logitCoeffs(1)-logitCoeffs(2).*uq));
        logitFitCiu = glmval(logitCci(:,1),uq,'logit');
        logitFitCil = glmval(logitCci(:,2),uq,'logit');
        
        [tested_gR,]=nanhistquantized(nluqEV,x1_gR,ones(1,length(y1_gR)));
        [frac_R_g05_gR,uq]=nanhistquantized(nluqEV,x1_gR,y1_gR);
        stst=statset('fitglm');
        stst.tolX=1e-20;
        logitMdl_gR= fitglm(x1_gR,y1_gR,'Distribution','binomial','Link','logit');
        logitCci_gR = coefCI(logitMdl_gR);
        logitFit_gR = glmval(logitMdl_gR.Coefficients.Estimate,uq,'logit'); % logitFit=1./(1+exp(-logitCoeffs(1)-logitCoeffs(2).*uq));
        logitFitCiu_gR = glmval(logitCci_gR(:,1),uq,'logit');
        logitFitCil_gR = glmval(logitCci_gR(:,2),uq,'logit');
        
        [tested_lR,]=nanhistquantized(nluqEV,x1_lR,ones(1,length(y1_lR)));
        [frac_R_g05_lR,uq]=nanhistquantized(nluqEV,x1_lR,y1_lR);
        stst=statset('fitglm');
        stst.tolX=1e-20;
        logitMdl_lR= fitglm(x1_lR,y1_lR,'Distribution','binomial','Link','logit');
        logitCci_lR = coefCI(logitMdl_lR);
        logitFit_lR = glmval(logitMdl_lR.Coefficients.Estimate,uq,'logit'); % logitFit=1./(1+exp(-logitCoeffs(1)-logitCoeffs(2).*uq));
        logitFitCiu_lR = glmval(logitCci_lR(:,1),uq,'logit');
        logitFitCil_lR = glmval(logitCci_lR(:,2),uq,'logit');
        
        
        plot(uq,frac_R_g05_gR./tested_gR,'.','color',[0 0 50]./255); hold on;
        plot([-1e3 1e3],[0 0],'k-',[0 0],[-10 10],'k-');
        plot(uq,logitFit_gR,'-','color',[0 0 50]./255);
        fill([uq -uq],[logitFitCiu_gR' flip(logitFitCil_gR')],[0 0 50]./255,'edgecolor','none','facealpha',.15);
        text(abs(max(x1))*.6,.75,sprintf(['\\beta_1 %1.2f' getpstars(logitMdl_gR.Coefficients.pValue(2)) ';'],...
            logitMdl_gR.Coefficients.Estimate(2)),'color',[0 0 50]./255);
        
        plot(uq,frac_R_g05./tested,'.','color',[0 0 150]./255); hold on;
        plot([-1e3 1e3],[0 0],'k-',[0 0],[-10 10],'k-');
        plot(uq,logitFit,'-','color',[0 0 150]./255);
        fill([uq -uq],[logitFitCiu' flip(logitFitCil')],[0 0 150]./255,'edgecolor','none','facealpha',.15);
        text(abs(max(x1))*.6,.65,sprintf(['\\beta_1 %1.2f' getpstars(logitMdl.Coefficients.pValue(2)) ';'],...
            logitMdl.Coefficients.Estimate(2)),'color',[0 0 150]./255);
        
        plot(uq,frac_R_g05_lR./tested_lR,'.','color',[0 0 250]./255); hold on;
        plot([-1e3 1e3],[0 0],'k-',[0 0],[-10 10],'k-');
        plot(uq,logitFit_lR,'-','color',[0 0 250]./255);
        fill([uq -uq],[logitFitCiu_lR' flip(logitFitCil_lR')],[0 0 250]./255,'edgecolor','none','facealpha',.15);
        text(abs(max(x1))*.6,.55,sprintf(['\\beta_1 %1.2f' getpstars(logitMdl_lR.Coefficients.pValue(2)) ';'],...
            logitMdl_lR.Coefficients.Estimate(2)),'color',[0 0 250]./255);
        
        plot([0 0],[0 1],'k-'); gch=gca; gch.XTick=-3:3;
        xlabel('EV_R - EV_L');
        xlim([-3 3]); ylim([0 1]); daspect([4 1 1]);
        title([twinNames{tw} ' time']);
        if tw==1; ylabel('average chR'); end
    end
    supertitle('Subjects 1:2, all sessions');
    
    print(hfig,[figures_folder '/Fig2A-Eyetrack-chR-vs-EVdiff.png'],'-dpng','-r600');
    saveas(hfig,[figures_folder '/Fig2A-Eyetrack-chR-vs-EVdiff.fig']);
    saveas(hfig,[figures_folder '/Fig2A-Eyetrack-chR-vs-EVdiff.svg']);
end


if plotopts.Fig2B
    %% Fig 2B - fR vs chR for EVR>EVL and EVR<EVL
    hfig=figure('units','normalized','position',[0 0 1 .2]);
    
    for trgrp=1:3 %trgrp 1=bestL; 2=all trials; 3=bestR
        for tw=1:8
            subplot(1,8,tw); hold on;
            s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
            s1cpx2=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
            s1cpy1=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
            s1cpy2=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
            s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
            s2cpx2=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
            s2cpy1=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
            s2cpy2=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
            
            chosenR_flag=double([s1gbi.choose12(s1gbi.ifLsR)==2; s1gbi.choose12(s1gbi.ifRsL)==2;...
                s2gbi.choose12(s2gbi.ifLsR)==2; s2gbi.choose12(s2gbi.ifRsL)==2]);
            offLev=double([s1gbi.offer1ev(s1gbi.ifLsR); s1gbi.offer2ev(s1gbi.ifRsL);...
                s2gbi.offer1ev(s2gbi.ifLsR); s2gbi.offer2ev(s2gbi.ifRsL)]);
            offRev=double([s1gbi.offer2ev(s1gbi.ifLsR); s1gbi.offer1ev(s1gbi.ifRsL);...
                s2gbi.offer2ev(s2gbi.ifLsR); s2gbi.offer1ev(s2gbi.ifRsL)]);
            
            nx11=size(s1cpx1,2); nx12=size(s1cpx2,2); nx21=size(s2cpx1,2); nx22=size(s2cpx2,2);
            nx=max([nx11 nx12 nx21 nx22]);
            ny11=size(s1cpy1,2); ny12=size(s1cpy2,2); ny21=size(s2cpy1,2); ny22=size(s2cpy2,2);
            ny=max([ny11 ny12 ny21 ny22]);
            cpx=[s1cpx1 nan(size(s1cpx1,1),nx-nx11); -s1cpx2 nan(size(s1cpx2,1),nx-nx12); ...
                s2cpx1 nan(size(s2cpx1,1),nx-nx21);  -s2cpx2 nan(size(s2cpx2,1),nx-nx22)];
            cpy=[s1cpy1 nan(size(s1cpy1,1),ny-ny11); s1cpy2 nan(size(s1cpy2,1),ny-ny12); ...
                s2cpy1 nan(size(s2cpy1,1),ny-ny21); s2cpy2 nan(size(s2cpy2,1),ny-ny22)];
            cpn=(isnan(cpx) | isnan(cpy)); cp0=(cpx==0 & cpy==0); cpx(cpn | cp0)=nan; cpy(cpn | cp0)=nan;
            if tw==1
                cpx=cpx(:,end-400+1:end);
                cpy=cpy(:,end-400+1:end);
            else
                cpx=cpx(:,mean(~isnan(cpx),1)>=.999);
                cpy=cpy(:,mean(~isnan(cpx),1)>=.999);
            end
            
            sum_L=sum(~isnan(rmnansandthreshold(-cpx,0)),2);
            sum_R=sum(~isnan(rmnansandthreshold(+cpx,0)),2);
            
            x1=((sum_R)-(sum_L))./((((sum_R)+(sum_L)))+eps);
            nluqEV=-1:.1:1;
            
            if trgrp==3
                x1(offRev>offLev)=[]; % bestR
                chosenR_flag(offRev>offLev)=[]; % bestR
            elseif trgrp==1
                x1(offLev>offRev)=[]; % bestL
                chosenR_flag(offLev>offRev)=[]; % bestL
            end
            
            [tested,]=nanhistquantized(nluqEV,x1,ones(1,length(chosenR_flag)));
            [chosen_R,uq]=nanhistquantized(nluqEV,x1,chosenR_flag);
            
            logitMdl= fitglm(x1,chosenR_flag,'Distribution','binomial','Link','logit');
            
            logitCci = coefCI(logitMdl);
            logitFit = glmval(logitMdl.Coefficients.Estimate,uq,'logit'); % logitFit=1./(1+exp(-logitCoeffs(1)-logitCoeffs(2).*uq));
            logitFitCiu = glmval(logitCci(:,1),uq,'logit');
            logitFitCil = glmval(logitCci(:,2),uq,'logit');
            
            plot(uq,chosen_R./tested,'.','color',(0.3.*trgrp.*[0 176 80]./255)); hold on;
            plot([-1e3 1e3],[0 0],'k-',[0 0],[-10 10],'k-');
            plot(uq,logitFit,'-','color',.3.*trgrp.*[0 176 80]./255);
            fill([uq -uq],[logitFitCiu' flip(logitFitCil')],.3.*trgrp.*[0 176 80]./255,'edgecolor','none','facealpha',.15);
            xlabel('t_R-t_L'); daspect([1.5 1 1]);
            xlim([-1 1].*abs(max(x1))); ylim([0 1]); %daspect([2*abs(max(x1)) 1 1]);
            title([twinNames{tw} ' time']); if tw>1; gch=gca; gch.XTick=[]; end
            text(abs(max(x1))*.6,1.45-.15.*trgrp,sprintf(['\\beta_1 %1.2f' getpstars(logitMdl.Coefficients.pValue(2)) ';'],...
                logitMdl.Coefficients.Estimate(2)),'color',.3.*trgrp.*[0 176 80]./255);
            
            if tw==1
                text(1.1,.15,sprintf('nchR = %.0d;',sum(chosenR_flag)));
                text(1.1,0,sprintf('chR = %.0f%%;',mean(chosenR_flag).*100));
                ylabel('P(choosing R)'); else;  gch=gca; gch.YTick=[];
            end
        end
        if trgrp==1;  supertitle('Subjects 1:2, all sessions, all trials'); end
    end
    
    print(hfig,[figures_folder '/Fig2B-Eyetrack-chR-vs-timediff.png'],'-dpng','-r600');
    saveas(hfig,[figures_folder '/Fig2B-Eyetrack-chR-vs-timediff.fig']);
    saveas(hfig,[figures_folder '/Fig2B-Eyetrack-chR-vs-timediff.svg']);
    
end

if plotopts.Fig2C
    %% Fig2C - Logistic GLM of choice
    
    for isweightnormalized=0%:1
        for focusaroundstimuli=0:1 %uncomment here for supp. fig. S4
            hfig=figure('units','normalized','outerposition',[0 0 1 .75]);
            for ss=1:3
                btw=nan(6,8);
                ptw=nan(6,8);
                bmax=[]; bmin=[];
                for tw=1:8 %7
                    
                    evL=[s1gbi.offerLev(s1gbi.ifLsR); s1gbi.offerLev(s1gbi.ifRsL);      s2gbi.offerLev(s2gbi.ifLsR); s2gbi.offerLev(s2gbi.ifRsL)];
                    evR=[s1gbi.offerRev(s1gbi.ifLsR); s1gbi.offerRev(s1gbi.ifRsL);      s2gbi.offerRev(s2gbi.ifLsR); s2gbi.offerRev(s2gbi.ifRsL)];
                    varL=[s1gbi.offerLvar(s1gbi.ifLsR); s1gbi.offerLvar(s1gbi.ifRsL);   s2gbi.offerLvar(s2gbi.ifLsR); s2gbi.offerLvar(s2gbi.ifRsL)];
                    varR=[s1gbi.offerRvar(s1gbi.ifLsR); s1gbi.offerRvar(s1gbi.ifRsL);   s2gbi.offerRvar(s2gbi.ifLsR); s2gbi.offerRvar(s2gbi.ifRsL)];
                    sideRL=[ones(length(s1gbi.ifLsR),1); -ones(length(s1gbi.ifRsL),1);  ones(length(s2gbi.ifLsR),1);  -ones(length(s2gbi.ifRsL),1)];
                    
                    s1sn1gbi=getBehavioralInfos(subjsData(1).sbdata,1); s1sn2gbi=getBehavioralInfos(subjsData(1).sbdata,2);
                    s1sn3gbi=getBehavioralInfos(subjsData(1).sbdata,3); s1sn4gbi=getBehavioralInfos(subjsData(1).sbdata,4);
                    s2sn1gbi=getBehavioralInfos(subjsData(2).sbdata,1); s2sn2gbi=getBehavioralInfos(subjsData(2).sbdata,2);
                    s2sn3gbi=getBehavioralInfos(subjsData(2).sbdata,3); s2sn4gbi=getBehavioralInfos(subjsData(2).sbdata,4);
                    
                    s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                    s1cpx2=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                    s1cpy1=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifLsR,:))']);
                    s1cpy2=eval(['(s1gem' num2str(tw) '.posY(s1gbi.ifRsL,:))']);
                    s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                    s2cpx2=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                    s2cpy1=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifLsR,:))']);
                    s2cpy2=eval(['(s2gem' num2str(tw) '.posY(s2gbi.ifRsL,:))']);
                    s1cpx=[s1cpx1; s1cpx2]; s1cpy=[s1cpy1; s1cpy2];
                    s2cpx=[s2cpx1; s2cpx2]; s2cpy=[s2cpy1; s2cpy2];
                    
                    
                    if tw==1
                        s1cpx=s1cpx(:,end-400+1:end);
                        s1cpy=s1cpy(:,end-400+1:end);
                        s2cpx=s2cpx(:,end-400+1:end);
                        s2cpy=s2cpy(:,end-400+1:end);
                    else
                        s1cpx=s1cpx(:,mean(~isnan(s1cpx),1)>=.999);
                        s2cpx=s2cpx(:,mean(~isnan(s2cpx),1)>=.999);
                        s1cpy=s1cpy(:,mean(~isnan(s1cpx),1)>=.999);
                        s2cpy=s2cpy(:,mean(~isnan(s2cpx),1)>=.999);
                    end
                    
                    if focusaroundstimuli
                        s1cpx_Rts=rmnansandwithin(s1cpx,[7.5e3,+10e3]);
                        s1cpy_Rts=nan(size(s1cpx_Rts)); s1cpy_Rts(~isnan(s1cpx_Rts))=s1cpy(~isnan(s1cpx_Rts));
                        s1cpy_Rts=rmnansandwithin(s1cpy_Rts,[-8e3 8e3]);
                        s1cpx_Rts(isnan(s1cpy_Rts))=s1cpy_Rts(isnan(s1cpy_Rts));
                        
                        s1cpx_Lts=rmnansandwithin(s1cpx,[-10e3,-7.5e3]);
                        s1cpy_Lts=nan(size(s1cpx_Lts)); s1cpy_Lts(~isnan(s1cpx_Lts))=s1cpy(~isnan(s1cpx_Lts));
                        s1cpy_Lts=rmnansandwithin(s1cpy_Lts,[-8e3 8e3]);
                        s1cpx_Lts(isnan(s1cpy_Lts))=s1cpy_Lts(isnan(s1cpy_Lts));
                        
                        s2cpx_Rts=rmnansandwithin(s2cpx,[7.5e3,+10e3]);
                        s2cpy_Rts=nan(size(s2cpx_Rts)); s2cpy_Rts(~isnan(s2cpx_Rts))=s2cpy(~isnan(s2cpx_Rts));
                        s2cpy_Rts=rmnansandwithin(s2cpy_Rts,[-8e3 8e3]);
                        s2cpx_Rts(isnan(s2cpy_Rts))=s2cpy_Rts(isnan(s2cpy_Rts));
                        
                        s2cpx_Lts=rmnansandwithin(s2cpx,[-10e3,-7.5e3]);
                        s2cpy_Lts=nan(size(s2cpx_Lts)); s2cpy_Lts(~isnan(s2cpx_Lts))=s2cpy(~isnan(s2cpx_Lts));
                        s2cpy_Lts=rmnansandwithin(s2cpy_Lts,[-8e3 8e3]);
                        s2cpx_Lts(isnan(s2cpy_Lts))=s2cpy_Lts(isnan(s2cpy_Lts));
                        
                        s1cpx=nan(size(s1cpx)); s1cpx(~isnan(s1cpx_Rts))=s1cpx_Rts(~isnan(s1cpx_Rts)); s1cpx(~isnan(s1cpx_Lts))=s1cpx_Lts(~isnan(s1cpx_Lts));
                        s1cpy=nan(size(s1cpy)); s1cpy(~isnan(s1cpy_Rts))=s1cpy_Rts(~isnan(s1cpy_Rts)); s1cpx(~isnan(s1cpy_Lts))=s1cpy_Lts(~isnan(s1cpy_Lts));
                        s2cpx=nan(size(s2cpx)); s2cpx(~isnan(s2cpx_Rts))=s2cpx_Rts(~isnan(s2cpx_Rts)); s2cpx(~isnan(s2cpx_Lts))=s2cpx_Lts(~isnan(s2cpx_Lts));
                        s2cpy=nan(size(s2cpy)); s2cpy(~isnan(s2cpy_Rts))=s2cpy_Rts(~isnan(s2cpy_Rts)); s2cpx(~isnan(s2cpy_Lts))=s2cpy_Lts(~isnan(s2cpy_Lts));
                    end
                    
                    s1cpn=(isnan(s1cpx) | isnan(s1cpy)); s1cp0=(s1cpx==0 & s1cpy==0); s1cpx(s1cpn | s1cp0)=nan; s1cpy(s1cpn | s1cp0)=nan;
                    s1sum_L=sum(~isnan(rmnansandthreshold(-s1cpx,0)),2);
                    s1sum_R=sum(~isnan(rmnansandthreshold(+s1cpx,0)),2);
                    s1frac_R=s1sum_R./(s1sum_R+s1sum_L+eps);
                    
                    s2cpn=(isnan(s2cpx) | isnan(s2cpy)); s2cp0=(s2cpx==0 & s2cpy==0); s2cpx(s2cpn | s2cp0)=nan; s2cpy(s2cpn | s2cp0)=nan;
                    s2sum_L=sum(~isnan(rmnansandthreshold(-s2cpx,0)),2);
                    s2sum_R=sum(~isnan(rmnansandthreshold(+s2cpx,0)),2);
                    s2frac_R=s2sum_R./(s2sum_R+s2sum_L+eps);
                    
                    %frac_R=[s1frac_R; s2frac_R];
                    frac_R=[mat2gray(s1frac_R); mat2gray(s2frac_R)];
                    
                    chosenR_flag=double([s1gbi.choose12(s1gbi.ifLsR)==2; s1gbi.choose12(s1gbi.ifRsL)==1;...
                        s2gbi.choose12(s2gbi.ifLsR)==2; s2gbi.choose12(s2gbi.ifRsL)==1]);
                    if isweightnormalized
                        %regressVars=[evL evR varL varR sideRL prevexprew frac_R];
                        regressVars=[evL evR varL varR sideRL frac_R];
                        %regressVars=[mat2gray(evL) mat2gray(evR) mat2gray(varL) mat2gray(varR) (sideRL) (prevexprew) mat2gray(frac_R)];
                        %remove first trials of each session they dont have a prev reward
                        %s12_rmtrs=[s1_rmtrs length(s1sorted)+s2_rmtrs];
                        %regressVars(s12_rmtrs(:),:)=[];
                        %chosenR_flag(s12_rmtrs(:))=[];
                    else
                        regressVars=[evL./max(evL) evR./max(evR) varL./max(varL) varR./max(varR) sideRL frac_R./max(frac_R)];
                        %regressVars=[mat2gray(evL) mat2gray(evR) mat2gray(varL) mat2gray(varR) mat2gray(frac_R)];
                    end
                    
                    if ss==1
                        %remove second subj to check just first
                        regressVars(size(s1cpx,1)+1:end,:)=[];
                        chosenR_flag(size(s1cpx,1)+1:end)=[];
                    elseif ss==2
                        %remove first subj to check just second
                        regressVars(1:size(s1cpx,1),:)=[];
                        chosenR_flag(1:size(s1cpx,1))=[];
                    end
                    
                    %[b,dev,stats]=glmfit(regressVars,chosenR_flag);
                    [b,dev,stats]=glmfit(regressVars,chosenR_flag,'binomial','link','logit');
                    btw(:,tw)=stats.beta(2:end);
                    ptw(:,tw)=stats.p(2:end);
                    % logit(p(choose R)) = beta*X;  % might cause singularities
                    %[b,dev,stats]=glmfit(regressVars,ratioChoiceR,'binomial','link','logit');
                    
                    subplot(3,8,(ss-1)*8+tw);
                    bar(1:length(b)-1,stats.beta(2:end));
                    bmax=max([bmax; b(2:end)]); bmin=min([bmin; b(2:end)]);
                    if ss<3; title(['Subj' num2str(ss) ' - ' twinNames{tw}]);
                    else; title(['Subjs1:2 - ' twinNames{tw}]); end
                    gch=gca; gch.XTick=1:length(b)-1;
                    gch.XTickLabel={'EV_L','EV_R','\sigma^2_L','\sigma^2_R','s_L_R','f_R'};
                    gch.XTickLabelRotation=90;
                    if tw==1; ylabel('glm weights'); else;  gch=gca; gch.YTick=[]; end
                    
                end
                [~,~,~,ptwfdr]=fdr_bh(ptw);
                for tw=1:8
                    subplot(3,8,(ss-1)*8+tw);
                    for jj=1:size(btw,1)
                        text(jj-.5,btw(jj,tw)+.1*(btw(jj,tw)>0)-0.1*(btw(jj,tw)<0),getpstars(ptwfdr(jj,tw)));
                    end
                    ylim([bmin bmax].*1.25);
                end
            end
            if isweightnormalized
                isspstr='nonnrm';
            else
                isspstr='maxnrm';
            end
            if focusaroundstimuli
                focstr='isfocusedonst';
            else
                focstr='nonfocusedonst';
            end
            supertitle(['Subjects 1, 2, 1:2, all sessions, focus: ' num2str(focusaroundstimuli) ', norm weights: ' num2str(isweightnormalized)]);
            print(hfig,[figures_folder '/Fig2C-Behavioral-Choice-GLM-' isspstr '-' focstr '.svg'],'-dsvg','-r600');
            print(hfig,[figures_folder '/Fig2C-Behavioral-Choice-GLM-' isspstr '-' focstr '.png'],'-dpng','-r600');
            saveas(hfig,[figures_folder '/Fig2C-Behavioral-Choice-GLM-' isspstr '-' focstr '.fig']);
        end
    end
end

if plotopts.Fig2D
    %% Fig 2D - Logistic GLM of choice with epochs combined
    hfig=figure('units','normalized','position',[0 0 1 1]);
    for inclsideLR=0:1
        if ~(trpool<3 && inclsideLR==1)
            for focusaroundstimuli=0:1
                
                fracR_off1_del1_off2_del2=nan(length([s1gbi.offer1ev; s2gbi.offer1ev]),4);
                
                for tw=2:5
                    s1cpx1=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifLsR,:))']);
                    s1cpx2=eval(['(s1gem' num2str(tw) '.posX(s1gbi.ifRsL,:))']);
                    s2cpx1=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifLsR,:))']);
                    s2cpx2=eval(['(s2gem' num2str(tw) '.posX(s2gbi.ifRsL,:))']);
                    
                    
                    s1cpx1=s1cpx1(:,mean(~isnan(s1cpx1),1)>=.999); s1cpx2=s1cpx2(:,mean(~isnan(s1cpx2),1)>=.999);
                    s2cpx1=s2cpx1(:,mean(~isnan(s2cpx1),1)>=.999); s2cpx2=s2cpx2(:,mean(~isnan(s2cpx2),1)>=.999);
                    
                    s1cpx=[s1cpx1(:,1:min(size(s1cpx1,2),size(s1cpx2,2))); -s1cpx2(:,1:min(size(s1cpx1,2),size(s1cpx2,2)))];
                    s2cpx=[s2cpx1(:,1:min(size(s2cpx1,2),size(s2cpx2,2))); -s2cpx2(:,1:min(size(s2cpx1,2),size(s2cpx2,2)))];
                    s1cpx=s1cpx(:,mean(~isnan(s1cpx),1)>=.999);
                    s2cpx=s2cpx(:,mean(~isnan(s2cpx),1)>=.999);
                    
                    %cpx=[s1cpx; s2cpx];
                    %if tw==2 || tw==4
                    %    %    cpx=[s1cpx(:,end-200+1:end); s2cpx(:,end-200+1:end)];
                    %    s1cpx=s1cpx(:,end-200+1:end);
                    %    s2cpx=s2cpx(:,end-200+1:end);
                    %else
                    %    %    cpx=[s1cpx(:,1:200); s2cpx(:,1:200)];
                    %    s1cpx=s1cpx(:,1:200);
                    %    s2cpx=s2cpx(:,1:200);
                    %end
                    
                    if focusaroundstimuli
                        s1cpx_Rts=rmnansandwithin(s1cpx,[7.5e3,+10e3]);
                        %s1cpy_Rts=nan(size(s1cpx_Rts)); s1cpy_Rts(~isnan(s1cpx_Rts))=s1cpy(~isnan(s1cpx_Rts));
                        %s1cpy_Rts=rmnansandwithin(s1cpy_Rts,[-8e3 8e3]);
                        %s1cpx_Rts(isnan(s1cpy_Rts))=s1cpy_Rts(isnan(s1cpy_Rts));
                        
                        s1cpx_Lts=rmnansandwithin(s1cpx,[-10e3,-7.5e3]);
                        %s1cpy_Lts=nan(size(s1cpx_Lts)); s1cpy_Lts(~isnan(s1cpx_Lts))=s1cpy(~isnan(s1cpx_Lts));
                        %s1cpy_Lts=rmnansandwithin(s1cpy_Lts,[-8e3 8e3]);
                        %s1cpx_Lts(isnan(s1cpy_Lts))=s1cpy_Lts(isnan(s1cpy_Lts));
                        
                        s2cpx_Rts=rmnansandwithin(s2cpx,[7.5e3,+10e3]);
                        %s2cpy_Rts=nan(size(s2cpx_Rts)); s2cpy_Rts(~isnan(s2cpx_Rts))=s2cpy(~isnan(s2cpx_Rts));
                        %s2cpy_Rts=rmnansandwithin(s2cpy_Rts,[-8e3 8e3]);
                        %s2cpx_Rts(isnan(s2cpy_Rts))=s2cpy_Rts(isnan(s2cpy_Rts));
                        
                        s2cpx_Lts=rmnansandwithin(s2cpx,[-10e3,-7.5e3]);
                        %s2cpy_Lts=nan(size(s2cpx_Lts)); s2cpy_Lts(~isnan(s2cpx_Lts))=s2cpy(~isnan(s2cpx_Lts));
                        %s2cpy_Lts=rmnansandwithin(s2cpy_Lts,[-8e3 8e3]);
                        %s2cpx_Lts(isnan(s2cpy_Lts))=s2cpy_Lts(isnan(s2cpy_Lts));
                        
                        s1cpx=nan(size(s1cpx)); s1cpx(~isnan(s1cpx_Rts))=s1cpx_Rts(~isnan(s1cpx_Rts)); s1cpx(~isnan(s1cpx_Lts))=s1cpx_Lts(~isnan(s1cpx_Lts));
                        %s1cpy=nan(size(s1cpy)); s1cpy(~isnan(s1cpy_Rts))=s1cpy_Rts(~isnan(s1cpy_Rts)); s1cpx(~isnan(s1cpy_Lts))=s1cpy_Lts(~isnan(s1cpy_Lts));
                        s2cpx=nan(size(s2cpx)); s2cpx(~isnan(s2cpx_Rts))=s2cpx_Rts(~isnan(s2cpx_Rts)); s2cpx(~isnan(s2cpx_Lts))=s2cpx_Lts(~isnan(s2cpx_Lts));
                        %s2cpy=nan(size(s2cpy)); s2cpy(~isnan(s2cpy_Rts))=s2cpy_Rts(~isnan(s2cpy_Rts)); s2cpx(~isnan(s2cpy_Lts))=s2cpy_Lts(~isnan(s2cpy_Lts));
                    end
                    
                    %     subplot(2,4,tw-1);
                    %     plotperc(1:size(cpx,2),cpx,[66 33; 20 80; 49.9 50.1],'k',.15);
                    %     plot(1:size(cpx,2),nanmean(cpx),'k','linewidth',1.5);
                    %     plothline(0,'k:'); ylim([-3 3]/2*1e4);
                    %
                    %     subplot(2,4,4+tw-1);
                    %     plotperc(1:size(cpx,2),[s1cpx1; s2cpx1],[66 33; 20 80],'b',.1);
                    %     plotperc(1:size(cpx,2),-[s1cpx2; s2cpx2],[66 33; 20 80],'r',.1);
                    %     plot(1:size(cpx,2),nanmean(-[s1cpx2; s2cpx2]),'r','linewidth',1.5);
                    %     plot(1:size(cpx,2),nanmean([s1cpx1; s2cpx1]),'b','linewidth',1.5);
                    %     plothline(0,'k:'); ylim([-3 3]/2*1e4);
                    %     if tw==2 || tw==4; plotvline(200,'k:'); elseif tw==3 || tw==5; plotvline(size(cpx,2)-200,'k:'); end
                    %
                    
                    %subplot(2,4,8+tw-1);
                    %plotperc(1:size(cpx,2),[s1cpx2; s2cpx2],[49 51; 66 33; 20 80],'k',.15);
                    %ylim([-3 3]*1e4);
                    
                    %s1cpn=(isnan(s1cpx) | isnan(s1cpy)); s1cp0=(s1cpx==0 & s1cpy==0); s1cpx(s1cpn | s1cp0)=nan; s1cpy(s1cpn | s1cp0)=nan;
                    s1cpn=(isnan(s1cpx)); s1cp0=(s1cpx==0); s1cpx(s1cpn | s1cp0)=nan;
                    s1sum_L=sum(~isnan(rmnansandthreshold(-s1cpx,0)),2);
                    s1sum_R=sum(~isnan(rmnansandthreshold(+s1cpx,0)),2);
                    s1frac_R=s1sum_R./(s1sum_R+s1sum_L+eps);
                    
                    %s2cpn=(isnan(s2cpx) | isnan(s2cpy)); s2cp0=(s2cpx==0 & s2cpy==0); s2cpx(s2cpn | s2cp0)=nan; s2cpy(s2cpn | s2cp0)=nan;
                    s2cpn=(isnan(s2cpx)); s2cp0=(s2cpx==0); s2cpx(s2cpn | s2cp0)=nan;
                    s2sum_L=sum(~isnan(rmnansandthreshold(-s2cpx,0)),2);
                    s2sum_R=sum(~isnan(rmnansandthreshold(+s2cpx,0)),2);
                    s2frac_R=s2sum_R./(s2sum_R+s2sum_L+eps);
                    
                    %frac_R=[s1frac_R; s2frac_R];
                    frac_R=[mat2gray(s1frac_R); mat2gray(s2frac_R)];
                    %frac_L=sum_L./(sum_L+sum_R+eps); frac_R=1-frac_L;
                    %frac_R=sum_R./(sum_L+sum_R+eps);
                    %sum_R=mat2gray(sum_R); sum_L=mat2gray(sum_L);
                    %frac_R=sum_R./(sum_R+sum_L); frac_L=sum_L./(sum_L+sum_R);
                    
                    fracR_off1_del1_off2_del2(:,tw-1)=frac_R;
                    %fracR_off1_del1_off2_del2(:,tw-1)=nanmean(double([s1cpx; s2cpx]>0),2)./(nanmean(double([s1cpx; s2cpx]>0),2)+nanmean(double([s1cpx; s2cpx]<0),2)+eps);
                end
                
                evL   =[s1gbi.offerLev(s1gbi.ifLsR);  s1gbi.offerLev(s1gbi.ifRsL);  s2gbi.offerLev(s2gbi.ifLsR);  s2gbi.offerLev(s2gbi.ifRsL)];
                evR   =[s1gbi.offerRev(s1gbi.ifLsR);  s1gbi.offerRev(s1gbi.ifRsL);  s2gbi.offerRev(s2gbi.ifLsR);  s2gbi.offerRev(s2gbi.ifRsL)];
                varL  =[s1gbi.offerLvar(s1gbi.ifLsR); s1gbi.offerLvar(s1gbi.ifRsL); s2gbi.offerLvar(s2gbi.ifLsR); s2gbi.offerLvar(s2gbi.ifRsL)];
                varR  =[s1gbi.offerRvar(s1gbi.ifLsR); s1gbi.offerRvar(s1gbi.ifRsL); s2gbi.offerRvar(s2gbi.ifLsR); s2gbi.offerRvar(s2gbi.ifRsL)];
                sideRL=[ones(length(s1gbi.ifLsR),1); -ones(length(s1gbi.ifRsL),1);  ones(length(s2gbi.ifLsR),1); -ones(length(s2gbi.ifRsL),1)];
                s1chsd1=gbi1.choose12(gbi1.ifLsR)==2; s1chsd2=gbi1.choose12(gbi1.ifRsL)==1;
                s2chsd1=gbi2.choose12(gbi2.ifLsR)==2; s2chsd2=gbi2.choose12(gbi2.ifRsL)==1;
                
                if inclsideLR
                    x1=[evL/max(evL) evR/max(evR) varL/max(varL) varR/max(varR) sideRL fracR_off1_del1_off2_del2];
                else
                    x1=[evL/max(evL) evR/max(evR) varL/max(varL) varR/max(varR) fracR_off1_del1_off2_del2];
                end
                %x1=fracR_off1_del1_off2_del2;
                y1=[s1chsd1; s1chsd2; s2chsd1; s2chsd2];
                
                logitMdl= fitglm(x1,y1,'Distribution','binomial','Link','logit');
                
                subplot(2,2,focusaroundstimuli+2*inclsideLR+1)
                bar(1:size(x1,2),logitMdl.Coefficients.Estimate(2:end));
                for jj=1:size(x1,2); text(jj-.15,7,getpstars(logitMdl.Coefficients.pValue(jj+1))); end
                ylim([-1 1].*8); ylabel('logistic model weights'); box off
                if inclsideLR
                    gch=gca; gch.XTickLabel={'EV_L','EV_R','\sigma^2_L','\sigma^2_R','s_L_R','f_R offer1','f_R delay1','f_R offer2','f_R delay2'};
                else
                    gch=gca; gch.XTickLabel={'EV_L','EV_R','\sigma^2_L','\sigma^2_R','f_R offer1','f_R delay1','f_R offer2','f_R delay2'};
                end
                
                rectangle('Position',[4.7+inclsideLR 5 .6 1.5],'facecolor','k');
                rectangle('Position',[5.7+inclsideLR 5 .6 1.5],'facecolor','k');
                rectangle('Position',[6.7+inclsideLR 5 .6 1.5],'facecolor','k');
                rectangle('Position',[7.7+inclsideLR 5 .6 1.5],'facecolor','k');
                
                rectangle('Position',[4.8+inclsideLR 5.2 .1 1],'facecolor',[0 .7 0]);  rectangle('Position',[4.8+inclsideLR 6 .1 .2],'facecolor',[.7 0 0]);
                rectangle('Position',[7.1+inclsideLR 5.2 .1 1],'facecolor',[0 0 .7]);  rectangle('Position',[7.1+inclsideLR 5.7 .1 .5],'facecolor',[.7 0 0]);
                %gch=gca; gch.XTickLabel={'f_R offer1','f_R delay1','f_R offer2','f_R delay2'};
                gch.XTickLabelRotation=30;
                if ~focusaroundstimuli && inclsideLR
                    title('logit model of choice - include side s_L_R')
                elseif ~focusaroundstimuli && ~inclsideLR
                    title('logit model of choice')
                elseif focusaroundstimuli && inclsideLR
                    title('logit model of choice - include side s_L_R - focus around stimuli')
                elseif focusaroundstimuli && ~inclsideLR
                    title('logit model of choice - focus around stimuli')
                end
            end
        end
    end
  supertitle('All trials');
  saveas(hfig,[figures_folder '/Fig2D-Logit_off1del1off2del2_ALL.png']);
  saveas(hfig,[figures_folder '/Fig2D-Logit_off1del1off2del2_ALL.fig']);
  saveas(hfig,[figures_folder '/Fig2D-Logit_off1del1off2del2_ALL.svg']);
    
    
end