romma.c

支持向量机工具箱

/*---------------------------------------------------------------------------
[Alpha,bias,sol,t,kercnt,margin,trnerr]= 
  romma(data,labels,eps,ker,arg,tmax,C) 

ROMMA Relaxed On-line Maximum Margin Algorithm.
 It solves the Support vector Machines problem with quadratic 
 cost function for classification violations.

 Inputs:
   data [dim x N] training patterns
   labels [1 x N] labels of training patterns
   eps [real] defines maximal possible violation of KKT-cond.;
   ker [string] kernel, see 'help kernel'.
   arg [...] argument of given kernel, see 'help kernel'.
   tmax [int] maximal number of iterations.
   C [real] trade-off between margin and training error.
  
 Outputs:
   Alpha [1xN] Lagrangians defining found decision rule.
   bias [real] bias (threshold) of found decision rule.
   sol [int] 1 solution is found
             0 algorithm stoped (t == tmax) before converged.
            -1 hyperplane with margin greater then epsilon 
               does not exist.
   t [int] number of iterations.
   kercnt [int] number of kernel evaluations.
   margin [real] margin between classes.
   trnerr [real] training error.

 See also SVM.

 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, 13.4.2000
 Modifications
 20-jun-2002, VF
 14-jun-2002, VF
-------------------------------------------------------------------- */

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

#include "kernel.h"

#define MINUS_INF INT_MIN
#define PLUS_INF  INT_MAX


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

#define MAX(A,B)   (((A) > (B)) ? (A) : (B) )
#define MIN(A,B)   (((A) < (B)) ? (A) : (B) )
#define ABS(A)   (((A) < (0)) ? (-A) : (A) )

#define y(A) (labels[A]*2-3)

double kadd;         /* diagonal additional term */

/*------------------------------------------------*/
double ckernel( long i, long j) 
{
  if( i!=j ) return( kernel(i,j)+1); else return(kernel(i,j)+kadd+1);
}

/* ==============================================================
 Main MEX function - interface to Matlab.
============================================================== */
void mexFunction( int nlhs, mxArray *plhs[],
		  int nrhs, const mxArray*prhs[] )
{
   char skernel[10];
   long t;              /* iteration number */
   long i, j;           /* loop variables */
   int sol;             /* solution: 1=found, 0=not found, -1=does not exist*/

   long inx1, inx2;     /* --//--              */
   double k;            /* --//--              */
   double *wx;
   double minwx1, maxwx2;
   long t1, t2;
   double ker11, ker12, ker22;
   double w2;
   double margin2;
 
   double *labels;      /* pointer at labels */
   long N;              /* number of training patterns */
   double eps;         /* stopping criterion */
   long tmax;           /* maximal number of iterations */
   double C;            /* trade-off constant */

   double *alpha;       /* Lagrangians */
   double *bias;        /* threshold of the learned indicator function */
   double margin;       /* margin in the original space */
   double trn_err;      /* training error */
   double dfun;         /* value of decision function */
   double ct, dt, tmp, kk;
 

  /* ---- CHECK INPUT ARGUMENTS  ----------------------- */

   if(nrhs < 7)
      mexErrMsgTxt("Not enough input arguments.");
   if(nlhs < 3)
      mexErrMsgTxt("Not enough output arguments.");


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

   /* labels [1 x N ] */
   if( !mxIsNumeric(prhs[1]) || !mxIsDouble(prhs[1]) ||
       mxIsEmpty(prhs[1])    || mxIsComplex(prhs[1]) )
      mexErrMsgTxt("Input I must be a real vector.");

   /*  stopping condition */
   if( !mxIsNumeric(prhs[2]) || !mxIsDouble(prhs[2]) ||
       mxIsEmpty(prhs[2])    || mxIsComplex(prhs[2]))
      mexErrMsgTxt("Input stop must be a real number.");

   /* a string as kernel identifier ('linear',poly','rbf' ) */
   if( mxIsChar(prhs[3]) != 1 || mxGetM(prhs[3]) != 1 )
      mexErrMsgTxt("Input ker must be a string");
   else {
       /* check which kernel  */
       mxGetString( prhs[3], 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[4]) || !mxIsDouble(prhs[4]) ||
         mxIsEmpty(prhs[4])    || mxIsComplex(prhs[4]) ||
         mxGetN(prhs[4]) != 1  || mxGetM(prhs[4]) != 1 )
         mexErrMsgTxt("Input arg must be a real scalar.");
      else {
         arg1 = mxGetScalar(prhs[4]);  /* take kernel argument */

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

   /*  tmax  */
   if( !mxIsNumeric(prhs[5]) || !mxIsDouble(prhs[5]) ||
       mxIsEmpty(prhs[5])    || mxIsComplex(prhs[5]) ||
       (mxGetN(prhs[5]) != 1  && mxGetM(prhs[5]) != 1 ))
      mexErrMsgTxt("Input tmax must be an integer.");

   /*  one or two real trade-off*/
   if( !mxIsNumeric(prhs[6]) || !mxIsDouble(prhs[6]) ||
       mxIsEmpty(prhs[6])    || mxIsComplex(prhs[6]) ||
       (mxGetN(prhs[6]) != 1  && mxGetM(prhs[6]) != 1 ))
      mexErrMsgTxt("Input C must be a real scalar.");


   /* ---- GET INPUT ARGUMENTS ------------------------------- */

   dataA = mxGetPr(prhs[0]);  /* pointer at patterns */
   dataB = mxGetPr(prhs[0]);  /* pointer at patterns */
   labels = mxGetPr(prhs[1]); /* pointer at labels */

   dim = mxGetM(prhs[0]);     /* data dimension */
   N = mxGetN(prhs[0]);       /* number of data */

   eps = mxGetScalar(prhs[2]);

   if( mxIsInf( mxGetScalar(prhs[5])) ) {
      tmax = INT_MAX;
   } else {
     tmax = (long)mxGetScalar(prhs[5]);
   }
   C = mxGetScalar(prhs[6]);

   // computes additional term to kernel value on the diagonal
   if( C != 0 ) kadd = 1/(2*C); else kadd = 0;
 
   /* create vector for Lagrangeians */
   plhs[0] = mxCreateDoubleMatrix(1,N,mxREAL);
   alpha = mxGetPr(plhs[0]);


   /*-- INICIALIZATION ------------------------------*/

   ker_cnt = 0;  /* counter for number of kernel evaluetions */

   // inicialization of cached values
   if( (wx = mxCalloc(N, sizeof(double))) == NULL) {
      mexErrMsgTxt("Not enough memory for error cache.");
   }

   /* takes two vectors as an initial solution */
   for( i=0; i < N; i++ ) {
        alpha[i] = 0;
   }

   alpha[0]=y(0)/ckernel(0,0);
   w2 = alpha[0]*alpha[0]*ckernel(0,0);

   //aKa=alpha(1)*alpha(1)*K(1,1);
   //xka=zeros(num_data,1);
   //for i=1:num_data,
   //  xka(i)=alpha(1)*K(i,1);
   //end

   /* inits cache values */
   //   mexPrintf("wx=");
   for( i=0; i < N; i++) {
     wx[i] = alpha[0]*ckernel(0,i);
     //     mexPrintf("%f ", wx[i]);
   }
   
   sol=0;
   t = 0;

   /* -- MAIN OPTIMIZATION CYCLE ------------------------ */
   while( sol == 0 && tmax > t )
   {

      t++;

      sol=1;
      for( i=0; i< N; i++ )
      {

        //      ywx = yt*K(i,:)*alpha;
    
        if( y(i)*wx[i] < 1-eps) 
        {

          kk=ckernel(i,i);
  
          //           ct= (K(i,i)*aKa-yt*xka(i))/(K(i,i)*aKa - xka(i)^2);
          ct=(kk*w2 -y(i)*wx[i])/(kk*w2 - wx[i]*wx[i]);

          //       dt= (aKa*(yt-xka(i)))/(K(i,i)*aKa - xka(i)^2);
          dt=(w2*(y(i)-wx[i]))/(kk*w2 - wx[i]*wx[i]);
      
          //      tmp=0;
          //      for j=1:num_data,
          //        tmp=tmp+alpha(j)*K(i,j);
          //        xka(j) = ct*xka(j) + dt*K(i,j);
          //      end
          //      aKa= ct^2 * aKa + 2*ct*dt*tmp + dt^2*K(i,i);

          for( tmp=0,j=0; j <N; j++ ) {
            tmp += alpha[j]*ckernel(i,j);
            wx[j] = ct*wx[j] + dt*ckernel(i,j);
          }
          w2=ct*ct*w2 + 2*ct*dt*tmp + dt*dt*kk;

          for( tmp=0,j=0; j <N; j++ ) {
            alpha[j] = ct*alpha[j];
          }
          alpha[i] += dt;
          //      alpha=ct*alpha;
          //      alpha(i)=alpha(i)+dt;
              
          sol=0;
        }
      }   
   }  // while(...)

   /* --- COMPUTATION OF OUTPUT VALUES ----------------------- */

   // threshold 
   plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
   bias = mxGetPr(plhs[1]);
   for(tmp=0, i=0; i < N; i++) {
     tmp += alpha[i];
     alpha[i]*=y(i);
   }
   *bias=-tmp;

   //b=sum(alpha);
   //alpha=(alpha.*y)';

   // 
   plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL);
   *(mxGetPr(plhs[2]))=sol;

   plhs[3] = mxCreateDoubleMatrix(1,1,mxREAL);
   *(mxGetPr(plhs[3]))=t;

   plhs[4] = mxCreateDoubleMatrix(1,1,mxREAL);
   *(mxGetPr(plhs[4]))=ker_cnt;


   // compute margin 
   if( nlhs >= 6 ) {
     margin = 0;

     margin = 0;
     for(i = 0; i < N; i++ ) { 
       for( j=0; j < N; j++ ) { 
         if( alpha[i] != 0 && alpha[j] != 0 ) {
            if( labels[i] == labels[j] )
               margin += alpha[i]*alpha[j]*kernel(i,j);  
            else
               margin -= alpha[i]*alpha[j]*kernel(i,j);  
          }
       }
     }

     margin = 1/sqrt(margin); 
     plhs[5] = mxCreateDoubleMatrix(1,1,mxREAL);
     (*mxGetPr(plhs[5])) = margin;
   }


   // training errors
   if( nlhs >= 7 ) 
   {

     trn_err = 0;
     for(i = 0; i < N; i++ ) {
       dfun = 0;

       for( j=0; j < N; j++ ) {
         if( alpha[j] != 0 ) {
            if( labels[j] == 1) 
              dfun += alpha[j]*kernel(i,j); 
            else
              dfun -= alpha[j]*kernel(i,j);
         }
       }

       if( (3-labels[i]*2)*(dfun + *bias) < 0) trn_err++;
     }

     plhs[6] = mxCreateDoubleMatrix(1,1,mxREAL);
     (*mxGetPr(plhs[6])) = trn_err/N;

   }


   /* ----- FREE MEMORY ----------------------- */
   mxFree( wx );
 
}