% localization example 3, step 2: analysis

if(1) % load the data
    addpath(genpath('../uncertainSPT/'))
    clear
    
    ROIwidth=9;
    px2nm=80;
    
    % file
    resFolder='EX3results';
    
    sgStr={'truth->MLE','truth->MAP','FRS->MLE','FRS->MAP'}
    agStr={'truth->MLE','truth->MAP','FRS->MLE','FRS->MAP'}
    
    psfStr={'514nm NA1.49','','639nm NA1.4','680nm NA1.40'};
    
    Rfiles=dir(fullfile(resFolder,'*RW*.mat'));
    NR=numel(Rfiles);
    
    vars={'Nphotons','bg','tE','dt','psfString','psfIndex','emTrj_px',...
        'SMopt','sigMin',...
        'NA','lambda','frsKept','frsDetected'...
        'PSFsymGauss', 'SG_coord','SG_xyErr','SG_bNSdr',   'SG_cov','SG_CRLB',...
        'PSFasymGauss', 'AG_coord','AG_xyErr','AG_bNS12vdr',   'AG_cov','SG_CRLB',...
        };
    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;
            
            % recompute CRLB 
            a=1;
            s2a=Rk.SG_bNSdr{m}(:,3).^2+a^2/12;
            N  =Rk.SG_bNSdr{m}(:,2);
            b  =Rk.SG_bNSdr{m}(:,1);
            tau =2*pi.*b.*s2a./N/a^2;
            F0=1+4*tau+sqrt(2*tau./(1+4*tau)); % approximation from Rieger&Stallinga        
            Rk.SG_CRLB{m} = 2*s2a./N.*F0;
        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
    ROIwidth=9;
    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
% plot
if(1)
    figVmax=10;
    ncols=4
    nrows=1
    nbst=3
    plotCSA=[1 1 1];
    
    % estimated RMS bins
    Lrms=150;
    NBmin=300; % minimum number of dots/bin
    dRMS=7.5;
    RMSedge=dRMS/2:dRMS:Lrms;
    Lxy=100;

    sLaplCol='m';
    markSize=4;
    
    r_fs=[1 3 2 4];
    r_fa=[1 3 2 4];
    mFRS=[3 4];
    
    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)(sgStr{r_fs(kr)});
    posFun=@(ax,kr)(1)%set(ax,'position',[0.35 0.24 0.63 0.75]));
    for kr=1:numel(r_fs)
        %% collect data
        fs=r_fs(kr);        
        fa=r_fa(kr);        
        
        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=[];
        
        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','frsDetected','frsKept',...
                                'SG_bNSdr','SG_cov','SG_CRLB','SG_xyErr','Sn',...
                                'AG_bNS12vdr','AG_cov','AG_xyErr','aSn'});                            
                            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);
                            
                            % detected by the FRS algorithm
                            if(ismember(fs,mFRS))
                                cfrs = A.frsKept==1;
                            else
                                cfrs=true(size(A.frsKept));
                            end
                            
                            % symmetric Gaussian
                            ccSG= c0SG & cRs & cUSG & cz & cSsg & cfrs;
                            % asymmetric Gaussian
                            ccAG= c0AG & cRa & cUAG & cz & cSag & cfrs;
                            
                            % 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)];
                            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))];
                            end
                            if(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
        % conditional normalizedf errors
        if(1)
            figure(figVmax+15)
            subplot(nrows,ncols,kr)
            hold off
            pause(0.1)
            markSize=4;
            % 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');
                hp(end+1)=plot(sqrt(XC)-dRMS/5,YM,'ks','markerface','w','markersize',markSize)
                leg{end+1}='CRLB';
            end
            if(plotCSA(2)) % RMSD | sym. Laplace
                XC=mean(xcSL,1,'omitnan');
                YM=median(yncSL,'omitnan');
                hp(end+1)=plot(sqrt(XC)+dRMS/5,YM,'co','markerface','w','markersize',markSize,'color',sLaplCol);
                leg{end+1}='sym.L.';
            end
            if(plotCSA(3)) % RMSD | asym. Laplace
                XC=mean(xcAL,1,'omitnan');
                YM=median(yncAL,'omitnan');
                hp(end+1)=plot(sqrt(XC)-dRMS/5,YM,'b^','markerface','w','markersize',markSize);
                leg{end+1}='asym.L.';
            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'))
            legend(hp,leg,'location','best')
            
            xlabel('est.RMSE [nm]')
            ylabel('true/est. cRMSE [nm]')
            grid off
            set(gca,'xscale','lin','yscale','lin')
            axis([dRMS/2 Lxy+9 0.55 1.85])
            posFun(gca,kr);
            %%%pp=get(gca,'position');pp([2 3 4])=[0.2 0.135 0.67];set(gca,'pos',pp)
            title(titleFun(kr))
            %set(gca,'xtick',get(gca,'ytick'))
        end
        % absolute errors
        if(1)
            figure(figVmax+5)
            subplot(nrows,ncols,kr)
            hold off
            markSize=4;
            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([5 Lxy 5 Lxy])
            legend(hp,leg,'location','best')
            title(titleFun(kr))
            posFun(gca,kr);
            %%%pp=get(gca,'position');pp([2 3 4])=[0.2 0.135 0.67];set(gca,'pos',pp)

            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         
        % 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);
            pp=get(gca,'position');pp([2 3 4])=[0.25 0.135 0.62];set(gca,'pos',pp)
            set(gca,'xtick',0:20:100)
        end
        pause(0.1)
    end
end
