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');
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
% (yields same final results but all at once is faster)
% 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=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);

% 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


%% 
chRb1=nan(Ntps_ds,Nccs); chRb1_shf=nan(Ntps_ds,Nccs,Nshf);
chRb2=nan(Ntps_ds,Nccs); chRb2_shf=nan(Ntps_ds,Nccs,Nshf);
chRb3=nan(Ntps_ds,Nccs); chRb3_shf=nan(Ntps_ds,Nccs,Nshf);
chRp1=nan(Ntps_ds,Nccs); chRp1_shf=nan(Ntps_ds,Nccs,Nshf);
chRp2=nan(Ntps_ds,Nccs); chRp2_shf=nan(Ntps_ds,Nccs,Nshf);
chRp3=nan(Ntps_ds,Nccs); chRp3_shf=nan(Ntps_ds,Nccs,Nshf);
chRr=nan(Ntps_ds,Nccs);  chRr_shf=nan(Ntps_ds,Nccs,Nshf);
chRa=nan(Ntps_ds,Nccs,4);%chRa_shf=nan(Ntps_ds,Nccs,Nshf);

chRb1_LL=nan(Ntps_ds,Nccs); chRb1_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRb2_LL=nan(Ntps_ds,Nccs); chRb2_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRb3_LL=nan(Ntps_ds,Nccs); chRb3_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRp1_LL=nan(Ntps_ds,Nccs); chRp1_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRp2_LL=nan(Ntps_ds,Nccs); chRp2_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRp3_LL=nan(Ntps_ds,Nccs); chRp3_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRr_LL=nan(Ntps_ds,Nccs);  chRr_LL_shf=nan(Ntps_ds,Nccs,Nshf);
chRa_LL=nan(Ntps_ds,Nccs,4);%chRa_LL_shf=nan(Ntps_ds,Nccs,Nshf);

chRb1_LR=nan(Ntps_ds,Nccs); chRb1_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRb2_LR=nan(Ntps_ds,Nccs); chRb2_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRb3_LR=nan(Ntps_ds,Nccs); chRb3_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRp1_LR=nan(Ntps_ds,Nccs); chRp1_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRp2_LR=nan(Ntps_ds,Nccs); chRp2_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRp3_LR=nan(Ntps_ds,Nccs); chRp3_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRr_LR=nan(Ntps_ds,Nccs);  chRr_LR_shf=nan(Ntps_ds,Nccs,Nshf);
chRa_LR=nan(Ntps_ds,Nccs,4);%chRa_LR_shf=nan(Ntps_ds,Nccs,Nshf);

if recompute_mfit
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));
    ccoffchR=squeeze(poolofferchR(:,cc));

    parfor tt=1:Ntps_ds
    %for tt=1%
        warning('off','all')
        if vtps_ds(tt)==1
            
            ctccspk=squeeze(ccspk(:,tt));
            ctccepx=squeeze(ccepx(:,tt));
            
            toskip=isnan(ccoffchR) | isnan(ctccspk) | isnan(ctccepx);
            
            vctccspk=ctccspk(~toskip);
            vctccepx=ctccepx(~toskip);
            vccoffchR=ccoffchR(~toskip);
            vccoffEVR=ccoffEVR(~toskip);
            vccoffEVL=ccoffEVL(~toskip);
            nicc=length(vctccspk);
            
            isLL=find(vctccepx<0); niLL=length(isLL);
            isLR=find(vctccepx>0); niLR=length(isLR);
            
            %ctwcclmR=fitglm(vctccspk,vccoffchR,'distribution','binomial','link','logit');
            %ctwcclmR_LL=fitglm(vctccspk(isLL),vccoffchR(isLL),'distribution','binomial','link','logit');
            %ctwcclmR_LR=fitglm(vctccspk(isLR),vccoffchR(isLR),'distribution','binomial','link','logit');
            
            
            %ctwcclmR=fitglm([vccoffEVR vccoffEVL vctccspk],vccoffchR,'distribution','binomial','link','logit');
            %ctwcclmR_LL=fitglm([vccoffEVR(isLL) vccoffEVL(isLL) vctccspk(isLL)],vccoffchR(isLL),'distribution','binomial','link','logit');
            %ctwcclmR_LR=fitglm([vccoffEVR(isLR) vccoffEVL(isLR) vctccspk(isLR)],vccoffchR(isLR),'distribution','binomial','link','logit');
            
            ctwcclmR=fitlm([vccoffEVL vccoffEVR vccoffchR], vctccspk);
            ctwcclmR_LL=fitlm([vccoffEVL(isLL) vccoffEVR(isLL) vccoffchR(isLL)], vctccspk(isLL));
            ctwcclmR_LR=fitlm([vccoffEVL(isLR) vccoffEVR(isLR) vccoffchR(isLR)], vctccspk(isLR));
                        
            % Main results
            %ccoff1ev=squeeze(pooloffer1ev(:,cc));
            %ccoff2ev=squeeze(pooloffer2ev(:,cc));
            %ctwcclm1=fitlm(ccoff1ev,ctccspk);
            %ctwcclm2=fitlm(ccoff2ev,ctccspk);
            
            chRb1(tt,cc)=ctwcclmR.Coefficients.Estimate(2); chRb2(tt,cc)=ctwcclmR.Coefficients.Estimate(3); chRb3(tt,cc)=ctwcclmR.Coefficients.Estimate(4);
            chRp1(tt,cc)=ctwcclmR.Coefficients.pValue(2);   chRp2(tt,cc)=ctwcclmR.Coefficients.pValue(3);   chRp3(tt,cc)=ctwcclmR.Coefficients.pValue(4);
            chRr(tt,cc)=ctwcclmR.Rsquared.Ordinary;
            
            
            chRb1_LL(tt,cc)=ctwcclmR_LL.Coefficients.Estimate(2); chRb2_LL(tt,cc)=ctwcclmR_LL.Coefficients.Estimate(3); chRb3_LL(tt,cc)=ctwcclmR_LL.Coefficients.Estimate(4);
            chRp1_LL(tt,cc)=ctwcclmR_LL.Coefficients.pValue(2);   chRp2_LL(tt,cc)=ctwcclmR_LL.Coefficients.pValue(3);   chRp3_LL(tt,cc)=ctwcclmR_LL.Coefficients.pValue(4);
            chRr_LL(tt,cc)=ctwcclmR_LL.Rsquared.Ordinary;
            
            
            chRb1_LR(tt,cc)=ctwcclmR_LR.Coefficients.Estimate(2); chRb2_LR(tt,cc)=ctwcclmR_LR.Coefficients.Estimate(3); chRb3_LR(tt,cc)=ctwcclmR_LR.Coefficients.Estimate(4);
            chRp1_LR(tt,cc)=ctwcclmR_LR.Coefficients.pValue(2);   chRp2_LR(tt,cc)=ctwcclmR_LR.Coefficients.pValue(3);   chRp3_LR(tt,cc)=ctwcclmR_LR.Coefficients.pValue(4);
            chRr_LR(tt,cc)=ctwcclmR_LR.Rsquared.Ordinary;

            
            for sh=1:Nshf
                isshf=randperm(nicc); i1shf=randperm(nicc); i2shf=randperm(nicc); icshf=randperm(nicc); %trial-shuffle
                %ctwcclmR_shf=fitglm(vccoffchR(i2shf),vctccspk(isshf),'distribution','binomial','link','logit');
                ctwcclmR_shf=fitlm([vccoffEVL(i2shf) vccoffEVR(i2shf) vccoffchR(i2shf)],vctccspk(isshf));
                %ctwcclmR_shf=fitglm([vccoffEVR(i2shf) vccoffEVL(i1shf) vctccspk(icshf)],vccoffchR(isshf),'distribution','binomial','link','logit');
                chRp1_shf(tt,cc,sh)=ctwcclmR_shf.Coefficients.pValue(2);   chRp2_shf(tt,cc,sh)=ctwcclmR_shf.Coefficients.pValue(3);   chRp3_shf(tt,cc,sh)=ctwcclmR_shf.Coefficients.pValue(4);   
                chRb1_shf(tt,cc,sh)=ctwcclmR_shf.Coefficients.Estimate(2); chRb2_shf(tt,cc,sh)=ctwcclmR_shf.Coefficients.Estimate(3); chRb3_shf(tt,cc,sh)=ctwcclmR_shf.Coefficients.Estimate(4);
                chRr_shf(tt,cc,sh)=ctwcclmR_shf.Rsquared.Ordinary;
                
                %LL randperm
                isLLshf=isLL(randperm(niLL)); i1LLshf=isLL(randperm(niLL)); i2LLshf=isLL(randperm(niLL)); icLLshf=isLL(randperm(niLL)); %trial-shuffle
                %ctwcclmR_LL_shf=fitglm(vccoffchR(i2LLshf),vctccspk(isLLshf),'distribution','binomial','link','logit');
                ctwcclmR_LL_shf=fitlm([vccoffEVL(i2LLshf) vccoffEVR(i2LLshf) vccoffchR(i2LLshf)],vctccspk(isLLshf));
                %ctwcclmR_LL_shf=fitglm([vccoffEVR(i2LLshf) vccoffEVL(i1LLshf) vctccspk(icLLshf)],vccoffchR(isLLshf),'distribution','binomial','link','logit');
                chRp1_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Coefficients.pValue(2);   chRp2_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Coefficients.pValue(3);   chRp3_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Coefficients.pValue(4);
                chRb1_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Coefficients.Estimate(2); chRb2_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Coefficients.Estimate(3); chRb3_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Coefficients.Estimate(4);
                chRr_LL_shf(tt,cc,sh)=ctwcclmR_LL_shf.Rsquared.Ordinary;
                
                %LR randperm
                isLRshf=isLR(randperm(niLR)); i1LRshf=isLR(randperm(niLR)); i2LRshf=isLR(randperm(niLR)); icLRshf=isLR(randperm(niLR)); %trial-shuffle
                %ctwcclmR_LR_shf=fitglm(vccoffchR(i2LRshf),vctccspk(isLRshf),'distribution','binomial','link','logit');
                ctwcclmR_LR_shf=fitlm([vccoffEVL(i2LRshf) vccoffEVR(i2LRshf) vccoffchR(i2LRshf)],vctccspk(isLRshf));
                %ctwcclmR_LR_shf=fitglm([vccoffEVR(i2LRshf) vccoffEVL(i1LRshf) vctccspk(icLRshf)],vccoffchR(isLRshf),'distribution','binomial','link','logit');
                chRp1_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Coefficients.pValue(2);   chRp2_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Coefficients.pValue(3);   chRp3_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Coefficients.pValue(4);
                chRb1_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Coefficients.Estimate(2); chRb2_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Coefficients.Estimate(3); chRb3_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Coefficients.Estimate(4);
                chRr_LR_shf(tt,cc,sh)=ctwcclmR_LR_shf.Rsquared.Ordinary;
            end
            
        end
        warning('on','all')
    end
    fprintf('%s - done: %d/%d, %2.1f%%\n',datetime('now','Format','HH:mm'),cc,Nccs,cc/Nccs.*100);
end
    save(['fitlmstats_chR_evR_evL_spk_span_' num2str(Nspan) 'ms_lag_' num2str(Nlag) 'ms_' date '.mat'],'chR*');
else
    load('fitlmstats_chR_evR_evL_spk_span_200ms_lag_10ms_16-Feb-2023');
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=[55 130 200]/255;
o3LRcol=[168 168 168]/255;



%% Fig 4A - 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(4,1,1); hold on; gch=gca; plotcmapdots([o1col; o2col; o3col]);
runlength(nanmean(chRp1(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp1_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o1col,.005);
runlength(nanmean(chRp2(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp2_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o2col,.01);
runlength(nanmean(chRp3(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp3_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o3col,.015);
ylim([0 .25]); for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('All trials'); ylabel('% signif. cells'); xlabel('time bins (10ms)'); legend('EV_L','EV_R','chR');

subplot(4,1,2); hold on; gch=gca; plotcmapdots([o1LLcol; o1LRcol]);
runlength(nanmean(chRp1_LL(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp1_LL_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o1LLcol,.005);
runlength(nanmean(chRp1_LR(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp1_LR_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o1LRcol,.005);
ylim([0 .25]); for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('EV_L'); ylabel('% signif. cells'); xlabel('time bins (10ms)'); legend('LookL','LookR');

subplot(4,1,3); hold on; gch=gca; plotcmapdots([o2LLcol; o2LRcol]);
runlength(nanmean(chRp2_LL(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp2_LL_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o2LLcol,.01);
runlength(nanmean(chRp2_LR(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp2_LR_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o2LRcol,.01);
ylim([0 .25]); for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('EV_R'); ylabel('% signif. cells'); xlabel('time bins (10ms)'); legend('LookL','LookR');

subplot(4,1,4); hold on; gch=gca; plotcmapdots([o3LLcol; o3LRcol]);
runlength(nanmean(chRp3_LL(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp3_LL_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o3LLcol,.015);
runlength(nanmean(chRp3_LR(vtps_ds==1,cc2p)<0.05,2),squeeze(nanmean(chRp3_LR_shf(vtps_ds==1,cc2p,:)<0.05,2)),1,rtmax_ds,o3LRcol,.015);
ylim([0 .25]); for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('Choice'); ylabel('% signif. cells'); xlabel('time bins (10ms)'); legend('LookL','LookR');

if ss<3
    supertitle(['Subject ' num2str(ss)]);
    saveas(hfig,[figures_folder '/SuppFigS10_sig_in_time_pv_chR_subj' num2str(ss) '.png']);
    saveas(hfig,[figures_folder '/SuppFigS10_sig_in_time_pv_chR_subj' num2str(ss) '.svg']);
else
    supertitle(['Subjects 1:2']);
    saveas(hfig,[figures_folder '/Fig4A_sig_in_time_pv_chR.png']);
    saveas(hfig,[figures_folder '/Fig4A_sig_in_time_pv_chR.svg']);
end
end


%% Fig 4B - Fraction of significant cells in time epochs

hfig=figure('units','normalized','position',[0 0 1 1]);
for tw=1:8
subplot(4,1,1); plotcmapdots([o1col; o2col; o3col]);
    nnm1=nanmean(chRp1_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2); 
    nnm2=nanmean(chRp2_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2);
    nnm3=nanmean(chRp3_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2);
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    prc053=prctile(nnm3,05,3); prc953=prctile(nnm3,95,3);
    alpha1=0.1+(mean(nanmean(chRp1(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc951)).*.9;
    alpha2=0.1+(mean(nanmean(chRp2(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc952)).*.9;
    alpha3=0.1+(mean(nanmean(chRp3(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc953)).*.9;
    barmsem(tw-.125,nanmean(chRp1(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o1col,.1,'none',1,alpha1);
    barmsem(tw+.000,nanmean(chRp2(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o2col,.1,'none',1,alpha2);    
    barmsem(tw+.125,nanmean(chRp3(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o3col,.1,'none',1,alpha3);    
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051(:)),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951(:)),'k:');
    plot(tw+.000+[-.1 .1],[1 1].*mean(prc052(:)),'k:'); plot(tw+.000+[-.1 .1],[1 1].*mean(prc952(:)),'k:');
    plot(tw+.125+[-.1 .1],[1 1].*mean(prc053(:)),'k:'); plot(tw+.125+[-.1 .1],[1 1].*mean(prc953(:)),'k:');
        ylim([0 .2]); legend('EV_L','EV_R','choice'); title('All trials');
        %text(tw-.1,0.14,getpstars(signrank(nanmean(lmELpv(:,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),nanmean(lmERpv(:,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2))));
        
    subplot(4,1,2); plotcmapdots([o1LLcol; o1LRcol]);
    nnm1=nanmean(chRp1_LL_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2); 
    nnm2=nanmean(chRp1_LR_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2); 
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    alpha1=0.1+(mean(nanmean(chRp1_LL(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc951)).*.9;
    alpha2=0.1+(mean(nanmean(chRp1_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc952)).*.9;
    barmsem(tw-.125,nanmean(chRp1_LL(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o1LLcol,.1,'none',1,alpha1);
    barmsem(tw+.000,nanmean(chRp1_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o1LRcol,.1,'none',1,alpha2);
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051(:)),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951(:)),'k:');
    plot(tw+.000+[-.1 .1],[1 1].*mean(prc052(:)),'k:'); plot(tw+.000+[-.1 .1],[1 1].*mean(prc952(:)),'k:');
    ylim([0 .2]); legend('LookL','LookR'); title('Encoding of EV_L');
    %    text(tw-.1,0.14,getpstars(signrank(nanmean(lmELpvLL(:,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),nanmean(lmELpvLR(:,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2))));
    
    subplot(4,1,3); plotcmapdots([o2LLcol; o2LRcol]);
    nnm1=nanmean(chRp2_LL_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2);
    nnm2=nanmean(chRp2_LR_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2);
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    alpha1=0.1+(mean(nanmean(chRp2_LL(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc951)).*.9;
    alpha2=0.1+(mean(nanmean(chRp2_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc952)).*.9;
    barmsem(tw-.125,nanmean(chRp2_LL(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o2LLcol,.1,'none',1,alpha1);
    barmsem(tw+.000,nanmean(chRp2_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o2LRcol,.1,'none',1,alpha2);
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051(:)),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951(:)),'k:');
    plot(tw+.000+[-.1 .1],[1 1].*mean(prc052(:)),'k:'); plot(tw+.000+[-.1 .1],[1 1].*mean(prc952(:)),'k:');
    ylim([0 .2]); legend('LookL','LookR'); title('Encoding of EV_R');
    
    subplot(4,1,4); plotcmapdots([o3LLcol; o3LRcol]);
    nnm1=nanmean(chRp3_LL_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2);
    nnm2=nanmean(chRp3_LR_shf(ntce_ds(tw,1):ntce_ds(tw,2),:,:)<0.05,2);
    prc051=prctile(nnm1,05,3); prc951=prctile(nnm1,95,3);
    prc052=prctile(nnm2,05,3); prc952=prctile(nnm2,95,3);
    alpha1=0.1+(mean(nanmean(chRp3_LL(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc951)).*.9;
    alpha2=0.1+(mean(nanmean(chRp3_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1))>mean(prc952)).*.9;
    barmsem(tw-.125,nanmean(chRp3_LL(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o3LLcol,.1,'none',1,alpha1);
    barmsem(tw+.000,nanmean(chRp3_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1),o3LRcol,.1,'none',1,alpha2);
    plot(tw-.125+[-.1 .1],[1 1].*mean(prc051(:)),'k:'); plot(tw-.125+[-.1 .1],[1 1].*mean(prc951(:)),'k:');
    plot(tw+.000+[-.1 .1],[1 1].*mean(prc052(:)),'k:'); plot(tw+.000+[-.1 .1],[1 1].*mean(prc952(:)),'k:');
    ylim([0 .2]); legend('LookL','LookR'); title('Encoding of choice');
    
    %tt=text(tw-.125,0.2,getpstars(signrank(nanmean(chRp1_LR(ntce_ds(tw,1):ntce_ds(tw,2),:)<0.05,1)));
    %text(tw-.1,0.14,getpstars(signrank(nanmean(lmERpvLL(:,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2),nanmean(lmERpvLR(:,ntce_ds(tw,1):ntce_ds(tw,2))<0.05,2))));

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


%% Plotting R2

hfig=figure('units','normalized','position',[0 0 1 1]);
subplot(3,1,1); hold on; gch=gca; plotcmapdots([0 0 0; .33 .33 .33; .66 .66 .66]);
runlength(nanmean(chRr(vtps_ds==1,:),2),squeeze(nanmean(chRr_shf(vtps_ds==1,:,:),2)),1,rtmax_ds,[0 0 0] + 0.0,.005);
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('All trials'); ylabel('R^2'); xlabel('time bins (10ms)'); legend('R^2');

subplot(3,1,2); hold on; gch=gca; plotcmapdots([0 0 1; 0 0 .66; 0 0 .33]);
runlength(nanmean(chRr_LL(vtps_ds==1,:),2),squeeze(nanmean(chRr_LL_shf(vtps_ds==1,:,:),2)),1,rtmax_ds,[0 0 1].* 1.0,.005);
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('LookL'); ylabel('R^2'); xlabel('time bins (10ms)'); legend('R^2');

subplot(3,1,3); hold on; gch=gca; plotcmapdots([1 0 0; .66 0 0; .33 0 0]);
runlength(nanmean(chRr_LR(vtps_ds==1,:),2),squeeze(nanmean(chRr_LR_shf(vtps_ds==1,:,:),2)),1,rtmax_ds,[1 0 0].* 1.0,.005);
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),gch.YLim,'k:'); text(sum(rtmes_ds(1:jj-1)),1.02.*gch.YLim(2),twLabels{jj}); end
title('LookR'); ylabel('R^2'); xlabel('time bins (10ms)'); legend('R^2');

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