%Script prova approssimazione minimi quadrati
%Scritto da Pavan Andrea - 08/07/2022
clear;
clc;


%% minimi quadrati 1D
clear;
clc;
fprintf('CASO 1D\n');

%dati
%x = [-1 0 1 1.5 2 0.3 0.2 0.2];
x = -2 + (5-(-2))*rand(1,50);
x = [x,x];
u = 0.2*x.^2 + 0.6*(rand(1,length(x))-0.5);

%calcolo coefficienti polinomio
n = length(x);
A = [n, sum(x), sum(x.^2);
    sum(x), sum(x.^2), sum(x.^3);
    sum(x.^2), sum(x.^3), sum(x.^4)];
b = [sum(u); sum(u.*x); sum(u.*x.*x)];
c = A\b;
% [Q,R] = qr(A);
% c = R\(Q\b);
uapprox = @(xeval) c(1) + c(2)*xeval + c(3)*xeval.^2;
fprintf('Condizionamento matrice A: %f\n',cond(A));
fprintf('\n');

%plot
plot(linspace(min(x),max(x)),uapprox(linspace(min(x),max(x))),'b-');
hold on;
plot(x,u,'ko');
title('Esempio minimi quadrati 1D');
xlabel('x');
ylabel('u(x)');
axis equal;
hold off;


%% minimi quadrati 2D
clear;
clc;
fprintf('CASO 2D\n');

%dati
%x = [1.2 1 -1 -1 0 -0.2 0.5 1 1];
%y = [1.2 -1 1 -1 0 0.2 -0.3 1 1.5];
x = -2 + (2-(-2))*rand(1,50);
y = -2 + (2-(-2))*rand(1,50);
u = 0.2*x.^2 + 0.5*y.^2 + 0.1*(rand(1,length(x))-0.5);

%calcolo coefficienti polinomio
T0 = clock;
n = length(x);
A = [n, sum(x), sum(y), sum(x.^2), sum(y.^2), sum(x.*y);
    sum(x), sum(x.^2), sum(y.*x), sum(x.^3), sum(y.^2.*x), sum(x.^2.*y);
    sum(y), sum(x.*y), sum(y.^2), sum(x.^2.*y), sum(y.^3), sum(x.*y.^2);
    sum(x.^2), sum(x.^3), sum(y.*x.^2), sum(x.^4), sum(y.^2.*x.^2), sum(x.^3.*y);
    sum(y.^2), sum(x.*y.^2), sum(y.^3), sum(x.^2.*y.^2), sum(y.^4), sum(x.*y.^3);
    sum(x.*y), sum(x.^2.*y), sum(y.^2.*x), sum(x.^3.*y), sum(y.^3.*x), sum(x.^2.*y.^2)];
b = [sum(u); sum(u.*x); sum(u.*y); sum(u.*x.^2); sum(u.*y.^2); sum(u.*x.*y)];
c = A\b;
timelinsys = etime(clock,T0);
fprintf('Tempo per risolvere Ac=b: %.3f s\n',timelinsys);

%prova mediante inversione matrice A
T0 = clock;
invA = inv(A);
timeinv = etime(clock,T0);
fprintf('Tempo per invertire A: %.3f s\n',timeinv);

uapprox = @(xeval,yeval) c(1) + c(2)*xeval + c(3)*yeval + c(4)*xeval.^2 + c(5)*yeval.^2 + c(6)*xeval.*yeval;
fprintf('Condizionamento matrice A: %f\n',cond(A));
fprintf('\n');

%plot
[xplot,yplot] = meshgrid(min(x):((max(x)-min(x))/50):max(x),min(y):((max(y)-min(y))/50):max(y));
figure();
surf(xplot,yplot,uapprox(xplot,yplot));
hold on;
plot3(x,y,u,'w.');
title('Esempio minimi quadrati 2D');
xlabel('x');
ylabel('y');
zlabel('u(x,y)');
axis equal;
hold off;


%% minimi quadrati 1D - derivata imposta
clear;
clc;
fprintf('CASO 1D\n');

%dati
%x = [-1 0 1 1.5 2 0.3 0.2 0.2];
dydxg = -1;      %valore derivata imposta in x=0
x = -2 + (5-(-2))*rand(1,50);
x = [x,x];
u = 0.2*x.^2 + dydxg*x + 0.2*(rand(1,length(x))-0.5);

%calcolo coefficienti polinomio
n = length(x);
A = [n, sum(x.^2);
    sum(x.^2), sum(x.^4)];
b = [sum(u)-dydxg*sum(x);
    sum(u.*x.*x)-dydxg*sum(x.^3)];
c = A\b;
% [Q,R] = qr(A);
% c = R\(Q\b);
uapprox = @(xeval) c(1) + dydxg*xeval + c(2)*xeval.^2;
fprintf('Condizionamento matrice A: %f\n',cond(A));
fprintf('\n');

%plot
plot(linspace(min(x),max(x)),uapprox(linspace(min(x),max(x))),'b-');
hold on;
plot(x,u,'ko');
title('Esempio minimi quadrati 1D');
xlabel('x');
ylabel('u(x)');
axis equal;
hold off;


%% minimi quadrati 2D
clear;
clc;
fprintf('CASO 2D\n');

%dati
%x = [1.2 1 -1 -1 0 -0.2 0.5 1 1];
%y = [1.2 -1 1 -1 0 0.2 -0.3 1 1.5];
x = -2 + (2-(-2))*rand(1,50);
y = -2 + (2-(-2))*rand(1,50);
u = 0.2*x.^2 + 0.5*y.^2 + 0.1*(rand(1,length(x))-0.5);

%calcolo coefficienti polinomio
T0 = clock;
n = length(x);
A = [n, sum(x), sum(y), sum(x.^2), sum(y.^2), sum(x.*y);
    sum(x), sum(x.^2), sum(y.*x), sum(x.^3), sum(y.^2.*x), sum(x.^2.*y);
    sum(y), sum(x.*y), sum(y.^2), sum(x.^2.*y), sum(y.^3), sum(x.*y.^2);
    sum(x.^2), sum(x.^3), sum(y.*x.^2), sum(x.^4), sum(y.^2.*x.^2), sum(x.^3.*y);
    sum(y.^2), sum(x.*y.^2), sum(y.^3), sum(x.^2.*y.^2), sum(y.^4), sum(x.*y.^3);
    sum(x.*y), sum(x.^2.*y), sum(y.^2.*x), sum(x.^3.*y), sum(y.^3.*x), sum(x.^2.*y.^2)];
b = [sum(u); sum(u.*x); sum(u.*y); sum(u.*x.^2); sum(u.*y.^2); sum(u.*x.*y)];
c = A\b;
timelinsys = etime(clock,T0);
fprintf('Tempo per risolvere Ac=b: %.3f s\n',timelinsys);

%prova mediante inversione matrice A
T0 = clock;
invA = inv(A);
timeinv = etime(clock,T0);
fprintf('Tempo per invertire A: %.3f s\n',timeinv);

uapprox = @(xeval,yeval) c(1) + c(2)*xeval + c(3)*yeval + c(4)*xeval.^2 + c(5)*yeval.^2 + c(6)*xeval.*yeval;
fprintf('Condizionamento matrice A: %f\n',cond(A));
fprintf('\n');

%plot
[xplot,yplot] = meshgrid(min(x):((max(x)-min(x))/50):max(x),min(y):((max(y)-min(y))/50):max(y));
figure();
surf(xplot,yplot,uapprox(xplot,yplot));
hold on;
plot3(x,y,u,'w.');
title('Esempio minimi quadrati 2D');
xlabel('x');
ylabel('y');
zlabel('u(x,y)');
axis equal;
hold off;

