% Determine the best theorized linear fits to one set of data at a time
%
% Requires 'Statistics and Machine Learning Toolbox' with R2017a or later
%

clear all;
close all;

show_figs = 1;

load('mednoise_quadSVM.mat');       % to use in the classificatoin procedure

%% Generate Data to Fit
all_ages = sort(randi([18 94],1100,1));     % for generating random set of ages, note the ages ought to be sorted
ages = unique(all_ages);                    % removes ages that are repeated
age_span = 18:94;                           % for specifying the 77 ages used in the classification procedure

noise_level                 = .75;           % 0.75 based on calculations from Fjell data; constant across ages; 0.075 as a 'highly reduced noise' version; .4 for middle value

Curve_Type                  =  3;           % specify which type of data to simulate

critical_age = randi([45 75],1,1);
switch Curve_Type
    case 0
        data = all_ages*-.04 + randn(size(all_ages))*noise_level;
        critical_age = 0;
    case 1
        data = .001*(all_ages-critical_age).^2 + randn(size(all_ages))*noise_level;
    case 2
        data = -.8./(1+exp(-1*(all_ages-critical_age))) + randn(size(all_ages))*noise_level;
    case 3
        data = all_ages*-.04 + ((all_ages-critical_age)*-.04)./(1+exp(-1*(all_ages-critical_age))) + randn(size(all_ages))*noise_level;
end
       
data_to_fit = [all_ages, data]';

%% Fit a Smoothing Spline to the Data
bf_p = bic_smoothing_spline_range(data_to_fit(1,:)', data_to_fit(2,:)');  % determine the best fitting spline, constrained by BIC
output_vals = bf_p;
[f, goodness, output] = fit(data_to_fit(1,:)', data_to_fit(2,:)','smoothingspline','SmoothingParam',output_vals); % generate the best fitting spline
bf_spline = feval(f,ages)';

%% Analyze Spline for Critical Ages
[fx fxx] = differentiate(f,ages);
fxxx = [fxx(2:end)-fxx(1:(end-1));(fxx(end)-fxx(end-1))];

first_0_pts = nonzeros((abs(conv(sign(fx),[-1 1],'same'))>1).*(1:size(fx))');
spline_num_first_guesses = length(ages(first_0_pts));
[junk, index] = min(abs(fxx(first_0_pts)));
temp_spline_age_guess_first = ages(first_0_pts(index));

second_0_pts = nonzeros((abs(conv(sign(fxx),[-1 1],'same'))>1).*(1:size(fx))');
spline_num_second_guesses = length(ages(second_0_pts));
[junk, index] = max((fx(second_0_pts)-mean(fx)).^2);
temp_spline_age_guess_second = ages(second_0_pts(index));

third_0_pts = nonzeros((abs(conv(sign(fxxx),[-1 1],'same'))>1).*(1:size(fx))');
spline_num_third_guesses = length(ages(third_0_pts));
[junk,index] = max(abs(fxx));
temp_spline_age_guess_third = ages(index);

if isempty(temp_spline_age_guess_first)
    first_deriv_guess = 0;
else
    first_deriv_guess = temp_spline_age_guess_first;
end
if isempty(temp_spline_age_guess_second)
    second_deriv_guess = 0;
else
    second_deriv_guess = temp_spline_age_guess_second;
end
if isempty(temp_spline_age_guess_third)
    third_deriv_guess = 0;
else
    third_deriv_guess = temp_spline_age_guess_third;
end

spline_critical_age_guess_per_derivative = [first_deriv_guess,second_deriv_guess,third_deriv_guess];

%% Classify Spline
fjell_spline_curve = feval(f,age_span);
[fx_curve fxx_curve] = differentiate(f,age_span);

spline_to_categorize = (fx_curve'-mean(fx_curve))/std(fx_curve); % z-transform

spline_stats = [mean(fx_curve), std(fx_curve)]; % for 78th classification value

X = [spline_to_categorize,spline_stats(2)<.0005];

spline_category = mednoise_quadSVM.predictFcn(array2table(X));

if spline_category==categorical(4)
    spline_category=categorical(0);
    categorized_spline_guess=0;
elseif spline_category==categorical(3)
    categorized_spline_guess=third_deriv_guess;
elseif spline_category==categorical(2)
    categorized_spline_guess=second_deriv_guess;
elseif spline_category==categorical(1)
    if first_deriv_guess < 40
        categorized_spline_guess=third_deriv_guess;
    else
        categorized_spline_guess=first_deriv_guess;
    end
end

%% Graph Results
if show_figs
    figure;
    extreme_max_value = max(data_to_fit(2,:));
    extreme_min_value = min(data_to_fit(2,:));
    
    subplot(211);
    hold on;
    plot(data_to_fit(1,:), data_to_fit(2,:),'ko');
    h=plot(ages, bf_spline, 'linewidth', 3); set(h,'Color',[.5 .5 .5]);
    plot([categorized_spline_guess,categorized_spline_guess],[extreme_min_value,extreme_max_value],'k--', 'linewidth', 2)
    axis([0, 100, extreme_min_value, extreme_max_value]);
    ylabel('Data Points');
    
    extreme_max_value = max(fx);
    extreme_min_value = min(fx);
    
    subplot(212);
    hold on;
    h=plot(ages, fx, 'linewidth', 3); set(h,'Color',[.5 .5 .5]);
    plot([categorized_spline_guess,categorized_spline_guess],[extreme_min_value,extreme_max_value],'k--', 'linewidth', 2)
    axis([0, 100, extreme_min_value, extreme_max_value]);
    xlabel('Age');
    ylabel('Maturation Rate')
end

%% Return Results
results = [ Curve_Type, double(spline_category);
            critical_age, categorized_spline_guess];
display(sprintf('\t\t\t Actual   Guess \n Curve Type \t%d\t\t%d\n Critical Age \t%d\t\t%d',results(1),results(3),results(2),results(4)));