function [T,P,U,W,Q,B,V,ssq,eingen]=pls2sim2(X,Y, ncomp)
%   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

    [n, p] = size(X);
    m = size(Y, 2);
    
    % Initialize matrices
    T = zeros(n, ncomp);
    U = zeros(n, ncomp);
    P = zeros(p, ncomp);
    Q = zeros(m, ncomp);
    W = zeros(p, ncomp);
    V = zeros(p, ncomp);
    
    % Calculate covariance matrix
    S = X' * Y;
    
    % Precompute XtX if beneficial for performance
    if n > 3 * p
        XtX = X' * X;
    else
        XtX = [];
    end
    
    ssqX = zeros(ncomp, 1);
    ssqY = zeros(ncomp, 1);
    
    totalSSX = sum(X(:).^2);
    totalSSY = sum(Y(:).^2);
    
    nm1 = n - 1;
    z = zeros(m, 1);
    
    % SIMPLS algorithm
    for a = 1:ncomp
        % Deflation for components after first
        if a > 1
            StS = StS - z * z';
        else
            StS = S' * S;
        end
        
        % Find dominant eigenvector
        [eigVecs, eigVals] = eig(StS);
        lambda = diag(eigVals);
        [lambda_max, idx] = max(lambda);
        q = eigVecs(:, idx);
        
        % Calculate weights
        r = S * q;
        t = X * r;
        
        % Calculate X loadings
        if isempty(XtX)
            p_vec = X' * t;
        else
            p_vec = XtX * r;
        end
        
        % Normalize
        if n > p
            d = sqrt(r' * p_vec / nm1);
        else
            d = sqrt(t' * t / nm1);
        end
        
        % Calculate Y scores
        u = Y * q;
        
        % SIMPLS deflation
        v = p_vec;
        if a > 1
            v = v - V(:, 1:a-1) * (V(:, 1:a-1)' * p_vec);
        end
        v = v / sqrt(v' * v);
        z = (v' * S)';
        S = S - v * z';
        
        % Store components
        V(:, a) = v;
        W(:, a) = r / d;                        % X weights
        P(:, a) = p_vec / (d * nm1);            % X loadings  
        T(:, a) = t / d;                        % X scores
        U(:, a) = u;                            % Y scores
        Q(:, a) = q * (lambda_max / (nm1 * d)); % Y loadings
        
        % Calculate explained variance
        ssqX(a) = 100 * (sum(sum((T(:, a) * P(:, a)').^2)) / totalSSX);
        ssqY(a) = 100 * (sum(sum((U(:, a) * Q(:, a)').^2)) / totalSSY);
    end
    
    % Calculate regression coefficients
    % B=W*Q';
    % B = W * (P' * W) \ Q';
    B = W * (P' * W) / Q;
    ssq={ssqX ssqY};
    
    % Calculate eigenvalues
    eingen = sum(T.^2, 1)';
end

