bsvm.cpp

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

				for (i=0;i<=prob->n;i++)
					w[i] = 0;
				for (i=0;i<l;i++)
				{
					if (y[i] == 1 && alpha[i] >= Cpj)
						alpha[i] = Cp;
					else 
						if (y[i] == -1 && alpha[i] >= Cnj)
							alpha[i] = Cn;
						else
							alpha[i] *= Cp/Cpj;
					double yalpha = y[i]*alpha[i];
					for (const svm_node *px = prob->x[i];px->index != -1;px++)
						w[px->index] += yalpha*px->value;
					w[0] += yalpha;
				}
			}
			else
			{
				for (i=0;i<l;i++)
					alpha[i] *= param->Cstep;
				for (i=0;i<=prob->n;i++)
					w[i] *= param->Cstep;
			}
			Cpj *= param->Cstep;
			Cnj *= param->Cstep;
		}
		totaliter += s.Solve(l, prob->x, minus_ones, y, alpha, w, Cp, Cn,
		param->eps, si, param->shrinking, param->qpsize);
		//info("\noptimization finished, #iter = %d\n",totaliter);
	}
	else
	{
		Solver_B s;
		s.Solve(l, BSVC_Q(*prob,*param,y), minus_ones, y, alpha, Cp, Cn, 
		param->eps, si, param->shrinking, param->qpsize);
	}

	double sum_alpha=0;
	for(i=0;i<l;i++)
		sum_alpha += alpha[i];

	//info("nu = %f\n", sum_alpha/(param->C*prob->l));

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

	delete[] minus_ones;
	delete[] y;
}

static void solve_epsilon_svr(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver_B::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;
	}

	if (param->kernel_type == LINEAR)
	{
		double *w = new double[prob->n+1];
		for (i=0;i<=prob->n;i++)
			w[i] = 0;
		struct svm_node **x = new svm_node*[2*l];
		for (i=0;i<l;i++)
			x[i] = x[i+l] = prob->x[i];
		Solver_B_linear s;
		int totaliter = 0;
		double Cj = param->Cbegin;
		while (Cj < param->C)
		{
			totaliter += s.Solve(2*l, x, linear_term, y, alpha2, w, 
			Cj, Cj, param->eps, si, param->shrinking, param->qpsize);
			if (Cj*param->Cstep >= param->C)
			{
				for (i=0;i<=prob->n;i++)
					w[i] = 0;
				for (i=0;i<2*l;i++)
				{
					if (alpha2[i] >= Cj)
						alpha2[i] = param->C;
					else 
						alpha2[i] *= param->C/Cj;
					double yalpha = y[i]*alpha2[i];
					for (const svm_node *px = x[i];px->index != -1;px++)
						w[px->index] += yalpha*px->value;
					w[0] += yalpha;
				}
			}
			else
			{
				for (i=0;i<2*l;i++)
					alpha2[i] *= param->Cstep;
				for (i=0;i<=prob->n;i++)
					w[i] *= param->Cstep;
			}
			Cj *= param->Cstep;
		}
		totaliter += s.Solve(2*l, x, linear_term, y, alpha2, w, param->C,
			param->C, param->eps, si, param->shrinking, param->qpsize);
		//info("\noptimization finished, #iter = %d\n",totaliter);
	}
	else
	{
		Solver_B s;
		s.Solve(2*l, BSVR_Q(*prob,*param), linear_term, y, alpha2, param->C,
			param->C, param->eps, si, param->shrinking, param->qpsize);
	}

	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[] y;
	delete[] alpha2;
	delete[] linear_term;
}

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

decision_function solve_spoc(const svm_problem *prob, 
	const svm_parameter *param, int nr_class, double *weighted_C)
{
	int i, m, l = prob->l;
	double *alpha = new double[l*nr_class];
	short *y = new short[l];
	
	for (i=0;i<l;i++)
	{
		for (m=0;m<nr_class;m++)
			alpha[i*nr_class+m] = 0;
		y[i] = (short) prob->y[i];	
	}
	
	Solver_SPOC s;
	s.Solve(l, ONE_CLASS_Q(*prob, *param), alpha, y, weighted_C,
	param->eps, param->shrinking, nr_class);
	
	delete[] y;
	decision_function f;
	f.alpha = alpha;
	return f;  	
}

decision_function solve_msvm(const svm_problem *prob, 
	const svm_parameter *param, int nr_class, double *weighted_C, int *count)
{
	Solver_B::SolutionInfo si;
	int i, l = prob->l*(nr_class - 1);
	double *alpha = Malloc(double, l);
	short *y = new short[l];
	
	for (i=0;i<l;i++)
		alpha[i] = 0;
	
	int p = 0;
	for (i=0;i<nr_class;i++)
	{
		int q = 0;
		for (int j=0;j<nr_class;j++)
			if (i != j)
				for (int k=0;k<count[j];k++,q++)
					y[p++] = (short) prob->y[q];
			else
				q += count[j];
	}

	Solver_MB s;
	s.Solve(l, BONE_CLASS_Q(*prob,*param), -2, alpha, y, weighted_C,
	2*param->eps, &si, param->shrinking, param->qpsize, nr_class, count);

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

        // output SVs

	int nSV = 0, nBSV = 0;
	p = 0;
	for (i=0;i<nr_class;i++)
	{
		int q = 0;
		for (int j=0;j<nr_class;j++)
			if (i != j)
				for (int k=0;k<count[j];k++,q++)
				{
					if (fabs(alpha[p]) > 0)
					{
						++nSV;
						if (fabs(alpha[p]) >= si.upper_bound[(int) prob->y[q]])
							++nBSV;
					}
					p++;
				}
			else
				q += count[j];
	}
	//info("nSV = %d, nBSV = %d\n",nSV,nBSV);

	delete[] y;
	delete[] si.upper_bound;
	decision_function f;
	f.alpha = alpha;
	return f;
} 


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

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

//
// Interface functions
//


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

	int svm_type = param->svm_type;
	if(svm_type != C_SVC &&
	   svm_type != EPSILON_SVR &&
	   svm_type != KBB &&
	   svm_type != SPOC)
		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)
		return "unknown kernel type";

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

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

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

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

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

	if(svm_type == C_SVC ||
	   svm_type == KBB ||
	   svm_type == SPOC)
		if(param->qpsize < 2)
			return "qpsize < 2";

	if(kernel_type == LINEAR)
		if (param->Cbegin <= 0)
			return "Cbegin <= 0";

	if(kernel_type == LINEAR)
		if (param->Cstep <= 1)
			return "Cstep <= 1";

	return NULL;
}



extern "C" {

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

  struct svm_node ** sparsifyb (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 lala3(const svm_problem *prob, const svm_parameter* param, 
		  double *alpha,  double *weighted_C, Solver_B::SolutionInfo* sii, int nr_class, int *count)
  {
    int l = prob->l;
    int i;
    double Cp = param->C;
    double Cn = param->C;
     
    switch(param->svm_type)
      {
      case C_SVC:
	{
	  double *alpha = new double[l];
	  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->kernel_type == LINEAR)
	{
		double *w = new double[prob->n+1];
		for (i=0;i<=prob->n;i++)
			w[i] = 0;
		Solver_B_linear s;
		int totaliter = 0;
		double Cpj = param->Cbegin, Cnj = param->Cbegin*Cn/Cp;
		while (Cpj < Cp)
		{
			totaliter += s.Solve(l, prob->x, minus_ones, y, alpha, w, 
			Cpj, Cnj, param->eps, sii, param->shrinking, param->qpsize);
			if (Cpj*param->Cstep >= Cp)
			{
			  for (i=0;i<=prob->n;i++)
					w[i] = 0;
				for (i=0;i<l;i++)
				{
					if (y[i] == 1 && alpha[i] >= Cpj)
						alpha[i] = Cp;
					else 
						if (y[i] == -1 && alpha[i] >= Cnj)
							alpha[i] = Cn;
						else
							alpha[i] *= Cp/Cpj;
					double yalpha = y[i]*alpha[i];
					for (const svm_node *px = prob->x[i];px->index != -1;px++)
						w[px->index] += yalpha*px->value;
					w[0] += yalpha;
				}
			}
			else
			{
				for (i=0;i<l;i++)
					alpha[i] *= param->Cstep;
				for (i=0;i<=prob->n;i++)
					w[i] *= param->Cstep;
			}
			Cpj *= param->Cstep;
			Cnj *= param->Cstep;
		}
		totaliter += s.Solve(l, prob->x, minus_ones, y, alpha, w, Cp, Cn,
				     param->eps, sii, param->shrinking, param->qpsize);
		//info("\noptimization finished, #iter = %d\n",totaliter);


	}
	}
	break;
      case EPSILON_SVR:
	{
	  double *alpha = new double[2*l];
	double *linear_term = new double[2*l];
	schar *y = new schar[2*l];

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

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

	if (param->kernel_type == LINEAR)
	{
		double *w = new double[prob->n+1];
		for (i=0;i<=prob->n;i++)
			w[i] = 0;
		struct svm_node **x = new svm_node*[2*l];
		for (i=0;i<l;i++)
			x[i] = x[i+l] = prob->x[i];
		Solver_B_linear s;
		int totaliter = 0;
		double Cj = param->Cbegin;
		while (Cj < param->C)
		{
			totaliter += s.Solve(2*l, x, linear_term, y, alpha, w, 
			Cj, Cj, param->eps, sii, param->shrinking, param->qpsize);
			if (Cj*param->Cstep >= param->C)
			{
				for (i=0;i<=prob->n;i++)
					w[i] = 0;
				for (i=0;i<2*l;i++)
				{
					if (alpha[i] >= Cj)
						alpha[i] = param->C;
					else 
						alpha[i] *= param->C/Cj;
					double yalpha = y[i]*alpha[i];
					for (const svm_node *px = x[i];px->index != -1;px++)
						w[px->index] += yalpha*px->value;
					w[0] += yalpha;
				}
			}
			else
			{
				for (i=0;i<2*l;i++)
					alpha[i] *= param->Cstep;
				for (i=0;i<=prob->n;i++)
					w[i] *= param->Cstep;
			}
			Cj *= param->Cstep;
		}
		totaliter += s.Solve(2*l, x, linear_term, y, alpha, w, param->C,
			param->C, param->eps, sii, param->shrinking, param->qpsize);
		//info("\noptimization finished, #iter = %d\n",totaliter);

	}
	}
	break;
      case KBB:
	{
	  Solver_B::SolutionInfo si;
	  int i, ll = l*(nr_class - 1);
	  //  double *alpha = Malloc(double, l);
	  short *y = new short[ll];
	  
	  for (i=0;i<ll;i++)
	    alpha[i] = 0;
	  
	  int p = 0;
	  for (i=0;i<nr_class;i++)
	    {
	      int q = 0;
	      for (int j=0;j<nr_class;j++)
		if (i != j)
		  for (int k=0;k<count[j];k++,q++)
		    y[p++] = (short) prob->y[q];
		else
		  q += count[j];
	    }
	  
	  	  Solver_MB s;
	   s.Solve(ll, BONE_CLASS_Q(*prob,*param), -2, alpha, y, weighted_C,
	  	  2*param->eps, &si, param->shrinking, param->qpsize, nr_class, count);
	  
	   //info("obj = %f, rho = %f\n",si.obj,0.0);
	  
	  // output SVs
	  
	  int nSV = 0, nBSV = 0;
	  p = 0;
	for (i=0;i<nr_class;i++)
	  {
	    int q = 0;
	    for (int j=0;j<nr_class;j++)
	      if (i != j)
		for (int k=0;k<count[j];k++,q++)
		  {
		    if (fabs(alpha[p]) > 0)
		      {
			++nSV;
			if (fabs(alpha[p]) >= si.upper_bound[(int) prob->y[q]])
			  ++nBSV;
		      }
		    p++;
		  }
			else
			  q += count[j];
	  }
	//info("nSV = %d, nBSV = %d\n",nSV,nBSV);

	delete[] y;
	delete[] si.upper_bound;
	}
	  break;
      case SPOC:
	{
	  int m;
	  short *y = new short[l];
	  for (i=0;i<l;i++)
	    {
	      for (m=0;m<nr_class;m++)
	  	alpha[i*nr_class+m] = 0;
	      y[i] = (short) prob->y[i];	
	    }

	  Solver_SPOC s;
	  s.Solve(l, ONE_CLASS_Q(*prob, *param), alpha, y, weighted_C,
		  param->eps, param->shrinking, nr_class);
	  free(weighted_C);
	  delete[] y;
	}
	break;
      }
  }	

  void lala2(svm_problem *prob, svm_parameter* param,  double *alpha,
	     double *weighted_C, Solver_B::SolutionInfo* si, int nr_class)
  {
    int l = prob->l;
    int i;
   
    	  int m;
	  //	  	  double *alpha = new double[l*nr_class];
	  short *y = new short[l];
	   for (i=0;i<l;i++)
	   {
	     for (m=0;m<nr_class;m++)
	  	alpha[i*nr_class+m] = 0;
	     y[i] = (short) prob->y[i];	
	   }
	  
	   Solver_SPOC s;
	   s.Solve(l, ONE_CLASS_Q(*prob, *param), alpha, y, weighted_C,
		   param->eps, param->shrinking, nr_class);
	   free(weighted_C);
	delete[] y;
  }

  SEXP tron_optim(SEXP x,
		  SEXP r, 
		  SEXP c, 
		  SEXP y, 
		  SEXP nclass,
		  SEXP countc,
		  SEXP kernel_type, 
		  SEXP svm_type, 
		  SEXP cost, 
		  SEXP eps, 
		  SEXP gamma, 
		  SEXP degree, 
		  SEXP coef0,
		  SEXP Cbegin,
		  SEXP Cstep,
		  SEXP weightlabels, 
		  SEXP weights, 
		  SEXP nweights,
		  SEXP weightedc,
		  SEXP cache,
		  SEXP epsilon, 
		  SEXP qpsize,
		  SEXP shrinking,
		  SEXP expr,
		  SEXP rho
		 )
  {
    
    struct svm_parameter param;
    struct svm_problem  prob;
    int i ,*count = NULL;
    double *alpha2 = NULL;
    SEXP alpha3 = NULL;
    int nr_class;
    const char* s;
    struct Solver_B::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.Cbegin      = *REAL(Cbegin);
    param.Cstep       = *REAL(Cstep);
    param.qpsize      = *INTEGER(qpsize);
    param.expr        = expr;
    param.rho         = rho;
    nr_class          = *INTEGER(nclass);
    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);
    memcpy(prob.y, REAL(y), prob.l*sizeof(double));
    prob.x = sparsifyb(REAL(x), *INTEGER(r), *INTEGER(c)); 
    s = svm_check_parameterb(&prob, &param);
   
    if (s) 
      printf("%s",s);
    else {
      double *weighted_C = Malloc(double, nr_class);
      memcpy(weighted_C, REAL(weightedc), nr_class*sizeof(double));
      
      if(param.svm_type == 2)
	{ 
	 alpha2 = (double *) malloc (sizeof(double) * prob.l*nr_class);
   	}
      //if(param.svm_type == 1)
      else 
	{ 
	  count = Malloc(int, nr_class);
	  memcpy(count, REAL(countc), nr_class*sizeof(int));
	  alpha2 = (double *) malloc (sizeof(double) * prob.l*(nr_class-1));
	}
      //           svm_model *l=svm_train(&prob,&param);
	          lala3(&prob, &param, alpha2, weighted_C , &si, nr_class, count); 
      //      lala2(&prob, &param, alpha2, weighted_C , &si, nr_class); 
    }
    
    /* clean up memory */
    if (param.nr_weight > 0) {
      free(param.weight);
      free(param.weight_label);
    }
   
     if(param.svm_type == 2)
      {  
	PROTECT(alpha3 = allocVector(REALSXP, (nr_class*prob.l)));
	UNPROTECT(1);  
	for (i = 0; i < prob.l; i++) 
	  free (prob.x[i]); 
	for (i = 0; i <nr_class*prob.l; i++) 
	  REAL(alpha3)[i] = *(alpha2+i); 
      }
     
     if(param.svm_type ==1)
       { 
	 PROTECT(alpha3 = allocVector(REALSXP, ((nr_class-1)*prob.l)));
	 UNPROTECT(1);   
	 free(count);
	 for (i = 0; i < prob.l; i++) 
	   free (prob.x[i]); 
	 for (i = 0; i <(nr_class-1)*prob.l; i++) 
	   REAL(alpha3)[i] = *(alpha2+i); 
       }
     if(param.svm_type == 0||param.svm_type== 3)
      { 
	PROTECT(alpha3 = allocVector(REALSXP, prob.l+1));
	UNPROTECT(1);  
	for (i = 0; i < prob.l; i++) {free (prob.x[i]); REAL(alpha3)[i] = *(alpha2+i); } 
      }
    free(prob.x);
    free(prob.y);
    free(alpha2);
    return alpha3;
  }
 }