

%=================
% LOAD PROXY DATA
%=================

disp('Loading data...')
addpath('../')

% ALL PROXY TYPES, LMR DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
% INDEX ONLY RECONSTRUCTION WITH FULL ENSEMBLES
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_16-Nov-2017_15:57:59.mat';


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


%varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
%sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
%lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
%lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');



%-------------------------
% LOAD SPATIAL FIELDS
%-------------------------

% Temperature
k=find(state_tp=='t');
Xvar=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
Xt2m=reshape(Xvar,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);

% k=find(state_tp=='n');i=3; % = Nino 3.4
% xa_nino=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
% 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);
% 
% k=find(state_tp=='z');i=3; % = Pacific (Schneider)
% xa_itcz_pac=Xa_m((xmeta{k}.id_X(1)+i-1),:);
% xb_itcz_pac=Xb_o((xmeta{k}.id_X(1)+i-1),:);
% i=8; % = South America
% xa_itcz_sa=Xa_m((xmeta{k}.id_X(1)+i-1),:);
% xb_itcz_sa=Xb_o((xmeta{k}.id_X(1)+i-1),:);
% i=1; % = Atlantic
% xa_itcz_atl=Xa_m((xmeta{k}.id_X(1)+i-1),:);
% xb_itcz_atl=Xb_o((xmeta{k}.id_X(1)+i-1),:);
% 
% k=find(state_tp=='o'); % AMO
% xa_amo=Xa_m((xmeta{k}.id_X(1)+i-1),:);

%clearvars -except xa_amo xb_itcz_sa xb_itcz_atl xb_itcz_pac xa_itcz_sa xa_itcz_atl xa_itcz_pac Xt2m xa_nino xb_nino

clearvars -except Xt2m

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


% What years to look at?
analyrs1=800:1600;
analyrs2=1601:1925;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
[~,~,ir2] = intersect(analyrs2,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;


% Expand out PDSI land-only fields from reconstruction
% NEED TO MAKE 'LNDIDX' AVAILABLE IN THE OUTPUT FILE!!! MAYBE LOAD IT AS A GLOBAL VARIABLE SOMEWHERE? OR JUST LOAD IT IN PARAMS SCRIPT?
%load(['/d1/nsteiger/climate-data/cesm-lme/cesm_lme_010_pdsi_output_biascorr_AWC_c_7_u2_gcm.mat'],'lndidx')
%varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
%sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);

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

lat=xmeta{k}.lat;
lon=xmeta{k}.lon;

% PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
%load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
%lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
%lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');
% NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% Find nearest lat/lon equal to or outside bounds
lonext=[235 255];latext=[31 42];
ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
nasw=zeros(size(sftlf));
nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);

nasw(nasw==0)=NaN;
Xv_msk=bsxfun(@times,Xpdsi,nasw);
A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));
Xa_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX

% Drought index over the analysis years
di=Xa_di(ir1);
dtyp=1; % 1 = 11 year smooth; 2 = 2S2E
disp(['Computing drought type as ' num2str(dtyp)])
ds0=mdrought(di,dtyp,'y'); % 'y' indicates that the mean is removed

%disp('Loading data...')
%load('md_som_anal_input.mat','ds0','r_o','r_f','Xt2m','lat','lon','analyrs1','analyrs2')



% Compute SOMs over the entire period of interest
analyrs0=800:1925;
[~,~,ir0] = intersect(analyrs0,r_o:r_f);

% Detrend data b/c of long term cooling trend
D0=reshape(Xt2m(:,:,ir0),length(lat)*length(lon),length(ir0));
D1=detrend(D0')';% Detrend the data, if desired

% Non-detrended data
%D=reshape(Xt2m(:,:,ir0),length(lat)*length(lon),length(ir0));


% Make weighting and ocean masking
omsk=sftlf(:);
oids=find(omsk==0);
W=cosd(repmat(lat,[1 length(lon)]));W=W(:);
wt=W(oids);
% Just take ocean values
D=D1(oids,:);

%==================
% COMPUTE THE SOMS
%==================


somRows=4;somCols=2;
disp(['SOMs = ' num2str(somRows*somCols)])
init='l'; % (l)inear or (r)andom initialization
%wt='y'; % (y)es or (n)o latitude weighting
% D[lat*lon,time]
disp('Computing SOM...')
tic
[ptrns0,bmus,rmsErr,pfrq,sMap,sD] = som_make_njs(D,somRows,somCols,init,wt);
disp('SOM completed');toc


% Expand out the ptrns
ptrns=NaN(somRows*somCols,length(lat)*length(lon));
ptrns(:,oids)=ptrns0;



%figure
%for i=1:somRows*somCols
%  subplot(somRows,somCols,i)
%  h=pcolor(reshape(ptrns(i,:),length(lat),length(lon)));
%  set(h,'edgecolor','none')
%end
%ht=suptitle(['' fn(1:enddd-13)]);
%set(ht,'interpreter','none')

timestp2=strrep(datestr(now),' ','_');

save(['output_som' num2str(somRows*somCols) '_d_' num2str(dtyp) '_' timestp2 '.mat'])

% %----------------------------
% % PLOT BMUS FOR EACH PATTERN
% %----------------------------
% 
% % BMU frequency during droughts (note droughts are determined analysis period relative)
% ds_pfrq=zeros(length(ds0),somRows*somCols);
% for i=1:length(ds0)
%   ds_pfrq(i,:)=histcounts(bmus(ds0(i,1):ds0(i,2)),0.5:1:(somRows*somCols+0.5))/length(ds0(i,1):ds0(i,2));
% end
% 
% % BMU freq for pre-1600
% [~,ay2,~] = intersect(analyrs0,analyrs2);
% a2_pfrq=histcounts(bmus(ay2),0.5:1:(somRows*somCols+0.5))/length(ay2);
% 
% 
% figure
% hold on
% % Mean freq during mega D
% bar(1:somRows*somCols,mean(ds_pfrq))
% % Mean freq 1601-1925
% plot(1:somRows*somCols,a2_pfrq,'k*')
% title('Ptrn. Freq. MDrought')
% xlabel('Pattern')
% ylabel('Freq.')
% 
% figure
% plot(analyrs0,bmus,'or')
% ylabel('Pattern')
% xlabel('Year')
% 
% 
% % Model the transition frequencies as Markov chain
% 
% % Compute transition matrix
% [T1] = transM(bmus,1);
% 
% % Markov chain object
% mc = dtmc(T1);
% 
% % Directed graph with color = probabilities
% figure
% graphplot(mc,'ColorEdges',true);



% figure
% hold on
% plot(mon_r(irm1_a:irm1_b),zscore(xa_nino(irm1_a:irm1_b)))
% plot(ir1,zscore(Xa_di(ir1)))
% legend('Nino','PDSI')
% y=zscore(Xa_di(ir1));
% yMin=min(y)-0.5;yMax=max(y)+0.5;
% 
% % ADD NEXT IN THE PLOT VERTICAL BARS TO SEE PERIODS OF DROUGHT
% % Add a patches for top five times from truth
% for i=1:length(ds)
%     p=patch(ir1(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],'k');
%     p.FaceAlpha=0.25;p.EdgeColor='none';
% end
% xlim([ir1(1) ir1(end)])
% 
% if prntfg=='y';print(['figs/nino_pdsi_megaD_highlight_prd1.png'],'-dpng');end











