%approssimazione ai minimi quadrati

% n : grado del polinomio
% m : numero di nodi
% a : vettore coefficienti polinomio
% x : nodi 
% y : f(x)
% xx: punti da valutare

clear all
close all
clc

n=1;
x=linspace(-1,1,5);
y=x.^4
[C] = SIST(x,n);
[Q,R] = QRC(C);
[a] = SolveTriuQr(Q,R,y);
xx=[1; 10; 40];
[h] = POLYVAL(xx,a)

f=@(xx) POLYVAL(xx,a);
fplot(f)
hold on
plot(x, y, '*')

function [C] = SIST(x,n)
% C:per definizione del sistema
m=length(x);
C=ones(m,n);

for i=1:m
    for j=2:n+1
        C(i,j)=x(i)*C(i,j-1);
    end
end

end

function [Q,R] = QRC(C)
%DECOMPOSIZIONE QR DI C, Qortog simm e R trapez sup
% Alla fine del ciclo Q in realtà è Q'

[l,q]=size(C);
Q=eye(l);
R=C;
for i=1:q
    v=C(i:l,i);
    u=zeros(l,1);
    if round(v(1),8)==0
        v(1)=0;
        u(i:l)= v +norm(v)*eye(l-i+1,1);
    else
    u(i:l)=v+sign(v(1))*norm(v)*eye(l-i+1,1);
    end

    Qi= eye(l)-2*(u*u')/(u'*u);
    Q=Qi*Q;
    C=Qi*C;
end
R=C;

end

function [a] = SolveTriuQr(Q,R,b)
% Risoluzione del sistema triangolare
[~,n]=size(R);
d=Q*b';
d1=d(1:n);
R1=R(1:n,1:n);
a=eye(n,1);
a(n)=d1(n)/R1(n,n);

for i=n-1:-1:1
    a(i) = (d1(i) - R1(i, i+1:n)*a(i+1:n))/(R1(i,i));
end
end

function [P] = POLYVAL(xx,a)
%valutazione con i minimi quadrati
%base convenzionale {1, xx, xx^2...,xx^n}
n=length(a)-1;
P=0;

for i=0:n
     P= P + a(i+1)*xx.^i;
end

end