% SVDt.m
% Singular Value Decomposition to determine the eigenfunctions
% and amplitudes from satellite imag time series in which the number of 
% spatial points is large.  Does not require calculation of covariance matrix.
% Adapted from code written by % T. Mavor 12 November 1997
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%		eigenfunction	amplitude
%		-------------   ---------
%  TEMPORAL	     U		    SV		<---------
%  SPATIAL	     V		    US
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For temporal case, eigenfunctions are mapped and amplitudes
% show the time dependence

clear

% julian day for image time series
% jd = [50 70 85 95 105 115 125 135 145 155 165 175];	
% jd = [60 90 110 130 150 170]; 
jd = [90 110 130 150 170]; 

%%%%%%%%%%%%%%%%%%%%%%%
% Load data
%%%%%%%%%%%%%%%%%%%%%%%
load feofinput.mat 		% temporal mean not removed; in gray values
			% saved eofinput ssidx meanregchl

%%% CHANGE to log(CHL) 
[nr,nc]=size(eofinput);
for C = 1:nc
 eofinput(:,C) = log(gray2chl(eofinput(:,C))) - log(gray2chl(meanregchl)); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SVD --> eigenvalues and eigenfunctions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[U,S,V]=svd(eofinput,0);	% calculate the Single Value Decomposition
				% for temporal SVD 
				% U contains the eigenfunctions
EGN=diag(S*S');			% the eigenvalues
energy=EGN./sum(EGN);		% and the energy associated with each mode

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine number of modes needed to capture threshold % variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thresh = 0.90;
idx = find(cumsum(energy) > thresh); nmode = min(idx);
disp([num2str(nmode), ' modes are needed to capture ',num2str(100*thresh),'% of the variance'])

%%%%%%%%%%%
% AMPLITUDE 
%%%%%%%%%%%
for mode = 1:nmode
 A(:,mode)=V(:,mode)*sqrt(EGN(mode));		% Amplitude 
end

U=U*(-1); A=A*(-1);		% For more intuitive representation

%%%%%%%%%%%%%%%%%%%%%%%
% RECONSTRUCT INPUT
%%%%%%%%%%%%%%%%%%%%%%%
%% Reconstruct the data from first N modes
%% Eigenfunctions are in U
[nr nc] = size(eofinput);
disp(['Calculating time series from first ' int2str(nmode) ' modes...'])
for T = 1:nc
 base = zeros(nr,1);
 %% For log-transformed, first reconstruct log(chl)-log(meanchl), i.e. input;
 for mode = 1:nmode
%  base = base + exp(U(:,mode)*A(T,mode));
   base = base + U(:,mode)*A(T,mode);
 end

%% Add back log(meanregchl -- which is in gray values)
 base = base + log(gray2chl(meanregchl));

%% Convert back to chl
 base = exp(base);

 % Convert back to gray
 recon(:,T) = chl2gray(base);
% recon(:,T) = chl2gray(base) + meanregchl;
end

save reconTS80 recon ssidx 
		
%%%%%%%%%%%%%%%%%%%%%%%
% DISPLAY 
%%%%%%%%%%%%%%%%%%%%%%%
%   Load Ancillary data 
load mabtopo
load neall
load NECtopo
load mlab_chlpal

save SVD

meanCZCS = ones(512,512) * NaN; meanCZCS(ssidx) = meanregchl;
eofaxis = ([min(nealon) -66 37 42]);
numc = 128;

% Plot the temporal mean (ensemble average)
figure(1); clf; whitebg white;
axes('position',[0 0 1 1]);
image(nealon,nealat,meanCZCS); 
colormap(chlpal); hold on
p = eofaxis; land_mask;
plot(iso100(:,1),iso100(:,2),'k--');
set(gca,'XTick',[],'YTick',[]);

for mode = 3:3
 eval(['figure(',int2str(mode+1),'); clf'])
 whitebg white
 cmode = U(:,mode);
 datlim = [min(cmode) max(cmode)]; datrange = datlim(2)-datlim(1);
 cax = [-.012 .012];
 cim = ones(512,512) * NaN; cim(ssidx) = cmode;
 axes('position',[.10 .37 .75 .54]);
 pcolor(nealon,nealat,flipud(cim)); shading flat;
 colormap(jet(numc)); hold on
% caxis(datalim);
 caxis(cax);
% image(nealon,nealat,cim);
% caxis(cax); colormap(jet(numc)); hold on
 contour(nealon,nealat,flipud(cim),[0 0],'k');
 p = eofaxis; land_mask; 
 iso=plot(iso100(:,1),iso100(:,2),'w--'); set(iso,'color',[.5 .5 .5],'linewidth',2);
% extcontour(nealon,nealat,flipud(NECtopo),[-50 -50],'label',8000,'k--');
 set(gca,'TickDir','out','fontweight','bold','fontsize',15)
 set(gca,'YTick',[37:1:42],'XTick',[-76:2:-66]);
 set(gca,'YTicklabels',['37 N';'    ',;'    ';'    ';'    ';'42 N']);
 set(gca,'XTicklabels',['76 W';'    ';'    ',;'    ';'    ';'66 W']);
% t = text(-69.5,42,'100m'); 
% set(t,'fontweight','bold','rotation',-45)
 axis(eofaxis);

 %% Plot colorbar
 axes('position',[.94 .38 .03 .52])
%  cbrange=[datlim(1):datrange/numc:datlim(2)]';
cbrange=[cax(1):datrange/numc:cax(2)]';
 clear scale; scale(:,1)=cbrange;scale(:,2)=cbrange;
 pcolor([0:1],cbrange,scale)
 colormap(jet(numc)); shading flat;caxis(cax);	%caxis(datlim);
 set(gca,'XTick',[],'fontweight','bold','fontsize',15,'box','off','YTick',[-.01 0 .01])

 %% Plot coefficient
 axes('position',[.1 .12 .8 .15]);
 set(gca,'fontweight','bold','fontsize',15);
 plot(jd,A(:,mode),'b*'); hold on; pf=plot(jd,A(:,mode),'b'); set(pf,'linewidth',2);
 set(gca,'fontweight','bold','fontsize',15);
 axis([min(jd)-1 max(jd)+1 min(A(:,mode))-.05*datrange max(A(:,mode))+.05*datrange]);
 xl = get(gca,'xlim');
 plot(xl,[0 0],'k--');

 %% Plot the energy
 % semilogy(1:size(energy),energy,'-',1:size(energy),energy,'go')
 % axis([0 100 10^(-5) 1])

 axes('position',[0 0 1 1]); set(gca,'visible','off');

 %% EOF Plot label
 str1 = ['MODE ',int2str(mode),' : (',num2str(energy(mode)*100), ' % of the variance)'];
 putit = ['t = text(.5,.95,''',str1,' ''',');']; eval(putit)
 set(t,'fontweight','bold','fontsize',15,'horizontalalignment','center')
 tstr=['stext(.5,.9,' '''\bold\15 ' tstr ',''', 'Horizontalalignment''',',','''Center''',')'];
stext(.5,.04,'\bold\15 Julian day','horizontalalignment','center');
 %%  Time axis label
% t = text(.21,.07,'March            April             May            June')
% set(t,'fontweight','bold')

 %% Coefficient Plot label
% str2 = ['MODE ',int2str(mode),' coefficient']; 
% putit = ['t = text(.3,.27,''',str2,' ''',')']; eval(putit)
% set(t,'fontweight','bold','fontsize',15)

keyboard
end

