%% simulate 4000 samples of CRU obs of temperature anomalies over 20 year
% intervals, beginning with  1852--1871, in the absence of any data. To check on the extremes. 

%Use median process-level parameters

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

% order of the scalar parameters
%1 = alpha, the AR(1) coefficient
%2 = mu, the constant par in the linear mean of the AR(1) process
%3 = sigma2, the partial sill in the spatial covariance matrix
%4 = phi, the range parameter in the spatial covariance matrix
%5 = tau2_I, the Inst measurement error

alpha=Parameters_pers3(1,2);
mu=Parameters_pers3(2,2);
sig2=Parameters_pers3(3,2);
phi=Parameters_pers3(4,2);
tau=sqrt(Parameters_pers3(5,2));

% get the locations:
Locs=BARCAST_INPUT.Master_Locs(BARCAST_INPUT.Inds_Central,:);
N_Locs=length(Locs(:,1));
% get the distance matrix:
All_DistMat=EarthChordDistances(Locs);

%spat cov:
spat_cov=sig2*exp(-phi*All_DistMat);
chol_spatcov=chol(spat_cov);


% set the current value of the TempField:
TempField=zeros(N_Locs,1)+mu;

% propagate forward in time a number Npre of times, to get to convergence:
Npre=100;
for kk=1:1:Npre    
    TempField=alpha*TempField + mu*(1-alpha) + chol_spatcov'*randn(N_Locs, 1);
end

%set the number of time points:
N_ys=20;

%number of iterations
N_its=10000;
%number of time steps to propagate forward between iterations:
%(to reduce correlation between iteraitons)
N_sp=10;

%% Save eveything as a strucutre. 
%Get a mask, want only those times and locations over the 20 year periods
%when/where there are CRU obs. 
timeline_data=BARCAST_INPUT.Data_timeline;

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

N_Ints=length(y_starts);


Temp_pars=struct();

% need to save:
%Inds_gg_1
%N_Locs_gg_1
%Inst_Obs_1_sub_mask
%inders_nonnan_1

for kk=1:1:N_Ints
    
    %indices of the years
    indsc=[find(timeline_data==y_starts(kk)), find(timeline_data==y_ends(kk))];
    Inst_Obs=BARCAST_INPUT.Inst_Data(BARCAST_INPUT.Inds_Central, indsc(1):1:indsc(2));
    
    % get the indices of locs where there are obs:
    N_Inst_Obs=N_ys-sum(isnan(Inst_Obs), 2);
    Inds_gg=find(N_Inst_Obs>0);
    N_Locs_gg=length(Inds_gg);
    Inst_Obs_sub=Inst_Obs(Inds_gg, :);
    
    %replace all the non-NaN's with ones, so we can use element-wise
    %multiplication below:
    Inst_Obs_sub_mask=Inst_Obs_sub;
    inders_nonnan=find(isnan(Inst_Obs_sub_mask)==0);
    Inst_Obs_sub_mask(inders_nonnan)=1;
    
    Temp_pars(kk).Inds_gg=Inds_gg;
    Temp_pars(kk).N_Locs_gg=N_Locs_gg;
    Temp_pars(kk).Inst_Obs_sub_mask=Inst_Obs_sub_mask;
    Temp_pars(kk).inders_nonnan=inders_nonnan;
    
    
    clear Inds_gg N_Locs_gg Inst_Obs_sub_mask inders_nonnan
    
end
%
%% set empty arrays to fill

Sim_CRU_max_20y=struct();

for kk=1:1:N_Ints
Sim_CRU_max_20y(kk).Max_vals=NaN(Temp_pars(kk).N_Locs_gg, N_its);
Sim_CRU_max_20y(kk).Top_Twenty=NaN(20, N_its);

end


%% now go through the iterations (outer) and times (inner)

% for each iteration, generate the N_Locs by N_Ys array, then add the inst
% obs noise, and then take max across time at each loc. 

%Means_All=NaN(N_its, 1);
for KK=1:1:N_its
    
    %Empty Array:
    Space_Time_T=NaN(N_Locs, N_ys);
    for JJ=1:1:N_ys
        
        %take the simulation forward in time by one:
        TempField=alpha*TempField + mu*(1-alpha) + chol_spatcov'*randn(N_Locs, 1);
        
        %fill the temp field:
        Space_Time_T(:, JJ)=TempField;
    end
    
    
    % get maxima stuff for each interval
    for PP=1:1:N_Ints

        %subset to the locations where there are obs, and add the obs error
        %noise
        Space_Time_T_sub=Space_Time_T(Temp_pars(PP).Inds_gg,:)+randn(Temp_pars(PP).N_Locs_gg, N_ys)*tau;
        %hit with the NaN/1 mask to mimic the CRU:
        Space_Time_T_sub=Space_Time_T_sub.*Temp_pars(PP).Inst_Obs_sub_mask;
        
        %do the max stuff:
        Sim_CRU_max_20y(PP).Max_vals(:,KK) =  max(Space_Time_T_sub,[],2);
        
        %also get the top X values:
        temp=sort(Space_Time_T_sub(Temp_pars(PP).inders_nonnan));
        Sim_CRU_max_20y(PP).Top_Twenty(:,KK)=temp(end-19:1:end);

        clear Space_Time_T_sub temp
        
    end
    
    for MM=1:1:N_sp
        TempField=alpha*TempField + mu*(1-alpha) + chol_spatcov'*randn(N_Locs, 1);
    end
    KK
end

%
save Sim_CRU_max_20y Sim_CRU_max_20y
