function [XfitResults,XvalResults,XtestResults] = NSIMCA(Method,FacMax,FacMin,Xnam,classid,alpha,Pret,ordPret,indblk)
%% This is main of NSIMCA calculation
%%
%% NSIMCA performs disjoint N-way components models for each category
%% (N.B. Samples are in MODE1)
%%
%% Classification rules ar ebased on distances (both original-SIMCA and alternative-SIMCA frameworks)
%% Scores Distance is based on Leverage (H); Orthogonal Distance on Residuals (Q)
%% Calculation of H and Q statistics according to several criteria for assessing statistical limits:
%% 
%% limF_Hfit:  Forina et al., Chemom. Intel. Lab. Syst. 96 (2009) 239?245.
%% limAP_h: Pomerantsev et al.,  J. Chemometrics 22 (2008) 601-609.
%% limAP_H_CV Pomerantsev in CV
%% HTrlim_CV: uses the 95% percentile
%%
%% Qfit_lim: Jackson -Mudholkar approximation
%% limAP_Q: Pomerantsev et al.
%% Qrlim_CV: uses the 95% percentile
%% limAP_QCV: uses Pomerantsev with CV residuals
%%
%% H and Q are also predicted for the objects belonging to the non modeled class
%%
%%  input
%%     Method  = decomposition Method, PARAFAC ('par') or TUCKER ('tuc')..
%%               for 2-way array use PARAFAC, that gives same results as PCA
%%     FacMax  = maximum number of factors to be tested, eg. for PARAFAC ([10; ...; ...]) 
%%               for Tucker ([5 7 12; 3 5 8; ...]) 
%%               indicized n. of classes x order of the array
%%     FacMin  = minimum number of factors to be tested same format as FacMax
%%     Xnam    = name of data array
%%     classid = name of class assignation vector, the vector has to have as many elements as objects in
%%              MODE1, eg. [ 1 1 2 1 1 2 2 2 ..]
%%     Pret = ['  '; '   '; '   '];
%%
%%     Pret = 'non'  no pretreatment is applied
%%     Pret = 'mnc'  mean centering across given MODE  (Default if Pret is empty is 'mnc' on first MODE) 
%%     Pret = 'psc'  pareto scaling, i.e. weight = 1/sqrt(st.dev)
%%     Pret = 'aut'  standard deviation scaling: i.e. weight = 1/st.dev 
%%     Pret = 'blk'  blockscaling
%%     indblk = is a strctured array indicized with respect to number of modes.
%%            It contains indblk as field (indblk.indblk{..}), that is a block assignation vector,
%%            as many elements as variables in given MODE; needed only if Pret='blk', else []
%%     
%%     ordPret = order of pretreatment es: [1 2 3] first pretreat Mode 1 then mode 2 then Mode 3
%%
%%     alpha, e.g. can be equal to 0.05 or 0.01
%%     alpha= significance level for T2 and Q distribution
%%
%%  output
%%     XfitResults
%%     XvalResults
%%     XtestResults
%%     structured arrays with fields: H, Q, Hlim, Qlim, rmsec, r2 for
%%     calibration (i.e. a given class), crossvalidation and prediction (
%%     i.e. the other classes) respectively.
%% 
%%    [XfitResults,XvalResults,XtestResults] = NSIMCA(Method,FacMax,FacMin,Xnam,classid,alpha,Pret,ordPret,indblk);

%%  $ Version 1 $ Feb 2013 $ Marina Cocchi $ Not compiled $
%%   please cite:  
%%   C. Durante,R. Bro,M. Cocchi, A classification tool for N-way array based on SIMCA methodology
%%   Chemom. Intell. Lab. Syst. 2011 (106) 73-85.
%%
%%
% set parameters for decomposition algorithms
it = 0;
maxit = 1000;   % maximum number of iterations
const= [2 2 2]; %
% read input
Xtot = getfield(load(Xnam), Xnam);
classid = getfield(load(classid), classid);
if exist('indblk')
   indblk = getfield(load(indblk), indblk);
else
    indblk=[];
end
% set default pretreatment as centering in MODE1
if ~exist('Pret')
    Pret=['mnc';'non';'non'];
    ordPret=[1 2 3];
end
%
%
DimX = size(Xtot);
ord = length(DimX);
%
Xreshape=reshape(Xtot,DimX(1),prod(DimX(2:end)));
[I,J] = size(Xreshape);
clear Xreshape % to save space
%

% Determine the number of classes and X_calibration and X_test
if length(classid)~=DimX(1)
   error('Input (classid) must be a vector with the same number  of elements as number of rows of (data).');
end

[class,ind]   = sort(classid);   % Sort the data by class
indtest=ind(class==0);
Xtest_test=Xtot(indtest,:,:); 
ind(class==0) = [];     %drop class "zero" if there
x = Xtot(ind,:,:); % 3D array sorted according to classes 
clear Xtot % to save space
classrid=classid(ind);
classes=unique(classrid);
noclass = length(classes);     %number of classes

%% choose CV criteria and split for each class
for ii=1:noclass

    I=length(find(class==classes(ii)));
    [SegmentsID xvalm]= MakeSegments(I);
    CVparam(ii).Segments=SegmentsID;
    CVparam(ii).Xvalmeth=xvalm;
    
    
    Segments_ncd = 11;
    out = ones(I,J);
    out(1:Segments_ncd:end)=NaN;
    out2 = ones(size(x));
    out(find(isnan(out2))) = NaN;
    while any(sum(isnan(out))==I) | any(sum(isnan(out'))==J)
        Segments_ncd=Segments_ncd+2;
        out = ones(I,J);
        out(1:Segments_ncd:end)=NaN;
        out2 = ones(size(x));
        out(find(isnan(out2))) = NaN;
    end

end
%%
% prepare a vector which is the Fixmode parameter to be used in PARAFAC
% prediction 
if ord <=2; Fixm=[0 1]; else Fixm=[0 1 1];end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:noclass   % loop on classes
    ctrain=find(classrid==classes(ii));
    ctest=find(classrid~=classes(ii));
    Xcal  = x(ctrain,:,:);
    Xtest = x(ctest,:,:);    
    DimXtest_test=size(Xtest_test);
    DimXcal=size(Xcal);
    DimXtest=size(Xtest);
    
    %% do pretreatment
    if ord<=2
        if Pret(1,1:3) =='mnc' & Pret(2,1:3)~='blk'
            [Xcal_p, Mx_cal, Sx_cal]=scal2D(Xcal,Pret,[],[]);
            if ~isempty(Xtest)
                Xtest_p=scal2D(Xtest,Pret, Mx_cal, Sx_cal);
            end
            if ~isempty(Xtest_test)
                Xtest_test_p=scal2D(Xtest_test,Pret,Mx_cal,Sx_cal);
            end
        elseif Pret(1,1:3)=='mnc' & Pret(2,1:3)=='blk'
            Mx_cal=missmean(Xcal)
            Xm=Xcal - ones(size(Xcal,1),1)*Mx_cal;
            ind_blk=indblk(2).indblk;
            nblock=length(unique(ind_blk));
            SStot=missss(Xm);
            for iblk=1:nblock
                indb=find(ind_blk==iblk);
                SSblock(iblk)=missss(Xm(:,indb));
                Sx_cal(1,indb)=1./(sqrt(SStot./(SSblock(iblk).*nblock)));
            end
            Xcal_p=Xm./(ones(size(Xcal_p,1),1)*Sx_cal);
            clear Xm
            if ~isempty(Xtest)
                Xtest_p=Xtest - ones(size(Xtest,1),1)*Mx_cal;
                Xtest_p=Xtest_p./(ones(size(Xtest_p,1),1)*Sx_cal);
            end
            if ~isempty(Xtest_test)
                Xtest_test_p=Xtest_test - ones(size(Xtest_test,1),1)*Mx_cal;
                Xtest_test_p=Xtest_test_p./(ones(size(Xtest_test_p,1),1)*Sx_cal);
            end            
        elseif Pret(1,1:3)=='non' & Pret(2,1:3)=='blk'
            ind_blk=indblk(2).indblk;
            nblock=length(unique(ind_blk));
            SStot=missss(X_cal);
            for iblk=1:nblock
                indb=find(ind_blk==iblk);
                SSblock(iblk)=missss(Xcal(:,indb));
                Sx_cal(1,indb)=1./(sqrt(SStot./(SSblock(iblk).*nblock)));
            end
            Xcal_p=Xcal./(ones(size(Xcal_p,1),1)*Sx_cal);            
            if ~isempty(Xtest)
                Xtest_p=Xtest./(ones(size(Xtest,1),1)*Sx_cal);
            end
            if ~isempty(Xtest_test)
                Xtest_test_p=Xtest_test./(ones(size(Xtest_test,1),1)*Sx_cal);
            end
            Mx_cal=zeros(1,DimXcal(2));
        elseif Pret(1,1:3)=='non' & Pret(2,1:3)=='non'
            Xcal_p=Xcal;
            if ~isempty(Xtest)
                Xtest_p=Xtest;
            end
            if ~isempty(Xtest_test)
                Xtest_test_p=Xtest_test;
            end
            Mx_cal=zeros(1,prod(DimXcal(2:end)));
            Sx_cal=ones(1,DimXcal(2));
        end
    else
        [Xcal_p, Mx_cal, Sx_cal]=pretX_miss_0(Xcal, DimXcal, Pret, ordPret, indblk);
        if ~isempty(Xtest)
            Xtest_p=pretX_miss_0(Xtest, DimXtest, Pret, ordPret, indblk, Mx_cal, Sx_cal, 1);
        end
        if ~isempty(Xtest_test)
            Xtest_test_p=pretX_miss_0(Xtest_test,DimXtest_test, Pret, ordPret, indblk, Mx_cal, Sx_cal,1);
        end
    end
   %% end of pretreatment
    
   
    ssX=sum(sum(sum((Xcal).^2))); % Total SS in original unit to be used for %V calculation
    if ~isempty(Xtest)
        ssXtest=sum(sum(sum((Xtest).^2)));
    end
    %save space
    %clear Xtestm Xm 
    %
    if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
        ftot =  FacMax(ii);
        finit = FacMin(ii);
        possibleCombs=[];
    else  %% Method TUCKER
        % Find all    %possibleCombs
        PossibleNumber = [min(FacMin(ii,:)):max(FacMax(ii,:))]'*ones(1,ord);
        possibleCombs = unique(nchoosek(PossibleNumber(:),ord),'rows');
        %% added Ott 29 2012 to take only wanted combinations of factors
        for pp=1:size(possibleCombs,1)
            for kk=1:ord
                rr=(intersect([FacMin(ii,kk):FacMax(ii,kk)],possibleCombs(pp,kk)));
                if isempty(rr)
                    aa(pp,kk)=0;
                else
                    aa(pp,kk)=1;
                end
            end
        end
        aaok=find(prod(aa,2)>0);
        possibleCombs=possibleCombs(aaok,:);
        %%%% 
        %remove useless
        f2 = [];
        for f1 = 1:size(possibleCombs,1)
            if (prod(possibleCombs(f1,:))/max(possibleCombs(f1,:)))<max(possibleCombs(f1,:)) % Check that the largest mode is larger than the product of the other
                f2 = [f2;f1];
            elseif any(possibleCombs(f1,:)>FacMax(ii,:))
                f2=[f2;f1];
            end
        end
        possibleCombs(f2,:)=[];
        [f1,f2]=sort(sum(possibleCombs'));
        possibleCombs = [possibleCombs(f2,:) f1'];
        ftot=size( possibleCombs ,1);
        finit=1;
        if ftot >= 100
            sprintf('The number of TUCKER3 models to test are:%i',ftot)
            warning('Consider to abort this run, and pre-screen combination to be tested, e.g. run nvalidate prior to run NSIMCA')
        end
        
    end
    XfitResults(1,ii).DimXcal=DimXcal;
    XfitResults(1,ii).finit=finit;
    XfitResults(1,ii).ftot=ftot;
    
    for f = finit:ftot  %loop on factors the results are saved for each explored dimensionality
        %%%
        %%% MODEL for CALIBRATION SET (i.e. class ii)
        %%%
        if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
            Factors_f=parafac(Xcal_p,f,[1e-5 10 0 2 NaN maxit],const);
            G=[];
            nlv1=f;
        else %% Method TUCKER
            [Factors_f,G] = tucker(Xcal_p,possibleCombs(f,1:end-1),[1e-5 0 0 2 NaN maxit]);
            nlv1=possibleCombs(f,1);
        end
        if ord<=2
            [scores_fit load1_fit]=fac2let(Factors_f);
            load2_fit=[];
        else
            [scores_fit load1_fit load2_fit]=fac2let(Factors_f);
        end
        Xfit_cal=nmodel(Factors_f, G); % fitted data
        
        %%%%%%%%%  rescale recalculated data to original unit %%%%%%%%%%%%%%%
    
        if ord<=2
            if Pret(1,1:3) =='mnc' & Pret(2,1:3)~='blk'
                Xfit_o=rescal2D(Xfit_cal,Mx_cal,Sx_cal);
            elseif Pret(1,1:3)=='non' & Pret(2,1:3)=='non'
                Xfit_o=Xfit_cal;
            elseif Pret(1,1:3)=='mnc' & Pret(2,1:3)=='blk'
                Xfit_o=Xfit_cal.*(ones(size(Xfit_cal,1),1)*Sx_cal);
                Xfit_o=Xfit_o + ones(size(Xfit_cal,1),1)*Mx_cal;
            elseif Pret(1,1:3)=='non' & Pret(2,1:3)=='blk'
                Xfit_o=Xfit_cal.*(ones(size(Xfit_cal,1),1)*Sx_cal);
            end
        else
            Xfit_o=pretX_miss_0(Xfit_cal, DimXcal, Pret, ordPret, indblk, Mx_cal, Sx_cal, -1);
        end
        
        clear Xfit_cal %% Save space
        %% end of rescale
        
        
        resid_fit=Xcal-Xfit_o ; % 3-way residuals array
        resid2_fit=(Xcal-Xfit_o).^2;  % squared model residuals
        if ord<=2
            fit1=resid2_fit';
            fit2=sum(fit1);
        else
            fit1=permute(resid2_fit,[3 2 1]);
            fit2=sum(sum(fit1));
        end
        
        Qfit=shiftdim(fit2); % Q or model residuals
        
        clear fit1 fit2 %% save space
        
        %%%% calculate Q-statistic using residuallimit PLS_toolbox routine
        %%%% Jackson & Mudholkar %%%%%%%%%%%%%%
        limit=residuallimit(resid_fit,1-alpha);
        
        %%%% calculate Dof for Q or OD according to A. Pomerantsev
        %%%% J.Chemom. 2008 by using IQR approximation modified by me to
        %%%% get reasonable numbers i.e. :
        %%%% df=(exp(4.36ln(1.24/(IQR/ho)))^0.72
        dfAP_Q=(exp(4.36*log(1.24/(iqr(Qfit)./mean(Qfit)))))^0.72;
        limAP_Q=(mean(Qfit)*chi2inv(1-alpha,dfAP_Q))./dfAP_Q;
        % save space
        clear resid_fit Xfit_o
        %
        Qlimit_fit=limit(1);
        
        % calculate leverage
        hT_fit=diag(scores_fit*(inv(scores_fit'*scores_fit))*scores_fit'); %calculate MODE1 leverage
        
        % H lim_fit  according to Forina Chemo. Lab. 2009 
        limF_hTfit=betainv(1-alpha,nlv1/2,(DimXcal(1)-(nlv1-1))/2)*((DimXcal(1)-1)/DimXcal(1));
        %
        %%%% calculate Dof for h or SD according to A. Pomerantsev
        %%%% J.Chemom. 2008 by using IQR approximation, modified by me to
        %%%% get reasonable numbers i.e. :
        %%%% df=(exp(4.36ln(1.24/(IQR/ho)))^0.72
        dfAP_h=(exp(4.36*log(1.24/(iqr(hT_fit)/mean(hT_fit)))))^0.72;
        limAP_h=(mean(hT_fit)*chi2inv(1-alpha,dfAP_h))./dfAP_h;
        

        % model fit parameters
        ress_fit=sum(sum(sum(resid2_fit)));
        rmsec=sqrt(ress_fit/(DimXcal(1)-1));
        r2fit=1-ress_fit/ssX;
        %%
        % save space
        clear resid2_fit ress_fit
        %
        %%% save model and parameters for CALIBRATION SET in XfitResults
        
        XfitResults(f,ii).factors_fit=Factors_f;  % store factor loads
        XfitResults(f,ii).coreG=G;

        XfitResults(f,ii).Qfit=Qfit;     % store Q 
        XfitResults(f,ii).hT_fit=hT_fit; % store leverage
       
        %% store D_fit  disattivated
        %% XfitResults(f,ii).D_fit=D_fit; 
        %%
        %% XfitResults(f,ii).Dlim_fit=Dlim_fit; % store D limit
        
        XfitResults(f,ii).Qlimit_fit=Qlimit_fit;  % store residual limit
        
        XfitResults(f,ii).limF_hTfit=limF_hTfit;  % store hT limit for train (Forina)
       
        
        XfitResults(f,ii).Hlim_fitAP=limAP_h; % store h limit Pomerantsev
        XfitResults(f,ii).Qlimit_fitAP=limAP_Q;  % store Q limit Pomerantsev
        XfitResults(f,ii).QdfAP=dfAP_Q;  % store Q dof Pomerantsev
        XfitResults(f,ii).HdfAP=dfAP_h;  % store h dof Pomerantsev
        
        
        XfitResults(f,ii).RMSEC=rmsec;
        XfitResults(f,ii).R2=r2fit;
        
        if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
             XfitResults(f,ii).PossibleCombs=[];
        else %% Method TUCKER
             XfitResults(f,ii).PossibleCombs=possibleCombs(f,:);
        end
        %% temporary saving
        save tempfit XfitResults f ii
        %%

        
        
        %% APPLY MODEL on TEST SET (i.e. classes different from ii) %%%
        
        if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
            
            if ~isempty(Xtest)
                Fcal=Factors_f;Fcal{1}=rand(DimXtest(1),f);
                
                fac_new=parafac(Xtest_p,f,[1e-5 10 0 2 NaN maxit],const,Fcal,Fixm);
                Xnew=nmodel(fac_new);
                scoresnew=fac_new{1};
            end
            if ~isempty(Xtest_test)
                %% prediction of the inner test set
                
                Fcal_test=Factors_f;Fcal_test{1}=rand(DimXtest_test(1),f);
                fac_testnew=parafac(Xtest_test_p,f,[1e-5 10 0 2 NaN maxit],const,Fcal_test,[0 1 1]);
                Xnew_test=nmodel(fac_testnew);
                scoresnew_test=fac_testnew{1};
            end
            %%
            
            % save space
            if ~isempty(Xtest)
                clear Fcal
            end
            if ~isempty(Xtest_test)
                clear Fcal_test
            end
            %
               
        else %% Method TUCKER
            if ~isempty(Xtest)
                Xtest2=reshape(Xtest_p,DimXtest(1),prod(DimXtest(2:end)));
                %DimG=size(G);
                %G2=reshape(G,DimG(1),prod(DimG(2:end)));
                %Z=G2*kron(Factors_f{3},Factors_f{2})'; %ask Rasmus how to generalize to order>3
                
                Z=createZ(G,Factors_f);
                
                scoresnew=Xtest2*Z'*pinv(Z*Z');
                fac_new=Factors_f;
                fac_new{1}=scoresnew;
                Xnew=nmodel(fac_new,G);
            end
            if ~isempty(Xtest_test)
                %% prediction of the inner test set
                Xtest2_test=reshape(Xtest_test_p,DimXtest_test(1),prod(DimXtest_test(2:end)));
                %DimG=size(G);
                %G2=reshape(G,DimG(1),prod(DimG(2:end)));
                %Z=G2*kron(Factors_f{3},Factors_f{2})'; %ask Rasmus how to generalize to order>3
                
                Z=createZ(G,Factors_f);
                
                scoresnew_test=Xtest2_test*Z'*pinv(Z*Z');

                fac_new_test=Factors_f;
                fac_new_test{1}=scoresnew_test;
                Xnew_test=nmodel(fac_new_test,G);
            end

           
            % save space
            clear  scor1 scor2 Xtest2 
            %
        end
        if ~isempty(Xtest)
            Xnew3=reshape(Xnew,DimXtest);
        end
        if ~isempty(Xtest_test)
            Xnew3test=reshape(Xnew_test,DimXtest_test);
        end
        
        %%% rescale to original unit %%%%%
        if ord<=2
            if Pret(1,1:3)=='mnc'& Pret(2,1:3)~='blk';
                if ~isempty(Xtest); Xnew_o=rescal2D(Xnew3,Mx_cal, Sx_cal);end
                if ~isempty(Xtest_test); Xnew_test_o= rescal2D(Xnew3test,Mx_cal, Sx_cal);end
            elseif Pret(1,1:3)=='non' & Pret(2,1:3)=='non'
                if ~isempty(Xtest);Xnew_o=Xnew3;end
                if ~isempty(Xtest_test);Xnew_test_o=Xnew3test;end
            elseif Pret(1,1:3)=='mnc' & Pret(2,1:3)=='blk'
                if ~isempty(Xtest); Xnew_o=Xtest.*(ones(size(Xtest,1),1)*Sx_cal); Xnew_o=Xnew_o + ones(size(Xtest,1),1)*Mx_cal;end
                if ~isempty(Xtest_test); Xnew_test_o=Xtest_test.*(ones(size(Xtest_test,1),1)*Sx_cal); Xnew_test_o=Xnew_test_o + ones(size(Xtest_test,1),1)*Mx_cal;end
            elseif Pret(1,1:3)=='non' & Pret(2,1:3)=='blk'
                if ~isempty(Xtest); Xnew_o=Xtest.*(ones(size(Xtest,1),1)*Sx_cal);end
                if ~isempty(Xtest_test); Xnew_test_o=Xtest_test.*(ones(size(Xtest_test,1),1)*Sx_cal);end
            end
        else
            if ~isempty(Xtest); Xnew_o=pretX_miss_0(Xnew3, DimXtest, Pret, ordPret,indblk,Mx_cal, Sx_cal,-1);end % rescale Xnew to original units
            if ~isempty(Xtest_test); Xnew_test_o=pretX_miss_0(Xnew3test, DimXtest_test, Pret, ordPret, indblk,Mx_cal, Sx_cal,-1); end % rescale Xnew to original units
        end
        
        if ~isempty(Xtest); clear Xnew3; end  
        if ~isempty(Xtest_test); clear Xnew3test; end%% to save space
        %% end of rescaling
        if ~isempty(Xtest)
            hT_new=diag(scoresnew*(inv(scores_fit'*scores_fit))*scoresnew'); %leverage acording to A. Smilde
        end
        if ~isempty(Xtest_test)
            hT_new_test=diag(scoresnew_test*(inv(scores_fit'*scores_fit))*scoresnew_test');
        end
        %%
        if ~isempty(Xtest)
            resid2_test=(Xnew_o-Xtest).^2;
            if ord<=2
                new1=resid2_test';
                new2=sum(new1);
            else
                new1=permute(resid2_test,[3 2 1]);
                new2=sum(sum(new1));
            end
            Qnew=shiftdim(new2); % model residual for objects of other classes
        end
        if ~isempty(Xtest_test)
            resid2_test_test=(Xnew_test_o-Xtest_test).^2;
            if ord<=2
                new1_test=resid2_test_test';
                new2_test=sum(new1_test);
            else
                new1_test=permute(resid2_test_test,[3 2 1]);
                new2_test=sum(sum(new1_test));
            end
            Qnew_test=shiftdim(new2_test);
        end
        
        % pred parameters
        if ~isempty(Xtest)
            press=sum(sum(sum(resid2_test)));
            rmsep=sqrt(press/DimXtest(1));
            r2test=1-press/ssXtest;
        end
        % save space
        if ~isempty(Xtest)
            clear resid2_test Xnew_test new1 new2 press Xnew_o Xnew
        end
        if ~isempty(Xtest_test) 
            clear new1_test new2_test Xnew_test_o
        end 
        %
        % save model and parameters for TEST SET in XtestResults
        if ~isempty(Xtest); XtestResults(f,ii).scoresnew=scoresnew;end
        if ~isempty(Xtest_test); XtestResults(f,ii).scoresnew_test=scoresnew_test; end
        if ~isempty(Xtest); XtestResults(f,ii).hT_test=hT_new; end% store T2 predicted
        if ~isempty(Xtest_test); XtestResults(f,ii).hT_test_test=hT_new_test; end 
        
        %%
        if ~isempty(Xtest)
            XtestResults(f,ii).Qtest=Qnew; % store predicted model residuals
            XtestResults(f,ii).RMSEP=rmsep;
            XtestResults(f,ii).R2test=r2test;
        end
        if ~isempty(Xtest_test); XtestResults(f,ii).Qtest_test=Qnew_test; end% store predicted model residuals
 
         %% temporary saving
        if ~isempty(XtestResults); save tempTest XtestResults f ii ; end
        %%
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% MAKE SUB-MODELS and do CROSSVALIDATION  %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % initialize scoresCV
        scoresCV=[];
        X_CVpred=[];
        %
        %% CV to obtain predicted scores and residuals ( leaving entire objects Out )%%%%%%%%%%%
        SegmentsID=CVparam(ii).Segments;
        for iseg=1:size(SegmentsID,2)
            In = find(~SegmentsID(:,iseg));
            Out = find(SegmentsID(:,iseg));
            Xin=Xcal(In,:,:);  %  In samples
            Xout=Xcal(Out,:,:); % Out samples
            DimXin=size(Xin);
            DimXout=size(Xout);
            
            %% pretreatment for CV %%
            if ord<=2;
                if Pret(1,1:3)=='mnc' & Pret(2,1:3)~='blk'; [Xin_p, Mx_in, Sx_in]=scal2D(Xin,Pret,[],[]); Xout_p=scal2D(Xout,Pret,Mx_in, Sx_in);end
                if Pret(1,1:3)=='mnc' & Pret(2,1:3)=='blk'
                    Mx_in=missmean(Xin)
                    Xm=Xin - ones(size(Xin,1),1)*Mx_in;
                    ind_blk=indblk(2).indblk;
                    nblock=length(unique(ind_blk));
                    SStot=missss(Xm);
                    for iblk=1:nblock
                        indb=find(ind_blk==iblk);
                        SSblock(iblk)=missss(Xm(:,indb));
                        Sx_in(1,indb)=1./(sqrt(SStot./(SSblock(iblk).*nblock)));
                    end
                    Xin_p=Xm./(ones(size(Xin,1),1)*Sx_in);
                    Xout_p=Xout - ones(size(Xout,1),1)*Mx_in;
                    Xout_p=Xout_p./(ones(size(Xout,1),1)*Sx_in);                   
                    clear Xm
                end
                if Pret(1,1:3)=='non' & Pret(2,1:3)=='blk'
                    ind_blk=indblk(2).indblk;
                    nblock=length(unique(ind_blk));
                    SStot=missss(Xin);
                    for iblk=1:nblock
                        indb=find(ind_blk==iblk);
                        SSblock(iblk)=missss(Xin(:,indb));
                        Sx_in(1,indb)=1./(sqrt(SStot./(SSblock(iblk).*nblock)));
                    end
                    Xin_p=Xin./(ones(size(Xcal_p,1),1)*Sx_in);
                    Xout_p=Xout./(ones(size(Xout,1),1)*Sx_in);
                    Mx_in=zeros(1,size(Xin,2));
                end
                if Pret(1,1:3)=='non' & Pret(2,1:3)=='non'; Xin_p=Xin;Xout_p=Xout;Mx_in=zeros(1,prod(DimXin(2:end)));
                    Sx_in=ones(1,DimXin(2));
                end
            else
                [Xin_p, Mx_in, Sx_in]=pretX_miss_0(Xin, DimXin, Pret, ordPret, indblk);
                Xout_p=pretX_miss_0(Xout, DimXout, Pret, ordPret, indblk, Mx_in, Sx_in, 1);
            end
            %% end of pretreatment
            
            if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
                factCV = parafac(Xin_p,f,[1e-5 10 0 2 NaN maxit],const); % PARAFAC model for iseg CV Segments
                G_CV=[];
            else %% Method TUCKER
                [factCV, G_CV] = tucker(Xin_p,possibleCombs(f,1:end-1),[1e-5 0 0 2 NaN maxit]);
            end
            
            
            
            %Predict left-out samples
            if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
                %% ricordati di dare fixmode=[0 1] per biway
                Fcv=factCV;Fcv{1}=rand(DimXout(1),f);
                predCV=parafac(Xout_p,f,[1e-5 10 0 2 NaN maxit],const,Fcv,Fixm); 
                Xpred_out=nmodel(predCV);
                out_scores=predCV{1};
                % save space
                clear Fcv
                %
                G_CV=[];
            else %% Method TUCKER
                %out_scores = pred_tuck_score1(Xout_p, factCV, G_CV);
                Xout2=reshape(Xout_p,DimXout(1),prod(DimXout(2:end)));
                %DimGCV=size(G_CV);
                %G2CV=reshape(G_CV,DimGCV(1),prod(DimGCV(2:end)));
                %Z=G2CV*kron(factCV{3},factCV{2})'; %ask Rasmus how to generalize to order>3
                
                Z = createZ(G_CV,factCV);
               
                out_scores=Xout2*Z'*pinv(Z*Z');

                predCV=factCV;
                predCV{1}=out_scores;
                Xpred_out=nmodel(predCV,G_CV);

                % save space
%                 clear scor1 scor2 Xout_p_unf prodG_CV
            end
            Xpred_out3D=reshape(Xpred_out,DimXout);
            %save space
             clear Xpred_out
             %
             scoresCV(Out,1:size(out_scores,2))=out_scores;

             if ord<=2
                 if Pret(1,1:3)=='mnc' & Pret(2,1:3)~='blk'; Xpred_out_o=rescal2D(Xpred_out3D,Mx_in, Sx_in);end
                 if Pret(1,1:3)=='mnc' & Pret(2,1:3)=='blk'; Xpred_out_o=Xpred_out3D - ones(size(Xpred_out3D,1),1)*Mx_in; Xpred_out_o=Xpred_out_o./(ones(size(Xpred_out3D,1),1)*Sx_in);end
                 if Pret(1,1:3)=='non' & Pret(2,1:3)=='blk'; Xpred_out_o=Xpred_out3D./(ones(size(Xpred_out3D,1),1)*Sx_in);end                 
                 if Pret(1,1:3)=='non' & Pret(2,1:3)=='non'; Xpred_out_o=Xpred_out3D; end
             else
                 Xpred_out_o=pretX_miss_0(Xpred_out3D, DimXout, Pret, ordPret, indblk, Mx_in, Sx_in, -1); % rescale to original unit
             end
             if ord <=2
                 X_CVpred(Out,:)=Xpred_out_o;
             else
                 X_CVpred(Out,:,:)=Xpred_out_o;
             end
         end % end of CV loop
         
         %%
         save temp1 scoresCV f ii out_scores
         %%
         
         resid_CV=(X_CVpred-Xcal);
         resid2_CV=resid_CV.^2;
         if ord <=2
             Q_CV=sum(resid2_CV');
         else
             cv1=permute(resid2_CV,[3 2 1]);
             cv2=sum(sum(cv1));
             Q_CV=shiftdim(cv2);
         end
         hT_CV=diag(scoresCV*(inv(scores_fit'*scores_fit))*scoresCV'); %% analogous to what is done for test samples
         
         
         limitCV=residuallimit(resid_CV,1-alpha);
         Qlimit_CV=limitCV(1);
         
         %dof and limit for QCV and hT_CV based on A. Pomerentsev paper
         dfAP_QCV=(exp(4.36*log(1.24/(iqr(Q_CV)./mean(Q_CV)))))^0.72;
         limAP_QCV=(mean(Q_CV)*chi2inv(1-alpha,dfAP_QCV))./dfAP_QCV;
         dfAP_hT_CV=(exp(4.36*log(1.24/(iqr(hT_CV)./mean(hT_CV)))))^0.72;
         limAP_hT_CV=(mean(hT_CV)*chi2inv(1-alpha,dfAP_hT_CV))./dfAP_hT_CV;
 
         
         %% hT_CVlim, Qlim and Dlim_CV according to empirical 95 percentile on
         %% hT_CV and D_CV respectively
         [hTrlim_CV out_indH] = rlim(hT_CV,alpha);
         [Qrlim_CV out_indQ]  = rlim(Q_CV,alpha); % added Feb07
         
         

         % save model and parameters in CV mode in XvalResults
         XvalResults(f,ii).scores_CV=scoresCV;
         XvalResults(f,ii).Q_CV=Q_CV;
         XvalResults(f,ii).hT_CV=hT_CV;
         
         
         XvalResults(f,ii).Qlimit_CV=Qlimit_CV;
         %
         XvalResults(f,ii).limAP_QCV=limAP_QCV;
         XvalResults(f,ii).limAP_hT_CV=limAP_hT_CV;
         XvalResults(f,ii).dofAP_QCV=dfAP_QCV;
         XvalResults(f,ii).dofAP_hT_CV=dfAP_hT_CV;
         
         %
         XvalResults(f,ii).indout_hTrlim=out_indH;
         XvalResults(f,ii).indout_Qrlim=out_indQ;
                  
         
         % store outlier index for rlim calculation      
         
         XvalResults(f,ii).hTrlim_CV=hTrlim_CV;
         XvalResults(f,ii).Qrlim_CV=Qrlim_CV;
       
         % clear of temporary array
         clear cv1 cv2 resid_CV resid2_CV pressCV Xpred_out3D 
         % save space
         clear  Xin Xin_p Xout        
         %%
         save tempval XvalResults
         %%
    end %% end loop factors
    
        
    %%%
end %% end loop classes

 save  XfitResults   XfitResults
 save  XvalResults   XvalResults
 save  XtestResults  XtestResults



%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [SegmentsID,XvalMeth] = MakeSegments(I);
  
XvalMeth=questdlg('Which type of validation do you want to perform (ENTER => full Xval)?','Choose validation','Full X-validation','Segmented','Prespecified','Full X-validation');

switch XvalMeth
case 'Full X-validation'
   SegmentsID = speye(I);
   
case 'Segmented'
   prompt={'Enter the number of segments:'};
   eval(['def={''',num2str(min(I,max(3,round(I/7)))),'''};']);
   dlgTitle='Number of segments';
   lineNo=1;
   answer=inputdlg(prompt,dlgTitle,lineNo,def);
   NumbSegm=eval(answer{1});
   
   % Make sure the number of segments is OK
   while NumbSegm<2|NumbSegm>I
      prompt={'INCONSISTENT NUMBER CHOSEN (must be > 1 and <= samples)'};
      eval(['def={''',num2str(min(I,max(3,round(I/7)))),'''};']);
      dlgTitle='Number of segments';
      lineNo=1;
      answer=inputdlg(prompt,dlgTitle,lineNo,def);
      NumbSegm=eval(answer{1});
      NumbSegm<2|NumbSegm>I;
  end
   
   XvalSegm=questdlg('How should segments be chosen?','Choose segmentation','111222333...','123123123...','Random','123123123...');
   switch XvalSegm
      
   case '111222333...'
      SegmentsID = sparse(I,NumbSegm);
      NumbInEachSegm = floor(I/NumbSegm);
      Additional = I-NumbInEachSegm*NumbSegm;
      currentsample = 1;
      for i=1:NumbSegm
         if i <=Additional
            add = NumbInEachSegm+1;
         elseif i<NumbSegm
            add = NumbInEachSegm;
         else
            add = I-currentsample+1;
        end
         SegmentsID(currentsample:currentsample+add-1,i)=1;
         currentsample = currentsample + add;
     end
   case '123123123...'
      SegmentsID = sparse(I,NumbSegm);
      NumbInEachSegm = floor(I/NumbSegm);
      for i=1:NumbSegm
         SegmentsID(i:NumbSegm:end,i)=1;
     end
   case 'Random'
      % Make nonrandom and then randomize order
      SegmentsID = sparse(I,NumbSegm);
      NumbInEachSegm = floor(I/NumbSegm);
      for i=1:NumbSegm
         SegmentsID(i:NumbSegm:end,i)=1;
     end
      rand('state',sum(100*clock)) %Randomize randomizer
      [a,b] = sort(rand(I,1));
      SegmentsID = SegmentsID(b,:);
end
   
case 'Prespecified' 
      prompt={'Enter the name of the file defining the subsets'};
      def={'SegmentsID'};
      dlgTitle='Import definition';
      lineNo=1;
      answer=inputdlg(prompt,dlgTitle,lineNo,def);
      SegmentsID=eval(answer{1});
end % switch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    