kgilbert.c

支持向量机工具箱

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

KGILBERT kernel Gilbert's 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
   stop [real] defines precision of the found hyperplane;
   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
 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) )


double kadd;         /* diagonal additional term */

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

/* ==============================================================
 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 *stop;         /* 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 */
 

  /* ---- 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 */

   stop = mxGetPr(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( inx1 = -1, inx2 = -1, i=0; i < N && (inx1==-1 || inx2==-1); i++ ) {
     if( labels[i] == 1 && inx1 == -1) {
        inx1 = i;
        alpha[i] = 1;
     } else if( labels[i] == 2 && inx2 == -1) {
        alpha[i] = -1; 
        inx2 = i;
     }
   }

   /* inits cache values */
   //   mexPrintf("wx=");
   for( i=0; i < N; i++) {
     wx[i] = ckernel(inx1,i) -ckernel(inx2,i);
     //     mexPrintf("%f ", wx[i]);
   }
   
   w2 = ckernel(inx1,inx1)+ckernel(inx2,inx2)-2*ckernel(inx1,inx2);
   //   mexPrintf("\nw2=%f\n", w2);
 
   sol=0;
   t = 0;

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

      t++;
   
      // -- compute auxciliary variables --   

//         [minw1x,t1]=min(wx(inx1));
//   [minw2x,t2]=max(wx(inx2));

     minwx1=PLUS_INF;
     maxwx2=MINUS_INF;
     for(i=0; i < N; i++ ) {
       if( labels[i] ==1 ) 
       {
         if( minwx1 > wx[i] ) {minwx1=wx[i]; t1 = i; }
       }
       else 
       {
         if( maxwx2 < wx[i] ) {maxwx2=wx[i]; t2 = i; }
       }           
     }

     //     mexPrintf("minwx1=%f, t1=%d, maxwx2=%f, t2=%d\n",minwx1,t1,maxwx2,t2);
      
     /* --- stoping condition for the 1st class ------ */
     if( (stop[0]==1 && (sqrt(w2)-(minwx1-maxwx2)/sqrt(w2)) >= stop[1] ) ||
         (stop[0]==2 && (1-(minwx1-maxwx2)/w2 >= stop[1]) )
       )
     {

//      k=min(1,abs((w2-wx(inx1(t1))+wx(inx2(t2)))/...
//          (w2-2*(wx(inx1(t1))-wx(inx2(t2))) + ...
//          (xt1'*xt1-2*xt1'*xt2 + xt2'*xt2))));

       ker11 = ckernel(t1,t1);
       ker12 = ckernel(t1,t2);
       ker22 = ckernel(t2,t2);
      
       k=MIN(1, ABS( (w2 - wx[t1] + wx[t2])/
                     (w2-2*(wx[t1]-wx[t2])+(ker11-2*ker12+ker22))));
       
                     
//      w2=k^2*(xt1'*xt1+xt2'*xt2-2*xt1'*xt2)+(1-k)^2*w2+...
//          2*k*(1-k)*(wx(inx1(t1))-wx(inx2(t2)));

       w2=k*k*(ker11+ker22-2*ker12)+(1-k)*(1-k)*w2+
          2*k*(1-k)*(wx[t1]-wx[t2]);
       //       mexPrintf("k=%f, w2=%f\n",k,w2);
      
//     %----------------------
//     for i=1:num_data,
//%       wx(i)=(1-k)*(wx(i)-a1*data(:,i)'*xt1-a2*data(:,i)'*xt2)...
//%           +((1-k)*a1+k)*data(:,i)'*xt1...
//%           +((1-k)*a2-k)*data(:,i)'*xt2;
//        wx(i)=k*(data(:,i)'*xt1-data(:,i)'*xt2)+(1-k)*wx(i);
//      end

//        alpha=alpha*(1-k);
//        alpha(inx1(t1))=alpha(inx1(t1))+k;
//        alpha(inx2(t2))=alpha(inx2(t2))-k;
      
       //       mexPrintf("wx=");
        for(i=0; i <N; i++ ) {
          wx[i] = k*(ckernel(i,t1)-ckernel(i,t2))+(1-k)*wx[i];

          //          mexPrintf("%f ", wx[i]);          
          alpha[i]=alpha[i]*(1-k);
        }
        //       mexPrintf("\n");
   
        alpha[t1]+=k;
        alpha[t2]-=k;

        //        mexPrintf("t=%d:alpha=",t);
        //        for(i=0;i<N;i++) mexPrintf("%f ",alpha[i]);
        //        mexPrintf("\n");
        
     }
     else
     {
        sol=1;
     }
  
     if( w2  <= 0 ) {
        // algorithm has converged to the zero vector --> classes overlap
        sol = -1;
     }

   }  // while(...)

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

   // sqared margin in transformed space
   margin2 = w2; 

   // threshold after normalization
   plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
   bias = mxGetPr(plhs[1]);
   *bias = -(minwx1 + maxwx2)/margin2;

   //  mexPrintf("0.5*(min <w,x1> - max<w,x2>)/|w|=%f\n", 0.5*(minwx1 - maxwx2)/sqrt(w2));

   // solution (normal vect. in the transformed space) after normalization
   for( i=0; i < N; i++ ) {
     if(labels[i]==1) alpha[i] *= 2/margin2; else alpha[i] *= -2/margin2;
   }

   // 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;

   }

   // number of kernel evaluations
   if( nlhs >= 5 ) {
     plhs[4] = mxCreateDoubleMatrix(1,1,mxREAL);
     (*mxGetPr(plhs[4])) = (double)ker_cnt;
   }

   // solution 1 (found), 0 (not found), -1 (does not exist)
   if( nlhs >= 3 ) {
     plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL);
     (*mxGetPr(plhs[2])) = (double)sol;
   }

   // number of iterations
   if( nlhs >= 4 ) {
     plhs[3] = mxCreateDoubleMatrix(1,1,mxREAL);
     (*mxGetPr(plhs[3])) = (double)t;
   }


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