svr.cpp

基于线性规划的回归支持向量机源程序

// SVC.cpp: implementation of the CSVR class.
//
//////////////////////////////////////////////////////////////////////

#include "SVR.h"
#include "matlib.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CSVR::CSVR()
{
	m_type=LINEAR;
}

CSVR::~CSVR()
{

}
double CSVR::kernel(vector& x1,vector& x2)
{
	double data=0;
	int m=x1();
	int i;
	switch(m_type)
	{
		case LINEAR:
			data=dot(x1,x2);
			break;
		case GAUSS:
			for(i=0;i<m;i++)
			{
				double d=x1[i]-x2[i];
				data+=d*d;
			}
			data/=(m_sigma*m_sigma);
			data*=-0.5;
			data=exp(data);
			break;
	}
	return data;
}
double CSVR::SolveByNorm8(matrix& x,vector& y, double epslon, double C)
{
	int N=y();
	int n=mcol(x);
	matrix A;
	matrix A1;
	vector b;
	vector c;
	vector c1;
	IntVector js;
	int i;
	double werr=0;
	if(m_type==LINEAR)
	{
		A.initial(N+n+2,3*N+2*n+1);
		A1.initial(N+n+2,3*N+2*n+1+N+n+2);
		A=0;
		A1=0;
		b.initial(N+n+2);
		c.initial(3*N+2*n+1);
		c1.initial(3*N+2*n+1+N+n+2);
		b=0;c=0;c1=0;
		js.initial(N+n+2);
		for(i=0;i<n;i++)
		{
			for(int j=0;j<N;j++)
			{
				double data=x(j,i);
				A(i+1,j)=data;
				A(i+1,j+N)=-data;
				A1(i+1,j)=data;
				A1(i+1,j+N)=-data;
			}
			A(0,2*N+i)=1;
			A(0,2*N+n+i)=1;
			A(i+1,2*N+i)=1;
			A(i+1,2*N+n+i)=-1;
			A1(0,2*N+i)=1;
			A1(0,2*N+n+i)=1;
			A1(i+1,2*N+i)=1;
			A1(i+1,2*N+n+i)=-1;

		}
		for(i=0;i<N;i++)
		{
			A(n+1,i)=1;
			A(n+1,N+i)=-1;
			A(n+2+i,i)=1;
			A(n+2+i,i+N)=1;
			A(n+2+i,i+2*N+2*n+1)=1;
			A1(n+1,i)=1;
			A1(n+1,N+i)=-1;
			A1(n+2+i,i)=1;
			A1(n+2+i,i+N)=1;
			A1(n+2+i,i+2*N+2*n+1)=1;
			b[n+2+i]=C;
			c[i]=epslon-y[i];
			c[N+i]=epslon+y[i];
		}
		A(0,2*N+2*n)=1;
		A1(0,2*N+2*n)=1;
		b[0]=1;
		
		for(i=0;i<js();i++)
		{
			A1(i,3*N+2*n+1+i)=1;
			c1[3*N+2*n+1+i]=1;
			js[i]=3*N+2*n+1+i;
		}
	ofstream ofile("debug.txt");
		vector v;
ofile<<"A1=\n"<<A1<<"b="<<b<<"\n c1="<<c1<<"\n js="<<js<<"\n dim of js="<<js()<<endl;
		if(lp_simplex(A1,b,c1,js,v)<=0)
		{
			cout<<"no soluntion for feasible solution"<<endl;
		}
		js.sort();
ofile<<"v="<<v<<"\n js="<<js<<endl;
		int mm=3*N+2*n+1;
		for(i=0;i<js();i++)
			if(js[i]>=mm)
			{
				cout<<"No limite solution"<<endl;
				return -1;
			}
		vector u;
		if(lp_simplex(A,b,c,js,u)<=0)
			{
				cout<<"no soluntion for optimal"<<endl;
				return -1;
			}
cout<<"u=\n"<<u<<endl;
		c*=-1;
		int m=0;
		mm=2*N+2*n;
		js.sort();
ofile<<"after compute dual\n A=\n"<<A<<"b="<<b<<"\n c="<<c<<"\n js="<<js<<"\n dim of js="<<js()<<endl;
		for(i=0;i<js();i++)
		{
			if(js[i]<mm)
			{
				m++;
			}
		}
ofile<<"2N+2n="<<mm<<", num of eqation="<<m<<endl;
		IntVector is(m),ks(N+n+2);
		matrix a(m,m);
		vector c1(m);
		v.initial(N+n+2);
		v=0;
		c1=0;
		ks=1;
		int s=0;
		for(i=0;i<js();i++)
		{
			if(js[i]<mm)
			{
				is[s]=js[i];
				s++;
			}
			else
			{
				if(js[i]==mm)
					ks[0]=0;
				else
					ks[js[i]-mm-1+n+2]=0;
			}
		}
		s=0;
ofile<<"is="<<is<<"\n ks="<<ks<<endl;
		for(i=0;i<ks();i++)
		{
			if(ks[i])
			{
				for(int j=0;j<m;j++)
				{
					a(j,s)=A(i,is[j]);
				}
				s++;
			}
		}
		for(i=0;i<m;i++)
			c1[i]=c[is[i]];
ofile<<"a=\n"<<a<<"b="<<c1<<endl;
matrix e;
e=pinv(a);
A=e*a*e;

ofile<<"pinv(a)=\n"<<e;
dcinv(a);
ofile<<"dcinv(a)=\n"<<e;
		u=a*c1;
ofile<<"a^-1*b="<<u<<endl;
		s=0;
		for(i=0;i<ks();i++)
		{
			if(ks[i])
			{
				v[i]=u[s];
				s++;
			}
		}
		m_b=v[n+1];
		m_alpha.initial(n);
		for(i=0;i<n;i++)
			m_alpha[i]=v[i+1];
		werr=0;
		for(i=n+2;i<v();i++)
			werr+=v[i];
		ofile<<"v="<<v<<endl;
		ofile<<"m_alpha="<<m_alpha<<endl;
		werr/=N;
	}
	else
	{
		A.initial(2*N+2,5*N+1);
		A1.initial(2*N+2,5*N+1+2*N+2);
		A=0;
		A1=0;
		b.initial(2*N+2);
		c.initial(5*N+1);
		c1.initial(5*N+1+2*N+2);
		b=0;c=0;c1=0;
		js.initial(2*N+2);
		for(i=0;i<N;i++)
		{
			for(int j=0;j<N;j++)
			{
				double d=1;
				if(j!=i) d=kernel(x[i],x[j]);
				A(i+1,j)=d;
				A(i+1,j+N)=-d;
			}
			A(0,2*N+i)=1;
			A(0,3*N+i)=1;
			A(i+1,2*N+i)=1;
			A(i+1,3*N+i)=-1;
			A(N+1,i)=1;
			A(N+1,N+i)=-1;
			A(N+2+i,i)=1;
			A(N+2+i,i+N)=1;
			A(N+2+i,i+4*N+1)=1;

			A1(0,2*N+i)=1;
			A1(0,3*N+i)=1;
			A1(i+1,2*N+i)=1;
			A1(i+1,3*N+i)=-1;
			A1(N+1,i)=1;
			A1(N+1,N+i)=-1;
			A1(N+2+i,i)=1;
			A1(N+2+i,i+N)=1;
			A1(N+2+i,i+4*N+1)=1;

			c[i]=epslon-y[i];
			c[N+i]=epslon+y[i];
			b[N+2+i]=C;
		}
		A(0,4*N)=1;
		A1(0,4*N)=1;
		b[0]=1;
		
		for(i=0;i<js();i++)
		{
			A1(i,5*N+1+i)=1;
			c1[5*N+1+i]=1;
			js[i]=5*N+1+i;
		}
		vector v;
		if(lp_simplex(A1,b,c1,js,v)<=0)
		{
			cout<<"no soluntion for feasible solution"<<endl;
		}
		js.sort();
		
		for(i=0;i<js();i++)
			if(js[i]>=4*N+1)
			{
				cout<<"No limite solution"<<endl;
				return 0;
			}
		vector u;
	ofstream ofile("debug.txt");
	ofile<<"A=\n"<<A<<"b="<<b<<"\n c="<<c<<"\n js="<<js<<"\n dim of js="<<js()<<endl;
		if(lp_simplex(A,b,c,js,u)<=0)
			{
				cout<<"no soluntion for optimal"<<endl;
			}
//cout<<"u=\n"<<u<<endl;
		c*=-1;
		int m=0;
		n=4*N;
		js.sort();
ofile<<"after compute dual\n A=\n"<<A<<"b="<<b<<"\n c="<<c<<"\n js="<<js<<"\n dim of js="<<js()<<endl;
		for(i=0;i<js();i++)
		{
			if(js[i]<n)
			{
				m++;
			}
		}
ofile<<"4*N="<<n<<", num of eqation="<<m<<endl;
		IntVector is(m),ks(2*N+2);
		matrix a(m,m);
		vector c1(m);
		v.initial(2*N+2);
		v=0;
		c1=0;
		ks=1;
		int s=0;
		for(i=0;i<js();i++)
		{
			if(js[i]<n)
			{
				is[s]=js[i];
				s++;
			}
			else
			{
				if(js[i]==n)
					ks[0]=0;
				else
					ks[js[i]-n-1+N+2]=0;
			}
		}
		s=0;
ofile<<"is="<<is<<"\n ks="<<ks<<endl;
		for(i=0;i<ks();i++)
		{
			if(ks[i])
			{
				for(int j=0;j<m;j++)
				{
					a(j,s)=A(i,is[j]);
				}
				s++;
			}
		}
		for(i=0;i<m;i++)
			c1[i]=c[is[i]];
ofile<<"a=\n"<<a<<"b="<<c1<<endl;
		dcinv(a);
		u=a*c1;
ofile<<"a^-1*b="<<u<<endl;
		s=0;
		for(i=0;i<ks();i++)
		{
			if(ks[i])
			{
				v[i]=u[s];
				s++;
			}
		}
ofile<<"v="<<v<<endl;
		m_b=v[N+1];
		int k=0;
		for(i=1;i<=N;i++)
		{
			if(ks[i])
				k++;
		}
		js.initial(k);
		m_x.initial(k,n);
		m_y.initial(k);
		m_alpha.initial(k);
		k=0;
		for(i=0;i<N;i++)
		{
			if(ks[i+1])
			{
				js[k]=i+1;
				m_alpha[k]=v[i+1];
				m_y[k]=y[i];
				for(int j=0;j<mcol(x);j++)
					m_x[k][j]=x[i][j];
				k++;
			}
		}
		werr=0;
		for(i=N+2;i<v();i++)
			werr+=v[i];
		werr/=N;
	}
//ofstream ofile("debug.txt");
//ofile<<"A=\n"<<A<<"b="<<b<<"c="<<c<<endl;
	return werr;
}

int CSVR::SolveByNorm1(matrix& x,vector& y, double C)
{
	int N=y();
	int n=mcol(x);
	matrix A;
	vector b;
	vector c;
	IntVector js;
	int i;
	if(m_type==LINEAR)
	{
		A.initial(N,2*n+2*N+2);
		A=0;
		b.initial(N);
		c.initial(2*n+2*N+2);
		b=1;c=0;
		js.initial(N);
		for(i=0;i<N;i++)
		{
			for(int j=0;j<n;j++)
			{
				double d=y[i]*x[i][j];
				A(i,j)=d;
				A(i,n+j)=-d;
			}
			A(i,2*n)=y[i];
			A(i,2*n+1)=-y[i];
			A(i,2*n+2+i)=1;
			A(i,2*n+N+2+i)=-1;
			c[2*n+2+i]=C;
			js[i]=2*n+2+i;
		}
		for(i=0;i<n;i++)
		{
			c[i]=1;
			c[n+i]=1;
		}
	}
	else
	{
		A.initial(N,4*N+2);
		b.initial(N);
		c.initial(4*N+2);
		js.initial(N);
		A=0;b=1;c=0;
		for(i=0;i<N;i++)
		{
			for(int j=0;j<N;j++)
			{
				double d=y[i];
				if(j!=i) d*=kernel(x[i],x[j]);
				A(i,j)=d;
				A(i,N+j)=-d;
			}
			A(i,2*N)=y[i];
			A(i,2*N+1)=-y[i];
			A(i,2*N+2+i)=1;
			A(i,3*N+2+i)=-1;
			c[i]=1;
			c[N+i]=1;
			c[2*N+2+i]=C;
			js[i]=2*N+2+i;
		}
	}
ofstream ofile("debug.txt");
ofile<<"A="<<endl;
ofile<<A<<"b="<<b;
	vector u;
	lp_simplex(A,b,c,js,u);
//	lp_affineinner(A,b,c,u);
//cout<<"u=\n"<<u<<endl;
	int k=0;
	switch(m_type)
	{
	case LINEAR:
		m_b=u[2*n]-u[2*n+1];
		m_alpha.initial(n);
		for(i=0;i<n;i++)
		{
			if((u[i]>0.0000000001)&&(u[n+i]>0.0000000001))
			{
				cout<<i<<"'s, no proper:"<<u[i]<<"-"<<u[n+i]<<endl;
			}
			m_alpha[i]=u[i]-u[n+i];
		}
		break;
	case GAUSS:
		m_b=u[2*N]-u[2*N+1];
		b.initial(N);
		k=N;
		for(i=0;i<N;i++)
		{
			if((u[i]>0.0000000001)&&(u[N+i]>0.0000000001))
			{
				cout<<i<<"'s, no proper:"<<u[i]<<"-"<<u[N+i]<<endl;
			}
			b[i]=u[i]-u[N+i];
			double data=fabs(b[i]);
			if(data<0.0000000001)
				k--;
		}
		js.initial(k);
		m_x.initial(k,n);
		m_y.initial(k);
		m_alpha.initial(k);
		k=0;
		for(i=0;i<N;i++)
		{
			double data=fabs(b[i]);
			if(data>=0.0000000001)
			{
				js[k]=i;
				m_alpha[k]=b[i];
				m_y[k]=y[i];
				for(int j=0;j<mcol(x);j++)
					m_x[k][j]=x[i][j];
				k++;
			}
		}
		break;
	}
	return 1;
}


int CSVR::SolveByNorm2(matrix& x,vector& y, double C)
{
	return 1;
}
double CSVR::ERR(matrix& x, vector& y)
{
	double derr=0;
	int N=y();
	int i;
	if(m_type==LINEAR)
	{
		vector v;
		v=x*m_alpha;
		v+=m_b;
		for(int i=0;i<v();i++)
		{
			double data=fabs(v[i]-y[i]);
			derr+=data;
		}
	}
	else
	{
		for(i=0;i<N;i++)
		{
			double data=simulate(x[i])-y[i];
			derr+=fabs(data);
		}
	}
	derr/=N;
	return derr;
}
double CSVR::simulate(vector&v)
{
	double data=0;
	int m=mrow(m_x);
	int i;
	if(m_type==LINEAR)
	{
		data=dot(v,m_alpha);
		data+=m_b;
	}
	else
	{
		for(i=0;i<m;i++)
		{
			double d;
			d=m_alpha[i]*kernel(v,m_x[i]);
			data+=d;
		}
		data+=m_b;
	}
	return data;
}