

load distinct_colors_20.mat

%-------------------
% KMEANS CLUSTERING
%-------------------

% Load pre-computed k-means for consistency (sometimes algorithm struggles to converge)
load('kmeans_clusters_8-May-2019.mat')

% Sort by amount of precip
[cent_s,is]=sort(cent);

%rnc{1}=nrn; % no rain
rnc{1}=find(cidx==is(1)); 
rnc{2}=find(cidx==is(2)); 
rnc{3}=find(cidx==is(3)); 
rnc{4}=find(cidx==is(4)); 


%==============
% SOM analysis
%==============

%load('output_som30_era5_mslp_26-Mar-2019_20:04:25.mat','bmus')
%load('output_som35_era5_mslp_26-Mar-2019_20:08:26.mat','bmus')
%load('output_som40_era5_mslp_26-Mar-2019_21:24:30.mat','bmus')

%load('output_som40_era5_mslp_23-May-2019_16:42:35.mat','bmus','somRows','somCols')
%load('output_som48_era5_mslp_23-May-2019_16:26:27.mat','bmus','somRows','somCols')
%load('output_som54_era5_mslp_23-May-2019_16:42:18.mat','bmus','somRows','somCols')

% standardized data, smaller area, more soms
%load('output_som48_era5_mslp_27-May-2019_21:40:32.mat','bmus','somRows','somCols','slat','slon','ptrns');
%load('output_som54_era5_mslp_27-May-2019_21:42:16.mat','bmus','somRows','somCols','slat','slon','ptrns');
%load('output_som63_era5_mslp_24-Jun-2019_20:23:04.mat','bmus','somRows','somCols','slat','slon','ptrns');
%load('output_som72_era5_mslp_24-Jun-2019_21:07:41.mat','bmus','somRows','somCols','slat','slon','ptrns','dt','xv2_ds','eic_lat','eic_lon');
load('output_som81_era5_mslp_25-Jun-2019_11:56:59.mat','bmus','somRows','somCols','slat','slon','ptrns','dt','xv2_ds','eic_lat','eic_lon');
%load('output_som90_era5_mslp_25-Jun-2019_14:50:15.mat','bmus','somRows','somCols','slat','slon','ptrns','dt','xv2_ds','eic_lat','eic_lon');

% histcounts for observed pattern 4
hp4=histcounts(bmus(rnc{4}),0.5:1:(somRows*somCols+0.5))/length(rnc{4});


% Randomly sample BMU distribution and find out what counts as > 2*sigma occurrences

rsz=1000;
sm_rd=zeros(rsz,somRows*somCols);
for i=1:rsz
   rs=randsample(bmus,length(rnc{4}),'true');
   sm_rd(i,:)=histcounts(rs,0.5:1:(somRows*somCols+0.5))/length(rs);
end

% Get percentiles of the above uncertainty estimates
sm_rdp=prctile(sm_rd,[2.5,97.5],1);

% % Find which values lay outside this range 
%outrng=[];j=1;
%for i=1:somRows*somCols
%   if hp4(i)<sm_rdp(1,i)||hp4(i)>sm_rdp(2,i)
%      outrng(j)=i;
%      j=j+1;
%   end
%end

% Values outside upper range only
outrng=[];j=1;
for i=1:somRows*somCols
   if hp4(i)>sm_rdp(2,i)
      outrng(j)=i;
      j=j+1;
   end
end

% SOM frequency vs SOMs during clusters
figure
histogram(bmus(rnc{4}),somRows*somCols,'normalization','pdf')
hold on
plot(1:somRows*somCols,sm_rdp','-o','color',clrs.grey) 
plot(outrng,0.065*ones(length(outrng)),'*','color',clrs.black)
title('Cluster 4')
ylim([0 0.07])
ylabel('Frequency (norm)')
xlabel('SOM')
set(gca,'fontsize',14)
%print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_hist_' num2str(somRows*somCols) '.png'],'-dpng','-r300');



% PLOT OUTLIER SOMS

C=load('./colormaps/NCV_blue_red.mat');cmap_br=C.cmap;
C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;

[X,Y] = meshgrid(slon,slat);

%----------------
% PLOT SOMS
%----------------

% Precip
xf=reshape(xv2_ds,length(slat),length(slon),length(dt));
%xf=reshape(xv2_ds,length(slat),length(slon),length(1:10*366));
xvmn=squeeze(mean(xf,3));


%figure
for i=1:length(outrng)
   figure
   subplot(2,1,1)
   h=contour(X,Y,reshape(ptrns(outrng(i),:),length(slat),length(slon)));
   caxis([-1.5 1.5])
%   caxis([-2000 2000])
   colormap(cmap_br)
   hold on
   plot(eic_lon,eic_lat,'kp')
   set(gca,'xtick',[],'ytick',[])
   %title(['MSLP ' num2str(outrng(i))])
   
   subplot(2,1,2)
   cind=find(bmus==outrng(i));
   %cmp=mean(xf(:,:,cind),3)-xvmn;
   cmp=mean(xf(:,:,cind),3);
   %h=contour(X,Y,cmp);
   h=pcolor(X,Y,cmp);
   set(h,'EdgeColor','none');
   caxis([-4e-3 4e-3]) % Make the axis a log-axis?
   colormap(cmap_bg)
   hold on
   plot(eic_lon,eic_lat,'kp')
   set(gca,'xtick',[],'ytick',[])
   suptitle(['SOM ' num2str(outrng(i))])

   fig = gcf;fig.PaperUnits = 'inches';
   fig.PaperPosition = [0 2 9 6.5]; % [left bottom width height]
%   print(['/home/nsteiger/Dropbox/easterisle/ei_soms_' num2str(somRows*somCols) '_unfrq_' num2str(outrng(i)) '_mslp_pcp.png'],'-dpng','-r300');
end
%suptitle('MSLP')
%fig = gcf;fig.PaperUnits = 'inches';
%fig.PaperPosition = [0 2 9 6.5]; % [left bottom width height]

% %timestp2=strrep(datestr(now),' ','_');
%print(['/home/nsteiger/Dropbox/easterisle/ei_soms' num2str(somRows*somCols) '_mslp_' timestp2 '.png'],'-dpng','-r300');

%---------------------------
% COMPOSITE OTHER VARIABLES
%---------------------------

% % Precip
% xf=reshape(xv2_ds,length(slat),length(slon),length(dt));
% %xf=reshape(xv2_ds,length(slat),length(slon),length(1:10*366));
% xvmn=squeeze(mean(xf,3));
% 
% %figure
% for i=1:length(outrng)
%    %subplot(somRows,somCols,i)
%    figure
%    cind=find(bmus==outrng(i));
%    %cmp=mean(xf(:,:,cind),3)-xvmn;
%    cmp=mean(xf(:,:,cind),3);
%    %h=contour(X,Y,cmp);
%    h=pcolor(X,Y,cmp);
%    set(h,'EdgeColor','none');
%    caxis([-4e-3 4e-3]) % Make the axis a log-axis?
%    colormap(cmap_bg)
%    hold on
%    plot(eic_lon,eic_lat,'kp')
%    set(gca,'xtick',[],'ytick',[])
%    title(['\Delta PCP ' num2str(outrng(i))])
% end
% %suptitle('\Delta PCP')
% 
% %fig = gcf;fig.PaperUnits = 'inches';
% %fig.PaperPosition = [0 2 9 6.5]; % [left bottom width height]
% %print(['/home/nsteiger/Dropbox/easterisle/ei_cmp' num2str(somRows*somCols) '_pcp_' timestp2 '.png'],'-dpng','-r300');




% What do the transitions into and out of key rain events look like?

% Compute transition matrix
[tM cS] = transM(bmus,1);
% Normalize it and convert to Matlab object
mc=dtmc(tM);
% Pij = prob of i-->j (along the row, which sums to 1)

% Find highest probability transition (other than self)

pbfaf=zeros(length(outrng),5);
pbfaf(:,3)=outrng;
for i=1:length(outrng)
   % After key node
   p0=mc.P(outrng(i),:);
   p0(outrng(i))=0; % not interested in self-probability
   [m0,im0]=max(p0); % step+1
   pbfaf(i,4)=im0;

   p1=mc.P(im0,:);
   p1(im0)=0;
   [m1,im1]=max(p1); % step+2
   pbfaf(i,5)=im1;

   % Before key node search columns
   p2=mc.P(:,outrng(i));
   p2(outrng(i))=0;
   [m2,im2]=max(p2); % step-1
   pbfaf(i,2)=im2;

   p3=mc.P(:,im2);
   p3(im2)=0;
   [m3,im3]=max(p3); % step-2
   pbfaf(i,1)=im3;
end


for i=1:length(outrng)
figure
for j=1:5
   subplot(5,1,j)
   h=contour(X,Y,reshape(ptrns(pbfaf(i,j),:),length(slat),length(slon)));
   caxis([-1.5 1.5])
%   caxis([-2000 2000])
   colormap(cmap_br)
   hold on
   plot(eic_lon,eic_lat,'kp')
   set(gca,'xtick',[],'ytick',[])
   %title(['MSLP ' num2str(outrng(i))])
   
   fig = gcf;fig.PaperUnits = 'inches';
   fig.PaperPosition = [0 2 9 6.5]; % [left bottom width height]
%   print(['/home/nsteiger/Dropbox/easterisle/ei_soms_' num2str(somRows*somCols) '_unfrq_' num2str(outrng(i)) '_mslp_pcp.png'],'-dpng','-r300');
end
end
