
tic

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

%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_r8002000_p998_state_tq_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_12-Jul-2018_14:50:03.mat';
%efl='sst_MCruns_ensemble_mean.nc';
%efl='sst_ensemble_mean_nino34_regressed_off_productionFinal.nc';
%efl='hgt500_MCruns_ensemble_mean.nc';

pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
%pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])
%pth='/home/nsteiger/MATLAB/steiger-scidata-2018/output-recon/';load([pth,efl])
%pth='/home/nsteiger/MATLAB/hydro-recon/lmr/';


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

efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';disp('Loading DJF data...')
%load([pth,efl],'state_tp','xmeta','Xa_m')
load([pth,efl])
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);
Xt2m=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);


% % LMR data
% %sst0=double(ncread([pth,efl],'sst'));
% x0=double(ncread([pth,efl],'hgt500'));
% lat=double(ncread([pth,efl],'lat'));
% lon=double(ncread([pth,efl],'lon'));
% % MC ensemble mean
% %sstmn=reshape(permute(squeeze(mean(sst0,3)),[2 1 3]),length(lat)*length(lon),size(sst0,4));
% xmn=reshape(permute(squeeze(mean(x0,3)),[2 1 3]),length(lat)*length(lon),size(x0,4));
% 
% %sstmn=reshape(permute(sst0,[2 1 3]),length(lat)*length(lon),size(sst0,3)); % for ENSO regressed out mean
% r_o=0;r_f=2000;

%clearvars -except Xt2m

toc


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

disp(['Analysis years = ' num2str(analyrs0(1)) '-' num2str(analyrs0(end))])

disp('Data preparation...')


% Detrend for regular ensemble mean
%D1=detrend(Xt2m(:,ir0)')';
%D1=detrend(sstmn(:,ir0)')';

% detrend by removing smoothed time series
D0=Xt2m(:,ir0);
D1=zeros(size(D0));
for i=1:size(D0,1)
   D1(i,:)=D0(i,:)'-smooth(D0(i,:),100,'lowess');
end

% don't detrend
%D1=Xt2m(:,ir0);

% FULL SPATIAL GRID, NOT JUST SSTS
%D=detrend(xmn(:,ir0)')';
%W=cosd(repmat(lat,[1 length(lon)]));wt=W(:);

% for full ensemble: Detrend data b/c of long term cooling trend
%D1=zeros(length(lat)*length(lon),length(ir0),sub_ens);
%for j=1:sub_ens
%end
%D1=reshape(D1,length(lat)*length(lon),length(ir0)*sub_ens);

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


disp('Just select tropical values')
% Make weighting and ocean masking
varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
sftlf0=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);

%itlat=find(lat<45 & lat>-45);
itlat=find(lat<60 & lat>-60); % full Pacific
itlon=find(lon>100 & lon<300); % corresponds to 130 E to 60 W
sftlf=sftlf0(itlat,itlon);

slat=lat(itlat);
slon=lon(itlon);

omsk=sftlf(:);
oids=find(omsk==0);

W=cosd(repmat(lat(itlat),[1 length(lon(itlon))]));W=W(:);
wt=W(oids);
% Just take ocean values
D1=reshape(D1,length(lat),length(lon),length(ir0));
D2=reshape(D1(itlat,itlon,:),length(itlat)*length(itlon),length(ir0));
D=D2(oids,:);


%--------------------------------------
% UNCOMMENT BELOW FOR FULL GLOBAL SSTS
%--------------------------------------
% % Make weighting and ocean masking
% 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]);
% 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,:);
% %clear D1 D0 Xt2m





% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% UNCOMMENT BELOW FOR SST VARIABLE
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% % Make weighting and ocean masking
% % SOMETHING IS WRONG WITH THE MASKING VARIABLE, IT'S MISSING LOCATIONS...
% %sftlf=permute(ncread([pth,'oceanmask_landmask_lmr.nc'],'landmask'),[2 1]);
% %omsk=sftlf(:);
% oids=find(~isnan(sstmn(:,1)));
% W=cosd(repmat(lat,[1 length(lon)]));W=W(:);
% wt=W(oids);
% % Just take ocean values
% D=D1(oids,:);


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

toc

syrs=analyrs0;

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...')
[ptrns0,bmus,rmsErr,pfrq,sMap,sD] = som_make_njs(D,somRows,somCols,init,wt);
disp('SOM completed');toc

%clear D

% Expand out the ptrns
%ptrns=NaN(somRows*somCols,length(lat)*length(lon));
ptrns1=NaN(somRows*somCols,length(itlat)*length(itlon));
ptrns=NaN(somRows*somCols,length(itlat),length(itlon));
ptrns1(:,oids)=ptrns0;
for i=1:somRows*somCols
   ptrns(i,:,:)=reshape(ptrns1(i,:),length(itlat),length(itlon));
end


%ptrns=ptrns0;disp('Expanding patterns for full spatial field only...')
mf='n';
if mf=='y'
% CHECK THAT SOMS AND BMUS MAKE SENSE
figure
for i=1:somRows*somCols
  subplot(somRows,somCols,i)
  h=pcolor(squeeze(ptrns(i,:,:)));
  set(h,'edgecolor','none')
  caxis([-1.5 1.5])
end
ht=suptitle('SOM Patterns');
set(ht,'interpreter','none')
print(['/home/nsteiger/Dropbox/enso-volcanoes/sompatterns_' num2str(somRows*somCols) '_phyda_' num2str(analyrs0(1)) '-' num2str(analyrs0(end)) '_smoothDT.png'],'-dpng','-r300');

figure
plot(syrs,bmus,'o')
title('SOM BMUs')
print(['/home/nsteiger/Dropbox/enso-volcanoes/sombmus_' num2str(somRows*somCols) '_phyda_' num2str(analyrs0(1)) '-' num2str(analyrs0(end)) '_smoothDT.png'],'-dpng','-r300');

end



disp('Saving output...')


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


%save(['output_phyda_djf_tropicalPlus_som' num2str(somRows*somCols) '_yrs' num2str(analyrs0(1)) num2str(analyrs0(end)) '_' timestp2 '.mat'],'-v7.3')
save(['output_phyda_djf_pacific_soms' num2str(somRows*somCols) '_row' num2str(somRows) '_col' num2str(somCols) '_yrs' num2str(analyrs0(1)) num2str(analyrs0(end)) '_' timestp2 '.mat'],'ptrns','bmus','somRows','somCols','slat','slon','syrs','-v7.3')

toc

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











