%% Data fitting examples
%
%% Polynomial interpolation
% Matlab function polyfit calculates coefficients and polyval evaluates at points.
%
%% Basic example

xd=[-1 -.96 -.65 0.1 0.4 1];
yd=[-1 -.1512 .386 .4802 .8838 1];

N=length(xd);
c=polyfit(xd,yd,N-1);   % N datapts -> deg(poly) = N-1

% Evaluate poly at many x-points.
x=linspace(-2,2);
y=polyval(c,x);
plot(xd,yd,'*',x,y);grid on;shg
xlim([-1.1 1.1])
title('Polynomial interpolation')

%
%% Problems: Runge's phenomenon
% Runge's phenomenon is a problem of oscillation at the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree over a set of equispaced interpolation points. Runge 1901: https://en.wikipedia.org/wiki/Runge%27s_phenomenon
%
% Consider
% g(x) = \frac{1}{1+25 x^2}
% at N equally spaced on [-1,1].

close all
g=@(x)1./(1+25*x.^2)
x=linspace(-1,1);
plot(x,g(x),'k');
hold on; axis([-1,1,-0.4,1.3]);grid on;shg

%
%

for n = [4,6,8,10,20]
  xd = linspace(-1,1,n);
  yd = g(xd);
  plot(x,polyval(polyfit(xd,yd,n-1),x));shg
  pause(1)
end
title('Polynomials of deg 4,6,8,10,20')
shg

%
%% Chebyshev points
% Next choose x-points as non-uniformly spaced Chebyshev-points
%
% x_j =  \cos(\frac{j \pi}{N}), j = 0\ldots N

chebyshevPts = @(n)( cos( (0:n)*pi/n ) )
close all
x=linspace(-1,1);
plot(x,g(x),'k')
hold on;grid on
for n = [4,6,8,10,20]
  xd = chebyshevPts(n);
  yd = g(xd);
  c=polyfit(xd,yd,n-1);
  plot(x,polyval(c,x));shg
  pause(1)
end
title('Interpolation on Chebysev points, deg 4,6,8,10,20')
shg

%
%% Chebysev points

x1 = chebyshevPts(9)
x2 = chebyshevPts(19);

figure

subplot(2,1,1)
plot(x1,0,'x')
title(['Cheb-points n=10'])
grid on
subplot(2,1,2)
plot(x2,0,'x')
title(['Cheb-points n=20'])
grid on
shg

%
%% LSQ fitting
% Polynomial fitting can be done using Vandermonde as seen already. polyfit does it also.
%
%% Example:
% Generate points x = 0,1,2,...10 and y = 1 + 2x + epsilon, where epsilon is nomally distributed random variable with mean 0 ja std 0.5. Fit an LSQ line to this data.

close all
clear
x = linspace(0,1,10)'; 
y=2*x+1 + randn(10,1)*0.5;
A = [ones(length(x),1) x]     % Vandermonde matrix
b=A\y
xx=linspace(min(x),max(x));
plot(x,y,'o',xx,b(1)+xx*b(2),'r')
hold on
plot(xx,2*xx+1,'b--')
% Different runs produce different random sequences, thus different
% results.

%
%% polyfit does the same

close all
x = linspace(0,1,10)'; 
y=2*x+1 + randn(10,1)*0.5;
c=polyfit(x,y,1)  % Now 1 means degree 1 polynomial.
xx=linspace(min(x),max(x));
p=polyval(c,xx);
plot(x,y,'o',xx,p,xx,2*xx+1,'b--');shg
%
% This was technically exactly same as polynomial interpolation,
% exept we chose the degree of polynomial this time to be 1.
% Needless to say, we get coefficients of LSQ-fit polynomial of deg k by
% c=polyfit(x,y,k)

%
%% Spline interpolation
% A cubic spline is a piecewise cubic polynomial. The pieces have 2^{nd} continuous derivatives at the joints ("knots"), so look just smooth. There remains 2 free parameters to choose, which are called boundary conditions. Here we choose Matlab default, s.k. "not-a-knot".
% Use of Matlab's spline-function is easy. It combines polyfit and polyval, so it is also 'smooth' for the user. Let's use the example data fom the previous:

hold off
plot(x,y,'o',xx,spline(x,y,xx))
title('Spline interpolation to previous data')
shg

%
%% Where polyfit doesn't work
% If you use other basic functions than monomials, you have to form the generalized Vandermonde:  \phi_j(t(i))
%
%% Example: Splines fitted to Runge function

close all
g=@(x)1./(1+25*x.^2)
x=linspace(-1,1);
plot(x,g(x),'k');
hold on; axis([-1,1,-0.4,1.3]);grid on;shg
n=9;
xd = linspace(-1,1,n);
yd = g(xd)
plot(x,spline(xd,yd,x));shg
title(['Spline to data of ',num2str(n), ' points'] )
legend('Runge function','Spline')
grid on;shg