%% changepoint screening function
%
% This file is called from within the "calculate_changepoints_usingsizer.m"
% routine. It uses the output obtained from sizer analysis to screen for
% the starting point of sustained, significant climate trends
%
% 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 further documentation.
%
%% function
function [R] = sizer_changepoints(R, max_year, filter_clip)

for j = 1:length(R);

%load sizer output data
[x,y,z0,z90,zall] = deal(R(j).sizer_x,R(j).sizer_y,R(j).sizer_z0,R(j).sizer_z90,R(j).sizer_zall);

%find points of change in trend significance
z90_ex = z90(:,2:end) - z90(:,1:end-1);

%adjust for any changes from NaNs (flagged by -999)
z90_ex(find(z90_ex == 1000)) = -2;
z90_ex(find(z90_ex == 1003)) = 1;

%screen matrices to assess onset of most recent sustained and
%significant trend
%Assess for recent warming or cooling trend

    for i = 1:size(z90,1); %screen across each filter length
        
        %step 1: determine if most recent sig trend is warming or cooling
        %cooling = 4, warming = 1
        sigc = find(zall(i,:) == 4,1,'last');
        sigw = find(zall(i,:) == 1,1,'last');
        
        if isempty(sigc) == 1
            ysigc = x(1,1);
        else
            ysigc = x(1,sigc);
        end
        
        
        if isempty(sigw) == 1
            ysigw = x(1,1);
        else
            ysigw = x(1,sigw);
        end
        
        %assign classification
        if ysigc < ysigw
            class = 1; %warming
        elseif ysigc > ysigw
            class = -1; %cooling
        else
            class = 0; %neutral
        end
            
        %Step 2; determine if sign of most recent significant trend is
        %continuous to present
        if class == -1
            sigy_w = max_year;
            
            w = find(zall(i,:) == 2,1,'last');
                if isempty(w) == 1
                yw = x(1,1);
                else
                yw = x(1,w);
                end
                
            if yw > ysigc
               sigy = max_year; %sigy = x(1,end); to assign to last year instead
               sigy_c = sigy;
            else
               sigc = find(z90_ex(i,:) == 1,1,'last');
               if isempty(sigc) == 1
               sigy = x(1,1);
               else
               sigy = x(1,sigc+1);
               end
               sigy_c = sigy;
            end
            
        elseif class == 1
            sigy_c = max_year;
            
            c = find(zall(i,:) == 3,1,'last');
                if isempty(c) == 1
                yc = x(1,1);
                else
                yc = x(1,c);
                end
                
            if yc > ysigw
               sigy = max_year;
               sigy_w = sigy;
            else
               sigw = find(z90_ex(i,:) == -2,1,'last');
               if isempty(sigw) == 1
               sigy = x(1,1);
               else
               sigy = x(1,sigw+1);
               end
               sigy_w = sigy;
            end
            
        else %class = 0
            sigy = max_year;
            sigy_w = sigy;
            sigy_c = sigy;
        end
        
        sig_year(i,1) = sigy;
        sig_class(i,1) = class;
        sigwarm_year(i,1) = sigy_w;
        sigcool_year(i,1) = sigy_c;
        
    end

%calculate medians across selected filter widths (11-50y), and save to array
 sigwarm_median = median(sigwarm_year(filter_clip:end,1));
 sigcool_median = median(sigcool_year(filter_clip:end,1));
 
 
 R(j).cp.sig_year = sig_year;
 R(j).cp.sig_class = sig_class;
 R(j).cp.sigwarm_year = sigwarm_year;
 R(j).cp.sigcool_year = sigcool_year;
 R(j).cp.sigwarm_median = sigwarm_median;
 R(j).cp.sigcool_median = sigcool_median;   
 R(j).cp.filter_clip = filter_clip;
 R(j).cp.max_year = max_year;   
end


end

