
% VERIFY GLOBAL MEAN TEMPREATURE RECONSTRUCTION (QUICK RECON VERIFICATION)

addpath('../')
addpath(genpath('boundedline-pkg/'))
load('./colormaps/distinct_colors_20.mat')

% ALL PROXY TYPES, LMR DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2604_swtchbld100_21-Dec-2017_16:05:26.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2598_swtchbld100_21-Dec-2017_16:19:04.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2588_swtchbld100_21-Dec-2017_16:30:40.mat';

% INDEX ONLY RECONSTRUCTION WITH FULL ENSEMBLES
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:11:15.mat';


% ALL PROXY TYPES, PAGES2K DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_3_prxtp_tca_430_swtchbld100_28-Sep-2017_20:59:59.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_3_prxtp_tca_430_swtchbld100_28-Sep-2017_20:55:04.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_3_prxtp_tca_412_swtchbld100_28-Sep-2017_20:51:19.mat';

% TEST RECONSTRUCTIONS WITH ONLY TREES AND PDSI AS THE PSM
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_t_2165_swtchbld100_24-Sep-2017_16:40:12.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_t_2090_swtchbld100_24-Sep-2017_16:47:39.mat';


%-----------------------------------------

% ALL PROXY TYPES, LMR DATABASE
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';

disp(' ')
disp(['Analyzing file: ' efl])


%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])
load(['./output-recon/',efl])

% new path for figures
prntpth='./figs/';

%load([outpath,fln])

% Save figures?
prntfg='n';

% Years to compute the skill metrics over?
%analyrs=1871:1919;
analyrs=1871:2000;
disp(['Analysis years: ' num2str(analyrs(1)) '-' num2str(analyrs(end))])

% Get the global mean and ensemble

if ismember('2',state_tp)
    k=find(state_tp=='2');
    Xa_gmt_m=Xa_m(xmeta{k}.id_X(1),:)';
   
    sigma_t=Xa_sigma(xmeta{k}.id_X(1),:);
 
    
    %================================================================
    %                       FIGURES
    %================================================================
    
    
    % MEAN TEMPERATURE OVER THE INSTRUMENTAL PERIOD
    load('/d1/nsteiger/climate-data/hadcrut/hadcrut4_4_monthly_global_mean_T.mat')
    tas_had=mon_to_ann(h_median,mon_avg_o,mon_avg_f);
    tas_u=mon_to_ann(h_upprbnd,mon_avg_o,mon_avg_f);
    tas_l=mon_to_ann(h_lwrbnd,mon_avg_o,mon_avg_f);
    
    
    % get years of overlap
    h_yrs=1850:(1850+length(tas_had)-1);
    %[~,ir,ih] = intersect(r_o:r_f,h_yrs);
    ovrlprng=rng_intr([r_o,r_f],[h_yrs(1),h_yrs(end)]);
   
    [~,~,ir] = intersect(analyrs,r_o:r_f); % indices of recon segment
    [~,~,ih] = intersect(analyrs,h_yrs); % indices of obs

    % Prior reference (for CRPSS)
    xb_gmt=Xb_o(xmeta{k}.id_X(1),:);
    xb_gmt=xb_gmt-mean(xb_gmt)+mean(tas_had(ih));

    % make recon have same reference mean as HadCRUT4
    xa_gmt=Xa_gmt_m-mean(Xa_gmt_m(ir));
    xa_gmt=xa_gmt+mean(tas_had(ih));
    %Xa_gmt_rm=Xa_gmt_rm+mean(tas_had);

    x_rec=r_o:r_f;x_rec=x_rec(ir);

    %=======================
    % PLOT THE LINE FIGURES
    %=======================

    %plot_line_verif(y_rec,x_rec,x_sig,xb_rec,y_obs,ttl,xlbl,ylbl,prntfg,prntpth,prntnm,cohyn,dt)    
    plot_line_verif(xa_gmt(ir),x_rec,sigma_t(ir),xb_gmt,tas_had(ih),...
       'Global mean temperature','Years','Temperature anom. (deg C)',...
       prntfg,prntpth,['gmt_' fln(1:end-4)],'y',1,'northwest')    
    
    % Compute and display CE
    if length(analyrs)>50
    [~,~,vidx] = intersect(1871:1919,h_yrs); % indices of obs
    CE=nanCEv2(tas_had(ih),xa_gmt(ir),vidx);
    disp(['T2m: CE = ' sprintf('%0.2f',CE)])
    %disp(['T2m: CE0 = ' sprintf('%0.2f',nan_ce(tas_had(ih),xa_gmt(ir)))])
    end
    disp(['T2m: r2 = ' sprintf('%0.2f',nancorr(tas_had(ih),xa_gmt(ir))^2)]);
    disp(['T2m: CE0 = ' sprintf('%0.2f',nan_ce(tas_had(ih),xa_gmt(ir)))])
end


if ismember('w',state_tp) % North American Southwest drought
    
    k=find(state_tp=='w');
    xa_gmt=Xa_m(xmeta{k}.id_X(1),:)';
    sigma_t=Xa_sigma(xmeta{k}.id_X(1),:);
    
    xa_gmt=xa_gmt-mean(xa_gmt);
    
    
    %================================================================
    %                       FIGURES
    %================================================================
    
    
    % North American Southwest region: lat = 25-->42.5, lon = 125-->105 W
    
    % Load data and manipulate
    % flip lat and time
    pdsi=flip(flip(permute(ncread('/d1/nsteiger/drought-atlas/old_drought_atlases/NADAv2-2008.nc','PDSI'),[2 1 3]),1),3);
    lat=flip(ncread('/d1/nsteiger/drought-atlas/old_drought_atlases/NADAv2-2008.nc','lat'));
    lon=ncread('/d1/nsteiger/drought-atlas/old_drought_atlases/NADAv2-2008.nc','lon');
    time=flip(ncread('/d1/nsteiger/drought-atlas/old_drought_atlases/NADAv2-2008.nc','time'));
    
    % Just take over the region
    naswpdsi=pdsi(4:11,6:14,:);
    A=cosd(repmat(lat(4:11),[1 length(6:14)]));
    % Compute the index
    nasw_pdsi_idx=wmean_a(naswpdsi,A);
    
    % make recon have same reference mean as HadCRUT4
    xa_gmt=xa_gmt+mean(nasw_pdsi_idx);
    %Xa_gmt_rm=Xa_gmt_rm+mean(tas_had);
    
    xb_gmt=Xb_o(xmeta{k}.id_X(1),:);
    xb_gmt=xb_gmt-mean(xb_gmt)+mean(nasw_pdsi_idx(ih));

    % get years of overlap
    %[~,ir,ih] = intersect(r_o:r_f,time);
    [~,~,ir] = intersect(analyrs,r_o:r_f); % indices of recon segment
    [~,~,ih] = intersect(analyrs,time); % indices of obs
    x_rec=r_o:r_f;x_rec=x_rec(ir);

    %=======================
    % PLOT THE LINE FIGURES
    %=======================

    %plot_line_verif(y_rec,x_rec,x_sig,xb_rec,y_obs,ttl,xlbl,ylbl,prntfg,prntpth,prntnm,cohyn,dt)    
    plot_line_verif(xa_gmt(ir),x_rec,sigma_t(ir),xb_gmt,nasw_pdsi_idx(ih),...
       'NASW PDSI','Years','JJA PDSI',...
       prntfg,prntpth,['nasw_' fln(1:end-4)],'y',1,'northwest')    
    
    % Compute and display CE
    if length(analyrs)>50
    [~,~,vidx] = intersect(1871:1919,time); % indices of obs
    CE=nanCEv2(nasw_pdsi_idx(ih),xa_gmt(ir),vidx);
    disp(['NASW PDSI: CE = ' sprintf('%0.2f',CE)])
    %disp(['NASW PDSI: CE0 = ' sprintf('%0.2f',nan_ce(nasw_pdsi_idx(ih),xa_gmt(ir)))])
    end
    
    disp(['NASW PDSI: CE0 = ' sprintf('%0.2f',nan_ce(nasw_pdsi_idx(ih),xa_gmt(ir)))])
end


if ismember('n',state_tp)
    
    k=find(state_tp=='n');

    % Load observational data
    %[X_obs,xlat,xlon,time0] = load_obs_Smon('soda_sst');
    [X_obs,xlat,xlon,time0] = load_obs_Smon('bearth');
    X_obs=reshape(X_obs,length(xlat),length(xlon),length(time0)*12);

   for i=1:length(xmeta{k}.indinfo)
    

       xa_nino=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
       %xa_nino=Xa_nino_m-mean(Xa_nino_m);
       %Xa_gmt_rm=Xa_gmt_ens-repmat(mean(Xa_gmt_ens,1),[reconYrs,1]);
       
       sigma_t=reshape(Xa_sigma((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
       
       
       % Historical Nino data
       %load('/d1/nsteiger/climate-data/climate-indices/nino34_1950_2015.mat')
       %time=nino34(:,1);
       %tnd=load('/d1/nsteiger/climate-data/climate-indices/nino34_NOAA_HadISST_1870_2016.mat');
       %time=tnd.nino34(:,1);
       %t_nino=reshape(tnd.nino34(:,2:end)',12*length(time),1);

       [t_nino0]=nino(X_obs,xlat,xlon,i,'n');
       % start at 1870 instead of begining of record	
       %si=find(time0==1871);
       si=find(time0==analyrs(1)); % use first year of analysis period
       time=time0(si:end);
       t_nino=t_nino0((si*12-11):end);

       %t_nino=reshape(nino34(:,2:end)',12*length(time),1);
       % Remove the seasonal cycle
       %t_nino1=bsxfun(@minus,nino34(:,2:end),mean(nino34(:,2:end),1));
       %t_nino=reshape(t_nino1',12*length(time),1);
       
       
       % make recon have same reference mean
       %t_nino=t_nino-mean(t_nino);
       %Xa_gmt_rm=Xa_gmt_rm+mean(tas_had);
       
       % monthly indices
       %ovrlprng=rng_intr([r_o,r_f],[time(1),time(end)]);
       ovrlprng=rng_intr([r_o,r_f],[analyrs(1),analyrs(end)]);
       mon_r=r_o:(1/12):(r_f+11/12);
       mon_t=time(1):(1/12):(time(end)+11/12);
       
       % get years of overlap
       ir1=find(ovrlprng(1)==mon_r);
       ir2=find((ovrlprng(2)-1/12)==mon_r)+12;
       ih1=find(ovrlprng(1)==mon_t);
       ih2=find((ovrlprng(2)-1/12)==mon_t)+12;
       
       % Just look at skill of anomalies
       xa_nino=xa_nino-mean(xa_nino)+mean(t_nino(ih1:ih2));
       
       xb_nino=reshape(Xb_o((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*length(p_yrs),1);
       xb_nino=xb_nino-mean(xb_nino)+mean(t_nino(ih1:ih2));

       x_rec=mon_r(ir1:ir2);

       %=======================
       % PLOT THE LINE FIGURES
       %=======================

       %plot_line_verif(y_rec,x_rec,x_sig,xb_rec,y_obs,ttl,xlbl,ylbl,prntfg,prntpth,prntnm,cohyn,dt)    
       plot_line_verif(xa_nino(ir1:ir2),x_rec,sigma_t(ir1:ir2),xb_nino,t_nino(ih1:ih2),...
          xmeta{k}.indinfo{i}.indxnm,'Years','Temperature (deg C)',...
	  prntfg,prntpth,['nino_index_' num2str(i) '_' fln(1:end-4)],'y',1/12,'northwest')    
       
       
       % Compute and display CE
       if length(analyrs)>50
       ovrlprng=rng_intr([r_o,r_f],[1871,1919]);
       ihv1=find(ovrlprng(1)==mon_t);
       ihv2=find((ovrlprng(2)-1/12)==mon_t)+12;
       CE=nanCEv2(t_nino(ih1:ih2),xa_nino(ir1:ir2),ihv1:ihv2);
       disp([xmeta{k}.indinfo{i}.indxnm ': CE = ' sprintf('%0.2f',CE)])
       %disp([xmeta{k}.indinfo{i}.indxnm ': CE0 = ' sprintf('%0.2f',nan_ce(t_nino(ih1:ih2),xa_nino(ir1:ir2)))])
       end
 
       disp([xmeta{k}.indinfo{i}.indxnm ': CE0 = ' sprintf('%0.2f',nan_ce(t_nino(ih1:ih2),xa_nino(ir1:ir2)))])
       disp([xmeta{k}.indinfo{i}.indxnm ': CE0 ann = ' sprintf('%0.2f',nan_ce(mon2ann(t_nino(ih1:ih2),1,12),...
         mon2ann(xa_nino(ir1:ir2),1,12)))])
       % NORMAL SPECTRA
       
       % [Pxy,F] = cpsd(x,y,window,noverlap,nfft,fs)
       % [pxx,f] = periodogram(x,window,f,fs)
       % [pxx,f] = pwelch(x,window,noverlap,f,fs)
       %[p_t,~]=periodogram(detrend(t_nino(ih)),hanning(length(t_nino(ih))),[],12);
       %[p_r,f]=periodogram(detrend(xa_nino(ir)),hanning(length(t_nino(ih))),[],12);
       
       % UNCOMMENT BELOW HERE FOR COMPUTING SPECTRA
       % [p_t,~]=pwelch(detrend(t_nino(ih1:ih2)),round(length(t_nino(ih1:ih2))/8),[],[],12);
       % [p_r,f]=pwelch(detrend(xa_nino(ir1:ir2)),round(length(t_nino(ih1:ih2))/8),[],[],12);
       % 
       % figure
       % hold on
       % plot(f,real(p_t),'--','linewidth',2,'color',[0.5 0.5 0.5])
       % plot(f,real(p_r),'linewidth',2,'color',[0.494 0.184 0.556])
       % legend('Obs','Recon')
       % xlim([0 0.5]) % limit view to 100 years
       % ax = gca;
       % ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
       % ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
       % fig = gcf;fig.PaperUnits = 'inches';
       % fig.PaperPosition = [0 2 8.5 4.5]; % [left bottom width height]
       % box on
       % title([xmeta{k}.indinfo{i}.indxnm ' index'],'fontsize',18,'fontweight','bold')
       % xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
       % ylabel('Power spectral density','fontsize',14)
       % set(gca,'fontsize',14)
       % hold off
       % 
       % if make_figs=='y';print([pth 'nino_spectra_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
    
   end    
    
end


% AMO VERIFICATION
if ismember('o',state_tp)
    
    k=find(state_tp=='o');
    Xa_gmt_m=Xa_m(xmeta{k}.id_X(1),:)';
   
    sigma_t=Xa_sigma(xmeta{k}.id_X(1),:);
    
    % Global mean from update step
    xa_gmt=Xa_gmt_m-mean(Xa_gmt_m);
    %Xa_gmt_rm=Xa_gmt_ens-repmat(mean(Xa_gmt_ens,1),[reconYrs,1]);
    
    
    %================================================================
    %                       FIGURES
    %================================================================
    
    
    
    % Load data and manipulate
    % flip lat and time
    load('/d1/nsteiger/climate-data/climate-indices/AMO_KaplanSST.mat','amo_mon','ayrs')
    amo_mon1=amo_mon';
    amo_mon1=amo_mon1(:); % convert from original format to monthly time series
    amo=mon2ann(amo_mon1,mon_avg_o,mon_avg_f);
    ayrs=ayrs(1:length(amo));
    
    % make recon have same reference mean as HadCRUT4
    xa_gmt=xa_gmt+mean(amo);
    %Xa_gmt_rm=Xa_gmt_rm+mean(tas_had);

    xb_gmt=Xb_o(xmeta{k}.id_X(1),:);
    xb_gmt=xb_gmt-mean(xb_gmt)+mean(amo);
    
    % get years of overlap
    %[~,ir,ih] = intersect(r_o:r_f,ayrs);
    ovrlprng=rng_intr([r_o,r_f],[ayrs(1),ayrs(end)]);
    [~,~,ir] = intersect(analyrs,r_o:r_f); % indices of recon segment
    [~,~,ih] = intersect(analyrs,ayrs); % indices of obs
    
    x_rec=r_o:r_f;x_rec=x_rec(ir);

    %=======================
    % PLOT THE LINE FIGURES
    %=======================

    %plot_line_verif(y_rec,x_rec,x_sig,xb_rec,y_obs,ttl,xlbl,ylbl,prntfg,prntpth,prntnm,cohyn,dt)    
    plot_line_verif(xa_gmt(ir),x_rec,sigma_t(ir),xb_gmt,amo(ih),...
       'AMO (NASST index)','Years','Temperature (deg C)',...
       prntfg,prntpth,['amo_' fln(1:end-4)],'y',1,'northwest')    

    % Compute and display CE
    if length(analyrs)>50
    [~,~,vidx] = intersect(1871:1919,ayrs); % indices of obs
    CE=nanCEv2(amo(ih),xa_gmt(ir),vidx);
    disp(['AMO: CE = ' sprintf('%0.2f',CE)])
    %disp(['AMO: CE0 = ' sprintf('%0.2f',nan_ce(amo(ih),xa_gmt(ir)))])
    end

    disp(['AMO: r2 = ' sprintf('%0.2f',nancorr(amo(ih),xa_gmt(ir))^2)]);
    disp(['AMO: CE0 = ' sprintf('%0.2f',nan_ce(amo(ih),xa_gmt(ir)))])
end



%----------- ITCZ verification --------------

if ismember('z',state_tp)
    
    
    
    %================================================================
    %                       FIGURES
    %================================================================
    
    
    % Load data and manipulate
    % flip lat and time
    Pr_mon=permute(ncread('/d1/nsteiger/climate-data/gpcp/gpcp_cdr_v23rB1_1979_2016.nc','precip'),[2 1 3]);
    prlat=double(ncread('/d1/nsteiger/climate-data/gpcp/gpcp_cdr_v23rB1_1979_2016.nc','latitude'));
    prlon=double(ncread('/d1/nsteiger/climate-data/gpcp/gpcp_cdr_v23rB1_1979_2016.nc','longitude'));
    pryrs=1979:2016;
    
    Pr=mon2ann(Pr_mon,mon_avg_o,mon_avg_f);
    pryrs=pryrs(1:size(Pr,3));% account for seasonal averaging
    
    k=find(state_tp=='z');
    for i=1:length(xmeta{k}.indinfo)
      xa_var_m=Xa_m(xmeta{k}.id_X(1)+i-1:xmeta{k}.id_X(1)+i-1,:);
      sigma_t=Xa_sigma(xmeta{k}.id_X(1)+i-1:xmeta{k}.id_X(1)+i-1,:);

      %------------------------------
      % Scaled reconstruction state
      %------------------------------

      % Remove time mean
      xa_var=bsxfun(@minus,xa_var_m,mean(xa_var_m,2));
      %Xa_gmt_rm=Xa_gmt_ens-repmat(mean(Xa_gmt_ens,1),[reconYrs,1]);

      % ITCZ location over the ocean regions
      %Atlantic = [311 350]
      %ITCZ(1,:) = itcz_lon(Pr,prlat,prlon,[311 350]);
      %Pacific = [128 278]
      %ITCZ(2,:) = itcz_lon(Pr,prlat,prlon,[128 278]);
      %Indian = [50 96]
      %ITCZ(3,:) = itcz_lon(Pr,prlat,prlon,[50 96]);

      %tvar=ITCZ;
      %itcznms={'Atlantic','Pacific','Indian'};

      tvar=itcz_lon(Pr,prlat,prlon,[xmeta{k}.indinfo{i}.lon1 xmeta{k}.indinfo{i}.lon2]);


      % make recon have same reference mean as HadCRUT4
      xa_var=bsxfun(@plus,xa_var,mean(tvar,2));
      %Xa_gmt_rm=Xa_gmt_rm+mean(tas_had);

      % get years of overlap
      [~,ir,ih] = intersect(r_o:r_f,pryrs);

      xb_gmt=Xb_o(xmeta{k}.id_X(1)+i-1,:);
      xb_gmt=xb_gmt-mean(xb_gmt)+mean(tvar);

       x_rec=r_o:r_f;x_rec=x_rec(ir);

       %=======================
       % PLOT THE LINE FIGURES
       %=======================

       %plot_line_verif(y_rec,x_rec,x_sig,xb_rec,y_obs,ttl,xlbl,ylbl,prntfg,prntpth,prntnm,cohyn,dt)    
       plot_line_verif(xa_var(ir),x_rec,sigma_t(ir),xb_gmt,tvar(ih),...
	  xmeta{k}.indinfo{i}.varnmlng,'Years','Latitude (deg)',...
	  prntfg,prntpth,['itcz_' num2str(i) '_' fln(1:end-4)],'n',1,'northwest')    

    end
end

































