function [Lambda,U,mu]=incremental_PCA3_step(xprim,p,alpha, U, Lambda, mu)
% One incremental step for dynamic PCA
% Reference: Y. Li, "On Incremental and Robust Subspace Learning", Pattern
% Recognition, 37, pp. 1509-1518, 2004.
% 
% Input: xprim   row vector of k attributes
%        p       number of reduced PCA dimensions
%        alpha   memory coefficient for the algorithm
%        U       matrix of eigenvectors from previous iteration
%        Lambda  (diagonal!)matrix of eigenvalues from previous iteration
%        mu      old values of means for each attribute   
% Output:
% Lambda: resulting matrix of eigenvalues (p most significant, not
% necessarily ordered)
% U: Contains new values of p most significant eigenvectors
% mu: new values of means for each attribute   

% D. Pokrajac, 7/7/2004
% 7/12/2004 >incremental_PCA Now, output is array of matrices during
% iterations
% 7/12/2004 >incremental_PCA3 performs only one STEP

    xprim=xprim'; %need be column vector
    k=length(xprim);
    x=xprim-mu; %equation 1
    mu=mu+(1-alpha)*x; %equation 2
    A=[U x]*[sqrt(alpha)*sqrt(Lambda) zeros(p,1);zeros(1,p) sqrt(1-alpha)];  %equations 3-5
    %Should compute correctly
    B=A'*A; %equation 9
    %Should we EXPLICITLY MAKE B such that only diagonal, last row and last
    %column are non-zero? We definitely can, and this is what I'm doing
    %here...
    Bb=diag(diag(B));
    Bb(:,p+1)=B(:,p+1);
    Bb(p+1,:)=B(p+1,:);
    B=Bb;
    [e_full,lambda_full]=eig(B); %equation 10 
    e_full=A*e_full; %equation 12 
    %I should take all but the smallest eigenvalue
    min(lambda_full);
    dd=diag(lambda_full);
    [dumb,index_min]=min(diag(lambda_full));
    Lambda=diag(dd([1:index_min-1 index_min+1:p+1])); %not necessarily sorted, now!
    %imagesc(Lambda);title(num2str(i));pause(0.1);
    U=e_full(:,[1:index_min-1 index_min+1:p+1]); %note: WHY this eigenvectors are not orthonormal???\
    
    U=U./repmat(sqrt(diag(U'*U))',k,1); %Let's try to orthonormalize...
    %Question: Do we need assumption of orthonormality???
    %Answer: It seems that we DO need to make U orthonormal. Otherwise, will not work, diverges!
    %MAybe could be corrected somehow (to embed nonorthonormal U above...)



return


