clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  In this demo, we'll design a synthetic data set with a standing wave
%  in it, and perform EOF analysis on that data set.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Start by defining the space and time dimensions:

[xx, tt] = meshgrid(0:10, 1:100);   %  11 points in the xx-dimension
                                    %  (state), 100 points in the
                                    %  tt-dimension (sampling). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define standing wave

data = sin(2*pi*xx/10) .* cos(2*pi*tt/20);

%  Let's have a look at the resulting data

figure(1);
orient tall;
wysiwyg

subplot(2,2,1); cla;

  %  Make positive contours solid, negative contours dashed, the zero
  %  contour thicker:
  [cp, pp] = contour(xx, tt, data, [0.2:.2:1], '-k');
  hold on;
  [cn, pn] = contour(xx, tt, data, [-1:.2:-.2], '--k');
  [cz, pz] = contour(xx, tt, data, [0 0], '-k');
  set(pz, 'linewidth', 2);
  hold off;
  
  axis([0 10 1 100]);
  t1 = title('Raw Data', 'fontsize', 12);
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Now, let's add some gaussian noise to this data

data2 = data + randn(100, 11);

%  Again, let's have a look at the resulting data
subplot(2,2,2); cla;

  %  Make positive contours solid, negative contours dashed, the zero
  %  contour thicker:
  [cp, pp] = contour(xx, tt, data2, [0.5:.5:3], '-k');
  hold on;
  [cn, pn] = contour(xx, tt, data2, [-3:.5:-.5], '--k');
  [cz, pz] = contour(xx, tt, data2, [0 0], '-k');
  set(pz, 'linewidth', 2);
  hold off;
  
  axis([0 10 1 100]);
  t1 = title('Noisy Data', 'fontsize', 12);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Now, let's apply EOF analysis to the noisy data.  We'll do this using
%  the SVD method (see eof_example.m for examples of both the SVD method
%  and the 'eigenvectors of the covariance matrix' method.

[u, s, v] = svd(data2);
lam = diag(s.^2)/99;
pcs = u * s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Plot the eigenspectrum

subplot(4,2,5);
  errorbar(1:11, lam, lam*sqrt(2/99), 'ok');
  set(gca, 'XTick', 1:11);
  t1 = title('Eigenspectrum', 'fontsize', 12);
  grid on;
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Plot the EOFs, scaled by the square root of the eigenvalue.  

eof1 = v(:,1)*sqrt(lam(1));
eof2 = v(:,2)*sqrt(lam(2));
eof3 = v(:,3)*sqrt(lam(3));

%  Note, we coud have done the exact same thing using the following
%  commands:
%
%  eofs = (s/sqrt(99)) * v';
%  eofs = eofs(1:11, 1:11)';
%
%  Now, the first column of 'eofs' is eof1, the second eof2, etc ...
%

%  Plot the EOFs
subplot(4,2,6); cla;
  pp = plot(0:10, eof1, '-k', 0:10, eof2, '-r', 0:10, eof3, '-b');
  %  Let's make EOF1 more visible:
  set(pp(1), 'linewidth', 2);
  
  %  For kicks, let's also plot our constructed data:
  hold on;
  pp2 = plot(0:10, sin(2*pi*[0:10]/10), '--k');
  set(pp2, 'linewidth', 2);
  hold off;
  
  axis([0 10 -1 1]);
  t1 = title('EOFs 1-3', 'fontsize', 12);
  grid on;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Finally, plot the standardized PC's
pc1 = pcs(:,1)/sqrt(lam(1));
pc2 = pcs(:,2)/sqrt(lam(2));
pc3 = pcs(:,3)/sqrt(lam(3));

%  How could we have accomplished the above set of commands using linear
%  algebra?

subplot(4,1,4); cla;
  pp = plot(1:100, pc1, 'k', 1:100, pc2, 'r', 1:100, pc3, 'b');
  %  Let's make pc1 more visible:
  set(pp(1), 'linewidth', 2); 
  
  %  Plot theoretical standardized (hence the sqrt(2)) PC1:
  hold on;
  pp2 = plot(1:100, cos(2*pi*[1:100]/20)*sqrt(2), '--k');
  set(pp2, 'linewidth', 2);
  hold off;
  
  axis([1 100 -3 3]);
  grid on;
  t1 = title('PCs 1-3', 'fontsize', 12);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Done!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Here are some other things to try:
%
%  1.  Try a propagating wave instead of a standing wave.
%  2.  Change the number of samples (time) to see what happens to the
%      separation between the eigenvalues, and the noisiness of the EOFs
%      and PCs.
%  3.  Plot the unit-length EOFs and the original PCs (so that var(pcs) =
%      lam).  What plotting convention is more useful in this case?
%  4.  Calculate the EOF's, eigenvalues, and PC's by calculating the
%      eigenvectors of the covariance matrix.  Make sure that the plots
%      you produce are identical to the ones produced here (note that if
%      you use a different noise realization, they won't be identical,
%      but the amplitude and such should be appropriate).