%lab2.m
graphics_toolkit("gnuplot")
echo on
disp("functions")
%functions write a simple function with 2 outputs and 2 parameters input
%returning the value fo f(x)=x^2*a a parameter default value 1
function [y,dy]=sq(x,a=1)
   y=x*x*a;
   dy=2*x*a;
   if(nargin<2)
    disp("Warning a not defined- a =1 by default");
   end 
   if(nargout<2)
    disp("Warning df lost - not called for ");
   end 
end   

[y,dy]=sq(2,3)

[y,dy]=sq(2)

y=sq(2)

y=sq(3,5)


%function with a parameter which is a function handle 
function y=rect(f,a,b)
 y=(b-a)*f(0.5*(a+b));
end

rect(@sin,-1,1)
rect(@(x) 1+x*x,0,1)
 
%global variable
global aa
aa=3
function y=f(x)
  global aa
  y=sin(aa*x);
end

rect(@f,0,1)  
aa=4

rect(@f,0,1)  

%linear algebra
%solving linear system
A=[1 2;3 4];f=A(:,1);x=A\f
A=[1 2; 1 2];f=A(:,1);x=A\f
%LLSP
A=[1 2;3 4;5 6];f=A(:,1);x=A\f
f=randn(3,1);x=A\f
norm(A*x-f)

%interpolation Lagrange - compute p polynomial such that 
%p(x(k))=y(k) k=0,..,n solving vandermonde system
x=[0 1 2 4]'
y=[0 0 1 0]'
A=vander(x);
p=A\y;
norm(polyval(p,x)-y,"inf")
z=linspace(0,4,200);
plot(z,polyval(p,z),x,y,"+");pause(1);

%taking more and more equaidistant points on [0,1] and interpolating

for k=1:10,
  n=2+5*k;
  x=linspace(0,1,n+1)';
  y=sin(2*x);
  A=vander(x);
  p=A\y;
  e(k)=norm(polyval(p,x)-y,"inf");
end  


%linear regression
%write a function whcih for a gvien two sets of m points in two vectors x and y
%find [a;b] such that they minimize \sum_k |ax(k)+b-y(k)|^2
%The function should print out a warning if there is no uniqueness...
%HINT: one may solve the LLSP obtaining first QR factors of A; then chack if numerically the A has 
%max rank, and then use this factor to compute the solution (actually backslash does that without checking 
%the max rank of A
function [a,b]=linregression(x,y)
%INPUT x,y vectors of length m
%OUTPUT a,b  - the solution of the linear regression problem  
 if(isrow(x)) x=x';end
 if(isrow(y)) y=y';end
 A=[x,ones(size(x))];
 s=A\y;
 a=s(1);b=s(2);
 if(max(x)-min(x) < 1e-12) 
  disp("Warning: the problem may have no unique solution!");
 end
endfunction

x=randn(10,1);y=34*x-17;
[a,b]=linregression(x,y)
x=2*ones(10,1);
y=23*x-7 +1e-1*randn(size(x));
[a,b]=linregression(x,y)
