/*----------------------------------------------------------------------- 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 ); }