%% Extremes: histograms for different sections, 20 years


cd ..
cd Inputs
load BARCAST_INPUT_vHNL_Jan2012
cd ..
load Parameters_pers3
cd Fig3_Files/

timeline_data=BARCAST_INPUT.Data_timeline;

%variance of the AR(1) process is:
var_AR1= Parameters_pers3(3,2) / (1-Parameters_pers3(1,2)^2);
mu=Parameters_pers3(2,2);

% 
xs=-9:.001:9;
ys=normpdf(xs, mu, sqrt(var_AR1));

%
y_starts=1852:20:1992;
y_ends=1871:20:2011;


%load the mean and var over all
load mean_var_all

% Load the structures
% for the data:
load Histos_20Ys

%for the simulated CRU:
load Sim_CRU_max_20y

%get the number of sims:
N_Sims=length(Sim_CRU_max_20y(1).Max_vals(1,:));

%% DECIDE WHICH 20-YEAR INTERVAL TO CONSIDER:

Int_Select=8; 
%Int_Select=1; %1852:1:1871
%Int_Select=2; %1872:1:1891
%Int_Select=3; %1892:1:1911
%Int_Select=4; %1912:1:1931
%Int_Select=5; %1932:1:1951
%Int_Select=6; %1952:1:1971
%Int_Select=7; %1972:1:1991
%Int_Select=8; %1992:1:2011



%get the relevant counts from the data-constrained ensemble:
cents_20ys=Histos_20Ys(Int_Select).Cents;
counts_20ys=Histos_20Ys(Int_Select).Counts;
%and the mean/var from the data constrained ensemble
mean_var_20ys_Dats=Histos_20Ys(Int_Select).Mean_Var;

%select out the correct subinterval of the simulations
Max_vals_20ys_NoData=Sim_CRU_max_20y(Int_Select).Max_vals;
Top_Twenty_20ys_NoData=Sim_CRU_max_20y(Int_Select).Top_Twenty;

%set the interval
int=y_starts(Int_Select):1:y_ends(Int_Select);
    

% USE THE MEANS ACROSS all ensemble/times/locs TO ADJUST MEANS.
% AFTER ALL, IN ALL OF THESE PLOTS, LOCATIONS, NOT AREAS, RECEIVE EQUAL
% WEIGHT.
Max_vals_20ys_NoData_shift=Max_vals_20ys_NoData-mean_var_all(1)+mean_var_20ys_Dats(1);
Top_Twenty_20ys_NoData_shift=Top_Twenty_20ys_NoData-mean_var_all(1)+mean_var_20ys_Dats(1);


% extract the section of actual CRU data:
indsc=[find(timeline_data==int(1)), find(timeline_data==int(end))];
N_ys=length(int);

CRU_20year_pre=BARCAST_INPUT.Inst_Data(BARCAST_INPUT.Inds_Central, indsc(1):1:indsc(2));
% get non all nan's:
inds_gg=find(sum(isnan(CRU_20year_pre),2)<N_ys);
CRU_20year=CRU_20year_pre(inds_gg, :);

% get the max at each location:
Maxs_CRU_20year=max(CRU_20year, [], 2);

%max of all:
inders_gg=find(isnan(CRU_20year)==0);
CRU_20year_S=sort(CRU_20year(inders_gg));
CRU_20year_Top20=CRU_20year_S(end-19:1:end);


% figure out the varaince of the maxima (at each location) according to the shifted process:
var_proc_max=var(Max_vals_20ys_NoData_shift(:));
%and for the data:
var_data_max=var(Maxs_CRU_20year);
%get the varaince of each ensemble member for shifted process:
var_proc_max_each=var(Max_vals_20ys_NoData_shift);
%find percentile of the data:
per_dat_max=length(find(var_proc_max_each<var_data_max))/length(var_proc_max_each);
%sort the process level values, get indices of relevant percentiles:
[aa,ss]=sort(var_proc_max_each);


% for each location, get the percentile of the CRU value with respect to
% the mean-shifted process level. 

pers_each_loc=NaN(length(Maxs_CRU_20year),1);
NN=length(Max_vals_20ys_NoData_shift(1,:));
for kk=1:1:length(Maxs_CRU_20year)
        pers_each_loc(kk)=length(find(Max_vals_20ys_NoData_shift(kk,:)>Maxs_CRU_20year(kk)))/NN;    
end

%get the mean of the data max, find percentile under the simulations:
mean_max=mean(Maxs_CRU_20year);
length(find(Max_vals_20ys_NoData_shift(:)<mean_max))/length(Max_vals_20ys_NoData_shift(:));

%get the mean of each ensemble member for shifted process:
mean_proc_max_each=mean(Max_vals_20ys_NoData_shift);
%find percentile of the data:
per_dat_max_by_mean=length(find(mean_proc_max_each<mean_max))/length(mean_proc_max_each);


%number of hottest to plot; 
PP=5;


ctrip=[.2,.2,1];
ctrip2=[1,0,1]*.85;
ctrip3=[1,1,1]*0;

ww=.85;
vv=.45;
yy=.15;

hot_trip=[1,1,1]*.5;
clear ax
axers=[-3 5.5 0, .55];

figure(30), clf
set(gcf, 'position', [10 10 1300 750], 'renderer', 'painters')


ax(1)=subplot(3,1,1);
set(gca, 'fontsize', 15)

bar(cents_20ys, counts_20ys/sum(counts_20ys)/diff(cents_20ys(1:2)), 'barwidth', ww, 'facecolor', ctrip, 'edgecolor', ctrip), hold on
plot(xs, ys, 'k-', 'linewidth', 5), hold on
plot(xs, normpdf(xs, mean_var_20ys_Dats(1), sqrt(var_AR1)), 'r--', 'linewidth', 4), hold on

ggg=legend(['Posterior (', num2str(int(1)),'-', num2str(int(end)), ')'],'Simulated anomalies', ['Mean-shifted simulated anomalies'], 'location', 'northeast')

for kk=1:1:PP
    plot([1,1]*CRU_20year_Top20(kk+(20-PP)), [0, 4], 'k--', 'linewidth', 1.25', 'color', [kk/PP,.3,1-kk/PP]), hold on
end

text(-2.85, 0.5, ['(a)'], 'fontsize', 17)

axis([axers])
ylabel('Probability density')


ax(2)=subplot(3,1,2);
set(gca, 'fontsize', 15)


spc=.125;
cents_m=-4+spc/2:spc:(7-spc/2);



Data_to_use=Maxs_CRU_20year(:);
[counts, cts]=hist(Data_to_use(:), cents_m);
bar(cents_m, counts/sum(counts)/diff(cents_m(1:2)), 'barwidth', ww, 'facecolor', ctrip, 'edgecolor', ctrip), hold on


Sims_to_use=Max_vals_20ys_NoData_shift(:, [1:1:N_Sims]);
[counts, cts]=hist(Sims_to_use(:), cents_m);
bar(cents_m, counts/sum(counts)/diff(cents_m(1:2)), 'barwidth', vv, 'facecolor', ctrip3, 'edgecolor', ctrip3), hold on

%plot sim that is closest to data by variance, but within .05 .95.
Sims_to_use=Max_vals_20ys_NoData_shift(:, ss(max(round(N_Sims*min(.95,max(per_dat_max, .05)))-0, 1)));
[counts, cts]=hist(Sims_to_use(:), cents_m);
barplot_tops(counts/sum(counts)/diff(cents_m(1:2)), cents_m, 2, 1, ctrip2), hold on


text(-2.85, .15, ['(b)'], 'fontsize', 17)

legend(['Instrumental obs. (', num2str(int(1)),'-', num2str(int(end)), ')'],'Shifted simulations of inst. obs. (10000 draws)', 'Shifted simulations of inst. obs. (1 draw)', 'location', 'northwest')


for kk=1:1:PP
    plot([1,1]*CRU_20year_Top20(kk+(20-PP)), [0, 4], 'k--', 'linewidth', 1.25', 'color', [kk/PP,.3,1-kk/PP]), hold on
end

axis([-3 5.5 0 1.25])
ylabel('Probability density')


%max across space as well.
ax(3)=subplot(3,1,3);
set(gca, 'fontsize', 15)

kk=1;
plot([1,1]*CRU_20year_Top20(kk+15), [0, 4], 'k--', 'linewidth', 1.25', 'color', [kk/5,.3,1-kk/5]), hold on
[f,xi] = ksdensity(Top_Twenty_20ys_NoData_shift(kk+15,:), 2.5:.005:5.5);
plot(xi, f, 'k-', 'linewidth', 1.5, 'color', [kk/5,.3,1-kk/5]), hold on

legend(['5 Largest inst. obs. (', num2str(int(1)),'-', num2str(int(end)), ')'],'Shifted simulations of inst. obs. : Distributions of 5 largest', 'location', 'northwest')



for kk=1:1:PP
    plot([1,1]*CRU_20year_Top20(kk+(20-PP)), [0, 4], 'k--', 'linewidth', 1.25', 'color', [kk/PP,.3,1-kk/PP]), hold on
end
for kk=1:1:PP
[f,xi] = ksdensity(Top_Twenty_20ys_NoData_shift(kk+(20-PP),:), 1.0:.005:5.5);
plot(xi, f, 'k-', 'linewidth', 1.5, 'color', [kk/PP,.3,1-kk/PP]), hold on
end


axis([-3 5.5 0 1.95])
linkaxes(ax, 'x')
xlabel('Temperature Anomaly in ^oC')

text(-2.85, .25, ['(c)'], 'fontsize', 17)
ylabel('Probability density')