function [Vectors,Values,Psi] = pc_evectors(A,numvecs)
%PC_EVECTORS Get the top numvecs eigenvectors of the covariance matrix
%            of A, using Turk and Pentland's trick for numrows >> numcols
%            Returns the eigenvectors as the colums of Vectors and a
%            vector of ALL the eigenvectors in Values.

% Matthew Dailey 2000

  % Check arguments

  if nargin ~= 2
    error('usage: pc_evectors(A,numvecs)');
  end;

  nexamp = size(A,2);

  % Now compute the eigenvectors of the covariance matrix, using
  % a little trick from Turk and Pentland 1991

  % Compute the "average" vector
  % mean(A) gives you a row vector containing the mean of each column of A

  fprintf(1,'Computing average vector and vector differences from avg...\n');
  Psi = mean(A')';

  % Compute difference with average for each vector

  for i = 1:nexamp
    A(:,i) = A(:,i) - Psi;
  end;

  % Get the patternwise (nexamp x nexamp) covariance matrix

  fprintf(1,'Calculating L=A''A\n');
  L = A'*A;

  % Get the eigenvectors (columns of Vectors) and eigenvalues (diag of Values)

  fprintf(1,'Calculating eigenvectors of L...\n');
  [Vectors,Values] = eig(L);

  % Sort the vectors/values according to size of eigenvalue

  fprintf(1,'Sorting evectors/values...\n');
  [Vectors,Values] = sortem(Vectors,Values);

  % Convert the eigenvectors of A'*A into eigenvectors of A*A'

  fprintf(1,'Computing eigenvectors of the real covariance matrix..\n');
  Vectors = A*Vectors;

  % Get the eigenvalues out of the diagonal matrix and
  % normalize them so the evalues are specifically for cov(A'), not A*A'.

  Values = diag(Values);
  Values = Values / (nexamp-1);

  % Normalize Vectors to unit length, kill vectors corr. to tiny evalues

  num_good = 0;
  for i = 1:nexamp
    Vectors(:,i) = Vectors(:,i)/norm(Vectors(:,i));
    if Values(i) < 0.00001
      % Set the vector to the 0 vector; set the value to 0.
      Values(i) = 0;
      Vectors(:,i) = zeros(size(Vectors,1),1);
    else
      num_good = num_good + 1;
    end;
  end;
  if (numvecs > num_good)
    fprintf(1,'Warning: numvecs is %d; only %d exist.\n',numvecs,num_good);
    numvecs = num_good;
  end;
  Vectors = Vectors(:,1:numvecs);

