%% Code to run sizer analysis for detecting climate changepoints
%
% This file is called from within the "calculate_changepoints_usingsizer.m"
% routine.
%
% prepared and archived by Nerilie Abram, June 2016.
% for Abram et al., and the PAGES 2k Consortium (2016), Nature.
% see "calculate_changepoints_usingsizer.m" for full documentation.

%% Set up parameters for SiZer analysis
% see /sizer/sizerSM.m for description of all possible parameters
%
param = struct('imovie',0,'idatovlay',0,'ishoweffwind',0);

%define filtering parameters used in Abram et al (2016)
%number of bins (x-increments)
param.nbin = p2k_cp(1,1)-p2k_cp(end,1)+1;
%number of filters applied for family and sizer plots
param.nfh = 46;
param.nsh = 46;
%minimum filter length
param.fhmin = 5;
param.shmin = 5;
%maximum filter length
param.fhmax = 50;
param.shmax = 50;
%significance level and method
param.alpha = 0.1; %90% confidence
param.simflag = 2; %0 for less conservative (pointwise CI), c.f. 2 (default) = row-wise CI

%nja note: output of family filters is defined in sizerSM_edit, line 1024


%define plot labels and save output
%param.savestr = 'IndianCPS_sizer_rawoutput'
param.xlabelstr = ['year CE'];
param.ylabelstr = ['Temperature anomaly (\circC)'];
param.famoltitle = ['Reconstruction: ' int2str(z)]; %looping through iterations of data using for z = ...
param.sizertitle = ['SiZer map of significant trends'];
param.labelfontsize = 14;
param.titlefontsize = 15;
param.hhighlight = 0; %->> 0 = don't highlight h

%% Run sizer (with 90% trend confidence setting)
cd sizer
sizerSM_edit(p2k_cp,param)

%% Get filter output and save 90% confidence figure

%pull up filter output
load 'output_1'
p2k_filters = mfam;
p2k_filteryear = fxgrid;

% sizer figure for 90% trend significance
sizer_sig90 = gcf;

cd ../working

save('filters','p2k_filters','p2k_filteryear')
saveas(sizer_sig90,'sizer_sig90.fig')

cd ../

%% Rerun sizer (0% trend confidence - i.e. sign of trend only)

%change significance level in parameter file
param.alpha = 1;

%Run sizer
cd sizer
sizerSM_edit(p2k_cp,param)

cd ../working

%save all trends 
sizer_sig0 = gcf;
saveas(sizer_sig0,'sizer_sig0.fig') % for extracting sizer data

close all
clear all

%% Extract trend vs filter data from sizer figures, and save to working folder

%Get data for sig90 sizer run
data=hgload('sizer_sig90.fig');
get(data);
ch=get(data,'Children');
l=get(ch(1,1),'Children');
x_sizer=get(l,'Xdata');
y_sizer=get(l,'Ydata');
z_sizer = getimage(l);
sizer_z90 = z_sizer;
close gcf;

%Get data for sig90 sizer run
data=hgload('sizer_sig0.fig');
get(data);
ch=get(data,'Children');
l=get(ch(1,1),'Children');
x_sizer=get(l,'Xdata');
y_sizer=get(l,'Ydata');
z_sizer = getimage(l);
sizer_z0 = z_sizer;
close gcf;

%Create new data matrix from 0% and 90% significance sizer runs for plotting
%combined sizer plot (shows all trends plus significant trends)

%key for z90 matrix:
% 1 = sig warming, 3 = neutral, 4 = sig cooling, 2 = NaN;

%key for z0 matrix:
% 1 = warming, 4 = cooling, 2 = NaN;

%replace NaNs
sizer_z0(find(sizer_z0 == 2)) = -999;
sizer_z90(find(sizer_z90 == 2)) = -999;

%combine into single matrix (zall)
sizer_zall = (sizer_z0 + sizer_z90)/2; %average together
sizer_zall(find(sizer_zall == 3.5)) = 3; 
sizer_zall(find(sizer_zall < 0)) = -999;

%Key to data in sizer zall output matrix: 
%NaN = -999
%1 = significant warming
%2 = warming
%3 = cooling
%4 = significant cooling

%Get x and y data corresponding to matrix (z)
sizer_x = [x_sizer(1):1:x_sizer(end)];
sizer_y = [y_sizer(1):1:y_sizer(end)];

save('sizer','sizer_*')

close all
clear all

cd ../