% Seasonal cycle analysis of CMIP3 models, reanalysis data and IMD rainfall
% CCA following Tippett(2003)
% April 2009

clear
%% IMD data
load /Users/andy/data/India/IMD/IMD-Daily1DegRain-1901-2004/IMD-Daily1DegRain-1901-2004.mat
nx=35; ny=33; ndays=dayN-day1+1;
A=reshape(A,ndays,ny*nx);
x_axis=X(1,:)'; y_axis=Y(:,1);
X=reshape(X,1,ny*nx); Y=reshape(Y,1,ny*nx);
dvec=datevec(day1:dayN);

% select peninsular India, land points
% xysel=find(Y<=30 & X<=88 & A(end,:)~=-999);
xysel=find(A(end,:)~=-999);
A1=A(:,xysel);
X1=Y(xysel); Y1=Y(xysel);
nxy=length(xysel);

% Monthly averages,averaged over years 1971-2000
nmon=12;
rain=NaN*ones(nmon,nxy);
imon=0; 
for mon=1:12 
  imon=imon+1;
  dsel=find(dvec(:,2)==mon & dvec(:,1)>=1971 & dvec(:,1)<=2000);
  rain(imon,:)=mean(A1(dsel,:));
end
clear A A1

%% GCM 
% ECHAM5-MPI model
ncload('/Users/andy/data/CMIP3/MPI/data/monthly/mpi_echam5.20c3m.run1.mo.hus.1970-2000.nc')
[nmo nlev nlat nlon]=size(hus);
hus=reshape(hus,12,30,nlat,nlon);
husm=squeeze(mean(hus,2));

% select lat-lon domain, and calendar months
lat1=1; lat2=34;
lon1=17; lon2=97; % (30-180E)
husms=husm(:,lat1:lat2,lon1:lon2);
lat=lat(lat1:lat2); lon=lon(lon1:lon2);
[nmon nlat nlon]=size(husms);
hus=reshape(husms,nmon,nlat*nlon);
 
%% predictor (X) and predictand (Y) matrix: Y=A*X
X=(hus-repmat(mean(hus),nmon,1))'; % remove time mean and tranpose
Y=(rain-repmat(mean(rain),nmon,1))'; % remove time mean and tranpose

%% EOFs: here U-EOFs; V-PCs
[Ux,Sx,Vx]=svd(X,'econ'); 
VARx=100*diag(Sx).^2/sum(diag(Sx).^2);

[Uy,Sy,Vy]=svd(Y,'econ'); 
VARy=100*diag(Sy).^2/sum(diag(Sy).^2);
disp(['PC variances of Y and X: '])
disp([(1:12)' VARy VARx])

%% truncate PCs
nY=4; nX=3;
Ux=Ux(:,1:nX); Vx=Vx(:,1:nX); Sx=Sx(1:nX,1:nX);
Uy=Uy(:,1:nY); Vy=Vy(:,1:nY); Sy=Sy(1:nY,1:nY);

%% correlation matrix of PCs: B=Vy'*Vx
B=Vy'*Vx;
% SVD: B=RMS'
[R,M,S]=svd(B,'econ');
nCCA=length(M);
disp(['CCA Time Series Correlations:'])
disp(diag(M))

%% CCA observation and predictor homogeneous covariance maps
Cy=Uy*Sy*R;
Cx=Ux*Sx*S;

%% CCA time series
Ty=Vy*R;
Tx=Vx*S;

%% plot CCA time series & maps
nCCA=3;
figure(1); clf
plot(Tx)
legend('CCA 1','CCA 2','CCA 3')
hold on
plot(Ty,'--')

%%
figure(2); clf
Cx=reshape(Cx,nlat,nlon,nCCA);
count=0;
for k=1:nCCA
  count=count+1;
  subplot(3,2,count) % IMD rainfall map
  map=NaN*ones(ny,nx);
  map(xysel)=Cy(:,k);
  plot_map2(x_axis,y_axis,map,[-20 20])
  title(['CCA ',num2str(k),': Rainfall'])
  count=count+1;
  subplot(3,2,count) % predictor  map
  sf=1;
  map=Cx(:,:,k)*sf;
  max(max(map))
  plot_map2(lon,lat,map,[min(min(map)) max(max(map))])
  title(['CCA ',num2str(k),': SpHum'])
end

%% prediction: Y=A*X with regression matrix A=Cy*M*Cx^(-1)



