function [T,P,U,R,W,Q,B,ssq]=pls2nip(X,Y, facts)    
    [ns,npx]=size(X);
    [~,npy]=size(Y);
    
    ssqx=zeros(facts,1);
    ssqy=zeros(facts,1);
    
    ss2x=sum(sum(X.^2));
    ss2y=sum(sum(Y.^2));

    % Set E0=X and F0=Y
    E=X; 
    F=Y; 
    W=zeros(npx,facts); T=zeros(ns,facts); P=zeros(npx,facts); C=zeros(1,facts);
    Q=zeros(npy,facts); U=T;   % Initialization of the matrices

    % Loop over the number of LVs
    for i=1:facts                         
        [~,b]=max(std(F));      % Initial estimate of u as the Y(F) column with the largest variance
        u_old=F(:,b);
        crit=1;                 % initializing the convergence criterion
        while crit>10^(-9)      % loop for the computation of the ith component
            w=E'*u_old;
            w=w./sqrt(w'*w);
            t=E*w;
            q=F'*t; q=q./norm(q);
            u=F*q;
            crit=norm(u-u_old)./norm(u_old);  % check for convergence.
            u_old=u;
        end
        
        p=E'*t./(t'*t); c=u'*t./(t'*t);
        W(:,i)=w; T(:,i)=t; P(:,i)=p; C(i)= c; Q(:,i)=q; U(:,i)=u;
        E=E-t*p'; F=F-c*t*q';  % deflation of X and Y
        
        
        ssqx(i)=100*(sum(sum((T(:,i)*P(:,i)').^2))./ss2x);
        ssqy(i)=100*(sum(sum((U(:,i)*Q(:,i)').^2))./ss2y);
        
    end

    % R weights (allow direct calculation of T from X); 
    R=W/(P'*W);    

    % Calculation of PLS regression coefficients
    B=R*diag(C)*Q';
    ssq={ssqx ssqy};
end
