
function [] = browseEigen(sizeIm, meanIm, basisIm, sigma, sampleIm)
  % [] = browseEigen(sizeIm, meanIm, basisIm, sigma, sampleIm)
  % Browse the eigen image basis.
    
  %%%%%%%
  % Plot the first 20 singular values: dQ plot
  %%%%%%%
  % Compute total scaled variance of training set.
  sigma2 = sigma .* sigma;
  totalVar = sum(sigma2);
  dQ = sigma2/totalVar;
  Q = cumsum(dQ);
  figure(1); clf;
  selectInd = 1:20;
  plot(selectInd, dQ(selectInd));
  hold on;
  plot(selectInd, dQ(selectInd), 'bo');
  hold off;
  xlabel('Singular value index, k');
  ylabel('Fraction of Variance, dQ');
  title('dQ(k): Variance Fraction Explained by one s.v.');
  % Note the rapid drop off in the singular values.  The first
  % singular direction explains about 20% of the total variance,
  % the second about 8 or 9% and so on... as follows:
  dQ(1:5) * 100
  fprintf(2, 'Press any key to continue...');
  pause; fprintf(2, 'ok\n');

  % Plot Q(k) = variance fraction explained by first k components
  figure(1); clf;
  selectInd = 1:20;
  plot(selectInd, Q(selectInd));
  hold on;
  plot(selectInd, Q(selectInd), 'bo');
  hold off;
  xlabel('Singular value index');
  ylabel('Fraction of Variance');
  title('Variance Fraction Explained by Subspace');
  % Note the first five singular directions account for almost half
  % the variance, as follows:
  Q(1:5) * 100
  fprintf(2, 'Press any key to continue...');
  pause; fprintf(2, 'ok\n');

  % Plot log(dQ(k)) over a large range of components
  % (all but the last two... which give artificially small
  % or zero variances).
  x = 1:(size(sigma,1)-2);
  figure(1); clf;
  plot(x, log10(sigma2(x)/totalVar));
  xlabel('Singular value index');
  ylabel('Log_{10}(Fraction of Variance)');
  title('Variance Fraction Explained by Subspace of Dimension k');
  % Note the singular values drop off rapidly at first, then
  % decrease roughly exponentially (linearly in this plot) for
  % a large range of k values.  This is typical for a wide
  % variety of natural stimuli.
  fprintf(2, 'Press any key to continue...');
  pause; fprintf(2, 'ok\n');

  %%%%%%%%%%%%%%
  % Plot the first nBasis basis images
  %%%%%%%%%%% 
  nBasis = 20;
  for k=1:nBasis
    figure(1); clf;
    showIm(reshape(basisIm(:,k), sizeIm));
    title(sprintf('Basis Image #%d',k));
    fprintf(2, 'Press any key to continue...');
    pause; fprintf(2, 'ok\n');
  end
  % Notice how the basis images for small values of k depict
  % larger scales, while the scale (of blobs of one sign, for example)
  % gets finer as k increases.


  %%%%%%%%%%%%%%
  % Orthogonality of basis images.
  %%%%%%%%%%% 
  % The following should be essentially zero:
  max(max(abs(basisIm' * basisIm - eye(size(basisIm,2)))))
  fprintf(2, 'Press any key to continue...');
  pause; fprintf(2, 'ok\n');

  %% Remove mean image from sample images
  nIm = size(sampleIm, 2);
  diffIm = sampleIm - repmat(meanIm(:), [1, nIm]);
  
  %%%%%%%%%%%%%%
  % Coefficient Histograms
  %%%%%%%%%%%%%%
  % Histograms of coefficients are very roughly Gaussian, mean zero,
  % and standard deviation sigma(k) for component k
  % A strong assumption of the detection model is that these
  % distributions can be reasonably described in terms of Gaussians.
  figure(1); clf;
  superimpose=false;  % superimpose=false; % shows individual histograms
  USE_LOG = true;
  for cDim=1:20
    coeffs = basisIm(:,cDim)' * diffIm/sigma(cDim);

    [hist bin] = histo(coeffs, 64);
    histDens = hist/nIm/(bin(2)-bin(1));
    minHistDens = 0.1/nIm/(bin(2)-bin(1));
    % This normalizes the density so that:
    %    sum( histDens ) * (bin(2)-bin(1))
    % is one.
    figure(1);
    if superimpose
      hold on;
    end
    if ~USE_LOG
      plot(bin, histDens);
      hold on;
      dx = 0.1;
      x = -6:dx:6;
      y = 1.0/sqrt(2.0*pi) * exp(-(x .* x)/2);
      plot(x, 1.0/sqrt(2.0*pi) * exp(-(x .* x)/2), 'g');
      % Here sum(y) * dx is one.
      lam = sqrt(2.0);
      plot(x, lam/2.0 * exp(-abs(x*lam)), 'r');
      axis([-8 8 0 0.8]);
      ylabel('Probability');
    else
      idx = histDens > 0;
      logHistDens = log10(minHistDens) * ones(size(histDens));
      logHistDens(idx) = log10(histDens(idx));
      plot(bin, logHistDens);
      hold on;
      dx = 0.1;
      x = -6:dx:6;
      y = 1.0/sqrt(2.0*pi) * exp(-(x .* x)/2);
      plot(x, log10(1.0/sqrt(2.0*pi) * exp(-(x .* x)/2)), 'g');
      % Here sum(y) * dx is one.
      lam = sqrt(2.0);
      plot(x, log10(lam/2.0 * exp(-abs(x*lam))), 'r');
      axis([-8 8 -4 1]);
      ylabel('log10 Probability');
    end
    xlabel('coeff(k)/sigma(k)');
    title(sprintf('Coeff histogram for comp. # %d',cDim));
    hold off;
    fprintf(2, 'Press any key to continue...');
    pause; fprintf(2, 'ok\n');
  end
  % From the plot the tails of the distribution appear to
  % be longer than the Gaussian model (green curve).  For example,
  % a Laplace distribution (red curve) might be at least as
  % appropriate.  Plotting log frequency instead of frequency
  % gives a better idea of the tails of the distributions: 
  %   Repaste above code using USE_LOG = true;


  % Show 2D scatter plots of the eigen image coefficients.
  useMean = true;
  for cDim = [2:10 20 40 80]
    figure(1); clf;
    plot(basisIm(:,1)' * diffIm, basisIm(:,cDim)' * diffIm, 'or'); 
    axis equal;
    title(sprintf('Regular eyes projected onto components 1 and %d', cDim));
    xlabel('Coeff. of 1st comp.');
    ylabel(sprintf('Coeff for comp. # %d', cDim));
    hold on;
    if useMean 
      theta = (-1:0.01:1)*pi;
    else
      theta = (-1/2:0.01:1/2)*pi;
    end
    x = sigma(1)*cos(theta);
    y = sigma(cDim) * sin(theta);
    plot(x,y,'g', 2*x, 2*y, 'b');
    grid on;
    hold off;
    cDim
    fprintf(2, 'Press any key to continue...');
    pause; fprintf(2, 'ok\n');
  end

  % Show reconstructions
  idx = randperm(nIm);
  nBasisTry = [5 20 50];
  for t = 1:10
    im = sampleIm(:,idx(t));
    figure(1); clf;
    subplot(2,2,1); showIm(reshape(im, sizeIm), [-3 3]);
    title(sprintf('Image %d', idx(t)));
    for j = 1:3
      nBasis = nBasisTry(j);
      im0 = im - meanIm;
      coeffs = basisIm(:,1:nBasis)' * im0;
      modelIm = meanIm + ...
                basisIm(:, 1:nBasis) * coeffs;
      figure(1); subplot(2,2,1+j);
      showIm(reshape(modelIm, sizeIm), [-3 3]);
      title(sprintf('Recon, nBasis=%d', nBasis));
    end
    fprintf(2, 'Press any key to continue...');
    pause; fprintf(2, 'ok\n');
  end



  % To get an idea of range of possible images that the PCA model
  % represents, we show some random images from the PCA model.  That
  % is, we chose nBasis Gaussian coefficients with the std devs given
  % by sigma.

  nBasisTry = [5 10 20 50 100];
  for t = 1:20
    % Generate enough random coefficients for all reconstructions
    coeffs = randn(nBasisTry(end), 1);
    coeffs = coeffs .* sigma(1:nBasisTry(end));
    figure(1); clf;
    for j = 1:5
      % Show the random image using just the first nBasis coeffs.
      nBasis = nBasisTry(j);
      modelIm = meanIm + basisIm(:, 1:nBasis) * coeffs(1:nBasis);
      figure(1); subplot(2,3,j);
      showIm(reshape(modelIm, sizeIm), [-3 3]);
      title(sprintf('Rand nBasis=%d', nBasis));
    end
    fprintf(2, 'Press any key to continue...');
    pause; fprintf(2, 'ok\n');
  end
end
