function model = pca_svd(X,ncomp,climit,opt)  
%   -------------------------------INPUT--------------------------------------
%   X        = A dataset with the structure: samples x variables        
%   ncomp    = Number of componets (spectra) to fit.        
%   -------------------------------OUTPUT-------------------------------------
%   Inside the structure "model" we have the following variables: 
%   T: scores
%   P: loadings
%   Var_Exp: percentage of explained variance
%   Cum_Var_Exp: percentage of explained variance accumulated
%   E: matriz de resíduos do modelo PCA
%   T2: T2 Hotelling [samples x variables]
%   Qres: Q residuals [sampels x variables]
% 
% 	References: 
%   Bro, R., & Smilde, A. K. (2014). Principal component analysis. 
%   Anal. Methods, 6(9), 2812–2831. doi:10.1039/c3ay41907j 
%
%   This is a part of the GNAT        
%   Copyright  2023  <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
%   mathias.nilsson@manchester.ac.uk
%%

[m,n] = size(X); 
[meanMtx, ~, ~]=prepfn(X, opt.prepr{1});

[u,s,v]= svd(meanMtx, 'econ');
Eigenvalue = diag(s).^2/(m-1);                                       
Var_Exp = 100*((Eigenvalue)/sum(Eigenvalue));  % Explained variance
Cum_Var_Exp = cumsum(Var_Exp);                 % Explained variance accumulated

T_all = u*s;                                   % Scores (T)
P_all = v;                                     % Loadings (P)
                                         
Eigenvaluencomp = Eigenvalue(1:ncomp);
Var_Exp = Var_Exp(1:ncomp);
Cum_Var_Exp = Cum_Var_Exp(1:ncomp);
T = T_all(:,1:ncomp);                                                
P = P_all(:,1:ncomp);                                                  

AxesPlot = {};
for i=1:ncomp
    AxesPlot{i,1} = sprintf('Scores on PC %d (%4g%%)',i,round((Var_Exp(i,1)),2));
    AxesPlot{i,2} = sprintf('Loading on PC %d (%4g%%)',i,round((Var_Exp(i,1)),3)); 
end

% Q-Residuals
Qres=sum((meanMtx-T*P').^2,2);
sres=sum(Qres/(m-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


% I = zeros(size(T,2),size(T,2));                         % Hotelling's T.^2 
% for i=1:size(T,2)
%     I(i,i) = Eigenvaluencomp(i);
% end
% 
% T2 = [];
% for i=1:size(T,1)
%     T2(i) = T(i,:)*((I)\T(i,:)');
% end

% Hotelling's T.^2                                                     
T2=diag((m-1)*T/(T'*T)*T'); 
T2lim = t2limit(m,ncomp,climit);

% Numerical labels
label = {};
for k=1:m
    label{k} = num2str(k);
end

a = ((Qlim - Qres))/Qlim;
OutQlim = sum(a(a>0));
InQlim = sum(a(a<0));

b = ((T2lim - T2))/T2lim;
OutT2 = sum(b(b>0));
InT2 = sum(b(b<0));

% Model structure
model.type = 'pca';                                     
model.samples = m;
model.variables = n;
model.mean_matrix = meanMtx;
model.var_exp = {Var_Exp};
model.var_cum = {Cum_Var_Exp};
model.Eigenvalue = Eigenvaluencomp;
model.Scores = T;
model.Loadings = P;
model.Qres = Qres';
model.Qlim = Qlim;
model.T2 = T2';
model.T2lim = T2lim;
model.AxesPlot = AxesPlot;
model.Label = label';
model.OutQlim = OutQlim;
model.InQlim = InQlim;
model.OutT2 = OutT2;
model.InT2 = InT2;
% model.R2 = R2;
% model.RMSE = RMSE;
