function [R,out,deni] = RTpca(X,nPCA,nit,T)
%% Inputs :
% X : Input Matrix, each variable is fed into one row and each single
% column denotes a single datapoint
% nPCA : Number of required PCs to estimate 
% nit : number of iterations
%% Outputs :
% R : EigenVectors of each PC, each row is a complete eigenvector.
% out : Estimated PCs in realtime; The las row is an outcome made by the
% wiegth of the first three PCs.
% deni : eigenvalues
%% Method
% The real time PCA (RTpca) is a brilliant approach to estimate the PCs in
% real time with the help of Hebian learning rule and iterative learning
% with a single neuron. This method is first inspired and motivated by Oja
% et al but the current algorithm is an out of the work by Principe et al.
% The convergence speed is determined by T. This function works perfect for 
% estimating the first 4 PCs.

%% Reference
% [1] H. Dongho; Y.N. Rao, J.C. Principe, K. Gugel.
% ‘’ Real-Time PCA(Principa1 Component Analysis)’’. IEEE International Joint Conference
% on Neural Networks (IEEE Cat. No.04CH37541) , 2004, p2159-2159, 1p. Publisher: IEEE

%% Author : Hooman Sedghamiz     hooman.sedghamiz@gmail.com
% Copyright August 2014

%%
if nargin < 4
    T = 0.90;                                                                  % default is 0.90
    if nargin < 3
       nit = 1;
       if nargin < 2
          nPCA = size(X,1); 
       end
    end
end

% get the dimensionality
[m,n] = size(X);


% ----------------- random initial weights if no input --------------- %
w = rand(m,nPCA);


% -------- T determines convergance (higher faster convergence)------- %
out = zeros(nPCA,n);
numi = zeros(m,nPCA);                                                      %first iteration always zero
deni = zeros(nPCA,1);                                                      %first iteration always zero
y = zeros(nPCA,1);
R = zeros(m,nPCA);

%dummi =[];

 for iter = 1:nit
  XX = X;       
  for ii = 1:n
        for r = 1 : nPCA 
          %% =================== Network ========================= %%
          y(r) = w(:,r)'*XX(:,ii); 
          % update network nominator
          numi(:,r) = (1-(1/ii))*numi(:,r) + (1/ii)*(XX(:,ii)*y(r)); 
          deni(r) = (1-(1/ii))*deni(r) + (1/ii)*y(r)^2;
          w(:,r) = (1 - T)*w(:,r) + T*(numi(:,r)/deni(r));  
          out(r,ii) = w(:,r)'*XX(:,ii);
          %% ========== Deflation to compute the remining PCs =========%%
          deflation_s = (w(:,r)*y(r));
          XX(:,ii) = XX(:,ii) - deflation_s;
          %% =============== Check for convergance =================== %%
          if R(:,r) == w(:,r)
             fprintf('Converged! in %d iteration and %d Sample',iter,ii);    
          end
          %% =============== Save the EigenVector ==================== %%
          R(:,r) = w(:,r);      
        end
        %dummi =[dummi deni(:)]; %length of eigenvectors
        %M = deni(1:end)./sum(deni(:));
        %out(r+1,ii) = M(3)*X(1,ii) + M(2)*X(2,ii) + M(1)*X(3,ii);
  end
  
 end
end