%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  EOF / PC analysis example:  EOF / PC analysis applied to 
%  a matrix image.
%
%  In this case, the bitmap-style image is analagous to our 3D data set,
%  with the column dimension representing time, the row dimension
%  representing one spatial dimension (e.g. lat), and the color dimension
%  representing the other spatial dimension (e.g. lon).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Clear the workspace
clear


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Part 1:  Load and format the data:

%  Load the data - in this case, a picture of a Motu off Muri beach in
%  Rarotonga, Cook Islands.
cd /Users/dvimont/matlab/Data
x = imread('cooks2.jpg');

%  Note that 'x' is an 8-bit integer.  Convert to double precision, so
%  MATLAB routines work:
y = double(x); 

%  Our data is three-dimensional:  column by row by color (3-elements).
%  Note that this is analagous to a 3D data set, with time, lat, lon as
%  the three dimensions.  We need to reshape so we have time (column) by
%  space (row):
[ntim,nlat,nlon]=size(y);
y=reshape(y,ntim,nlat*nlon);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Part 2:  Calculate EOF's of the data:

%  Calculate EOF's of the data.  Note, we have not removed the mean
%  (ordinarily, you'd want to remove the mean), and we're calculating the
%  covariance matrix across the 'spatial' dimension (why do we do that,
%  rather than the 'temporal' dimension?)
[pcs,lam]=eig(y*y'/(nlat*nlon - 1));

%  Now, reorder the eigenvalues and PC's in decreasing amplitude (of the
%  eigenvalues):
lam=fliplr(flipud(lam));
pcs=fliplr(pcs);

%  Calculate percent variance explained:
per=diag(lam)/trace(lam);

%  Calculate EOF's by regressing data onto the PC's:
lds=(pcs'*y)';

%  Now, let's rescale all of this so that var(pcs) = lam, and lds are
%  orthonormal:
wgt = diag(1./(sqrt(diag(lam)*(nlat*nlon-1))));
lds = lds*wgt;

wgt = diag(sqrt(diag(lam)*(nlat*nlon-1)));
pcs = pcs*wgt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Part 3:  Display the results in a meaningful way

%  Set up a figure window
figure_tall(1); clf;

%  Set up the subplot geometry:
global_axes(6.4*.6, 4.8*.6, .5, .75, 1.5);

%  Plot the original picture:
subplot2(1,1);
  image(uint8(reshape(y, ntim, nlat, nlon)));
  t1 = title('\bfOriginal Picture', 'fontsize', 14);


%%%%%%%%%%%%%%%%%%%% Plot reconstructed picture %%%%%%%%%%%%%%%%%%
%  Reconstruct the picture, using a subset of the PCs:
last_eof = 1;  %  This is the last EOF to use in the reconstruction.  You
               %  should change this around.
ind = 1:last_eof;
y2 = pcs(:,ind)*lds(:,ind)';          %  Reconstruct by summing over EOF's
                                     %  and PC's
y2 = reshape(y2, ntim, nlat, nlon);  %  Reshape back to picture format

%  Plot the reconstructed data
global_axes(6.4*.6, 4.8*.6, .5, .75, 1.5);
subplot2(1,2)
image(uint8(y2));
title(['\bfModes 1 through ' num2str(last_eof) ':  ' ...
       num2str(round(10000*sum(per(ind)))/100) '% Var Expl.'], ...
      'fontsize', 14);
  
%  Plot the Principal Components:
global_axes(1.25, 4.8*.6, 6.4*.6+.5, .75, 1.5);
subplot2(2,4); cla;
pp = plot(pcs(:,last_eof), ntim:-1:1, 'k');

%  If we're plotting the first PC, make it look good (remember, we
%  haven't removed the mean)
if last_eof == 1; set(gca, 'XLim', [0 8000]); end;

%  Make it pretty:
set(pp, 'linewidth', 2);
set(gca, 'YTickLabel', [], 'YLim', [1 ntim]);
grid on;
t1 = title(['\bfPC' num2str(last_eof)], 'fontsize', 14);

%  And, plot the EOF's.  Remember that the EOF's are 2-dimensional:  row
%  dimension by color dimension.  So, plot three lines that correspond to
%  the three colors:
global_axes(6.4*.6, 1.25, .5, .5, 2*.6*4.8+0.5+0.5+1.5);
subplot2(1,1); cla;
tem = squeeze(reshape(lds(:,last_eof), nlat, nlon));
pp = plot(1:nlat, tem);
set(pp(1), 'color', 'r');
set(pp(2), 'color', 'g');
set(pp(3), 'color', 'b');

%  If we're plotting the first EOF, make it look good (remember, we
%  haven't removed the mean)
if last_eof == 1; set(gca, 'YLim', [0 0.03]); end;

%  Make it pretty:
set(pp, 'linewidth', 2);
set(gca, 'XLim', [1 nlat]);
grid on
xl = xlabel(['\bfEOF' num2str(last_eof) ...
     ', ' num2str(round(10000*per(last_eof))/100) '% Var Expl.'], ...
     'fontsize', 14);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Note:  I've included all of the plotting lines above in a .m script,
%  so we don't have to cut and paste over and over again.  Instead, just
%  type the following line in MATLAB (after downloading the appropriate
%  .m file into an appropriate directory):

cd ~/Class/aos575/matlab
last_eof = 2; plot_island_demo;