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

addpath('./support_functions');

clearvars -except subjsData poolstack* pooltw*
close all; clc; %clear all;

recompute_mfit=0; 
reload_eye_spk_data=0;
Nlag=10; Nspan=200; Nshf=10; 
%Nlag=10; Nspan=200; Nshf=1000; %Set default
normstr='minus';  %normstr='no'; %normstr='div';

poissfitflag=0; 
if poissfitflag && strcmpi(normstr,'div'); error('please combine poissfitflag with normstr=no'); end

figures_folder='../figures';

 if ~exist('subjsData','var')
     subjsData(1).sbdata=getSubjectData(1);
     subjsData(2).sbdata=getSubjectData(2);
 end

if ~exist('poolstackvars','var')
    poolstackvars=getPoolStackVars(subjsData);
end

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

nccs=arrayfun(@(sn) poolstackvars(sn).ncells, 1:8);
ntrs=arrayfun(@(sn) size(poolstackvars(sn).offer1sd,1),1:8);

Nccs=sum(nccs);
Ntrs12=max(ntrs);

pooloffer1sd=[]; poolofferLev=[]; poolofferRev=[]; pooloffer1ev=[]; pooloffer2ev=[]; 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)]));
    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

if reload_eye_spk_data
    if ~exist('pooltweyepos','var')
        pooltweyepos=getTimewinsPoolEyePos();
    end
    if ~exist('pooltwspksq','var')
        pooltwspksq=getTimewinsPoolSpikeSeqs();
    end
    
    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 pooltwspksq pooltweyepos % to reduce memory usage keep pooled version
    
    fchg=find(poolcovhist>=.999); fchl=find(poolcovhist<.999); % use 99.9% coverage
    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
    
    %if sparsity==1 it means 1spike/ms that is 1000 spikes/s
    %nanmean(vectorise((nansum(poolspksqmat,3))./(nansum(poolspksqmat>=0,3))))*1000
    %nanstd(vectorise((nansum(poolspksqmat,3))./(nansum(poolspksqmat>=0,3))))*1000
    
    % smooth by movmean filtering, downsample to 20ms, 50Hz resolution
    Ntps_ds=ceil(sum(Ntps12)/Nlag);
    tmax_ds=nanmovmean(1:size(pooleposxmat,3),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 % to reduce memory usage keep downsampled version
    
    save('../data/downsampled_eye_spk_data.mat','Ntps_ds','tmax_ds','vtps','ntc*','fch*','Ntps12',...
        'poolcovhist','poolspksqmat_ds','pooleposxmat_ds','-v7.3');
else
    load('../data/downsampled_eye_spk_data.mat');
end

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)]'];

% Apply movmean independently across time windows
% checktps_ds=zeros(1,length(poolcovhist_ds));
% for tw=1:8
%     tpvec=zeros(1,length(poolcovhist));
%     tpvec(ntce(tw,1):ntce(tw,2))=1;
%     tpvec_ds=zeros(1,length(poolcovhist_ds));
%     tpvec_ds(ntce_ds(tw,1):ntce_ds(tw,2))=1;
%     for cc=1:Nccs
%         nnspktemp=nanmovmean(squeeze(poolspksqmat(:,cc,tpvec==1))',Nlag,[0 Nspan])';
%         poolspksqmat_ds(:,cc,tpvec_ds==1)=nnspktemp(:,1:sum(tpvec_ds==1));
%         nnepxtemp=nanmovmean(squeeze(pooleposxmat(:,cc,tpvec==1))',Nlag,Nlag)';
%         pooleposxmat_ds(:,cc,tpvec_ds==1)=nnepxtemp(:,1:sum(tpvec_ds==1));
%     end
%     checktps_ds(tpvec_ds==1)=checktps_ds(tpvec_ds==1)+tw;
% end

vtmax_ds=(tmax_ds.*vtps_ds); 
vtmax_ds(isnan(vtmax_ds))=[];
rtmax_ds=1:length(vtmax_ds);
rtmes_ds=arrayfun(@(jj) length(ntce_ds(jj,1):ntce_ds(jj,2)), 1:8);

% baseline normalization
if strcmpi(normstr,'no')
    poolspksqmat_ds_bln=poolspksqmat_ds;
    %pooleposxmat_ds_bln=pooleposxmat_ds;
elseif strcmpi(normstr,'div')
    poolspksqmat_ds_bln=poolspksqmat_ds./repmat(nanmean(poolspksqmat_ds(:,:,1:floor(400/Nlag)),3)+eps,[1 1 Ntps_ds]);
    %pooleposxmat_ds_bln=pooleposxmat_ds./repmat(nanmean(pooleposxmat_ds(:,:,1:floor(400/Nlag)),3)+eps,[1 1 Ntps_ds]);
elseif strcmpi(normstr,'minus')
    poolspksqmat_ds_bln=poolspksqmat_ds-repmat(nanmean(poolspksqmat_ds(:,:,1:floor(400/Nlag)),3),[1 1 Ntps_ds]);
    %pooleposxmat_ds_bln=pooleposxmat_ds-repmat(nanmean(pooleposxmat_ds(:,:,1:floor(400/Nlag)),3),[1 1 Ntps_ds]);
end

clearvars poolspksqmat_ds %pooleposxmat_ds % to reduce memory usage

%% Regress EV vs Spiking
if recompute_mfit
lmELpv=nan(Nccs,Ntps_ds);   lmELpv_shf=nan(Nccs,Ntps_ds,Nshf);   lmERpv=nan(Nccs,Ntps_ds);   lmERpv_shf=nan(Nccs,Ntps_ds,Nshf);
lmELpvLL=nan(Nccs,Ntps_ds); lmELpvLL_shf=nan(Nccs,Ntps_ds,Nshf); lmERpvLL=nan(Nccs,Ntps_ds); lmERpvLL_shf=nan(Nccs,Ntps_ds,Nshf);
lmELpvLR=nan(Nccs,Ntps_ds); lmELpvLR_shf=nan(Nccs,Ntps_ds,Nshf); lmERpvLR=nan(Nccs,Ntps_ds); lmERpvLR_shf=nan(Nccs,Ntps_ds,Nshf);
lmELR2=nan(Nccs,Ntps_ds);   lmELR2_shf=nan(Nccs,Ntps_ds,Nshf);   lmERR2=nan(Nccs,Ntps_ds);   lmERR2_shf=nan(Nccs,Ntps_ds,Nshf);
lmELR2LL=nan(Nccs,Ntps_ds); lmELR2LL_shf=nan(Nccs,Ntps_ds,Nshf); lmERR2LL=nan(Nccs,Ntps_ds); lmERR2LL_shf=nan(Nccs,Ntps_ds,Nshf);
lmELR2LR=nan(Nccs,Ntps_ds); lmELR2LR_shf=nan(Nccs,Ntps_ds,Nshf); lmERR2LR=nan(Nccs,Ntps_ds); lmERR2LR_shf=nan(Nccs,Ntps_ds,Nshf);
lmELb0=nan(Nccs,Ntps_ds);   lmELb0_shf=nan(Nccs,Ntps_ds,Nshf);   lmERb0=nan(Nccs,Ntps_ds);   lmERb0_shf=nan(Nccs,Ntps_ds,Nshf);
lmELb0LL=nan(Nccs,Ntps_ds); lmELb0LL_shf=nan(Nccs,Ntps_ds,Nshf); lmERb0LL=nan(Nccs,Ntps_ds); lmERb0LL_shf=nan(Nccs,Ntps_ds,Nshf);
lmELb0LR=nan(Nccs,Ntps_ds); lmELb0LR_shf=nan(Nccs,Ntps_ds,Nshf); lmERb0LR=nan(Nccs,Ntps_ds); lmERb0LR_shf=nan(Nccs,Ntps_ds,Nshf);
lmELb1=nan(Nccs,Ntps_ds);   lmELb1_shf=nan(Nccs,Ntps_ds,Nshf);   lmERb1=nan(Nccs,Ntps_ds);   lmERb1_shf=nan(Nccs,Ntps_ds,Nshf);
lmELb1LL=nan(Nccs,Ntps_ds); lmELb1LL_shf=nan(Nccs,Ntps_ds,Nshf); lmERb1LL=nan(Nccs,Ntps_ds); lmERb1LL_shf=nan(Nccs,Ntps_ds,Nshf);
lmELb1LR=nan(Nccs,Ntps_ds); lmELb1LR_shf=nan(Nccs,Ntps_ds,Nshf); lmERb1LR=nan(Nccs,Ntps_ds); lmERb1LR_shf=nan(Nccs,Ntps_ds,Nshf);
lmntrsLL=nan(Nccs,Ntps_ds); lmntrsLR=nan(Nccs,Ntps_ds);
for cc=1:Nccs
    ccoff1ev=pooloffer1ev(:,cc);
    ccoff2ev=pooloffer2ev(:,cc);
%    ccoff1sd=pooloffer1sd(:,cc);
    ccspk=squeeze(poolspksqmat_ds_bln(:,cc,:));   % For linear fit
%    ccspk=squeeze(poolspksqmat_ds(:,cc,:));      % Force NO BL correction
    ccepx=squeeze(pooleposxmat_ds(:,cc,:));      % Force NO BL correction
    rmallnans=isnan(ccoff1ev);
    rmtrials=rmallnans;%(rmallnans | ccoff1rw~=1 | ccoff2rw~=1);
%    rmtrials=(rmallnans | ccoff1sd==1); % runs only first Left
%    rmtrials=(rmallnans | ccoff1sd==0); % runs only first Right
    ccspk(rmtrials,:)=[];
    ccepx(rmtrials,:)=[];
    ccoff1ev(rmtrials)=[];
    ccoff2ev(rmtrials)=[];    
    ccntrs=size(ccspk,1);
        
    parfor tt=1:Ntps_ds
    %for tt=1:Ntps_ds
        warning('off','all')
        if vtps_ds(tt)==1 % discards low trials times
            ctspk=squeeze(ccspk(:,tt));
            ctepx=squeeze(ccepx(:,tt));

            if ~poissfitflag
                ctlm1=fitlm(ccoff1ev,ctspk);
                ctlm2=fitlm(ccoff2ev,ctspk);
            else
                ctlm1=fitglm(ccoff1ev,ctspk,'linear','Distribution','poisson');
                ctlm2=fitglm(ccoff2ev,ctspk,'linear','Distribution','poisson');
            end
            lmELpv(cc,tt)=ctlm1.Coefficients.pValue(2);
            lmERpv(cc,tt)=ctlm2.Coefficients.pValue(2);
            lmELR2(cc,tt)=ctlm1.Rsquared.Ordinary;
            lmERR2(cc,tt)=ctlm2.Rsquared.Ordinary;
            lmELb0(cc,tt)=ctlm1.Coefficients.Estimate(1);
            lmERb0(cc,tt)=ctlm2.Coefficients.Estimate(1);
            lmELb1(cc,tt)=ctlm1.Coefficients.Estimate(2);
            lmERb1(cc,tt)=ctlm2.Coefficients.Estimate(2);
            
            isLL=find(ctepx<0);
            niLL=length(isLL);
            if ~poissfitflag
                ctlm1LL=fitlm(ccoff1ev(isLL),ctspk(isLL));
                ctlm2LL=fitlm(ccoff2ev(isLL),ctspk(isLL));
            else
                ctlm1LL=fitglm(ccoff1ev(isLL),ctspk(isLL),'linear','Distribution','poisson');
                ctlm2LL=fitglm(ccoff2ev(isLL),ctspk(isLL),'linear','Distribution','poisson');
            end
            lmELpvLL(cc,tt)=ctlm1LL.Coefficients.pValue(2);
            lmERpvLL(cc,tt)=ctlm2LL.Coefficients.pValue(2);
            lmELR2LL(cc,tt)=ctlm1LL.Rsquared.Ordinary;
            lmERR2LL(cc,tt)=ctlm2LL.Rsquared.Ordinary;
            lmELb0LL(cc,tt)=ctlm1LL.Coefficients.Estimate(1);
            lmERb0LL(cc,tt)=ctlm2LL.Coefficients.Estimate(1);
            lmELb1LL(cc,tt)=ctlm1LL.Coefficients.Estimate(2);
            lmERb1LL(cc,tt)=ctlm2LL.Coefficients.Estimate(2);
            lmntrsLL(cc,tt)=niLL;
            
            isLR=find(ctepx>0);
            niLR=length(isLR);
            if ~poissfitflag
                ctlm1LR=fitlm(ccoff1ev(isLR),ctspk(isLR));
                ctlm2LR=fitlm(ccoff2ev(isLR),ctspk(isLR));
            else
                ctlm1LR=fitglm(ccoff1ev(isLR),ctspk(isLR),'linear','Distribution','poisson');
                ctlm2LR=fitglm(ccoff2ev(isLR),ctspk(isLR),'linear','Distribution','poisson');
            end
            lmELpvLR(cc,tt)=ctlm1LR.Coefficients.pValue(2);
            lmERpvLR(cc,tt)=ctlm2LR.Coefficients.pValue(2);
            lmELR2LR(cc,tt)=ctlm1LR.Rsquared.Ordinary;
            lmERR2LR(cc,tt)=ctlm2LR.Rsquared.Ordinary;
            lmELb0LR(cc,tt)=ctlm1LR.Coefficients.Estimate(1);
            lmERb0LR(cc,tt)=ctlm2LR.Coefficients.Estimate(1);
            lmELb1LR(cc,tt)=ctlm1LR.Coefficients.Estimate(2);
            lmERb1LR(cc,tt)=ctlm2LR.Coefficients.Estimate(2);
            lmntrsLR(cc,tt)=niLR;
            
            for shf=1:Nshf
                i1shf=randperm(ccntrs);                                   %trial-shuffle
                i2shf=randperm(ccntrs);                                   %trial-shuffle
                isshf=randperm(ccntrs);                                   %trial-shuffle
                if ~poissfitflag
                    ctlm1_shf=fitlm(ccoff1ev(i1shf),ctspk(isshf));            %trial-shuffle
                    ctlm2_shf=fitlm(ccoff2ev(i2shf),ctspk(isshf));            %trial-shuffle
                else
                    ctlm1_shf=fitglm(ccoff1ev(i1shf),ctspk(isshf),'linear','Distribution','poisson');
                    ctlm2_shf=fitglm(ccoff2ev(i2shf),ctspk(isshf),'linear','Distribution','poisson');
                end
                %ctlm1_shf=fitlm(ccoff1ev,squeeze(ccspk_shf(shf,:,tt))); %time-shuffle
                %ctlm2_shf=fitlm(ccoff2ev,squeeze(ccspk_shf(shf,:,tt))); %time-shuffle
                lmELpv_shf(cc,tt,shf)=ctlm1_shf.Coefficients.pValue(2);
                lmERpv_shf(cc,tt,shf)=ctlm2_shf.Coefficients.pValue(2);
                lmELR2_shf(cc,tt,shf)=ctlm1_shf.Rsquared.Ordinary;
                lmERR2_shf(cc,tt,shf)=ctlm2_shf.Rsquared.Ordinary;
                lmELb0_shf(cc,tt,shf)=ctlm1_shf.Coefficients.Estimate(1);
                lmERb0_shf(cc,tt,shf)=ctlm2_shf.Coefficients.Estimate(1);
                lmELb1_shf(cc,tt,shf)=ctlm1_shf.Coefficients.Estimate(2);
                lmERb1_shf(cc,tt,shf)=ctlm2_shf.Coefficients.Estimate(2);
                
                i1LLshf=isLL(randperm(niLL));                                        %trial-shuffle
                i2LLshf=isLL(randperm(niLL));                                        %trial-shuffle
                isLLshf=isLL(randperm(niLL));                                        %trial-shuffle 
                if ~poissfitflag
                    ctlm1LL_shf=fitlm(ccoff1ev(i1LLshf),ctspk(isLLshf));                 %trial-shuffle
                    ctlm2LL_shf=fitlm(ccoff2ev(i2LLshf),ctspk(isLLshf));                 %trial-shuffle
                else
                    ctlm1LL_shf=fitglm(ccoff1ev(i1LLshf),ctspk(isLLshf),'linear','Distribution','poisson');
                    ctlm2LL_shf=fitglm(ccoff2ev(i2LLshf),ctspk(isLLshf),'linear','Distribution','poisson');
                end
                %ctlm1LL_shf=fitlm(ccoff1ev(isLL),squeeze(ccspk_shf(shf,isLL,tt))); %time-shuffle
                %ctlm2LL_shf=fitlm(ccoff2ev(isLL),squeeze(ccspk_shf(shf,isLL,tt))); %time-shuffle
                lmELpvLL_shf(cc,tt,shf)=ctlm1LL_shf.Coefficients.pValue(2);
                lmERpvLL_shf(cc,tt,shf)=ctlm2LL_shf.Coefficients.pValue(2);
                lmELR2LL_shf(cc,tt,shf)=ctlm1LL_shf.Rsquared.Ordinary;
                lmERR2LL_shf(cc,tt,shf)=ctlm2LL_shf.Rsquared.Ordinary;
                lmELb0LL_shf(cc,tt,shf)=ctlm1LL_shf.Coefficients.Estimate(1);
                lmERb0LL_shf(cc,tt,shf)=ctlm2LL_shf.Coefficients.Estimate(1);
                lmELb1LL_shf(cc,tt,shf)=ctlm1LL_shf.Coefficients.Estimate(2);
                lmERb1LL_shf(cc,tt,shf)=ctlm2LL_shf.Coefficients.Estimate(2);
                
                i1LRshf=isLR(randperm(niLR));                                        %trial-shuffle
                i2LRshf=isLR(randperm(niLR));                                        %trial-shuffle
                isLRshf=isLR(randperm(niLR));                                        %trial-shuffle
                if ~poissfitflag
                    ctlm1LR_shf=fitlm(ccoff1ev(i1LRshf),ctspk(isLRshf));                 %trial-shuffle
                    ctlm2LR_shf=fitlm(ccoff2ev(i2LRshf),ctspk(isLRshf));                 %trial-shuffle
                else
                    ctlm1LR_shf=fitglm(ccoff1ev(i1LRshf),ctspk(isLRshf),'linear','Distribution','poisson');
                    ctlm2LR_shf=fitglm(ccoff2ev(i2LRshf),ctspk(isLRshf),'linear','Distribution','poisson');
                end
                %ctlm1LR_shf=fitlm(ccoff1ev(isLR),squeeze(ccspk_shf(shf,isLR,tt))); %time-shuffle
                %ctlm2LR_shf=fitlm(ccoff2ev(isLR),squeeze(ccspk_shf(shf,isLR,tt))); %time-shuffle
                lmELpvLR_shf(cc,tt,shf)=ctlm1LR_shf.Coefficients.pValue(2);
                lmERpvLR_shf(cc,tt,shf)=ctlm2LR_shf.Coefficients.pValue(2);
                lmELR2LR_shf(cc,tt,shf)=ctlm1LR_shf.Rsquared.Ordinary;
                lmERR2LR_shf(cc,tt,shf)=ctlm2LR_shf.Rsquared.Ordinary;
                lmELb0LR_shf(cc,tt,shf)=ctlm1LR_shf.Coefficients.Estimate(1);
                lmERb0LR_shf(cc,tt,shf)=ctlm2LR_shf.Coefficients.Estimate(1);
                lmELb1LR_shf(cc,tt,shf)=ctlm1LR_shf.Coefficients.Estimate(2);
                lmERb1LR_shf(cc,tt,shf)=ctlm2LR_shf.Coefficients.Estimate(2);
            end
        end
        warning('on','all')
    end
    disp([sprintf('%s',datetime('now','Format','HH:mm')) ' - done: ' num2str(cc) '/' num2str(Nccs) ', ' num2str(cc/Nccs.*100) '%']);
end
    save(['fitlmstats_span_' num2str(Nspan) 'ms_lag_' num2str(Nlag) 'ms_' date '.mat'],'lmE*','lmntr*');
else
    load(['fitlmstats_span_' num2str(Nspan) 'ms_lag_' num2str(Nlag) 'ms_16-Feb-2023.mat']);
end


%% Plot colors
o1col=[55 188 52]/255;
o2col=[55 131 200]/255;
o3col=[121 121 121]/255;

o1LLcol=[50 148 48]/255;
o2LLcol=[55 97 135]/255;
o3LLcol=[43 43 43]/255;

o1LRcol=[135 215 133]/255;
o2LRcol=[144 187 226]/255;
o3LRcol=[168 168 168]/255;

%% Fig 3B - Time-resolved fraction of significant cells
for ss=1:3
if ss==1; cc2p=1:sum(nccs(1:4)); elseif ss==2; cc2p=(sum(nccs(1:4))+1):Nccs; else; cc2p=1:Nccs; end
hfig=figure('units','normalized','position',[0 0 1 1]);
subplot(3,1,1); hold on;
plotcmapdots([o1col; o2col]);
runlength(mean(lmELpv(cc2p,vtps_ds==1)<0.05,1),squeeze(mean(lmELpv_shf(cc2p,vtps_ds==1,:)<0.05,1)),1,rtmax_ds,o1col,.005); 
runlength(mean(lmERpv(cc2p,vtps_ds==1)<0.05,1),squeeze(mean(lmERpv_shf(cc2p,vtps_ds==1,:)<0.05,1)),1,rtmax_ds,o2col,.010); 
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),[0 0.2],'k:'); text(sum(rtmes_ds(1:jj-1)),.2,twLabels{jj}); end
hl=legend({'E(L)','E(R)'}); ylim([0 0.2]); box(hl,'on');

subplot(3,1,2);
plotcmapdots([o1LLcol; o1LRcol]);
runlength(mean(lmELpvLL(cc2p,vtps_ds==1)<0.05,1),squeeze(mean(lmELpvLL_shf(cc2p,vtps_ds==1,:)<0.05,1)),1,rtmax_ds,o1LLcol,.005); 
runlength(mean(lmELpvLR(cc2p,vtps_ds==1)<0.05,1),squeeze(mean(lmELpvLR_shf(cc2p,vtps_ds==1,:)<0.05,1)),1,rtmax_ds,o1LRcol,.010); 
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),[0 0.2],'k:'); text(sum(rtmes_ds(1:jj-1)),.2,twLabels{jj}); end
hl=legend({'E(L), Look L','E(L), Look R'}); box(hl,'on'); ylim([0 0.2]);

subplot(3,1,3);
plotcmapdots([o2LLcol; o2LRcol]);
runlength(mean(lmERpvLL(cc2p,vtps_ds==1)<0.05,1),squeeze(mean(lmERpvLL_shf(cc2p,vtps_ds==1,:)<0.05,1)),1,rtmax_ds,o2LLcol,.005); 
runlength(mean(lmERpvLR(cc2p,vtps_ds==1)<0.05,1),squeeze(mean(lmERpvLR_shf(cc2p,vtps_ds==1,:)<0.05,1)),1,rtmax_ds,o2LRcol,.010); 
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),[0 0.2],'k:'); text(sum(rtmes_ds(1:jj-1)),.2,twLabels{jj}); end
hl=legend({'E(R), Look L','E(R), Look R'}); box(hl,'on'); ylim([0 0.2]);

if ss<3; supertitle(['Subject ' num2str(ss)]); else; supertitle('Both Subjects'); end
saveas(hfig,[figures_folder '/perc_sig_cells_pv_s' num2str(ss) '.png']);
saveas(hfig,[figures_folder '/perc_sig_cells_pv_s' num2str(ss) '.svg']);

end

%% Fig 3C - Time-averaged fractions of significant cells
for ss=1:3
if ss==1; cc2p=1:sum(nccs(1:4)); elseif ss==2; cc2p=(sum(nccs(1:4))+1):Nccs; else; cc2p=1:Nccs; end
hfig=figure('units','normalized','position',[0 0 1 1]);
for tw=1:8
subplot(3,1,1); plotcmapdots([o1col; o2col]); ylim([0 .15]);
    nnm1=nanmean(lmELpv_shf(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1); 
    nnm2=nanmean(lmERpv_shf(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1);
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    alpha1=0.05+(mean(nanmean(lmELpv(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,2))>mean(prc951)).*.95;
    alpha2=0.05+(mean(nanmean(lmERpv(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,2))>mean(prc952)).*.95; 
    barmsem(tw-.125,nanmean(lmELpv(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),o1col,.1,'none',1,alpha1);
    barmsem(tw+.125,nanmean(lmERpv(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),o2col,.1,'none',1,alpha2); 
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051,2),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951,2),'k:');
    plot(tw+.125+[-.1 .1],[1 1].*mean(prc052,2),'k:'); plot(tw+.125+[-.1 .1],[1 1].*mean(prc952,2),'k:');
    text(tw-.1,0.14,getpstars(signrank(nanmean(lmELpv(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),...
                                       nanmean(lmERpv(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2))));
        
subplot(3,1,2); plotcmapdots([o1LLcol; o1LRcol]); ylim([0 .15]);
    nnm1=nanmean(lmELpvLL_shf(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05); 
    nnm2=nanmean(lmELpvLR_shf(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05);
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    alpha1=0.05+(mean(nanmean(lmELpvLL(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,2))>mean(prc951)).*.95;
    alpha2=0.05+(mean(nanmean(lmELpvLR(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,2))>mean(prc952)).*.95;  
    barmsem(tw-.125,nanmean(lmELpvLL(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),o1LLcol,.1,'none',1,alpha1);
    barmsem(tw+.125,nanmean(lmELpvLR(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),o1LRcol,.1,'none',1,alpha2);    
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051,2),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951,2),'k:');
    plot(tw+.125+[-.1 .1],[1 1].*mean(prc052,2),'k:'); plot(tw+.125+[-.1 .1],[1 1].*mean(prc952,2),'k:');
    text(tw-.1,0.14,getpstars(signrank(nanmean(lmELpvLL(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),...
                                       nanmean(lmELpvLR(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2))));
    
subplot(3,1,3); plotcmapdots([o2LLcol; o2LRcol]); ylim([0 .15]);
    nnm1=nanmean(lmERpvLL_shf(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05); 
    nnm2=nanmean(lmERpvLR_shf(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05);
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    alpha1=0.05+(mean(nanmean(lmERpvLL(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,2))>mean(prc951)).*.95;
    alpha2=0.05+(mean(nanmean(lmERpvLR(cc2p,ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,2))>mean(prc952)).*.95;  
    barmsem(tw-.125,nanmean(lmERpvLL(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),o2LLcol,.1,'none',1,alpha1);
    barmsem(tw+.125,nanmean(lmERpvLR(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),o2LRcol,.1,'none',1,alpha2);    
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051,2),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951,2),'k:');
    plot(tw+.125+[-.1 .1],[1 1].*mean(prc052,2),'k:'); plot(tw+.125+[-.1 .1],[1 1].*mean(prc952,2),'k:');
    text(tw-.1,0.14,getpstars(signrank(nanmean(lmERpvLL(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),...
                                       nanmean(lmERpvLR(cc2p,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2))));

end
if ss<3; supertitle(['Subject ' num2str(ss)]); else supertitle('Both subjects'); end
saveas(hfig,[figures_folder '/timewins_sig_in_time_pv_subj' num2str(ss) '.png']);
saveas(hfig,[figures_folder '/timewins_sig_in_time_pv_subj' num2str(ss) '.svg']);
end


%% Supp. Fig. 8A - R2 for the linear regression

cc2p=1:Nccs;
hfig=figure('units','normalized','position',[0 0 1 1]);
subplot(3,1,1); plotcmapdots([o1col; o2col]);
runlength(nanmean(lmELR2(cc2p,vtps_ds==1),1),squeeze(nanmean(lmELR2_shf(cc2p,vtps_ds==1,:),1)),1,rtmax_ds*10,o1col,.0002); 
runlength(nanmean(lmERR2(cc2p,vtps_ds==1),1),squeeze(nanmean(lmERR2_shf(cc2p,vtps_ds==1,:),1)),1,rtmax_ds*10,o2col,.0004); 
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj))*10,[0 0.01],'k:'); text(sum(rtmes_ds(1:jj-1))*10,.005,twLabels{jj}); end
hl=legend({'E(L)','E(R)'}); box(hl,'on'); ylim([0 0.005]);

subplot(3,1,2); plotcmapdots([o1LLcol; o1LRcol]);
runlength(nanmean(lmELR2LR(cc2p,vtps_ds==1),1),squeeze(nanmean(lmELR2LR_shf(cc2p,vtps_ds==1,:),1)),1,rtmax_ds*10,o1LRcol,.0002); 
runlength(nanmean(lmELR2LL(cc2p,vtps_ds==1),1),squeeze(nanmean(lmELR2LL_shf(cc2p,vtps_ds==1,:),1)),1,rtmax_ds*10,o1LLcol,.0004); 
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj))*10,[0 0.02],'k:'); text(sum(rtmes_ds(1:jj-1))*10,.02,twLabels{jj}); end
hl=legend({'E(L), Look L','E(L), Look R'}); box(hl,'on'); ylim([0 0.025]);

subplot(3,1,3); plotcmapdots([o2LLcol; o2LRcol]); 
runlength(nanmean(lmERR2LR(cc2p,vtps_ds==1),1),squeeze(nanmean(lmERR2LR_shf(cc2p,vtps_ds==1,:),1)),1,rtmax_ds*10,o2LRcol,.0002); 
runlength(nanmean(lmERR2LL(cc2p,vtps_ds==1),1),squeeze(nanmean(lmERR2LL_shf(cc2p,vtps_ds==1,:),1)),1,rtmax_ds*10,o2LLcol,.0004); 
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj))*10,[0 0.02],'k:'); text(sum(rtmes_ds(1:jj-1))*10,.02,twLabels{jj}); end
hl=legend({'E(R), Look L','E(R), Look R'}); box(hl,'on'); ylim([0 0.025]);
xlabel('task time (ms)');

saveas(hfig,[figures_folder '/R2_vs_R2shf_ELER_sig_in_time.png']);
saveas(hfig,[figures_folder '/R2_vs_R2shf_ELER_sig_in_time.svg']);