% Plot function
function [] = draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true)  
    figure;
    plot(x_plot, f, 'b');
    hold on;
    plot(x, y, 'ro');
    hold on;
    plot(x_pred, y_pred, 'bo', 'LineWidth', 1.5);
    hold on;
    plot(x_pred, y_true, 'ro', 'LineWidth', 1.5);
    hold on;
    plot(x, p_y, 'bo');
    title('Curve fitting');
    xlabel('Days');
    ylabel('Drug released (ug)');
end

% Residual calculate function
function [] = res(x, p_y, y)
    residual = zeros(1, length(y));
    for i = 1:length(y)
        residual(i) = p_y(i) - y(i);
    end
    fprintf("Residual: ");
    disp(residual);
    sum = 0;
    for i = 1:length(residual)
        sum = sum + residual(i) * residual(i);
    end
    fprintf("Sum of residual squares: %.2f\n", sum);

    figure;
    bar(x, residual);
    title('Residual values');
    xlabel('x');
    ylabel('y');
end

% Print function
function [] = print(p, p_y)    
    fprintf("Parameters:\n");
    disp(p);
    fprintf("Model values:\n ");
    disp(p_y);
end

% Read the data
x = [];
y = [];
fin = fopen("data.txt", 'r');
while ~feof(fin)
    line = fgetl(fin);
    cells = strsplit(line, ',');
    x = [x, str2double(cells{1})];
    y = [y, str2double(cells{2})];
end

% Curve fitting
x_plot = linspace(min(x) - 30, max(x) + 30, 100);
x_pred = [12 16 21];
y_true = [28.4 28.5 29.5];
for j = 1:3
    if j == 1
        p = polyfit(x, y, 1); % Estimate function parameters
        p_y = polyval(p, x); % Estimate function value
        y_pred = polyval(p, x_pred);
        f = p(1)*x_plot + p(2); % Estimate function
        draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true);
        print(p, p_y);
        res(x, p_y, y);
    elseif j == 2
        p = polyfit(x, y, 2); % Estimate function parameters
        p_y = polyval(p, x); % Estimate function value
        y_pred = polyval(p, x_pred);
        f = p(1)*x_plot.^2 + p(2)*x_plot + p(3); % Estimate function
        draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true);
        print(p, p_y);
        res(x, p_y, y);
    else
        p = polyfit(x, y, 3); % Estimate function parameters
        p_y = polyval(p, x); % Estimate function value
        y_pred = polyval(p, x_pred);
        f = p(1)*x_plot.^3 + p(2)*x_plot.^2 + p(3)*x_plot + p(4); % Estimate function
        draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true);
        print(p, p_y);
        res(x, p_y, y);
    end
end

% Fminsearch
% Initial guess
x0 = [10, 2];

% The square of sum of residual
func = @(params) sum((params(1)*(1 - exp(-params(2)*x)) - y).^2);

% fminsearch implement
optimal_params = fminsearch(func, x0);

% Create the estimate function with the result parameter by the fminsearch
optimal_func = optimal_params(1)*(1 - exp(-optimal_params(2)*x_plot));
optimal_func_pred = optimal_params(1)*(1 - exp(-optimal_params(2)*x_pred));

% Calculate the values with respect to the result parameter
p_y = optimal_params(1)*(1 - exp(-optimal_params(2)*x));

%Plot the graph
draw(x, x_plot, y, optimal_func, p_y, x_pred, optimal_func_pred, y_true);

%Print out value and the residual
print(optimal_params, p_y);
res(x, p_y, y);

% Terminate section
fclose(fin);