clear

%use predictable components or EOFS or something from CFS

%loaddap('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/.MONTHLY/.t2m/Y/-20/20/RANGE/S/%281%20Oct%29VALUES/S/%281%20Oct%201960%29%281%20Oct%202001%29RANGE/L/3/runningAverage/L/3.5/VALUE/dods')

loaddap('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/.MONTHLY/.prec/Y/-20/20/RANGE/S/%281%20Oct%29VALUES/S/%281%20Oct%201960%29%281%20Oct%202001%29RANGE/L/3/runningAverage/L/3.5/VALUE/86400000/mul/dods')

model = squeeze(prec.prec);
%model = squeeze(t2m.t2m);

[p,q,bst,st,nn,nm,np] = sn_eofs_function_2(model,Y);
for im=1:np
  subplot(2,1,1)
  my_pcolor(X,Y,p(:,:,im))
  cbar('vert',[-1.5:0.5:1.5]/2)
  pacific_noresize_at_all
  subplot(2,1,2)
  bar((S+L)/12+1960,bst(:,im))
  pause
end

%eofs

%loaddap('http://iridl.ldeo.columbia.edu/expert/IRIONLY/ds:/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/.MONTHLY/.t2m/Y/-15/15/RANGE/S/%281%20Oct%29VALUES/S/%281%20Oct%201960%29%281%20Oct%202001%29RANGE%5BM%5Daverage/L/3/runningAverage/L/3.5/VALUE/home/.tippett/.scratch/.na/.hindcasts/.obs/.t2m/S/%281%20Oct%29VALUES/S/%281%20Oct%201960%29%281%20Oct%202001%29RANGE/L/3.5/VALUE/%28name%29%28obs%29def/:ds/dods')

loaddap('http://iridl.ldeo.columbia.edu/expert/IRIONLY/ds:/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/.MONTHLY/.t2m/Y/-15/15/RANGE/S/%281%20Oct%29VALUES/S/%281%20Oct%201960%29%281%20Oct%202001%29RANGE%5BM%5Daverage/L/3/runningAverage/L/3.5/VALUE/home/.tippett/.scratch/.na/.hindcasts/.obs/.prcp/S/%281%20Oct%29VALUES/S/%281%20Oct%201960%29%281%20Oct%202001%29RANGE/L/3.5/VALUE/%28name%29%28obs%29def/:ds/dods')


obs = squeeze(obs.obs);
[my,mx,nt] = size(obs);
y = reshape(obs,my*mx,nt);


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

t2m_mean = mean(t2m,3);
t2m_std = std(t2m,[],3);
t2ma = (t2m - t2m_mean(:,:,ones(1,nt)))./t2m_std(:,:,ones(1,nt));;

[u,s,v] = svd(reshape(t2ma,ny*nx,nt),'econ');




nmodes=6

for im=1:nmodes
  subplot(2,1,1)
  my_pcolor(X_01,Y_01,t2m_std.*s(im,im).*reshape(u(:,im),ny,nx))
  cbar
  pacific_noresize_at_all
  subplot(2,1,2)
  bar((S+L)/12+1960,v(:,im))
  pause
end
clf


%x = [ones(1,nt);v(:,1:nmodes)'];
x = [ones(1,nt);bst'];




m=size(x,1);

a = y*x'*inv(x*x');
for n=1:m
  ind = 1:n;
  sse(:,n) = sum((y-a(:,ind)*x(ind,:)).^2,2);
end

wrk = repmat([1:m],[my*mx,1]);

%mallow's cp
cp = sse./sse(:,m*ones(1,m))*nt - nt + 2*wrk;
[dum,imin] = min(cp,[],2);

cp_trunc = reshape(imin,my,mx);
cp_trunc(isnan(dum))=NaN;

C = colormap;


my_pcolor(X,Y,cp_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('Mallow''s C_p')
pause

%aic
aic = nt*log(sse) + 2*wrk;
[dum,imin] = min(aic,[],2);
aic_trunc = reshape(imin,my,mx);
aic_trunc(isnan(dum))=NaN;

my_pcolor(X,Y,aic_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('AIC')
pause


%corrected aic
aicc = nt*log(sse) + 2*wrk + 2*wrk.*(wrk+1)./(nt-wrk-1);
[dum,imin] = min(aicc,[],2);
aicc_trunc = reshape(imin,my,mx);
aicc_trunc(isnan(dum))=NaN;

my_pcolor(X,Y,aicc_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('corrected AIC')
pause



%bic
bic = nt*log(sse) + log(nt)*wrk;
[dum,imin] = min(bic,[],2);
bic_trunc = reshape(imin,my,mx);
bic_trunc(isnan(dum))=NaN;

my_pcolor(X,Y,bic_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('BIC')
pause


%f-test
f1 = (sse(:,ones(1,m))-sse)./(wrk-1);
f2 = sse./(nt-wrk);

f = f1./f2;
f(:,1)=NaN;
p = 1 - fcdf(f,wrk-1,nt-wrk);

[dum,imin] = min(p,[],2);
imin(dum>0.05)=1;
f_trunc = reshape(imin,my,mx);
f_trunc(isnan(dum))=NaN;

my_pcolor(X,Y,f_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('F-test')
pause


%cross-validation
for n=1:m
  ind=1:n;
  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
  sse_cv(:,n) = sum(res_cv.^2,2);
end

[dum,imin] = min(sse_cv,[],2);
cv_trunc = reshape(imin,my,mx);
cv_trunc(isnan(dum))=NaN;

my_pcolor(X,Y,cv_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('Cross-validation')
pause


%split data

indt = 1:(nt/2);
a = y(:,indt)*x(:,indt)'*inv(x(:,indt)*x(:,indt)');
for n=1:m
  ind = 1:n;
  indt = (nt/2+1):nt;
  sse2(:,n) = sum((y(:,indt)-a(:,ind)*x(ind,indt)).^2,2);
end


[dum,imin] = min(sse2,[],2);
split_trunc = reshape(imin,my,mx);
split_trunc(isnan(dum))=NaN;


my_pcolor(X,Y,split_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('Split 1')
pause


indt =(nt/2+1):nt;
a = y(:,indt)*x(:,indt)'*inv(x(:,indt)*x(:,indt)');
for n=1:m
  ind = 1:n;
  indt = 1:(nt/2);
  sse2(:,n) = sum((y(:,indt)-a(:,ind)*x(ind,indt)).^2,2);
end

[dum,imin] = min(sse2,[],2);
split_trunc = reshape(imin,my,mx);
split_trunc(isnan(dum))=NaN;


my_pcolor(X,Y,split_trunc')   
cbar('vert',[1:m])
pacific_noresize_at_all
title('Split 2')
pause


iok = ~isnan(y(:,1));

[mean(aic_trunc(iok)<bic_trunc(iok)), ...
mean(aic_trunc(iok)>bic_trunc(iok)), ...
mean(aic_trunc(iok)==bic_trunc(iok))]

[mean(aic_trunc(iok)<cv_trunc(iok)), ...
mean(aic_trunc(iok)>cv_trunc(iok)), ...
mean(aic_trunc(iok)==cv_trunc(iok))]

[mean(aic_trunc(iok)<cp_trunc(iok)),...
mean(aic_trunc(iok)>cp_trunc(iok)),...
mean(aic_trunc(iok)==cp_trunc(iok))]

[mean(aic_trunc(iok)<aicc_trunc(iok)),...
mean(aic_trunc(iok)>aicc_trunc(iok)),...
mean(aic_trunc(iok)==aicc_trunc(iok))]