npa.c~

一个工具包

/*-----------------------------------------------------------------------
 NPA algorithm for separable single-class SVM problem.

 int single_npa(TKerFun ker,
             long num_data,
             long tmax,
             double tolabs,
             double tolrel,
             double *Alpha,
             double *UB,
             double *LB,
             long *t,
             double *History
             int verb)

 tmax, tolabs, tolrel ... Define stopping conditions: 

    UB-LB <= tolabs           ->  exit_flag = 1   Abs. tolerance.
    (UB-LB)/(LB+1) <= tolrel  ->  exit_flag = 2   Relative tolerance.
    t >= tmax                 ->  exit_flag = 0   Number of iterations.
 
 Alpha ... Lagrangians defining found decision rule. 
 UB ... Achieved upper bound on the optimal solution.
 LB ... Achieved lower bound on the optimal solution.
 t  ... Number of iterations.
 History ... Value of LB and UB wrt. number of iterations.
 verb .... if 1 then print info about training.

 Modifications:
 16-Npv-2004, VF
 29-Feb-2004, VF
 23-Jan-2004, VF

-------------------------------------------------------------------- */


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

#define HISTORY_BUF 1000000

#define MINUS_INF INT_MIN
#define PLUS_INF  INT_MAX

#define ABS(A) ((A >= 0) ? A : -A)
#define MIN(A,B) ((A < B) ? A : B)
#define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)

/* ==============================================================
 Kernel NPA algorithm.
============================================================== */

int single_npa(const double (*kernel_fce)(long, long),
            long num_data, 
            long tmax,
            double tolabs,
            double tolrel,
            double **out_Alpha,
            double *out_UB,
            double *out_LB,
            long  *out_t,
               double **out_History,
            int verb )
{
  double *Alpha;
  double *History;
  double LB;
  double UB2;
  double *ProjX;
  double *K_Diag;
  double *tmp_ptr;
  long min_inx;
  long max_inx;
  long i;
  long t;
  long u, v;
  double kuu, kvv, kuv;
  double x11, x22, x33, x12, x23, x13;
  double f1, f2;
  double e11, e22, e12;
  double den;
  double lambda, omega;
  double lambda12, lambda13, lambda23;
  double tmp1;
  double UB123, UB12, UB13, UB23;
  double UB_buf[3];
  long nearest_segment;
  long rule1, rule2, rule3, rule4;
  long History_size;
  int exitflag;

  /* allocate memory */
  Alpha = mxCalloc(num_data, sizeof(double));
  if( Alpha == NULL ) mexErrMsgTxt("Not enough memory.");

  ProjX = mxCalloc(num_data, sizeof(double));
  if( ProjX == NULL ) mexErrMsgTxt("Not enough memory.");

  K_Diag = mxCalloc(num_data, sizeof(double));
  if( K_Diag == NULL ) mexErrMsgTxt("Not enough memory.");

  History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
  History = mxCalloc(History_size*2,sizeof(double));
  if( History == NULL ) mexErrMsgTxt("Not enough memory.");

  /* == Inicialization == */
  for( LB =  PLUS_INF, i = 0; i < num_data; i++ ) 
  {
    Alpha[i] = 0;

    ProjX[i] = kernel_fce( 0, i );

    K_Diag[i] = kernel_fce( i, i );

    if( ProjX[i] < LB ) {
      LB = ProjX[i];
      min_inx = i;
    }
  }

  max_inx = 0;
  UB2 = K_Diag[0];
  LB = LB/sqrt(UB2);
  Alpha[0] = 1;
  t = 0;
  rule1=0; rule2=0; rule3=0; rule4=0;
  History[INDEX(0,0,2)] = LB;
  History[INDEX(1,0,2)] = sqrt(UB2);

  if( sqrt(UB2) <= tolabs ) exitflag = 1;
  else if((sqrt(UB2)-LB)/(ABS(LB)+1) <= tolrel ) exitflag = 2;
  else exitflag = -1;

  /* == Main cycle == */

  while( exitflag == -1 ) 
  {
    t++;     

    if(verb)  mexPrintf("%d: ", t);

    u = max_inx;
    v = min_inx;
    kuu = K_Diag[u];
    kvv = K_Diag[v];
    kuv = kernel_fce(u,v);


    x11 = UB2;
    x12 = ProjX[v];
    x22 = K_Diag[v];
    x13 = UB2 + Alpha[u]*(ProjX[v]-ProjX[u]);
    x23 = ProjX[v] + Alpha[u]*(kvv-kuv);
    x33 = UB2 + 2*Alpha[u]*(ProjX[v]-ProjX[u]) + 
          Alpha[u]*Alpha[u]*(kvv - 2*kuv+kuu);

    /* Pocitej vzdalenost od trojuhelniku x1,x2,x3 */
    f1 = x11 - x12;
    f2 = x11 - x13;
    e11 = x22 - 2*x12 + x11;
    e22 = x33 - 2*x13 + x11;
    e12 = x23 - x13 - x12 + x11;
    den = e11*e22 - e12*e12;

    /* Trojuhelnik */
    if( den > 0 ) 
    {
       lambda = (f1*e22 - f2*e12)/den;
       omega = (e11*f2 - f1*e12)/den;
  
       if( (lambda+omega <= 1) && (lambda >= 0) && (omega >=0 ) )
       {

          UB123 = x11*(1-lambda-omega)*(1-lambda-omega) + 
                  x22*lambda*lambda + x33*omega*omega + 
                  2*x12*lambda*(1-lambda-omega) + 
                  2*x13*omega*(1-lambda-omega) +
                  2*x23*lambda*omega;

          tmp1 = Alpha[u]*omega;
          for( i = 0; i < num_data; i++ )
          {
             ProjX[i] = ProjX[i]*(1-lambda) + 
                        (lambda+tmp1)*kernel_fce(i,v) -
                        tmp1*kernel_fce(i,u);

             Alpha[i] = Alpha[i]*(1-lambda);
          }
       
          Alpha[v] = Alpha[v] + lambda + tmp1;
          Alpha[u] = Alpha[u] - tmp1;

          UB2 = UB123;
          rule4 = rule4 + 1;
/*          mexPrintf("rule 4, lambda = %f, omega = %f, ", lambda, omega);*/
       } 
       else
       {
          UB123 = PLUS_INF;
       }
    } 
    else
    {
       UB123 = PLUS_INF;
    }

    /* Zkus line segmenty 1-2, 1-3 a 2-3 */
    if( UB123 == PLUS_INF)
    {
       /* Pocitej vzdalenost od segmentu x1 - x2 */
       lambda12 = MIN( (x11 - x12)/(x11 - 2*x12 + x22) , 1 );
       UB12 = x11*(1-lambda12)*(1-lambda12)+2*lambda12*(1-lambda12)*x12
              +lambda12*lambda12*x22;

       /* Pocitej vzdalenost od segmentu x1 - x3 */
       lambda13 = MIN( (x11 - x13)/(x11 - 2*x13 + x33) , 1);
       UB13 = x11*(1-lambda13)*(1-lambda13)+2*lambda13*(1-lambda13)*x13
              +lambda13*lambda13*x33;
       
       /* Pocitej vzdalenost od segmentu x2 - x3 */
       den = (x22 - 2*x23 + x33);
       if( den > 0 )
       {
          lambda23 = MIN((x22 - x23)/den , 1 );
          UB23 = x22*(1-lambda23)*(1-lambda23)+2*lambda23*(1-lambda23)*x23
                 +lambda23*lambda23*x33;
       }
       else
       {
          UB23 = PLUS_INF;
       }
    
       /* Vyber nejblizsi segment */
/*       mexPrintf("UB12=%f,UB13=%f,UB23=%f, ", UB12, UB13, UB23);*/
       UB_buf[0] = UB12; UB_buf[1] = UB13; UB_buf[2] = UB23;
       for( UB2 = PLUS_INF, i = 0; i < 3; i++ ) 
       {
          if( UB2 > UB_buf[i] )
          {
             UB2 = UB_buf[i];
             nearest_segment = i;
          }
       }

       /* Pouzij nejblizsi segment */
       switch( nearest_segment )
       {
          case 0: 
            for( i = 0; i < num_data; i++ )
            {
               ProjX[i] = ProjX[i]*(1-lambda12) + kernel_fce(i,v)*lambda12;

               Alpha[i] = Alpha[i]*(1 - lambda12);
            }
            Alpha[v] = Alpha[v] + lambda12;

            rule1 = rule1+1;
/*            mexPrintf("rule 1, lambda12=%f, ", lambda12);*/
            break;
 
          case 1:
            tmp1 = lambda13*Alpha[u];
            for( i = 0; i < num_data; i++ )
            {
               ProjX[i] = ProjX[i] + tmp1*(kernel_fce(i,v) - kernel_fce(i,u));
            }

            Alpha[v] = Alpha[v] + tmp1;
            Alpha[u] = Alpha[u] - tmp1;

            rule2 = rule2+1;
/*            mexPrintf("rule 2, lambda13=%f, ", lambda13);*/
            break;

          case 2:
            tmp1 = Alpha[u]*lambda23;
            for( i = 0; i < num_data; i++ )
            {
               ProjX[i] = lambda23*ProjX[i] 
                          + (1-lambda23+tmp1)*kernel_fce(i,v) 
                          - tmp1*kernel_fce(i,u);

               Alpha[i] = Alpha[i]*lambda23;
            }

            Alpha[v] = Alpha[v] + 1-lambda23 + tmp1;
            Alpha[u] = Alpha[u] - tmp1;
            rule3 = rule3+1;
/*            mexPrintf("rule 3, lambda23=%f, ", lambda23);*/
            break;
 
       }
    }

    /* Pocita dolni mez, min_inx a max_inx */
    for( LB = PLUS_INF, tmp1 = MINUS_INF, i = 0; i < num_data; i++ ) 
    {
       if( LB > ProjX[i] ) 
       {
          LB = ProjX[i];
          min_inx = i;
       }

       if( (Alpha[i] > 0) && (tmp1 < ProjX[i] ))
       {
          tmp1 = ProjX[i];
          max_inx = i;
       }
    }

    LB = LB/sqrt(UB2);

    if( verb ) mexPrintf("LB=%f, UB=%f, tolabs=%f, tolrel=%f\n", 
       LB, sqrt(UB2), sqrt(UB2)-LB, (sqrt(UB2)-LB)/(ABS(LB)+1));

/*    for(i = 0; i < num_data; i++ ) {*/
/*      if( Alpha[i]>0 ) mexPrintf("%d:%f ", i+1,ProjX[i]);*/
/*    }*/
/*    mexPrintf("\n");*/

    /* Stopping conditions */
    if( sqrt(UB2)-LB <= tolabs ) exitflag = 1; 
    else if( ((sqrt(UB2)-LB)/(ABS(LB)+1)) <= tolrel ) exitflag = 2; 
    else if(t >= tmax) exitflag = 0; 

    /* Store selected values */
    if( t < History_size ) {
      History[INDEX(0,t,2)] = LB;
      History[INDEX(1,t,2)] = sqrt(UB2);
    }
    else {
      tmp_ptr = mxCalloc((History_size+HISTORY_BUF)*2,sizeof(double));
      if( tmp_ptr == NULL ) mexErrMsgTxt("Not enough memory.");
      for( i = 0; i < History_size; i++ ) {
        tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
        tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
      }
      tmp_ptr[INDEX(0,t,2)] = LB;
      tmp_ptr[INDEX(1,t,2)] = sqrt(UB2);
      
      History_size += HISTORY_BUF;
      mxFree( History );
      History = tmp_ptr;
    }

/*    if( t >= 100) exitflag=1; */

  } /* while( exitflag == -1 )  */

/*  mexPrintf("rule1=%d, rule2=%d, rule3=%d, rule4=%d\n", */
/*      rule1, rule2, rule3, rule4);*/

  /* transform Alphas to obtain canonical hyperplane representation */
  for( i = 0; i < num_data; i++ ) {
    Alpha[i] = Alpha[i] / (LB*sqrt(UB2));
  }

  /* outputs */
  (*out_Alpha) = Alpha;
  (*out_UB) = sqrt(UB2);
  (*out_LB) = LB;
  (*out_t) = t;
  (*out_History) = History;

  /**/
  mxFree( ProjX );
  mxFree( K_Diag );

  return( exitflag ); 
}