covernn.cpp

thank you for using this software,which is fugai method,when you can t catch it,please ask me for my

	while(k <= Dim_in && k<allpt)
	{
		GetProjection(pt_proj,cen.center,P_B,k,Ei,times);//求cen.center到P_B的垂足pt_proj

		d = 0;
		bool new_x = false;
		for(i=0;i<Num;i++)
		{
			if(NoCov[i] && Target[i] != cen.kind)//样本未被覆盖且与pt不同类
			{
				bool tag = true;
				for(j=0;j<k;j++)//再排除样本属于P_B集合的情况
				{
					if(P_B[j] == i)
					{
						tag = false;
					}
				}
				if(tag)
				{
					double d_x;
					for(j=0;j<=Dim_in;j++)
					{
						temp[j] = Data[P_B[0]][j] - Data[i][j];//temp = c - x (c是P_B中任意一点)
					}
					d_x = GetInnerProduct(cen.center,temp,Dim_in+1) / GetInnerProduct(pt_proj,temp,Dim_in+1);

					if(d_x == 0.0)
					{
						new_x = false;
						P_B[k] = i;//add this point to the set P_B
						k++;
						break;//break for(i)
					}
					else
					{
						if(new_x)
						{
							if(d>d_x)
							{
								d = d_x;
								P_B[k] = i;
							}
						}
						else
						{
							d = d_x;
							P_B[k] = i;
						}
						new_x = true;
					}
				}//------end of if(tag)
			}//------end of if(Covered[i] && Target[i])
		}//-----end of for(i)
		if(new_x)
		{
			k++;
			double module = 0.0;
			for(j=0;j<=Dim_in;j++)
			{
				temp[j] = cen.center[j] - d * pt_proj[j];
				module += temp[j] * temp[j];
			}
			module = R / sqrt(module);//这里“模”已经改变了它本来的意义，便于简化下面的计算
			
			for(j=0;j<=Dim_in;j++)//将向量temp投影到球面上
			{
				cen.center[j] = temp[j] * module;
			}
		}
	}//-----end of while 

	delete[] Ei;
	delete[] temp;
	delete[] pt_proj;
	delete[] P_B;
	return 0;
}

///////////求与cen距离最近的非本类的样本点集，返回值为样本点的个数/////////////
int CoverNN::GetSet_minDist(const CovCell &cen, int* P_B)
{
	int counter = 0;
	int i;
	double tmp;
	double dist = 2*R;//距离初始化一个最大值，等于投影球面的直径

	for(i=0;i<Num;i++)
	{
		if(NoCov[i] && Target[i] != cen.kind)//样本未被覆盖且与cen不同类
		{
			tmp = GetDistance(cen.center,Data[i],Dim_in+1);//使用欧氏距离计算
			
			if(dist>tmp)
			{
				counter=1;
				dist = tmp;
				P_B[0] = i;//0 == counter-1
			}
			else if(dist==tmp)
			{
				if(counter < Dim_in+1)
				{
					P_B[counter] = i;//P_B有Dim_in+1个位置，多余的就没有必要再放了
				}
				counter++;
			}
		}
	}

	return counter;
}

////////////////////求投影点(方法1、2的效果是等价的)/////////////////
int CoverNN::GetProjection(double *pt, const double *a, const int *set, const int k, double **Ei, int &times)
{
	int i,j,l;
	double * d = new double[Dim_in+1];
	double tmp;

	//------------先求归一化的正交基---------------
	if(times == 0)//第一次调用该函数
	{
		for(i=1;i<k;i++)
		{
			for(j=0;j<=Dim_in;j++)
			{
				d[j] = Data[set[i]][j] - Data[set[0]][j];//求得di
			}
			for(j=1;j<i;j++)
			{
				tmp = GetInnerProduct(d,Ei[j],Dim_in+1);
				for(l=0;l<=Dim_in;l++)
				{
					d[l] -= tmp * Ei[j][l];
				}
			}
			tmp = GetModule(d,Dim_in+1);
			for(j=0;j<=Dim_in;j++)
			{
				Ei[i][j] = d[j] / tmp;//归一化求得Ei[i]
			}
		}
		//-----------再求投影------------
/*1		for(i=0;i<=Dim_in;i++)
		{
			d[i] = a[i] - Data[set[0]][i];//求得a - c0
			pt[i] = Data[set[0]][i];//初始化赋值pt=c0
		}
		for(i=1;i<k;i++)
		{
			tmp = GetInnerProduct(d,Ei[i],Dim_in+1);
			for(j=0;j<=Dim_in;j++)
			{
				pt[j] += tmp * Ei[i][j];
			}
		}
*/		//---------------求得投影点pt----------------
	}
	else//已经不是第一次调用该函数，k肯定大于1
	{
		//------------先求归一化的正交基---------------
		for(i=0;i<=Dim_in;i++)//求得d（即算法中的dk）
		{
			d[i] = Data[set[k-1]][i] - Data[set[0]][i];
		}
		for(i=1;i<k-1;i++)//求得未归一化的e(k-1)，存于d
		{
			tmp = GetInnerProduct(d,Ei[i],Dim_in+1);
			for(j=0;j<=Dim_in;j++)
			{
				d[j] -= tmp * Ei[i][j];
			}
		}
		tmp = GetModule(d,Dim_in+1);
		for(i=0;i<=Dim_in;i++)
		{
			Ei[k-1][i] = d[i] / tmp;//归一化求得Ei[k-1]
		}
		//-----------再求投影(之前所求的pt依然有效)------------
/*1		for(i=0;i<=Dim_in;i++)
		{
			d[i] = a[i] - pt[i];//求得a - c0
		}
		for(i=1;i<k;i++)
		{
			tmp = GetInnerProduct(d,Ei[i],Dim_in+1);
			for(j=0;j<=Dim_in;j++)
			{
				pt[j] += tmp * Ei[i][j];
			}
		}
*/		//---------------求得投影点pt----------------
	}
	//-----------再求投影(因a改变，之前所求的pt已经无效)------------
	for(i=0;i<=Dim_in;i++)
	{
		d[i] = a[i] - Data[set[0]][i];//求得a - c0
		pt[i] = Data[set[0]][i];//初始化赋值pt=c0
	}
	for(i=1;i<k;i++)
	{
		tmp = GetInnerProduct(d,Ei[i],Dim_in+1);
		for(j=0;j<=Dim_in;j++)
		{
			pt[j] += tmp * Ei[i][j];
		}
	}
/*2*/	//---------------求得投影点pt----------------
	times++;

	delete[] d;
	return 0;
}

//////////计算向量的模///////////////
inline double CoverNN::GetModule(double *vec, int dim)
{
	double module = 0;
	for(int i=0;i<dim;i++)
	{
		module += vec[i] * vec[i];
	}
	return sqrt(module);
}

//////////应用网络进行识别/////////////
double CoverNN::Recognize()
{
	CovCell * p;
	double tmp;
	int i,j,counter;
	int * result = new int[Num];//记录识别结果，用于和Target比较

	//------------将数据投影到球面上---------------
	for(i=0;i<Num;i++)
	{
		tmp = 0;
		for(j=0;j<Dim_in;j++)
		{
			tmp += Data[i][j]*Data[i][j];
		}
		Data[i][Dim_in] = sqrt(R*R - tmp);
	}
	//------------识别测试-------------
	for(i=0;i<Num;i++)
	{
		p = Covers->next;//p初始化为第一个覆盖单元
		while(p != NULL)
		{
			if( IsCovered(*p, Data[i]) )
			{
				result[i] = p->kind;
				break;
			}
			else
			{
				p = p->next;
			}
		}
	}
	//--------得到识别结果result,进行识别率统计----------------
	ofstream out("Notrecog.txt");
	counter = 0;
	for(i=0;i<Num;i++)
	{
		if(result[i] == Target[i])
		{
			counter++;
		}
		else
		{
			out<<Data[i][0]<<","<<Data[i][1]<<" ";//专为CAD中画图所写文件的格式
		}
	}
	out.close();

	tmp = (double)counter / (double)Num;//求得识别正确率

	delete[] result;

	return tmp;
}

/////////////////////////////////////////////////
/////////////神经网络数据文件格式
//LINE1: InnPro Num_C Dim_in R
//LINEi: CovCell.center CovCell.threshold CovCell.kind
/////////////////////////////////////////////////
/////////导出神经网络数据////////////
int CoverNN::ExportNNData()
{
	CovCell * p = Covers->next;
	ofstream out("ann.txt");

	out<<InnPro<<" "<<Num_C<<" "<<Dim_in<<" "<<R<<endl;
	while(p!=NULL)
	{
		for(int i=0;i<=Dim_in;i++)
		{
			out<<p->center[i]<<" ";
		}
		out<<p->threshold<<" "<<p->kind<<endl;//

		p = p->next;
	}

	out.close();

	return 0;
}

//////////载入神经网络数据///////////
int CoverNN::LoadNNData()
{
	CovCell * p = Covers;
	CovCell * q;
	ifstream inf("ann.txt");

	//载入数据前先执行清理
	if(p->next != NULL)
	{
		q = p->next;
		p->next = NULL;
		p = q;
		while(q != NULL)
		{
			p = p->next;
			delete[] q->center;
			delete q;
			q = p;
		}
		p = Covers;//依然需要为p重新赋值
	}

	int inn;
	inf>>inn>>Num_C>>Dim_in>>R;
	if(inn == 0)
		InnPro = false;
	else
		InnPro = true;
	for(int i=0;i<Num_C;i++)
	{
		q = new CovCell;
		q->center = new double[Dim_in+1];

		for(int j=0;j<=Dim_in;j++)
		{
			inf>>q->center[j];
		}
		inf>>q->threshold>>q->kind;

		p->next = q;
		p = p->next;
	}

	inf.close();

	return 0;
}

/////////////求一系列样本点中，最远的两点的中点(用内积计算)/////////////
int CoverNN::GetMidpoint(double *midpt, const bool *sfx)
{
	double tmp, dist = R * R;
	int endpt[2];
	int i,j;
	bool found = false;

	for(i=0;i<Num;i++)
	{
		if(sfx[i])
		{
			for(j=i+1;j<Num;j++)
			{
				if(sfx[j])
				{
					tmp = GetInnerProduct(Data[i],Data[j],Dim_in+1);
					if(dist>tmp)
					{
						endpt[0] = i;
						endpt[1] = j;
						dist = tmp;
						found = true;
					}
				}
			}
		}
	}
	if(found)
	{
		tmp = 0;
		//这种投影方式使其尽可能靠近螺旋线中心（有点耍赖）
		for(i=0;i<Dim_in;i++)
		{
			midpt[i] = ( Data[endpt[0]][i] + Data[endpt[1]][i] ) / 2;
			tmp += pow( midpt[i] , 2 );
		}
		midpt[Dim_in] = sqrt( R*R - tmp );

		return 0;//正常返回中心
	}
	else
	{
		return 1;//未达到预期结果返回，保持原“重心”
	}
}

/////////专用于平面螺旋线（范围：[PI/2 , 6*PI] ）///////////////
int CoverNN::ProductTestData(int n, int num_each, double disturb, unsigned int seed)
{
	const double PI = 3.1415926535;

	double x, rnd, dif = 2*PI / n;
	double cosx, sinx;
	int i, j, s, counter = 0;
	int num;
	
	num = 0;
	srand(seed);
	for(i=0;i<n;i++)
	{
		for(j=0;j<num_each;j++)
		{
			x = PI/2 + (j*11*PI) / (2*num_each);//缺点：没有保证在整条线上点的分布是均匀的（越靠外越稀疏）
			cosx = cos(x - i * dif);
			sinx = sin(x - i * dif);
			//---------添加随机扰动---------
			s = rand()%100;
			if(s%2 == 0)
			{
				rnd = disturb * s * PI / (n*100);
			}
			else
			{
				rnd = - disturb * s * PI / (n*100);
			}
			Data[num][0] = (x + rnd) * cosx;
			Data[num][1] = (x + rnd) * sinx;
			Data[num][2] = sqrt(R*R - Data[num][0]*Data[num][0] - Data[num][1]*Data[num][1]);
			Target[num] = i;
			num++;
		}
	}

	return 0;
}

int CoverNN::ProductTrainData(int n)
{
	const double PI = 3.1415926535;
	const int n_max = 78;//每条螺旋线的训练样本数

	ofstream out("train.txt");
	int i, j, k;
	double x, x1, y1, dif = 2*PI / n;//样本产生过程中角度的增量

	out<<n*n_max<<" "<<2<<endl;//写入样本总数和维数

	for(i=1;i<12;i++)
	{
		x = i * PI / 2;
		for(k=0;k<n;k++)
		{
			x1 = x * cos(x-k*dif);
			y1 = x * sin(x-k*dif);
			out<<x1<<" "<<y1<<" "<<k<<endl;
		}
		for(j=0;j<i;j++)
		{
			x = x + PI / (2 * (i+1));
			for(k=0;k<n;k++)
			{
				x1 = x * cos(x-k*dif);
				y1 = x * sin(x-k*dif);
				out<<x1<<" "<<y1<<" "<<k<<endl;
			}
		}
	}
	x = 6 * PI;
	for(i=0;i<n;i++)
	{
		x1 = x * cos(x-i*dif);
		y1 = x * sin(x-i*dif);
		out<<x1<<" "<<y1<<" "<<i<<endl;
	}
	out.close();

	return 0;
}

double CoverNN::Test(int n, int num_each, double disturb, unsigned int seed)
{
	double rate;
	Num = n*num_each;
	Dim_in = 2;

	Data = new double*[Num];
	Target = new int[Num];
	for(int i=0;i<Num;i++)
	{
		Data[i] = new double[3];//即Dim_in+1，因为专为平面螺旋线
	}

	ProductTestData(n, num_each, disturb, seed);
	rate = Recognize();

	//------------清理测试数据所占的内存--------------
	delete[] Data;
	delete[] Target;

	return rate;
}