clear
close all
figure
set(gcf,'color',[1,1,1])
%read obs
yearstr = '%281971%29%282007%29';
%yearstr = '%281971%29%282005%29';

loaddods(['http://iridl.ldeo.columbia.edu/expert/SOURCES/.WCRP/.GCOS/.GPCC/.FDP/.version4/.0p5/.prcp/X/%28115E%29%28130E%29RANGE/Y/%285N%29%2820N%29RANGE/T/3/runningAverage/T/%28Apr-Jun%29VALUES/T/',yearstr,'RANGE/dup/SOURCES/.Features/.Political/.World/.Countries/.the_geom/objectid/185/VALUE/rasterize/0/flaggt/0/maskle/mul/12/mul/365.25/div/(name)/(prcp)/def/dods'])
obs = prcp.prcp;
[my,mx,nt] = size(obs);
Xo=X;
Yo=Y;

loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/T/3/runningAverage/T/%28Apr-Jun%29VALUES/T/%281971%29%282007%29RANGE/dods')
nino_sim=NINO34.NINO34;

c_nino_sim = correlation2D1D(obs,nino_sim);
c_nino_sim(abs(c_nino_sim)<2/sqrt(nt))=NaN;
my_pcolor(X,Y,c_nino_sim)
cbar('vert',[-1:0.2:1])
pacific

loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/T/3/runningAverage/T/%28Jan-Mar%29VALUES/T/%281971%29%282007%29RANGE/dods')
nino_1=NINO34.NINO34;

c_nino_1 = correlation2D1D(obs,nino_1);
c_nino_1(abs(c_nino_1)<2/sqrt(nt))=NaN;
my_pcolor(X,Y,c_nino_1)
cbar('vert',[-1:0.2:1])
pacific

%preceeding SST
loaddods(['http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCDC/.ERSST/.version2/.SST/X/100./220./RANGEEDGES/Y/-25/25/RANGEEDGES/T/3/runningAverage/T/%28Jan-Mar%29VALUES/T/',yearstr,'RANGE/%28name%29%28pred%29def%5B/X/Y%5DREORDER/dods'])

pred=pred.pred;
[ny,nx,nt] = size(pred);


my_pcolor(X,Y,pred(:,:,1)-mean(pred,3))
cbar('vert',[-1.6:0.4:1.6])
title('SST anomaly JFM 1971')
pacific
filename='sst_jfm_1971';
export_fig('-pdf',filename)

my_pcolor(Xo,Yo,obs(:,:,1)-mean(obs,3))
cbar('vert',[-1.6:0.4:1.6]*5)
title('Rainfall anomalies AMJ 1971')
pacific
green_brown_cm2
filename='prcp_amj_1971';
export_fig('-pdf',filename)



x = reshape(pred,ny*nx,nt);
xm = mean(x,2);
msk = ~isnan(xm);
xa = x - xm(:,ones(1,nt)); 
xa(~msk,:)=0;

[u,s,v] = svd(xa,'econ');
u(~msk,:)=NaN;
colormap('default')
fac=[1,.5,0.5];
varexpl=diag(s.^2)./trace(s.^2);

for ii=1:3
  subplot(4,1,1:3)
  my_pcolor(X,Y,reshape(s(ii,ii)*u(:,ii),ny,nx)/sqrt(nt))
  cbar('vert',[-1:0.2:1]*fac(ii))
  pacific
  title(['EOF ',int2str(ii),' variance explained =',num2str(varexpl(ii)*100,2),'%'])
  subplot(4,1,4)
  bar(T/12+1960,v(:,ii)*sqrt(nt))
  grid
  axis('tight')
  pause
  filename=['sst_jfm_eof',int2str(ii)];
  export_fig('-pdf',filename)

end
clf

for ii=1:5
  c = correlation2D1D(obs,v(:,ii)');
  c(abs(c)<2/sqrt(nt))=NaN;
  my_pcolor(Xo,Yo,c)
  cbar('vert',[-1:0.2:1])
  pacific
  pause
end


%pcr
x = [ones(nt,1),v(:,1:3)]';
y = reshape(obs,my*mx,nt);

inds={[1],[1,2]};
m=length(inds);
sumse = zeros(my*mx,m);
nvars=zeros(1,m);
for n=1:m
  ind = inds{n};
  a = y*x(ind,:)'*inv(x(ind,:)*x(ind,:)');
  sumse(:,n) = sum((y-a*x(ind,:)).^2,2);
  nvars(n) = nnz(ind);
end
wrk = repmat(nvars,[my*mx,1]);
%aic
aic = nt*log(sumse) + 2*wrk + 2*wrk.*(wrk+1)./(nt-wrk-1);
[dum,imin] = min(aic,[],2);
aic_trunc = reshape(imin,my,mx);
aic_trunc(isnan(dum))=NaN;

colormap(jet(m))
my_pcolor(Xo,Yo,aic_trunc)
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',{'CL','ENSO','2','3','4'})
pacific
title('AIC selected model')
pause
filename='aic_trunc_ph';
export_fig('-pdf',filename)

%model selection is not included in the cross-validation
colormap('default')
yhat = zeros(size(y))+NaN;
for ii=1:m
  ind = inds{ii};
  indx = reshape(aic_trunc==ii,my*mx,1);
  %ind=[1,2];
  %indx=1:(my*mx);
  V = x(ind,:)'*inv(x(ind,:)*x(ind,:)')*x(ind,:);%HAT matrix
  res = y*(eye(nt)-V);
  wrk2 = ones(my*mx,1)*(1./(1-diag(V)))';
  res_cv = res.*wrk2;%cv residual w/o cv
  yhat(indx,:) = y(indx,:) - res_cv(indx,:);
end

clf
c = reshape(correlation(y,yhat),my,mx);
cpw = c;
c(abs(c)<2/sqrt(nt))=NaN;
my_pcolor(Xo,Yo,c)
pacific
cbar('vert',[-1:0.2:1])
title('correlation')
filename='corr_ph';
export_fig('-pdf',filename)

pause

ym = mean(y,2);
sterror_pw = 1-reshape(mean((y-yhat).^2,2)./mean((y-ym(:,ones(1,nt))).^2,2),my,mx);
my_pcolor(Xo,Yo,sterror_pw)
pacific
cbar('vert',[0:0.2:1]*0.5)
title('1-error var/clim var')
filename='explvar_ph';
export_fig('-pdf',filename)

%model selection IS included in the cross-validation
sumse = zeros(my*mx,m);

colormap('default')
yhat = zeros(size(y))+NaN;
for it=1:nt
  indt = setdiff([1:nt],it);
  for n=1:m
    ind = inds{n};
    a = y(:,indt)*x(ind,indt)'*inv(x(ind,indt)*x(ind,indt)');
    sumse(:,n) = sum((y(:,indt)-a*x(ind,indt)).^2,2);
    nvars(n) = nnz(ind);
  end
  wrk = repmat(nvars,[my*mx,1]);
  %aic
  aic = (nt-1)*log(sumse) + 2*wrk + 2*wrk.*(wrk+1)./(nt-wrk-1);
  [dum,imin] = min(aic,[],2);
  aic_trunc = reshape(imin,my,mx);
  aic_trunc(isnan(dum))=NaN;
  AIC_trunc(:,:,it) = aic_trunc;
  for ii=1:m
    ind = inds{ii};
    indx = reshape(aic_trunc==ii,my*mx,1);
    %ind=[1,2];
    %indx=1:(my*mx);
    V = x(ind,:)'*inv(x(ind,:)*x(ind,:)')*x(ind,:);%HAT matrix
    res = y*(eye(nt)-V);
    wrk2 = ones(my*mx,1)*(1./(1-diag(V)))';
    res_cv = res.*wrk2;%cv residual w/o cv
    yhat(indx,it) = y(indx,it) - res_cv(indx,it);
  end
end


colormap(jet(m))
my_pcolor(Xo,Yo,max(AIC_trunc,[],3)) 
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',{'CL','ENSO','2','3','4'})
pacific
title('max AIC selected model')
pause
filename='aic_trunc_max_ph';
export_fig('-pdf',filename)

my_pcolor(Xo,Yo,min(AIC_trunc,[],3)) 
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',{'CL','ENSO','2','3','4'})
pacific
title('min AIC selected model')
pause
filename='aic_trunc_min_ph';
export_fig('-pdf',filename)

colormap('default')

clf
c = reshape(correlation(y,yhat),my,mx);
cpw2 = c;
c(abs(c)<2/sqrt(nt))=NaN;



plot(cpw(:),cpw2(:),'.')
hold on
plot([-1,1],[-1,1],'k')
hold off
grid
axis('equal')
set(gca,'xlim',[-1,1],'ylim',[-1,1])
xlabel('w/o cv pred. sel.')
ylabel('w/ cv pred. sel.')
title('correlation')
filename='delta_corr_scatter_ph';
export_fig('-pdf',filename)


ym = mean(y,2);
sterror_pw2 = 1-reshape(mean((y-yhat).^2,2)./mean((y-ym(:,ones(1,nt))).^2,2),my,mx);

my_pcolor(Xo,Yo,sterror_pw2-sterror_pw)
pacific
cbar('vert',[-1:0.2:1]*0.2)
title('Change in 1-error var/clim var')
filename='delta_explvar_ph';
export_fig('-pdf',filename)

my_pcolor(Xo,Yo,cpw2-cpw)
pacific
cbar('vert',[-1:0.2:1]/2)
title('Change in correlation')
filename='delta_corr_ph';
export_fig('-pdf',filename)



plot(sterror_pw(:),sterror_pw2(:),'.')
hold on
plot([-1.3,1],[-1.3,1],'k')
hold off
grid
axis('equal')
set(gca,'xlim',[-.3,0.5],'ylim',[-.3,0.5])
xlabel('w/o cv pred. sel.')
ylabel('w/ cv pred. sel.')
title('1-error var/clim var')
filename='delta_explvar_scatter_ph';
export_fig('-pdf',filename)



%-----------------------------------------------------------------------------
%check if pcr is better than pattter/pattern
% 2 obs eof + 1 pred eof

y = reshape(obs,my*mx,nt);
msk = ~isnan(y(:,1));
y = y(msk,:);
ym = mean(y,2);
[uy,sy,vy]=svd(y-ym(:,ones(1,nt)),'econ');
[ux,sx,vx]=svd(xa,'econ');

vy(:,1:2)'*vx(:,1)*inv(vx(:,1)'*vx(:,1))

corr2(vy(:,1),vx(:,1))

%pick eof truncation


inds={[1],[1,2]};
m=length(inds);
sumse = zeros(nt,m);
nvars=zeros(1,m);
for n=1:m
  ind = inds{n};
  a = vy'*x(ind,:)'*inv(x(ind,:)*x(ind,:)');
  sumse(:,n) = sum((vy'-a*x(ind,:)).^2,2);
  nvars(n) = nnz(ind);
end
wrk = repmat(nvars,[nt,1]);
%aic
aic = nt*log(sumse) + 2*wrk + 2*wrk.*(wrk+1)./(nt-wrk-1);
[dum,imin] = min(aic,[],2);
aic_trunc = imin;

subplot(2,1,1)
bar(aic_trunc(1:(nt-1)))
set(gca,'ytick',[1,2],'yticklabel',{'CL','ENSO'},'xtick',[5:5:35],'xlim',[0,37],'ylim',[0,2.2])
xlabel('PCs of y')
subplot(2,1,2)
bar(aic(1:(nt-1),1)-aic(1:(nt-1),2))
set(gca,'xtick',[5:5:35],'xlim',[0,37])
xlabel('PCs of y')
axis('tight')
grid
ylabel('\Delta AIC')

filename='pca_y_aic_trunc_ph';
export_fig('-pdf',filename)

clf

im=1;
in=1:2;
clear yhat yhat2
for ii=1:nt
  ind = setdiff([1:nt],[ii]);%<-PITFALL
  xx = sx(im,im)*vx(ind,im)';

  ym = mean(y,2);
  ya = y - ym(:,ones(1,nt));
  [uy,sy,vy] = svd(ya,'econ');
  yy = sy(in,in)*vy(ind,in)';
  A = yy*xx'*inv(xx*xx');
  yhat(:,ii) = ym+uy(:,in)*A*[sx(im,im)*vx(ii,im)'];
  
  ym = mean(y(:,ind),2);
  ya = y - ym(:,ones(1,nt));
  [uy,sy,vy] = svd(ya(:,ind),'econ');
  yy = sy(in,in)*vy(:,in)';
  A = yy*xx'*inv(xx*xx');
  yhat2(:,ii) = ym+uy(:,in)*A*[sx(im,im)*vx(ii,im)'];

end

wrk = zeros(my*mx,nt)+NaN;
wrk(msk,:)=yhat;
yhat=wrk;

wrk = zeros(my*mx,nt)+NaN;
wrk(msk,:)=yhat2;
yhat2=wrk;

y = reshape(obs,my*mx,nt);
msk = ~isnan(y(:,1));
y = y(msk,:);
ym = mean(y,2);
ya = y - ym(:,ones(1,nt));
[uy,sy,vy] = svd(ya,'econ');
green_brown_cm2
for ii=1:2
  wrk = zeros(size(msk))+NaN;
  wrk(msk,:)=uy(:,ii);
  subplot(4,1,1:3)
  my_pcolor(Xo,Yo,sy(ii,ii)*reshape(wrk,my,mx)/sqrt(nt))
  pacific       
  cbar('vert',[-4:4])
  title(['EOF ',int2str(ii),' of prcp'])
  subplot(4,1,4)
  bar(T/12+1960,vy(:,ii)*sqrt(nt))
  grid
  axis('tight')

  pause
  filename=['prcp_amj_eof',int2str(ii)];
  export_fig('-pdf',filename)

end
uu=[0.8591   -0.5118]'
wrk = zeros(size(msk))+NaN;
wrk(msk,:)=uy(:,1:2)*sy(1:2,1:2)*uu;
subplot(4,1,1:3)
my_pcolor(Xo,Yo,reshape(wrk,my,mx)/sqrt(nt))
pacific       
cbar('vert',[-4:4])
title(['Pattern of prcp'])
subplot(4,1,4)
bar(T/12+1960,vy(:,1:2)*uu*sqrt(nt))
grid
axis('tight')

pause
filename=['prcp_amj_cca'];
export_fig('-pdf',filename)


colormap('default')
clf
y = reshape(obs,my*mx,nt);
c = reshape(correlation(y,yhat),my,mx);
ceof=c;
c(abs(c)<2/sqrt(nt))=NaN;
my_pcolor(Xo,Yo,c)
pacific
cbar('vert',[-1:0.2:1])
title('correlation')
pause

c = reshape(correlation(y,yhat2),my,mx);
ceof2=c;
c(abs(c)<2/sqrt(nt))=NaN;
my_pcolor(Xo,Yo,c)
pacific
cbar('vert',[-1:0.2:1])
title('correlation')
pause
filename='corr_eof2_ph';
export_fig('-pdf',filename)


pause

ym = mean(y,2);
sterror_eof = 1-reshape(mean((y-yhat).^2,2)./mean((y-ym(:,ones(1,nt))).^2,2),my,mx);
sterror_eof2 = 1-reshape(mean((y-yhat2).^2,2)./mean((y-ym(:,ones(1,nt))).^2,2),my,mx);
my_pcolor(Xo,Yo,sterror_eof2)
pacific
cbar('vert',[0:0.2:1]*0.5)
title('1-error var/clim var')
filename='explvar_eof2_ph';
export_fig('-pdf',filename)

clf

plot(cpw(:),ceof(:),'.')
grid
axis('equal')
set(gca,'xlim',[0,1])
set(gca,'ylim',[0,1])
hold on
plot([0,1],[0,1])
hold off
xlabel('pointwise')
ylabel('pattern')
title('correlation')
filename='pw_vs_pat_corr_ph';
%set(gcf,'PaperPosition',[0.25 2.5 4 3])
export_fig('-pdf',filename)

plot(cpw(:),ceof2(:),'.')
grid
axis('equal')
set(gca,'xlim',[0,1])
set(gca,'ylim',[0,1])
hold on
plot([0,1],[0,1])
hold off
xlabel('pointwise')
ylabel('pattern')
title('correlation')
filename='pw_vs_pat_corr_cv_ph';
%set(gcf,'PaperPosition',[0.25 2.5 4 3])
export_fig('-pdf',filename)



plot(sterror_pw(:),sterror_eof(:),'.')
grid
axis('equal')
set(gca,'xlim',[-0.1,0.6])
set(gca,'ylim',[-0.1,0.6])
hold on
plot([0,1],[0,1])
hold off
xlabel('pointwise')
ylabel('pattern')
title('1-error var/cl var')
filename='pw_vs_pat_error_ph';
%set(gcf,'PaperPosition',[0.25 2.5 4 3])
export_fig('-pdf',filename)

plot(sterror_pw(:),sterror_eof2(:),'.')
grid
axis('equal')
set(gca,'xlim',[-0.1,0.6])
set(gca,'ylim',[-0.1,0.6])
hold on
plot([0,1],[0,1])
hold off
xlabel('pointwise')
ylabel('pattern')
title('1-error var/cl var')
filename='pw_vs_pat_error2_ph';
%set(gcf,'PaperPosition',[0.25 2.5 4 3])
export_fig('-pdf',filename)

%%%%%%%%%%%%not quite cv'd
