% localization example 1, step 2:
% this script compiles data generated by calls to step1_run.m, and creates
% some plots.

% Requirements: 
% - matlab R2015a or later, with 
%   image_toolbox
%   optimization_toolbox
%   statistics_toolbox
% - uncertainSPT, with all folders and subfolders on your matlab path. Try,
%   e.g. addpath(genpath(fullfile(pwd,'..','..','..','uncertainSPT'))); 

if(1) % collect data
    clear    
    % file
    resFolder='EX1results';
    
    sgStr={'MLE','MAP'};
    agStr={'MLE','MAP'};
    % asymmetric Gauss MLE
    % hard limit on S1,S2, otherwise MLE fit
    % hard limit on S1,S2, lnB~N(0,ln30), (dS1+dS2)/2/S0~exp(1), ln(S1/S2)~N(0,log(3))
    % hard limit on S1,S2, lnB~N(0,ln30), (dS1+dS2)/2/S0~exp(1), ln(S1/S2)~N(0,log(2))
    % hard limit on S1,S2, lnB~N(0,ln30), (dS1+dS2)/2/S0~exp(1), ln(S1/S2)~N(0,log(1.5))
    
    psfStr={'514nm NA1.49','639nm NA1.4','680nm NA1.49'};    
    Rfiles=dir(fullfile(resFolder,'*RW*.mat'));
    NR=numel(Rfiles);
    
    vars={'Nphotons','bg','tE','dt','psfString','psfIndex','emTrj_px',...
        'SMopt','sigMin'...
        'NA','lambda',...
        'PSFsymGauss', 'SG_coord','SG_bNSdr',   'SG_cov','SG_CRLB',...
        'PSFasymGauss','AG_coord','AG_bNS12vdr','AG_cov'};
    disp('              ')
    for k=1:numel(Rfiles)
        Rk=load(fullfile(resFolder,Rfiles(k).name),vars{:});
        
        % compute error and normalized point properties, symGauss
        for m=1:numel(Rk.SG_coord)
            Rk.SG_xyErr{m}=Rk.SG_coord{m}-Rk.emTrj_px(:,1:2);
            Rk.bn{m}=Rk.SG_bNSdr{m}(:,1)/Rk.bg;
            Rk.Nn{m}=Rk.SG_bNSdr{m}(:,2)/Rk.Nphotons;
            Rk.Sn{m}=Rk.SG_bNSdr{m}(:,3)/Rk.sigMin;
        end
        % compute error and normalized point properties, asymGauss
        for m=1:numel(Rk.AG_coord)
            Rk.AG_xyErr{m}=Rk.AG_coord{m}-Rk.emTrj_px(:,1:2);
            Rk.abn{m}=Rk.AG_bNS12vdr{m}(:,1)/Rk.bg;
            Rk.aNn{m}=Rk.AG_bNS12vdr{m}(:,2)/Rk.Nphotons;
            Rk.aSn{m}=Rk.AG_bNS12vdr{m}(:,3:4)/Rk.sigMin;
        end
        
        R(k)=Rk;
        fprintf('\b\b\b\b\b\b\b\b\b\b%4d %5d',k,numel(Rfiles))
    end
    fprintf('\n')
    
    vars=fieldnames(R(1));
    % selection variables
    S=struct;
    S.Nphotons=[R(:).Nphotons];
    S.bg=[R(:).bg];
    S.tE=[R(:).tE];
    S.psfIndex=[R(:).psfIndex];
    S.lambda=[R(:).lambda];
    S.NA=[R(:).NA];
    S.sigMin=[R(:).sigMin];
    
    % label variables
    L=S;
    fa=fieldnames(S);
    for k=1:numel(fa)
        L.(fa{k})=sort(unique(S.(fa{k})));
    end
    % global parameters and trucated std curve
    px2nm=80;
    nbst=2; % overall bootstrap level. >100 gives decent estimates, , set 3 to maximize execution
    
    % calibration curve for truncated RMS estimates
    dx=randn(1,1e6);
    dy=randn(1,1e6);
    drMax=[0.01 0.1:0.1:5 10 20 10000];
    dxStd=0*drMax;
    for k=1:numel(drMax)
        dxStd(k)=std(dx( sqrt(dx.^2+dy.^2)<=drMax(k)));
    end
    trStdFun=@(std0,rMax)(std0.*interp1(drMax,dxStd,rMax./std0,'linear'));
    clear dx dy drMax dxStd
end
% single-panel production plots
if(1)
    figVmax=150;
    ncols=1
    nrows=1
    nbst=300
    % fig 2:
     plotCSA=[1 1 1];
    % fig 3:
    %plotCSA=[0 1 1];
    
    % estimated RMS bins
    Lrms=150;
    NBmin=300; % minimum number of dots/bin
    dRMS=7.5;
    RMSedge=dRMS/2:dRMS:Lrms;
    Lxy=80;

    %sLaplCol=[0 0.5 1];
    sLaplCol='m';
    markSize=3;
    
    %fig 2a,b
    fs=1;fa=1;
    % fig 2e,f, fig 3
    fs=2;fa=2;

    kr=1;
    sMin=0.0;
    sMax=inf;
    sMaxPX=9;
    zMin=00;
    zMax=inf;
    drMax=4;
    varMax=16^2;
    
    r_kE=1:numel(L.tE)    
    r_kN=1:numel(L.Nphotons)
    r_bg=1:numel(L.bg)
    r_kP=2;
    
    titleFun=@(kr)([]);%sprintf(agStr{r_fa(kr)}));%sprintf('MLE fit, s_{PSF} > %.1f s_{min}',r_sMin(kr)));
    posFun=@(ax,kr)(set(ax,'position',[0.35 0.24 0.63 0.75]));
    %% collect data
    clear c0Frac czFrac xErr xLap xCRLB rLap rCRLB dxErr dxLap dxCRLB drLap drCRLB fKeptSG dfKeptSG
    clear c0Frac czFrac xaErr xaLap raLap dxaErr dxaLap draLap fKeptAG dfKeptAG
    clear xAsErr dxAsErr xaSErr dxaSErr
    clear cc1 ccN
    clear xcCRLB ycCRLB yncCRLB xcSL ycSL yncSL xcAL ycAL yncAL
    
    sBNS=[];aBNS12=[];
    snBNS=[];anBNS12=[];
    
    np=0; % parameter counter
    ns=0;
    na=0;
    nsa=0;
    
    neSG =[]; % [xErr xLap xCRLB ]
    neAG =[]; % [xErr xLap]
    
    for kN=r_kN
        for kb=r_bg
            for kE=r_kE
                for kP=r_kP
                    Aind=S.Nphotons==L.Nphotons(kN) & S.bg==L.bg(kb) ...
                        & S.tE==L.tE(kE) & S.psfIndex==kP;
                    if(sum(Aind)>0)
                        A=vecAggregate(R,Aind,[],...
                            {'emTrj_px',...
                            'SG_bNSdr','SG_cov','SG_CRLB','SG_xyErr','Sn','bn','Nn',...
                            'AG_bNS12vdr','AG_cov','AG_xyErr','aSn','abn','aNn'});
                        
                        np=np+1;
                        
                        % c0 : convergence with second derivative, inside ROI mid-circle
                        c0SG=isfinite(A.SG_cov{fs}(:,1)) & A.SG_cov{fs}(:,3) >0;
                        c0AG=isfinite(A.AG_cov{fa}(:,1)) & A.AG_cov{fa}(:,3) >0;
                        
                        % cE : upper bounds on CRLB & LaplCov
                        cUSG =   A.SG_cov{fs}(:,1)<varMax & A.SG_cov{fs}(:,3)<varMax;% & A.SG_CRLB{fs}(:,1)<varMax;
                        cUAG =	 A.AG_cov{fa}(:,1)<varMax & A.AG_cov{fa}(:,3)<varMax;
                        
                        % range of PSF widths
                        cSsg =  A.Sn{fs}(:,1)>= sMin & A.Sn{fs}(:,1)<= sMax & A.SG_bNSdr{fs}(:,3) <= sMaxPX;
                        cSag =  A.aSn{fa}(:,1)>=sMin & A.aSn{fa}(:,2)>=sMin & min(A.AG_bNS12vdr{fa}(:,3:4),[],2) <= sMaxPX ...
                            & A.aSn{fa}(:,1)<=sMax & A.aSn{fa}(:,2)<=sMax ;
                        
                        % maximum displacement
                        cRs = A.SG_bNSdr{fs}(:,4) < drMax;
                        cRa = A.AG_bNS12vdr{fa}(:,6) < drMax;
                        
                        cz = abs(A.emTrj_px(:,3)*px2nm) <= zMax & zMin <= abs(A.emTrj_px(:,3)*px2nm);
                        
                        % symmetric Gaussian
                        ccSG= c0SG & cRs & cUSG & cz & cSsg;
                        % asymmetric Gaussian
                        ccAG= c0AG & cRa & cUAG & cz & cSag;
                        
                        % compare them ?
                        ccSA=ccSG & ccAG;
                        
                        fprintf('%s, N=%d ph., bg=%d ph./px, tE=%.1fms, %s\n',...
                            sgStr{fs},L.Nphotons(kN),L.bg(kb),L.tE(kE)*1e3,psfStr{kP})
                        fprintf(' symGauss: kept %d of %d, fraction %6.3f\n',sum(ccSG),numel(ccSG),sum(ccSG)/numel(ccSG))
                        fprintf('asymGauss: kept %d of %d, fraction %6.3f\n',sum(ccAG),numel(ccAG),sum(ccAG)/numel(ccAG))
                        disp('---')
                        
                        % symm Gauss [ error vLaplace vCRLB]
                        neSG0=[A.SG_xyErr{fs}(ccSG,1)*px2nm A.SG_cov{fs}(ccSG,1)*px2nm^2 A.SG_CRLB{fs}(ccSG,1)*px2nm^2];
                        neSG=[neSG; neSG0];
                        
                        % asymm Gauss [ error vLaplace]
                        neAG0=[A.AG_xyErr{fa}(ccAG,1)*px2nm A.AG_cov{fa}(ccAG,1)*px2nm^2 ];
                        neAG=[neAG; neAG0];
                        
                        % conditional errors in each data set
                        [ycCRLB(np,:),xcCRLB(np,:),~,~]= condMeanHist(neSG0(:,3),neSG0(:,1),(RMSedge).^2,@std,NBmin);
                        [ycSL(np,:),  xcSL(np,:),  ~,~]= condMeanHist(neSG0(:,2),neSG0(:,1),(RMSedge).^2,@std,NBmin);
                        [ycAL(np,:),  xcAL(np,:),  ~,~]= condMeanHist(neAG0(:,2),neAG0(:,1),(RMSedge).^2,@std,NBmin);
                        
                        % normalized conditional errors in each data set
                        [yncCRLB(np,:),xcCRLB(np,:),~,~]= condMeanHist(neSG0(:,3),neSG0(:,1)./sqrt(neSG0(:,3)),(RMSedge-0*dRMS/5).^2,@std,NBmin);
                        [yncSL(np,:),  xcSL(np,:),  ~,~]= condMeanHist(neSG0(:,2),neSG0(:,1)./sqrt(neSG0(:,2)),(RMSedge+0*dRMS/5).^2,@std,NBmin);
                        [yncAL(np,:),  xcAL(np,:),  ~,~]= condMeanHist(neAG0(:,2),neAG0(:,1)./sqrt(neAG0(:,2)),(RMSedge-0*dRMS/5).^2,@std,NBmin);
                        
                        % averages over this parameter set
                        if(sum(ccSG>0)>100)
                            ns=ns+1;
                            % Laplace err , symm Gauss
                            x=A.SG_xyErr{fs}(ccSG,1)*px2nm;
                            xErr(ns) =std(x);
                            dxErr(ns)=std(bootstrp(nbst,@std,x));
                            
                            x=A.SG_cov{fs}(ccSG,1)*px2nm^2;
                            xLap(ns) =sqrt(mean(x));
                            dxLap(ns)=std(bootstrp(nbst,@(y)(sqrt(mean(y))),x));
                            
                            x=[A.SG_cov{fs}(ccSG,1)*px2nm^2 A.SG_xyErr{fs}(ccSG,1)*px2nm];
                            rLap(ns)=std(x(:,2))/sqrt(mean(x(:,1)));
                            drLap(ns)=std(bootstrp(nbst,@(y)(std(y(:,2))/sqrt(mean(y(:,1)))),x));
                            
                            % CRLB, sym Gauss
                            x=A.SG_CRLB{fs}(ccSG,1)*px2nm^2;
                            xCRLB(ns) =sqrt(mean(x));
                            dxCRLB(ns)=std(bootstrp(nbst,@(y)(sqrt(mean(y))),x));
                            
                            x=[A.SG_CRLB{fs}(ccSG,1)*px2nm^2 A.SG_xyErr{fs}(ccSG,1)*px2nm];
                            rCRLB(ns)=std(x(:,2))/sqrt(mean(x(:,1)));
                            drCRLB(ns)=std(bootstrp(nbst,@(y)(sqrt(mean(y(:,1)))/std(y(:,2))),x));
                            
                            x=ccSG;
                            fKeptSG(ns) =mean(x);
                            dfKeptSG(ns)=std(bootstrp(nbst,@mean,x));
                            
                            cc1SG(ns,fa)=sum(ccSG==1);
                            ccNSG(ns,fa)=numel(ccSG);
                            
                            % dot parameters
                            sBNS=[sBNS ; A.SG_bNSdr{fs}(ccSG,1:3)];

                            % normalized dot parameters
                            snBNS=[snBNS ; [A.bn{fs}(ccSG,1) A.Nn{fs}(ccSG,1) A.Sn{fs}(ccSG,1)]];
                        end
                        if(sum(ccAG>0)>100)
                            na=na+1;
                            % Laplace err, asym Gauss
                            x=A.AG_xyErr{fa}(ccAG,1)*px2nm;
                            xaErr(na) =std(x);
                            dxaErr(na)=std(bootstrp(nbst,@std,x));
                            
                            x=A.AG_cov{fa}(ccAG,1)*px2nm^2;
                            xaLap(na) =sqrt(mean(x));
                            dxaLap(na)=std(bootstrp(nbst,@(y)(sqrt(mean(y))),x));
                            
                            x=[A.AG_cov{fa}(ccAG,1)*px2nm^2 A.AG_xyErr{fa}(ccAG,1)*px2nm];
                            raLap(na)=std(x(:,2))/sqrt(mean(x(:,1)));
                            draLap(na)=std(bootstrp(nbst,@(y)(std(y(:,2))/sqrt(mean(y(:,1)))),x));
                            
                            x=ccAG;
                            fKeptAG(na) =mean(x);
                            dfKeptAG(na)=std(bootstrp(nbst,@mean,x));
                            
                            cc1AG(na,fa)=sum(ccAG==1);
                            ccNAG(na,fa)=numel(ccAG);
                            
                            % dot parameters
                            aBNS12=[aBNS12 ; (A.AG_bNS12vdr{fa}(ccSG,1:4))];

                            % normalized dot parameters
                            anBNS12=[anBNS12 ; [A.abn{fa}(ccAG,1) A.aNn{fa}(ccAG,1) A.aSn{fa}(ccAG,1:2)]];
                            
                        end
                        if(false && sum(ccSA>0)>100)
                            nsa=nsa+1;
                            % Laplace err, sym Gauss
                            x=A.SG_xyErr{fs}(ccSA,1)*px2nm;
                            xaSErr(nsa) =std(x);
                            dxaSErr(nsa)=std(bootstrp(nbst,@std,x));
                            
                            % Laplace err, asym Gauss
                            x=A.AG_xyErr{fa}(ccSA,1)*px2nm;
                            xAsErr(nsa) =std(x);
                            dxAsErr(nsa)=std(bootstrp(nbst,@std,x));
                            
                        end
                    end
                end
            end
        end
    end
    %% plot
    % MLboxplot of normalized conditional RMSD|est RMSD
    if(1)
        figure(figVmax+15)
        subplot(nrows,ncols,kr)
        hold off
        
        % conditional truncated std
        std0=linspace(0,2*Lxy,300); % nm
        plot(std0,trStdFun(std0,drMax*px2nm)./std0,'k--')
        hold on
        plot(std0,0.9*trStdFun(std0,drMax*px2nm)./std0,'k:')
        plot(std0,1.1*trStdFun(std0,drMax*px2nm)./std0,'k:')
        hp=[];
        leg={};
        
        % box plots
        if(plotCSA(1)) % RMSD | CRLB, averaged over all points
            leg{end+1}='CRLB';
            XC=mean(xcCRLB,1,'omitnan'); % predicted precision (variance)
            boxplot(yncCRLB,'position',sqrt(XC)-dRMS/5,'plotStyle','traditional',...
                'colors','k','symbol','k+','outliersize',3,'boxstyle','filled',...
                'whisker',inf)
            %YM=median(yncCRLB,'omitnan');
            %plot(sqrt(XC),YM,'ks','markerface','w','markersize',markSize)
        end
        if(plotCSA(2)) % RMSD | sym. Laplace
            leg{end+1}='Lapl.';
            XC=mean(xcSL,1,'omitnan');
            boxplot(yncSL,'position',(sqrt(XC)+dRMS/5),'plotStyle','traditional',...
                'colors',sLaplCol,'symbol','c+','outliersize',3,'boxstyle','filled',...
                'whisker',inf)
            %YM=median(yncSL,'omitnan');
            %plot(sqrt(XC),YM,'ro','markerface','w','markersize',markSize)
        end
        if(plotCSA(3)) % RMSD | asym. Laplace
            leg{end+1}='asym.L.';
            XC=mean(xcAL,1,'omitnan');
            boxplot(yncAL,'position',(sqrt(XC)-dRMS/5),'plotStyle','traditional',...
                'colors','b','symbol','k+','outliersize',3,'boxstyle','filled',...
                'whisker',inf)
            %YM=median(yncAL,'omitnan');
            %plot(sqrt(XC),YM,'b^','markerface','w','markersize',markSize)
        end
        % medians on top
        if(plotCSA(1)) % RMSD | CRLB, averaged over all points
            XC=mean(xcCRLB,1,'omitnan'); % predicted precision (variance)
            YM=median(yncCRLB,'omitnan');
            plot(sqrt(XC)-dRMS/5,YM,'ks','markerface','w','markersize',markSize)
        end
        if(plotCSA(2)) % RMSD | sym. Laplace
            XC=mean(xcSL,1,'omitnan');
            YM=median(yncSL,'omitnan');
            plot(sqrt(XC)+dRMS/5,YM,'co','markerface','w','markersize',markSize,'color',sLaplCol)
        end
        if(plotCSA(3)) % RMSD | asym. Laplace
            XC=mean(xcAL,1,'omitnan');
            YM=median(yncAL,'omitnan');
            plot(sqrt(XC)-dRMS/5,YM,'b^','markerface','w','markersize',markSize)
        end
        
        % a spurious boxplot to get reasonable x ticks
        xTick=0:20:Lrms;
        xTickData=-1e6+rand(3,length(xTick));
        boxplot(xTickData,'position',xTick,'Labels',num2str(xTick'))
        
        
        xlabel('est.RMSE [nm]')
        ylabel('true/est. cRMSE [nm]')
        grid off
        set(gca,'xscale','lin','yscale','lin')
        % fig 2b
        axis([dRMS/2 Lxy+9 0.85 2.35])
        % fig 2f, 3b
        %axis([dRMS/2 Lxy+9 0.85 1.35])
        posFun(gca,kr);
        title(titleFun(kr))
        %set(gca,'xtick',get(gca,'ytick'))
        pause(0.1)
    end
    % absolute errors
    if(1)
        figure(figVmax+5)
        subplot(nrows,ncols,kr)
        hold off
        
        hp=[];
        leg={};
        if(plotCSA(1))
            hp(end+1)=errorbar(xCRLB,xErr,dxErr,dxErr,dxCRLB,dxCRLB,'ks','markersize',markSize,'markerface','k','capsize',0);
            leg{end+1}='CRLB';
            hold on
        end
        if(plotCSA(2))
            hp(end+1)=errorbar(xLap,xErr,dxErr,dxErr,dxLap,dxLap,'mo','markersize',markSize,'markerface','w','capsize',0,'color',sLaplCol);
            leg{end+1}='Lapl.';
            hold on
        end
        if(plotCSA(3))
            hp(end+1)=errorbar(xaLap,xaErr,dxaErr,dxaErr,dxaLap,dxaLap,'b^','markersize',markSize,'markerface','b','capsize',0);
            leg{end+1}='asym.L.';
            hold on
        end
        std0=linspace(1e-3,Lxy,300); % nm
        plot(std0,trStdFun(std0,drMax*px2nm),'k--')
        plot(std0,1.1*trStdFun(std0,drMax*px2nm),'k:')
        plot(std0,0.9*trStdFun(std0,drMax*px2nm),'k:')
        
        set(gca,'xscale','log','yscale','log')
        %axis equal
        axis([10 Lxy 10 Lxy])
        legend(hp,leg,'location','best')
        title(titleFun(kr))
        posFun(gca,kr);
        xlabel('est. RMSE [nm]')
        ylabel('true RMSE [nm]')
        box on
        grid off
        
        set(gca,'ytick',[15 30 60])
        set(gca,'xtick',get(gca,'ytick'))
        set(gca,'ytick',get(gca,'xtick'))
    end
    % Gaussian probability plots
    if(1)
        std0=0;
        std1=3*80;
        figure(figVmax+3)
        subplot(nrows,ncols,kr)
        hold off
        dXY=5;
        hp=[];
                
        % CRLB
        if(false && plotCSA(1))
            ind= std0 <= neSG(:,3) & neSG(:,3) <= std1;
            [x,y,Q,yTick,qTick]=QQplot_N01(neSG(ind,1)./sqrt(neSG(ind,3)));
            hp(end+1)=plot(x,y,'k-');
            hold on
        end
        % sym. Laplace
        if(true && plotCSA(2))
            ind= std0 <= neSG(:,2) & neSG(:,2) <= std1;
            [x,y,Q,yTick,qTick]=QQplot_N01(neSG(ind,1)./sqrt(neSG(ind,2)));
            hp(end+1)=plot(x,y,'m-');
            hold on
        end
        if(plotCSA(3))
            % asymmetric Laplace
            ind= std0 <= neAG(:,2) & neAG(:,2) <= std1;
            [x,y,Q,yTick,qTick]=QQplot_N01(neAG(ind,1)./sqrt(neAG(ind,2)));
            hp(end+1)=plot(x,y,'r-');
            hold on
        end
        set(gca,'ytick',yTick,'yticklabel',num2str(qTick',4),'xlim',[-5 5])
        
        % truncated Gaussian from asymmetric Laplace std
        % simulate truncated Gaussian distribution with same mixture of stds
        yStd=sqrt(neSG(ind,2));
        yErr=zeros(size(yStd));
        ind=(1:numel(yErr))';
        it=0;
        while(~isempty(ind))
            dx=yStd(ind).*randn(size(ind));
            dy=yStd(ind).*randn(size(ind));
            i2= (sqrt(dx.^2+dy.^2)<drMax*px2nm);
            yErr(ind(i2))=dx(i2);
            ind=find(~i2);
            it=it+1
        end
        [x,y,Q,yTick,qTick]=QQplot_N01(yErr./yStd);
        hp(end+1)=plot(x,y,'k--');
        %set(gca,'ytick',ytick','yticklabel',num2str(qTick',2))
        
        legend(hp,['sCRLB ' sgStr{fs}],['sLapl. ' sgStr{fs}],['aLapl. ' agStr{fa}],'location','best')
        %legend(hp,'CRLB','Lapl.','asym.L.','location','best')
        %legend(hp,'CRLB','Laplace','asym. Lapl.','trunc. Gauss','location','southeast')
        posFun(gca,kr);
        title(titleFun(kr))
        xlabel('true/est. RMS err.')
        ylabel('probability')
        axis equal
        axis(dXY*[-1 1 -1 1])
        set(gca,'position',[0.35 0.3 0.61 0.68])
        grid on
        box on
    end
    % errors of symmetric vs asymmetric Gauss PSF
    if(1)
        figure(figVmax+6)
        subplot(nrows,ncols,kr)
        hold off
        plot(2*[0 Lxy],2*[0 Lxy],'k--')
        hold on
        
        errorbar(xErr,xaErr,dxaErr,dxaErr,dxErr,dxErr,'k.','capsize',0)
        
        xlabel('true RMSE sym  [nm]')
        ylabel('true RMSE asym [nm]')
        legend(['asym: ' agStr{fa}])
        
        box on
        grid on
        title(titleFun(kr))
        posFun(gca,kr);
        axis(1.2*Lxy*[0 1 0 1])
        
        keptFractionSG=sum(cc1SG,1)./sum(ccNSG,1)
        %keptFractionAG=sum(cc1AG,1)./sum(ccNAG,1)
    end
    % dot convergence and retention frequencies
    if(1)
        figure(figVmax+7)
        subplot(nrows,ncols,kr)
        hold off
        hp=[];
        set(plot(xErr,fKeptSG,'k.'),'marker','.','linestyle','none')
        if(plotCSA(3))
            hold on
            set(plot(xaErr,fKeptAG,'r.'),'marker','+','linestyle','none')
        end
        box on
        grid on
        xlabel('true RMSE [nm]')
        ylabel('fraction kept dots')
        legend(hp,['sLapl. ' sgStr{fs}],['aLapl. ' agStr{fa}],'location','best')
        %legend(hp,'sym.','asym.','location','best')
        %legend(hp,'symmetric Gauss','asymmetric Gauss','location','best')
        %legend(hp,'symmetric Gauss','location','best')
        axis([0 Lxy 0.0 1.05])
        title(titleFun(kr))
        posFun(gca,kr);
        set(gca,'xtick',0:20:100)
    end
    % PSF width distributions
    if(1)
        figure(figVmax+8)
        subplot(nrows,ncols,kr)
        hold off
        hp=[];
        
        S0=S.sigMin(kr)*px2nm;%0.21*639/1.4;
        dS=5;
        sBin=dS/2:dS:1000;
        [a,b]=hist(sBNS(:,3)*px2nm,sBin);
        a=a([1:end end])/sum(a)/dS;
        b=[b-dS/2 b(end)+dS/2];

        hp=plot(S0*[1 1],[0 5*max(a)],'k--');
        hold on        

        ss=linspace(0,5*S0,1000);
        yy=1/S0*exp(-(ss-S0)/S0);
        yy(ss<S0)=0;
        hp(end+1)=plot(ss,yy,'-c');
        
        hp(end+1)=stairs(b,a,'k-');

        [a,b]=hist(min(aBNS12(:,3:4),[],2)*px2nm,sBin);
        a=a([1:end end])/sum(a)/dS;
        b=[b-dS/2 b(end)+dS/2];
        hp(end+1)=stairs(b,a,'b-');
        
        [a,b]=hist(max(aBNS12(:,3:4),[],2)*px2nm,sBin);
        a=a([1:end end])/sum(a)/dS;
        b=[b-dS/2 b(end)+dS/2];
        hp(end+1)=stairs(b,a,'r-');
        legend(hp,'0.21\lambda/NA','prior','sym.','asym.min','asym.max','location','best')        
        
        axis([0 390 0 0.05])
        xlabel('PSF width [nm]')
        ylabel('freq.')
    end
    % normalized PSF amplitude distribution
    if(1)
        figure(figVmax+9)
        subplot(nrows,ncols,kr)
        hold off
        hp=[];
        leg={};
        
        nBin=linspace(0,4,301);
        dn=mean(diff(nBin));
        [a,b]=hist(snBNS(:,2),nBin);
        a=a/sum(a)/dn;
        
        hp(end+1)=plot(b,a,'k-');
        leg{end+1}='sym';
        hold on
        
        [a,b]=hist(anBNS12(:,2),nBin);
        a=a/sum(a)/dn;
        hp(end+1)=plot(b,a,'r-');
        leg{end+1}='asym';
        
        legend(hp,leg,'location','best')
        
        axis([0 3 0 2.8])
        xlabel('N / N_{true}')
        ylabel('freq.')
    end
    % normalized background distribution
    if(1)
        figure(figVmax+10)
        subplot(nrows,ncols,kr)
        hold off
        hp=[];
        leg={};
        
        nBin=linspace(0,4,301);
        dn=mean(diff(nBin));
        [a,b]=hist(snBNS(:,1),nBin);
        a=a/sum(a)/dn;
        
        hp(end+1)=plot(b,a,'k-');
        leg{end+1}='sym';
        hold on
        
        [a,b]=hist(anBNS12(:,1),nBin);
        a=a/sum(a)/dn;
        hp(end+1)=plot(b,a,'r-');
        leg{end+1}='asym';
        
        legend(hp,leg,'location','best')
        
        axis([-0.01 3 0 2.8])
        xlabel('b / b_{true}')
        ylabel('freq.')
    end
    pause(0.1)
end
