


%===========================================================================
% 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;
%analyrs3=800: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;



%=====================================
% 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;


% %--------------------------------------
% % Calculate CPS reconstruction
% %--------------------------------------
% 
% % DIDN'T REALLY SEEM TO WORK OUT WITHOUT VARIANCE MATCHING THE WHOLE RECON THROUGH TIME...
% 
% %##########################################################################
% % create proxy data just as done in PaiCo reconstruction
% efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
% pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
% load([pth,efl],'proxy','xmeta')
% clat=xmeta{1}.lat;
% clon=xmeta{1}.lon;
% % Find proxies over southern South America
% p_lat=zeros(length(proxy),1);
% p_lon=zeros(length(proxy),1);
% for i=1:length(proxy)
%    p_lat(i)=proxy{i}.location(1);
%    p_lon(i)=proxy{i}.location(2);
% end
% % Find which proxies are over south america land areas
% %pi=find(p_lat<0 & p_lat>-57 & p_lon>270 & p_lon<320);
% % Find which proxies are in the north-central Chile region
% pi=find(p_lat<-26 & p_lat>-45 & p_lon>270 & p_lon<293);
% % Which proxies to use?
% prxytyp='p'; % PDSI-like trees
% %prxytyp='a'; % all proxies
% j=1;
% data=struct();
% for i=1:length(pi)
%    if strcmp(prxytyp,'p')
%       if strcmp(proxy{pi(i)}.ptype,'Tree Rings')
% 	 if strcmp(proxy{pi(i)}.psm_mthd,'p')
% 	 %if strcmp(proxy{pi(i)}.psm_mthd,'t')
% 	    % Proxy locations
% 	    data.proxy{j}.locations=[proxy{pi(i)}.location(1) proxy{pi(i)}.location(2)];
%             % Proxy data, make sure the sign of the proxy data lines up with positive correlation with PDSI
% 	    %data.proxy{j}.data=(proxy{pi(i)}.data)*sign(proxy{pi(i)}.pcr);
%             data.proxy{j}.data=proxy{pi(i)}.data; % don't manually adjust sign
% 	    % Proxy years
% 	    data.proxy{j}.times=proxy{pi(i)}.age_yng';
%             % Proxy correlation with obs
% 	    data.proxy{j}.corr=proxy{pi(i)}.pcr;
% 	    j=j+1;
% 	 end
%       end
%    elseif strcmp(prxytyp,'a')
%       % Proxy locations
%       data.proxy{j}.locations=[proxy{pi(i)}.location(1) proxy{pi(i)}.location(2)];
%       % Proxy data
%       data.proxy{j}.data=proxy{pi(i)}.data;
%       % Proxy years
%       data.proxy{j}.times=proxy{pi(i)}.age_yng';
%       % Proxy correlation with obs
%       data.proxy{j}.corr=proxy{pi(i)}.pcr;
%       j=j+1;
%    end
% end
% %##########################################################################
% 
% % Get instrumental data
% pth2='/home/nsteiger/MATLAB/PaiCo/phyda/';
% fln2='phyda_data_chile_pdsilikePrxs_nosignadj.mat';
% d=load([pth2,fln2],'data');
% % vandershier PDSI DJF over Chile
% instd=d.data.instrumental.data;
% instt=d.data.instrumental.times;
% 
% % Align positively with obs and compute correlation for weighting
% pd=[];
% pwt=zeros(length(data.proxy),1);
% for i=1:length(data.proxy)
%    [~,i_i,i_p] = intersect(instt,data.proxy{i}.times); 
%    pd0=data.proxy{i}.data;
%    pc0=nancorr(instd(i_i),pd0(i_p));
%    pd(i,:)=pd0*sign(pc0); % orient data according to correlation
%    pwt(i)=pc0^2; % weight is related to correlation
% end
% 
% disp('Disregard weights...')
% pwt=ones(size(pwt));
% 
% % Bootstrap composite reconstruction
% enz=100;
% pwmne=zeros(size(pd,2),enz);
% % overlap period between instrumental and proxy data
% [~,i_i,i_p] = intersect(instt,data.proxy{1}.times); 
% inst_ovrlp=instd(i_i);
% inst_std=nanstd(inst_ovrlp);
% inst_mn=nanmean(inst_ovrlp);
% for i=1:enz
%    % Random sample of proxies
%    ri=randsample(size(pd,1),round(size(pd,1)-1));
%    %----------------------------------------------------
%    % Compute the composite
%    rcn0=nansum(bsxfun(@times,pwt(ri),pd(ri,:)),1)/sum(pwt(ri));
%    rcn0(rcn0==0)=NaN;
%    %----------------------------------------------------
%    % Method-of-moments aka Mean-variance matching
%    rcn_ovrlp=rcn0(i_p);
%    rcn_ovrlp_mn=nanmean(rcn_ovrlp);
%    mul=inst_std/nanstd(rcn_ovrlp);
%    pwmne(:,i)=(rcn0-rcn_ovrlp_mn)*mul+inst_mn;
% 
% end
% 
% 
% figure
% plot(pwmne)
% title('Composite + scale over obs period (no weights)')
% xlabel('Year','fontsize',18)
% ylabel(['PDSI'],'fontsize',18)
% xlim([analyrs3(1) analyrs3(end)])
% 
% %print(['/home/nsteiger/Dropbox/sa-phyda/phyda_cps_chile.png'],'-dpng','-r300');
% print(['/home/nsteiger/Dropbox/sa-phyda/phyda_cpsNoWght_chile.png'],'-dpng','-r300');



%--------------------------------------
% Load Paico-based PDSI reconstruction
%--------------------------------------

%fln='chile_phyda_r.1000_1995_prxyN.47_27-Sep-2019_14:08:51.mat';
%fln='chile_phyda_r.1000_1995_prxyN.47_27-Sep-2019_14:25:00.mat';
fln='chile_phyda_r.800_1995_prxyN.47_03-Oct-2019_16:36:59.mat';% all PDSI-like proxies in Chile box
%fln='chile_phyda_r.800_1995_prxyN.77_07-Jan-2020_12:15:30.mat';% all proxies in Chile box
pth='/home/nsteiger/MATLAB/PaiCo/phyda/results/';
load([pth,fln])

%figure;plot(result.times,result.resample.signals')
%hold on;plot(result.times,mean(result.resample.signals),'k','linewidth',2)


% vandershier PDSI DJF over Chile --> get from PaiCo recon
vdsp=data.instrumental.data;
vdst=data.instrumental.times;



load distinct_colors_20.mat

phydatm=xa_di_sa(ir3);
% remove mean?
%phydatm=phydatm-mean(phydatm);
phyoff=mean(phydatm);

[~,~,ip3] = intersect(analyrs3,result.times); % indices of recon segment
paicotm=mean(result.resample.signals(:,ip3));
paicotm=paicotm-mean(paicotm)+mean(phydatm);

cphy=clrs.maroon;
cpai=clrs.blue;

% PHYDA VS PAICO
figure
hold on
%h1=plot(analyrs3,xa_di_na(ir3),'linewidth',1,'color',[clrs.blue 0.5]); % extra number on the end adds transparency
h3=plot(result.times(ip3),paicotm-phyoff,'linewidth',1,'color',[cpai 0.35]);
h5=plot(result.times(ip3),smooth(paicotm,11)-phyoff,'color',cpai,'linewidth',2.1);

h2=plot(analyrs3,phydatm-phyoff,'linewidth',1,'color',[cphy 0.35]);
h4=plot(analyrs3,smooth(phydatm,11)-phyoff,'linewidth',2.1,'color',cphy);
%plot(analyrs3,smooth(DA_di(id3),11),'linewidth',2,'color',clrs.red)
%h3=plot(analyrs3,smooth(xa_di_na(ir3),11),'linewidth',2,'color',clrs.blue);

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])
ylim([-6 6])

legend([h4,h5],{'SASW PHYDA','SASW PaiCo'},'location','northwest')
hold off


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


disp('Skill of phyda and paico recons over chile')

% validation years
vyrs=1903:1995; % can't go beyond 1995 b/c PaiCo doesn't go that recent b/c of data lack
[~,~,pvi] = intersect(vyrs,result.times); % years in paico
[~,~,dvi] = intersect(vyrs,r_o:r_f); % years in phyda
[~,~,ovi] = intersect(vyrs,vdst); % years in obs

% common mean
phyda_v=xa_di_sa(dvi)-mean(xa_di_sa(dvi));
paico_v=mean(result.resample.signals(:,pvi))-mean(mean(result.resample.signals(:,pvi)));
vobs_v=vdsp(ovi)-mean(vdsp(ovi));

% compute correlation
disp(['corr(phyda,vandrsr) = ' num2str(nancorr(phyda_v(:),vobs_v(:)))])
disp(['corr(paico,vandrsr) = ' num2str(nancorr(paico_v(:),vobs_v(:)))])

% compute mae
disp(['mae(phyda,vandrsr) = ' num2str(mae(vobs_v(:),phyda_v(:)))])
disp(['mae(paico,vandrsr) = ' num2str(mae(vobs_v(:),paico_v(:)))])

% compute nash-sutcliffe efficiency (aka CE)
disp(['nse(phyda,vandrsr) = ' num2str(nse(vobs_v(:),phyda_v(:)))])
disp(['nse(paico,vandrsr) = ' num2str(nse(vobs_v(:),paico_v(:)))])

% compute crpss
%[mCRPS] = crps_m(f_mu,f_sig,y)
phyda_sig=std(xa_di_ens_sa(dvi,:),0,2);
paico_sig=std(result.resample.signals(:,pvi),0,1);

crps_r=crps_m(phyda_v,phyda_sig,vobs_v);
% resample largest std for recon std with no mean (mean removed above from the time series)
std0=repmat(max(std(result.resample.signals,0,1)),length(dvi),1);
%std0=repmat(max(std(xa_di_ens_sa,0,2)),length(dvi),1);
crps_0=crps_m(zeros(length(dvi),1),std0,vobs_v);
crpss=1-crps_r/crps_0;

disp(['crps(phyda,vandrsr) = ' num2str(crps_r)])
disp(['crpss(phyda,vandrsr) = ' num2str(crpss)])

crps_r=crps_m(paico_v,paico_sig,vobs_v);
% resample largest std for recon std with no mean (mean removed above from the time series)
% use same std0 from above
%std0=repmat(max(std(result.resample.signals,0,1)),length(dvi),1);
%crps_0=crps_m(zeros(length(dvi),1),std0,vobv_v);
crpss=1-crps_r/crps_0;

disp(['crps(paico,vandrsr) = ' num2str(crps_r)])
disp(['crpss(paico,vandrsr) = ' num2str(crpss)])


addpath('./boundedline-pkg/Inpaint_nans/')

figure
hold on
[hl1,hp1]=boundedline(vyrs,phyda_v,2*phyda_sig,'cmap',clrs.maroon,'alpha');
set(hl1,'linewidth',2)
%set(hp1,'alpha',0.35)
%hl1.LineWidth=2;
[hl2,hp2]=boundedline(vyrs,paico_v,2*paico_sig,'cmap',clrs.blue,'alpha');
set(hl2,'linewidth',2)
%set(hp2,'alpha',0.35)
%hl2.LineWidth=2;

%plot(vyrs,phyda_v,'linewidth',2)
%plot(vyrs,paico_v,'linewidth',2)
yt=plot(vyrs,vobs_v,'k','linewidth',2.5);
legend([yt,hl2,hl1],{'VDS','PaiCo','PHYDA'})
xlabel('Years','fontsize',18)
ylabel('PDSI','fontsize',18)
set(gca,'fontsize',16)
ylim([-8,8])
xlim([vyrs(1),vyrs(end)])
hold off
box on

%print(['/home/nsteiger/Dropbox/sa-na-coupled-mega/phyda_paico_chile_verif.png'],'-dpng','-r300');


phydate=prctile(xa_di_ens_sa(ir3,:),[5 95],2);


[~,~,ip3] = intersect(analyrs3,result.times); % indices of recon segment
paicote=prctile(result.resample.signals(:,ip3),[5 95],1);
paicotem=paicote-mean(paicote,2)+mean(phydate)';

cphy=clrs.maroon;
cpai=clrs.blue;

% ANNUAL RECONSTRUCTIONS WITH UNCERTAINTIES
% PHYDA VS PAICO
figure
hold on
%h1=plot(analyrs3,xa_di_na(ir3),'linewidth',1,'color',[clrs.blue 0.5]); % extra number on the end adds transparency

%h3=plot(result.times(ip3),paicotem,'linewidth',1,'color',[cpai 0.5]);
%h5=plot(result.times(ip3),paicotm,'color',cpai,'linewidth',1.25);

%h2=plot(analyrs3,phydate,'linewidth',1,'color',[cphy 0.5]);
%h4=plot(analyrs3,phydatm,'linewidth',1.25,'color',cphy);

s0=5;
h1=plot(result.times(ip3),smooth(paicotem(1,:),s0,'lowess')-phyoff,'linewidth',1,'color',[cpai 0.35]);
h3=plot(result.times(ip3),smooth(paicotem(2,:),s0,'lowess')-phyoff,'linewidth',1,'color',[cpai 0.35]);
h5=plot(result.times(ip3),smooth(paicotm,s0,'lowess')-phyoff,'color',cpai,'linewidth',1.75);

h0=plot(analyrs3,smooth(phydate(:,1),s0,'lowess')-phyoff,'linewidth',1,'color',[cphy 0.35]);
h2=plot(analyrs3,smooth(phydate(:,2),s0,'lowess')-phyoff,'linewidth',1,'color',[cphy 0.35]);
h4=plot(analyrs3,smooth(phydatm,s0,'lowess')-phyoff,'linewidth',1.75,'color',cphy);

plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([1400 analyrs3(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 8 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])
ylim([-7 7])% match Fig 1

legend([h4,h2(1),h5,h3(1)],{'SASW PHYDA','5th, 95th','SASW PaiCo','5th, 95th'},...
'Orientation','horizontal','location','south','fontsize',12)
%legend([h2,h3],{'Chile PHYDA','Chile PaiCo'},'location','northwest')
hold off



print(['/home/nsteiger/Dropbox/phyda_paico_chile_ann14002000_smth' num2str(s0) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/phyda_paico_chile_ann14002000_smth' num2str(s0) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/phyda_paicoAllPrxs_chile.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/phyda_paico_argentina.png'],'-dpng','-r300');










return

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


% 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,:)=[];


%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],clrs.blue);
    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],clrs.purple);
    p.FaceAlpha=0.35;p.EdgeColor='none';
end

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

%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');





