% Demonstrate Gaussian RBF networks on 1D toy datasets.

% Dataset for 1D regression
rng(1778);
N = 100; Nv = 150; 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

% RBF network
ll = 10.^([-7 -5 -3]);				% Regularisation parameter
ss = [0.05 0.5 1];				% Width
HH = [3 10 30 100];				% # BFs
g = cell(length(ll),length(ss),length(HH));
ET = zeros(length(ll),length(ss),length(HH)); EV = ET;
x = linspace(ax(1),ax(2),1000)'; fx = f(x);
gx = zeros(length(x),length(ll),length(ss),length(HH));
for k=1:length(HH)
  for i=1:length(ll)
    for j=1:length(ss)
      % Reset the seed of the RNG so every RBF network uses the same centres
      rng(77); g{i,j,k} = rbftrain(X,Y,HH(k),ss(j),ll(i));
      e = Y - rbf(X,g{i,j,k}); ET(i,j,k) = e(:)'*e(:)/2;
      ev = Yv - rbf(Xv,g{i,j,k}); EV(i,j,k) = ev(:)'*ev(:)/2;
      gx(:,i,j,k) = rbf(x,g{i,j,k});
    end
  end
end
ef = ef/N; ET = ET/N; EV = EV/Nv;


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

% Training & validation errors
figure(1); clf;
set(gcf,'Name','1D dataset: RBF networks'' training & validation error');
for i=1:length(ll)
  subplot(1,length(ll),i); hold on;
  plot(1:length(HH),squeeze(ET(i,:,:)),'o-','LineWidth',2);
  plot(1:length(HH),squeeze(EV(i,:,:)),'o--',[1 length(HH)],[ef ef],'k-');
  hold off;
  axis([1 length(HH) 0 4*ef]); box on;
  set(gca,'XTick',1:length(HH),'XTickLabel',HH);
  xlabel('H (= # basis functions)');
  ylabel('RBF network error: training (-) & validation (--)');
  title(['\lambda = ' num2str(ll(i),'%4.0e')]);
end
legend([repmat('\sigma=',length(ss),1) num2str(ss')]);
text(length(HH),ef,' \leftarrow true fcn');

% Dataset and regression function estimates
figure(2);
set(gcf,'Name','1D dataset: RBF networks for different hyperparameter values');
for i=1:length(ll)
  for j=1:length(ss)
    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,j,:)),'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(['RBF network: \lambda = ' num2str(ll(i),'%4.0e') ',  \sigma = '...
           num2str(ss(j))]);
    pause
  end
end

