% curvefit_demo.m

% David Pan, UAH

% Example curve: y = 1 + 2x + 3x^2

N = 4
% Generate 4 data points for training
rng(1);
x = 10*rand(1, N);

Z = 1 + 2*x + 3*x.^2;  % Target values

% Formulate the input data matrix
Fx = zeros(3,N);
for i = 1:N
    Fx(:,i) = [1, x(i), x(i)^2];
end

W = Z / Fx
W2 =  Z * pinv(Fx)


[p,s] = polyfit(x, Z, 2);   % notice the reversed order
wrev(p)   % Show weights from low to high orders


% Now with noise added
N = 4
rng(1);
x = 10*rand(1, N);

rng(1);
Z = 1 + 2*x + 3*x.^2 + randn(1, N);

% Formulat the input data matrix
Fx = zeros(3,N);
for i = 1:N
    Fx(:,i) = [1, x(i), x(i)^2];
end

W = Z / Fx

% Mimic the Matlab implemenation
% The Vondermonde matrix 
V = fliplr(Fx');    % Transpose followed by left to right flip
[Q,R] = qr(V,0);   % Economy-size QR decomposition
% Orthornormal columns
Q'*Q

p2 = mldivide(R, Q'*Z');
p2'
% Compared wit the structure returned by polyfit( )
[p,s] = polyfit(x, Z, 2);
p
s.R
s.normr
% Compare the norm of the residues

normr2 = norm(W*Fx - Z)
normr3 = norm(V*p2 - Z')

% With more training data with much worse noise added
N = 100;
rng(1);
x = 10*rand(1,N);
rng(1);
Z = 1 + 2*x + 3*x.^2 + 20*randn(1, N);

% Export the dataset to a csv file (to be imported by sklearn)
data_label = horzcat(x', Z');
writematrix(data_label, 'regression_dataset.csv');


Fx = zeros(3,N);
for i = 1:N
    Fx(:,i) = [1, x(i), x(i)^2];
end

scatter(x,Z);grid

W = Z / Fx
W2 =  Z * pinv(Fx)
p = polyfit(x, Z, 2);   % notice the reversed order
wrev(p) 


hold on;
xx = min(x):0.01:max(x);
plot(xx, W(1) + W(2)*xx + W(3)*xx.^2);

[p,s] = polyfit(x, Z, 2);
p
s.R
s.normr

% Compared with the added noise norm
rng(1);
noise = 20*randn(1, N);
norm(noise)

% SVD method as a comparison with sklearn
[A,S,B]= svd(Fx','econ');
S
p_svd = B*S^(-1)*A'*Z'