svm.cpp

这是核学习的一个基础软件包

		else
		{
			alpha[i] = min(1.0,sum_neg);
			sum_neg -= alpha[i];
		}

	double *zeros = new double[l];

	for(i=0;i<l;i++)
		zeros[i] = 0;

	Solver_NU s;
	s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
		alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
	double r = si->r;

	//info("C = %f\n",1/r);

	for(i=0;i<l;i++)
		alpha[i] *= y[i]/r;

	si->rho /= r;
	si->obj /= (r*r);
	si->upper_bound_p = 1/r;
	si->upper_bound_n = 1/r;

	delete[] y;
	delete[] zeros;
}

static void solve_one_class(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver::SolutionInfo* si)
{
	int l = prob->l;
	double *zeros = new double[l];
	schar *ones = new schar[l];
	int i;

	int n = (int)(param->nu*prob->l);	// # of alpha's at upper bound

	for(i=0;i<n;i++)
		alpha[i] = 1;
	alpha[n] = param->nu * prob->l - n;
	for(i=n+1;i<l;i++)
		alpha[i] = 0;

	for(i=0;i<l;i++)
	{
		zeros[i] = 0;
		ones[i] = 1;
	}

	Solver s;
	s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
		alpha, 1.0, 1.0, param->eps, si, param->shrinking);

	delete[] zeros;
	delete[] ones;
}

static void solve_epsilon_svr(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver::SolutionInfo* si)
{
	int l = prob->l;
	double *alpha2 = new double[2*l];
	double *linear_term = new double[2*l];
	schar *y = new schar[2*l];
	int i;

	for(i=0;i<l;i++)
	{
		alpha2[i] = 0;
		linear_term[i] = param->p - prob->y[i];
		y[i] = 1;

		alpha2[i+l] = 0;
		linear_term[i+l] = param->p + prob->y[i];
		y[i+l] = -1;
	}

	Solver s;
	s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
		alpha2, param->C, param->C, param->eps, si, param->shrinking);

	double sum_alpha = 0;
	for(i=0;i<l;i++)
	{
		alpha[i] = alpha2[i] - alpha2[i+l];
		sum_alpha += fabs(alpha[i]);
	}
	//info("nu = %f\n",sum_alpha/(param->C*l));

	delete[] alpha2;
	delete[] linear_term;
	delete[] y;
}

static void solve_nu_svr(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver::SolutionInfo* si)
{
	int l = prob->l;
	double C = param->C;
	double *alpha2 = new double[2*l];
	double *linear_term = new double[2*l];
	schar *y = new schar[2*l];
	int i;

	double sum = C * param->nu * l / 2;
	for(i=0;i<l;i++)
	{
		alpha2[i] = alpha2[i+l] = min(sum,C);
		sum -= alpha2[i];

		linear_term[i] = - prob->y[i];
		y[i] = 1;

		linear_term[i+l] = prob->y[i];
		y[i+l] = -1;
	}

	Solver_NU s;
	s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
		alpha2, C, C, param->eps, si, param->shrinking);

	//info("epsilon = %f\n",-si->r);

	for(i=0;i<l;i++)
		alpha[i] = alpha2[i] - alpha2[i+l];

	delete[] alpha2;
	delete[] linear_term;
	delete[] y;
}

//
// decision_function
//
struct decision_function
{
	double *alpha;
	double rho;	
};

decision_function svm_train_one(
	const svm_problem *prob, const svm_parameter *param,
	double Cp, double Cn)
{
	double *alpha = Malloc(double,prob->l);
	Solver::SolutionInfo si;
	switch(param->svm_type)
	{
		case C_SVC:
			solve_c_svc(prob,param,alpha,&si,Cp,Cn);
			break;
		case NU_SVC:
			solve_nu_svc(prob,param,alpha,&si);
			break;
		case ONE_CLASS:
			solve_one_class(prob,param,alpha,&si);
			break;
		case EPSILON_SVR:
			solve_epsilon_svr(prob,param,alpha,&si);
			break;
		case NU_SVR:
			solve_nu_svr(prob,param,alpha,&si);
			break;
	}

	//info("obj = %f, rho = %f\n",si.obj,si.rho);

	// output SVs

// 	int nSV = 0;
// 	int nBSV = 0;
// 	for(int i=0;i<prob->l;i++)
// 	{
// 		if(fabs(alpha[i]) > 0)
// 		{
// 			++nSV;
// 			if(prob->y[i] > 0)
// 			{
// 				if(fabs(alpha[i]) >= si.upper_bound_p)
// 					++nBSV;
// 			}
// 			else
// 			{
// 				if(fabs(alpha[i]) >= si.upper_bound_n)
// 					++nBSV;
// 			}
// 		}
// 	}

// 	info("nSV = %d, nBSV = %d\n",nSV,nBSV);

	decision_function f;
	f.alpha = alpha;
	f.rho = si.rho;
	return f;
}

//
// svm_model
//
struct svm_model
{
	svm_parameter param;	// parameter
	int nr_class;		// number of classes, = 2 in regression/one class svm
	int l;			// total #SV
	svm_node **SV;		// SVs (SV[l])
	double **sv_coef;	// coefficients for SVs in decision functions (sv_coef[n-1][l])
	double *rho;		// constants in decision functions (rho[n*(n-1)/2])

	// for classification only

	int *label;		// label of each class (label[n])
	int *nSV;		// number of SVs for each class (nSV[n])
				// nSV[0] + nSV[1] + ... + nSV[n-1] = l
	// XXX
	int free_sv;		// 1 if svm_model is created by svm_load_model
				// 0 if svm_model is created by svm_train
};



const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
{
	// svm_type

	int svm_type = param->svm_type;
	if(svm_type != C_SVC &&
	   svm_type != NU_SVC &&
	   svm_type != ONE_CLASS &&
	   svm_type != EPSILON_SVR &&
	   svm_type != NU_SVR)
		return "unknown svm type";
	
	// kernel_type
	
	int kernel_type = param->kernel_type;
	if(kernel_type != LINEAR &&
	   kernel_type != POLY &&
	   kernel_type != RBF &&
	   kernel_type != SIGMOID &&
	   kernel_type != R)
		return "unknown kernel type";

	// cache_size,eps,C,nu,p,shrinking

	if(param->cache_size <= 0)
		return "cache_size <= 0";

	if(param->eps <= 0)
		return "eps <= 0";

	if(svm_type == C_SVC ||
	   svm_type == EPSILON_SVR ||
	   svm_type == NU_SVR)
		if(param->C <= 0)
			return "C <= 0";

	if(svm_type == NU_SVC ||
	   svm_type == ONE_CLASS ||
	   svm_type == NU_SVR)
		if(param->nu < 0 || param->nu > 1)
			return "nu < 0 or nu > 1";

	if(svm_type == EPSILON_SVR)
		if(param->p < 0)
			return "p < 0";

	if(param->shrinking != 0 &&
	   param->shrinking != 1)
		return "shrinking != 0 and shrinking != 1";


	// check whether nu-svc is feasible
	
	if(svm_type == NU_SVC)
	{
		int l = prob->l;
		int max_nr_class = 16;
		int nr_class = 0;
		int *label = Malloc(int,max_nr_class);
		int *count = Malloc(int,max_nr_class);

		int i;
		for(i=0;i<l;i++)
		{
			int this_label = (int)prob->y[i];
			int j;
			for(j=0;j<nr_class;j++)
				if(this_label == label[j])
				{
					++count[j];
					break;
				}
			if(j == nr_class)
			{
				if(nr_class == max_nr_class)
				{
					max_nr_class *= 2;
					label = (int *)realloc(label,max_nr_class*sizeof(int));
					count = (int *)realloc(count,max_nr_class*sizeof(int));
				}
				label[nr_class] = this_label;
				count[nr_class] = 1;
				++nr_class;
			}
		}
	
		for(i=0;i<nr_class;i++)
		{
			int n1 = count[i];
			for(int j=i+1;j<nr_class;j++)
			{
				int n2 = count[j];
				if(param->nu*(n1+n2)/2 > min(n1,n2))
				{
					free(label);
					free(count);
					return "specified nu is infeasible";
				}
			}
		}
	}

	return NULL; 
}


extern "C" {

#include <R.h>
#include <Rinternals.h>
#include <stdio.h>  

  struct svm_node ** sparsify (double *x, int r, int c)
  {
    struct svm_node** sparse;
    int         i, ii, count;
    
    sparse = (struct svm_node **) malloc (r * sizeof(struct svm_node *));
    for (i = 0; i < r; i++) {
      /* determine nr. of non-zero elements */
      for (count = ii = 0; ii < c; ii++)
	if (x[i * c + ii] != 0) count++;
      
      /* allocate memory for column elements */
      sparse[i] = (struct svm_node *) malloc ((count + 1) * sizeof(struct svm_node));
      
      /* set column elements */
      for (count = ii = 0; ii < c; ii++)
	if (x[i * c + ii] != 0) {
	  sparse[i][count].index = ii;
	  sparse[i][count].value = x[i * c + ii];
	  count++;
	}
      
      /* set termination element */
      sparse[i][count].index = -1;
    }
    
    return sparse;
  }
  
  void solve_smo(const svm_problem *prob, const svm_parameter* param,
		 double *alpha, Solver::SolutionInfo* si, double C, double *linear_term)
  {
    int l = prob->l;
    int i;

    switch(param->svm_type)
      {
      case C_SVC:
	{ double Cp,Cn;
	  double *minus_ones = new double[l];
	  schar *y = new schar[l];
	  for(i=0;i<l;i++)
	    {
	      alpha[i] = 0;
	      minus_ones[i] = -1;
	      if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;
	    }
	  if(param->nr_weight > 0)
	    {
	      Cp = C*param->weight[0];
	      Cn = C*param->weight[1];
	    }
	  else 
	    Cp = Cn = C;
	  Solver s; //have to weight cost parameter for multiclass. problems 
	  s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
		  alpha, Cp, Cn, param->eps, si, param->shrinking);
	  delete[] minus_ones;
	  delete[] y;
	}
	break;
      case NU_SVC:
	{
	  schar *y = new schar[l];
	  double nu = param->nu;
	  double sum_pos = nu*l/2;
	  double sum_neg = nu*l/2;
	  for(i=0;i<l;i++)
	    if(prob->y[i]>0) 
	      { 
		y[i] = +1;  
		alpha[i] = min(1.0,sum_pos);
		sum_pos -= alpha[i];
	      }
	    else {
	      y[i] = -1;
	      alpha[i] = min(1.0,sum_neg);
	      sum_neg -= alpha[i];
	  }
	  double *zeros = new double[l];
	  for(i=0;i<l;i++)
	    zeros[i] = 0;
	  Solver_NU s;
	  s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
		  alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
	  double r = si->r;
	  //info("C = %f\n",1/r);
	  for(i=0;i<l;i++)
	    alpha[i] *= y[i]/r;
	  si->rho /= r;
	  si->obj /= (r*r);
	  si->upper_bound_p = 1/r;
	  si->upper_bound_n = 1/r;
	  delete[] y;
	  delete[] zeros;
	}
	break;
      case ONE_CLASS:
	{
	  double *zeros = new double[l];
	  schar *ones = new schar[l];
	  int n = (int)(param->nu*l);	// # of alpha's at upper bound
	  // set initial alpha probably usefull for smo 
	  for(i=0;i<n;i++)
	    alpha[i] = 1;
	  alpha[n] = param->nu * l - n;
	  for(i=n+1;i<l;i++)
	    alpha[i] = 0;
	  for(i=0;i<l;i++)
	    {
	      zeros[i] = 0;
	      ones[i] = 1;
	    }
	  
	  Solver s;
	  s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
		  alpha, 1.0, 1.0, param->eps, si, param->shrinking);

	  delete[] zeros;
	  delete[] ones;
	}
	  break;
      case EPSILON_SVR:
	{
	  double *alpha2 = new double[2*l];
	  double *linear_term = new double[2*l];
	  schar *y = new schar[2*l];
	  
	  for(i=0;i<l;i++)
	    {
	      alpha2[i] = 0;
	      linear_term[i] = param->p - prob->y[i];
	      y[i] = 1;
	      alpha2[i+l] = 0;
	      linear_term[i+l] = param->p + prob->y[i];
	      y[i+l] = -1;
	    }
	  Solver s;
	  s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
		  alpha2, param->C, param->C, param->eps, si, param->shrinking);
	  double sum_alpha = 0;
	  for(i=0;i<l;i++)
	    {
	      alpha[i] = alpha2[i] - alpha2[i+l];
	      sum_alpha += fabs(alpha[i]);
	    }
	  //info("nu = %f\n",sum_alpha/(param->C*l));
	  
	  delete[] alpha2;
	  delete[] linear_term;
	  delete[] y; 
	}
	break;
		case NU_SVR:
		  {
		    double C = param->C;
		    double *alpha2 = new double[2*l];
		    double *linear_term = new double[2*l];
		    schar *y = new schar[2*l];
		    double sum = C * param->nu * l / 2;
		    for(i=0;i<l;i++)
		      {
			alpha2[i] = alpha2[i+l] = min(sum,C);
			sum -= alpha2[i];
			
			linear_term[i] = - prob->y[i];
			y[i] = 1;
			
			linear_term[i+l] = prob->y[i];
			y[i+l] = -1;
		      }
		    
		    Solver_NU s;
		    s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
			    alpha2, C, C, param->eps, si, param->shrinking);
		    
		    //info("epsilon = %f\n",-si->r);
		    
		    for(i=0;i<l;i++)
		      alpha[i] = alpha2[i] - alpha2[i+l];
		    
		    delete[] alpha2;
		    delete[] linear_term;
		    delete[] y;
		  }
			break;
      }
  }

  SEXP smo_optim(SEXP x,
		 SEXP r, 
		 SEXP c, 
		 SEXP y, 
		 SEXP linear_term, 
		 SEXP kernel_type, 
		 SEXP svm_type, 
		 SEXP cost, 
		 SEXP nu, 
		 SEXP eps, 
		 SEXP gamma, 
		 SEXP degree, 
		 SEXP coef0, 
		 SEXP weightlabels, 
		 SEXP weights, 
		 SEXP nweights, 
		 SEXP cache,
		 SEXP epsilon, 
		 SEXP shrinking,
		 SEXP expr,
		 SEXP rho
		 )
  {
    
    SEXP alpha;
    double *alpha2;
    struct svm_parameter param;
    struct svm_problem  prob;
    int i;  
    const char* s;
    struct Solver::SolutionInfo si;
    param.svm_type    = *INTEGER(svm_type);
    param.kernel_type = *INTEGER(kernel_type); 
    param.degree      = *REAL(degree); 
    param.gamma       = *REAL(gamma);
    param.coef0       = *REAL(coef0);
    param.cache_size  = *REAL(cache);
    param.eps         = *REAL(epsilon);
    param.C           = *REAL(cost);
    param.nu          = *REAL(nu);
    param.expr        = expr;
    param.rho         = rho;
    param.nr_weight   = *INTEGER(nweights);
    if (param.nr_weight > 0) {
      param.weight      = (double *) malloc (sizeof(double) * param.nr_weight);
      memcpy (param.weight, REAL(weights), param.nr_weight * sizeof(double));
      param.weight_label = (int *) malloc (sizeof(int) * param.nr_weight);
      memcpy (param.weight_label, INTEGER(weightlabels), param.nr_weight * sizeof(int));
    }
    param.p           = *REAL(eps);
    param.shrinking   = *INTEGER(shrinking);

    /* set problem */
    prob.l = *INTEGER(r);
    // prob.y =  (double *) malloc (sizeof(double) * prob.l);
    prob.y = REAL(y);

    prob.x = sparsify(REAL(x), *INTEGER(r), *INTEGER(c)); 
    alpha2      = (double *) malloc (sizeof(double) * prob.l);
    s = svm_check_parameter(&prob, &param);
    if (s) {
      printf("%s",s);
    } else {

            solve_smo(&prob, &param, alpha2, &si, *REAL(cost), REAL(linear_term));
    }
    PROTECT(alpha = allocVector(REALSXP, prob.l+1));
    /* clean up memory */
    if (param.nr_weight > 0) {
      free(param.weight);
      free(param.weight_label);
    }
    for (i = 0; i < prob.l; i++) {free (prob.x[i]); REAL(alpha)[i] = *(alpha2+i); } 
    free (prob.x);
    REAL(alpha)[prob.l]=si.rho;
    free(alpha2); 
    UNPROTECT(1);  
    return alpha;
  }

  void smo_setup(double *x,int *r, int *c, double *y, double *linear_term, int *kernel_type, 
		 int *svm_type, double *cost, double *nu, double *gamma,double *eps,int degree,
		 double *coef0, int *weightlabels, double *weights, int *nweights, 
		 double *cache, double *epsilon, int *shrinking, double *alpha, double *beta, char  **error) 
  {
    struct svm_parameter param;
    struct svm_problem  prob;
    int i;  const char* s;
    struct Solver::SolutionInfo si;
    /* set parameters */
    param.svm_type    = *svm_type;
    param.kernel_type = *kernel_type;
    param.degree      = degree;
    param.gamma       = *gamma;
    param.coef0       = *coef0;
    param.cache_size  = *cache;
    param.eps         = *eps;
    param.C           = *cost;
    param.nu          = *nu;
    param.nr_weight   = *nweights;
    if (param.nr_weight > 0) {
      param.weight      = (double *) malloc (sizeof(double) * param.nr_weight);
      memcpy (param.weight, weights, param.nr_weight * sizeof(double));
      param.weight_label = (int *) malloc (sizeof(int) * param.nr_weight);
      memcpy (param.weight_label, weightlabels, param.nr_weight * sizeof(int));
    }
    param.p           = *epsilon;
    param.shrinking   = *shrinking;
    /* set problem */
    prob.l = *r;
    prob.y = y;
    prob.x = sparsify(x, *r, *c);
     s = svm_check_parameter(&prob, &param);
    if (s) {
	strcpy(*error, s);
    } else {
	/* call smo  */
    solve_smo(&prob, &param, alpha, &si,*cost, linear_term);
    }

    *beta = si.rho;
    /* clean up memory */
    if (param.nr_weight > 0) {
      free(param.weight);
      free(param.weight_label);
    }
    for (i = 0; i < *r; i++) free (prob.x[i]);
    free (prob.x);
  }
}