m=length(inds);
colormap(jet(m))
C=colormap;
C(1,:) = [1,1,1]*0.8;
colormap(C);

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]);

%mallow's cp
cp = sse./sse(:,m*ones(1,m))*(nt-m) - 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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('Mallow''s C_p')

pause
if make_pic
  filename=[casename,'_pred_sel_cp'];
  export_fig('-pdf',filename)
end

%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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('AIC')

pause
if make_pic
  filename=[casename,'_pred_sel_aic'];
  export_fig('-pdf',filename)
end



%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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('corrected AIC')

pause
if make_pic
  filename=[casename,'_pred_sel_aicc'];
  export_fig('-pdf',filename)
end




%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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('BIC')

pause
if make_pic
  filename=[casename,'_pred_sel_bic'];
  export_fig('-pdf',filename)
end



%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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('F-test')

pause
if make_pic
  filename=[casename,'_pred_sel_ftest'];
  export_fig('-pdf',filename)
end


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

my_pcolor(X,Y,fc_trunc')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('F-test corrected')

pause
if make_pic
  filename=[casename,'_pred_sel_ftest_corr'];
  export_fig('-pdf',filename)
end



%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
  sse_cv(:,n) = mean(res_cv.^2,2);
  sse_cv_std(:,n) = std(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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('Cross-validation')

pause
if make_pic
  filename=[casename,'_pred_sel_cv'];
  export_fig('-pdf',filename)
end


%cv + interval = hard
for ii=1:my*mx
  i1 = imin(ii);
  if (sse_cv(ii,1) - sse_cv_std(ii,1)/sqrt(nt)) < sse_cv(ii,i1) ; 
    imin(ii) = 1;
  end
end

%my_pcolor(X,Y,reshape(imin,my,mx)')   
%cbar('vert',[0.5:0.5:(m+0.5)])
%pacific
%title('Cross-validation + interval')
%pause
%if make_pic
%  filename=[casename,'_pred_sel_cp'];
%  export_fig('-pdf',filename)
%end


%split data

indt1 = 1:(nt/2);
indt2 = (nt/2+1):nt;
for n=1:m
  ind = inds{n};
  a = y(:,indt1)*x(ind,indt1)'*inv(x(ind,indt1)*x(ind,indt1)');
  sse2(:,n) = sum((y(:,indt2)-a*x(ind,indt2)).^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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('predict second half')

pause
if make_pic
  filename=[casename,'_pred_sel_split2'];
  export_fig('-pdf',filename)
end


indt2 = 1:(nt/2);
indt1 = (nt/2+1):nt;
for n=1:m
  ind = inds{n};
  a = y(:,indt1)*x(ind,indt1)'*inv(x(ind,indt1)*x(ind,indt1)');
  sse2(:,n) = sum((y(:,indt2)-a*x(ind,indt2)).^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')   
h=cbar('vert',[0.5:0.5:(m+0.5)]);
set(h,'ytick',[1:m],'yticklabel',lbs,'fontsize',8)
pacific
title('predict first half')

pause
if make_pic
  filename=[casename,'_pred_sel_split1'];
  export_fig('-pdf',filename)
end



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

[mean(aicc_trunc(iok)<bic_trunc(iok)), ...
mean(aicc_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))]

hist([aic_trunc(:),cp_trunc(:),aicc_trunc(:),cv_trunc(:),bic_trunc(:),f_trunc(:),fc_trunc(:)])
legend('aic','cp','aicc','cv','bic','F','Fc','location','best')
colormap(lines(7))
set(gca,'xtick',[1:m])
set(gca,'xticklabel',lbs)
pause
if make_pic
  filename=[casename,'_pred_sel_histogram'];
  export_fig('-pdf',filename)
end


mean((bic_trunc(iok)==f_trunc(iok)))

mean((aic_trunc(iok)==aicc_trunc(iok))&(aic_trunc(iok)==cv_trunc(iok))&(aic_trunc(iok)==cp_trunc(iok)))

