gmnplib.c

很好的matlab模式识别工具箱

        case 2: 
          aHa = tmp_aHa2;
          ac = tmp_ac2;
          /*
            Ha = Ha + gamma2*alpha(v)*(H(:,u)-H(:,v));
            alpha = alpha*(1-gamma2)+gamma2*alpha2;
          */
          tmp = alpha[v]*gamma2;
          for( i = 0; i < dim; i++ ) {
            Ha[i] = Ha[i] + tmp*(col_u[i] - col_v[i]);
          }
          alpha[u] = alpha[u] + tmp;
          alpha[v] = alpha[v] - tmp;
          break;

        case 3: 
          aHa = tmp_aHa3;
          ac = tmp_ac3;
          /* 
            Ha = gamma3*Ha + H(:,u)*(1-gamma3+gamma3*alpha(v)) - ...
                 gamma3*alpha(v)*H(:,v);
            alpha = alpha1*(1-gamma3)+gamma3*alpha2;
          */
          tmp = alpha[v]*gamma3;
          for( i = 0; i < dim; i++ ) {
            Ha[i] = gamma3*Ha[i] + col_u[i]*(1-gamma3+tmp) - tmp*col_v[i];
            alpha[i] = alpha[i]*gamma3;
          }          
          alpha[u] = alpha[u] + 1 - gamma3 + tmp;
          alpha[v] = alpha[v] - tmp;

          break;
      }  
    }

    UB = 0.5*aHa + ac;

    min_beta = PLUS_INF; 
    max_beta = MINUS_INF;
    for( i = 0; i < dim; i++ ) 
    {
       beta = Ha[i]+ vector_c[i];

       if( alpha[i] !=0 && max_beta < beta ) 
       {
         v = i;
         max_beta = beta;
       }

       if( beta < min_beta )
       { 
         u = i;
         min_beta = beta;
       }
    }    

    LB = min_beta - 0.5*aHa; 

    /* Stopping conditions */
    if( UB-LB <= tolabs ) exitflag = 1; 
    else if( UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
    else if(LB > th ) exitflag = 3;
    else if(t >= tmax) exitflag = 0; 

    /* print info */
    if(verb && (t % verb) == 0 ) {
      mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
        t, UB, LB, UB-LB,(UB-LB)/UB);
    }  

    /* Store selected values */
    if( t < History_size ) {
      History[INDEX(0,t,2)] = LB;
      History[INDEX(1,t,2)] = UB;
    }
    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)] = UB;
      
      History_size += HISTORY_BUF;
      mxFree( History );
      History = tmp_ptr;
    }
  }

  /* print info about last iteration*/
  if(verb && (t % verb) ) {
    mexPrintf("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
      UB, LB, UB-LB,(UB-LB)/UB);
  }  

  /*------------------------------------------------------- */
  /* Set outputs                                            */
  /*------------------------------------------------------- */
  (*ptr_t) = t;
  (*ptr_History) = History;

  /* Free memory */
  mxFree( Ha );
  
  return( exitflag ); 
}


/* --------------------------------------------------------------
 GMNP solver based on the Kowalczyk's algorithm 
 (Maximal Margin Perceptron).

 Usage: exitflag = gmnp_kowalczyk( &get_col, diag_H, vector_c, dim,  
                  tmax, tolabs, tolrel, th, &alpha, &t, &History ); 
-------------------------------------------------------------- */
int gmnp_kowalczyk(const void* (*get_col)(long,long),
            double *diag_H,
            double *vector_c,
            long dim, 
            long tmax,
            double tolabs,
            double tolrel,
            double th,
            double *alpha,
            long  *ptr_t,
            double **ptr_History,
            long verb)
{
  double LB;
  double UB;
  double aHa, ac;
  double tmp, tmp1;
  double tmp_aHa, tmp_ac, tmp_UB, tmp_gamma;
  double min_beta, beta;
  double x10, x20, x11, x12, x22;
  double delta;
  double gamma;
  double *History;
  double *Ha;
  double *tmp_ptr;
  double *col_inx;
  long i;
  long inx;
  long t;
  long History_size;
  int exitflag;

  /* ------------------------------------------------------------ */
  /* Initialization                                               */
  /* ------------------------------------------------------------ */

  Ha = mxCalloc(dim, sizeof(double));
  if( Ha == 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.");

  /* inx = argmin(0.5*diag_H + vector_c ); */
  for( tmp1 =  PLUS_INF, i = 0; i < dim; i++ ) {
    tmp = 0.5*diag_H[i] + vector_c[i];
    if( tmp1 > tmp) {
      tmp1 = tmp;
      inx = i;
    }
  }

  col_inx = (double*)get_col(inx,-1);

  for( min_beta = PLUS_INF, i = 0; i < dim; i++ ) 
  {
    alpha[i] = 0;
    Ha[i] = col_inx[i];

    beta = Ha[i] + vector_c[i];
    if( beta < min_beta ) {
      min_beta = beta;
    }
  }

  alpha[inx] = 1;
  aHa = diag_H[inx];
  ac = vector_c[inx];

  UB = 0.5*aHa + ac;
  LB = min_beta - 0.5*aHa;
  t = 0;
  History[INDEX(0,0,2)] = LB;
  History[INDEX(1,0,2)] = UB;

  if( verb ) {
    mexPrintf("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
      UB, LB, UB-LB,(UB-LB)/UB);
  }  

  /* Stopping conditions */
  if( UB-LB <= tolabs ) exitflag = 1;
  else if(UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
  else if(LB > th ) exitflag = 3;
  else exitflag = -1;

  /* ------------------------------------------------------------ */
  /* Main optimization loop                                       */
  /* ------------------------------------------------------------ */

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

    x11 = aHa;
    x10 = ac;
    /* searches for rule which yileds the biggest improvement */
    for( i = 0; i < dim; i++ ) 
    {
       delta = Ha[i] + vector_c[i] - aHa - ac;
    
       tmp_UB = PLUS_INF;
       if( delta < 0 ) 
       {
          /* Kozinec rule */
          x12 = Ha[i];
          x20 = vector_c[i];
          x22 = diag_H[i];
      
          tmp_gamma = (x11 - x12 + x10 - x20)/(x11 - 2*x12 + x22);
          tmp_gamma = MIN(1,tmp_gamma);

          tmp_aHa = (1-tmp_gamma)*(1-tmp_gamma)*x11+2*(1-tmp_gamma)*tmp_gamma*x12
                  +tmp_gamma*tmp_gamma*x22;
          tmp_ac = (1-tmp_gamma)*x10+tmp_gamma*x20;
          tmp_UB = 0.5*tmp_aHa + tmp_ac;
       }
       else if( delta > 0 && alpha[i] < 1 && alpha[i] > 0)
       {
          x12 = (x11-alpha[i]*Ha[i])/(1-alpha[i]);
          x22 = (x11-2*alpha[i]*Ha[i]+alpha[i]*alpha[i]*diag_H[i])/
                 ((1-alpha[i])*(1-alpha[i]));
          x20 = (x10-alpha[i]*vector_c[i])/(1-alpha[i]);

          tmp_gamma = (x11 - x12 + x10 - x20)/(x11 - 2*x12 + x22);
          tmp_gamma = MIN(1,tmp_gamma);
  
          tmp_aHa = (1-tmp_gamma)*(1-tmp_gamma)*x11+2*(1-tmp_gamma)*tmp_gamma*x12
                  + tmp_gamma*tmp_gamma*x22;
          tmp_ac = (1-tmp_gamma)*x10+tmp_gamma*x20;
          tmp_UB = 0.5*tmp_aHa + tmp_ac;
       }

       if( tmp_UB < UB ) 
       {
         UB = tmp_UB;
         gamma = tmp_gamma;
         aHa = tmp_aHa;
         ac = tmp_ac;
         inx = i;
       }
    }

    col_inx = (double*)get_col(inx,-1);

    /* Use the update with biggest improvement */
    delta = Ha[inx] + vector_c[inx] - x11 - x10;
    if( delta < 0 ) 
    {
       /* Kozinec rule */
      for(i = 0; i < dim; i++ ) {
        Ha[i] = Ha[i]*(1-gamma) + gamma*col_inx[i];
        alpha[i] = alpha[i]*(1-gamma);
      }
      alpha[inx] = alpha[inx] + gamma;
    }
    else
    {
      /* Inverse Kozinec rule */
      tmp = gamma*alpha[inx];
      tmp1 = 1-alpha[inx];
      for(i = 0; i < dim; i++ ) {
        Ha[i] = (Ha[i]*(tmp1+tmp) - tmp*col_inx[i])/tmp1; 

        alpha[i] = alpha[i]*(1-gamma) + gamma*alpha[i]/tmp1;
      }
      alpha[inx] = alpha[inx] - tmp/tmp1;
    }

    min_beta = PLUS_INF; 
    for( i = 0; i < dim; i++ ) 
    {
       beta = Ha[i]+ vector_c[i];
       if( beta < min_beta ) { min_beta = beta; }
    }    

    LB = min_beta - 0.5*aHa; 

    /* Stopping conditions */
    if( UB-LB <= tolabs ) exitflag = 1; 
    else if( UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
    else if( LB > th ) exitflag = 3;
    else if(t >= tmax) exitflag = 0; 

    /* print info */
    if(verb && (t % verb) == 0 ) {
      mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
        t, UB, LB, UB-LB,(UB-LB)/UB);
    }  

    /* Store selected values */
    if( t < History_size ) {
      History[INDEX(0,t,2)] = LB;
      History[INDEX(1,t,2)] = UB;
    }
    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)] = UB;
      
      History_size += HISTORY_BUF;
      mxFree( History );
      History = tmp_ptr;
    }
  }

  /* print info about last iteration*/
  if(verb && (t % verb) ) {
    mexPrintf("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
      UB, LB, UB-LB,(UB-LB)/UB);
  }  

  /*------------------------------------------------------- */
  /* Set outputs                                            */
  /*------------------------------------------------------- */
  (*ptr_t) = t;
  (*ptr_History) = History;

  /* Free memory */
  mxFree( Ha );
  
  return( exitflag ); 
}

/* --------------------------------------------------------------
 GMNP solver based on the Kozinec algorithm.

 Usage: exitflag = gmnp_kozinec( &get_col, diag_H, vector_c, dim,  
                  tmax, tolabs, tolrel, th, &alpha, &t, &History );
-------------------------------------------------------------- */
int gmnp_kozinec(const void* (*get_col)(long,long),
            double *diag_H,
            double *vector_c,
            long dim, 
            long tmax,
            double tolabs,
            double tolrel,
            double th,
            double *alpha,
            long  *ptr_t,
            double **ptr_History,
            long verb)
{
  double LB;
  double UB;
  double aHa, ac;
  double tmp, tmp1;
  double min_beta, beta;
  double lambda;
  double *History;
  double *Ha;
  double *tmp_ptr;
  double *col_u;
  long u, inx;
  long new_u;
  long i;
  long t;
  long History_size;
  int exitflag;

  /* ------------------------------------------------------------ */
  /* Initialization                                               */
  /* ------------------------------------------------------------ */

  Ha = mxCalloc(dim, sizeof(double));
  if( Ha == 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.");

  /* inx = argmin(0.5*diag_H + vector_c ); */
  for( tmp1 =  PLUS_INF, i = 0; i < dim; i++ ) {
    tmp = 0.5*diag_H[i] + vector_c[i];
    if( tmp1 > tmp) {
      tmp1 = tmp;
      inx = i;
    }
  }
  col_u = (double*)get_col(inx,-1);

  for( min_beta = PLUS_INF, i = 0; i < dim; i++ ) 
  {
    alpha[i] = 0;
    Ha[i] = col_u[i];

    beta = Ha[i] + vector_c[i];
    if( beta < min_beta ) {
      min_beta = beta;
      u = i;
    }
  }

  alpha[inx] = 1;
  aHa = diag_H[inx];
  ac = vector_c[inx];

  UB = 0.5*aHa + ac;
  LB = min_beta - 0.5*aHa;
  t = 0;
  History[INDEX(0,0,2)] = LB;
  History[INDEX(1,0,2)] = UB;

  if( verb ) {
    mexPrintf("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
      UB, LB, UB-LB,(UB-LB)/UB);
  }  

  /* Stopping conditions */
  if( UB-LB <= tolabs ) exitflag = 1;
  else if(UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
  else if(LB > th ) exitflag = 3;
  else exitflag = -1;

  /* ------------------------------------------------------------ */
  /* Main optimization loop                                       */
  /* ------------------------------------------------------------ */

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

    col_u = (double*)get_col(u,-1);

    /* Adaptation rule and update */
    lambda = (aHa - Ha[u] + ac - vector_c[u])/(aHa - 2*Ha[u] + diag_H[u]);
    lambda = MIN(1,lambda);

    aHa = aHa*(1-lambda)*(1-lambda) + 2*lambda*(1-lambda)*Ha[u] 
         + diag_H[u]*lambda*lambda;
    ac = ac*(1-lambda) + lambda*vector_c[u];

    min_beta = PLUS_INF;
    for( i = 0; i < dim; i++ ) {
      alpha[i] = alpha[i]*(1-lambda);
      Ha[i] = Ha[i]*(1-lambda) + lambda*col_u[i];

      beta = Ha[i] + vector_c[i];
      if( min_beta > beta ) { 
        min_beta = beta; 
        new_u = i; 
      }
    }
    alpha[u] = alpha[u] + lambda;

    UB = 0.5*aHa + ac;
    LB = min_beta - 0.5*aHa; 
    u = new_u;    


    /* Stopping conditions */
    if( UB-LB <= tolabs ) exitflag = 1; 
    else if( UB-LB <= ABS(UB)*tolrel ) exitflag = 2;
    else if( LB > th ) exitflag = 3;
    else if(t >= tmax) exitflag = 0; 

    /* print info about last iteration*/
    if(verb && (t % verb) == 0 ) {
      mexPrintf("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
        t, UB, LB, UB-LB,(UB-LB)/UB);
    }  

    /* Store selected values */
    if( t < History_size ) {
      History[INDEX(0,t,2)] = LB;
      History[INDEX(1,t,2)] = UB;
    }
    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)] = UB;
      
      History_size += HISTORY_BUF;
      mxFree( History );
      History = tmp_ptr;
    }
  }

  /* print info about last iteration*/
  if(verb && (t % verb) ) {
    mexPrintf("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
      UB, LB, UB-LB,(UB-LB)/UB);
  }  

  /*------------------------------------------------------- */
  /* Set outputs                                            */
  /*------------------------------------------------------- */
  (*ptr_t) = t;
  (*ptr_History) = History;

  /* Free memory */
  mxFree( Ha );
  
  return( exitflag ); 
}