function Z = rbf(X, Y, sigma)
%RBF Radial basis function.
%
%  Syntax
%
%    Z = rbf(X, Y, sigma)
%
%
%  Description
%
%    RBF computes the radial basis function 
%         Z = exp(-dist(X, Y).^2/(2*sigma^2))
%    where dist(X, Y) = sqrt(||X||.^2 + ||Y||.^2 - 2*X.*Y)
%
%    Z = RBF(X,Y,SIGMA) takes,
%      X     - m inputs of dimension n (m-by-n matrix).
%      Y     - k inputs of dimension n (k-by-n matrix).
%      SIGMA - Sigma used within the radial basis function, default = 1.
%    and returns,
%      Z     - RBF kernel matrix (m-by-k matrix).

if nargin < 3
    sigma = 1;
end

if sigma < eps
    error('Sigma must be positive.');
end

[m,n] = size(X);
k = size(Y, 1);

if n ~= size(Y, 2) || m < 1 || k < 1
    error('First and/or second argument have wrong dimensions.');
end


s = 0;
switch s
    case 2
        % two for loops 
        D = zeros(m, k);
        for i = 1:m
            for j = 1:k
                D(i, j) = norm(X(i,:) - Y(j,:))^2;
            end
        end
    case 1
        % one for loop 
        D = zeros(m, k);
        for i = 1:m
            D(i, :) = sum((Y - repmat(X(i,:), k, 1))'.^2);
        end
    case 0
        % no for loops
        XX = sum(X'.*X');
        YY = sum(Y'.*Y');
        XY = X*Y'; 
        D = repmat(XX', 1, k) + repmat(YY, m, 1) - 2*XY;
        % X^2 - 2XY + Y^2   
end

Z = exp(-D/(2*sigma^2));
