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     - n inputs of dimension m (m-by-n matrix).
%      Y     - k inputs of dimension m (m-by-k matrix).
%      SIGMA - Sigma used within the radial basis function, default = 1.
%    and returns,
%      Z     - RBF kernel matrix (n-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, 2);

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

% Quadrierte Distanzmatrix berechnen
s = 0;
switch s
    case 2
        % Quadrierte Distanzmatrix mit zwei for-Schleifen berechnen. 
        D = zeros(n, k); % Speicher reservieren
        for i = 1:n
            for j = 1:k
                D(i, j) = norm(X(:,i) - Y(:,j))^2;
            end
        end
    case 1
        % Quadrierte Distanzmatrix mit einer for-Schleife berechnen. 
        D = zeros(n, k); % Speicher reservieren
        for i = 1:n
            D(i, :) = sum((Y - repmat(X(:,i), 1, k)).^2);
        end
    case 0
        % Quadrierte Distanzmatrix ohne for-Schleifen berechnen. 
        XX = sum(X.*X);
        YY = sum(Y.*Y);
        XY = X'*Y; 
        D = repmat(XX', 1, k) + repmat(YY, n, 1) - 2*XY;
        % X^2 - 2XY - Y^2
    otherwise
        % Interne Dist-Funktion
        D = dist(X', Y).^2;
end

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