svm.cpp

采用c++编写的有关svm交叉验证的程序。

		svm_destroy_model(submodel);
		free(subprob.x);
		free(subprob.y);
	}		
	free(fold_start);
	free(perm);	
}


int svm_get_svm_type(const svm_model *model)
{
	return model->param.svm_type;
}

int svm_get_nr_class(const svm_model *model)
{
	return model->nr_class;
}

void svm_get_labels(const svm_model *model, int* label)
{
	if (model->label != NULL)
		for(int i=0;i<model->nr_class;i++)
			label[i] = model->label[i];
}

double svm_get_svr_probability(const svm_model *model)
{
	if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
	    model->probA!=NULL)
		return model->probA[0];
	else
	{
		//info("Model doesn't contain information for SVR probability inference\n");
		info+="模型不包括ＳＶＲ概率推论的信息\n";
		return 0;
	}
}

void svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
{
	if(model->param.svm_type == ONE_CLASS ||
	   model->param.svm_type == EPSILON_SVR ||
	   model->param.svm_type == NU_SVR)
	{
		double *sv_coef = model->sv_coef[0];
		double sum = 0;
		for(int i=0;i<model->l;i++)
			sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
		sum -= model->rho[0];
		*dec_values = sum;
	}
	else
	{
		int i;
		int nr_class = model->nr_class;
		int l = model->l;
		
		double *kvalue = Malloc(double,l);
		for(i=0;i<l;i++)
			kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);

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

		int p=0;
		for(i=0;i<nr_class;i++)
			for(int j=i+1;j<nr_class;j++)
			{
				double sum = 0;
				int si = start[i];
				int sj = start[j];
				int ci = model->nSV[i];
				int cj = model->nSV[j];
				
				int k;
				double *coef1 = model->sv_coef[j-1];
				double *coef2 = model->sv_coef[i];
				for(k=0;k<ci;k++)
					sum += coef1[si+k] * kvalue[si+k];
				for(k=0;k<cj;k++)
					sum += coef2[sj+k] * kvalue[sj+k];
				sum -= model->rho[p];
				dec_values[p] = sum;
				p++;
			}

		free(kvalue);
		free(start);
	}
}

double svm_predict(const svm_model *model, const svm_node *x)
{
	if(model->param.svm_type == ONE_CLASS ||
	   model->param.svm_type == EPSILON_SVR ||
	   model->param.svm_type == NU_SVR)
	{
		double res;
		svm_predict_values(model, x, &res);
		
		if(model->param.svm_type == ONE_CLASS)
			return (res>0)?1:-1;
		else
			return res;
	}
	else
	{
		int i;
		int nr_class = model->nr_class;
		double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
		svm_predict_values(model, x, dec_values);

		int *vote = Malloc(int,nr_class);
		for(i=0;i<nr_class;i++)
			vote[i] = 0;
		int pos=0;
		for(i=0;i<nr_class;i++)
			for(int j=i+1;j<nr_class;j++)
			{
				if(dec_values[pos++] > 0)
					++vote[i];
				else
					++vote[j];
			}

		int vote_max_idx = 0;
		for(i=1;i<nr_class;i++)
			if(vote[i] > vote[vote_max_idx])
				vote_max_idx = i;
		free(vote);
		free(dec_values);
		return model->label[vote_max_idx];
	}
}

double svm_predict_probability(
	const svm_model *model, const svm_node *x, double *prob_estimates)
{
	if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
	    model->probA!=NULL && model->probB!=NULL)
	{
		int i;
		int nr_class = model->nr_class;
		double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
		svm_predict_values(model, x, dec_values);

		double min_prob=1e-7;
		double **pairwise_prob=Malloc(double *,nr_class);
		for(i=0;i<nr_class;i++)
			pairwise_prob[i]=Malloc(double,nr_class);
		int k=0;
		for(i=0;i<nr_class;i++)
			for(int j=i+1;j<nr_class;j++)
			{
				pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
				pairwise_prob[j][i]=1-pairwise_prob[i][j];
				k++;
			}
		multiclass_probability(nr_class,pairwise_prob,prob_estimates);

		int prob_max_idx = 0;
		for(i=1;i<nr_class;i++)
			if(prob_estimates[i] > prob_estimates[prob_max_idx])
				prob_max_idx = i;
		for(i=0;i<nr_class;i++)
			free(pairwise_prob[i]);
		free(dec_values);
                free(pairwise_prob);	     
		return model->label[prob_max_idx];
	}
	else 
		return svm_predict(model, x);
}

const char *svm_type_table[] =
{
	"c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
};

const char *kernel_type_table[]=
{
	"linear","polynomial","rbf","sigmoid","precomputed",NULL
};

int svm_save_model(const char *model_file_name, const svm_model *model)
{
	FILE *fp = fopen(model_file_name,"w");
	if(fp==NULL) return -1;
    //CString str,str1;//添加
	const svm_parameter& param = model->param;

	fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
	//str1.Format("svm类型 %s\n", svm_type_table[param.svm_type]);//添加
	//str+=str1;//添加
	fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);
	//str1.Format("核函数类型 %s\n", kernel_type_table[param.kernel_type]);//添加
	//str+=str1;//添加
	if(param.kernel_type == POLY)
	{
		fprintf(fp,"degree %d\n", param.degree);
		//str1.Format("核函数度 %d\n", param.degree);//添加
	    //str+=str1;//添加
	}
	if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
	{
		fprintf(fp,"gamma %g\n", param.gamma);
		//str1.Format("核函数γ %g\n", param.gamma);//添加
	    //str+=str1;//添加
	}
	if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
	{
		fprintf(fp,"coef0 %g\n", param.coef0);
		//str1.Format("核函数coef0 %g\n", param.coef0);//添加
	    //str+=str1;//添加
	}
	int nr_class = model->nr_class;
	int l = model->l;
	fprintf(fp, "nr_class %d\n", nr_class);
	//str1.Format("类的数目 %d\n", nr_class);//添加
	//str+=str1;//添加
	fprintf(fp, "total_sv %d\n",l);
	//str1.Format("支持向量总个数 %d\n",l);//添加
	//str+=str1;//添加	

	{
		fprintf(fp, "rho");
		//str+="判决函数的常数项b";//添加
		for(int i=0;i<nr_class*(nr_class-1)/2;i++)
		{
			fprintf(fp," %g",model->rho[i]);
		    //str1.Format(" %g",model->rho[i]);//添加
	        //str+=str1;//添加
		}
		fprintf(fp, "\n");
		//str+="\n";//添加
	}
	
	if(model->label)
	{
		fprintf(fp, "label");
		//str+="类";//添加
		for(int i=0;i<nr_class;i++)
		{
			fprintf(fp," %d",model->label[i]);
		    //str1.Format(" %d",model->label[i]);//添加
	        //str+=str1;//添加
	    }
		fprintf(fp, "\n");
		//str+="\n";//添加
	}

	if(model->probA) // regression has probA only
	{
		fprintf(fp, "probA");
		//str+="对偶A概率";//添加
		for(int i=0;i<nr_class*(nr_class-1)/2;i++)
		{
			fprintf(fp," %g",model->probA[i]);
		    //str1.Format(" %g",model->probA[i]);//添加
	        //str+=str1;//添加
		}
		fprintf(fp, "\n");
		//str+="\n";//添加
	}
	if(model->probB)
	{
		fprintf(fp, "probB");
		//str+="对偶B概率";//添加
		for(int i=0;i<nr_class*(nr_class-1)/2;i++)
		{
			fprintf(fp," %g",model->probB[i]);
		    //str1.Format(" %g",model->probB[i]);//添加
	        //str+=str1;//添加
		}
		fprintf(fp, "\n");
		//str+="\n";//添加
	}

	if(model->nSV)
	{
		fprintf(fp, "nr_sv");
		//str+="支持向量数目";//添加
		for(int i=0;i<nr_class;i++)
		{
			fprintf(fp," %d",model->nSV[i]);
		    //str1.Format(" %d",model->nSV[i]);//添加
	        //str+=str1;//添加
		}
		fprintf(fp, "\n");
		//str+="\n";//添加
	}

	fprintf(fp, "SV\n");
	//str+="支持向量\n";//添加
	const double * const *sv_coef = model->sv_coef;
	const svm_node * const *SV = model->SV;

	for(int i=0;i<l;i++)
	{
		for(int j=0;j<nr_class-1;j++)
		{
			fprintf(fp, "%.16g ",sv_coef[j][i]);
		    //str1.Format("%.16g ",sv_coef[j][i]);//添加
	        //str+=str1;//添加
		}
		const svm_node *p = SV[i];

		if(param.kernel_type == PRECOMPUTED)
		{
			fprintf(fp,"0:%d ",(int)(p->value));
		    //str1.Format("0:%d ",(int)(p->value));//添加
	        //str+=str1;//添加
		}
		else
			while(p->index != -1)
			{
				fprintf(fp,"%d:%.8g ",p->index,p->value);
				//str1.Format("%d:%.8g ",p->index,p->value);//添加
				//str+=str1;//添加
				p++;
			}
		fprintf(fp, "\n");
		//str+="\n";//添加
	}

	if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
	else return 0;
	//AfxMessageBox(str);
}

svm_model *svm_load_model(const char *model_file_name)
{
	FILE *fp = fopen(model_file_name,"rb");
	if(fp==NULL) return NULL;
	
	// read parameters

	svm_model *model = Malloc(svm_model,1);
	svm_parameter& param = model->param;
	model->rho = NULL;
	model->probA = NULL;
	model->probB = NULL;
	model->label = NULL;
	model->nSV = NULL;

	char cmd[81];
	while(1)
	{
		fscanf(fp,"%80s",cmd);

		if(strcmp(cmd,"svm_type")==0)
		{
			fscanf(fp,"%80s",cmd);
			int i;
			for(i=0;svm_type_table[i];i++)
			{
				if(strcmp(svm_type_table[i],cmd)==0)
				{
					param.svm_type=i;
					break;
				}
			}
			if(svm_type_table[i] == NULL)
			{
				//fprintf(stderr,"unknown svm type.\n");
				AfxMessageBox("未知的svm类型\n");
				free(model->rho);
				free(model->label);
				free(model->nSV);
				free(model);
				return NULL;
			}
		}
		else if(strcmp(cmd,"kernel_type")==0)
		{		
			fscanf(fp,"%80s",cmd);
			int i;
			for(i=0;kernel_type_table[i];i++)
			{
				if(strcmp(kernel_type_table[i],cmd)==0)
				{
					param.kernel_type=i;
					break;
				}
			}
			if(kernel_type_table[i] == NULL)
			{
				//fprintf(stderr,"unknown kernel function.\n");
				AfxMessageBox("未知的未知的核函数\n");
				free(model->rho);
				free(model->label);
				free(model->nSV);
				free(model);
				return NULL;
			}
		}
		else if(strcmp(cmd,"degree")==0)
			fscanf(fp,"%d",&param.degree);
		else if(strcmp(cmd,"gamma")==0)
			fscanf(fp,"%lf",&param.gamma);
		else if(strcmp(cmd,"coef0")==0)
			fscanf(fp,"%lf",&param.coef0);
		else if(strcmp(cmd,"nr_class")==0)
			fscanf(fp,"%d",&model->nr_class);
		else if(strcmp(cmd,"total_sv")==0)
			fscanf(fp,"%d",&model->l);
		else if(strcmp(cmd,"rho")==0)
		{
			int n = model->nr_class * (model->nr_class-1)/2;
			model->rho = Malloc(double,n);
			for(int i=0;i<n;i++)
				fscanf(fp,"%lf",&model->rho[i]);
		}
		else if(strcmp(cmd,"label")==0)
		{
			int n = model->nr_class;
			model->label = Malloc(int,n);
			for(int i=0;i<n;i++)
				fscanf(fp,"%d",&model->label[i]);
		}
		else if(strcmp(cmd,"probA")==0)
		{
			int n = model->nr_class * (model->nr_class-1)/2;
			model->probA = Malloc(double,n);
			for(int i=0;i<n;i++)
				fscanf(fp,"%lf",&model->probA[i]);
		}
		else if(strcmp(cmd,"probB")==0)
		{
			int n = model->nr_class * (model->nr_class-1)/2;
			model->probB = Malloc(double,n);
			for(int i=0;i<n;i++)
				fscanf(fp,"%lf",&model->probB[i]);
		}
		else if(strcmp(cmd,"nr_sv")==0)
		{
			int n = model->nr_class;
			model->nSV = Malloc(int,n);
			for(int i=0;i<n;i++)
				fscanf(fp,"%d",&model->nSV[i]);
		}
		else if(strcmp(cmd,"SV")==0)
		{
			while(1)
			{
				int c = getc(fp);
				if(c==EOF || c=='\n') break;	
			}
			break;
		}
		//else if(strcmp(cmd,"#INF")==0)//添加//什么也不作，往后移动
		//{
		//	while(1)
		//	{
		//		int c = getc(fp);
		//		if(c==EOF || c=='\n') break;	
		//	}
		//	break;
		//}

		else
		{
			//fprintf(stderr,"unknown text in model file\n");
			AfxMessageBox("无法识别的字符在model文件\n");
			free(model->rho);
			free(model->label);
			free(model->nSV);
			free(model);
			return NULL;
		}
	}

	// read sv_coef and SV

	int elements = 0;
	long pos = ftell(fp);

	while(1)
	{
		int c = fgetc(fp);
		switch(c)
		{
			case '\n':
				// count the '-1' element
			case ':':
				++elements;
				break;
			case EOF:
				goto out;
			default:
				;
		}
	}
out:
	fseek(fp,pos,SEEK_SET);

	int m = model->nr_class - 1;
	int l = model->l;
	model->sv_coef = Malloc(double *,m);
	int i;
	for(i=0;i<m;i++)
		model->sv_coef[i] = Malloc(double,l);
	model->SV = Malloc(svm_node*,l);
	svm_node *x_space=NULL;
	if(l>0) x_space = Malloc(svm_node,elements);

	int j=0;
	for(i=0;i<l;i++)
	{
		model->SV[i] = &x_space[j];
		for(int k=0;k<m;k++)
			fscanf(fp,"%lf",&model->sv_coef[k][i]);
		while(1)
		{
			int c;
			do {
				c = getc(fp);
				//if(c=='#'|| c=='I'|| c=='N' || c=='F') 
				//	c='0';
				if(c=='\n') 
					goto out2;
			} while(isspace(c));
			ungetc(c,fp);
			fscanf(fp,"%d:%lf",&(x_space[j].index),&(x_space[j].value));
			++j;
		}	
out2:
		x_space[j++].index = -1;
	}

	if (ferror(fp) != 0 || fclose(fp) != 0) return NULL;

	model->free_sv = 1;	// XXX
	return model;
}

void svm_destroy_model(svm_model* model)
{
	if(model->free_sv && model->l > 0)
		free((void *)(model->SV[0]));
	for(int i=0;i<model->nr_class-1;i++)
		free(model->sv_coef[i]);
	free(model->SV);
	free(model->sv_coef);
	free(model->rho);
	free(model->label);
	free(model->probA);
	free(model->probB);
	free(model->nSV);
	free(model);
}

void svm_destroy_param(svm_parameter* param)
{
	free(param->weight_label);
	free(param->weight);
}

double 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 "未知的svm类型";//unknown svm type
	
	// kernel_type, degree
	
	int kernel_type = param->kernel_type;
	if(kernel_type != LINEAR &&
	   kernel_type != POLY &&
	   kernel_type != RBF &&
	   kernel_type != SIGMOID &&
	   kernel_type != PRECOMPUTED)
		return "未知的核";

	if(param->degree < 0)
		return "多项式核的度不能小于0";//degree of polynomial kernel < 0

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

	if(param->cache_size <= 0)
		return "缓存要大于0M";//cache_size <= 0

	if(param->eps <= 0)
		return "运算的终止判据要大于0";//eps <= 0

	if(svm_type == C_SVC ||
	   svm_type == EPSILON_SVR ||
	   svm_type == NU_SVR)
		if(param->C <= 0)
			return "惩罚因子要大于0";//C <= 0

	if(svm_type == NU_SVC ||
	   svm_type == ONE_CLASS ||
	   svm_type == NU_SVR)
		if(param->nu <= 0 || param->nu > 1)
			return "ν-SVC，一类SVM和ν- SVR的参数必须在0～1之间";//nu <= 0 or nu > 1

	if(svm_type == EPSILON_SVR)
		if(param->p < 0)
			return "ε-SVR 中损失函数ε的值要大于0";//p < 0

	if(param->shrinking != 0 &&
	   param->shrinking != 1)
		return "使用启发式，只能取0或1";//shrinking != 0 and shrinking != 1

	if(param->probability != 0 &&
	   param->probability != 1)
		return "是否需要预报概率，只能取0或1";//probability != 0 and probability != 1

	if(param->probability == 1 &&
	   svm_type == ONE_CLASS)
		return "一类SVM不支持概率预测";//one-class SVM probability output not supported yet

*/