clear

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/Y/-60/80/RANGE/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);


loaddap('http://iridl.ldeo.columbia.edu/expert/ds:/SOURCES/.Indices/.nino/.EXTENDED/a:/.NINO34/T/3/runningAverage/T/%28Sep-Nov%29VALUES/T/%281960%29%282001%29RANGE/:a:/.NINO4/T/3/runningAverage/T/%28Sep-Nov%29VALUES/T/%281960%29%282001%29RANGE/:a:/.NINO3/T/3/runningAverage/T/%28Sep-Nov%29VALUES/T/%281960%29%282001%29RANGE/:a:/.NINO12/T/3/runningAverage/T/%28Sep-Nov%29VALUES/T/%281960%29%282001%29RANGE/:a/:ds/dods')

n12 = NINO12.NINO12';
n34 = NINO34.NINO34';
n4 = NINO4.NINO4';
n3 = NINO3.NINO3';

ninos = [n12-mean(n12) ; n3-mean(n3); n34-mean(n34) ; n4-mean(n4);];
ninos = [n34-mean(n34) ];

x = [ones(1,nt);ninos;S'];
inds={[1],[1,2],[1,3],[1,2,3]};

m=length(inds);
colormap(jet(m))


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

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

%cross-validation
for n=1:m
  ind = inds{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
  yhat = y - res_cv;
  cor_cv(:,n) = correlation(yhat,y);
  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',[0.5:0.5:(m+0.5)])
pacific_noresize_at_all
title('Cross-validation')
pause

for ii=1:my*mx
  cor_best(ii,1) = cor_cv(ii,imin(ii));
end


my_pcolor(X,Y,1-reshape(sqrt(dum/nt)./std(y,[],2),my,mx).^2')   
cbar
pacific_noresize_at_all
title('Cross-validation rms')
pause


my_pcolor(X,Y,reshape(cor_best,my,mx)')   
cbar
pacific_noresize_at_all
title('Cross-validation correlation')
pause



%double cross-validation
for it =1:nt
  indt = setdiff(1:nt,it);
  for n=1:m
    ind = inds{n};
    V = x(ind,indt)'*inv(x(ind,indt)*x(ind,indt)')*x(ind,indt);%HAT matrix
    res = y(:,indt)*(eye(nt-1)-V);
    wrk2 = ones(my*mx,1)*(1./(1-diag(V)))';
    res_cv = res.*wrk2;%cv residual w/o cv
    sse_cv2(:,n) = sum(res_cv.^2,2);
  end
  [dum,imin] = min(sse_cv2,[],2);
  cv_trunc = reshape(imin,my,mx);
  cv_trunc(isnan(dum))=NaN;
  for n=1:m
    ind = inds{n};
    a = y(:,indt)*x(ind,indt)'*inv(x(ind,indt)*x(ind,indt)');
    yall(:,n) = a*x(ind,it);
  end
  for ii=1:my*mx
    yhat_cv2(ii,it) = yall(ii,imin(ii));
  end
end

error_cv2 = sum((y-yhat_cv2).^2,2);
my_pcolor(X,Y,1-reshape(sqrt(error_cv2/nt)./std(y,[],2),my,mx).^2')   
cbar
pacific_noresize_at_all
title('Cross-validation rms')

cor_cv2 = correlation(y,yhat_cv2);
my_pcolor(X,Y,reshape(cor_cv2,my,mx)')   
cbar
pacific_noresize_at_all
title('Cross-validation correlation')
pause


plot(cor_best(:),cor_cv2(:),'.') 
grid
pause

plot(dum(:)./std(y,[],2).^2/nt,error_cv2(:)./std(y,[],2).^2/nt,'.') 
grid
pause

%cross-validation/aic
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)');
    sse_cva(:,n) = sum((y(:,indt)-a*x(ind,indt)).^2,2);
    nvars(n) = nnz(ind);
  end
  aic = nt*log(sse_cva) + 2*wrk;
  [dum,imin] = min(aic,[],2);
  aic_trunc = reshape(imin,my,mx);
  aic_trunc(isnan(dum))=NaN;
  for n=1:m
    ind = inds{n};
    a = y(:,indt)*x(ind,indt)'*inv(x(ind,indt)*x(ind,indt)');
    yall(:,n) = a*x(ind,it);
  end
  for ii=1:my*mx
    yhat_cva(ii,it) = yall(ii,imin(ii));
  end
end

error_cva = sum((y-yhat_cva).^2,2);
my_pcolor(X,Y,1-reshape(sqrt(error_cva/nt)./std(y,[],2),my,mx).^2')   
cbar
pacific_noresize_at_all
title('Cross-validation rms')

cor_cva = correlation(y,yhat_cva);
my_pcolor(X,Y,reshape(cor_cva,my,mx)')   
cbar
pacific_noresize_at_all
title('Cross-validation correlation')
pause


plot(cor_best(:),cor_cva(:),'.') 
grid
pause

plot(dum(:)./std(y,[],2).^2/nt,error_cva(:)./std(y,[],2).^2/nt,'.') 
grid
pause




return


%split data

indt = 1:(nt/2);
a = y(:,indt)*x(:,indt)'*inv(x(:,indt)*x(:,indt)');
for n=1:m
  ind = inds{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',[0.5:0.5:(m+0.5)])
pacific_noresize_at_all
title('predict 1981-2001')
pause


indt =(nt/2+1):nt;
a = y(:,indt)*x(:,indt)'*inv(x(:,indt)*x(:,indt)');
for n=1:m
  ind = inds{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',[0.5:0.5:(m+0.5)])
pacific_noresize_at_all
title('Predict 1960-1980')
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))]