function  Nclass(XfitResults,XvalResults,XtestResults, Method, rtype)
% this function reads output from NSIMCA and:
% 
% - make class assignments by using different classification rules
% - calculate SENS and SPEC for each class and each factor
% - also for number of factor corresponding to max R2CV according to
%   ncrossdecomp
%
%% INPUT
%% XfitResults, XvalResults,XtestResults  structured array(factors,
%% class) from NSIMCA
%% Method (par or tuc)
%% rtype says type of range for scores to be used in SIMCA original
%% rtype=1 (paper original Wold tmax+1/2*t_student*std(t))
%% rtype=2 normal range, i.e. tmax
%% rtype=3 ( tmax + std(t))
%%
%**** Nclass(XfitResults,XvalResults,XtestResults, Method,rtype)

%% OUTPUT
%%
%% classResults: 
%% contains SENSf (sensitivity), SPECf (specificity) and classtestf (1 accepted; 0 rejected) 
%% 
%% Each column corresponds to a different classification rule as follow:
%%
%%  1. H_fit/Q_fit (limF_hTfit Qlim_fit)  % FORINA changed May 09
%%  2. H_fit/Q_fit (Hrlim_CV Qrlim_CV)
%%  3. H_fit/Q_fit (hT_lim_fitAP Qlim_fitAP) % POMERENTSEV fit added May 09
%%  4. H_fit/Q_fit (limAP_hT_CV limAP_QCV)    % POMERENTSEV CV changed May 09

%%  5. simca Wold fit
%%  6. simca Wold CV
%%  7. simca Wold fp  (not defined for 2D data)

%% Each row to a category; third dimension is the corrspondinf model dimensionality 
%%

%% Furthermore SENS and SPEC with the different criteria are also calculated for CV predicted samples 
%% and saved as SENS_CV and SPEC_CV  with columns corresponding to: 

%%  1. H_CV/Q_CV (limF_hTfit Qlim_fit)
%%  2. H_CV/Q_CV (Hrlim_CV Qrlim_CV)  
%%  3. H_CV/Q_CV (hT_lim_fitAP Qlim_fitAP) 
%%  4. H_CV/Q_CV (limAP_hT_CV limAP_QCV)    

%get number of factors and number of classes
[nf nc]=size(XvalResults);

% loop over classes
for i=1:nc
    fmin(i)=XfitResults(1,i).finit;
    Fmax(i)=XfitResults(1,i).ftot;
    DimXcal=XfitResults(1,i).DimXcal;
    
    
    
    for ff=fmin(i):Fmax(i)
        [SENS, SPEC, classtest, SENS_CV, SPEC_CV]=class_eval(nc,ff,i,Method, rtype, DimXcal,XfitResults,XvalResults,XtestResults);
        SENSf(i,:,ff)=SENS;
        SENSf_CV(i,:,ff)=SENS_CV;
        if ~isempty(SPEC_CV); SPECf_CV(i,:,ff)=SPEC_CV;else SPECf_CV(i,:,ff)=0;end
        if ~isempty(SPEC); SPECf(i,:,ff)=SPEC;else SPECf(i,:,ff)=0;end
        classtestf(ff,i).cassTest=classtest;        
        
    end
    %%
end
save classResults SENSf SPECf classtestf
save CVResults SENSf_CV SPECf_CV



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [SENS, SPEC, CLASStest,  SENS_CV, SPEC_CV]=class_eval(nc,indf,indc,Method, rtype, DimXcal,XfitResults,XvalResults,XtestResults,Ftab,Ftablefp)
    %% statistics corresponding to minum PRESS LV's %%%%%
    
    %% fit %%    
    H_fit=[XfitResults(indf,indc).hT_fit];    
    Q_fit=[XfitResults(indf,indc).Qfit];  
    Qlimit_fit=[XfitResults(indf,indc).Qlimit_fit];  

    Qlimit_fitAP=[XfitResults(indf,indc).Qlimit_fitAP];  %_added May 09 Pomerntsev 
    Hlim_fitAP=[XfitResults(indf,indc).Hlim_fitAP];  
    limF_hTfit=[XfitResults(indf,indc).limF_hTfit];      % added May 09 Forina

    %% CrossValidated
    H_CV=[XvalResults(indf,indc).hT_CV];    
    Q_CV=[XvalResults(indf,indc).Q_CV];      
    Qlimit_CV=[XvalResults(indf,indc).Qlimit_CV];
    
    limAP_QCV=[XvalResults(indf,indc).limAP_QCV];     %_added May 09 Pomerntsev 
    limAP_hT_CV=[XvalResults(indf,indc).limAP_hT_CV];  
 
    % 95 Percentile
    Hrlim_CV = [XvalResults(indf,indc).hTrlim_CV];
    Qrlim_CV = [XvalResults(indf,indc).Qrlim_CV]; % added Feb 07
    
    
    %% other  classes
    if nc > 1
        H_others=[XtestResults(indf,indc).hT_test];
        Q_others=[XtestResults(indf,indc).Qtest];
    else
        H_others=[];
        Q_others=[];
    end
    %% external test set 
    if isfield(XtestResults,'hT_test_test')
        H_test_est=[XtestResults(indf,indc).hT_test_test];
        Q_test_est=[XtestResults(indf,indc).Qtest_test];
    else
        H_test_est=[];
        Q_test_est=[];
    end
    %%% SENS SPEC according to the differnt classification rules
    %% PLS toolbox philosophy (equal weight to T2 and Q for the moment) 
    
    %% 1- leverage/Qfit  Forina beta distr. is used
    [SENS_1, SPEC_1, asstest_1]=calc_sensspec(H_fit, Q_fit, H_others, Q_others, H_test_est, Q_test_est, limF_hTfit, Qlimit_fit);

    [SENSCV_1, SPECCV_1]=calc_sensspec(H_CV, Q_CV, H_others, Q_others, H_test_est, Q_test_est, limF_hTfit, Qlimit_fit);

    %% 2- leverageCV /QCV  with rlim (95% percentile)
    [SENS_2, SPEC_2, asstest_2]=calc_sensspec(H_fit, Q_fit, H_others, Q_others, H_test_est, Q_test_est, Hrlim_CV, Qrlim_CV);     

    [SENSCV_2, SPECCV_2]=calc_sensspec(H_CV, Q_CV, H_others, Q_others, H_test_est, Q_test_est, Hrlim_CV, Qrlim_CV);     
    

    %% 3- leverage/Qfit Pomentsev chi2 distr. is used sep. for H and Q
    [SENS_3, SPEC_3, asstest_3]=calc_sensspec(H_fit, Q_fit, H_others, Q_others, H_test_est, Q_test_est, Hlim_fitAP, Qlimit_fitAP);
    
    [SENSCV_3, SPECCV_3]=calc_sensspec(H_CV, Q_CV, H_others, Q_others, H_test_est, Q_test_est, Hlim_fitAP, Qlimit_fitAP);

    %% 4- leverageCV Pomerentsev
    [SENS_4, SPEC_4, asstest_4]=calc_sensspec(H_fit, Q_fit, H_others, Q_others, H_test_est, Q_test_est, limAP_hT_CV, limAP_QCV); 

    [SENSCV_4, SPECCV_4]=calc_sensspec(H_CV, Q_CV, H_others, Q_others, H_test_est, Q_test_est, limAP_hT_CV, limAP_QCV); 
    
 
    %% SIMCA orig 

    %% computation of degree of freedom for OD according to G.R. Flaten et al. J.Chemometrics 2004 (18) 173-182
    %%  N.B. original SIMCA do not use dof nor distribution for SD since it is combined  with OD before 
    %%  invoking any distribution, the F-test is used not because of Hotelling T2 but because it is a an 
    %%  F-test on variances!!
    
    %%%% in fit and CV
    
    NumFac=indf; %factors for model of class i
        
    scores_fit=XfitResults(indf,indc).factors_fit{1};
    Nobj=size(scores_fit,1);
    scores_CV=XvalResults(indf,indc).scores_CV;
    if nc >1; 
        scores_others=XtestResults(indf,indc).scoresnew; 
    else scores_others=[]; 
    end
    if isfield(XtestResults,'scoresnew_test')
        Scoreest=XtestResults(indf,indc).scoresnew_test;
    else
        Scoreest=[];
    end
    Stheta=missstd(scores_fit); % score standard deviation as in Forina etc.. (i.e. divided by n-1)
    SthetaCV=sqrt(var(scores_CV).*((Nobj-1)/Nobj));      
    thetamax=max(scores_fit);
    thetamin=min(scores_fit);
    thetamaxCV=max(scores_CV);
    thetaminCV=min(scores_CV);
    tt=tinv(0.025,size(scores_fit,1));
    
     %  Equation 23 of Wold and Sjostrom:
     if rtype==1
         MaxLimtheta=thetamax+(1./2)*tt*Stheta;
         MaxLimthetaCV=thetamaxCV+(1./2)*tt*SthetaCV;
         %  Equation 24 of Wold and Sjostrom:
         MinLimtheta=thetamin-(1./2)*tt*Stheta;
         MinLimthetaCV=thetaminCV-(1./2)*tt*SthetaCV;
     elseif rtype==2
         MaxLimtheta=thetamax;
         MaxLimthetaCV=thetamaxCV;
         MinLimtheta=thetamin;
         MinLimthetaCV=thetaminCV;
     elseif rtype==3
         MaxLimtheta=thetamax+ Stheta;
         MaxLimthetaCV=thetamaxCV+ SthetaCV;
         MinLimtheta=thetamin- Stheta;
         MinLimthetaCV=thetaminCV-SthetaCV;
     end
    
    % class orthogonal distance
    ord=size(DimXcal,2);
    if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
        Nfac(1:ord)=NumFac;
    else
        PossibleCombs=XfitResults(indf,indc).PossibleCombs;
        Nfac=PossibleCombs;
    end
    prodM=1;
    for i=1:ord
        if i==1
            Nf=Nfac(i)+1;
        else
            Nf=Nfac(i);
        end
        prodM=prodM.*(DimXcal(i)-Nf);
    end
    inv_dof=(DimXcal(1)-1)./prodM;                  % according to Flaten training objs
    inv_dof_CV=1./(prodM./(DimXcal(1)-Nfac(1)-1));  % according to Flaten cv, others, test objs
    dof_CV=prodM./(DimXcal(1)-Nfac(1)-1);
    
    %% dof estimated from the number of free parameters Kroonenberg to be
    %% implemented May 09
    if ord > 2
        if lower(Method)=='tuc'
            nfp=calcnfp(PossibleCombs,DimXcal,Method); % number of free parametrs
        else
            nfp=calcnfp(Nfac,DimXcal,Method);
        end
        
        dof_fp=prod(DimXcal)-nfp ; % training objs
        dof_CV_fp=dof_fp + DimXcal(1)*Nfac(1) ; % cv, ohers, test objs but I am not sure of this
    end
    %%
    So2=sum(Q_fit).*inv_dof; % train
    % with nfp
    if ord > 2; So2_fp=sum(Q_fit)./dof_fp; else So2_fp=[]; end % train
    %
    So2CV=sum(Q_CV).*inv_dof_CV;  %  test
    Sp2train=Q_fit.*inv_dof;   % this is the part of So2 due to a single training object CONTROLLARE
    %
    if ord > 2; So2CV_fp=sum(Q_CV)./dof_CV_fp; else So2CV_fp=[]; end % test
    if ord > 2; Sp2train_fp=Q_fit./dof_fp; else Sp2train_fp=[];end % train a single obj
    %% Sp2trainCV_fp=Q_CV./dof_CV_fp; %train CV a single 
    if ord >2; Sp2_fp=Q_others./dof_CV_fp; else Sp2_fp=[];end
    if ord >2; Sp2est_fp=Q_test_est./dof_CV_fp; else Sp2est_fp=[];end

    %
    Sp2trainCV=Q_CV.*inv_dof_CV;  % this is the part of So2CV due to asingle cross-validated object CONTROLLARE
    Sp2=Q_others.*inv_dof_CV;
    Sp2est=Q_test_est.*inv_dof_CV;
    
    % Equation 26 from Wold and Sjostrom:
    %  Note missprint in Wold and Sjostrom, use So instead of Sp
    phia=sqrt(So2)*(ones(1,Nfac(1))./Stheta);
    %
    if ord > 2; phia_fp=sqrt(So2_fp)*(ones(1,Nfac(1))./Stheta);else phia_fp=[];end

    %
    phiaCV=sqrt(So2CV)*(ones(1,Nfac(1))./SthetaCV);
   
    % class score distance
    if nc > 1; 
        [Dp2]=scoreDist(scores_others, Nfac(1), MaxLimtheta, MinLimtheta, phia, Sp2);
    else Dp2=[]; 
    end
    if ~isempty(Scoreest)
        [Dp2est]=scoreDist(Scoreest, Nfac(1), MaxLimtheta, MinLimtheta, phia, Sp2est);
    end
    
    %%
    if ord > 2
        if nc >1;
            [Dp2_fp]=scoreDist(scores_others, Nfac(1), MaxLimtheta, MinLimtheta, phia_fp, Sp2_fp);
        else Dp2_fp=[];
        end
        [Dp2est_fp]=scoreDist(Scoreest, Nfac(1), MaxLimtheta, MinLimtheta, phia_fp, Sp2est_fp);
    else
        Dp2_fp=[]; 
        
        Dp2est_fp=[];
        
    end
    %%
    %% CV
    if nc >1
        [Dp2_CV]=scoreDist(scores_others, Nfac(1), MaxLimtheta, MinLimtheta, phiaCV, Sp2);
    else
        Dp2_CV=[];
    end
    if ~isempty(Scoreest)
        [Dp2estCV]=scoreDist(Scoreest, Nfac(1), MaxLimtheta, MinLimtheta, phiaCV, Sp2est);
    end

    %%%
      
    
    Fcalc_train=Sp2train/So2;
    Fcalc_trainCV=Sp2trainCV/So2CV;
    %
    if ord > 2; Fcalc_train_fp=Sp2train_fp/So2_fp;else Fcalc_train_fp=[]; end
    %
    
    % critical F value
    if ord > 2; Ftable_fp=ftest(0.05,dof_CV_fp, dof_fp);else Ftable_fp=[];end
    Ftable=ftest(0.05, dof_CV, prodM);
    

    if ~isempty(Scoreest)
        [SPEC_orig, asstest_orig, SENS_orig]=calc_simca(Ftable, Dp2, So2, Nobj, size(H_others,1),Fcalc_train, Dp2est, size(H_test_est,1));
        [SPEC_origCV, asstest_origCV, SENS_origCV]=calc_simca(Ftable, Dp2_CV, So2CV, Nobj, size(H_others,1), Fcalc_trainCV, Dp2estCV, size(H_test_est,1));
    else
        [SPEC_orig, asstest_orig, SENS_orig]=calc_simca(Ftable, Dp2, So2, Nobj, size(H_others,1),Fcalc_train);
        [SPEC_origCV, asstest_origCV, SENS_origCV]=calc_simca(Ftable, Dp2_CV, So2CV, Nobj, size(H_others,1), Fcalc_trainCV);
    end
    %%
    if ord > 2
        if ~isempty(Scoreest)
            [SPEC_orig_fp, asstest_orig_fp, SENS_orig_fp]=calc_simca(Ftable_fp, Dp2_fp, So2_fp, Nobj, size(H_others,1), Fcalc_train_fp, Dp2est_fp,size(H_test_est,1));
        else
            [SPEC_orig_fp, asstest_orig_fp, SENS_orig_fp]=calc_simca(Ftable_fp, Dp2_fp, So2_fp, Nobj, size(H_others,1), Fcalc_train_fp);
        end
    else
      SPEC_orig_fp=[]; asstest_orig_fp=[]; SENS_orig_fp=[];
    end
    %%
    %
    SENS = [SENS_1 SENS_2 SENS_3 SENS_4 SENS_orig SENS_origCV SENS_orig_fp];

    SENS_CV = [SENSCV_1 SENSCV_2 SENSCV_3  SENSCV_4];    

    SPEC =  [SPEC_1 SPEC_2 SPEC_3 SPEC_4 SPEC_orig SPEC_origCV SPEC_orig_fp ];
    
    SPEC_CV =  [SPECCV_1 SPECCV_2 SPECCV_3 SPECCV_4];

    CLASStest = [asstest_1 asstest_2 asstest_3 asstest_4 asstest_orig asstest_origCV asstest_orig_fp ];

    %%

    
 %%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%
 
 function [sensibility, specificity, assign_test]=calc_sensspec(scoreD, orthD, scoreDother, orthDother, scoreDtest, orthDtest, scoreDlim, orthDlim)
 %%
 Naccepted=length(find(sqrt((scoreD./scoreDlim).^2 + (orthD./orthDlim).^2)<=sqrt(2)));
 if ~isempty(scoreDother); Nrejected=length(find(sqrt((scoreDother./scoreDlim).^2 + (orthDother./orthDlim).^2)>sqrt(2))); end
 sensibility=100*(Naccepted./size(scoreD,1));
 if ~isempty(scoreDother); specificity=100*(Nrejected./size(scoreDother,1));else specificity=[]; end
%% test prediction
if ~isempty(scoreDtest)
    assign_test=zeros(size(scoreDtest,1),1);
    iass=find(sqrt((scoreDtest./scoreDlim).^2 + (orthDtest./orthDlim).^2)<=sqrt(2));
    assign_test(iass)=1;
else
    assign_test=[];
end
%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function [Dp2]=scoreDist(scores, NumFac, maxl, minl, phi, Sp2)
     for ll=1:size(scores,1)
         for zz=1:NumFac
             if scores (ll,zz) > maxl(zz)
                 ScoreDiff(ll,zz)=scores(ll,zz)-maxl(zz);
             elseif scores(ll,zz) < minl(zz)
                 ScoreDiff(ll,zz)=scores(ll,zz)-minl(zz);
             else
                 ScoreDiff(ll,zz)=0;
             end %if
         end %for zz
         ScoreDiff2(ll,1)=sum((phi.^2).*(ScoreDiff(ll,:).^2));
         Dp2(ll,1)=Sp2(ll,1)+ScoreDiff2(ll,1);
     end %for ll
    
    %%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%
    function [SPEC_sim, asstest_sim, SENS_sim]=calc_simca(Ftable, dist2, So2, Nobj, Nother, Fcal_train, dist2est, Ntest)
    if Fcal_train ~=0
        Naccepted=length(find(Fcal_train <= Ftable));
        SENS_sim=100*(Naccepted./Nobj);
    end
    if Nother >0
        Fcal=dist2/So2;
        Nrejected=length(find(Fcal >= Ftable));
        SPEC_sim=100*(Nrejected./Nother);
    else
        SPEC_sim=[];
    end
   if nargin > 6
       Fcalest=dist2est/So2;
       asstest_sim=zeros(Ntest,1);
       iass=find(Fcalest <=Ftable);
       asstest_sim(iass)=1;
   else
       asstest_sim=[];
   end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   