function out = get_opt_poly(X, Y, poly_range, win)
%
% This function was built to optimize the polynomial order used for 
% f(P) & k2_f(P) polynomials derived for MBARI pH calibrations in the
% Matlab pH calibration GUI
%
% INPUTS:
%   X          - Independant variable (Pressure array in calibration data)
%   Y          - Dependent variable (k2 or f(P) arrays derived from
%                calibration data)
%   poly_range - an integer array of polynomial orders (i.e. [0 1 2 3 4 5 6])
%   win        - filter window size for median filter. Used to get
%                anomalies to get error limits for BIC estimate
%                k2   default = 5
%                f(P) default = 7

% ***********************************************************************
% TESTING
% load('C:\temp\k2fP_test_data.mat'); %TESTING

% K2 f(P)
%X = d.BinK2_data(:,1);
% X = d.BinK2_data(:,1);
% Y = d.BinK2_data(:,2);
% poly_range = 0:6;
% win = 5;

% f(P)
% X = d.binfP_data(:,1);
% Y = d.binfP_data(:,2);
% poly_range = 0:8;
% win = 7;
% ***********************************************************************

% SET DEFAULT VALUES IF NO DATA
out.hdr  = {'poly order' 'NSSE' 'BIC'};
out.data = [NaN NaN NaN];
out.BIC  = NaN;
out.poly_order = NaN;
out.X = X;
out.Y = Y;
out.range = NaN;
out.tstr = sprintf('No data to fit');

% Check that data exist to fit in data subset
if isempty(X) || isempty(Y) % no data to fit
    fprintf('At least one input is empty!\n')
    return
end

% Check polynomial integration range vs data points
if max(poly_range) >= max(size(X))
    fprintf(['Max poly order (%0.0f) >= than # of data points (%0.0f).', ...
        'poly order will be truncated to max # of data points-1.\m'], ...
        max(poly_range), max(size(X)));
    poly_range(poly_range >  max(size(X))) = [];
end

% USE A MEDIAN FILTER TO ASSES MEASUREMENT NOISE (ERROR LIMIT FOR BIC)
anom = Y - movmedian(Y,win);
error_limit = std(anom);

% NOW STEP THROUGH POLY_ORDERS & CALC BIC
rpoly = max(size(poly_range));
ppY = ones(size(Y,1), rpoly)*NaN; % A matrix of polynomil fits w increasing order
for ct = 1:max(size(poly_range))
    pp = polyfit(X,Y, poly_range(ct));
    ppY(:,ct) = polyval(pp,X);
end

N     = max(size(Y)); % sample ct
anoms = Y - ppY;
SSE   = sum(anoms.^2, 1,'omitnan');
KK    = (poly_range + 1) .*2;
%BIC = log(SSE./N) + KK.*log(N)./N;

% error_limit.^2 is an additional non standard term add to the tradtional BIC
% equation to force lower polynomial selection
BIC = log(SSE./N + error_limit.^2) + KK.*log(N)./N;
% BICa = log(SSE./N + error_limit.^2) ;
% BICb = KK.*log(N)./N;

% https://surveillance.cancer.gov/help/joinpoint/tech-help/frequently-asked-questions/bayesian-information-criterion-bic-details
% The equation for computing the BIC for a k-joinpoint model is:
% BIC(k) = ln{SSE(k)/#Obs} + {#Parm(k) /#Obs} * ln(#Obs),
% where SSE(k) is the sum of squared errors of the k-joinpoint regression model, 
% #Parm(k)=2*(k+1) is the number of parameters of the k-joinpoint model 
% and #Obs is the number of observations.
% The k-joinpoint model with the minimum value of BIC(k) is selected as the final model.

% ASSES BIC (Round to 0.X) - below that is really just noise
round_bic = round(BIC,1);
t1 = round_bic == min(round_bic,[],2,'omitnan');
min_poly  = min(poly_range(t1));
idx       = find(poly_range == min_poly, 1);

%[~,idx] = min(BIC,[],2,'omitnan');

sY = sort(Y);
%rY = (sY(end) - sY(1))./ mean(sY,'omitnan')*100; % for title strings

out.hdr  = {'poly order' 'NSSE' 'BIC'};
out.data = [poly_range; SSE./N; BIC]';
out.BIC  = BIC(idx);
out.poly_order = poly_range(idx);
out.X = X;
out.Y = Y;
%out.range = abs(rY);

out.range = sY(end) - sY(1);


out.tstr = sprintf('opt. poly order = %0.0f   (BIC = %0.1f)',...
    poly_range(idx),BIC(idx));





