clear all;

% Paremeters
N = 100000;
S = 25;

% Generate 201 numbers grid from -1 to 1
u = linspace(-1, 1, 201);

% Generate N random numbers using the given function
u_k = 2 * rand(1,N) - 1;

% Initialize lists
ai_List = zeros(1, N);
bi_List = zeros(1, N);
ai_Mean_List = zeros(1, S);
bi_Mean_List = zeros(1, S);

u_dash_List = zeros(1, length(u));
u_dash_num = zeros(1,length(u));
u_dash_den = zeros(1,length(u));

% Generate legendre values for u_k
P_k = legendre(S,u_k,'norm');

for s = 1:S
    for k = 1:N
        % Trigonometric fi as basis function 
        ai_List(k) = fi(s,u_k(k));
        y_k = traingle_func(u_k(k)) + (rand() - 0.5);
        bi_List(k) = y_k * fi(s,u_k(k));

        % Legendre as basis function
        %ai_List(k) = P_k(s+1,k);
        %y_k = traingle_func(u_k(k)) + (rand() - 0.5);
        %bi_List(k) = y_k * P_k(s+1,k);
    end
    ai_Mean_List(s) = mean(ai_List);
    bi_Mean_List(s) = mean(bi_List);
end

% Generate legendre values for u
P = legendre(S,u,'norm');

for x = 1:length(u)
    for i = 1:S
        % Trigonometric fi as basis function 
        u_dash_num(x) = u_dash_num(x) +  bi_Mean_List(i) * fi(i,u(x));
        u_dash_den(x) = u_dash_den(x) +  ai_Mean_List(i) * fi(i,u(x));
        
        % Legendre as basis function
        %u_dash_num(x) = u_dash_num(x) +  bi_Mean_List(i) * P(i+1,x);
        %u_dash_den(x) = u_dash_den(x) +  ai_Mean_List(i) * P(i+1,x);
    end
    u_dash_List(x) = u_dash_num(x)/u_dash_den(x);
end

% Plot the theoretical and estimated Regressions
figure;
plot(u, arrayfun(@traingle_func, u), 'LineWidth', 4, 'DisplayName', 'Theoretical Regression', 'Color', [0, 255, 255] / 255);
hold on;
plot(u, u_dash_List,'LineWidth', 2, 'DisplayName', 'Estimated Regression', 'Color', [255, 105, 180] / 255);
grid on;
legend('show',Location='best');
title(['Estimated and Theoretical Regressions for S=', num2str(S)]);
hold off;

% Function Definition
function fi_result = fi(i, x)
    if i == 1
        fi_result = 1 / sqrt(2);
    elseif mod(i, 2) == 0
        fi_result = sin(i/2 * pi * x);
    else
        fi_result = cos((i-1)/2 * pi * x);
    end
end

function y = traingle_func(x)
    if x >= -1 && x < -0.5
        y = -2*x-2;
    elseif x >= -0.5 && x <= 0.5
        y = 2*x;
    elseif x >= 0.5 && x <= 1
        y = -2*x+2;
    else
        y = 0;
    end
end