kmatrix.c

this a SVM toolbox,it is very useful for someone who just learn SVM.In order to be undestood easily,

/* --------------------------------------------------------------------
 [K] = kmatrix( data, ker, arg )
 [K] = kmatrix( dataA, dataB, ker, arg )
 Computes kernel matrix of given data.

  1) K(i,j) = kernel( data(:,i), data(:,j)) i,j=1,...,N  
 or 
  2) K(i,j) = kernel( dataA(:,i), dataB(:,j)) i=1,...,N1, j=1,...,N2

 Example: for 'linear' kernel it returns 1) data'*data or 2) dataA'*dataB.

 Make executable file by a command 'mex kmatrix.c kernel.c'.
 
 Inputs:
  1) data [D x N ] matrix of N training D-dimensional patterns.
 or 
  2) dataA [D x N1 ] matrix of N1 training D-dimensional patterns.
     dataB [D x N2 ] matrix of N2 training D-dimensional patterns.

  ker [string] kernel identifier.
  arg [real] kernel argument.

    ker      arg      Kernel function
   -----------------------------------------------------------------
    'linear' []       Linear kernel: k(a,b) = a'*b
    'poly'   d [int]  Polynom: k(a,b)=(a'*b +1)^d
    'rbf'    s [real] Radial Basis Functions: k(a,b)=exp^(0.5*||a-b||^2/s^2)

 Outputs:
   1) K [N x N] kernel matrix.
  or
   2) K [N1 x N2] kernel matrix.

 Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
 (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
 Written Vojtech Franc (diploma thesis) 02.11.1999

 Modifications:
 13-sep-2002, VF
 21-October-2001, V.Franc.
 30-September-2001, V.Franc, created.
 -------------------------------------------------------------------- */

#include "mex.h"
#include "matrix.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "kernel.h"

/* case insensitive string comparision */
#ifdef __BORLANDC__
   #define STR_COMPARE(A,B,C)      strncmpi(A,B,C)  /* Borland */
#else
  #define STR_COMPARE(A,B,C)      strncasecmp(A,B,C) /* Linux gcc */
#endif

/* ==============================================================
 Main MEX function - interface to Matlab.
============================================================== */
void mexFunction( int nlhs, mxArray *plhs[],
		  int nrhs, const mxArray*prhs[] )
{
   long i, j, N1, N2;
   double tmp;
   double *K;
   char skernel[10];

   /* ---- check input arguments  ----------------------- */

   if( nrhs == 3) 
   {
      /* data matrix [dim x N1 ] */
      if( !mxIsNumeric(prhs[0]) || !mxIsDouble(prhs[0]) ||
        mxIsEmpty(prhs[0])    || mxIsComplex(prhs[0]) )
        mexErrMsgTxt("Input data must be a real matrix.");

      /* a string as kernel identifier ('linear',poly','rbf' ) */
      if( mxIsChar(prhs[1]) != 1 || mxGetM(prhs[1]) != 1 )
        mexErrMsgTxt("Input ker must be a string");
      else {
         /* check which kernel  */
         mxGetString( prhs[1], skernel, 10 );
       if( STR_COMPARE( skernel, "linear", 6) == 0 ) {
          ker = 0;
       } else if( STR_COMPARE( skernel, "poly", 4) == 0 ) {
          ker = 1;
       } else if( STR_COMPARE( skernel, "rbf", 3) == 0 ) {
          ker = 2;
       } else
          mexErrMsgTxt("Unknown kernel identifier.");
     }

     /*  real input argument for polynomial and rbf kernel   */
     if( ker == 1 || ker == 2) {
       if( !mxIsNumeric(prhs[2]) || !mxIsDouble(prhs[2]) ||
         mxIsEmpty(prhs[2])    || mxIsComplex(prhs[2]) ||
         mxGetN(prhs[2]) != 1  || mxGetM(prhs[2]) != 1 )
         mexErrMsgTxt("Input arg must be a real scalar.");
       else {
         arg1 = mxGetScalar(prhs[2]);  /* take kernel argument */

         /* if kernel is RBF than recompute its argument */
         if( ker == 2) arg1 = -2*arg1*arg1;
      }
     }

     dataA = mxGetPr(prhs[0]);    /* pointer at patterns */
     dataB = dataA;
     dim = mxGetM(prhs[0]);       /* data dimension */
     N1 = mxGetN(prhs[0]);         /* number of data */

     /* creates output kernel matrix. */
     plhs[0] = mxCreateDoubleMatrix(N1,N1,mxREAL);
     K = mxGetPr(plhs[0]);

     /* computes kenrel matrix. */
     for( i = 0; i < N1; i++ ) {
        for( j = i; j < N1; j++ ) {
           tmp = kernel( i, j );
           K[i*N1+j] = tmp;
           K[j*N1+i] = tmp;
        }
     }
   } 
   else if( nrhs == 4)
   {
      /* data matrix [dim x N1 ] */
      if( !mxIsNumeric(prhs[0]) || !mxIsDouble(prhs[0]) ||
        mxIsEmpty(prhs[0])    || mxIsComplex(prhs[0]) )
        mexErrMsgTxt("Input dataA must be a real matrix.");

      /* data matrix [dim x N1 ] */
      if( !mxIsNumeric(prhs[1]) || !mxIsDouble(prhs[1]) ||
        mxIsEmpty(prhs[1])    || mxIsComplex(prhs[1]) )
        mexErrMsgTxt("Input dataB must be a real matrix.");

      /* a string as kernel identifier ('linear',poly','rbf' ) */
      if( mxIsChar(prhs[2]) != 1 || mxGetM(prhs[2]) != 1 )
        mexErrMsgTxt("Input ker must be a string");
      else {
         /* check which kernel  */
         mxGetString( prhs[2], skernel, 10 );
       if( STR_COMPARE( skernel, "linear", 6) == 0 ) {
          ker = 0;
       } else if( STR_COMPARE( skernel, "poly", 4) == 0 ) {
          ker = 1;
       } else if( STR_COMPARE( skernel, "rbf", 3) == 0 ) {
          ker = 2;
       } else
          mexErrMsgTxt("Unknown kernel identifier.");
     }

     /*  real input argument for polynomial and rbf kernel   */
     if( ker == 1 || ker == 2) {
       if( !mxIsNumeric(prhs[3]) || !mxIsDouble(prhs[3]) ||
         mxIsEmpty(prhs[3])    || mxIsComplex(prhs[3]) ||
         mxGetN(prhs[3]) != 1  || mxGetM(prhs[3]) != 1 )
         mexErrMsgTxt("Input arg must be a real scalar.");
       else {
         arg1 = mxGetScalar(prhs[3]);  /* take kernel argument */

         /* if kernel is RBF than recompute its argument */
         if( ker == 2) arg1 = -2*arg1*arg1;
      }
     }

     dataA = mxGetPr(prhs[0]);    /* pointer at patterns */
     dataB = mxGetPr(prhs[1]);    /* pointer at patterns */
     dim = mxGetM(prhs[0]);       /* data dimension */
     N1 = mxGetN(prhs[0]);         /* number of data */
     N2 = mxGetN(prhs[1]);         /* number of data */

     /* creates output kernel matrix. */
     plhs[0] = mxCreateDoubleMatrix(N1,N2,mxREAL);
     K = mxGetPr(plhs[0]);

     /* computes kenrel matrix. */
     for( i = 0; i < N1; i++ ) {
        for( j = 0; j < N2; j++ ) {
           K[j*N1+i] = kernel( i, j );
        }
     }
   }
   else
   {
      mexErrMsgTxt("Wrong number of input arguments.");
   }
}