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

addpath(genpath('../../../../my-tools'));
addpath(genpath('../../../../external-tools'));

clearvars -except subjsData poolstack* pooltw*

rng(0)
curr_rng=rng; % saves random seed as the results may vary

recomputelm=1;
Nlag=10; Nspan=200; 
Nshf=0; selectdifftrials=1;

if ~exist('poolstackvars','var')
    load '../../poolstackvars_270923.mat';
end
if ~exist('pooltweyepos','var')
    load '../latestpooltweyepos_14feb22.mat'
end
if ~exist('pooltwspksq','var')
    load '../latestpooltwspks_14feb22.mat'
end

twLabels={'preoffer1','offer1','delay1','offer2','delay2','refix','ch-go','ch-hold'};

nccs=[ 51 59 24 29 18 32  9 26];
ntrs=arrayfun(@(xx) size(poolstackvars(xx).offer1sd,1),1:8);

Nccs=sum(nccs);
Ntrs12=max(ntrs);
Ntps12=arrayfun(@(tw) size(pooltweyepos(tw).posX,3), 1:8);
Ntps12(1)=400;

pooleposxmat=cat(3,pooltweyepos(1).posX,pooltweyepos(2).posX,pooltweyepos(3).posX,pooltweyepos(4).posX,...
    pooltweyepos(5).posX,pooltweyepos(6).posX,pooltweyepos(7).posX,pooltweyepos(8).posX);
poolspksqmat=cat(3,pooltwspksq(1).spkseqs,pooltwspksq(2).spkseqs,pooltwspksq(3).spkseqs,pooltwspksq(4).spkseqs,...
    pooltwspksq(5).spkseqs,pooltwspksq(6).spkseqs,pooltwspksq(7).spkseqs,pooltwspksq(8).spkseqs);
poolcovhist=cat(2,pooltwspksq(1).covhist,pooltwspksq(2).covhist,pooltwspksq(3).covhist,pooltwspksq(4).covhist,...
    pooltwspksq(5).covhist,pooltwspksq(6).covhist,pooltwspksq(7).covhist,pooltwspksq(8).covhist);

clearvars pooltweyepos pooltwspksq


pooloffer1sd=[]; poolofferLev=[]; poolofferRev=[]; pooloffer1ev=[]; pooloffer2ev=[]; pooloffer1var=[]; pooloffer2var=[]; 
poolofferchR=[]; pooloffer1pr=[]; pooloffer2pr=[]; pooloffer1rw=[]; pooloffer2rw=[];
for sn=1:8
    pooloffer1sd=cat(2, pooloffer1sd, repmat(cat(1,poolstackvars(sn).offer1sd,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    poolofferLev=cat(2, poolofferLev, repmat(cat(1,poolstackvars(sn).offerLev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    poolofferRev=cat(2, poolofferRev, repmat(cat(1,poolstackvars(sn).offerRev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer1ev=cat(2, pooloffer1ev, repmat(cat(1,poolstackvars(sn).offer1ev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer2ev=cat(2, pooloffer2ev, repmat(cat(1,poolstackvars(sn).offer2ev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer1var=cat(2, pooloffer1var, repmat(cat(1,poolstackvars(sn).offer1var,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer2var=cat(2, pooloffer2var, repmat(cat(1,poolstackvars(sn).offer2var,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    poolofferchR=cat(2, poolofferchR, repmat(cat(1,poolstackvars(sn).chooseR, nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer1pr=cat(2, pooloffer1pr, repmat(cat(1,poolstackvars(sn).offer1pr,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer2pr=cat(2, pooloffer2pr, repmat(cat(1,poolstackvars(sn).offer2pr,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer1rw=cat(2, pooloffer1rw, repmat(cat(1,poolstackvars(sn).offer1rw,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer2rw=cat(2, pooloffer2rw, repmat(cat(1,poolstackvars(sn).offer2rw,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
end

%clearvars pooltwspksq pooltweyepos

%fchg=find(poolcovhist>=.8); fchl=find(poolcovhist<.8);
fchg=find(poolcovhist>=.999); fchl=find(poolcovhist<.999);
ntcv=[fchg(diff([fchg sum(Ntps12)])>1)' fchl([diff(fchl) sum(Ntps12)]>1)'];
ntce=[[1; 401; ntcv(1:end-1,2)] [400; ntcv(:,1)]];
vtps=ones(1,sum(Ntps12)); for jj=1:7; vtps(ntcv(jj,1):ntcv(jj,2))=nan; end

%% smooth by movmean filtering, downsample to 20ms, 50Hz resolution
Ntps_ds=ceil(sum(Ntps12)/Nlag);
%tmax_ds=nanmovmean(1:size(pooleposxmat,3),Nlag,Nlag);
tmax_ds=nanmovmean(1:sum(Ntps12),Nlag,Nlag);
poolspksqmat_ds=nan(Ntrs12,Nccs,length(tmax_ds));
pooleposxmat_ds=nan(Ntrs12,Nccs,length(tmax_ds));

for cc=1:Nccs
    poolspksqmat_ds(:,cc,:)=nanmovmean(squeeze(poolspksqmat(:,cc,:))',Nlag,[0 Nspan])';
    pooleposxmat_ds(:,cc,:)=nanmovmean(squeeze(pooleposxmat(:,cc,:))',Nlag,Nlag)';
end

clearvars poolspksqmat pooleposxmat

poolcovhist_ds=nanmovmean(poolcovhist,Nlag,Nlag);
%fcgh_ds=find(poolcovhist_ds>=.8); fchl_ds=find(poolcovhist_ds<.8);
fcgh_ds=find(poolcovhist_ds>=.999); fchl_ds=find(poolcovhist_ds<.999);
ntcv_ds=[fcgh_ds(diff([fcgh_ds Ntps_ds])>1)' fchl_ds([diff(fchl_ds) Ntps_ds]>1)'];
vtps_ds=ones(1,Ntps_ds);  for jj=1:7; vtps_ds(ntcv_ds(jj,1):ntcv_ds(jj,2))=nan; end
ntce_ds=[[1 ceil(400/Nlag)+1 find(vtps_ds(2:end)+isnan(vtps_ds(1:end-1))==2)+1]' [ceil(400/Nlag) find(vtps_ds(1:end-1)+isnan(vtps_ds(2:end))==2)]'];


vtmax_ds=rmnans(tmax_ds.*vtps_ds);%find(rmnans(vtps_ds));
rtmax_ds=1:length(vtmax_ds);
rtmes_ds=arrayfun(@(jj) length(ntce_ds(jj,1):ntce_ds(jj,2)), 1:8);


%No baseline normalization
poolspksqmat_ds_bln=poolspksqmat_ds;
clearvars poolspksqmat_ds

%%
del2tpoi=ntce_ds(5,1):ntce_ds(5,2);
Ntps_ds=length(del2tpoi);

%beh_weights_mat=nan(Ntps_ds,Nccs,7);
%beh_pvalues_mat=nan(Ntps_ds,Nccs,7);
%zetas_ss1_mat=nan(Ntrs12,Nccs,Ntps_ds);     zetas_ss2_mat=nan(Ntrs12,Nccs,Ntps_ds);
%zetasLL_ss1_mat=nan(Ntrs12,Nccs,Ntps_ds);   zetasLL_ss2_mat=nan(Ntrs12,Nccs,Ntps_ds);
%zetasLR_ss1_mat=nan(Ntrs12,Nccs,Ntps_ds);   zetasLR_ss2_mat=nan(Ntrs12,Nccs,Ntps_ds);

%zetas_ch_bw1_ss1LL_mat=nan(Ntrs12,Nccs,Ntps_ds);  zetas_ch_bw1_ss1LR_mat=nan(Ntrs12,Nccs,Ntps_ds); 
%zetas_ch_bw1_ss2LL_mat=nan(Ntrs12,Nccs,Ntps_ds);  zetas_ch_bw1_ss2LR_mat=nan(Ntrs12,Nccs,Ntps_ds); 
%zetas_ch_bw2_ss1LL_mat=nan(Ntrs12,Nccs,Ntps_ds);  zetas_ch_bw2_ss1LR_mat=nan(Ntrs12,Nccs,Ntps_ds); 
%zetas_ch_bw2_ss2LL_mat=nan(Ntrs12,Nccs,Ntps_ds);  zetas_ch_bw2_ss2LR_mat=nan(Ntrs12,Nccs,Ntps_ds); 

Nsubg=1;

zetas_noch_bw0_ss1LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw0_ss1LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 
zetas_noch_bw0_ss2LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw0_ss2LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 
zetas_noch_bw1_ss1LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw1_ss1LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 
zetas_noch_bw1_ss2LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw1_ss2LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 
zetas_noch_bw2_ss1LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw2_ss1LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 
zetas_noch_bw2_ss2LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw2_ss2LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 

zetas_noch_bw0_ss0LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);  zetas_noch_bw0_ss0LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw0_ss0LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw0_ss0LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 
pvals_noch_bw0_ss0LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  pvals_noch_bw0_ss0LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 

betas_ch_bw0_ss1LL_mat=nan(4,Nccs,Ntps_ds,Nsubg);  betas_ch_bw0_ss1LR_mat=nan(4,Nccs,Ntps_ds,Nsubg); 
betas_ch_bw0_ss2LL_mat=nan(4,Nccs,Ntps_ds,Nsubg);  betas_ch_bw0_ss2LR_mat=nan(4,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw0_ss1LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw0_ss1LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw0_ss2LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw0_ss2LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 
betas_ch_bw1_ss1LL_mat=nan(4,Nccs,Ntps_ds,Nsubg);  betas_ch_bw1_ss1LR_mat=nan(4,Nccs,Ntps_ds,Nsubg); 
betas_ch_bw1_ss2LL_mat=nan(4,Nccs,Ntps_ds,Nsubg);  betas_ch_bw1_ss2LR_mat=nan(4,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw1_ss1LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw1_ss1LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw1_ss2LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw1_ss2LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 
betas_ch_bw2_ss1LL_mat=nan(4,Nccs,Ntps_ds,Nsubg);  betas_ch_bw2_ss1LR_mat=nan(4,Nccs,Ntps_ds,Nsubg); 
betas_ch_bw2_ss2LL_mat=nan(4,Nccs,Ntps_ds,Nsubg);  betas_ch_bw2_ss2LR_mat=nan(4,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw2_ss1LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw2_ss1LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 
betas_noch_bw2_ss2LL_mat=nan(3,Nccs,Ntps_ds,Nsubg);  betas_noch_bw2_ss2LR_mat=nan(3,Nccs,Ntps_ds,Nsubg); 

useinterceptflag=0;
bbetas_ch_bw0_ss1LL_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg);  bbetas_ch_bw0_ss1LR_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg); 
bbetas_ch_bw0_ss2LL_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg);  bbetas_ch_bw0_ss2LR_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg); 
bbetas_ch_bw1_ss1LL_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg);  bbetas_ch_bw1_ss1LR_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg); 
bbetas_ch_bw1_ss2LL_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg);  bbetas_ch_bw1_ss2LR_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg); 
bbetas_ch_bw2_ss1LL_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg);  bbetas_ch_bw2_ss1LR_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg); 
bbetas_ch_bw2_ss2LL_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg);  bbetas_ch_bw2_ss2LR_mat=nan(1+useinterceptflag,Nccs,Ntps_ds,Nsubg); 

pvals_ch_bw0_ss1LL_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg));  pvals_ch_bw0_ss1LR_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg)); 
pvals_ch_bw0_ss2LL_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg));  pvals_ch_bw0_ss2LR_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg)); 
pvals_noch_bw0_ss1LL_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg));  pvals_noch_bw0_ss1LR_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg)); 
pvals_noch_bw0_ss2LL_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg));  pvals_noch_bw0_ss2LR_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg)); 
pvals_ch_bw1_ss1LL_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg));  pvals_ch_bw1_ss1LR_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg)); 
pvals_ch_bw1_ss2LL_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg));  pvals_ch_bw1_ss2LR_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg)); 
pvals_noch_bw1_ss1LL_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg));  pvals_noch_bw1_ss1LR_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg)); 
pvals_noch_bw1_ss2LL_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg));  pvals_noch_bw1_ss2LR_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg)); 
pvals_ch_bw2_ss1LL_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg));  pvals_ch_bw2_ss1LR_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg)); 
pvals_ch_bw2_ss2LL_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg));  pvals_ch_bw2_ss2LR_mat=uint8(zeros(4,Nccs,Ntps_ds,Nsubg)); 
pvals_noch_bw2_ss1LL_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg));  pvals_noch_bw2_ss1LR_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg)); 
pvals_noch_bw2_ss2LL_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg));  pvals_noch_bw2_ss2LR_mat=uint8(zeros(3,Nccs,Ntps_ds,Nsubg)); 

ppvals_ch_bw0_ss1LL_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg));  ppvals_ch_bw0_ss1LR_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg)); 
ppvals_ch_bw0_ss2LL_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg));  ppvals_ch_bw0_ss2LR_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg)); 
ppvals_ch_bw1_ss1LL_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg));  ppvals_ch_bw1_ss1LR_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg)); 
ppvals_ch_bw1_ss2LL_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg));  ppvals_ch_bw1_ss2LR_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg)); 
ppvals_ch_bw2_ss1LL_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg));  ppvals_ch_bw2_ss1LR_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg)); 
ppvals_ch_bw2_ss2LL_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg));  ppvals_ch_bw2_ss2LR_mat=uint8(zeros(1+useinterceptflag,Nccs,Ntps_ds,Nsubg)); 

offch_ss0LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   offch_ss0LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);
offch_ss1LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   offch_ss1LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   
offch_ss2LL_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   offch_ss2LR_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);
offch_ss0LL_noflip_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   offch_ss0LR_noflip_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);
offch_ss1LL_noflip_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   offch_ss1LR_noflip_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   
offch_ss2LL_noflip_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);   offch_ss2LR_noflip_mat=nan(Ntrs12,Nccs,Ntps_ds,Nsubg);


% beta0s_ss1_mat=nan(Nccs,Ntps_ds);           beta0s_ss2_mat=nan(Nccs,Ntps_ds);
% beta1s_ss1_mat=nan(Nccs,Ntps_ds);           beta1s_ss2_mat=nan(Nccs,Ntps_ds);
% beta2s_ss1_mat=nan(Nccs,Ntps_ds);           beta2s_ss2_mat=nan(Nccs,Ntps_ds);
% beta0sLL_ss1_mat=nan(Nccs,Ntps_ds);         beta0sLL_ss2_mat=nan(Nccs,Ntps_ds);
% beta1sLL_ss1_mat=nan(Nccs,Ntps_ds);         beta1sLL_ss2_mat=nan(Nccs,Ntps_ds);
% beta2sLL_ss1_mat=nan(Nccs,Ntps_ds);         beta2sLL_ss2_mat=nan(Nccs,Ntps_ds);
% beta0sLR_ss1_mat=nan(Nccs,Ntps_ds);         beta0sLR_ss2_mat=nan(Nccs,Ntps_ds);
% beta1sLR_ss1_mat=nan(Nccs,Ntps_ds);         beta1sLR_ss2_mat=nan(Nccs,Ntps_ds);
% beta2sLR_ss1_mat=nan(Nccs,Ntps_ds);         beta2sLR_ss2_mat=nan(Nccs,Ntps_ds);
% 
% beta00s_ss1_mat=nan(Nccs,Ntps_ds);          beta00s_ss2_mat=nan(Nccs,Ntps_ds);
% beta01s_ss1_mat=nan(Nccs,Ntps_ds);          beta01s_ss2_mat=nan(Nccs,Ntps_ds);
% beta00sLL_ss1_mat=nan(Nccs,Ntps_ds);        beta00sLL_ss2_mat=nan(Nccs,Ntps_ds);
% beta01sLL_ss1_mat=nan(Nccs,Ntps_ds);        beta01sLL_ss2_mat=nan(Nccs,Ntps_ds);
% beta00sLR_ss1_mat=nan(Nccs,Ntps_ds);        beta00sLR_ss2_mat=nan(Nccs,Ntps_ds);
% beta01sLR_ss1_mat=nan(Nccs,Ntps_ds);        beta01sLR_ss2_mat=nan(Nccs,Ntps_ds);

if recomputelm
for sg=1:Nsubg
for cc=1:Nccs
    %ccspk=squeeze(poolspksqmat_ds(:,cc,:));
    ccspk=squeeze(poolspksqmat_ds_bln(:,cc,:));
    ccepx=squeeze(pooleposxmat_ds(:,cc,:));
    ccoffEVL=squeeze(pooloffer1ev(:,cc));
    ccoffEVR=squeeze(pooloffer2ev(:,cc));
    ccoffPrL=squeeze(pooloffer1pr(:,cc));
    ccoffPrR=squeeze(pooloffer2pr(:,cc));
    ccoffRwL=squeeze(pooloffer1rw(:,cc));
    ccoffRwR=squeeze(pooloffer2rw(:,cc));
    ccoffvarL=squeeze(pooloffer1var(:,cc));
    ccoffvarR=squeeze(pooloffer2var(:,cc));
    ccoffchR=squeeze(poolofferchR(:,cc));
    ccoff1sd=squeeze(pooloffer1sd(:,cc));

    toskip=isnan(ccoffchR);
    if selectdifftrials
        toskip=isnan(ccoffchR)+double(abs(ccoffEVL-ccoffEVR)>1);
    end

    vccoff1sd=ccoff1sd(~toskip);  % 1stL=1, 1stR=0
    vccoffchR=ccoffchR(~toskip);
    vccoffEVR=ccoffEVR(~toskip);
    vccoffEVL=ccoffEVL(~toskip);
    %max(abs(vccoffEVL-vccoffEVR));
    %vccoffPrR=ccoffPrR(~toskip);
    %vccoffPrL=ccoffPrL(~toskip);
    vccoffRwR=ccoffRwR(~toskip);
    vccoffRwL=ccoffRwL(~toskip);
    vccoffvarR=ccoffvarR(~toskip);
    vccoffvarL=ccoffvarL(~toskip);

    % FOR SVs WE USE chR as 0,1
    vccoffchR_flip=vccoffchR;
    vccoffchR_flip(vccoff1sd==0)=1-vccoffchR_flip(vccoff1sd==0);
    %vccoffchR_flip(vccoffchR_flip==0)=-1;
    % FOR THE LINEAR MODEL WE USE chR AS -1 +1
    vccoffchR(vccoffchR==0)=-1;

    nicc=length(vccoffchR);
    ss1=randperm(nicc,round(nicc/2));
    ss2=setdiff(1:nicc,ss1);
    
    behvars=[-vccoffRwL -vccoffEVL -vccoffvarL vccoffRwR vccoffEVR vccoffvarR];
    behchglm_bw0=fitglm(behvars,vccoffchR_flip,'Distribution','binomial','link','logit');
    behchglm_weights_bw0=behchglm_bw0.Coefficients.Estimate(2:end);
    behchglm_bw1=fitglm(behvars(ss1,:),vccoffchR_flip(ss1),'Distribution','binomial','link','logit');
    behchglm_weights_bw1=behchglm_bw1.Coefficients.Estimate(2:end);
    behchglm_bw2=fitglm(behvars(ss2,:),vccoffchR_flip(ss2),'Distribution','binomial','link','logit');
    behchglm_weights_bw2=behchglm_bw2.Coefficients.Estimate(2:end);

    vccoffSVL_bw0=behchglm_weights_bw0(1).*vccoffRwL+behchglm_weights_bw0(2).*vccoffEVL+behchglm_weights_bw0(3).*vccoffvarL;
    vccoffSVR_bw0=behchglm_weights_bw0(4).*vccoffRwR+behchglm_weights_bw0(5).*vccoffEVR+behchglm_weights_bw0(6).*vccoffvarR;
    vccoffSVL_bw1=behchglm_weights_bw1(1).*vccoffRwL+behchglm_weights_bw1(2).*vccoffEVL+behchglm_weights_bw1(3).*vccoffvarL;
    vccoffSVR_bw1=behchglm_weights_bw1(4).*vccoffRwR+behchglm_weights_bw1(5).*vccoffEVR+behchglm_weights_bw1(6).*vccoffvarR;
    vccoffSVL_bw2=behchglm_weights_bw2(1).*vccoffRwL+behchglm_weights_bw2(2).*vccoffEVL+behchglm_weights_bw2(3).*vccoffvarL;
    vccoffSVR_bw2=behchglm_weights_bw2(4).*vccoffRwR+behchglm_weights_bw2(5).*vccoffEVR+behchglm_weights_bw2(6).*vccoffvarR;

    parfor tt=1:Ntps_ds
        %warning('off','all')
        if vtps_ds(del2tpoi(tt))==1
            
             ctccspk=squeeze(ccspk(:,del2tpoi(tt)));
             ctccepx=squeeze(ccepx(:,del2tpoi(tt)));
            %toskip=isnan(ccoffchR); % from outer loop
            %toskip=isnan(ccoffchR) | isnan(ctccspk) | isnan(ctccepx);
             vctccspk=ctccspk(~toskip);
             vctccepx=ctccepx(~toskip);
             ss0LL=find(vctccepx<0); niLL=length(ss0LL);
             ss0LR=find(vctccepx>0); niLR=length(ss0LR);
            % one randomization for all bins
             ss1LL=intersect(ss0LL,ss1);%randperm(niLL,round(niLL/2));
             ss1LR=intersect(ss0LR,ss1);%randperm(niLR,round(niLR/2));
             ss2LL=intersect(ss0LL,ss2);%setdiff(1:niLL,ss1LL);
             ss2LR=intersect(ss0LR,ss2);%setdiff(1:niLR,ss1LR);
            % bin-by-bin randomization
            %ss1LL=randperm(niLL,round(niLL/2));
            %ss1LR=randperm(niLR,round(niLR/2));
            %ss2LL=setdiff(1:niLL,ss1LL);
            %ss2LR=setdiff(1:niLR,ss1LR);

            % FITTING LM FOR SV: NOCH LL/LR (SS0, BW0)
            ctwcclm_noch_bw0_ss0LL=fitlm([vccoffSVL_bw0(ss0LL) vccoffSVR_bw0(ss0LL)], vctccspk(ss0LL));
            zetas_noch_bw0_ss0LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw0_ss0LL.Residuals.Raw; nan(Ntrs12-length(ss0LL),1)];
            betas_noch_bw0_ss0LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw0_ss0LL.Coefficients.Estimate;
            pvals_noch_bw0_ss0LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw0_ss0LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw0_ss0LR=fitlm([vccoffSVL_bw0(ss0LR) vccoffSVR_bw0(ss0LR)], vctccspk(ss0LR));
            zetas_noch_bw0_ss0LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw0_ss0LR.Residuals.Raw; nan(Ntrs12-length(ss0LR),1)];
            betas_noch_bw0_ss0LR_mat(:,cc,tt,sg)=ctwcclm_noch_bw0_ss0LR.Coefficients.Estimate;
            pvals_noch_bw0_ss0LR_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw0_ss0LR.Coefficients.pValue<0.05);


            % %% Quick check about the sign of residuals
            % eta=vctccspk;
            % svL=vccoffSVL_bw0;
            % svR=vccoffSVR_bw0;
            % lmL=fitlm(svL,eta);
            % lmR=fitlm(svL,eta);
            % zetaL=lmL.Residuals.Raw;
            % betaL=lmL.Coefficients.Estimate;
            % zetaR=lmR.Residuals.Raw;
            % betaR=lmR.Coefficients.Estimate;
            % 
            % fullfig([0 0 .4 .4]);
            % subplot(2,2,1);
            % plot(svL,eta,'k.'); hold on
            % plot(sort(svL),beta(1)+sort(svL).*beta(2),'-');
            % xlabel('SV_L'); ylabel('\eta');
            % subplot(2,2,2);
            % plot(svR,eta,'k.'); hold on
            % plot(sort(svR),beta(1)+sort(svR).*beta(2),'-');
            % xlabel('SV_R'); ylabel('\eta');
            % 
            % subplot(2,2,3);
            % plot(svL,zeta,'r.'); hold on
            % xlabel('SV_L'); ylabel('residuals');
            % subplot(2,2,4);
            % plot(svR,zeta,'r.'); hold on
            % xlabel('SV_R'); ylabel('residuals');
            %
            % %%

            % FITTING LM FOR SV: CH/NOCH LL/LR (SS1, BW0) 
            ctwcclm_ch_bw0_ss1LL=fitlm([vccoffSVL_bw0(ss1LL) vccoffSVR_bw0(ss1LL) vccoffchR(ss1LL)], vctccspk(ss1LL));
            %zetas_ch_bw0_ss1LL_mat(:,cc,tt,sg)=[ctwcclm_ch_bw0_ss1LL.Residuals.Raw; nan(Ntrs12-length(ss1LL),1)];
            betas_ch_bw0_ss1LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss1LL.Coefficients.Estimate;
            pvals_ch_bw0_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss1LL.Coefficients.pValue<0.05);
            ctwcclm_ch_bw0_ss1LR=fitlm([vccoffSVL_bw0(ss1LR) vccoffSVR_bw0(ss1LR) vccoffchR(ss1LR)], vctccspk(ss1LR));
            %zetas_ch_bw0_ss1LR_mat(:,cc,tt,sg)=[ctwcclm_ch_bw0_ss1LR.Residuals.Raw; nan(Ntrs12-length(ss1LR),1)];
            betas_ch_bw0_ss1LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss1LR.Coefficients.Estimate;
            pvals_ch_bw0_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss1LR.Coefficients.pValue<0.05);

            ctwcclm_noch_bw0_ss1LL=fitlm([vccoffSVL_bw0(ss1LL) vccoffSVR_bw0(ss1LL)], vctccspk(ss1LL));
            zetas_noch_bw0_ss1LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw0_ss1LL.Residuals.Raw; nan(Ntrs12-length(ss1LL),1)];
            betas_noch_bw0_ss1LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw0_ss1LL.Coefficients.Estimate;
            pvals_noch_bw0_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw0_ss1LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw0_ss1LR=fitlm([vccoffSVL_bw0(ss1LR) vccoffSVR_bw0(ss1LR)], vctccspk(ss1LR));
            zetas_noch_bw0_ss1LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw0_ss1LR.Residuals.Raw; nan(Ntrs12-length(ss1LR),1)];
            betas_noch_bw0_ss1LR_mat(:,cc,tt,sg)=ctwcclm_noch_bw0_ss1LR.Coefficients.Estimate;
            pvals_noch_bw0_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw0_ss1LR.Coefficients.pValue<0.05);

            ctwcclm_ch_bw0_ss1LL_res=fitlm(vccoffchR(ss1LL),ctwcclm_noch_bw0_ss1LL.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw0_ss1LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss1LL_res.Coefficients.Estimate;
            ppvals_ch_bw0_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss1LL_res.Coefficients.pValue<0.05);
            ctwcclm_ch_bw0_ss1LR_res=fitlm(vccoffchR(ss1LR),ctwcclm_noch_bw0_ss1LR.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw0_ss1LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss1LR_res.Coefficients.Estimate;
            ppvals_ch_bw0_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss1LR_res.Coefficients.pValue<0.05);

            % FITTING LM FOR SV: CH/NOCH LL/LR (SS2, BW0)
            ctwcclm_ch_bw0_ss2LL=fitlm([vccoffSVL_bw0(ss2LL) vccoffSVR_bw0(ss2LL) vccoffchR(ss2LL)], vctccspk(ss2LL));
            %zetas_ch_bw0_ss2LL_mat(:,cc,tt,sg)=[ctwcclm_ch_bw0_ss2LL.Residuals.Raw; nan(Ntrs12-length(ss2LL),1)];
            betas_ch_bw0_ss2LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss2LL.Coefficients.Estimate;
            pvals_ch_bw0_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss2LL.Coefficients.pValue<0.05);
            ctwcclm_ch_bw0_ss2LR=fitlm([vccoffSVL_bw0(ss2LR) vccoffSVR_bw0(ss2LR) vccoffchR(ss2LR)], vctccspk(ss2LR));
            %zetas_ch_bw0_ss2LR_mat(:,cc,tt,sg)=[ctwcclm_ch_bw0_ss2LR.Residuals.Raw; nan(Ntrs12-length(ss2LR),1)];
            betas_ch_bw0_ss2LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss2LR.Coefficients.Estimate;
            pvals_ch_bw0_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss2LR.Coefficients.pValue<0.05);

            ctwcclm_noch_bw0_ss2LL=fitlm([vccoffSVL_bw0(ss2LL) vccoffSVR_bw0(ss2LL)], vctccspk(ss2LL));
            zetas_noch_bw0_ss2LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw0_ss2LL.Residuals.Raw; nan(Ntrs12-length(ss2LL),1)];
            betas_noch_bw0_ss2LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw0_ss2LL.Coefficients.Estimate;
            pvals_noch_bw0_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw0_ss2LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw0_ss2LR=fitlm([vccoffSVL_bw0(ss2LR) vccoffSVR_bw0(ss2LR)], vctccspk(ss2LR));
            zetas_noch_bw0_ss2LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw0_ss2LR.Residuals.Raw; nan(Ntrs12-length(ss2LR),1)];
            pvals_noch_bw0_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw0_ss2LR.Coefficients.pValue<0.05);

            ctwcclm_ch_bw0_ss2LL_res=fitlm(vccoffchR(ss2LL),ctwcclm_noch_bw0_ss2LL.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw0_ss2LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss2LL_res.Coefficients.Estimate;
            ppvals_ch_bw0_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss2LL_res.Coefficients.pValue<0.05);
            ctwcclm_ch_bw0_ss2LR_res=fitlm(vccoffchR(ss2LR),ctwcclm_noch_bw0_ss2LR.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw0_ss2LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw0_ss2LR_res.Coefficients.Estimate;
            ppvals_ch_bw0_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw0_ss2LR_res.Coefficients.pValue<0.05);


            % FITTING LM FOR SV: CH/NOCH LL/LR (SS1, BW1)
            ctwcclm_ch_bw1_ss1LL=fitlm([vccoffSVL_bw1(ss1LL) vccoffSVR_bw1(ss1LL) vccoffchR(ss1LL)], vctccspk(ss1LL));
            %zetas_ch_bw1_ss1LL_mat(:,cc,tt,sg)=[ctwcclm_ch_bw1_ss1LL.Residuals.Raw; nan(Ntrs12-length(ss1LL),1)];
            betas_ch_bw1_ss1LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss1LL.Coefficients.Estimate;
            pvals_ch_bw1_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss1LL.Coefficients.pValue<0.05);
            ctwcclm_ch_bw1_ss1LR=fitlm([vccoffSVL_bw1(ss1LR) vccoffSVR_bw1(ss1LR) vccoffchR(ss1LR)], vctccspk(ss1LR));
            %zetas_ch_bw1_ss1LR_mat(:,cc,tt,sg)=[ctwcclm_ch_bw1_ss1LR.Residuals.Raw; nan(Ntrs12-length(ss1LR),1)];
            betas_ch_bw1_ss1LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss1LR.Coefficients.Estimate;
            pvals_ch_bw1_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss1LR.Coefficients.pValue<0.05);

            ctwcclm_noch_bw1_ss1LL=fitlm([vccoffSVL_bw1(ss1LL) vccoffSVR_bw1(ss1LL)], vctccspk(ss1LL));
            zetas_noch_bw1_ss1LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw1_ss1LL.Residuals.Raw; nan(Ntrs12-length(ss1LL),1)];
            betas_noch_bw1_ss1LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw1_ss1LL.Coefficients.Estimate;
            pvals_noch_bw1_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw1_ss1LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw1_ss1LR=fitlm([vccoffSVL_bw1(ss1LR) vccoffSVR_bw1(ss1LR)], vctccspk(ss1LR));
            zetas_noch_bw1_ss1LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw1_ss1LR.Residuals.Raw; nan(Ntrs12-length(ss1LR),1)];
            pvals_noch_bw1_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw1_ss1LR.Coefficients.pValue<0.05);

            ctwcclm_ch_bw1_ss1LL_res=fitlm(vccoffchR(ss1LL),ctwcclm_noch_bw1_ss1LL.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw1_ss1LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss1LL_res.Coefficients.Estimate;
            ppvals_ch_bw1_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss1LL_res.Coefficients.pValue<0.05);
            ctwcclm_ch_bw1_ss1LR_res=fitlm(vccoffchR(ss1LR),ctwcclm_noch_bw1_ss1LR.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw1_ss1LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss1LR_res.Coefficients.Estimate;
            ppvals_ch_bw1_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss1LR_res.Coefficients.pValue<0.05);

            % FITTING LM FOR SV: CH/NOCH LL/LR (SS2, BW1)
            ctwcclm_ch_bw1_ss2LL=fitlm([vccoffSVL_bw1(ss2LL) vccoffSVR_bw1(ss2LL) vccoffchR(ss2LL)], vctccspk(ss2LL));
            %zetas_ch_bw1_ss2LL_mat(:,cc,tt,sg)=[ctwcclm_ch_bw1_ss2LL.Residuals.Raw; nan(Ntrs12-length(ss2LL),1)];
            betas_ch_bw1_ss2LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss2LL.Coefficients.Estimate;
            pvals_ch_bw1_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss2LL.Coefficients.pValue<0.05);
            ctwcclm_ch_bw1_ss2LR=fitlm([vccoffSVL_bw1(ss2LR) vccoffSVR_bw1(ss2LR) vccoffchR(ss2LR)], vctccspk(ss2LR));
            %zetas_ch_bw1_ss2LR_mat(:,cc,tt,sg)=[ctwcclm_ch_bw1_ss2LR.Residuals.Raw; nan(Ntrs12-length(ss2LR),1)];
            betas_ch_bw1_ss2LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss2LR.Coefficients.Estimate;
            pvals_ch_bw1_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss2LR.Coefficients.pValue<0.05);

            ctwcclm_noch_bw1_ss2LL=fitlm([vccoffSVL_bw1(ss2LL) vccoffSVR_bw1(ss2LL)], vctccspk(ss2LL));
            zetas_noch_bw1_ss2LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw1_ss2LL.Residuals.Raw; nan(Ntrs12-length(ss2LL),1)];
            betas_noch_bw1_ss2LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw1_ss2LL.Coefficients.Estimate;
            pvals_noch_bw1_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw1_ss2LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw1_ss2LR=fitlm([vccoffSVL_bw1(ss2LR) vccoffSVR_bw1(ss2LR)], vctccspk(ss2LR));
            zetas_noch_bw1_ss2LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw1_ss2LR.Residuals.Raw; nan(Ntrs12-length(ss2LR),1)];
            betas_noch_bw1_ss2LR_mat(:,cc,tt,sg)=ctwcclm_noch_bw1_ss2LR.Coefficients.Estimate;
            pvals_noch_bw1_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw1_ss2LR.Coefficients.pValue<0.05);

            ctwcclm_ch_bw1_ss2LL_res=fitlm(vccoffchR(ss2LL),ctwcclm_noch_bw1_ss2LL.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw1_ss2LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss2LL_res.Coefficients.Estimate;
            ppvals_ch_bw1_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss2LL_res.Coefficients.pValue<0.05);
            ctwcclm_ch_bw1_ss2LR_res=fitlm(vccoffchR(ss2LR),ctwcclm_noch_bw1_ss2LR.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw1_ss2LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw1_ss2LR_res.Coefficients.Estimate;
            ppvals_ch_bw1_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw1_ss2LR_res.Coefficients.pValue<0.05);

            % FITTING LM FOR SV: CH/NOCH LL/LR (SS1, BW2)
            ctwcclm_ch_bw2_ss1LL=fitlm([vccoffSVL_bw2(ss1LL) vccoffSVR_bw2(ss1LL) vccoffchR(ss1LL)], vctccspk(ss1LL));
            %zetas_ch_bw2_ss1LL_mat(:,cc,tt,sg)=[ctwcclm_ch_bw2_ss1LL.Residuals.Raw; nan(Ntrs12-length(ss1LL),1)];
            betas_ch_bw2_ss1LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss1LL.Coefficients.Estimate;
            pvals_ch_bw2_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss1LL.Coefficients.pValue<0.05);
            ctwcclm_ch_bw2_ss1LR=fitlm([vccoffSVL_bw2(ss1LR) vccoffSVR_bw2(ss1LR) vccoffchR(ss1LR)], vctccspk(ss1LR));
            %zetas_ch_bw2_ss1LR_mat(:,cc,tt,sg)=[ctwcclm_ch_bw2_ss1LR.Residuals.Raw; nan(Ntrs12-length(ss1LR),1)];
            betas_ch_bw2_ss1LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss1LR.Coefficients.Estimate;
            pvals_ch_bw2_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss1LR.Coefficients.pValue<0.05);

            ctwcclm_noch_bw2_ss1LL=fitlm([vccoffSVL_bw2(ss1LL) vccoffSVR_bw2(ss1LL)], vctccspk(ss1LL));
            zetas_noch_bw2_ss1LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw2_ss1LL.Residuals.Raw; nan(Ntrs12-length(ss1LL),1)];
            betas_noch_bw2_ss1LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw2_ss1LL.Coefficients.Estimate;
            pvals_noch_bw2_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw2_ss1LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw2_ss1LR=fitlm([vccoffSVL_bw2(ss1LR) vccoffSVR_bw2(ss1LR)], vctccspk(ss1LR));
            zetas_noch_bw2_ss1LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw2_ss1LR.Residuals.Raw; nan(Ntrs12-length(ss1LR),1)];
            betas_noch_bw2_ss1LR_mat(:,cc,tt,sg)=ctwcclm_noch_bw2_ss1LR.Coefficients.Estimate;
            pvals_noch_bw2_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw2_ss1LR.Coefficients.pValue<0.05);

            ctwcclm_ch_bw2_ss1LL_res=fitlm(vccoffchR(ss1LL),ctwcclm_noch_bw2_ss1LL.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw2_ss1LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss1LL_res.Coefficients.Estimate;
            ppvals_ch_bw2_ss1LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss1LL_res.Coefficients.pValue<0.05);
            ctwcclm_ch_bw2_ss1LR_res=fitlm(vccoffchR(ss1LR),ctwcclm_noch_bw2_ss1LR.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw2_ss1LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss1LR_res.Coefficients.Estimate;
            ppvals_ch_bw2_ss1LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss1LR_res.Coefficients.pValue<0.05);

            % FITTING LM FOR SV: CH/NOCH LL/LR (SS2, BW2)
            ctwcclm_ch_bw2_ss2LL=fitlm([vccoffSVL_bw2(ss2LL) vccoffSVR_bw2(ss2LL) vccoffchR(ss2LL)], vctccspk(ss2LL));
            %zetas_ch_bw2_ss2LL_mat(:,cc,tt,sg)=[ctwcclm_ch_bw2_ss2LL.Residuals.Raw; nan(Ntrs12-length(ss2LL),1)];
            betas_ch_bw2_ss2LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss2LL.Coefficients.Estimate;
            pvals_ch_bw2_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss2LL.Coefficients.pValue<0.05);
            ctwcclm_ch_bw2_ss2LR=fitlm([vccoffSVL_bw2(ss2LR) vccoffSVR_bw2(ss2LR) vccoffchR(ss2LR)], vctccspk(ss2LR));
            %zetas_ch_bw2_ss2LR_mat(:,cc,tt,sg)=[ctwcclm_ch_bw2_ss2LR.Residuals.Raw; nan(Ntrs12-length(ss2LR),1)];
            betas_ch_bw2_ss2LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss2LR.Coefficients.Estimate;
            pvals_ch_bw2_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss2LR.Coefficients.pValue<0.05);

            ctwcclm_noch_bw2_ss2LL=fitlm([vccoffSVL_bw2(ss2LL) vccoffSVR_bw2(ss2LL)], vctccspk(ss2LL));
            zetas_noch_bw2_ss2LL_mat(:,cc,tt,sg)=[ctwcclm_noch_bw2_ss2LL.Residuals.Raw; nan(Ntrs12-length(ss2LL),1)];
            betas_noch_bw2_ss2LL_mat(:,cc,tt,sg)=ctwcclm_noch_bw2_ss2LL.Coefficients.Estimate;
            pvals_noch_bw2_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_noch_bw2_ss2LL.Coefficients.pValue<0.05);
            ctwcclm_noch_bw2_ss2LR=fitlm([vccoffSVL_bw2(ss2LR) vccoffSVR_bw2(ss2LR)], vctccspk(ss2LR));
            zetas_noch_bw2_ss2LR_mat(:,cc,tt,sg)=[ctwcclm_noch_bw2_ss2LR.Residuals.Raw; nan(Ntrs12-length(ss2LR),1)];
            betas_noch_bw2_ss2LR_mat(:,cc,tt,sg)=ctwcclm_noch_bw2_ss2LR.Coefficients.Estimate;
            pvals_noch_bw2_ss2LR_mat(:,cc,tt,sg)=ctwcclm_noch_bw2_ss2LR.Coefficients.pValue;

            ctwcclm_ch_bw2_ss2LL_res=fitlm(vccoffchR(ss2LL),ctwcclm_noch_bw2_ss2LL.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw2_ss2LL_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss2LL_res.Coefficients.Estimate;
            ppvals_ch_bw2_ss2LL_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss2LL_res.Coefficients.pValue<0.05);
            ctwcclm_ch_bw2_ss2LR_res=fitlm(vccoffchR(ss2LR),ctwcclm_noch_bw2_ss2LR.Residuals.Raw,'Intercept',useinterceptflag==1);
            bbetas_ch_bw2_ss2LR_mat(:,cc,tt,sg)=ctwcclm_ch_bw2_ss2LR_res.Coefficients.Estimate;
            ppvals_ch_bw2_ss2LR_mat(:,cc,tt,sg)=uint8(ctwcclm_ch_bw2_ss2LR_res.Coefficients.pValue<0.05);

            % SAVING CHOICES TO BE USED FOR PREFERRED/NON-PREF. CHOICE
            offch_ss0LL_mat(:,cc,tt,sg)=[vccoffchR_flip(ss0LL); nan(Ntrs12-length(ss0LL),1)];
            offch_ss0LR_mat(:,cc,tt,sg)=[vccoffchR_flip(ss0LR); nan(Ntrs12-length(ss0LR),1)];
            offch_ss1LL_mat(:,cc,tt,sg)=[vccoffchR_flip(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offch_ss1LR_mat(:,cc,tt,sg)=[vccoffchR_flip(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offch_ss2LL_mat(:,cc,tt,sg)=[vccoffchR_flip(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offch_ss2LR_mat(:,cc,tt,sg)=[vccoffchR_flip(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            offch_ss0LL_noflip_mat(:,cc,tt,sg)=[vccoffchR(ss0LL); nan(Ntrs12-length(ss0LL),1)];
            offch_ss0LR_noflip_mat(:,cc,tt,sg)=[vccoffchR(ss0LR); nan(Ntrs12-length(ss0LR),1)];
            offch_ss1LL_noflip_mat(:,cc,tt,sg)=[vccoffchR(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offch_ss1LR_noflip_mat(:,cc,tt,sg)=[vccoffchR(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offch_ss2LL_noflip_mat(:,cc,tt,sg)=[vccoffchR(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offch_ss2LR_noflip_mat(:,cc,tt,sg)=[vccoffchR(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            
            offEVL_ss0LL_mat(:,cc,tt,sg)=[vccoffEVL(ss0LL);  nan(Ntrs12-length(ss0LL),1)];
            offEVL_ss0LR_mat(:,cc,tt,sg)=[vccoffEVL(ss0LR);  nan(Ntrs12-length(ss0LR),1)];
            offEVL_ss1LL_mat(:,cc,tt,sg)=[vccoffEVL(ss1LL);  nan(Ntrs12-length(ss1LL),1)];
            offEVL_ss1LR_mat(:,cc,tt,sg)=[vccoffEVL(ss1LR);  nan(Ntrs12-length(ss1LR),1)];
            offEVL_ss2LL_mat(:,cc,tt,sg)=[vccoffEVL(ss2LL);  nan(Ntrs12-length(ss2LL),1)];
            offEVL_ss2LR_mat(:,cc,tt,sg)=[vccoffEVL(ss2LR);  nan(Ntrs12-length(ss2LR),1)];
            
            offEVR_ss0LL_mat(:,cc,tt,sg)=[vccoffEVR(ss0LL);  nan(Ntrs12-length(ss0LL),1)];
            offEVR_ss0LR_mat(:,cc,tt,sg)=[vccoffEVR(ss0LR);  nan(Ntrs12-length(ss0LR),1)];
            offEVR_ss1LL_mat(:,cc,tt,sg)=[vccoffEVR(ss1LL);  nan(Ntrs12-length(ss1LL),1)];
            offEVR_ss1LR_mat(:,cc,tt,sg)=[vccoffEVR(ss1LR);  nan(Ntrs12-length(ss1LR),1)];
            offEVR_ss2LL_mat(:,cc,tt,sg)=[vccoffEVR(ss2LL);  nan(Ntrs12-length(ss2LL),1)];
            offEVR_ss2LR_mat(:,cc,tt,sg)=[vccoffEVR(ss2LR);  nan(Ntrs12-length(ss2LR),1)];
            

            offSVL_bw0ss0LL_mat(:,cc,tt,sg)=[vccoffSVL_bw0(ss0LL); nan(Ntrs12-length(ss0LL),1)];
            offSVL_bw0ss0LR_mat(:,cc,tt,sg)=[vccoffSVL_bw0(ss0LR); nan(Ntrs12-length(ss0LR),1)];
            offSVL_bw0ss1LL_mat(:,cc,tt,sg)=[vccoffSVL_bw0(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offSVL_bw0ss1LR_mat(:,cc,tt,sg)=[vccoffSVL_bw0(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offSVL_bw0ss2LL_mat(:,cc,tt,sg)=[vccoffSVL_bw0(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offSVL_bw0ss2LR_mat(:,cc,tt,sg)=[vccoffSVL_bw0(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            offSVL_bw1ss1LL_mat(:,cc,tt,sg)=[vccoffSVL_bw1(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offSVL_bw1ss1LR_mat(:,cc,tt,sg)=[vccoffSVL_bw1(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offSVL_bw1ss2LL_mat(:,cc,tt,sg)=[vccoffSVL_bw1(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offSVL_bw1ss2LR_mat(:,cc,tt,sg)=[vccoffSVL_bw1(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            offSVL_bw2ss1LL_mat(:,cc,tt,sg)=[vccoffSVL_bw2(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offSVL_bw2ss1LR_mat(:,cc,tt,sg)=[vccoffSVL_bw2(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offSVL_bw2ss2LL_mat(:,cc,tt,sg)=[vccoffSVL_bw2(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offSVL_bw2ss2LR_mat(:,cc,tt,sg)=[vccoffSVL_bw2(ss2LR); nan(Ntrs12-length(ss2LR),1)];

            offSVR_bw0ss0LL_mat(:,cc,tt,sg)=[vccoffSVR_bw0(ss0LL); nan(Ntrs12-length(ss0LL),1)];
            offSVR_bw0ss0LR_mat(:,cc,tt,sg)=[vccoffSVR_bw0(ss0LR); nan(Ntrs12-length(ss0LR),1)];
            offSVR_bw0ss1LL_mat(:,cc,tt,sg)=[vccoffSVR_bw0(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offSVR_bw0ss1LR_mat(:,cc,tt,sg)=[vccoffSVR_bw0(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offSVR_bw0ss2LL_mat(:,cc,tt,sg)=[vccoffSVR_bw0(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offSVR_bw0ss2LR_mat(:,cc,tt,sg)=[vccoffSVR_bw0(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            offSVR_bw1ss1LL_mat(:,cc,tt,sg)=[vccoffSVR_bw1(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offSVR_bw1ss1LR_mat(:,cc,tt,sg)=[vccoffSVR_bw1(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offSVR_bw1ss2LL_mat(:,cc,tt,sg)=[vccoffSVR_bw1(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offSVR_bw1ss2LR_mat(:,cc,tt,sg)=[vccoffSVR_bw1(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            offSVR_bw2ss1LL_mat(:,cc,tt,sg)=[vccoffSVR_bw2(ss1LL); nan(Ntrs12-length(ss1LL),1)];
            offSVR_bw2ss1LR_mat(:,cc,tt,sg)=[vccoffSVR_bw2(ss1LR); nan(Ntrs12-length(ss1LR),1)];
            offSVR_bw2ss2LL_mat(:,cc,tt,sg)=[vccoffSVR_bw2(ss2LL); nan(Ntrs12-length(ss2LL),1)];
            offSVR_bw2ss2LR_mat(:,cc,tt,sg)=[vccoffSVR_bw2(ss2LR); nan(Ntrs12-length(ss2LR),1)];
            
        end
        warning('on','all')
    end
    fprintf('%s - done: %d/%d, %2.1f%%\n',datetime('now','Format','HH:mm'),cc,Nccs,cc/Nccs.*100);
end
fprintf('Subgroup generation %d out of %d done.\n',sg,Nsubg)
end
    save(['fitlmstats_chR_evR_evL_spk_span_' num2str(Nspan) 'ms_lag_' num2str(Nlag) 'ms_' date '.mat'],'curr_rng','offch*','zeta*','beta*','bbeta*','pval*','ppval*','-v7.3');
else
    load('fitlmstats_chR_evR_evL_spk_span_200ms_lag_10ms_01-Mar-2024.mat');
end

disp('done')
error('force stop here')
%% FIG4E - CORRELATION BETWEEN CHOICE AND SV WEIGHTS
close all;

tbins_to_plot=1:58;%21:58;
subg_to_plot=1;
plot_sigbetas=1;

% bw0_ss1
sig_betas_ch_bw0_ss1LL_mat=betas_ch_bw0_ss1LL_mat; 
sig_bbetas_ch_bw0_ss1LL_mat=bbetas_ch_bw0_ss1LL_mat;
sig_betas_ch_bw0_ss1LR_mat=betas_ch_bw0_ss1LR_mat;
sig_bbetas_ch_bw0_ss1LR_mat=bbetas_ch_bw0_ss1LR_mat;
% bw0_ss2
sig_betas_ch_bw0_ss2LL_mat=betas_ch_bw0_ss2LL_mat;
sig_bbetas_ch_bw0_ss2LL_mat=bbetas_ch_bw0_ss2LL_mat;
sig_betas_ch_bw0_ss2LR_mat=betas_ch_bw0_ss2LR_mat;
sig_bbetas_ch_bw0_ss2LR_mat=bbetas_ch_bw0_ss2LR_mat;
% bw1_ss1
sig_betas_ch_bw1_ss1LL_mat=betas_ch_bw1_ss1LL_mat;
sig_bbetas_ch_bw1_ss1LL_mat=bbetas_ch_bw1_ss1LL_mat;
sig_betas_ch_bw1_ss1LR_mat=betas_ch_bw1_ss1LR_mat;
sig_bbetas_ch_bw1_ss1LR_mat=bbetas_ch_bw1_ss1LR_mat;
% bw2_ss2
sig_betas_ch_bw2_ss2LL_mat=betas_ch_bw2_ss2LL_mat;
sig_bbetas_ch_bw2_ss2LL_mat=bbetas_ch_bw2_ss2LL_mat;
sig_betas_ch_bw2_ss2LR_mat=betas_ch_bw2_ss2LR_mat;
sig_bbetas_ch_bw2_ss2LR_mat=bbetas_ch_bw2_ss2LR_mat;
% bw1_ss2
sig_betas_ch_bw1_ss2LL_mat=betas_ch_bw1_ss2LL_mat;
sig_bbetas_ch_bw1_ss2LL_mat=bbetas_ch_bw1_ss2LL_mat;
sig_betas_ch_bw1_ss2LR_mat=betas_ch_bw1_ss2LR_mat;
sig_bbetas_ch_bw1_ss2LR_mat=bbetas_ch_bw1_ss2LR_mat;
% bw2_ss1
sig_betas_ch_bw2_ss1LL_mat=betas_ch_bw2_ss1LL_mat;
sig_bbetas_ch_bw2_ss1LL_mat=bbetas_ch_bw2_ss1LL_mat;
sig_betas_ch_bw2_ss1LR_mat=betas_ch_bw2_ss1LR_mat;
sig_bbetas_ch_bw2_ss1LR_mat=bbetas_ch_bw2_ss1LR_mat;

if plot_sigbetas
    % bw0_ss1
    sig_betas_ch_bw0_ss1LL_mat(pvals_ch_bw0_ss1LL_mat==0)=nan;
    sig_bbetas_ch_bw0_ss1LL_mat(ppvals_ch_bw0_ss1LL_mat==0)=nan;
    sig_betas_ch_bw0_ss1LR_mat(pvals_ch_bw0_ss1LR_mat==0)=nan;
    sig_bbetas_ch_bw0_ss1LR_mat(ppvals_ch_bw0_ss1LR_mat==0)=nan;
    % bw0_ss2
    sig_betas_ch_bw0_ss2LL_mat(pvals_ch_bw0_ss2LL_mat==0)=nan;
    sig_bbetas_ch_bw0_ss2LL_mat(ppvals_ch_bw0_ss2LL_mat==0)=nan;
    sig_betas_ch_bw0_ss2LR_mat(pvals_ch_bw0_ss2LR_mat==0)=nan;
    sig_bbetas_ch_bw0_ss2LR_mat(ppvals_ch_bw0_ss2LR_mat==0)=nan;
    % bw1_ss1
    sig_betas_ch_bw1_ss1LL_mat(pvals_ch_bw1_ss1LL_mat==0)=nan;
    sig_bbetas_ch_bw1_ss1LL_mat(ppvals_ch_bw1_ss1LL_mat==0)=nan;
    sig_betas_ch_bw1_ss1LR_mat(pvals_ch_bw1_ss1LR_mat==0)=nan;
    sig_bbetas_ch_bw1_ss1LR_mat(ppvals_ch_bw1_ss1LR_mat==0)=nan;
    % bw2_ss2
    sig_betas_ch_bw2_ss2LL_mat(pvals_ch_bw2_ss2LL_mat==0)=nan;
    sig_bbetas_ch_bw2_ss2LL_mat(ppvals_ch_bw2_ss2LL_mat==0)=nan;
    sig_betas_ch_bw2_ss2LR_mat(pvals_ch_bw2_ss2LR_mat==0)=nan;
    sig_bbetas_ch_bw2_ss2LR_mat(ppvals_ch_bw2_ss2LR_mat==0)=nan;
    % bw1_ss2
    sig_betas_ch_bw1_ss2LL_mat(pvals_ch_bw1_ss2LL_mat==0)=nan;
    sig_bbetas_ch_bw1_ss2LL_mat(ppvals_ch_bw1_ss2LL_mat==0)=nan;
    sig_betas_ch_bw1_ss2LR_mat(pvals_ch_bw1_ss2LR_mat==0)=nan;
    sig_bbetas_ch_bw1_ss2LR_mat(ppvals_ch_bw1_ss2LR_mat==0)=nan;
    % bw2_ss1
    sig_betas_ch_bw2_ss1LL_mat(pvals_ch_bw2_ss1LL_mat==0)=nan;
    sig_bbetas_ch_bw2_ss1LL_mat(ppvals_ch_bw2_ss1LL_mat==0)=nan;
    sig_betas_ch_bw2_ss1LR_mat(pvals_ch_bw2_ss1LR_mat==0)=nan;
    sig_bbetas_ch_bw2_ss1LR_mat(ppvals_ch_bw2_ss1LR_mat==0)=nan;
end

% bw0_ss1
beta1_bw0_ss1LL=nanmean(squeeze(sig_betas_ch_bw0_ss1LL_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw0_ss1LL=nanmean(squeeze(sig_betas_ch_bw0_ss1LL_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw0_ss1LL=nanmean(squeeze(sig_betas_ch_bw0_ss1LL_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw0_ss1LL=nanmean(squeeze(sig_bbetas_ch_bw0_ss1LL_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
beta1_bw0_ss1LR=nanmean(squeeze(sig_betas_ch_bw0_ss1LR_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw0_ss1LR=nanmean(squeeze(sig_betas_ch_bw0_ss1LR_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw0_ss1LR=nanmean(squeeze(sig_betas_ch_bw0_ss1LR_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw0_ss1LR=nanmean(squeeze(sig_bbetas_ch_bw0_ss1LR_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
% bw0_ss2
beta1_bw0_ss2LL=nanmean(squeeze(sig_betas_ch_bw0_ss2LL_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw0_ss2LL=nanmean(squeeze(sig_betas_ch_bw0_ss2LL_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw0_ss2LL=nanmean(squeeze(sig_betas_ch_bw0_ss2LL_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw0_ss2LL=nanmean(squeeze(sig_bbetas_ch_bw0_ss2LL_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
beta1_bw0_ss2LR=nanmean(squeeze(sig_betas_ch_bw0_ss2LR_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw0_ss2LR=nanmean(squeeze(sig_betas_ch_bw0_ss2LR_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw0_ss2LR=nanmean(squeeze(sig_betas_ch_bw0_ss2LR_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw0_ss2LR=nanmean(squeeze(sig_bbetas_ch_bw0_ss2LR_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
% bw1_ss1
beta1_bw1_ss1LL=nanmean(squeeze(sig_betas_ch_bw1_ss1LL_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw1_ss1LL=nanmean(squeeze(sig_betas_ch_bw1_ss1LL_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw1_ss1LL=nanmean(squeeze(sig_betas_ch_bw1_ss1LL_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw1_ss1LL=nanmean(squeeze(sig_bbetas_ch_bw1_ss1LL_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
beta1_bw1_ss1LR=nanmean(squeeze(sig_betas_ch_bw1_ss1LR_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw1_ss1LR=nanmean(squeeze(sig_betas_ch_bw1_ss1LR_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw1_ss1LR=nanmean(squeeze(sig_betas_ch_bw1_ss1LR_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw1_ss1LR=nanmean(squeeze(sig_bbetas_ch_bw1_ss1LR_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
% bw2_ss2
beta1_bw2_ss2LL=nanmean(squeeze(sig_betas_ch_bw2_ss2LL_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw2_ss2LL=nanmean(squeeze(sig_betas_ch_bw2_ss2LL_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw2_ss2LL=nanmean(squeeze(sig_betas_ch_bw2_ss2LL_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw2_ss2LL=nanmean(squeeze(sig_bbetas_ch_bw2_ss2LL_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
beta1_bw2_ss2LR=nanmean(squeeze(sig_betas_ch_bw2_ss2LR_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw2_ss2LR=nanmean(squeeze(sig_betas_ch_bw2_ss2LR_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw2_ss2LR=nanmean(squeeze(sig_betas_ch_bw2_ss2LR_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw2_ss2LR=nanmean(squeeze(sig_bbetas_ch_bw2_ss2LR_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
% bw1_ss2
beta1_bw1_ss2LL=nanmean(squeeze(sig_betas_ch_bw1_ss2LL_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw1_ss2LL=nanmean(squeeze(sig_betas_ch_bw1_ss2LL_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw1_ss2LL=nanmean(squeeze(sig_betas_ch_bw1_ss2LL_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw1_ss2LL=nanmean(squeeze(sig_bbetas_ch_bw1_ss2LL_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
beta1_bw1_ss2LR=nanmean(squeeze(sig_betas_ch_bw1_ss2LR_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw1_ss2LR=nanmean(squeeze(sig_betas_ch_bw1_ss2LR_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw1_ss2LR=nanmean(squeeze(sig_betas_ch_bw1_ss2LR_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw1_ss2LR=nanmean(squeeze(sig_bbetas_ch_bw1_ss2LR_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
% bw2_ss1
beta1_bw2_ss1LL=nanmean(squeeze(sig_betas_ch_bw2_ss1LL_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw2_ss1LL=nanmean(squeeze(sig_betas_ch_bw2_ss1LL_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw2_ss1LL=nanmean(squeeze(sig_betas_ch_bw2_ss1LL_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw2_ss1LL=nanmean(squeeze(sig_bbetas_ch_bw2_ss1LL_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);
beta1_bw2_ss1LR=nanmean(squeeze(sig_betas_ch_bw2_ss1LR_mat(2,:,tbins_to_plot,subg_to_plot)),2:3);
beta2_bw2_ss1LR=nanmean(squeeze(sig_betas_ch_bw2_ss1LR_mat(3,:,tbins_to_plot,subg_to_plot)),2:3);
beta3_bw2_ss1LR=nanmean(squeeze(sig_betas_ch_bw2_ss1LR_mat(4,:,tbins_to_plot,subg_to_plot)),2:3);
bbeta3_bw2_ss1LR=nanmean(squeeze(sig_bbetas_ch_bw2_ss1LR_mat(1,:,tbins_to_plot,subg_to_plot)),2:3);

corrmethstr={'Pearson','Spearman'};

for mm=1:2
    corrmeth=corrmethstr{mm};

    hf=fullfig([0 0 1/4 1/3]+(mm-1)*[.5 0 0 0]);
    subplot(2,2,1);
    correlplot([[beta3_bw0_ss2LL; beta3_bw0_ss1LL],[beta1_bw0_ss1LL; beta1_bw0_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_1');
    subplot(2,2,2);
    correlplot([[beta3_bw0_ss2LR; beta3_bw0_ss1LR],[beta1_bw0_ss1LR; beta1_bw0_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_1');
    subplot(2,2,3);
    correlplot([[beta3_bw0_ss2LL; beta3_bw0_ss1LL],[beta2_bw0_ss1LL; beta2_bw0_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_2');
    subplot(2,2,4);
    correlplot([[beta3_bw0_ss2LR; beta3_bw0_ss1LR],[beta2_bw0_ss1LR; beta2_bw0_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_2');
    supertitle('ss1-bw0/ss2-bw0',8);

    hf=fullfig([1/4 0 1/4 1/3]+(mm-1)*[.5 0 0 0]);
    subplot(2,2,1);
    correlplot([[bbeta3_bw0_ss2LL; bbeta3_bw0_ss1LL],[beta1_bw0_ss1LL; beta1_bw0_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_1');
    subplot(2,2,2);
    correlplot([[bbeta3_bw0_ss2LR; bbeta3_bw0_ss1LR],[beta1_bw0_ss1LR; beta1_bw0_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_1');
    subplot(2,2,3);
    correlplot([[bbeta3_bw0_ss2LL; bbeta3_bw0_ss1LL],[beta2_bw0_ss1LL; beta2_bw0_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_2');
    subplot(2,2,4);
    correlplot([[bbeta3_bw0_ss2LR; bbeta3_bw0_ss1LR],[beta2_bw0_ss1LR; beta2_bw0_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_2');
    supertitle('ss1-bw0/ss2-bw0',8);

    hf=fullfig([0 1/3 1/4 1/3]+(mm-1)*[.5 0 0 0]);
    subplot(2,2,1);
    correlplot([[beta3_bw1_ss1LL; beta3_bw2_ss2LL],[beta1_bw2_ss2LL; beta1_bw1_ss1LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_1');
    subplot(2,2,2);
    correlplot([[beta3_bw1_ss1LR; beta3_bw2_ss2LR],[beta1_bw2_ss2LR; beta1_bw1_ss1LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_1');
    subplot(2,2,3);
    correlplot([[beta3_bw1_ss1LL; beta3_bw2_ss2LL],[beta2_bw2_ss2LL; beta2_bw1_ss1LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_2');
    subplot(2,2,4);
    correlplot([[beta3_bw1_ss1LR; beta3_bw2_ss2LR],[beta2_bw2_ss2LR; beta2_bw1_ss1LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_2');
    supertitle('ss1-bw1/ss2-bw2',8);

    hf=fullfig([1/4 1/3 1/4 1/3]+(mm-1)*[.5 0 0 0]);
    subplot(2,2,1);
    correlplot([[bbeta3_bw1_ss1LL; bbeta3_bw2_ss2LL],[beta1_bw2_ss2LL; beta1_bw1_ss1LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_1');
    subplot(2,2,2);
    correlplot([[bbeta3_bw1_ss1LR; bbeta3_bw2_ss2LR],[beta1_bw2_ss2LR; beta1_bw1_ss1LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_1');
    subplot(2,2,3);
    correlplot([[bbeta3_bw1_ss1LL; bbeta3_bw2_ss2LL],[beta2_bw2_ss2LL; beta2_bw1_ss1LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_2');
    subplot(2,2,4);
    correlplot([[bbeta3_bw1_ss1LR; bbeta3_bw2_ss2LR],[beta2_bw2_ss2LR; beta2_bw1_ss1LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_2');
    supertitle('ss1-bw1/ss2-bw2',8);

    hf=fullfig([0 2/3 1/4 1/3]+(mm-1)*[.5 0 0 0]);
    subplot(2,2,1);
    correlplot([[beta3_bw1_ss2LL; beta3_bw2_ss1LL],[beta1_bw2_ss1LL; beta1_bw1_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_1');
    subplot(2,2,2);
    correlplot([[beta3_bw1_ss2LR; beta3_bw2_ss1LR],[beta1_bw2_ss1LR; beta1_bw1_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_1');
    subplot(2,2,3);
    correlplot([[beta3_bw1_ss2LL; beta3_bw2_ss1LL],[beta2_bw2_ss1LL; beta2_bw1_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_2');
    subplot(2,2,4);
    correlplot([[beta3_bw1_ss2LR; beta3_bw2_ss1LR],[beta2_bw2_ss1LR; beta2_bw1_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_3'); ylabel('\beta_2');
    supertitle('ss2-bw1/ss2-bw2',8);

    hf=fullfig([1/4 2/3 1/4 1/3]+(mm-1)*[.5 0 0 0]);
    subplot(2,2,1);
    correlplot([[bbeta3_bw1_ss2LL; bbeta3_bw2_ss1LL],[beta1_bw2_ss1LL; beta1_bw1_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_1');
    subplot(2,2,2);
    correlplot([[bbeta3_bw1_ss2LR; bbeta3_bw2_ss1LR],[beta1_bw2_ss1LR; beta1_bw1_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_1');
    subplot(2,2,3);
    correlplot([[bbeta3_bw1_ss2LL; bbeta3_bw2_ss1LL],[beta2_bw2_ss1LL; beta2_bw1_ss2LL]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookL subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_2');
    subplot(2,2,4);
    correlplot([[bbeta3_bw1_ss2LR; bbeta3_bw2_ss1LR],[beta2_bw2_ss1LR; beta2_bw1_ss2LR]],corrmeth,[.1 .1 .1],[-2.5e-3 2.5e-3],[-2.5e-3 2.5e-3]);
    title('LookR subsets 1:2'); xlabel('\beta_1\prime'); ylabel('\beta_2');
    supertitle('ss2-bw1/ss1-bw2',8);
end


%% FIG 4F - RESIDUALS MAGNITUDES
tbins_to_plot=1:58;%51:58;
Ntb2p=length(tbins_to_plot);
subg_to_plot=1;
addpath('./support_functions')
clc

%% BW0 SS0
close all
clearvars Zstruct*
plot_sigbetas=0;
Zstruct0_trwise=get_bw_ss('bw0','ss0',plot_sigbetas,tbins_to_plot,subg_to_plot);
plot_residuals_comparison_scatter(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0/',['_bw0_ss0_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
%plot_residuals_comparison_distributions(Zstruct0_trwise, ...
%                                  './zetas_comparisons/',['_bw0_ss0_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0/',['_bw0_ss0_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

plot_sigbetas=1;
Zstruct0_trwise=get_bw_ss('bw0','ss0',plot_sigbetas,tbins_to_plot,subg_to_plot);
plot_residuals_comparison_scatter(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0_sig/',['_bw0_ss0_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
%plot_residuals_comparison_distributions(Zstruct0_trwise, ...
%                                  './zetas_comparisons/',['_bw0_ss0_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0_sig/',['_bw0_ss0_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);


%% BW0 SS1-SS2
addpath('./support_functions')
close all
clearvars Zstruct*
plot_sigbetas=0;
Zstruct1_trwise=get_bw_ss('bw0','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct2_trwise=get_bw_ss('bw0','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);

% sum(~isnan(nanmean(Zstruct1_trwise.Z_b1_LL_trwise(:,1,:,1),3)))+...
% sum(~isnan(nanmean(Zstruct2_trwise.Z_b1_LL_trwise(:,1,:,1),3))) % b1>0 chL
% sum(~isnan(nanmean(Zstruct1_trwise.Z_b1_LL_trwise(:,1,:,3),3)))+...
% sum(~isnan(nanmean(Zstruct2_trwise.Z_b1_LL_trwise(:,1,:,3),3))) % b1>0 chR
% 
% sum(~isnan(nanmean(Zstruct0_trwise.Z_b1_LL_trwise(:,1,:,1),3))) % b1>0 chL
% sum(~isnan(nanmean(Zstruct0_trwise.Z_b1_LL_trwise(:,1,:,3),3))) % b1>0 chR

plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss12/',['_bw0_ss12_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss12/',['_bw0_ss12_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

plot_sigbetas=1;
Zstruct1_trwise=get_bw_ss('bw0','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct2_trwise=get_bw_ss('bw0','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);

plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss12_sig/',['_bw0_ss12_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss12_sig/',['_bw0_ss12_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

%% BW1SS1-BW2SS2
addpath('./support_functions')
close all
clearvars Zstruct*
plot_sigbetas=0;
Zstruct1_trwise=get_bw_ss('bw1','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct2_trwise=get_bw_ss('bw2','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);

plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss12/',['_bw12_ss12_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss12/',['_bw12_ss12_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);


plot_sigbetas=1;
Zstruct1_trwise=get_bw_ss('bw1','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct2_trwise=get_bw_ss('bw2','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);

plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss12_sig/',['_bw12_ss12_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss12_sig/',['_bw12_ss12_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);


%% BW1SS2-BW1SS2
addpath('support_functions')
close all
clearvars Zstruct*
plot_sigbetas=0;
Zstruct1_trwise=get_bw_ss('bw1','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct2_trwise=get_bw_ss('bw2','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);

plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21/',['_bw12_ss21_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

[~,~,~,ZLA]=plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21/',['_bw12_ss21_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

%%
plot_sigbetas=1;
Zstruct1_trwise=get_bw_ss('bw1','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct2_trwise=get_bw_ss('bw2','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);

plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21_sig/',['_bw12_ss21_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21_sig/',['_bw12_ss21_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);


%% BW0 SS0 DIFFICULT
close all
clearvars Zstruct*
plot_sigbetas=0; subselectEVdiff=1;
Zstruct0_trwise=get_bw_ss('bw0','ss0',plot_sigbetas,tbins_to_plot,subg_to_plot,subselectEVdiff);
plot_residuals_comparison_scatter(Zstruct0_trwise,...
                                  './zetaplot_comparisons_choice_prob_cellwiseplot_comparisons_choice_prob_cellwises_comparisons_cameraready/bw0_ss0_diff/',['_bw0_ss0_diff_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0_diff/',['_bw0_ss0_diff_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

plot_sigbetas=1;
Zstruct0_trwise=get_bw_ss('bw0','ss0',plot_sigbetas,tbins_to_plot,subg_to_plot);
plot_residuals_comparison_scatter(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0_diff_sig/',['_bw0_ss0_diff_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct0_trwise,...
                                  './zetas_comparisons_cameraready/bw0_ss0_diff_sig/',['_bw0_ss0_diff_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

%% BW1SS2-BW2SS1 DIFFICULT
close all
clearvars Zstruct*
plot_sigbetas=0; subselectEVdiff=1;
Zstruct1_trwise=get_bw_ss('bw1','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot,subselectEVdiff);
Zstruct2_trwise=get_bw_ss('bw2','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot,subselectEVdiff);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);
plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21_diff/',['_bw12_ss21_diff_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
[~,~,~,ZLD]=plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21_diff/',['_bw12_ss21_diff_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
%%
plot_sigbetas=1; subselectEVdiff=1;
Zstruct1_trwise=get_bw_ss('bw1','ss2',plot_sigbetas,tbins_to_plot,subg_to_plot,subselectEVdiff);
Zstruct2_trwise=get_bw_ss('bw2','ss1',plot_sigbetas,tbins_to_plot,subg_to_plot,subselectEVdiff);
Zstruct12_trwise = combine_subsets(Zstruct1_trwise,Zstruct2_trwise);
plot_residuals_comparison_scatter(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21_diff_sig/',['_bw12_ss21_diff_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);
plot_comparisons_choice_prob_cellwise(Zstruct12_trwise,...
                                  './zetas_comparisons_cameraready/bw12_ss21_diff_sig/',['_bw12_ss21_diff_sig_' num2str(tbins_to_plot(1)) '-' num2str(tbins_to_plot(end))]);

%%
load('fitlmstats_chR_evR_evL_spk_span_200ms_lag_10ms_11-Mar-2024.mat');
nAll =sum(~isnan(zetas_noch_bw1_ss2LL_mat(:)))+sum(~isnan(zetas_noch_bw1_ss2LR_mat(:)))+...
      sum(~isnan(zetas_noch_bw2_ss1LL_mat(:)))+sum(~isnan(zetas_noch_bw2_ss1LL_mat(:)));
load('fitlmstats_chR_evR_evL_spk_span_200ms_lag_10ms_22-Mar-2024.mat');
nDiff=sum(~isnan(zetas_noch_bw1_ss2LL_mat(:)))+sum(~isnan(zetas_noch_bw1_ss2LR_mat(:)))+...
      sum(~isnan(zetas_noch_bw2_ss1LL_mat(:)))+sum(~isnan(zetas_noch_bw2_ss1LL_mat(:)));

fprintf('Difficult trials are: %.3f%% of total trials.',nDiff/nAll*100)


%% SUPPORT FUNCTIONS BELOW THIS POINT

%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

%%  Reshape betas based on beta signs

function [zr_mat_trwise_split,zr_mat_trwise_merged,zr_mat_travg] = reshapez_choice_betas(...
                cv_mat,bv_mat,zv_mat,tbins_to_plot,subg_to_plot,plot_sigbetas,pv_mat)
if nargin<6
    plot_sigbetas=0;
end
[Ntrs12, Nccs, Ntb2p, Nsubg]=size(zv_mat(:,:,tbins_to_plot,subg_to_plot));
cmv0=cv_mat(:,:,tbins_to_plot,subg_to_plot);%initially ch=1 means chooseR
bmv0=repmat(reshape(squeeze(bv_mat(1,:,tbins_to_plot,subg_to_plot)),[1 Nccs Ntb2p Nsubg]),Ntrs12,1,1);
zmv0=zv_mat(:,:,tbins_to_plot,subg_to_plot);
bmv0=repmat(nanmean(bmv0,3),[1 1 Ntb2p Nsubg]); %time-average betas (optional)
bg0_nm=nan(size(bmv0)); bg0_nm(bmv0>0)=1;
bl0_nm=nan(size(bmv0)); bl0_nm(bmv0<0)=1;
cL0_nm=nan(size(cmv0));  cL0_nm(cmv0==0)=1;
cR0_nm=nan(size(cmv0));  cR0_nm(cmv0==1)=1;
if plot_sigbetas
    ps0=double(repmat(reshape(squeeze(pv_mat(1,:,tbins_to_plot,subg_to_plot)),[1 Nccs Ntb2p Nsubg]),Ntrs12,1,1)==1);
    ps0(ps0==0)=nan;
    zmv0=zmv0.*ps0;
end
%sum(~isnan(nanmedian(zmv0,[1 3])),2)
%sum(sum(~isnan(bg0_nm(1,:,:)),3)==0) %there is some time bins where no beta is >0
pref1=bg0_nm.*cL0_nm;
pref2=bl0_nm.*cL0_nm;
pref3=bg0_nm.*cR0_nm;
pref4=bl0_nm.*cR0_nm;

z0_bpcL=zmv0.*pref1;  z_bpcL=nanmean(z0_bpcL,[1 3]);%nanmedian(z0_bpcL,[1 3]);
z0_bncL=zmv0.*pref2;  z_bncL=nanmean(z0_bncL,[1 3]);%nanmedian(z0_bncL,[1 3]);
z0_bpcR=zmv0.*pref3;  z_bpcR=nanmean(z0_bpcR,[1 3]);%nanmedian(z0_bpcR,[1 3]);
z0_bncR=zmv0.*pref4;  z_bncR=nanmean(z0_bncR,[1 3]);%nanmedian(z0_bncR,[1 3]);

%z0_bpcL=zmv0.*pref1;  z_bpcL=nanmedian(z0_bpcL,[1 3]);
%z0_bncL=zmv0.*pref2;  z_bncL=nanmedian(z0_bncL,[1 3]);
%z0_bpcR=zmv0.*pref3;  z_bpcR=nanmedian(z0_bpcR,[1 3]);
%z0_bncR=zmv0.*pref4;  z_bncR=nanmedian(z0_bncR,[1 3]);
%(sum(~isnan([z_bpcL; z_bncR; z_bncL; z_bpcR]),2))
%z_bpcL=z_bpcL(:); z_bncL=z_bncL(:); z_bpcR=z_bpcR(:); z_bncR=z_bncR(:);
%Nz=[size(z_bpcL,1)+size(z_bncR,1) size(z_bpcR,1)+size(z_bncL,1)];
%zr_mat_travg=[[z_bpcL; z_bncR; nan(max(Nz)-Nz(1),1)] [z_bncL; z_bpcR; nan(max(Nz)-Nz(2),1)]];
%size([z_bpcL z_bpcR])
%zr_mat_travg=[[z_bpcL' z_bpcR']; [z_bncR' z_bncL']]'; % [2*Nccs, 2]
zr_mat_travg=[z_bpcL z_bncR; z_bpcR, z_bncL]; % [2*Nccs, 2]
%zr_mat_travg=cat(2,cat(1,z_bpcL,z_bncR),cat(1,z_bpcR,z_bncL));
%z_bpcL=z0_bpcL;%nanmedian(z0_bpcL,3);
%z_bncL=z0_bncL;%nanmedian(z0_bncL,3);
%z_bpcR=z0_bpcR;%nanmedian(z0_bpcR,3);
%z_bncR=z0_bncR;%nanmedian(z0_bncR,3);
%Nz=[size(z_bpcL,2)+size(z_bncR,2) size(z_bpcR,2)+size(z_bncL,2)];
%size(z0_bpcL)
zr_mat_trwise_merged=cat(4,cat(2,z0_bpcL,z0_bncR),cat(2,z0_bpcR,z0_bncL)); %[Ntrs12 2*Nccs Ntb2p 2]
zr_mat_trwise_split=cat(4,z0_bpcL,z0_bncR,z0_bpcR,z0_bncL); %[Ntrs12 Nccs Ntb2p 4]
%zr_mat_trwise=cat(4,[z_bpcL z_bncR; nan(max(Nz)-Nz(1),1,N)],[z_bncL z_bpcR; nan(max(Nz)-Nz(2),1)]);
%sum(~isnan(zr_mat),1)
end



%% Combine disjoint subsets into pools of betas

function  Zstruct_trwise = combine_subsets(Zstructi,Zstructj)

Z_b1_LL_bw_ssi_trwise=Zstructi.Z_b1_LL_trwise;
Z_b1_LR_bw_ssi_trwise=Zstructi.Z_b1_LR_trwise;
Z_b2_LL_bw_ssi_trwise=Zstructi.Z_b2_LL_trwise;
Z_b2_LR_bw_ssi_trwise=Zstructi.Z_b2_LR_trwise;

Z_b1_LL_bw_ssj_trwise=Zstructj.Z_b1_LL_trwise;
Z_b1_LR_bw_ssj_trwise=Zstructj.Z_b1_LL_trwise;
Z_b2_LL_bw_ssj_trwise=Zstructj.Z_b2_LL_trwise;
Z_b2_LR_bw_ssj_trwise=Zstructj.Z_b2_LL_trwise;

Z_b1_LL_trwise_ij=cat(4,Z_b1_LL_bw_ssi_trwise(:,:,:,1),Z_b1_LL_bw_ssj_trwise(:,:,:,2),...
                        Z_b1_LL_bw_ssj_trwise(:,:,:,3),Z_b1_LL_bw_ssi_trwise(:,:,:,4));
Z_b1_LR_trwise_ij=cat(4,Z_b1_LR_bw_ssi_trwise(:,:,:,1),Z_b1_LR_bw_ssj_trwise(:,:,:,2),...
                        Z_b1_LR_bw_ssj_trwise(:,:,:,3),Z_b1_LR_bw_ssi_trwise(:,:,:,4));
Z_b2_LL_trwise_ij=cat(4,Z_b2_LL_bw_ssi_trwise(:,:,:,1),Z_b2_LL_bw_ssj_trwise(:,:,:,2),...
                        Z_b2_LL_bw_ssj_trwise(:,:,:,3),Z_b2_LL_bw_ssi_trwise(:,:,:,4));
Z_b2_LR_trwise_ij=cat(4,Z_b2_LR_bw_ssi_trwise(:,:,:,1),Z_b2_LR_bw_ssj_trwise(:,:,:,2),...
                        Z_b2_LR_bw_ssj_trwise(:,:,:,3),Z_b2_LR_bw_ssi_trwise(:,:,:,4));

Z_b1_LL_trwise_ji=cat(4,Z_b1_LL_bw_ssj_trwise(:,:,:,1),Z_b1_LL_bw_ssi_trwise(:,:,:,2),...
                        Z_b1_LL_bw_ssi_trwise(:,:,:,3),Z_b1_LL_bw_ssj_trwise(:,:,:,4));
Z_b1_LR_trwise_ji=cat(4,Z_b1_LR_bw_ssj_trwise(:,:,:,1),Z_b1_LR_bw_ssi_trwise(:,:,:,2),...
                        Z_b1_LR_bw_ssi_trwise(:,:,:,3),Z_b1_LR_bw_ssj_trwise(:,:,:,4));
Z_b2_LL_trwise_ji=cat(4,Z_b2_LL_bw_ssj_trwise(:,:,:,1),Z_b2_LL_bw_ssi_trwise(:,:,:,2),...
                        Z_b2_LL_bw_ssi_trwise(:,:,:,3),Z_b2_LL_bw_ssj_trwise(:,:,:,4));
Z_b2_LR_trwise_ji=cat(4,Z_b2_LR_bw_ssj_trwise(:,:,:,1),Z_b2_LR_bw_ssi_trwise(:,:,:,2),...
                        Z_b2_LR_bw_ssi_trwise(:,:,:,3),Z_b2_LR_bw_ssj_trwise(:,:,:,4));

Zstruct_trwise.Z_b1_LL_trwise=cat(1,Z_b1_LL_trwise_ij,Z_b1_LL_trwise_ji);
Zstruct_trwise.Z_b1_LR_trwise=cat(1,Z_b1_LR_trwise_ij,Z_b1_LR_trwise_ji);
Zstruct_trwise.Z_b2_LL_trwise=cat(1,Z_b2_LL_trwise_ij,Z_b2_LL_trwise_ji);
Zstruct_trwise.Z_b2_LR_trwise=cat(1,Z_b2_LR_trwise_ij,Z_b2_LR_trwise_ji);

end

%% Reshape residuals based on sign of betas

function Zstruct_trwise=get_bw_ss(bwstr,subsetstr,plot_sigbetas,tbins_to_plot,subg_to_plot,select_diff_trials)
if nargin<6
select_diff_trials=0;
end
% bwstr = 'bw0', 'bw1', 'bw2'
% subsetstr = 'ss0'
if select_diff_trials
    disp('loading residuals for diff trials')
    load('fitlmstats_chR_evR_evL_spk_span_200ms_lag_10ms_22-Mar-2024.mat');
else
    disp('loading residuals for all trials')
    load('fitlmstats_chR_evR_evL_spk_span_200ms_lag_10ms_11-Mar-2024.mat');
end
betasLL_mat=eval(['betas_noch_' bwstr '_' subsetstr 'LL_mat']);
betasLR_mat=eval(['betas_noch_' bwstr '_' subsetstr 'LR_mat']);
zetasLL_mat=eval(['zetas_noch_' bwstr '_' subsetstr 'LL_mat']);
zetasLR_mat=eval(['zetas_noch_' bwstr '_' subsetstr 'LR_mat']);
pvalsLL_mat=eval(['pvals_noch_' bwstr '_' subsetstr 'LL_mat']);
pvalsLR_mat=eval(['pvals_noch_' bwstr '_' subsetstr 'LR_mat']);
offchLL_mat=eval(['offch_' subsetstr 'LL_mat']);
offchLR_mat=eval(['offch_' subsetstr 'LR_mat']);

if select_diff_trials
    load 'EVSVsubsets.mat';
    %SVL_LL_mat=eval(['offSVL_' bwstr swstr 'LL_mat']);
    %SVL_LR_mat=eval(['offSVL_' bwstr swstr 'LR_mat']);
    %SVR_LL_mat=eval(['offSVR_' bwstr swstr 'LL_mat']);
    %SVR_LR_mat=eval(['offSVR_' bwstr swstr 'LR_mat']);
    EVL_LL_mat=eval(['offEVL_' subsetstr 'LL_mat']);
    EVL_LR_mat=eval(['offEVL_' subsetstr 'LR_mat']);
    EVR_LL_mat=eval(['offEVR_' subsetstr 'LL_mat']);
    EVR_LR_mat=eval(['offEVR_' subsetstr 'LR_mat']);
    EVDiffLL_mask=nan(size(EVL_LL_mat)); 
    EVDiffLL_mask(abs(EVL_LL_mat-EVR_LL_mat)<1)=1;
    EVDiffLR_mask=nan(size(EVL_LR_mat)); 
    EVDiffLR_mask(abs(EVL_LR_mat-EVR_LR_mat)<1)=1;
    offchLL_mat=offchLL_mat.*EVDiffLL_mask;
    offchLR_mat=offchLR_mat.*EVDiffLR_mask;
end

Z_b1_LL_trwise=reshapez_choice_betas(offchLL_mat,betasLL_mat(2,:,:),zetasLL_mat,...
    tbins_to_plot,subg_to_plot,plot_sigbetas,pvalsLL_mat(2,:,:));
Z_b1_LR_trwise=reshapez_choice_betas(offchLR_mat,betasLR_mat(2,:,:),zetasLR_mat,...
    tbins_to_plot,subg_to_plot,plot_sigbetas,pvalsLR_mat(2,:,:));
Z_b2_LL_trwise=reshapez_choice_betas(offchLL_mat,betasLL_mat(3,:,:),zetasLL_mat,...
    tbins_to_plot,subg_to_plot,plot_sigbetas,pvalsLL_mat(3,:,:));
Z_b2_LR_trwise=reshapez_choice_betas(offchLR_mat,betasLR_mat(3,:,:),zetasLR_mat,...
    tbins_to_plot,subg_to_plot,plot_sigbetas,pvalsLR_mat(3,:,:));

Zstruct_trwise.Z_b1_LL_trwise=Z_b1_LL_trwise;
Zstruct_trwise.Z_b1_LR_trwise=Z_b1_LR_trwise;
Zstruct_trwise.Z_b2_LL_trwise=Z_b2_LL_trwise;
Zstruct_trwise.Z_b2_LR_trwise=Z_b2_LR_trwise;

end

%% Compute ROC and AUC

function [FPR,TPR,AUC,FP,TP,FN,TN]=get_fpr_tpr_auc(X1,X2,Y,plot_option)
if nargin<4
    if nargin<3
    Y=0:.001:1;
    end
    plot_option=0;
end
    X1=X1(:)';
    X2=X2(:)';
    FP=sum(X2>Y',2);
    TP=sum(X1>Y',2);
    FN=sum(X1<Y',2);
    TN=sum(X2<Y',2);
    FPR=FP./(TN+FP);
    TPR=TP./(FN+TP);
    AUC=max(cumtrapz(flip(FPR),flip(TPR)));
    if plot_option
    figure();
    plot(FPR,TPR); hold on; plot([0 1],[0 1]);
    text(.3,.1,sprintf('AUC: %.3f',AUC));
    end
end