function [T,P,U,R,C,B,V,ssq,eingen]=pls2sim(Xs,Ys, facts)
%   Function:
%   [T,P,U,R,C,B,V,ssq,eingen]=pls2sim(Xs,Ys,facts)
%   --------------------------------INPUT-------------------------------------------------
%   Xs    - matrix (n,p), predictor matrix (assumed to be center)
%   Ys    - matrix (n,m), predictand (assumed to be center)
%   facts - number of components to be tested
%
%   --------------------------------OUTPUT------------------------------------------------
%   B    - matrix (p,m), regression coefficients
%
%   --------------------------------Example-----------------------------------------------
%   [T,P,U,R,C,B,V,ssq,eingen]=pls2sim(Xs,Ys, 10)
%
%   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
%
%   Hugo da Silva Rocha, PhD Student
%   School of Chemistry, University of Manchester,
%   hugo.rocha@postgrad.manchester.ac.uk

    [ns,nx]=size(Xs);
    ny=size(Ys,2);
    
    S=(Ys'*Xs)';
    % when XtX=[] and X is very "tall", the booster XtX is calculated
    if ns>3*nx
        XtX=Xs'*Xs;
    else
        XtX=[];
    end			
    
	% initialization of variables
    T=zeros(ns,facts);
    U=T;					
    R=zeros(nx,facts);
    P=R;
    V=R;
    C=zeros(ny,facts);
    z=zeros(ny,1);
    
    ssqx=zeros(facts,1);
    ssqy=zeros(facts,1);
    
    ss2x=sum(sum(Xs.^2));
    ss2y=sum(sum(Ys.^2));
        
    StS=S'*S;				    
    nm1=ns-1;

    % SIMPLS algorithm
    for a=1:facts
        StS=StS-z*z';
        [Q,LAMBDA]=eig(StS);
        [lambda,j]=max(diag(LAMBDA));
        q=Q(:,j(1));
        r=S*q;
        t=Xs*r;
        
        if isempty(XtX)
            p=(t'*Xs)';
        else
            p=XtX*r;
        end
        
        if ns>nx
            d=sqrt(r'*p/nm1);
        else
            d=sqrt(t'*t/nm1);
        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; 						    % X weights
        P(:,a)=p/(d*nm1); 						% X loadings
        T(:,a)=t/d;							    % X scores
        U(:,a)=Ys*q;							% Y scores
        C(:,a)=q*(lambda(1)/(nm1*d)); 			% Y loadings
        eingen(:,a) = lambda;

        ssqx(a)=100*(sum(sum((T(:,a)*P(:,a)').^2))./ss2x);
        ssqy(a)=100*(sum(sum((U(:,a)*C(:,a)').^2))./ss2y);
        
    end
    B=R*C';
    
    ssq={ssqx ssqy};
end
