function [B,C,P,T,U,R,R2X,R2Y]=simpls(X,Y,A,S,XtX)

%   Function:
%   [B]=simpls(X,Y,h,S,XtX)
%   --------------------------------INPUT-------------------------------------------------
%   X    - matrix (n,p), predictor matrix (assumed to be center)
%   Y    - matrix (n,m), predictand (assumed to be center)
%   A    - scalar, number of PLS factors (Latent Variables)
%   S    - matrix (n,m), S=X'*Y
%   XtX  - matrix (n,n), XtX=X'*X (boosts speed for tall X matrices in which n>>p)
%
%   --------------------------------OUTPUT------------------------------------------------
%   B    - matrix (p,m), regression coefficients
%
%   --------------------------------REFERENCE---------------------------------------------
%
%
%
%   This is a part of the GNAT
%   Copyright  2017  <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, Undergrad Studant
%   Institute of Chemistry, Univeristy of Brasilia,
%   POBOX 70910-900, Asa Norte - Brasília (DF), Brazil
%   rocha.hugo@aluno.unb.br
%
%   Lattest updated 29/10/2021 HSR

[~,p] = size(X); 
[n,m] = size(Y);

if nargin<5
    S = []; 
end

if isempty(S)
    S=(Y'*X)'; 
end	

if nargin<4
    XtX=[]; 
end	

if isempty(XtX) && n>3*p
    XtX = X'*X; 
end	

if nargin<3
    A=10; 
end

A = min([A p n-1]);
T = zeros(n ,A); U = T;				
R = zeros(p,A); P = R; V = R;
C = zeros(m ,A); 
R2Y = zeros(1,A);
z = zeros(m,1); 
v = zeros(p,1);

if n>p
    S0 = S;
end
StS = S'*S;
nm1 = n-1;
tol = 0;

for a = 1:A
  StS = StS-z*z'; 
  [Q,LAMBDA] = eig(StS); 
  [lambda,j] = max(diag(LAMBDA)); 
  q = Q(:,j(1));
  r = S*q;
  t = X*r;
  if isempty(XtX)
      p = (t'*X)';
  else
      p = XtX*r;
  end
  if n>p
      d = sqrt(r'*p/nm1); 
  else 
      d = sqrt(t'*t/nm1); 
  end

  if d<tol
	disp(' ')
        disp('WARNING: the required number of factors (A) is too high')
	disp('Fewer PLS factors were extracted') 
	disp(' ')
	break,
  else 
      tol=max(tol,d/1e5);
  end

  v = p-V(:,1:max(1,a-1))*(p'*V(:,1:max(1,a-1)))'; v = v/sqrt(v'*v); 
  z = (v'*S)'; 
  S = S-v*z'; 
							
  V(:,a) = v;
  R(:,a) = r/d; 						
  P(:,a) = p/(d*nm1); 					
  T(:,a) = t/d;							
  U(:,a) = Y*q;						
  C(:,a) = q*(lambda(1)/(nm1*d)); 			
  R2Y(1,a) =  lambda(1)/d;			
end

clear StS V LAMBDA Q p q r t v z;

if d<tol
    A=a-1; 
    a=A; 
    T=T(:,1:A); 
    U=U(:,1:A); 
    R=R(:,1:A); 
    P=P(:,1:A); 
    C=C(:,1:A);
end

while a>1
    U(:,a) = U(:,a)-T(:,1:a-1)*(U(:,a)'*T(:,1:a-1)/nm1)'; 
    a=a-1; 
end

B = R*C'; % Regression Coeficients (Matrix)

if isempty(XtX) 
    sumX2=sum(X.^2); 
else 
    sumX2 = sum(diag(XtX)); 
end

R2X = 100*nm1/sum(sumX2)*cumsum(sum(P.^2)); 
R2Y = 100/nm1/sum(sum(Y.^2))*cumsum(R2Y(1:A).^2);



