function [T, W, P, B, Q, a] = nipals(X,y,h)

%-----------------------------------------------------------------
%E0 = X; F0 = y;
[n,m] = size(X);
[E0, F0] = scale(X,y);
%----------------------------------------------------------------
u(:,1) = F0;

W(:,1) = (E0'*u(:,1))/(u(:,1)'*u(:,1));
W(:,1) = W(:,1)/norm(W(:,1));

T(:,1) = (E0*W(:,1));

C(:,1) = (F0'*T(:,1))/(T(:,1)'*T(:,1));
C(:,1) = C(:,1)/norm(C(:,1));

u(:,1) = (F0*C(:,1))/(C(:,1)'*C(:,1));

P(:,1) = (E0'*T(:,1))/(T(:,1)'*T(:,1));
Q(:,1) = (F0'*T(:,1))/(u(:,1)'*u(:,1));

b(1) = (u(:,1)'*u(:,1))/(T(:,1)'*T(:,1));

E(:,:,1) = E0-T(:,1)*P(:,1)';
F(:,1) = F0-b(1)*T(:,1)*Q(:,1)';

for i = 2:h
    u(:,i) = F(:,(i-1));
    
    W(:,i) = (E(:,:,(i-1))'*u(:,i))/(u(:,i)'*u(:,i));
    W(:,i) = W(:,i)/norm(W(:,i));
    
    T(:,i) = (E(:,:,(i-1))*W(:,i));
    
    C(:,i) = (F(:,i-1)'*T(:,1))/(T(:,i)'*T(:,i));
    C(:,i) = C(:,i)/norm(C(:,i));
    
    u(:,i) = (F(:,(i-1))*C(:,i))/(C(:,i)'*C(:,i));

    P(:,i) = (E(:,:,(i-1))'*T(:,i))/(T(:,i)'*T(:,i));
    Q(:,i) = (F(:,i-1)'*T(:,i))/(u(:,i)'*u(:,i));
    
    b(i) = (u(:,i)'*u(:,i))/(T(:,i)'*T(:,i));

    E(:,:,i) = E(:,:,(i-1))-T(:,i)*P(:,i)';
    F(:,i) = F(:,(i-1))-b(i)*T(:,i)*Q(:,i)';
end

B = diag(b);

%-------------------------------------------------------------------

Ws(:,1) = W(:,1);

for j = 2:h
    w = eye(m);
    for k = 1:(j-1)
        w = w*(eye(m)-W(:,k)*P(:,k)');
    end
    
    Ws(:,j) = w*W(:,j);
end

a = Ws*B*Q';

%--------------------------------------------------------------------