inflectionPointTestLinearMdl

%Testing code that calculates inflection points for a synthetic fitlm (linear) model

% 1. Generate synthetic data with known inflection points
x = linspace(-10, 10, 500)'; % Column vector as input for the model
y = exp(1.5 * x); % Define a function with inflection points

% 2. Fit a polynomial regression model to the data
degree = 7; % Choose a polynomial degree that can capture the inflection points
mdl = fitlm(x, y, sprintf('poly%d', degree)); % Fit the model

% 3. Predict y values using the model
x_highres = linspace(-10, 10, 5000)'; % Higher resolution for better inflection detection
y_pred = predict(mdl, x_highres); % Predicted values

% 4. Apply inflection point detection
% Interpolation and smoothing
y_interp = interp1(x_highres, y_pred, x_highres, 'pchip'); % Interpolation
y_smooth = smoothdata(y_interp, 'gaussian', 50); % Gaussian smoothing

% 5. Derivatives
dy_dx = gradient(y_smooth) ./ gradient(x_highres); % First derivative
d2y_dx2 = gradient(dy_dx) ./ gradient(x_highres); % Second derivative

% 6. Detect zero-crossings in the second derivative
sign_changes = find(diff(sign(d2y_dx2)) ~= 0); % Detect strict sign changes

% 7. Initialize arrays for inflection points
inflection_points_x = [];
inflection_points_y = [];

for i = 1:length(sign_changes)
    inflection_x = x_highres(sign_changes(i));
    inflection_y = y_smooth(sign_changes(i));
    
    % Ignore points near the edges, e.g., within 5% of each end of the x range
    if inflection_x > x_highres(round(0.05 * length(x_highres))) && ...
       inflection_x < x_highres(round(0.95 * length(x_highres)))
    
        % Ensure unique inflection points based on proximity
        if isempty(inflection_points_x) || all(abs(inflection_x - inflection_points_x) > 1e-3)
            inflection_points_x = [inflection_points_x, inflection_x];
            inflection_points_y = [inflection_points_y, inflection_y];
            
            % Display each unique inflection point
            disp(['Inflection point: (x, y) = (' num2str(inflection_x) ', ' num2str(inflection_y) ')']);
        end
    end
end

% 8. If no inflection points are found, display a message
if isempty(inflection_points_x)
    disp('No inflection points found.');
else
    % Plotting the results to visualize inflection points
    figure;
    plot(x_highres, y_pred, 'b-', 'LineWidth', 1.5); % Plot predicted values
    hold on;
    scatter(inflection_points_x, inflection_points_y, 100, 'r', 'filled'); % Mark inflection points
    title('Predicted Model with Inflection Points');
    xlabel('x');
    ylabel('y');
    legend('Model Prediction', 'Inflection Points');
    hold off;
end