function model = doPLSDA(X,Y,Factor,opt,Classes)
%   % Function to perform PLS-DA calculation
%   Function:
%   model = doPLSDA(X,Y,Factor,opt)
%   --------------------------------INPUT-------------------------------------------------
%   X           Matrix (n,p), predictor matrix (assumed to be center)
%   Y           Matrix (n,m), predictand (assumed to be center)
%   Factor      number of Latent variables used for PLS-DA model calculation
%   opt         options for PLS-DA model calculation 
%        opt.prepr - Preprocessing method for the X and Y method:
%               {'none'};
%               {'mean'};
%               {'auto'};
%               {'pareto'};
%        opt.algorithm - PLS-DA calculation method:
%               'simpls'
%               'nipals'
%        opt.detailedoutput - compute the T2 and Q-residual values:
%               'on'
%               'off'
%        opt.classcriterion - defines how classification is accomplished based on the values of the predicted Y. 
%               'maxY'   (samples are assigned based on the highest component of the predicted Y) 
%               'lda'    (samples are assigned after LDA on predicted Y or X scores.       
%        opt.ldax - parameters for the lda prediction method
%               'predy'  (LDA is made on PLS predicted Y after removale of one column)
%               'scores' (LDA is made on PLS scores T)
%        opt.prior - 
%               'uniform' (sequential order)
%               'frequency' (prior based on the class frequencies in the training set')
%   Classes  
%              Vector array with relation of classes for the calibration set 
%   --------------------------------OUTPUT------------------------------------------------
%   model       main structure with information of the PLS-DA model
%
%   --------------------------------References--------------------------------------------
% 
%    [1] de Jong, S. (1993) "SIMPLS: an alternative approach to partial least squares
%        regression", Chemometrics and Intelligent Laboratory Systems, 18:251-263.
%    [2] Rosipal, R. and N. Kramer (2006) "Overview and Recent Advances in Partial
%        Least Squares", in Subspace, Latent Structure and Feature Selection:
%        Statistical and Optimization Perspectives Workshop (SLSFS 2005),
%        Revised Selected Papers (Lecture Notes in Computer Science 3940), C.
%        Saunders et al. (Eds.) pp. 34-51, Springer.
%
%   This is a part of the GNAT
%   Copyright  2024  <Mathias Nilsson>%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%
%   Hugo da Silva Rocha, PhD Student
%   School of Chemistry, University of Manchester,
%   hugo.rocha@postgrad.manchester.ac.uk

    [~,ycl]=max(Y,[],2);
    nclass=size(Y,2);
%     Y = vet_matrix(Y,size(Y,2)); % array of classes

    % Preprocessing step for X and Y
    [Xp, mx, sx]=prepfn(X, opt.prepr{1});    %Preprocessing of X matrix
    [Yp, my, sy]=prepfn(Y, opt.prepr{1});    %Preprocessing of Y matrix
    [ns,nx]=size(X);
    facts=min([Factor nx ns-1]);  

    % Method for PLS-DA calculation
    switch opt.algorithm
        case 'nipals'
            [T,P,U,R,W,Q,B,ssq]=pls2nip(Xp,Yp, facts);
            model.scores={T U};
            model.loadings={P Q};
            model.weights=W;
            model.Rweights=R;
            model.explvar=ssq;
            model.regcoef=B;
            [vip, viptot]=myvip2(T,W,R,Q,B);  
%             vip_scores = calculate_vip(Xp, Yp, W, T, Q);
        case 'simpls'
            [T,P,U,R,Q,B,V,ssq,lambda]=pls2sim(Xp,Yp, facts);          
            model.scores={T U};
            model.loadings={P Q};
            model.Rweights=R;
            model.basis=V;
            model.explvar=ssq;
            model.regcoef=B;
            model.eigvalue=lambda;
            [vip, viptot]=myvip2(T,[],R,Q,B);  
%             vip_scores = calculate_vip(Xp, Yp, R, T, Q);
    end
    
    AxesXPlot = {};
    for i=1:Factor
        AxesXPlot{i,1} = sprintf('Scores on LV %d',i);
    end
    model.AxesXPlot= AxesXPlot;
    Yhat=Xp*B;
    Ypred=unprepfn(Yhat,opt.prepr{1},my,sy);

    % Backscale the Loadings and VIP score when using Autoscale
    switch opt.prepr{1}
        case 'auto'
            viptot=viptot.*sx';
            model.loadings{1,1}=model.loadings{1,1}.*sx';
    end

    model.classes=Y;
    model.preds=Ypred;
    model.YResiduals=Y-Ypred;
    model.bias=mean(model.YResiduals);
    model.RMSE=sqrt(sum(model.YResiduals.^2)/size(Ypred,1));

    % R2 Calculation
    % Training R2
    RSS=sum(sum(model.YResiduals.^2)); % Sum of Squared Errors for Test Set
    TSS=sum(sum((Y-repmat(mean(Y), size(Y,1), 1)).^2)); % Total Sum of Squares for Test Set
    model.R2=1-(RSS./TSS(1));
            
    model.vip=vip;
    model.TotVip=viptot;

    model.raw_matrix = X;
    model.mean_matrix = Xp;
        
    if strcmp(opt.detailedoutput,'on')        
        q=sum((Xp-T*P').^2,2);
        t2=diag((ns-1)*T/(T'*T)*T'); %tsquare
        t2lim=facts*(ns-1)*finv(0.95,facts,ns-facts)./(ns-facts);           
        sres=sum(q/(ns-1));        
        % Computation of qlim using J-M approx
        theta1 = sum(sres);
        theta2 = sum(sres.^2);
        theta3 = sum(sres.^3);
        if theta1==0
            qlim = 0;
        else
            h0     = 1-2*theta1*theta3/3/(theta2.^2);
            if h0<0.001
                h0 = 0.001;
            end
            ca    = sqrt(2)*erfinv(2*0.95-1);
            h1    = ca*sqrt(2*theta2*h0.^2)/theta1;
            h2    = theta2*h0*(h0-1)/(theta1.^2);
            qlim = theta1*(1+h1+h2).^(1/h0);
        end        
        model.tsq=t2;
        model.qres=q;
        model.tsqlim=t2lim;
        model.qlim=qlim;
        model.tsqr=t2/t2lim; 
        model.qr=q/qlim;           
    end
            
    model.preprpars={mx sx; my sy};
    model.options=opt;

    %folding the classes into a vector
    ycalV = vet_matrix(Y,size(unique(Classes),1));

    % class probabilities, Chemometrics and Intelligent Laboratory Systems 95 (2013) 122-128.
    for g=1:max(ycalV)
        mc(g) = mean(Ypred(find(ycalV==g),g));
        sc(g) = std(Ypred(find(ycalV==g),g));
        mnc(g) = mean(Ypred(find(ycalV~=g),g));
        snc(g) = std(Ypred(find(ycalV~=g),g));
        for i=1:size(X,1)
            Pc = 1./(sqrt(2*pi)*sc(g)) * exp(-0.5*((Ypred(i,g) - mc(g))/sc(g)).^2);
            Pnc = 1./(sqrt(2*pi)*snc(g)) * exp(-0.5*((Ypred(i,g) - mnc(g))/snc(g)).^2);
            prob(i,g) = Pc/(Pc + Pnc);
        end
    end

    % class evaluation
    [~,class_true] = max(Y');
    resthr = plsdafindthr(Ypred,class_true');
    switch opt.classcriterion
        case 'maxY'
            % assigns on the maximum calculated response
            [~,assigned_class] = max(prob');
        case 'threshold'
            % assigns on the bayesian discrimination threshold
            assigned_class = plsdafindclass(Ypred,resthr.class_thr);
        case 'lda'
            if strcmp(opt.ldax, 'scores')
                xx=model.scores{1};
            elseif strcmp(opt.ldax, 'predy')
                xx=model.preds;
                xx=xx(:,1:end-1);
            end
            
            mlda = ldaf(xx,xx,ycl, opt.prior);
            model.classification.PredClass=mlda.PredClass;
            model.classification.probability=mlda.PosteriorProb;
    end
    class_param = calc_class_param(assigned_class',class_true');

    model.classification.TrueClass=ycl;
    switch opt.classcriterion
        case 'maxY'
            [~,YPredClass]=max(model.preds,[],2);
            model.classification.PredClass=YPredClass;
        case 'lda'
            if strcmp(opt.ldax, 'scores')
                xx=model.scores{1};
            elseif strcmp(opt.ldax, 'predy')
                xx=model.preds;
                xx=xx(:,1:end-1);
            end
            
            mlda = ldaf(xx,xx,ycl, opt.prior);
            model.classification.PredClass=mlda.PredClass;
            model.classification.probability=mlda.PosteriorProb;
    end
    
    model.threshold = resthr;
    model.class_param = class_param;
    cmatrix=cmatrixcalc(model.classification.PredClass, ycl, nclass);
    model.classification.ConfusionMatrix=cmatrix;
    ccrate(1:size(cmatrix,2))=100*diag(cmatrix)./sum(cmatrix,2);
    model.classification.CorrClassRate=ccrate;
    model.classification.CorrClassTot=100*sum(diag(cmatrix))/length(ycl);        
    model.options=opt;             
end