% Demonstrate Gaussian RBF networks on MNIST.

load MNIST.mat

I = find(train_gnd==1); Iv = find(test_gnd==1);			% only 1s
%I = 1:size(train_gnd,1); Iv = 1:size(test_gnd,1);		% all types
N = 1000; Nv = 1000;
X = double(train_fea(I(randperm(length(I),N)),:))/255;
X = X + 1e-3*randn(size(X));
Xv = double(test_fea(Iv(randperm(length(Iv),Nv)),:))/255;
Xv = Xv + 1e-3*randn(size(Xv));
[N,D] = size(X); Nv = size(Xv,1);

% Ground-truth linear mapping and noise
F.type = 'linf';
% - Random linear mapping
% E = 10; F.W = randn(E,D)/10+0.02; F.w = randn(E,1)/10; s = 0.5;
% - Flip & clip
ox = 5; oy = 8; s = 0.001;
Z = reshape(1:D,sqrt(D),sqrt(D)); Z = Z(10:10+ox-1,15:15+oy-1); Z = Z(:);
E = length(Z); F.W=zeros(E,D); F.W(sub2ind([E D],1:E,Z'))=1; F.w=zeros(E,1);

% Construct output labels
Y = linf(X,F); Y = Y + s*randn(size(Y));
Yv = linf(Xv,F); Yv = Yv + s*randn(size(Yv));

% Linear regression, exactly by solving the normal equations
f = linfexact(X,Y);
e = Y - linf(X,f); efT = e(:)'*e(:)/2;		% Training error
ev = Yv - linf(Xv,f); efV = ev(:)'*ev(:)/2;	% Validation error

% RBF network, optimised approximately by running k-means on the centres
ll = 10.^([-14 -12 -10 -8 -7 -6 -5 -4 -3 -1]);	% Regularisation parameter
ss = [1 2 5 10 20 50 100 500 1000];		% Width
H = 500;					% # BFs
Eold = Inf; ET = zeros(length(ss),length(ll)); EV = ET;
warning('off','MATLAB:nearlySingularMatrix');
for i=1:length(ss)
  for j=1:length(ll)
    % Reset the seed of the RNG so every RBF network uses the same centres
    rng(77); gg = rbftrain(X,Y,H,ss(i),ll(j));
    e = Y - rbf(X,gg); ET(i,j) = e(:)'*e(:)/2;
    ev = Yv - rbf(Xv,gg); EV(i,j) = ev(:)'*ev(:)/2;
    if EV(i,j) < Eold Eold = EV(i,j); g = gg; end;	% keep best
  end
end
warning('on','MATLAB:nearlySingularMatrix');
efT = efT/N; efV = efV/Nv; ET = ET/N; EV = EV/Nv;

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

set(groot,'defaultAxesColorOrder',...
          [get(groot,'defaultAxesColorOrder');1 1 0;0 1 0;0 0 1;1 0 0]);

% Training & validation errors (curves)
figure(1); clf;
set(gcf,'Name','MNIST dataset: RBF networks'' training & validation error');
subplot(2,1,1); hold on;
plot([1 length(ss)],[efT efT],'k-');
h = plot(1:length(ss),ET,'o-','LineWidth',2);
hold off; axis([1 length(ss) 0 min(ET(:))*40]); box on;
set(gca,'XTick',1:length(ss),'XTickLabel',ss);
xlabel('\sigma'); ylabel('training error');
legend(h,[repmat('\lambda=',length(ll),1) num2str(ll','%4.0e')],...
       'Location','NorthWest');
text(length(ss),efT,' \leftarrow true fcn');
title(['RBF network with H = ' num2str(H) ' basis functions']);
subplot(2,1,2); hold on;
plot([1 length(ss)],[efV efV],'k-');
h = plot(1:length(ss),EV,'o--');
hold off; axis([1 length(ss) 0 min(EV(:))*40]); box on;
set(gca,'XTick',1:length(ss),'XTickLabel',ss);
xlabel('\sigma'); ylabel('validation error');
legend(h,[repmat('\lambda=',length(ll),1) num2str(ll','%4.0e')],...
       'Location','NorthWest');
text(length(ss),efV,' \leftarrow true fcn');
title(['RBF network with H = ' num2str(H) ' basis functions']);

% Validation errors (contour plot)
figure(2); clf;
set(gcf,'Name','MNIST dataset: RBF networks'' validation error');
[x,y] = meshgrid(linspace(log10(ss(1)),log10(ss(end)),2*length(ss)),...
                   linspace(log10(ll(1)),log10(ll(end)),2*length(ll)));
[SS,LL] = meshgrid(log10(ss),log10(ll));
EV1 = log10(EV'); EV1 = EV1 - min(EV1(:));
z = griddata(SS(:),LL(:),EV1(:),x,y,'cubic');
contour(x,y,z,100);
lll=length(ll); t1=[SS NaN(lll,1)]'; t2=[LL NaN(lll,1)]'; t3=[EV1 NaN(lll,1)]';
hold on; plot3(t1(:),t2(:),t3(:),'b:',SS(:),LL(:),EV1(:),'bo'); hold off;
xlabel('\sigma'); ylabel('\lambda'); title('log(validation error)');
set(gca,'XTick',log10(ss),'XTickLabel',ss,'YTick',log10(ll),'YTickLabel',ll);
[AZ,EL] = view;
for az=linspace(0,360,100) view(az,30+20*sin(az*pi/90)); pause(0.1); end;
view(AZ,EL);

% Plot of sample images and their outputs, using the best RBF network
idx = 1:10; lidx = length(idx); XX = Xv(idx,:); YY = Yv(idx,:);	% first 10 img
FXX = linf(XX,F); fXX = linf(XX,f); gXX = rbf(XX,g);
figure(3); clf; colormap(gray(256));
set(gcf,'Name','MNIST dataset: sample images and their outputs');
for i=1:lidx
  subplot(4,lidx,i);
  imagesc(reshape(XX(i,:),sqrt(D),sqrt(D)),[0 1]);
  set(gca,'XTick',[],'YTick',[]); axis image;
  if i==1 ylabel('x'); end
  subplot(4,lidx,i+lidx); imagesc(reshape(FXX(i,:),ox,oy),[0 1]);
  set(gca,'XTick',[],'YTick',[]); axis image;
  if i==1 ylabel('true f(x)'); end
  subplot(4,lidx,i+2*lidx); imagesc(reshape(fXX(i,:),ox,oy),[0 1]);
  set(gca,'XTick',[],'YTick',[]); axis image;
  if i==1 ylabel('optimal f(x)'); end
  subplot(4,lidx,i+3*lidx); imagesc(reshape(gXX(i,:),ox,oy),[0 1]);
  set(gca,'XTick',[],'YTick',[]); axis image;
  if i==1 ylabel('RBF f(x)'); end
end

% Plot of RBF centres, using the best RBF network
figure(4); clf; colormap(gray(256));
set(gcf,'Name','MNIST dataset: RBF network centres (subset)');
K2 = ceil(sqrt(min(H,100))); K1 = ceil(min(H,100)/K2);
for i=1:min(H,100)				% Subset of centres
  subplot(K1,K2,i); imagesc(reshape(g.C(i,:),28,28),[0 1]);
  set(gca,'XTick',[],'YTick',[]); axis image;
end


% Suggestions of things to try:
% - Try datasets containing different digits.
% - Vary the number of basis functions H.
% - Try using as centres a random subset vs the result of k-means.

