%% Detecting the onset of sustained,significant trends in paleoclimate data
%
% SiZer changepoints: version 1.0
%
% Code archive for publication:
% Abram, N.J., McGregor, H.V., Tierney, J.E., Evans, M.N., McKay, N.P., 
% Kaufmann, D. and the PAGES 2k Consortium. 2016. Early onset of 
% industrial-era warming across the oceans and continents. Nature.
%
% Please see publication for description of method.
%
% Code prepared by Nerilie Abram, The Australian National University
% Operating with MATLAB r2014a software, and MacOS 10.10 (Yosemite)
% Archived June 2016.
%
%
% Based on adaption of code by Steve Marron:
% SIZERSM, SIgnificance of ZERo crossings of derivatives,
%   Steve Marron's matlab function
%     For determining which features in a smooth 
%     (density estimate, or nonparametric regression)
%     are statistically significant.
%
% Original code available at: http://www.unc.edu/marron/marron_software.html
% 
% Original publication: Chaudhuri, P. & Marron, J. S. 
% SiZer for exploration of structure in curves. 
% Journal of the American Statistical Association 94, 807-823 (1999)
%
% *** Please cite original reference as well as Abram et al., (2016) when
% using this code package ***
%
%
% NOTE: parameters within this file, as well as "sizer_changepoints" and
% "sizer_run" may need to be changed if analysis needs are different to
% those used in Abram et al., (2016). Files are set up to screen for the
% onset of sustained, significant warming or cooling by analysing annually
% resolved paleoclimate reconstructions and model output since 1500CE. 
% Changepoints are assessed across filter widths from 15-50y and where 
% sustained recent trends are not detected these years are set to 2015.
% Input data files need to be in annual or pseudo-annual resolution for
% current analysis set up.
%
% Please contact Nerilie Abram (nerilie.abram@anu.edu.au) for advice on 
% using this code.

%% Data descriptions

% INPUT FILE:
% Required input data is: 
%   "y": a n by 1 matrix for time (with a timestep of 1), and 
%   "rec": a n by m matrix for reconstruction(s) data
% Input data is placed in the "input" folder. A simple example is given, 
% based on the PAGES2k and Ocean2k regional temperature reconstructions 
% since 1500CE. 
% The full input datasets used in Abram et al., 2016 (archived as
% Supplementary Data 2) are set up as MATLAB structures. Because of the
% different shapes of these structures (between regional reconstructions, 
% site-level records and model output) some modification of the data preparation 
% and looping/saving routines below is required to run this code package
% on the large data structures used in Abram et al., 2016. Please contact
% nerilie.abram@anu.edu.au for advice on code modifications.

% OUTPUT FILES:
% Output includes:
%   1. sizer plots for each reconstruction member (numbered according to
%   column m in input matrix "rec"), and for sizer screening of all trends
%   and for trends significant at p<0.1 level. Plots are saved in a
%   subfolder of the "output" folder
%
%   2. An output data file "cp_output", saved in the "output folder".
%   cp_output is a matlab structure (1-by-m size). Its contents for each
%   reconstruction member are:
%       ens = ensemble member (m of the n-by-m input "rec" file).
%
%       rec = input data used for SiZer analysis, with year in column 1 and
%       climate data in column 2.
%
%       filter_y = year for the sizer filters.
%
%       filters = the family of filters of rec based on the filter
%       widths screened in the sizer analysis.
%
%       sizer_x = x-dimension of sizer output (year).
%
%       sizer_y = y-dimension of sizer output (filter width).
%
%       sizer_z0 = sizer output matrix for occurence of warming and cooling
%       trends according to year (x) and filter width (y). 
%       Key for z0 matrix: 1 = warming, 4 = cooling, 2 = NaN;
%
%       sizer_z90 = sizer output matrix for occurence of significant warming and cooling
%       trends according to year (x) and filter width (y). 
%       Key for z90 matrix: 1 = sig warming, 3 = neutral, 4 = sig cooling, 2 = NaN;
%
%       sizer_zall = combination of sizer_z0 and sizer_z90 to show trends
%       and their significance according to year (x) and filter width (y).
%       Key for zall matrix: 1 = significant warming, 2 = warming, 3 = 
%       cooling, 4 = significant cooling, NaN = -999. The zall matix is
%       used for plotting the SiZer trend maps shown in Extended Data
%       Figure 2 and Supplementarty Figure 1 of Abram et al., 2016.
%
%       cp = a sub-structure containing output from the changepoint
%       screening
%
%
%   The "cp" sub-structure contains for each reconstruction:
%
%       sig_year = the year when the most recent sustained, significant
%       trend began for each filter width (set to screen across 5-50y).
%
%       sig_class = denotes whether the most recent sustained, significant
%       trend was a warming (1) or cooling (-1) trend for each filter
%       width.
%
%       sigwarm_year = the year when the most recent sustained, significant
%       warming trend began for each filter width (max_year used to flag
%       filters where recent sustained, significant warming is not
%       detected; set to 2015 in Abram et al., 2016 study).
%
%       sigcool_year = the year when the most recent sustained, significant
%       cooling trend began for each filter width (max_year used to flag
%       filters where recent sustained, significant cooling is not
%       detected; set to 2015 in Abram et al., 2016 study).
%
%       sigwarm_median = median year for the onset of sustained, significant
%       warming across chosen filter widths (15-50y in Abram et al., 2016).
%
%       sigcool_median = median year for the onset of sustained, significant
%       cooling across chosen filter widths (15-50y in Abram et al., 2016).
%
%       filter_clip = first first-width used for assessing the median onset
%       of warming (e.g. 11 corresponds to the 15y filter-width).
%
%       max_year = year used to indicate no recent significant trend (set
%       to 2015 for Abram et al., 2016 study).
%

%% Data preparation

close all
clear all

%load input data

cd input
load example
cd ../

%save input data to working folder
cd working
save('input')

%clear existing output file in working save file 
R = [];
save('R','R')
clear R
cd ../

%% run sizer analysis on climate reconstructions

%loop sizer analysis through a set of reconstructions, stored in input file
%'rec', with associated input file 'y'.
for z = 1:size(rec,2); 

%save iteration number and crop year in working folder
cd working
save('it','z')
cd ../ 

%set up data
data = rec(:,z);

p2k_cp = [y, data];
p2k_cp=p2k_cp(~any(isnan(p2k_cp),2),:); %remove any rows with NaNs

%check orientation of data, if needed flip to run young to old.
if p2k_cp(2,1)-p2k_cp(1,1) > 0;
p2k_cp = flipud(p2k_cp); 
else
end

cd working
save('rec','p2k_cp')
cd ../

%run sizer routine
sizer_run

%load sizer data from working folder
cd working
load filters 
load sizer
load param
load it
load rec
load R

%move sizer output into save array
R(z).ens = z;
R(z).rec = p2k_cp;
R(z).filter_y = p2k_filteryear;
R(z).filters = p2k_filters;
R(z).sizer_x = sizer_x;
R(z).sizer_y = sizer_y;
R(z).sizer_z0 = sizer_z0;
R(z).sizer_z90 = sizer_z90;
R(z).sizer_zall = sizer_zall;

save('R','R')

%save figures to output file 

open sizer_sig0.fig
cd ../output/plots_example % Note: change save folder name here if needed
print(gcf,'-djpeg',['example_rec' int2str(z) '_sig0']); % Note: change save figure name here if needed

cd ../../working
open sizer_sig90.fig
cd ../output/plots_example % Note: change save folder name here if needed
print(gcf,'-djpeg',['example_rec' int2str(z) '_sig90']); % Note: change save figure name here if needed

close all
clear all

%reload files to allow loop to continue
cd ../../working
load it
load input
load crop

cd ../

end

close all
clear all

%save compiled sizer output to changepoints folder
cd working
load R
cd ../output

cp_output = R;
save('example_output','cp_output') % Note: change save file name here if needed

close all
clear all
cd ../

%% run changepoint calculation
%
% Uses sizer output to screen for onset of sustained, significant trends
%

cd output 
load example_output % Note: change save file name here if needed
cd ../

%parameters
max_year = 2015; %year to put if no recent significant change is observed
filter_clip = 11; %filter cut off (filters set in sizer parameters to screen between 5 years (filter 1) and 50 years (filter 46). Set at 11 so that only filters from 15y and higher are considered in determining median changepoint

%put data into temporary file
R = cp_output;

%run changeoint function
R_cp = sizer_changepoints(R, max_year, filter_clip); 

%save data back to file
cd output
cp_output = R_cp; %save updated datafile
clear R_cp

save('example_output','cp_output') % Note: change save file name here if needed

close all
clear all

load example_output % Note: change save file name here if needed

cd ../

