bsvm.cpp

一个support vector machine的matlab工具箱。

		swap_index(start1[i], start1[i-1]);
	t = nr_class*nr_class;
	for (i=s+1;i<=t;i++)
		start1[i]++;
	for (i=0;i<=s;i++)
		start2[i]++;
}

//
// Q matrices for various formulations
//
class BSVC_Q: public Kernel
{ 
public:
	BSVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
	:Kernel(prob.l, prob.x, param)
	{
		clone(y,y_,prob.l);
		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int start;
		if((start = cache->get_data(i,&data,len)) < len)
		{
			for(int j=start;j<len;j++)
				data[j] = (Qfloat)y[i]*y[j]*((this->*kernel_function)(i,j) + 1);
		}
		return data;
	}

	void swap_index(int i, int j) const
	{
		cache->swap_index(i,j);
		Kernel::swap_index(i,j);
		swap(y[i],y[j]);
	}

	~BSVC_Q()
	{
		delete[] y;
		delete cache;
	}
private:
	schar *y;
	Cache *cache;
};

class ONE_CLASS_Q: public Kernel
{
public:
	ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
	:Kernel(prob.l, prob.x, param)
	{
		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int start;
		if((start = cache->get_data(i,&data,len)) < len)
		{
			for(int j=start;j<len;j++)
				data[j] = (Qfloat)(this->*kernel_function)(i,j);
		}
		return data;
	}

	void swap_index(int i, int j) const
	{
		cache->swap_index(i,j);
		Kernel::swap_index(i,j);
	}

	~ONE_CLASS_Q()
	{
		delete cache;
	}
private:
	Cache *cache;
};

class BONE_CLASS_Q: public Kernel
{
public:
	BONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
	:Kernel(prob.l, prob.x, param)
	{
		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int start;
		if((start = cache->get_data(i,&data,len)) < len)
		{
			for(int j=start;j<len;j++)
				data[j] = (Qfloat)(this->*kernel_function)(i,j) + 1;
		}
		return data;
	}

	~BONE_CLASS_Q()
	{
		delete cache;
	}
private:
	Cache *cache;
};

class BSVR_Q: public Kernel
{ 
public:
	BSVR_Q(const svm_problem& prob, const svm_parameter& param)
	:Kernel(prob.l, prob.x, param)
	{
		l = prob.l;
		cache = new Cache(l,(int)(param.cache_size*(1<<20)));
		sign = new schar[2*l];
		index = new int[2*l];
		for(int k=0;k<l;k++)
		{
			sign[k] = 1;
			sign[k+l] = -1;
			index[k] = k;
			index[k+l] = k;
		}
		q = param.qpsize;
		buffer = new Qfloat*[q];
		for (int i=0;i<q;i++)
			buffer[i] = new Qfloat[2*l];	
		next_buffer = 0;
	}

	void swap_index(int i, int j) const
	{
		swap(sign[i],sign[j]);
		swap(index[i],index[j]);
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int real_i = index[i];
		if(cache->get_data(real_i,&data,l) < l)
		{
			for(int j=0;j<l;j++)
				data[j] = (Qfloat)(this->*kernel_function)(real_i,j) + 1;
		}

		// reorder and copy
		Qfloat *buf = buffer[next_buffer];
		next_buffer = (next_buffer+1)%q;
		schar si = sign[i];
		for(int j=0;j<len;j++)
			buf[j] = si * sign[j] * data[index[j]];
		return buf;
	}

	~BSVR_Q()
	{
		delete cache;
		delete[] sign;
		delete[] index;
		for (int i=0;i<q;i++)
			delete[] buffer[i];
		delete[] buffer;
	}
private:
	int l, q;
	Cache *cache;
	schar *sign;
	int *index;
	mutable int next_buffer;
	Qfloat** buffer;
};

//
// construct and solve various formulations
//
static void solve_c_svc(
	const svm_problem *prob, const svm_parameter* param,
	double *alpha, Solver_B::SolutionInfo* si, double Cp, double Cn)
{
	int l = prob->l;
	double *minus_ones = new double[l];
	schar *y = new schar[l];

	int i;

	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, si, 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, 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;
} 

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

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

	// 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[0])
					++nBSV;
			}
			else
			{
				if (fabs(alpha[i]) >= si.upper_bound[1])
					++nBSV;
			}
		}
	}
	info("nSV = %d, nBSV = %d\n",nSV,nBSV);

	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
//
svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
{
	svm_model *model = Malloc(svm_model,1);
	model->param = *param;
	model->free_sv = 0;	// XXX

	if (param->svm_type == EPSILON_SVR)
	{

		// regression
		model->nr_class = 2;
		model->label = NULL;
		model->nSV = NULL;
		model->sv_coef = Malloc(double *,1);
		decision_function f = svm_train_one(prob,param,0,0);

		int nSV = 0;
		int i;
		for(i=0;i<prob->l;i++)
			if(fabs(f.alpha[i]) > 0) ++nSV;
		model->l = nSV;
		model->SV = Malloc(svm_node *,nSV);
		model->sv_coef[0] = Malloc(double,nSV);
		int j = 0;
		for(i=0;i<prob->l;i++)
			if(fabs(f.alpha[i]) > 0)
			{
				model->SV[j] = prob->x[i];
				model->sv_coef[0][j] = f.alpha[i];
				++j;
			}		

		free(f.alpha);

	}
	else
	{

		// classification
		// find out the number of classes
		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 *index = Malloc(int,l);

		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;
				}
			index[i] = j;
			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;
			}
		}

		// group training data of the same class

		int *start = Malloc(int,nr_class);
		start[0] = 0;
		for(i=1;i<nr_class;i++)
			start[i] = start[i-1]+count[i-1];

		svm_node **x = Malloc(svm_node *,l);
		
		for(i=0;i<l;i++)
		{
			x[start[index[i]]] = prob->x[i];
			++start[index[i]];
		}
		
		start[0] = 0;
		for(i=1;i<nr_class;i++)
			start[i] = start[i-1]+count[i-1];

		// calculate weighted C

		double *weighted_C = Malloc(double, nr_class);
		for(i=0;i<nr_class;i++)
			weighted_C[i] = param->C;
		for(i=0;i<param->nr_weight;i++)
		{	
			int j;
			for(j=0;j<nr_class;j++)
				if(param->weight_label[i] == label[j])
					break;
			if(j == nr_class)
				fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
			else
				weighted_C[j] *= param->weight[i];
		}

		bool *nonzero = Malloc(bool,l);
		for(i=0;i<l;i++)
			nonzero[i] = false;

	if (param->svm_type == C_SVC)
	{	

		// train n*(n-1)/2 models
		
		decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);

		int p = 0;
		for(i=0;i<nr_class;i++)
			for(int j=i+1;j<nr_class;j++)
			{
				svm_problem sub_prob;
				int si = start[i], sj = start[j];
				int ci = count[i], cj = count[j];
				sub_prob.l = ci+cj;