// Konvencija je da skripte dobivaju ekstenziju .sce, a 
// funkcije .sci
///////////////////////////////////////////////////////////
// Metoda najmanjih kvadrata.
// Zadanu funkciju aproksimiramo polinomom stupnja n
// u smislu najmanjih kvadrata. Radi se o diskretnom
// problemu najbolje aproksimacije.
///////////////////////////////////////////////////////////
//
// skripta se izvršava naredbom:
//   exec("ime_skripte");

//n=x_dialog("Stupanj aproksimativnog polinoma : ");
//N=x_dialog("Broj tocaka na itervalu [0,4] : ");

clear;

n=input("Stupanj aproksimativnog polinoma : ");  // unos podataka
N=input("Broj tocaka na itervalu [0,4] : ");

if n > N  then
   error("Broj tocaka mora biti veci od stupnja polinoma.");
end

//  Točke segementa [0,4] u kojima su zadane funkcijske vrijednosti
x=linspace(0,4,N)';

// Funkcijske vrijednosti ( f(x) = cos(x^2) )
y=cos(x.^2);

// Konstrukcija matrice (Vandermondeova matrica)

A=ones(N,n+1);

for i=1:n
  A(:,i+1) = A(:,i).*x;
// ili 
// A(:,i+1)=x.^i;
end

// Izračunaj rang matice
Arang=rank(A);

// U smisli najmanjih kvadrata riješava A coef = y, tj
//  || A coef -y ||^2 --> min

coef=lsq(A,y);

// Konstrukcija polinoma pomoću koeficijenata. Varijabla se zove s.
p=poly(coef,"s","coef");

// Novo polje x-ova za crtanje rješenja i polinoma

xx=(0:0.01:4)';
yy=cos(xx.^2);
xbasc()    // očisti ekran
plot2d(xx,[yy horner(p,xx)],leg="cos(x^2)@lsq-aproksimacija")
xtitle("Aproksimacija u smislu najmanjih kvadrata","x","y")

printf("Greska aproksimacije = || A coef - y || = %f\n",norm(A*coef-y));
printf("Rang matrice A = %d\n",Arang);

//  Napomena. Sa 128  točaka polinomi stupnja n=11 i n=12 daju najbolju
//            aproksimaciju; greška je ok 1.9


//  Riješimo sada sustav eksplicitno koristeći QR faktorizaciju, 
//  umjesto funkcije lsq()

printf("\nRjesenje metodom najmanjih kvadrata putem QR faktorizacije (funkcija qr)\n");
pause

[Q,R]=qr(A);  //  QR-faktorizacija (potpuna)

d=Q'*y;   //  

coef1 = R(1:n+1,:)\d(1:n+1);   //  Rhat coef1 = dhat  (hat =  educiran)
//norm(coef-coef1) koeficijenti se razlikuju jer se matrica R razlikuje
p1=poly(coef1,"s","coef");
xx=(0:0.05:4)';
plot2d(xx,horner(p1,xx),style=-2)
printf("Greska aproksimacije = || d(n+2:N) || = %f\n",norm(d(n+2:N)));
printf("Greska aproksimacije = || A coef1 - y || = %f\n",norm(A*coef1-y));


// Riješimo isti problem SVD faktorizacijom
printf("Rjesenje metodom najmanjih kvadrata putem SVD faktorizacije (funkcija svd)\n");
pause

[U,S,V]=svd(A);  // size(U) = N x N, size(S)= N x (n+1),  size(V) = (n+1) x (n+1)

d2 = U'*y;  // N x 1
coef2 = S(1:n+1,:)\d2(1:n+1);  // S coef2 = U' y  (prvih n+1 jednadžbi)
coef2=V*coef2;

p2=poly(coef2,"s","coef");
plot2d(xx,horner(p2,xx),style=-4)

printf("Greska aproksimacije = || d2(n+2:N) || = %f\n",norm(d2(n+2:N)));
printf("Greska aproksimacije = || A coef2 - y || = %f\n",norm(A*coef2-y));