function rbf
% RBF test program for a Gaussian RBF network trained either with least squares or
% regularised least squares. The training data consist of two clouds of
% two-dimensional data points each created randomly using a Gaussian distribution
% with random covariance

% constants
clmeans = [[1, 1]; [0, -1];];           % cluster means
n = 20;									% number of data points per cluster
max_sigma = 0.8;						% maximal cluster variance
n_test = 50;                            % test grid side length
lambda = 0.1;                           % regularisation parameter
rho = 0.3;                              % width of Gaussian basis function

% create data cloud of 2 Gaussian clusters
sigma = max_sigma*rand(2,1); % class omega 1 with random covariance in x and y direction
om1(:, 1) = repmat(clmeans(1, 1), n, 1) + sigma(1)*randn(n, 1); % x-component
om1(:, 2) = repmat(clmeans(1, 2), n, 1) + sigma(2)*randn(n, 1); % y-component
sigma = max_sigma*rand(2,1); % class omega 2 random covariance in x and y direction
om2(:, 1) = repmat(clmeans(2, 1), n, 1) + sigma(1)*randn(n, 1); % x-component
om2(:, 2) = repmat(clmeans(2, 2), n, 1) + sigma(2)*randn(n, 1); % y-component
input = [om1; om2;]; % concatenate both inpu matrices

% Train RBF
Y = zeros(2*n, 2*n);        % alloc design matrix
b = ones(2*n, 1);           % alloc vector containing margins = 1 
b(n+1:end) = -b(n+1:end);   % set all margins of class omega2 to -1
for i=1:2*n                 % fill design matrix
    for j=1:2*n             % each row contains one training exmaple mapped into feature space
        Y(i,j) = gauss(input(i,:)', input(j,:)', rho);  % determine phi(x_i - x)
    end
end
%a = inv(Y)*b;                          % uncomment for least squares solution
a = inv(Y'*Y + lambda*eye(2*n))*Y'*b;   % uncomment for regularised least squares solution

% make a grid of test points around training data
xmin = min(min(input(:,1)));
xmax = max(max(input(:,1)));
ymin = min(min(input(:,2)));
ymax = max(max(input(:,2)));
minval = min([xmin, ymin]);
maxval = max([ymax, xmax]);
test = 1.2*(minval:(maxval - minval)/(n_test - 1):maxval);

% compute discriminant function on all grid points
tinp = zeros(n_test^2, 2);
phi = zeros(2*n, 1);
cnt = 1;
for i=1:n_test
  for j=1:n_test
    tinp(cnt, 1) = test(i);             
    tinp(cnt, 2) = test(j);
    for k=1:2*n
        phi(k) = gauss(input(k,:)', tinp(cnt,:)', rho);
    end
    df(cnt) = a'*phi;
    cnt = cnt + 1;
  end
end


% show discriminant function
figure(1) % plot data points
prmesh = reshape(df, n_test, n_test);
shading interp;
pcolor(test, test, prmesh);
tmp = get(gca, 'children');
set(tmp(1), 'edgecolor', 'none');
hold on
plot(om1(:,1),om1(:,2),'k+');
plot(om2(:,1),om2(:,2),'ko');
hold off
end

% Gauss-RBF
function y = gauss(x1, x2, rho)
    y = exp(-0.5/rho^2*(x1 - x2)'*(x1 - x2));
    return
end