function make_Coomans_simca(XfitResults,XvalResults,XtestResults, Method, rtype, numF)
% this function reads output from NSIMCA and calculate orig SIMCA class
% distance and make Coomans plots
% 
% make_Coomans_simca(XfitResults,XvalResults,XtestResults, Method, rtype, numF)

%
%% INPUT
%% XfitResults, XvalResults,XtestResults  structured array(factors,
%% class) from NSIMCA_0
%% 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))
%% 
%%  numF model dimensionality to be cosidered
%%

%% OUTPUT
%%  plot and save simca distances and F for F-test in files:
%%
%%  distsimca & Fsimca
%% (structured array dimensioned 1 x number of classes
%% 

%get number of factors and number of classes
[nf nc]=size(XvalResults);
% loop over classes
for indc=1:nc

    fmin =XfitResults(1,indc).finit;
    Fmax =XfitResults(1,indc).ftot;
    DimXcal=XfitResults(1,indc).DimXcal;
    ord=size(DimXcal,2);    
    NumFac=numF(indc);
    if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
        Nfac(1:ord)=NumFac;
    else
        PossibleCombs=XfitResults(NumFac,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
    dof_CV = prodM./(DimXcal(1)-Nfac(1)-1);
    %dof_fp=XftiResults.
    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
    
    [Sp2train Sp2trainCV Sp2train_fp So2 So2_fp So2CV Dp2 Dp2_fp Dp2est Dp2est_fp Dp2_CV Dp2estCV]=class_eval(numF(indc),indc,Method, rtype, DimXcal,XfitResults,XvalResults,XtestResults);
    distsimca(indc).Sp2train{1}=Sp2train;
    distsimca(indc).Sp2train{2}=Sp2trainCV;
    
    if ord > 2; distsimca(indc).Sp2train{3}=Sp2train_fp;end
    
    distsimca(indc).Sp2_o{1}=Dp2;
    distsimca(indc).Sp2_o{2}=Dp2_CV;
    if ord > 2; distsimca(indc).Sp2_o{3}=Dp2_fp;end
    
    distsimca(indc).Sp2_est{1}=Dp2est;
    distsimca(indc).Sp2_est{2}=Dp2estCV;    
    if ord >2; distsimca(indc).Sp2_est{3}=Dp2est_fp; end

    save distsimca  distsimca

    Fcalc_train=Sp2train./So2;
    Fcalc_trainCV=Sp2trainCV./So2CV;
    %
    if ord > 2; Fcalc_train_fp=Sp2train_fp./So2_fp; end
    %
    
    % critical F value
    if ord > 2;
        Ftable_fp=ftest(0.05,dof_CV_fp, dof_fp);
        if isnan(Ftable_fp)
            disp('WARNING: Fcrit = NaN, compute outside matlab')
        end
    end
    Ftable=ftest(0.05, dof_CV, prodM);
    if isnan(Ftable)
        disp('WARNING: Fcrit = NaN, compute outside matlab')
    end
    
    Flimtr = Ftable.*So2;
    FlimCV = Ftable.*So2CV;
    if ord >2; Flimfp = Ftable_fp.*So2_fp; end
    
    Fsimca(indc).Flimtr=Flimtr;
    Fsimca(indc).FlimCV=FlimCV;
    if ord >2; Fsimca(indc).Flimfp=Flimfp; end
    
    
    save Fsimca Fsimca


    %%
end


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

function [Sp2train Sp2trainCV Sp2train_fp So2 So2_fp So2CV Dp2 Dp2_fp  Dp2est Dp2est_fp Dp2_CV Dp2estCV]=class_eval(indf,indc,Method, rtype, DimXcal,XfitResults,XvalResults,XtestResults)
    
    %% 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
    H_others=[XtestResults(indf,indc).hT_test];     
    Q_others=[XtestResults(indf,indc).Qtest];   
    
    %% 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];
    end
        
    %% 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 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;
    scores_others=XtestResults(indf,indc).scoresnew;
    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(:,1:3);
    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;  % test
        Sp2train_fp=Q_fit./dof_fp;  % train a single obj
    
        %% Sp2trainCV_fp=Q_CV./dof_CV_fp; %train CV a single
        Sp2_fp=Q_others./dof_CV_fp;
        if isfield(XtestResults,'Qtest_test')
            Sp2est_fp=Q_test_est./dof_CV_fp;
        else
            Sp2est_fp=[];
        end
    else
        So2CV_fp=[];
        Sp2train_fp=[];
        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;
    if isfield(XtestResults,'Qtest_test')
        Sp2est=Q_test_est.*inv_dof_CV;
    else
        Sp2est=[];
    end
    % 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);end

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

    %%
    %% CV
    [Dp2_CV]=scoreDist(scores_others, Nfac, MaxLimtheta, MinLimtheta, phiaCV, Sp2);
    if isfield(XtestResults,'Qtest_test')
        [Dp2estCV]=scoreDist(Scoreest, Nfac, MaxLimtheta, MinLimtheta,  phiaCV, Sp2est);
    else
        Dp2estCV=[];
    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
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
   