function [binary_nngt_comp,nnmd_diff_comp,nnmn_diff_comp,Zlengths] = plot_comparisons_choice_prob(...
    Zstruct_trwise,save_folder_path,filename_suffix)
addpath(genpath('../../../../my-tools'));

Nccs=size(Zstruct_trwise.Z_b1_LL_trwise,2);
sidestr={'Left','Right'};
casestr={'_b1_LL','_b1_LR','_b2_LL','_b2_LR'};
ttlab={'SV_L delay2 LookL','SV_L delay2 LookR','SV_R delay2 LookL','SV_R delay2 LookR'};

binary_nngt_comp=nan(4,Nccs); nnmd_diff_comp=nan(4,Nccs); nnmn_diff_comp=nan(4,Nccs);
Zlengths=nan(4,Nccs,4);
for cases=1:4
    Z1rnm=nanmean(eval(['Zstruct_trwise.Z' casestr{cases} '_trwise(:,:,:,1)']),3);
    Z2rnm=nanmean(eval(['Zstruct_trwise.Z' casestr{cases} '_trwise(:,:,:,2)']),3);
    Z3rnm=nanmean(eval(['Zstruct_trwise.Z' casestr{cases} '_trwise(:,:,:,3)']),3);
    Z4rnm=nanmean(eval(['Zstruct_trwise.Z' casestr{cases} '_trwise(:,:,:,4)']),3);   
    parfor cc=1:Nccs
        Z1rnm_cc=rmnans(Z1rnm(:,cc)); Z2rnm_cc=rmnans(Z2rnm(:,cc)); 
        Z3rnm_cc=rmnans(Z3rnm(:,cc)); Z4rnm_cc=rmnans(Z4rnm(:,cc));
        ZD11=nan; ZD21=nan; ZD31=nan;
        ZD12=nan; ZD22=nan; ZD32=nan;
        csw=(cases<=2)-(cases>2); % switches LookL>LookR in SVL to LookR>LookL in SVR
        Zlengths(cases,cc,:)=[length(Z1rnm_cc) length(Z2rnm_cc) length(Z3rnm_cc) length(Z4rnm_cc)];
        if ~isempty(Z1rnm_cc) && ~isempty(Z3rnm) % beta1 is positive: increasing p(chL), hence chL > chR
            [ZD11,ZD21,ZD31,~,ZD51]=vecdiffsubsamp(csw*Z1rnm_cc,csw.*Z3rnm_cc,100);
        end
        if ~isempty(Z2rnm_cc) && ~isempty(Z4rnm_cc) % beta1 is negative: increasing p(chR), hence chL < chR
            [ZD12,ZD22,ZD32,~,ZD52]=vecdiffsubsamp(csw.*Z4rnm_cc,csw.*Z2rnm_cc,100);
        end
        if cases==1 || cases==4
            binary_nngt_comp(cases,cc)=nanmax([ZD11 ZD12]);%nanmean(ZDiff>0,'all');
        else
            binary_nngt_comp(cases,cc)=nanmin([ZD11 ZD12]);%nanmean(ZDiff>0,'all');
        end
        nnmd_diff_comp(cases,cc)=nanmean([ZD21 ZD22]);%nanmedian(ZDiff,'all');
        nnmn_diff_comp(cases,cc)=nanmean([ZD31 ZD32]);%nanmedian(ZDiff,'all');
        %binary_nngt_comp2D(cases,cc,:)=[ZD11 ZD12];
        %clearvars ZD* Z1 Z2 Z3 Z4
    end
    fprintf('computing case:%d/4 %s\n',cases,ttlab{cases})
end

%% Quick safety checks
%save('FTPN.mat','TP','FP','TN','FN')
%FPR=nansum(FP,2)./(nansum(TN,2)+nansum(FP,2))
%TPR=nansum(TP,2)./(nansum(FN,2)+nansum(TP,2))
%nanmean(AUC1)
%nanmean(FPR1)
%nanmean(TPR1)
%plot(nanmean(FPR1),nanmean(TPR1),'k','linewidth',1.5); hold on
%plot([0 1],[0 1]); hold on
%error

hff=fullfig([0 0 .5 .5]); clo=colororder();
for cs=1:4
    subplot(2,2,cs); hold on
    if cs==2
        subplot(2,2,3);
    elseif cs==3
        subplot(2,2,2);
    end
    if cs==1 || cs==4; pvtail='right'; else pvtail='left'; end
    xvnm=(nnmn_diff_comp(cs,:)); pvnm=ranksum(xvnm,zeros(size(xvnm)),'tail',pvtail);
    xvnd=(nnmd_diff_comp(cs,:)); pvnd=ranksum(xvnd,zeros(size(xvnd)),'tail',pvtail);
    histogram(xvnm,-1e-3:2e-5:1e-3,'edgecolor','none','normalization','count'); hold on;
    histogram(xvnd,-1e-3:2e-5:1e-3,'edgecolor',[.5 .5 .5],'normalization','count','facecolor','none'); hold on; % ylim([0 20]);
    plot(nanmedian(xvnm),max(get(gca,'YLim')),'v','markeredgecolor',clo(1,:)); plotvline(0,'k-');
    plot(nanmedian(xvnd),max(get(gca,'YLim')),'v','markeredgecolor',[.5 .5 .5]); plotvline(0,'k-');
    text(min(get(gca,'XLim')).*.9,max(get(gca,'YLim')).*.9,sprintf('p=%.2i %s',pvnm,getpstars(pvnm)),'color',clo(1,:));
    text(min(get(gca,'XLim')).*.9,max(get(gca,'YLim')).*.8,sprintf('p=%.2i %s',pvnd,getpstars(pvnd)),'color',[.5 .5 .5]);
    xlabel('resudials difference'); ylabel('cell count'); box off;
    title(ttlab{cs});  legend({'mean','median'},'box','off')
end
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-4p' filename_suffix '.fig']);
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-4p' filename_suffix '.svg']);
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-4p' filename_suffix '.png']);

hff=fullfig([.5 .25 .5 .25]);
for cs=1:2
    subplot(1,2,cs); hold on
    xvnm=(nnmn_diff_comp(2*(cs-1)+1,:))-(nnmn_diff_comp(2*(cs-1)+2,:));
    if cs==2; xvnm=-xvnm; end
    pvnm=ranksum(xvnm,zeros(size(xvnm)),'tail','right');
    histogram(xvnm,-2e-3:5e-5:2e-3,'facecolor',clo(1,:),'edgecolor','none','normalization','count'); hold on; ylim([0 28]);
    plot(nanmedian(xvnm),max(get(gca,'YLim')),'rv','markeredgecolor',clo(1,:)); plotvline(0,'k-');

    xvnd=nnmd_diff_comp(2*(cs-1)+1,:)-nnmd_diff_comp(2*(cs-1)+2,:);
    if cs==2; xvnd=-xvnd; end
    pvnd=ranksum(xvnd,zeros(size(xvnd)),'tail','right');
    histogram(xvnd,-2e-3:5e-5:2e-3,'facecolor','none','edgecolor',[.3 .3 .3],'normalization','count'); hold on; ylim([0 28]);
    plot(nanmedian(xvnd),max(get(gca,'YLim')),'rv','markeredgecolor',[.3 .3 .3]); plotvline(0,'k-');

    text(-1.8e-3,max(get(gca,'YLim')).*.9,sprintf('p=%.2i %s',pvnm,getpstars(pvnm)),'color',clo(1,:));
    text(-1.8e-3,max(get(gca,'YLim')).*.8,sprintf('p=%.2i %s',pvnd,getpstars(pvnd)),'color',[.3 .3 .3]);
    xlabel('resudials difference'); ylabel('cell count');
    if cs==1
        title([ttlab{2*(cs-1)+1}(1:end-5) '(LookL - LookR)']);
    else
        title([ttlab{2*(cs-1)+1}(1:end-5) '(LookR - LookL)']);
    end
end
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-2p' filename_suffix '.fig']);
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-2p' filename_suffix '.svg']);
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-2p' filename_suffix '.png']);

hff=fullfig([.5 0 .5 .25]); comb=[1 4; 2 3];
for cs=1:2
    subplot(1,2,cs); hold on
    xvnm=[nnmn_diff_comp(comb(cs,1),:) nnmn_diff_comp(comb(cs,2),:)];
    %if cs==2; xvnm=-xvnm; end
    histogram(xvnm,-2e-3:2e-5:2e-3,'facecolor',[.5 .5 .5],'edgecolor','none','normalization','count'); hold on; ylim([0 38]);
    plot(nanmedian(xvnm),max(get(gca,'YLim')),'rv','markeredgecolor',[.5 .5 .5]); plotvline(0,'k-');

    xvnd=[nnmd_diff_comp(comb(cs,1),:) nnmd_diff_comp(comb(cs,2),:)];
    %if cs==2; xvnd=-xvnd; end
    histogram(xvnd,-2e-3:2e-5:2e-3,'facecolor',clo(1,:),'edgecolor','none','normalization','count'); hold on; ylim([0 38]);
    plot(nanmedian(xvnd),max(get(gca,'YLim')),'rv','markeredgecolor',clo(1,:)); plotvline(0,'k-');
    if cs==1
        pvnm=ranksum(xvnm,zeros(size(xvnm)),'tail','right');
        pvnd=ranksum(xvnd,zeros(size(xvnd)),'tail','right');
    else
        pvnm=ranksum(xvnm,zeros(size(xvnm)),'tail','left');
        pvnd=ranksum(xvnd,zeros(size(xvnd)),'tail','left');
    end
    xlim([-1.5e-3 1.5e-3]);
    text(min(get(gca,'XLim'))*.9,max(get(gca,'YLim')).*.9,sprintf('p=%.2i %s',pvnm,getpstars(pvnm)),'color',[.5 .5 .5]);
    text(min(get(gca,'XLim'))*.9,max(get(gca,'YLim')).*.8,sprintf('p=%.2i %s',pvnd,getpstars(pvnd)),'color',clo(1,:));
    xlabel('resudials difference (ch. ipsil. - ch. contral.)'); ylabel('cell count');
    if cs==1
        title('delay2 encoding (SV_L LookL, SV_R LookR)');
    else
        title('delay2 encoding (SV_L LookR, SV_R LookL)');
    end
end
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-2p-sides-' filename_suffix '.fig']);
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-2p-sides-' filename_suffix '.svg']);
saveas(hff,[save_folder_path 'ChoiceProb-SV-Look-2p-sides-' filename_suffix '.png']);


hff=fullfig([0 .5 .5 .5]); clo=colororder();  clo(1:2,:)=flipud(clo(1:2,:));
for cs=1:4
    subplot(2,2,cs);
    if cs==2; subplot(2,2,3); elseif cs==3; subplot(2,2,2); end
    xv=binary_nngt_comp(cs,:);
    if mod(cs,3)==1
        pv=ranksum(xv,.5*ones(size(xv)),'tail','right');
        text(0,14,sprintf('>.5 p=%.2i %s',pv,getpstars(pv)));
    else
        pv=ranksum(xv,.5*ones(size(xv)),'tail','left');
        text(0,14,sprintf('<.5 p=%.2i %s',pv,getpstars(pv)));
    end
    text(.25,25,sprintf('n=%2.0d',sum(xv<.5)),'color',clo(1,:)); hold on
    text(.75,25,sprintf('n=%2.0d',sum(xv>=.5)),'color',clo(2,:));
    histogram(xv(xv<.5),0:.01:1,'edgecolor','none','normalization','count','facecolor',clo(1,:)); hold on
    histogram(xv(xv>=.5),0:.01:1,'edgecolor','none','normalization','count','facecolor',clo(2,:)); hold on
    ylim([0 30]); plotvline(0.5,'k-'); ylim([0 30]);
    title(ttlab{cs}); xlabel([sidestr{ceil(cs/2)} ' choice probability']); ylabel('cell count'); box off
end
saveas(hff,[save_folder_path 'ChoicePref-SV-Look' filename_suffix '.fig']);
saveas(hff,[save_folder_path 'ChoicePref-SV-Look' filename_suffix '.svg']);
saveas(hff,[save_folder_path 'ChoicePref-SV-Look' filename_suffix '.png']);

hff=fullfig([.5 0 .5 .25]); clo=colororder();  clo(1:2,:)=flipud(clo(1:2,:));
for cs=1:2
    subplot(1,2,cs); hold on; plotcmapdots(clo);
    histogram(binary_nngt_comp(2*(cs-1)+1,:),0:7e-3:1,'edgecolor','none','normalization','count','facecolor',clo(1,:)); hold on
    histogram(binary_nngt_comp(2*(cs-1)+2,:),0:7e-3:1,'edgecolor','none','normalization','count','facecolor',clo(2,:)); hold on
    ylim([0 20]);
    nm1=nanmean(binary_nngt_comp(2*(cs-1)+1,:)); nsd1=nanstd(binary_nngt_comp(2*(cs-1)+1,:));
    nm2=nanmean(binary_nngt_comp(2*(cs-1)+2,:)); nsd2=nanstd(binary_nngt_comp(2*(cs-1)+2,:));
    plot(0:7e-3:1,Nccs*normpdf(0:7e-3:1,nm1,nsd1)*7e-3,'color',clo(1,:));
    plot(0:7e-3:1,Nccs*normpdf(0:7e-3:1,nm2,nsd2)*7e-3,'color',clo(2,:));
    plot(nanmedian(binary_nngt_comp(2*(cs-1)+1,:)),max(get(gca,'YLim')),'v','markeredgecolor','none','markerfacecolor',clo(1,:));
    plot([1 1].*nanmedian(binary_nngt_comp(2*(cs-1)+1,:)),get(gca,'YLim'),':','color',clo(1,:));
    plot(nanmedian(binary_nngt_comp(2*(cs-1)+2,:)),max(get(gca,'YLim')),'v','markeredgecolor','none','markerfacecolor',clo(2,:));
    plot([1 1].*nanmedian(binary_nngt_comp(2*(cs-1)+2,:)),get(gca,'YLim'),':','color',clo(2,:)); plotvline(.5,'k:')
    if cs==1
        pv=ranksum(binary_nngt_comp(2*(cs-1)+1,:),binary_nngt_comp(2*(cs-1)+2,:),'tail','right');
    else
        pv=ranksum(binary_nngt_comp(2*(cs-1)+1,:),binary_nngt_comp(2*(cs-1)+2,:),'tail','left');
    end
    %plotvline(0.5,'k-');
    text(0,14,sprintf('p=%.2i %s',pv,getpstars(pv)));
    title(ttlab{2*(cs-1)+1}(1:end-6)); xlabel([sidestr{cs} ' choice probability']); ylabel('cell count'); box off
    legend({'LookL','LookR'},'Location','northeast','box','off')
end
saveas(hff,[save_folder_path 'ChoicePref-SV-Look-comparison' filename_suffix '.fig']);
saveas(hff,[save_folder_path 'ChoicePref-SV-Look-comparison' filename_suffix '.svg']);
saveas(hff,[save_folder_path 'ChoicePref-SV-Look-comparison' filename_suffix '.png']);

hff=fullfig([.5 .25 .5 .25]); xorder=[1 2 4 3];
for cs=1:2
    subplot(1,2,cs);
    X1=binary_nngt_comp(xorder((cs-1)*2+1),:);
    X2=binary_nngt_comp(xorder((cs-1)*2+2),:);
    Y=0:.001:1;
    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);
    plot(FPR,TPR); hold on; plot([0 1],[0 1]); daspect([1 1 1]);
    xlabel('False positive rate'); ylabel('True positive rate');
    text(.3,.1,sprintf('AUC = %.3f',max(cumtrapz(flip(FPR),flip(TPR)))));
    title(['ROC ch=' sidestr{cs}(1) ', ' ttlab{xorder((cs-1)*2+1)} ' vs ' ttlab{xorder((cs-1)*2+2)}(end-4:end)]);
end
saveas(hff,[save_folder_path 'ChoicePref-SV-Look-comparison-ROC' filename_suffix '.fig']);
saveas(hff,[save_folder_path 'ChoicePref-SV-Look-comparison-ROC' filename_suffix '.svg']);
saveas(hff,[save_folder_path 'ChoicePref-SV-Look-comparison-ROC' filename_suffix '.png']);

end