clear
loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/.monthly/a:/.sst/:a:/.mask/1/masklt/add/:a/T/%28Jan%201982%29%28Dec%202009%29RANGE/yearly-anomalies/Y/-30/30/RANGE/X/%28110E%29180/%2860W%29RANGE/dods')
asum=asum.asum;
x=asum;mask = mean(asum,3)*0.0;
x(isnan(x))=0;

loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/.monthly/a:/.sst/:a:/.mask/1/masklt/add/:a/T/%28Jul%201982%29%28Jun%202010%29RANGE/yearly-anomalies/Y/-30/30/RANGE/X/%28110E%29180/%2860W%29RANGE/dods')

asum=asum.asum;
y=asum;
[ny,nx,nt]=size(y);
msk = isnan(sum(y,3));
msk = reshape(msk,ny*nx,1);
y(isnan(y))=0;
[uy,sy,vy]=svd(reshape(y,ny*nx,nt),'econ');
[ux,sx,vx]=svd(reshape(x,ny*nx,nt),'econ');
yvar = trace(sy*sy')/nt;
xvar = trace(sx*sx')/nt;


%pick eof truncation

mmax = 75;
yeof = sy*vy';
clear er
for m = 1:mmax;
  im=1:m;
  yhat = zeros(nt,nt);
  for ii=1:nt
    %ind = setdiff([1:nt],[(ii-6):(ii+6)]);
    %ind = setdiff([1:nt],[(ii-1):(ii+1)]);
    ind = setdiff([1:nt],[ii]);%<-PITFALL
    xx = sx(im,im)*vx(ind,im)';
    yy = sy(im,im)*vy(ind,im)';
    %yy = sy*vy(ind,:)';
    A = yy*xx'*inv(xx*xx');
    yhat(1:m,ii) = A*sx(im,im)*vx(ii,im)';
    %yhat(:,ii) = A*sx(im,im)*vx(ii,im)';
  end
  er(m) = mean((yeof(:)-yhat(:)).^2);
  %c(:,m) = correlation(uy*yeof,uy*yhat);
  %find(er==min(er));
end

im = 1:find(er==min(er));
yhat = zeros(nt,nt);
for ii=1:nt
  %ind = setdiff([1:nt],[(ii-6):(ii+6)]);
  %ind = setdiff([1:nt],[(ii-1):(ii+1)]);
  ind = setdiff([1:nt],[ii]);%<-PITFALL
  xx = sx(im,im)*vx(ind,im)';
  yy = sy(im,im)*vy(ind,im)';
  %yy = sy*vy(ind,:)';
  A = yy*xx'*inv(xx*xx');
  yhat(im,ii) = A*sx(im,im)*vx(ii,im)';
  %yhat(:,ii) = A*sx(im,im)*vx(ii,im)';
end
er(m) = mean((yeof(:)-yhat(:)).^2);
c = correlation(uy*yeof,uy*yhat);



my_pcolor(X,Y,reshape(c,ny,nx))
cbar('vert',[-1:0.2:1])
title(['mean correlation = ',num2str(nanmean(c(:)),2)])
pacific
filename='sst_cor_out1';
export_fig('-pdf',filename)

mmax = 75;
yeof = sy*vy';
clear er
for m = 1:mmax;
  im=1:m;
  yhat = zeros(nt,nt);
  for ii=1:nt
    ind = setdiff([1:nt],[(ii-6):(ii+6)]);
    %ind = setdiff([1:nt],[(ii-1):(ii+1)]);
    %ind = setdiff([1:nt],[ii]);%<-PITFALL
    xx = sx(im,im)*vx(ind,im)';
    yy = sy(im,im)*vy(ind,im)';
    %yy = sy*vy(ind,:)';
    A = yy*xx'*inv(xx*xx');
    yhat(1:m,ii) = A*sx(im,im)*vx(ii,im)';
    %yhat(:,ii) = A*sx(im,im)*vx(ii,im)';
  end
  er(m) = mean((yeof(:)-yhat(:)).^2);
  %c(:,m) = correlation(uy*yeof,uy*yhat);
  %find(er==min(er));
end

im = 1:find(er==min(er));
yhat = zeros(nt,nt);
for ii=1:nt
  ind = setdiff([1:nt],[(ii-6):(ii+6)]);
  %ind = setdiff([1:nt],[(ii-1):(ii+1)]);
  %ind = setdiff([1:nt],[ii]);%<-PITFALL
  xx = sx(im,im)*vx(ind,im)';
  yy = sy(im,im)*vy(ind,im)';
  %yy = sy*vy(ind,:)';
  A = yy*xx'*inv(xx*xx');
  yhat(im,ii) = A*sx(im,im)*vx(ii,im)';
  %yhat(:,ii) = A*sx(im,im)*vx(ii,im)';
end
er(m) = mean((yeof(:)-yhat(:)).^2)
c = correlation(uy*yeof,uy*yhat);

my_pcolor(X,Y,reshape(c,ny,nx))
cbar('vert',[-1:0.2:1])
title(['mean correlation = ',num2str(nanmean(c(:)),2)])
pacific
filename='sst_cor_out13';
export_fig('-pdf',filename)

return
%cca
xx = vx(:,im)';
yy = vy(:,im)';
A = yy*xx'*inv(xx*xx');
[u,s,v]=svd(A);

py = uy(:,im)*sy(im,im)*u ;
py(msk,:)=NaN;

px = ux(:,im)*sx(im,im)*v ;
px(msk,:)=NaN;

ty = u'*vy(:,im)';
tx = v'*vx(:,im)';

for ii=im
  subplot(3,1,1)
  my_pcolor(X,Y,reshape(px(:,ii),ny,nx)/sqrt(nt))
  cbar('vert',[-1:0.2:1]/2)
  subplot(3,1,2)
  my_pcolor(X,Y,s(ii,ii)*reshape(py(:,ii),ny,nx)/sqrt(nt))
  cbar('vert',[-1:0.2:1]/2)
  subplot(3,1,3)
  plot(T/12+1960,sqrt(nt)*[tx(ii,:);ty(ii,:)])
  title(corr2(ty(ii,:),tx(ii,:)))
  pause
end


%
xx = sx(im,im)*vx(:,im)';
yy = vy(:,im)';
A = yy*xx'*inv(xx*xx');
svd(A)

%
xx = vx(:,im)';
yy = sy(im,im)*vy(:,im)';
A = yy*xx'*inv(xx*xx');
svd(A)




