function H = kernel(ker,x1,x2,par)
% Function that computes the kernel matrix.
%
%  H = kernel(ker,x1,x2,par)
%
%  Parameters:
%
%              ker: kernel type. String.
%                  'linear'  - Linear kernel
%                  'poly_h'  - Homegenious polynomial kernel. par is degree of polynomial. 
%                  'poly_i'  - Inhomogenious polynomial kernel. par is degree of polynomial.
%                  'rbf'     - RBF kernel. par is width of rbfs (sigma).
%
%              x1: matrix of input vectors. Each row of x1 represents a d-dimensional vector.
%              x2: matrix of input vectors. Each row of x2 represents a d-dimensional vector.
%              par: parameter of the kernel.
%
%  Outputs:   
%
%              H: Matrix of dimension size(x1,1)*size(x2,1) and in each entry (H)i,j we have the kernel of x1(i,:) with x2(j,:).

  if (nargin < 3) % check correct number of arguments
     help kernel
  else

    % could check for correct number of args in here
    % but will slow things down further
    switch lower(ker)
    case 'linear'
       H = x1*x2';
    case 'poly_h'
       H = (x1*x2').^par;       
    case 'poly_i'
       H = (x1*x2' +1).^par;       
    case 'rbf'
       if(size(x1,2)==1)
          N1=size(x1,1);
          N2=size(x2,1);
          H = zeros(N1,N2);
          for i=1:N1
             H(i,:) = (exp(-(x2-ones(N2,1)*x1(i,:))'.*(x2-ones(N2,1)*x1(i,:))'/(2*par^2)));
          end
       else
          N1=size(x1,1);
          N2=size(x2,1);
          H = zeros(N1,N2);
          for i=1:N1
             H(i,:) = exp(-sum((x2-ones(N2,1)*x1(i,:))'.*(x2-ones(N2,1)*x1(i,:))')/(2*par^2));
          end
       end 
    otherwise
       H = x1*x2';
    end        
 end
