smo.cpp

支持向量机的SMO算法实现数据拟合训练的基本程序

#include <iostream.h>
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <stdlib.h>


#define MAX(x,y) (x>y?x:y)
#define MIN(x,y) (x<y?x:y)
#define MOD(x) (x-2*(int)(x/2)) 
#define ABS(x) (x>0?x:(-x))

#define N 6
#define K 50

double x[50][6],target[50];
double *alpha,*alpha_star;
double epsilon;
int k;
double C;
double b;

double kernel(double *x,double *y)
{
	int i,n;
	double result;

	n=N;

	result=0;
	for(i=0;i<n;i++)
		result+=((*x++)*(*y++));
	result=(result+1)*(result+1);
//	result=(result+1);
	return result;
}

int numalpha()
{
	int i;
	int k;
	int number;
	k=K;
	number=0;
	for(i=0;i<k;i++)
	{
		if(alpha[i]!=0&&alpha[i]!=C)
			number++;
	}
	return number;
}

double calculate_b(double *p1,double *p2,int l,int m)
{
	double b12,b1,b2;
	double y1,y2;
	int i,i1,i2;
	double *p;

	double kernel(double *x,double *y);

	p=NULL;
	
	i1=l;
	i2=m;

	y1=target[i1];
	y2=target[i2];
	b1=0;
	b2=0;

	for(i=0;i<K;i++)
	{
		p=&x[i][0];//p=x+N*i;
		b1+=((alpha[i]-alpha_star[i])*kernel(p,p1));
		b2+=((alpha[i]-alpha_star[i])*kernel(p,p2));
	}
	b1=y1-b1;
	b2=y2-b2;
	b12=(b1+b2)/2;
	return b12;
}
double svm(double *p1)
{
	int i;
	int n,k;
	double *p;
	double result;

	double kernel(double *x,double *y);
	
	n=N;
	k=K;
	p=NULL;
	p=p1;
	
	result=b;
	for(i=0;i<k;i++)
		result+=(alpha[i]-alpha_star[i])*kernel(p,&x[i][0]);
	return result;
}


int takestep(int l,int m)
{
	int i1,i2,n;
	double alpha1,alpha1_star,alpha2,alpha2_star,y1,y2;
	double alpha1old,alpha1old_star,alpha2old,alpha2old_star;

	double L,H;

	double phi1,phi2;
	double delta_phi;

	double k11,k12,k22;
	double eta;
	double gamma;
	double object1,object2;

	double a1,a2;
	double *p,*p1,*p2;

	int case1,case2,case3,case4,finished;
	n=N;

	p=NULL;
	p1=NULL;
	p2=NULL;

	double svm(double *x);
	double kernel(double *x,double *y);
	double calculate_b(double *p1,double *p2,int l,int m);


	i1=l;
	i2=m;

	cout<<i1<<"   "<<i2<<"\n";
	if (i1==i2) return 0;

	alpha1=alpha[i1];
	alpha1_star=alpha_star[i1];
	alpha2=alpha[i2];
	alpha2_star=alpha_star[i2];
	
	y1=target[i1];
	y2=target[i2];
	
	p=&x[i1][0];
	phi1=svm(p)-y1;
	cout<<"phi1="<<phi1<<"\n";
	
	p=&x[i2][0];
	phi2=svm(p)-y2;

	p1=&x[i1][0];
	p2=&x[i2][0];

	k11=kernel(p1,p1);
	k12=kernel(p1,p2);
	k22=kernel(p2,p2);

	eta=k11+k22+2*k12;
	gamma=alpha1-alpha1_star+alpha2-alpha2_star;

	cout<<"eta="<<eta<<"\n";
	cout<<"gamma="<<gamma<<"\n";

	alpha1old=alpha1;alpha1old_star=alpha1_star;
	alpha2old=alpha2;alpha2old_star=alpha2_star;

	delta_phi=phi1-phi2;
	cout<<"delta_phi="<<delta_phi<<"\n";

	case1=0;case2=0;case3=0;case4=0;finished=0;

	while(!finished)
	{
		bool aa;
		aa=((case1==0)&& 
			(alpha1>0||(alpha1_star==0 && delta_phi>0))&&
			(alpha2>0||(alpha2_star==0 &&delta_phi<0)));
		cout<<"aa="<<aa<<"\n";
		if(aa)
		{
//			cout<<"//compute L,H";

			L=MAX(0,gamma-C);
			H=MIN(C,gamma);
			cout<<"L="<<L<<"\n"<<"H="<<H<<"\n";
			if(L<H)
			{
				if(eta>0)
				{
					a2=alpha2-delta_phi/eta;
					a2=MIN(a2,H);
					a2=MAX(L,a2);
					a1=alpha1-(a2-alpha2);
				}
				else
				{
					a2=L;
					a1=alpha1-(a2-alpha2);
					object1=-0.5*a1*a1*eta+a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
					a2=H;
					a1=alpha1-(a2-alpha2);
					object2=-0.5*a1*a1*eta+a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
					if(object1>object2) a2=L;
					else a2=H;
					a1=alpha1-(a2-alpha2);
				}
				if(ABS((delta_phi-eta*(a1-alpha1)))>epsilon)
				{
					alpha1=a1;
					alpha2=a2;
				}
			
		}
			else
				finished=1;
			case1=1;
		}
		else
		{
			aa=((case2==0)&&
				(alpha1>0||(alpha1_star==0&&delta_phi>2*epsilon))&&
				(alpha2_star>0||(alpha2==0&&delta_phi<2*epsilon)));
			cout<<"aa2="<<aa<<"\n";
			if(aa)//((case2==0)&&(alpha1>0||(alpha1_star==0&&delta_phi>2*epsilon))&&(alpha2_star>0||(alpha2==0&&delta_phi<2*epsilon)))
			{
				//compute L,H;
				L=MAX(0,-gamma);
				H=MIN(C,C-gamma);
				//L=MAX(0,gamma);
				//H=MIN(C,C+gamma);
				cout<<"L="<<L<<"\n"<<"H="<<H<<"\n";

				if(L<H)
				{
					if(eta>0)
					{
						a2=alpha2_star+(delta_phi-2*epsilon)/eta;
						a2=MIN(a2,H);
						a2=MAX(L,a2);
						a1=alpha1+(a2-alpha2_star);
				
					}
					else
					{
						a2=L;
						a1=alpha1+(a2-alpha2_star);
						object1=-0.5*a1*a1*eta-2*epsilon*a1+a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
						a2=H;
						a1=alpha1+(a2-alpha2);
						object1=-0.5*a1*a1*eta-2*epsilon*a1+a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
						if(object1>object2) a2=L;
						else a2=H;
						a1=alpha1+(a2-alpha2_star);
					}
					if(ABS((delta_phi-eta*(a1-alpha1)))>epsilon)
					{
						cout<<"alpha1_star="<<a1<<"\n"<<"alpha2_star="<<a2<<"\n";
						alpha1=a1;
						alpha2_star=a2;
					}
					
				}
				else
					finished=1;
				case2=1;
			}
			else
		{
			aa=((case3==0)&&
			(alpha1_star>0||(alpha1==0&&delta_phi<2*epsilon))&&
			(alpha2>0||(alpha2_star==0&&delta_phi<2*epsilon)));
			cout<<"aa3="<<aa<<"\n";
			if(aa)//((case3==0)&&(alpha1_star>0||(alpha1==0&&delta_phi<2*epsilon))&&(alpha2>0||(alpha2_star==0&&delta_phi<2*epsilon)))
			{
				//compute L,H;
				L=MAX(0,gamma);
				H=MIN(C,C+gamma);
			//	L=MAX(0,-gamma);
			//	H=MIN(C,C-gamma);
				cout<<"L="<<L<<"\n"<<"H="<<H<<"\n";
				if(L<H)
				{
					if(eta>0)
					{
						a2=alpha2-(delta_phi+2*epsilon)/eta;//a2=alpha2-(delta_phi-2*epsilon)/eta;
						a2=MIN(a2,H);
						a2=MAX(L,a2);
						a1=alpha1_star+(a2-alpha2);
					}
					else
					{
						a2=L;
						a1=alpha1_star+(a2-alpha2);
						object1=-0.5*a1*a1*eta-2*epsilon*a1-a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
						a2=H;
						a1=alpha1_star+(a2-alpha2);
						object2=-0.5*a1*a1*eta-2*epsilon*a1-a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
						if(object1>object2) a2=L;
						else a2=H;
						a1=alpha1_star+(a2-alpha2);
					}
					if(ABS((delta_phi-eta*(-a1+alpha1_star)))>epsilon)
					{
						cout<<"alpha1_star="<<a1<<"\n"<<"alpha2="<<a2<<"\n";
						alpha1_star=a1;
						alpha2=a2;
					}


				}
				else
					finished=1;
				case3=1;
			}
		else
		{
			aa=((case4==0)&&
			(alpha1_star>0||(alpha1==0&&delta_phi<0))&&
			(alpha2_star>0||(alpha2==0&&delta_phi>0)));
			cout<<"aa4="<<aa<<"\n";
			if(aa)//((case4==0)&&(alpha1_star>0||(alpha1==0&&delta_phi<0))&&(alpha2_star>0||(alpha2==0&&delta_phi>0)))
			{
				//compute L,H;
				L=MAX(0,-gamma-C);
				H=MIN(C,-gamma);
				cout<<"L="<<L<<"\n"<<"H="<<H<<"\n";
				if(L<H)
				{
					if(eta>0)
					{
						a2=alpha2_star+delta_phi/eta;
						a2=MIN(a2,H);
						a2=MAX(L,a2);
						a1=alpha1_star-(a2-alpha2_star);
				
					}
					else
					{
						a2=L;
						a1=alpha1_star-(a2-alpha2_star);
						object1=-0.5*a1*a1*eta-a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
						a2=H;
						a1=alpha1_star-(a2-alpha2_star);
						object1=-0.5*a1*a1*eta-a1*(delta_phi+(alpha1old-alpha1old_star)*eta);
						if(object1>object2) a2=L;
						else a2=H;
						a1=alpha1_star-(a2-alpha2_star);
					}
					if(ABS((delta_phi-eta*(-a1+alpha1_star)))>epsilon)
					{
						cout<<"alpha1_star="<<a1<<"\n"<<"alpha2_star="<<a2<<"\n";
						alpha1_star=a1;
						alpha2_star=a2;
					}
				}
				else
					finished=1;
				case4=1;
			}
		else
			finished=1;
		}
		}
		}
		delta_phi=(delta_phi-eta*((alpha1-alpha1_star)-(alpha1old-alpha1old_star)));
	}
	
	alpha[i1]=alpha1;
	alpha_star[i1]=alpha1_star;
	alpha[i2]=alpha2;
	alpha_star[i2]=alpha2_star;

	for(int i=0;i<K;i++)
	{
		cout<<"alpha["<<i<<"]="<<alpha[i]<<"	"<<"alpha_star["<<i<<"]="<<alpha_star[i]<<"\n";
	}
//	update thresold

	b=calculate_b(p1,p2,i1,i2);
	cout<<"b="<<b<<"\n";

	if(ABS(delta_phi)>epsilon)
		return 1;
	else
		return 0;
}


int examineexample(int m)
{
	int i;
	int i1,i2;
	double y1,y2;
	double alpha2,alpha2_star;
	double phi1,phi2,delta_phi;

	double *p;
	int n;
	n=N;
	p=NULL;
	
	double svm(double *p1);
	int takestep (int i1,int i2);
	int numalpha();

	i2=m;
	y2=target[i2];
	alpha2=alpha[i2];
	alpha2_star=alpha_star[i2];

	p=&x[i2][0];
	phi2=svm(p)-y2;
	cout<<"phi2="<<phi2<<"\n";

	bool bb;
	bb=((phi2>epsilon&&alpha2_star<C)||
		(phi2<epsilon&&alpha2_star>0)||
		(-phi2>epsilon&&alpha2<C)||
		(-phi2>epsilon&&alpha2>0));
	if(bb)
	{
		if(numalpha()>1)
		{
			delta_phi=0;
			i1=0;
			for(i=0;i<K;i++)
			{
				p=&x[i][0];
				y1=target[i];
				phi1=svm(p)-y1;
				if(ABS(phi1-phi2)>delta_phi) 
				{
					delta_phi=ABS(phi1-phi2) ;
					i1=i;
					cout<<"according to the max(deltaphi):i1="<<i1<<"\n";
				}
			}
			if (takestep(i1,i2)) return 1;
		}

		for(i=0;i<K;i++)
		{
			i1=rand()%K;
			if (alpha[i1]!=0 &&alpha[i1]!=C)
			{	
				if (takestep(i1,i2)) return 1;
			}
		}
		for(i=0;i<K;i++)
		{
			i1=rand()%K;
			if (takestep(i1,i2)) return 1;
		}

	}
	return 0;
}
	
void smomain()
{
	int i,j;
	int numsamples;
	double sigfig;//epsilon,
	int numchanged,examineall,loopcounter,minimumnumchanged;

//	double C;
	double pobject,dobject;
	double *p1,*p2;
	
	int examineexample(int i);

	p1=NULL;
	p2=NULL;

	k=K;
	numsamples=k;

	C=0.009;
	b=0;

	alpha=(double *)malloc(k*sizeof(double));
	alpha_star=(double *)malloc(k*sizeof(double));

	for(i=0;i<k;i++)
	{
		alpha[i]=0;
		alpha_star[i]=0;
	}
	epsilon=0.001;
	numchanged=0;
	examineall=1;
	sigfig=-100;
	loopcounter=0;

	while((((numchanged>0||examineall))||(sigfig<=-3))&&(loopcounter<10))
	{
		cout<<"loopcounter="<<loopcounter<<"\n";
		loopcounter++;
		numchanged=0;

		if(examineall)
			for(i=0;i<k;i++)
			{
				numchanged+=examineexample(i);
				cout<<"numchanged="<<numchanged<<"\n";
			}
		else
			for(i=0;i<k;i++)
			{
				if(alpha[i]!=0 && alpha[i]!=C)
					numchanged+=examineexample(i);
				cout<<"numchanged="<<numchanged<<"\n";
			}

		if(MOD(loopcounter)==0)
			minimumnumchanged=(int)MAX(1,0.1*numsamples);
		else
			minimumnumchanged=1;
		
		if(examineall==1)
			examineall=0;
		else
			if(numchanged<minimumnumchanged)
				examineall=1;
	
		dobject=0;
		pobject=0;
		for(i=0;i<K;i++)
		{
			p1=&x[i][0];//p1=x+N*i;
//			for(j=0;j<K;j++)
//			{
//				p2=&x[j][0];//p2=x+N*j;
//				pobject+=(-0.5*(alpha[i]-alpha_star[i])*(alpha[j]-alpha_star[j])*kernel(p1,p2));
//			}
//			pobject+=(target[i]*(alpha[i]-alpha_star[i])-epsilon*(alpha[i]+alpha_star[i]));
			dobject+=(MAX(0,svm(p1)-target[i]-epsilon)*(C-alpha_star[i]));
			dobject-=(MIN(0,svm(p1)-target[i]-epsilon)*alpha_star[i]);
			dobject+=(MIN(0,target[i]-epsilon-svm(p1))*(C-alpha[i]));
			dobject-=(MAX(0,target[i]-epsilon-svm(p1))*alpha[i]);
		}

		for(i=0;i<K;i++)
		{
			p1=&x[i][0];
			pobject+=(0.5*(alpha[i]-alpha_star[i])*(svm(p1)-b));
			pobject-=(epsilon*(alpha[i]+alpha_star[i]));
			pobject+=(target[i]*(alpha[i]-alpha_star[i]));
		}
		pobject+=dobject;
		
		cout<<"dobject="<<dobject<<"pobject="<<ABS(pobject)<<"\n";

		sigfig=log10(dobject/(ABS(pobject)+1));
		cout<<"SigFig="<<sigfig<<"\n";

		FILE *fp;
		fp=fopen("d:\\zhb\\regression\\sigfig.txt","a");
		fprintf(fp,"%f\n",sigfig);
		fclose(fp);
//		char qq;
//		cin>>qq;
	
	}
}

void main()
{
	int i,j,q;//,n,k;
	double x1[50][6],u[50];
	double x2[6];
	double *p;
	int n,k;
	p=NULL;
	
	FILE *filename;
	char temp;
	float value;
	filename=fopen("d:\\zhb\\regression\\data\\table5_9.txt","r");
//	fscanf(filename,"d% %d",&N,&K);
//	temp=fgetc(filename);
	for(i=0;i<K;i++)
	{
		for(j=0;j<N;j++)
		{
			fscanf(filename,"%f",&value);
			x1[i][j]=value;
			temp=fgetc(filename);
		}
		fscanf(filename,"%f",&value);
		u[i]=value;
	}
	fclose(filename);

	
	for(i=0;i<K;i++)
	{
		target[i]=u[i];
		cout<<target[i]<<"\n";
		for(j=0;j<N;j++)
		{
			x[i][j]=x1[i][j];
		}
	}

	smomain();
	
	for(i=0;i<K;i++)
		cout<<"alpha["<<i<<"]="<<alpha[i]<<"   "<<"alpha_star["<<i<<"]="<<alpha_star[i]<<"\n";
	
	filename=fopen("d:\\zhb\\regression\\data\\check5_9.txt","r");
	fscanf(filename,"%d",&q);
	temp=fgetc(filename);
	for(i=0;i<q;i++)
	{
		for(j=0;j<N;j++)
		{
			fscanf(filename,"%f",&value);
			x2[j]=value;
			temp=fgetc(filename);
		}

		p=&x2[0];
		double result;
		result=svm(p);
		FILE *fp1;
		fp1=fopen("d:\\zhb\\regression\\result1.txt","a");
		fprintf(fp1,"%f\n",result);
		fclose(fp1);
		cout<<"result="<<result<<"\n";
	}
	fclose(filename);
	free(alpha);
	free(alpha_star);
}