% Demonstrate polynomial regression on 1D toy datasets.

% Dataset for 1D regression
rng(1778);
%N = 6; Nv = 10; D = 1; DD = 1;
N = 100; Nv = 200; D = 1; DD = 1;
a = 4; b = -0.5; f = @(x) sin(a*x) + b; s = 0.1;	% Nonlinear fcn & noise
X = 3*rand(N,D); Y = f(X) + s*randn(N,DD);		% Training set
Xv = 3*rand(Nv,D); Yv = f(Xv) + s*randn(Nv,DD);		% Validation set
ax = [-1 4 -2 1]; AX = [-1.2 0.2;-0.2 1];		% Ranges for plots

% Error of the true function (> optimal error because of the noise added)
e = Y - f(X); ef = e(:)'*e(:)/2;			% Optimal error

% Polynomial f(x) = w0 + w1.x + ... + wH.x^H
ll = 10.^([-10 -2 0 5]);			% Regularisation parameter
HH = [1 2 3 4 6 2*N];				% Polynomial degree + 1
g = cell(length(ll),length(HH));
ET = zeros(length(ll),length(HH)); EV = ET;
x = linspace(ax(1),ax(2),10000)'; fx = f(x);
gx = zeros(length(x),length(ll),length(HH));
warning('off','MATLAB:nearlySingularMatrix');
for i=1:length(ll)
  for k=1:length(HH)
    g{i,k} = polytrain(X,Y,HH(k),ll(i));
    e = Y - polyf(X,g{i,k}); ET(i,k) = e(:)'*e(:)/2;
    ev = Yv - polyf(Xv,g{i,k}); EV(i,k) = ev(:)'*ev(:)/2;
    gx(:,i,k) = polyf(x,g{i,k});
  end
end
warning('on','MATLAB:nearlySingularMatrix');
ef = ef/N; ET = ET/N; EV = EV/Nv;


% ---------------------------------- Plots ----------------------------------

% Training & validation errors
figure(1); clf;
set(gcf,'Name','1D dataset: polynomials'' training & validation error');
plot(1:length(HH),ET,'o-','LineWidth',2);
legend([repmat('\lambda=',length(ll),1) num2str(ll','%4.0e')]);
hold on; plot(1:length(HH),EV,'o--',[1 length(HH)],[ef ef],'k-'); hold off;
axis([1 length(HH) 0 max(ET(:,1))]); box on;
set(gca,'XTick',1:length(HH),'XTickLabel',HH);
xlabel('H (= polynomial degree + 1)');
ylabel('Polynomial error: training (-) & validation (--)');
text(length(HH),ef,' \leftarrow true fcn');

% Dataset and regression function estimates
figure(2);
set(gcf,'Name','1D dataset: polynomials for different hyperparameter values');
for i=1:length(ll)
  set(0,'CurrentFigure',2);
  clf; plot(X,Y,'bo','MarkerSize',6);
  hold on; plot(Xv,Yv,'g+',x,fx,'k-');
  h1 = plot(x,squeeze(gx(:,i,:)),'LineWidth',2); hold off;
  axis(ax); daspect([1 1 1]); box on; xlabel('x'); ylabel('y');
  legend(h1,[repmat('H=',length(HH),1) num2str(HH')]);
  title(['Polynomial: \lambda = ' num2str(ll(i),'%4.0e')]);
  pause
end


% Suggestions of things to try:
% - Change the random seed to rng(1999);
% - Change the dataset (number of points, true function).
%   If you use few points (say N=6) you can see some polynomials interpolating
%   all the points (going through them exactly).
% - Change the hyperparameter: degree of polynomial, weight decay l.

