% Copyright (c) 2016 by Miguel A. Carreira-Perpinan
% for use in CSE176 Introduction to Machine Learning at UC Merced

% Demonstrate ensemble learning with bagging for regression.

% Dataset for 1D regression containing N training points (inputs X, outputs Y)
rng(18);
%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
et = Y - f(X); eft = et(:)'*et(:)/2/N;			% on the training set
ev = Yv - f(Xv); efv = ev(:)'*ev(:)/2/Nv;		% on the validation set

% Generate L bootstrap samples from the training set X:
% XX(:,i), YY(:,i) are the N training points (inputs, outputs) for bootstrap
% sample i (for =1:L)
L = 10; idx = randi(N,[L N]); XX = X(idx)'; YY = Y(idx)';

% Train L learners in the bootstrap samples
l = 1e-8;					% Regularisation parameter
HH = [1 2 3 4 6 10 20 30];			% Polynomial degree + 1
g = cell(length(HH),L);
ET = zeros(L,length(HH)); EV = ET;

% Outputs of the individual learners, for computing the ensemble output:
gt = zeros(N,L,length(HH));		% on the training set
gv = zeros(Nv,L,length(HH));		% on the validation set
x = linspace(ax(1),ax(2),10000)'; fx = f(x);
gx = zeros(length(x),L,length(HH));	% on an X-axis grid, for plotting

warning('off','MATLAB:nearlySingularMatrix');
for k=1:length(HH)
  for i=1:L
    % Train ith learner (on the ith bootstrap sample)
    g{i,k} = polytrain(XX(:,i),YY(:,i),HH(k),l);
    % Training set error
    gt(:,i,k) = polyf(X,g{i,k}); et=Y-gt(:,i,k); ET(i,k)=et(:)'*et(:)/2/N;
    % Validation set error
    gv(:,i,k) = polyf(Xv,g{i,k}); ev=Yv-gv(:,i,k); EV(i,k)=ev(:)'*ev(:)/2/Nv;
    % (x,y) values of the learned regressor for the plot below
    gx(:,i,k) = polyf(x,g{i,k});
  end
end
warning('on','MATLAB:nearlySingularMatrix');

% Ensemble result.
% The ensemble prediction is the mean of the L learners' predictions.
gt_ens = squeeze(mean(gt,2));		% on the training set
ET_ens = sum(bsxfun(@minus,Y,gt_ens).^2,1)/2/N;
gv_ens = squeeze(mean(gv,2));		% on the validation set
EV_ens = sum(bsxfun(@minus,Yv,gv_ens).^2,1)/2/Nv;
gx_ens = squeeze(mean(gx,2));		% on the X-axis grid, for plotting

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

% Training & validation errors, as a function of the hyperparameter
figure(1); clf;
set(gcf,'Name','Ensemble of polynomials: training & validation error');
subplot(1,2,1); cla; hold on;
h1 = plot(1:length(HH),ET_ens,'ro-','LineWidth',2);
h2 = plot(1:length(HH),ET,'bo');
h3 = plot([1 length(HH)],[eft eft],'k-');
hold off; axis([1 length(HH) 0 max([ET(:);EV(:)])]); box on;
set(gca,'XTick',1:length(HH),'XTickLabel',HH);
legend([h1;h2(1);h3],'ensemble','learners','true fcn');
xlabel('H (= polynomial degree + 1)'); ylabel('training error');
text(length(HH),eft,' \leftarrow true fcn');
subplot(1,2,2); cla; hold on;
h1 = plot(1:length(HH),EV_ens,'ro-','LineWidth',2);
h2 = plot(1:length(HH),EV,'bo');
h3 = plot([1 length(HH)],[efv efv],'k-');
hold off; axis([1 length(HH) 0 max([ET(:);EV(:)])]); box on;
set(gca,'XTick',1:length(HH),'XTickLabel',HH);
legend([h1;h2(1);h3],'ensemble','learners','true fcn');
xlabel('H (= polynomial degree + 1)'); ylabel('validation error');
text(length(HH),efv,' \leftarrow true fcn');

% Dataset and regression function estimates
figure(2);
set(gcf,'Name','Ensemble of polynomials: different hyperparameter values');
for k=1:length(HH)
  set(0,'CurrentFigure',2);
  clf; plot(X,Y,'bo','MarkerSize',12);
  hold on; plot(Xv,Yv,'g+'); h3 = plot(x,fx,'k-');
  h2 = plot(x,gx(:,:,k),'b-');
  h1 = plot(x,gx_ens(:,k),'r-','LineWidth',2); hold off;
  axis(ax); daspect([1 1 1]); box on; xlabel('x'); ylabel('y');
  legend([h1;h2(1);h3],'ensemble','learners','true fcn');
  title(['Learners: polynomials of degree ' num2str(HH(k)-1)...
         ', \lambda = ' num2str(l,'%4.0e')]);
  pause
end


% Suggestions of things to try:
% - Level of noise in the dataset.
% - Other datasets.
% - Other types of base learners beyond polynomials.

