
addpath('../')

%===========================================================================
% LOAD CLIMATE INDICES

% Load full ensembles for the nino and amo variables
efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:24:35.mat';
% LOAD JJA FOR AMO ANALYSIS
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:11:15.mat';

pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
%load([pth,efl],'r_o','r_f','state_tp','xmeta','Xa_m','Xa_ens','reconYrs','p_yrs')
load([pth,efl])


% What years to look at?
analyrs1=1000:1925;
%analyrs2=1801:1925;
analyrs3=1000:2000;
%analyrs4=1000:1925;
%analyrs4=1000:1800;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
%[~,~,ir2] = intersect(analyrs2,r_o:r_f); % indices of recon segment
[~,~,ir3] = intersect(analyrs3,r_o:r_f); % indices of recon segment
%[~,~,ir4] = intersect(analyrs4,r_o:r_f); % indices of recon segment
% monthly indices
ovrlprng=rng_intr([r_o,r_f],[analyrs1(1),analyrs1(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm1_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm1_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

%ovrlprng=rng_intr([r_o,r_f],[analyrs2(1),analyrs2(end)]);
%mon_r=r_o:(1/12):(r_f+11/12);
%irm2_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
%irm2_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

ovrlprng=rng_intr([r_o,r_f],[analyrs3(1),analyrs3(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm3_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm3_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;



k=find(state_tp=='o');
xa_amo_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);
% Ensemble mean of AMO
xa_amo_mn=mean(xa_amo_ens,2);
amo_off=mean(mean(xa_amo_ens(ir1,:),2));


k=find(state_tp=='2');
xa_gmt_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);

k=find(state_tp=='n'); % specific nino variable
if strcmp(xmeta{k}.indinfo{3}.indxnm,'Nino 3.4'); i=3; else; error('Wrong Nino specification'); end
xa_nino_ens=reshape(permute(Xa_ens((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:,:),[1 3 2]),12*reconYrs,length(p_yrs));
xa_nino_mn=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
% Nino offset with respect to reference period
nino_off=mean(mean(xa_nino_ens(irm1_a:irm1_b,:),2));

% Annual nino from monthly data (leaves year ends not quite correct)
xa_nino_ann=mon2ann(xa_nino_mn,4,3);
% Annual nino from spatial field
%xa_nino_ann=nino(Xt2m,lat,lon,3,'n');



% Get annual Nino 3.4 value from annual spatial field
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';
load([pth,efl],'state_tp','xmeta','Xa_m')
k=find(state_tp=='t');
lat=xmeta{k}.lat;
lon=xmeta{k}.lon;
Xt2m=reshape(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:),length(lat),length(lon),reconYrs);

% Load mask for Chile/Argentina
%msk=load('ChlArg_mask.mat');
%msk=load('Pampas_mask.mat');
%lndmsk=msk.sftlf;
%lndmsk(lndmsk==0)=NaN;

% Land mask over Argentina
%msk=load('argentina_mask.mat');
%lndmsk=msk.argntn;
% Land mask over Chile
msk=load('chile_mask.mat');
lndmsk=msk.chl;
% Land mask over Chile + Argentina
%msk=load('chile+argentina_mask.mat');
%lndmsk=msk.chlarg;




% Latitude weighting
A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));

% Get the ensemble mean; identical to the original reconstruction
Xv_msk=bsxfun(@times,Xt2m,lndmsk);
t2m_casa=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX




% Drought atlases use JJA instead of annual
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
%pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])
%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])

% PDSI ONLY, ENSEMBLE RECONSTRUCTION
%efl='cesm_lme010_r12000_p998_state_q_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_06-Aug-2018_17:32:52.mat';
%pth='./output-recon/';load([pth,efl],'xmeta')


% Load DJF PDSI
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
%load([pth,efl],'state_tp','xmeta','Xa_m')

% Load A2M PDSI
k=find(state_tp=='q');
Xvar=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
Xpdsi(xmeta{k}.lndidx,:)=Xvar;
Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);



%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

% Load both south american and north america drought indices
%load('svmd_ens_90pagree_sa.mat')
%load('svmd_ens_90pagree_pampas.mat')
%load('pampas_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('argentina_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
load('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chlarg_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa=Xa_di;
xa_di_ens_sa=Xa_di_ens;
ds_sa=svmd;
mdagr_sa=ens_drght;
sa_yrs=analyrs;

%load('svmd_ens_50pagree.mat')
%load('svmd_ens_90pagree.mat')
load('nasw_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_na=Xa_di;
xa_di_ens_na=Xa_di_ens;
ds_na=svmd;
mdagr_na=ens_drght;
na_yrs=analyrs;

%ds=svmd;


%===========================================================================
% LOAD FORCED CLIMATE RESPONSE ESTIMATES

% Full forcing
%load('../megad-forcing/output_twoboxCM_full_06-Aug-2018_15:09:13.mat','gmt_fr','myrs')
load('../ebm-held/output_twoboxCM_full_21-Sep-2018_13:29:49.mat','gmt_fr','myrs')


% Plot mean forced response
gmt_est_mn=mean(gmt_fr,2);
[~,~,ig2] = intersect(analyrs1,myrs); 
%[~,~,ig3] = intersect(analyrs3,myrs); 
%[~,~,ig4] = intersect(analyrs4,myrs); 

% remove mean from ref
gmt_est_ens=gmt_fr-mean(gmt_est_mn(ig2));
gmt_est_f=gmt_est_mn-mean(gmt_est_mn(ig2));

% solar 
load('../ebm-held/output_twoboxCM_sol_21-Sep-2018_13:29:03.mat','gmt_fr')
%load('../megad-forcing/output_twoboxCM_solar_08-Aug-2018_13:47:28.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_s=gmt_est_mn-mean(gmt_est_mn(ig2));

% volcanic 
%load('../megad-forcing/output_twoboxCM_volc_08-Aug-2018_13:38:26.mat','gmt_fr')
load('../ebm-held/output_twoboxCM_volc_21-Sep-2018_13:35:31.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_v=gmt_est_mn-mean(gmt_est_mn(ig2));

% ghg
%load('../megad-forcing/output_twoboxCM_ghg_08-Aug-2018_14:05:19.mat','gmt_fr')
load('../ebm-held/output_twoboxCM_ghg_21-Sep-2018_13:32:35.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_g=gmt_est_mn-mean(gmt_est_mn(ig2));

% ghg + solar
load('../ebm-held/output_twoboxCM_solghg_26-Sep-2018_14:06:47.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_sg=gmt_est_mn-mean(gmt_est_mn(ig2));

% volc + solar
load('../ebm-held/output_twoboxCM_volcsol_26-Sep-2018_15:37:39.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_vs=gmt_est_mn-mean(gmt_est_mn(ig2));

%===========================================================================



load distinct_colors_20.mat

clrna=clrs.teal;
clrsa=clrs.maroon;


% PHYDA AND HIGHLIGHT MEGADROUGHTS

% Droughts in the real data mean
ds0_n=mdrought(xa_di_na(analyrs),1,'y');
ds0_s=mdrought(xa_di_sa(analyrs),1,'y');

% Droughts too frequent, throw out based on relative threshold
%dcut=67;
dcut=50;
np67=prctile(ds0_n(:,3),dcut);
rmv=find(ds0_n(:,3)>np67);
ds0_n(rmv,:)=[];

sp67=prctile(ds0_s(:,3),dcut);
rmv=find(ds0_s(:,3)>sp67);
ds0_s(rmv,:)=[];

% compute some drought stats
disp(['nasw drought severity - median = ' num2str(median(ds0_n(:,3)))])
disp(['chile drought severity - median = ' num2str(median(ds0_s(:,3)))])
disp(' ')
disp(['nasw drought severity - mean = ' num2str(mean(ds0_n(:,3)))])
disp(['chile drought severity - mean = ' num2str(mean(ds0_s(:,3)))])
disp(' ')
disp(['nasw drought length - median = ' num2str(median(ds0_n(:,2)-ds0_n(:,1)))])
disp(['nasw drought length - mean = ' num2str(mean(ds0_n(:,2)-ds0_n(:,1)))])
disp(' ')
disp(['chile drought length - median = ' num2str(median(ds0_s(:,2)-ds0_s(:,1)))])
disp(['chile drought length - mean = ' num2str(mean(ds0_s(:,2)-ds0_s(:,1)))])


figure
hold on

yMin=-7;
yMax=7;

%d=ds_na(5:end,:);% only plot droughts that begin after 1000
%d=ds_na;
d=ds0_n;
% Add a patches for top five times from truth
for i=1:length(d)
    p=patch(na_yrs(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrna);
    p.FaceAlpha=0.35;p.EdgeColor='none';
end

%d=ds_sa;
d=ds0_s;
% Add a patches for top five times from truth
for i=1:length(d)
    p=patch(sa_yrs(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrsa);
    p.FaceAlpha=0.35;p.EdgeColor='none';
end

naoff=mean(xa_di_na(ir3));
saoff=mean(xa_di_sa(ir3));


h1=plot(analyrs3,xa_di_na(ir3)-naoff,'linewidth',1,'color',[clrna 0.5]); % extra number on the end adds transparency
h2=plot(analyrs3,xa_di_sa(ir3)-saoff,'linewidth',1,'color',[clrsa 0.5]);
%plot(analyrs3,smooth(DA_di(id3),11),'linewidth',2,'color',clrs.red)
h3=plot(analyrs3,smooth(xa_di_na(ir3),11)-naoff,'linewidth',2,'color',clrna);
h4=plot(analyrs3,smooth(xa_di_sa(ir3),11)-saoff,'linewidth',2,'color',clrsa);
plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs3(1) analyrs3(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%y=Xa_di(ir3); % using 'l' from above
%yMin=min(y)-0.5;yMax=max(y)+0.5;
%ax=gca;yMin=ax.YLim(1);yMax=ax.YLim(2);
ylim([yMin yMax])


% % Plot patches over high drought agreement periods
% dyrs_sa=find(mdagr_sa>=80);
% dyrs_na=find(mdagr_na>=80);
% 
% % Find continguous years
% [~,~,irn] = intersect(na_yrs,r_o:r_f); % indices of recon segment
% [ds_na]=contig(dyrs_na,irn);
% 
% [~,~,irs] = intersect(sa_yrs,r_o:r_f); % indices of recon segment
% [ds_sa]=contig(dyrs_sa,irs);





legend([h3,h4],{'NASW','SASW'},'location','southeast')
%legend([h3,h4],{'NASW','Argentina'},'location','southeast')
%legend([h3,h4],{'NASW','ChlArg'},'location','southeast')
hold off

print(['/home/nsteiger/Dropbox/drght_nasw_chile_phyda_11mean_tealmaroon.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_na_chile_phyda_11mean_mean_p' num2str(dcut) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_na_argentina_phyda_11mean_mean_p' num2str(dcut) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_na_chlarg_phyda_11mean_mean_p' num2str(dcut) '.png'],'-dpng','-r300');


return


% AR null for mean
addpath('../ar-fit/')
ntr=1000;
d0=xa_di_na(analyrs);
[w,A,C,~,~,~]=arfit(d0,1,5); % explore up to 5, but it looks like all are 2nd order
na_ar_ens=squeeze(arsim(w,A,C,[length(analyrs),ntr]));% generate base ensemble AR noise time series 

d0=xa_di_sa(analyrs);
[w,A,C,~,~,~]=arfit(d0,1,5); % explore up to 5, but it looks like all are 2nd order
sa_ar_ens=squeeze(arsim(w,A,C,[length(analyrs),ntr]));% generate base ensemble AR noise time series 

% Ensemble estimate of drought periods
dsens_n=struct();
dsens_s=struct();
for i=1:ntr
   % Compute the megadroughts
   dsens_n.ds{i}=mdrought(na_ar_ens(:,i),1,'y');
   dsens_s.ds{i}=mdrought(sa_ar_ens(:,i),1,'y');
   
   np67_rn=prctile(dsens_n.ds{i}(:,3),dcut);% Throw out droughts using noise data (self-consistent)
   rmv=find(dsens_n.ds{i}(:,3)>np67_rn);
   %rmv=find(dsens_n.ds{i}(:,3)>np67);% Throw out droughts below Xth percentile from real data
   dsens_n.ds{i}(rmv,:)=[];
   
   sp67_rn=prctile(dsens_s.ds{i}(:,3),dcut);% Throw out droughts using noise data (self-consistent)
   rmv=find(dsens_s.ds{i}(:,3)>sp67_rn);
   %rmv=find(dsens_s.ds{i}(:,3)>sp67);% Throw out droughts below Xth percentile from real data
   dsens_s.ds{i}(rmv,:)=[];
end


% Find which years there is agreement
pd_nthr=zeros(ntr,1);
pd_na=zeros(ntr,1);
pd_sa=zeros(ntr,1);
pd_bh=zeros(ntr,1);
for i=1:ntr

   % North America
   dn0=dsens_n.ds{i}(:,1:2); % pull out dates of droughts
   dn0=dn0';dn0=dn0(:);

   % Fill out all years of droughts in vector
   k=1;dn1=[];
   while k<length(dn0)
      dn1=cat(2,dn1,dn0(k):dn0(k+1));
      k=k+2;
   end

   % South America
   ds0=dsens_s.ds{i}(:,1:2); % pull out dates of droughts
   ds0=ds0';ds0=ds0(:);

   % Fill out all years of droughts in vector
   k=1;ds1=[];
   while k<length(ds0)
      ds1=cat(2,ds1,ds0(k):ds0(k+1));
      k=k+2;
   end

   % Figure out indices of neither, NA, SA, both
   inthr=setdiff(1:length(analyrs),union(dn1,ds1)); % neither
   [~,ina,isa]=setxor(dn1,ds1); % indices of NA only and SA only
   ibth=intersect(dn1,ds1); % both

   % Probabilities for differention possibilities
   pd_nthr(i)=length(inthr)/length(analyrs);
   pd_na(i)=length(ina)/length(analyrs);
   pd_sa(i)=length(isa)/length(analyrs);
   pd_bh(i)=length(ibth)/length(analyrs);

end


%----------------------------
% Probabilites for real data
%----------------------------

% North America
dn0=ds0_n(:,1:2)';dn0=dn0(:);

% Fill out all years of droughts in vector
k=1;dn1=[];
while k<length(dn0)
   dn1=cat(2,dn1,dn0(k):dn0(k+1));
   k=k+2;
end

% South America
ds0=ds0_s(:,1:2)';ds0=ds0(:);

% Fill out all years of droughts in vector
k=1;ds1=[];
while k<length(ds0)
   ds1=cat(2,ds1,ds0(k):ds0(k+1));
   k=k+2;
end


% Figure out indices of neither, NA, SA, both
inthr=setdiff(1:length(analyrs),union(dn1,ds1)); % neither
[~,ina,isa]=setxor(dn1,ds1); % indices of NA only and SA only
ibth=intersect(dn1,ds1); % both

% Probabilities for differention possibilities
pd_nthr0=length(inthr)/length(analyrs);
pd_na0=length(ina)/length(analyrs);
pd_sa0=length(isa)/length(analyrs);
pd_bh0=length(ibth)/length(analyrs);


% Combine noise null
pnl(:,1)=pd_bh;
pnl(:,2)=pd_nthr;
pnl(:,3)=pd_na;
pnl(:,4)=pd_sa;

%c0={'Both','None','NASW','Pampas'};
c0={'Both','Neither','NASW','SASW'};
%c0={'Both','None','NASW','Arg'};
%c0={'Both','None','NASW','ChlArg'};
c=categorical(c0,c0);

% Colors for boxplots
%cm2=cat(1,clrs.red,clrs.grey,clrs.blue,clrs.purple);
cm2=cat(1,clrs.orange,clrs.grey);


figure
hold on
%hl=boxplot(pnl,c,'symbol','','colors',cm2,'whisker',1); % whisker=1 => 2 sigma coverage
hl=boxplot(pnl(:,1:2),c(1:2),'symbol','','colors',cm2,'whisker',1); % whisker=1 => 2 sigma coverage
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
ylim([0 0.7])
%line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
ylabel('Probability','fontsize',18)
set(gca,'FontSize',16,'XTickLabelRotation',45)

h1=plot(1,pd_bh0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
h2=plot(2,pd_nthr0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
% only show first two options "both" and "none"
%h3=plot(3,pd_na0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
%h4=plot(4,pd_sa0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);

%legend([h1,h2,h3,h4],c0,'location','northeast')
%legend([h1],'PHYDA','location','northeast')

fig = gcf;fig.PaperUnits = 'inches';
%fig.PaperPosition = [0 2 4 4]; % [left bottom width height]
fig.PaperPosition = [0 2 2 4]; % [left bottom width height]


%print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_null_na_chile_phyda_p' num2str(dcut) '_orange.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_null_na_argentina_phyda_p' num2str(dcut) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_null_na_chlarg_phyda_p' num2str(dcut) '.png'],'-dpng','-r300');


return

%-------------------------------------------------
%  COMPARE PDFS DURING MEGADROUGHTS
%-------------------------------------------------

error('Need to update below with correct accounting of types...')

% PDFs for different kinds of years 
enso_na_pdf=xa_nino_ann();
enso_sa_pdf=xa_nino_ann(ds1);
enso_nne_pdf=xa_nino_ann(dnone);
enso_bth_pdf=xa_nino_ann(dbth);

x=(min(xa_nino_ann(analyrs1))-1):0.01:(max(xa_nino_ann(analyrs1))+1);

pk1=fitdist(enso_na_pdf(:),'kernel');
Y_na=pdf(pk1,x); 

pk1=fitdist(enso_sa_pdf(:),'kernel');
Y_sa=pdf(pk1,x); 

pk1=fitdist(enso_nne_pdf(:),'kernel');
Y_nne=pdf(pk1,x); 

pk1=fitdist(enso_bth_pdf(:),'kernel');
Y_bth=pdf(pk1,x); 

figure
hold on
plot(x,Y_na,'color',clrs.blue,'linewidth',3)
plot(x,Y_sa,'color',clrs.purple,'linewidth',3)
plot(x,Y_nne,'color',clrs.black,'linewidth',3)
plot(x,Y_bth,'color',clrs.red,'linewidth',3)
box on

%legend('NASW','Chile','None','Both')
%legend('NASW','Arg','None','Both')
legend('NASW','ChlArg','None','Both')

%print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_pdfs_na_chlarg_phyda_p' num2str(dcut) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_pdfs_na_argentina_phyda_p' num2str(dcut) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_pdfs_na_chile_phyda_p' num2str(dcut) '.png'],'-dpng','-r300');

return




% Get modern pdf and mean offset value
enso_mdrn_pdf=xa_nino_ens(irm1_a:irm1_b,:);
enp=struct([]);

%ds=ds1;
ird_all=[];
ird_ann=[];
figure; hold on;

% PDF of the reference period
enso_mdrn_pdf=enso_mdrn_pdf(:)-nino_off;
x=(min(enso_mdrn_pdf(:))-1):0.01:(max(enso_mdrn_pdf(:))+1);
pk1=fitdist(enso_mdrn_pdf(:),'kernel');
Yk1=pdf(pk1,x); 
%plot(x,Yk1,'color',clrs.blue,'linewidth',2.5)
%area(x,Yk1,'facecolor',clrs.cyan,'linestyle','none')

dn1
ds1
dnone=find(dtm==0); % no drought will still be zero
dbth=ibth; 

%for i=1:size(ds,1)
for i=1:4

  ovrlprng=rng_intr([r_o,r_f],[ds(i,1)+analyrs1(1)-1,ds(i,2)+analyrs1(1)-1]);
  irm_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
  irm_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

  % Collect all months of drought
  ird_all=cat(2,ird_all,irm_a:irm_b);

  [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),r_o:r_f); 
  % Collect all years of drought
  ird_ann=cat(1,ird_ann,irab);
  
  % Pull out distributions during this drought
  enso_pdf=xa_nino_ens(irm_a:irm_b,:)-nino_off; 
  enp{i}.pdf=enso_pdf(:);

  x2=(min(enso_pdf(:))-1):0.01:(max(enso_pdf(:))+1);
  pk2=fitdist(enso_pdf(:),'kernel');
  Yk2=pdf(pk2,x2); 
  %plot(x,Yk1,'color',[0.5 0.5 0.5],'linewidth',1.75)
  plot(x2,Yk2,'color',clrs.orange,'linewidth',1.75)
  title('Nino3.4')
  %xlabel('Temp. anom. (deg C)','fontsize',18)
  %ylabel('Probability density','fontsize',18)
  %set(gca,'fontsize',18)
  box on



end


% PDF of the non-drought years
irw = setdiff(irm1_a:irm1_b,ird_all);
ensow_pdf=xa_nino_ens(irw,:)-nino_off;
x_w=(min(ensow_pdf(:))-1):0.01:(max(ensow_pdf(:))+1);
pk_w=fitdist(ensow_pdf(:),'kernel');
Yk_w=pdf(pk_w,x_w); 
plot(x_w,Yk_w,'k','linewidth',3)

% Plot contemporary after so visible
plot(x,Yk1,'color',clrs.teal,'linewidth',3)

% Remove the labels, squash down the figure
xlim([-4.75 4.75])
set(gca,'xticklabel',[])
set(gca,'yticklabel',[])
set(gca,'XDir','reverse'); % Reverse direction for flipping on its side
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 4 2.5]; % [left bottom width height]

%print(['./figs/enso_droughtpdf_shifts.png'],'-dpng','-r300');


% Take analysis years
Xt2ma=Xt2m(:,:,analyrs);
Xpdsia=Xpdsi(:,:,analyrs);

% Composites over the possibilities
Xt2maM=bsxfun(@minus,Xt2ma,mean(Xt2ma,3));
X1=zeros(length(lat),length(lon),4);
X1(:,:,1)=mean(Xt2maM(:,:,df0),3);
X1(:,:,2)=mean(Xt2maM(:,:,dn1),3);
X1(:,:,3)=mean(Xt2maM(:,:,ds1),3);
X1(:,:,4)=mean(Xt2maM(:,:,ibth),3);

XpaM=bsxfun(@minus,Xpdsia,mean(Xpdsia,3));
X2=zeros(length(lat),length(lon),4);
X2(:,:,1)=mean(XpaM(:,:,df0),3);
X2(:,:,2)=mean(XpaM(:,:,dn1),3);
X2(:,:,3)=mean(XpaM(:,:,ds1),3);
X2(:,:,4)=mean(XpaM(:,:,ibth),3);

cnm={'none','na','sa','both'};

X2(isnan(X2))=1e20;






















return
















%------------------------------
% ESTIMATE NULL WITH ENSEMBLE 
% (TURNS OUT BEING AN OPTIMIZATION PROBLEM AND NOT ROBUST...)

% Generate AR noise time series similar to the ensemble PDSI reconstructions
%addpath('../ar-fit/')
% v[time,vars]
% [w,A,C,SBC,FPE,th]=ARFIT(v,pmin,pmax)
ntr=100;

di_rens_na=zeros(length(analyrs),size(xa_di_ens_na,2),ntr);
di_rens_sa=zeros(length(analyrs),size(xa_di_ens_na,2),ntr);
for i=1:ntr
   d0=xa_di_ens_na(analyrs,:);
   dn=d0+0.25*randn(size(d0)); % have to add noise because adding AR makes the correlations go up 
   [w,A,C,~,~,~]=arfit(dn,1,2); % explore up to 5, but it looks like all are 2nd order
   a0=arsim(w,A,C,[length(analyrs),1]);% generate base ensemble AR noise time series 
   di_rens_na(:,:,i)=bsxfun(@minus,a0,mean(a0,1))+mean(xa_di_ens_na,1); % add ensemble offset

   d0=xa_di_ens_sa(analyrs,:);
   dn=d0+0.25*randn(size(d0)); % have to add noise because adding AR makes the correlations go up 
   [w,A,C,~,~,~]=arfit(dn,1,2); % explore up to 5, but it looks like all are 2nd order
   a0=arsim(w,A,C,[length(analyrs),1]);% generate base ensemble AR noise time series 
   di_rens_sa(:,:,i)=bsxfun(@minus,a0,mean(a0,1))+mean(xa_di_ens_sa,1); % add ensemble offset
end



% Generate AR noise time series similar to the ensemble PDSI reconstructions
addpath('../ar-fit/')
% v[time,vars]
% [w,A,C,SBC,FPE,th]=ARFIT(v,pmin,pmax)
ntr=500;
di_rens_na=zeros(length(analyrs),ntr,size(xa_di_ens_na,2));
di_rens_sa=zeros(length(analyrs),ntr,size(xa_di_ens_na,2));
parfor i=1:size(xa_di_ens_na,2)
   [w,A,C,~,~,~]=arfit(xa_di_ens_na(analyrs,i),1,5); % explore up to 5, but it looks like all are 2nd order
   di_rens_na(:,:,i)=arsim(w,A,C,[length(analyrs),ntr]); % generate ntr realizations of length n
end
parfor i=1:size(xa_di_ens_na,2)
   [w,A,C,~,~,~]=arfit(xa_di_ens_sa(analyrs,i),1,5); % explore up to 5, but it looks like all are 2nd order
   di_rens_sa(:,:,i)=arsim(w,A,C,[length(analyrs),ntr]); % generate ntr realizations of length n
end

% Generate AR noise time series similar to the ensemble PDSI reconstructions
addpath('../ar-fit/')
% v[time,vars]
% [w,A,C,SBC,FPE,th]=ARFIT(v,pmin,pmax)
ntr=100;
di_rens_na=zeros(length(analyrs),size(xa_di_ens_na,2),ntr);
Cm_na=corr(xa_di_ens_na); % correlation matrix over analysis period
di_rens_sa=zeros(length(analyrs),size(xa_di_ens_na,2),ntr);
Cm_sa=corr(xa_di_ens_sa); % correlation matrix over analysis period
parfor i=1:ntr
   [w,A,C,~,~,~]=arfit(xa_di_ens_na(analyrs,randsample(size(xa_di_ens_na,2),1)),1,5); % explore up to 5, but it looks like all are 2nd order
   a0=arsim(w,A,C,[length(analyrs),size(xa_di_ens_na,2)]);% generate base ensemble AR noise time series 
   % Make a set of data with a specified correlation matrix
   di_rens_na(:,:,i)=squeeze(a0)*chol(Cm_na);
end
parfor i=1:ntr
   [w,A,C,~,~,~]=arfit(xa_di_ens_sa(analyrs,randsample(size(xa_di_ens_sa,2),1)),1,5); % explore up to 5, but it looks like all are 2nd order
   a0=arsim(w,A,C,[length(analyrs),size(xa_di_ens_sa,2)]);% generate base ensemble AR noise time series 
   % Make a set of data with a specified correlation matrix
   di_rens_sa(:,:,i)=squeeze(a0)*chol(Cm_sa);
end



% Run the noise time series through the same algorithm to get drought probability
mdsi_na=zeros(length(analyrs),ntr);
mdsi_sa=zeros(length(analyrs),ntr);
sub_ens=100;
for k=1:ntr
   % Ensemble estimate of drought periods
   dsens_na=struct();
   dsens_sa=struct();
   for i=1:sub_ens
      dsens_na.ds{i}=mdrought(squeeze(di_rens_na(:,i,k)),1,'y');
      dsens_sa.ds{i}=mdrought(squeeze(di_rens_sa(:,i,k)),1,'y');
   end

   mds=zeros(length(analyrs),1);
   % Find which years there is agreement
   for i=1:sub_ens
      for j=1:length(dsens_na.ds{i})
	 m0=zeros(length(analyrs),1);
	 m0(dsens_na.ds{i}(j,1):dsens_na.ds{i}(j,2))=1;
	 mds=mds+m0; % add up droughts
      end
   end
   mdsi_na(:,k)=mds;

   mds=zeros(length(analyrs),1);
   % Find which years there is agreement
   for i=1:sub_ens
      for j=1:length(dsens_sa.ds{i})
	 m0=zeros(length(analyrs),1);
	 m0(dsens_sa.ds{i}(j,1):dsens_sa.ds{i}(j,2))=1;
	 mds=mds+m0; % add up droughts
      end
   end
   mdsi_sa(:,k)=mds;
end


% Compare the null with the actual: threshold >0.5 or >0.66? 
%===============
%   JOINT PDF 
%===============

mna1=mdsi_na(:,1:100);
msa1=mdsi_sa(:,1:100);

X1=cat(2,mna1(:),msa1(:));

%bwn=[0.0422,0.0097];
%[f1,x1] = ksdensity(X1,'bandwidth',bwn);
[f1,x1] = ksdensity(X1);

xl=[0,100];yl=[0,100];
%[xq,yq] = meshgrid(linspace(-4.75,4.75),linspace(-1,1));
[xq,yq] = meshgrid(linspace(xl(1),xl(2),500),linspace(yl(1),yl(2),500));
z1 = griddata(x1(:,1),x1(:,2),f1,xq,yq);
%z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);

%clb.Label.HorizontalAlignment='right';
%[c2,h2]=contour(xq,yq,z2,'showtext','on','linecolor',clrs.red,'linestyle','-.');
%[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.red,'linestyle','-.');
%clabel(c2,h2,'Color',clrs.red)



figure
hold on
[c1,h1]=contourf(xq,yq,z1);%,0.1:0.2:1.5,'showtext','on');
colormap(flip(gray))
%clb=colorbar;
%clb.Label.String='Probability Density';
%clb.Label.FontSize=16;


% PLOT EACH MEGADROUGHT PDF INDIVIDUALLY
X2=cat(2,mdagr_na(:),mdagr_sa(:));
[f2,x2] = ksdensity(X2);
z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);

%[c2,h2]=contour(xq,yq,z2,0.1:0.2:1.5,'showtext','on','linecolor',clrs.red,'linestyle','-','linewidth',1.25);
[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.blue,'linestyle','-','linewidth',1.25);
clabel(c2,h2,'Color',clrs.blue)


scatter(mdagr_na(:),mdagr_sa(:),'jitter','on','jitterAmount',0.25,'MarkerEdgeColor',clrs.navy)

set(gca,'fontsize',16)
xlabel('NASW drought prob.','fontsize',18)
ylabel('Pampas drought prob.','fontsize',18)
box on

%print(['/home/nsteiger/Dropbox/sa-phyda/sa_na_droughtprob_null.png'],'-dpng','-r300');


figure
plot(analyrs,mdagr_na,'o','color',clrs.red)
hold on
plot(analyrs,mdagr_sa,'o','color',clrs.blue)

set(gca,'fontsize',16)
xlim([analyrs(1) analyrs(end)])
xlabel('Year','fontsize',18)
ylabel('Percent agreement','fontsize',18)
box on
legend('NASW','Pampas')
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%print(['/home/nsteiger/Dropbox/sa-phyda/sa_na_drought_agree.png'],'-dpng','-r300');




return % Stop here

%-----------
% RUNLENGTH
%-----------

% Max allowable length
%ml=max(dst(:,2)-dst(:,1));
%ml=ceil(mean(dst(:,2)-dst(:,1)));
ml=ceil(mean(ds_sa(:,2)-ds_sa(:,1)));

% Compute the run length
[nl,ns,ri] = runlength(dtm,ml);
ri(end+1)=length(dtm);


% Select out segments for binning and regression
xD=zeros(length(nl),2);
w0=zeros(length(nl),1);
y0=zeros(length(nl),1);

% Binning 
for i=1:length(nl)
   [~,~,irab] = intersect((ri(i)+analyrs4(1)-1):(ri(i+1)-1+analyrs4(1)-1),r_o:r_f); 
   w0(i)=length(irab); % bin size 

   % Mean of bins
   xD(i,1)=mean(xa_nino_ann(irab));
   % xD(i,2)=mean(xa_amo_mn(irab));
   xD(i,2)=mean(gmt_est_f(irab)); % myrs and r_o start in the same year so this indexing is OK
   %xD(i,3)=mean(gmt_est_sg(irab)); % myrs and r_o start in the same year so this indexing is OK
   %xD(i,3)=mean(gmt_est_vs(irab)); % myrs and r_o start in the same year so this indexing is OK
   %xD(i,4)=mean(gmt_est_v(irab));
   %xD(i,5)=mean(gmt_est_s(irab));


end

%disp('Add gaussian random variable')
%xD(:,3)=randn(size(dwdd,1),1);

% Test sensitivy of removing large volcanoes
%disp('Removing outlier volcanic segments')
xDnv=xD;
%volyr=[21,24,26,33];
%volyr=[24,33];
% Remove value that is more than three scaled median absolute deviations (MAD) away from the median
volyr=find(isoutlier(xD(:,2))==1); % i=[24,25,33,50], which are two MD and two non-MD 
%volyr=[];

xDnv(volyr,:)=NaN; % remove all volcanic year info
%xDnv(volyr,3)=NaN; % remove just forcing values
%xDnv(volyr,4)=NaN;


% Standardize data to compare regression coefficients
%Xd=zscore(xD);
X0 = bsxfun(@minus,xD,nanmean(xD,1));
Xd = bsxfun(@rdivide,X0,nanstd(X0,[],1));
X0 = bsxfun(@minus,xDnv,nanmean(xDnv,1));
Xdnv = bsxfun(@rdivide,X0,nanstd(X0,[],1));
% don't standardize (for getting actual coefficient values)
%Xd=xD;
%Xdnv=xDnv;

% Normalize weights?
wD=w0;
%wD=w0/max(w0);
%wD=w0/length(analyrs4)*100;
%wD=ones(size(w0));


% Make the Y variable
Yd=categorical(ns(:));
%error('Specific matlab function does not allow both weighted regression and mulitnomial regression...')
% PERHAPS USE SEPERATE LOGISTIC REGRESSIONS?

%figure;autocorr(yD);title('Is the drought time series autocorrelated?')




%--------------------
%    Regression
%--------------------

% Fit all varieties of models: What can explain PDSI in the American southwest?
% USE SAME AS ABOVE....
%xm={1,2,3,4,5,[1,2],[1,3],[1,4],[1,5],[1,2,3],[1,2,4],[1,2,5]};
%c0={'Nino3.4','AMO','Forc','Vol','Sol','N+A','N+F','N+V','N+S','N+A+F','N+A+V','N+A+S'};
%xm={1,2,3,[1,2],[1,3],[2,3],[1,2,3],4};
%c0={'Nino3.4','AMO','Forc','N+A','N+F','A+F','N+A+F','R'};
%xm={1,2,[1,2],3};
%c0={'Nino3.4','Forc','N+F','R'};
xm={1,2,[1,2]};
c0={'Nino3.4','Forc','N+F'};
c=categorical(c0,c0);

% Jackknife sampling for regression coefficients
%js=999;
js=300;

aic_j=zeros(length(xm),js,1);
bic_j=zeros(length(xm),js,1);
b_j=NaN(length(xm),size(xD,2)+1,js);
%b_j=NaN(length(xm),5,js); % adding noise var
%p_j=NaN(length(xm),4,js);

aic_jnv=zeros(length(xm),js,1);
bic_jnv=zeros(length(xm),js,1);
b_jnv=NaN(length(xm),size(xD,2)+1,js);
%b_jnv=NaN(length(xm),5,js); % adding noise var
%p_jnv=NaN(length(xm),4,js);

% Weighted logistic regression
for i=1:js

   % Jacknife sample
   %rsp=randsample(length(Yd),length(Yd)-1);
   % Bootstrap
   rsp=randsample(length(Yd),length(Yd),true);

   % Random variable
   %Xd(:,3)=zscore(randn(length(nl),1));
   %Xdnv(:,3)=zscore(randn(length(nl),1));

   for ii=1:length(xm)
      % Regression
      mdl = fitglm(Xd(rsp,xm{ii}),Yd(rsp),'Distribution','binomial','Weights',wD(rsp));
      % Grab stats
      b_j(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.Estimate;
      %p_j(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.pValue;
      aic_j(ii,i)=mdl.ModelCriterion.AICc;
      bic_j(ii,i)=mdl.ModelCriterion.BIC;

      % Regression WITH VOLCANOES REMOVED
      mdl = fitglm(Xdnv(rsp,xm{ii}),Yd(rsp),'Distribution','binomial','Weights',wD(rsp));
      % Grab stats
      b_jnv(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.Estimate;
      %p_jnv(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.pValue;
      aic_jnv(ii,i)=mdl.ModelCriterion.AICc;
      bic_jnv(ii,i)=mdl.ModelCriterion.BIC;
   end
end


clear b1
b1(:,1)=squeeze(b_j(1,2,:));
b1(:,2)=squeeze(b_j(2,3,:));
%b1(:,3)=squeeze(b_j(3,4,:));
%b1(:,4)=squeeze(b_j(4,5,:));
%b1(:,5)=squeeze(b_j(5,6,:));

clear b01
b01(:,1)=squeeze(b_j(length(xm)-1,2,:));
b01(:,2)=squeeze(b_j(length(xm)-1,3,:));
%b01(:,3)=squeeze(b_j(7,4,:));

clear b2
b2(:,1)=squeeze(b_jnv(1,2,:));
b2(:,2)=squeeze(b_jnv(2,3,:));
%b2(:,3)=squeeze(b_jnv(3,4,:));

clear b02
b02(:,1)=squeeze(b_jnv(length(xm)-1,2,:));
b02(:,2)=squeeze(b_jnv(length(xm)-1,3,:));
%b02(:,3)=squeeze(b_jnv(7,4,:));

%clear br
%br(:,1)=squeeze(b_j(length(xm),5,:));
%br(:,2)=squeeze(b_jnv(length(xm),5,:));
%br(:,1)=squeeze(b_j(length(xm),4,:));
%br(:,2)=squeeze(b_jnv(length(xm),4,:));


%clear p1
%p1(:,1)=squeeze(p_j(1,2,:));
%p1(:,2)=squeeze(p_j(2,3,:));
%p1(:,3)=squeeze(p_j(3,4,:));
%p1(:,4)=squeeze(p_j(4,5,:));
%p1(:,5)=squeeze(p_j(5,6,:));

% Model coefficients

% Restructure the data
bg(:,1)=b1(:,1);
bg(:,2)=b2(:,1);
bg(:,3)=b1(:,2);
bg(:,4)=b2(:,2);
%bg(:,5)=b1(:,3);
%bg(:,6)=b2(:,3);

% New categories
%c01={'Nino3.4-All','Nino3.4-NLV','AMO-All','AMO-NLV','Forc-All','Forc-NLV'};
c01={'Nino3.4-All','Nino3.4-NLV','Forc-All','Forc-NLV'};
c1=categorical(c01,c01);

% Color matrix
%cm=cat(1,clrs.green,clrs.green,clrs.blue,clrs.blue,clrs.red,clrs.red);
cm=cat(1,clrs.green,clrs.green,clrs.red,clrs.red);

figure
hold on
hl=boxplot(bg,c1,'Colors',cm,'Symbol','');
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
ax.YLim=[-4 2];
line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
ax.YGrid = 'on';
%ax.GridLineStyle = '--';
% Add percentage of values below:above 0
for j=1:size(bg,2);
   bpc=length(find(bg(:,j)<0))/length(bg(:,j))*100;
   text(j,ax.YLim(2)-0.5,[sprintf('%0.0f',bpc) ':' sprintf('%0.0f',100-bpc)],'HorizontalAlignment','center','fontsize',14)
end
ylabel('\beta_1','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 5 4]; % [left bottom width height]

%print(['figs/sa_logress_betas_bootstrap_V_NV.png'],'-dpng','-r300');


%figure
%boxplot(br)



% Best model via AIC and BIC?
[~,m_aic]=min(aic_j); % minimum of the AIC
[~,m_bic]=min(bic_j); % minimum of the AIC
hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
hcab0=cat(1,hc_a,hc_b);

% NO VOLCANOES
[~,m_aic]=min(aic_jnv); % minimum of the AIC
[~,m_bic]=min(bic_jnv); % minimum of the AIC
hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
hcab1=cat(1,hc_a,hc_b);

hcab=cat(1,hcab0,hcab1);

figure
b=bar(c,hcab','facecolor','flat','edgecolor','flat');
b(1).CData = clrs.white;
b(1).EdgeColor = clrs.grey;
b(1).LineWidth = 2;
b(2).CData = clrs.white;
b(2).EdgeColor = clrs.black;
b(2).LineWidth = 2;
b(3).CData = clrs.grey;
b(4).CData = clrs.black;
legend('AIC-All','BIC-All','AIC-NLV','BIC-NLV','location','northwest')
ylabel('Best model (%)','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]

%print(['figs/sa_logress_AICBIC_bootstrap_V_NV.png'],'-dpng','-r300');














% %---------------------------------------------
% % HOW WELL DO FORCINGS EXPLAIN NINO AND AMO?
% 
% 
% %xm={1,2,3};
% %c0={'Forc','Vol','Sol'};
% c0={'Forc'};
% c=categorical(c0,c0);
% 
% % Averaging length
% %al=10;
% 
% bs=200;
% %al=[8 9 10 11 12];
% al=[9 10 11];
% rsqrd_n=zeros(bs,length(al));
% rsqrd_a=zeros(bs,length(al));
% beta_n=NaN(2,bs,length(al));
% beta_a=NaN(2,bs,length(al));
% 
% k=1;
% for j=al 
%    % estimate NASW PDSI
%    clear X4 y4n y4a
%    y4n0=detrend(ann2xt(xa_nino_ann(ir4),j));
%    y4a0=detrend(ann2xt(xa_amo_mn(ir4),j));
%    X4=detrend(ann2xt(gmt_est_f(ig4),j));
%    %X4(:,2)=detrend(ann2xt(gmt_est_v(ig4),j)); 
%    %X4(:,3)=detrend(ann2xt(gmt_est_s(ig4),j)); 
%    % Standardize data to compare regression coefficients
%    y4n=(y4n0-nanmean(y4n0))./nanstd(y4n0);
%    y4a=(y4a0-nanmean(y4a0))./nanstd(y4a0);
%    X0 = bsxfun(@minus,X4,nanmean(X4,1));
%    X4 = bsxfun(@rdivide,X0,nanstd(X0,[],1));
% 
% 
%    for ii=1:bs 
%       % Bootstrap sample
%       rsp=randsample(length(y4n),length(y4n),true);
%       % Jacknife sample
% %         rsp=randsample(length(y4),length(y4)-1);
% 
%       % ESTIMATE AMO
%       mdl=fitlm(X4(rsp),y4a(rsp));
%       rsqrd_a(ii,k)=mdl.Rsquared.Adjusted;
%       beta_a(:,ii,k)=mdl.Coefficients.Estimate;
% 
%       % ESTIMATE NINO
%       mdl=fitlm(X4(rsp,xm{i}),y4n(rsp));
%       rsqrd_n(ii,k)=mdl.Rsquared.Adjusted;
%       beta_n(:,ii,k)=mdl.Coefficients.Estimate;
% 
%    end
%    k=k+1;
% end
% 
% r2a=reshape(rsqrd_a,bs*length(al));
% r2n=reshape(rsqrd_n,bs*length(al));
% b0a=reshape(beta_a,4,bs*length(al));
% b0n=reshape(beta_n,4,bs*length(al));
% 
% % Just take values for individual models
% clear b1a b1n
% b1a(:,1)=squeeze(b0a(1,2,:));
% %b1a(:,2)=squeeze(b0a(2,3,:));
% %b1a(:,3)=squeeze(b0a(3,4,:));
% 
% b1n(:,1)=squeeze(b0n(1,2,:));
% %b1n(:,2)=squeeze(b0n(2,3,:));
% %b1n(:,3)=squeeze(b0n(3,4,:));
% 
% 
% 
% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2a',c)
% plot(1:length(c),zeros(length(c),1),'--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% title('AMO')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcamo_r2_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcamo_r2_bootstrap.png'],'-dpng','-r300');
% 
% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2n',c)
% plot(1:length(c),zeros(length(c),1),'--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% title('Nino3.4')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcnino_r2_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcnino_r2_bootstrap.png'],'-dpng','-r300');
% 
% % Model coefficients
% figure
% hold on
% boxplot(b1a,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% title('AMO')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcamo_betas_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcamo_betas_bootstrap.png'],'-dpng','-r300');
% 
% % Model coefficients
% figure
% hold on
% boxplot(b1n,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% title('Nino3.4')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcnino_betas_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcnino_betas_bootstrap.png'],'-dpng','-r300');


%---------------------------------------
% DIFFERENT TIME PERIODS FOR REGRESSION

% % estimate NASW PDSI
% y1=detrend(ann2xt(Xa_di(ir1),al));
% X1(:,1)=detrend(ann2xt(xa_nino_ann(ir1),al));
% X1(:,2)=detrend(ann2xt(xa_amo_mn(ir1),al));
% X1(:,3)=detrend(ann2xt(gmt_est_f(ig1),al)); % full forcing
% X1(:,4)=detrend(ann2xt(gmt_est_v(ig1),al)); % full forcing
% X1(:,5)=detrend(ann2xt(gmt_est_s(ig1),al)); % full forcing
% 
% % estimate NASW PDSI
% y2=detrend(ann2xt(Xa_di(ir2),al));
% X2(:,1)=detrend(ann2xt(xa_nino_ann(ir2),al));
% X2(:,2)=detrend(ann2xt(xa_amo_mn(ir2),al));
% X2(:,3)=detrend(ann2xt(gmt_est_f(ig2),al)); % full forcing
% X2(:,4)=detrend(ann2xt(gmt_est_v(ig2),al)); % full forcing
% X2(:,5)=detrend(ann2xt(gmt_est_s(ig2),al)); % full forcing
% 
% % estimate NASW PDSI
% y4=detrend(ann2xt(Xa_di(ir4),al));
% X4(:,1)=detrend(ann2xt(xa_nino_ann(ir4),al));
% X4(:,2)=detrend(ann2xt(xa_amo_mn(ir4),al));
% X4(:,3)=detrend(ann2xt(gmt_est_f(ig4),al)); % full forcing
% X4(:,4)=detrend(ann2xt(gmt_est_v(ig4),al)); % full forcing
% X4(:,5)=detrend(ann2xt(gmt_est_s(ig4),al)); % full forcing
% 
% % y4o=ann2xt(Xa_di(ir4),al);
% % X4o(:,1)=ann2xt(xa_nino_ann(ir4),al);
% % X4o(:,2)=ann2xt(xa_amo_mn(ir4),al);
% % X4o(:,3)=ann2xt(gmt_est_f(ig4),al); % full forcing
% % X4o(:,4)=ann2xt(gmt_est_v(ig4),al); % full forcing
% % X4o(:,5)=ann2xt(gmt_est_s(ig4),al); % full forcing
% 
% rsqrd_avg1=zeros(length(xm),1);
% rsqrd_avg2=zeros(length(xm),1);
% rsqrd=zeros(length(xm),1);
% %rsqrdo=zeros(length(xm),1);
% aic_bs_avg1=zeros(length(xm),1);
% aic_bs_avg2=zeros(length(xm),1);
% aic_bs=zeros(length(xm),1);
% dwp_avg1=zeros(length(xm),1);
% dwp_avg2=zeros(length(xm),1);
% dwp_avg4=zeros(length(xm),1);
% %dwp_avg4o=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X1(:,xm{i}),y1);
%    rsqrd_avg1(i)=mdl.Rsquared.Adjusted;
%    dwp_avg1(i)=dwtest(mdl);
%    aic_bs_avg1(i)=mdl.ModelCriterion.AICc;   
% 
%    mdl=fitlm(X2(:,xm{i}),y2);
%    rsqrd_avg2(i)=mdl.Rsquared.Adjusted;
%    dwp_avg2(i)=dwtest(mdl);
%    aic_bs_avg2(i)=mdl.ModelCriterion.AICc;   
% 
%    mdl=fitlm(X4(:,xm{i}),y4);
%    rsqrd(i)=mdl.Rsquared.Adjusted;
%    dwp_avg4(i)=dwtest(mdl);
%    aic_bs(i)=mdl.ModelCriterion.AICc;   
% 
%    %mdl=fitlm(X4o(:,xm{i}),y4o);
%    %rsqrdo(i)=mdl.Rsquared.Adjusted;
%    %dwp_avg4o(i)=dwtest(mdl);
% 
% end
% 
% % horizontal bar plot
% figure
% % ax1=subplot(3,1,1);
% % barh(c,rsqrd_avg1)
% % title('800-1600')
% % ax1.XGrid='on';
% % xlabel(['Adjusted r^2'])
% % 
% % ax2=subplot(3,1,2);
% % barh(c,rsqrd_avg2)
% % title('1601-1925')
% % ax2.XGrid='on';
% % xlabel(['Adjusted r^2'])
% 
% %ax3=subplot(3,1,3);
% barh(c,rsqrd)
% title('800-1925')
% %ax3.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% %linkaxes([ax3,ax2,ax1],'x');
% 
% %set(gca,'fontsize',16)
% 
% %fig = gcf;fig.PaperUnits = 'inches';
% %fig.PaperPosition = [0 2 5 7.5]; % [left bottom width height]
% 
% %print(['figs/linearfit_comparison_r2_10yr.png'],'-dpng','-r300');
% 




%------------------------------------
% LINEAR REGRESSION WITH ANNUAL DATA
%------------------------------------


% % estimate NASW PDSI
% y1=detrend(Xa_di(ir1));
% X1(:,1)=detrend(xa_nino_ann(ir1));
% X1(:,2)=detrend(xa_amo_mn(ir1));
% X1(:,3)=detrend(gmt_est_f(ig1)); % full forcing
% X1(:,4)=detrend(gmt_est_v(ig1)); % full forcing
% X1(:,5)=detrend(gmt_est_s(ig1)); % full forcing
% 
% % estimate NASW PDSI
% y2=detrend(Xa_di(ir2));
% X2(:,1)=detrend(xa_nino_ann(ir2));
% X2(:,2)=detrend(xa_amo_mn(ir2));
% X2(:,3)=detrend(gmt_est_f(ig2)); % full forcing
% X2(:,4)=detrend(gmt_est_v(ig2)); % full forcing
% X2(:,5)=detrend(gmt_est_s(ig2)); % full forcing
% 
% % estimate NASW PDSI
% y4=detrend(Xa_di(ir4));
% X4(:,1)=detrend(xa_nino_ann(ir4));
% X4(:,2)=detrend(xa_amo_mn(ir4));
% X4(:,3)=detrend(gmt_est_f(ig4)); % full forcing
% X4(:,4)=detrend(gmt_est_v(ig4)); % full forcing
% X4(:,5)=detrend(gmt_est_s(ig4)); % full forcing
% 
% % ACCESS LINEAR MODEL FIT OUTPUT:
% % http://www.mathworks.com/help/stats/linearmodel-class.html
% 
% rsqrd_ann1=zeros(length(xm),1);
% rsqrd_ann2=zeros(length(xm),1);
% rsqrd_ann4=zeros(length(xm),1);
% aic_bs_ann1=zeros(length(xm),1);
% aic_bs_ann2=zeros(length(xm),1);
% aic_bs_ann4=zeros(length(xm),1);
% dwp_ann1=zeros(length(xm),1);
% dwp_ann2=zeros(length(xm),1);
% dwp_ann4=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X1(:,xm{i}),y1);
%    rsqrd_ann1(i)=mdl.Rsquared.Adjusted;
%    aic_bs_ann1(i)=mdl.ModelCriterion.AICc;   
%    dwp_ann1(i)=dwtest(mdl);
% 
%    mdl=fitlm(X2(:,xm{i}),y2);
%    rsqrd_ann2(i)=mdl.Rsquared.Adjusted;
%    aic_bs_ann2(i)=mdl.ModelCriterion.AICc;   
%    dwp_ann2(i)=dwtest(mdl);
%    
%    mdl=fitlm(X4(:,xm{i}),y4);
%    rsqrd_ann4(i)=mdl.Rsquared.Adjusted;
%    aic_bs_ann4(i)=mdl.ModelCriterion.AICc;   
%    dwp_ann4(i)=dwtest(mdl);
% end
% 
% % horizontal bar plot
% figure
% ax1=subplot(3,1,1);
% barh(c,rsqrd_ann1)
% title('800-1600')
% ax1.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% ax2=subplot(3,1,2);
% barh(c,rsqrd_ann2)
% title('1601-1925')
% ax2.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% ax3=subplot(3,1,3);
% barh(c,rsqrd_ann4)
% title('800-1925')
% ax3.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% linkaxes([ax3,ax2,ax1],'x');
% 
% %set(gca,'fontsize',16)
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 7.5]; % [left bottom width height]
% 
% %print(['figs/linearfit_comparison_r2_ann.png'],'-dpng','-r300');
% 
% 
% % horizontal bar plot
% figure
% ax1=subplot(3,1,1);
% barh(c,1e-3.*aic_bs_ann1)
% title('800-1600')
% ax1.XGrid='on';
% xlabel(['AIC'])
% 
% ax2=subplot(3,1,2);
% barh(c,1e-3.*aic_bs_ann2)
% title('1601-1925')
% ax2.XGrid='on';
% xlabel(['AIC'])
% 
% ax3=subplot(3,1,3);
% barh(c,1e-3.*aic_bs_ann4)
% title('800-1925')
% ax3.XGrid='on';
% xlabel(['AIC'])
% 
% linkaxes([ax3,ax2,ax1],'x');
% 
% %set(gca,'fontsize',16)
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 7.5]; % [left bottom width height]
% 
% %print(['figs/linearfit_comparison_r2_ann.png'],'-dpng','-r300');


% is proportional to the probability that the ith model minimizes the (estimated) information loss
% exp((AICmin − AICi)/2)


% % Can any forcings explain Nino or AMO?
% xm={3,4,5};
% rsqrd_n=zeros(length(xm),1);
% rsqrd_a=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X(:,xm{i}),y2);rsqrd_n(i)=mdl.Rsquared.Adjusted;
%    mdl=fitlm(X(:,xm{i}),y3);rsqrd_a(i)=mdl.Rsquared.Adjusted;
% end
% 
% 
% % horizontal bar plot
% c0={'Forc','Vol','Sol'};
% c=categorical(c0,c0);
% figure
% barh(c,rsqrd_n)
% title('Forcing explain Nino3.4?')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_nino_forcings_r2.png'],'-dpng','-r300');
% 
% figure
% barh(c,rsqrd_a)
% title('Forcing explain AMO?')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_amo_forcings_r2.png'],'-dpng','-r300');










% % Averaged PDSI
% y0=Xa_di(ir4)-mean(Xa_di(ir2));
% 
% % Break up time series by 'bl' year segments; find how many are drought years and bin accordingly
% bl=10; % bin length; 21 = smallest bin length with no autocorrelation in the period 800-1600
% segs=fix(length(y0)/bl);
% fn=length(y0)-mod(length(y0),bl);
% 
% ynD=sum(reshape(y0(1:fn),bl,segs))./bl;
% 
% % Mean climate and forcing values in each bin
% X=zeros(segs,6);
% X0=xa_nino_ann(ir4)-mean(xa_nino_ann(ir2));
% X(1:segs,1)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,1)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=xa_amo_mn(ir4)-mean(xa_amo_mn(ir2));
% X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_f(ig4)-mean(gmt_est_f(ig2));
% X(1:segs,3)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,3)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_v(ig4)-mean(gmt_est_v(ig2));
% X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_s(ig4)-mean(gmt_est_s(ig2));
% X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_g(ig4)-mean(gmt_est_g(ig2));
% X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs));
% 
% 
% 
% % Fit all varieties of models: What can explain PDSI in the American southwest?
% xm={1,2,3,4,5,6,[1,2],[1,3],[1,4],[1,5],[1,6],[1,2,3],[1,2,4]};
% rsqrd=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X(:,xm{i}),ynD);rsqrd(i)=mdl.Rsquared.Adjusted;
% end
% 
% % horizontal bar plot
% c0={'Nino','AMO','Forc','Vol','Sol','GHG','N+A','N+F','N+V','N+S','N+G','N+A+F','N+A+V'};
% c=categorical(c0,c0);
% figure
% barh(c,rsqrd)
% grid on
% xlabel(['Adjusted r^2'],'fontsize',16)
% set(gca,'fontsize',14)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% title(['Linear regression, ' num2str(bl) ' yr binned data'])
% 
% 
% 
% % Can any forcings explain Nino or AMO?
% xm={3,4,5,6};
% rsqrd_n=zeros(length(xm),1);
% rsqrd_a=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X(:,xm{i}),X(:,1));rsqrd_n(i)=mdl.Rsquared.Adjusted;
%    mdl=fitlm(X(:,xm{i}),X(:,2));rsqrd_a(i)=mdl.Rsquared.Adjusted;
% end
% 
% 
% % horizontal bar plot
% c0={'Forc','Vol','Sol','GHG'};
% c=categorical(c0,c0);
% figure
% barh(c,rsqrd_n)
% title('Forcing explain Nino3.4? (avg)')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_nino_forcings_r2.png'],'-dpng','-r300');
% 
% figure
% barh(c,rsqrd_a)
% title('Forcing explain AMO? (avg)')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_amo_forcings_r2.png'],'-dpng','-r300');




% %------------------------------------------
% % LOGISTIC REGRESSION ANALYSIS
% %------------------------------------------
% 
% % Bin the time series and do a logistric regression 
% 
% %ds=svmd; % only use megadroughts
% %ds=ds4; % use all decadal droughts
% ird_ann=[];
% for i=1:size(ds,1)
%   [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),r_o:r_f); 
%   % Collect all years of drought
%   ird_ann=cat(1,ird_ann,irab);
% end
% 
% % Break up time series into single megadrought events and non-events
% 
% 
% y0=zeros(reconYrs,1);
% y0(ird_ann)=1;
% yd=y0(ir4);
% 
% % Break up time series by 'bl' year segments; find how many are drought years and bin accordingly
% %bl=10; % bin length; 21 = smallest bin length with no autocorrelation in the period 800-1600
% segs=fix(length(yd)/bl);
% fn=length(yd)-mod(length(yd),bl);
% ydr=reshape(yd(1:fn),bl,segs);
% ydi=find(sum(ydr)>=round(bl/2));
% yD=zeros(segs,1); % bins with not enough drought = 0
% yD(ydi)=1; % bins with more drought = 1
% 
% figure
% autocorr(yD)
% title('Is the drought time series autocorrelated?')
% 
% % Mean climate and forcing values in each bin
% X=zeros(segs,6);
% X0=xa_nino_ann(ir4)-mean(xa_nino_ann(ir2));
% X(1:segs,1)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,1)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=xa_amo_mn(ir4)-mean(xa_amo_mn(ir2));
% X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_f(ig4)-mean(gmt_est_f(ig2));
% X(1:segs,3)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,3)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_v(ig4)-mean(gmt_est_v(ig2));
% X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_s(ig4)-mean(gmt_est_s(ig2));
% X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_g(ig4)-mean(gmt_est_g(ig2));
% X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs));
% 
% 
% Y=categorical(yD);
% %disp('logistic reg: Only using the ENSO and AMO as predictors')
% %[B,dev,stats] = mnrfit(X(:,1:2),Y);
% [B,dev,stats] = mnrfit(X,Y);
% 
% %close all
% 
% figure
% plot(1:length(stats.p),stats.p,'o')
% hold on
% plot(1:length(stats.p),0.05*ones(7,1))
% title('Logistic regression p-values')
% xticklabels({'Int.','Nino','AMO','Forc','Vol','Sol','GHG'});
% 
% %print(['figs/logregss_p_bsues.png'],'-dpng','-r300');
% 
% figure
% plot(stats.resid,'.')
% title('Logistic regression residuals')



% % Explore the predicted probabilities using the range of ENSO and AMO
% 
% % Define a grid of all possible Nino and AMO values
% %ns=xa_nino_ann(ir4,:)-nino_off;as=xa_amo_ens(ir4,:)-amo_off;
% % Just pick the mean values over which the model is built
% %ns=mean(xa_nino_ann(ir4,:)-nino_off,2);as=mean(xa_amo_ens(ir4,:)-amo_off,2);
% ns=X(:,1);as=X(:,2);
% [mx,my] = meshgrid(min(ns(:)):0.01:max(ns(:)),min(as(:)):0.01:max(as(:)));
% Xm(:,1)=mx(:);
% Xm(:,2)=my(:);
% 
% % Get probability for random observations
% pihat = mnrval(B,Xm);
% [pihat,dlow,dhi] = mnrval(B,Xm,stats);% get uncertainties
% % NOTE THAT UPPER AND LOWER BOUNDS ARE THE SAME!!
% 
% 
% % Reshape to match the grid
% pM=reshape(pihat(:,2),size(mx));
% pL=reshape(dlow(:,2),size(mx));
% %pU=reshape(dhi(:,2),size(mx));
% 
% 
% % Alternative plot with both probability and uncertainty
% figure
% hold on
% [C1,h1]=contour(mx,my,pM,'fill','on');
% set(gca,'tickdir','both')
% colormap(flip(gray))
% c=colorbar;
% c.Label.String = 'Est. decadal drought prob.';
% c.Label.FontSize = 16;
% [Cc,h]=contour(mx,my,pL,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% clabel(Cc,h,'Color',clrs.red)
% set(gca,'fontsize',16)
% xlabel('Nino3.4 (10 yr anom. deg. C)','fontsize',18)
% ylabel('AMO (10 yr anom. deg. C)','fontsize',18)
% box on

%if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');end
%print(['./figs/logreg_decdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');










% % FORCING AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(analyrs3,gmt_est_f(ig3),'linewidth',2,'color',clrs.black);
% h2=plot(analyrs3,gmt_est_v(ig3),'linewidth',2,'color',clrs.red);
% h3=plot(analyrs3,gmt_est_s(ig3),'linewidth',2,'color',clrs.blue);
% h4=plot(analyrs3,gmt_est_g(ig3),'linewidth',2,'color',clrs.green);
% %plot(analyrs2,smooth(Xa_di(ir2),10,'lowess'),'linewidth',2,'color',clrs.red)
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.grey)
% set(gca,'fontsize',16)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Forced response (anom. ' char(176) 'C)'],'fontsize',18)
% box on
% 
% 
% ax=gca;
% yMin=ax.YLim(1)+0.5;yMax=ax.YLim(2)-0.5;
% %y=gmt_est_mn(ig3); % using 'l' from above
% %yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for megadroughts
% for i=1:length(ds)
%     p=patch(analyrs3(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],clrs.orange);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% % Add patches for all other decadal droughts
% % for i=10:size(ds,1)
% %     p=patch(analyrs3(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],clrs.olive);
% %     p.FaceAlpha=0.25;p.EdgeColor='none';
% % end
% 
% legend([h1,h2,h3,h4],{'All','Volc.','Solar','GHG'},'location','southwest')
% hold off
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% %print(['figs/forcings_alldroughts_800_2000.png'],'-dpng','-r300');
% 
% % FORCING, NINO, AMO, AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(mon_r(irm3_a:irm3_b),xa_nino_mn(irm3_a:irm3_b)-nino_off,'linewidth',1.5,'color',clrs.blue);
% h3=plot(analyrs3,xa_amo_mn(ir3)-amo_off,'linewidth',2,'color',clrs.green);
% h2=plot(analyrs3,gmt_est_f(ig3),'linewidth',2,'color',clrs.red);
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
% set(gca,'fontsize',12)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',13)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',13)
% box on
% 
% ax=gca;
% yMin=ax.YLim(1)+0.5;yMax=ax.YLim(2)-0.5;
% %y=gmt_est_mn(ig3); % using 'l' from above
% %yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for megadroughts
% for i=1:length(ds)
%     p=patch(analyrs3(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],clrs.orange);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% %legend([h1,h3,h2],{'Nino3.4','AMO','Forcing'},'location','southeast')
% hold off
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8 3]; % [left bottom width height]
% 
% %print(['figs/mdrought_nino_forcing_amo_800_2000.png'],'-dpng','-r300');

















% % Bivariate logistic regression with interaction
% %Xi(:,1)=mean(xa_nino_ann(ir4,:)-nino_off,2);
% %Xi(:,2)=mean(xa_amo_ens(ir4,:)-amo_off,2);
% %Xi(:,3)=Xi(:,1).*Xi(:,2);
% %[Bi,devi,statsi] = mnrfit(Xi,Y);
% 
% 
% 
% % Explore the predicted probabilities using the range of ENSO and AMO
% 
% % Define a grid of all possible Nino and AMO values
% %ns=xa_nino_ann(ir4,:)-nino_off;as=xa_amo_ens(ir4,:)-amo_off;
% % Just pick the mean values over which the model is built
% ns=mean(xa_nino_ann(ir4,:)-nino_off,2);as=mean(xa_amo_ens(ir4,:)-amo_off,2);
% [mx,my] = meshgrid(min(ns(:)):0.01:max(ns(:)),min(as(:)):0.01:max(as(:)));
% Xm(:,1)=mx(:);
% Xm(:,2)=my(:);
% 
% % Get probability for random observations
% pihat = mnrval(B,Xm);
% [pihat,dlow,dhi] = mnrval(B,Xm,stats);% get uncertainties
% % NOTE THAT UPPER AND LOWER BOUNDS ARE THE SAME!!
% 
% 
% % Reshape to match the grid
% pM=reshape(pihat(:,2),size(mx));
% pL=reshape(dlow(:,2),size(mx));
% %pU=reshape(dhi(:,2),size(mx));
% 
% % figure
% % [Cc,h]=contour(mx,my,pM,'ShowText','on','fill','on','linecolor','r');
% % colormap(flip(gray))
% % clabel(Cc,h,'FontSize',14,'Color','red')
% % title('Predicted megadrought probability','fontsize',20)
% % set(gca,'fontsize',16)
% % xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% % ylabel('AMO (anom. deg. C)','fontsize',18)
% % 
% % if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo.png'],'-dpng','-r300');end
% % 
% % 
% % figure
% % [Cc,h]=contour(mx,my,pL,'ShowText','on','linecolor','k');
% % clabel(Cc,h,'FontSize',14,'Color','black')
% % title('Probability 95% conf.','fontsize',20)
% % set(gca,'fontsize',16)
% % xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% % ylabel('AMO (anom. deg. C)','fontsize',18)
% % 
% % if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_unc.png'],'-dpng','-r300');end
% 
% 
% % Alternative plot with both probability and uncertainty
% figure
% hold on
% [C1,h1]=contour(mx,my,pM,'fill','on');
% set(gca,'tickdir','both')
% colormap(flip(gray))
% c=colorbar;
% c.Label.String = 'Est. megadrought probability';
% c.Label.FontSize = 16;
% [Cc,h]=contour(mx,my,pL,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% clabel(Cc,h,'Color',clrs.red)
% set(gca,'fontsize',16)
% xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% ylabel('AMO (anom. deg. C)','fontsize',18)
% box on
% 
% %if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');end
% 
















% % AMO offset
% amo_off=mean(mean(xa_amo_ens(ir2,:),2));
% xa_amo_mn=mean(xa_amo_ens,2);
% 
% % FORCING AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(analyrs3,xa_amo_mn(ir3)-amo_off,'linewidth',2,'color',[0 0 0]);
% plot(analyrs3,gmt_est_mn(ig3),'linewidth',2,'color',clrs.red)
% %plot(analyrs2,smooth(Xa_di(ir2),10,'lowess'),'linewidth',2,'color',clrs.red)
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
% set(gca,'fontsize',16)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',18)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% y=Xa_di(ir3); % using 'l' from above
% yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for top five times from truth
% for i=1:9
%     p=patch(analyrs3(1)+[ds4(i,1) ds4(i,2) ds4(i,2) ds4(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% %legend({'NADA','NADA-11'},'location','southeast')
% hold off
% 
% %print(['figs/gmtest_drght_AMO_cmp.png'],'-dpng','-r300');
% 
% 
% 
% 
% 
% nino_off=mean(mean(xa_nino_ens(irm2_a:irm2_b,:),2));
% 
% 
% % FORCING AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(mon_r(irm3_a:irm3_b),xa_nino_mn(irm3_a:irm3_b)-nino_off,'linewidth',2,'color',[0 0 0]);
% plot(analyrs3,gmt_est_mn(ig3),'linewidth',2,'color',clrs.red)
% %plot(analyrs2,smooth(Xa_di(ir2),10,'lowess'),'linewidth',2,'color',clrs.red)
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
% set(gca,'fontsize',16)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',18)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% y=Xa_di(ir3); % using 'l' from above
% yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for top five times from truth
% for i=1:9
%     p=patch(analyrs3(1)+[ds4(i,1) ds4(i,2) ds4(i,2) ds4(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% %legend({'NADA','NADA-11'},'location','southeast')
% hold off
% 
% %print(['figs/gmtest_drght_Nino_cmp.png'],'-dpng','-r300');
% 
% 









