function [theta,lambda] = RegrHAH(X,N,epochs,gamma0,gamma00)

%   [theta,lambda] = RegrHAH(X,N,epochs,gamma0,gamma00)
%   [theta,lambda] = RegrHAH(X,N,epochs,gamma)
%   [theta,lambda] = RegrHAH(X,N)
%
% Hebbian - Anti-Hebbian Regression demonstration.
%
% Input parameters:
%  - X: Input data block (k x n)
%  - N: Number of latent variables to be extracted
%  - epochs: Number of iterations (default epochs=100)
%  - gamma0: Initial step size (default gamma0=0.1)
%  - gamma00: Final step size (default gamma00=gamma0)
% Return parameters:
%  - theta: Eigenvectors
%  - lambda: Sequence of eigenvalue vectors
%
% Heikki Hyotyniemi June 20, 2003


[k,n] = size(X);

if nargin < 3
   epochs = 100;
end
if nargin < 4
   gamma0 = 0.1;
end
if nargin < 5
   gamma00 = gamma0;
end
gamma = logspace(log10(gamma0),log10(gamma00),epochs);

[THETA,LAMBDA] = RegrPCA(X,n);

Rzz = eye(N);
Rzx = randn(N,n);
RT = zeros(N,N);
for i = 1:N
    RT(i,1:i) = ones(1,i);
end

lambda = NaN*ones(N,epochs);

for i = 1:epochs
   theta = (inv(Rzz)*Rzx)';
   Z = X*theta;
   Rzx = (1-gamma(i))*Rzx + gamma(i)*Z'*X/k;
   Rzz = (1-gamma(i))*Rzz + gamma(i)*Z'*Z/k;
   
   Rzz = RT.*Rzz;

   lambda(:,i) = diag(Rzz);
   clf
   hold on
   plot([1:i]',ones(i,1)*LAMBDA');
   plot([1:i]',lambda(:,1:i)');   
   plot([1:i]',lambda(:,1:i)','*');   
   title(['Behavior of the eigenvalues - epoch ',num2str(i)]);
   drawnow
end