面向对象的模拟退火编程技术.cpp

由于K－均值聚类算法局部最优的特点

void SAA::NewCenterPlus(ClusterType *p1,int t,double *p2,int dim)
{
	int i;
	for(i=0;i<dim;i++)
		p1[t].center[i]=p1[t].center[i]+(p2[i]-p1[t].center[i])/(p1[t].sonnum);
}


void SAA::NewCenterReduce(ClusterType *p1,int t,double *p2,int dim)
{
	int i;
	for(i=0;i<dim;i++)
		p1[t].center[i]=p1[t].center[i]+(p1[t].center[i]-p2[i])/(p1[t].sonnum);
}


void SAA::GetDataset(DataType *p1,double *p2,int datanum,int dim)
{
	int i,j;
	for(i=0;i<datanum;i++)
	{
		
		for(j=0;j<dim;j++)
			p1[i].data[j]=p2[i*dim+j];
	}
}

void SAA::GetValue(double *str1,double *str2,int dim)
{
	int i;
	for(i=0;i<dim;i++)		
		str1[i]=str2[i];
}

int  SAA::FindFather(double *p,int k)
{
	int i,N=0;
	double min=30000;

	for(i=0;i<k;i++)	
		if(p[i]<min)
		{
			min=p[i];
			N=i;
		}
	return N;
}

double SAA::SquareDistance(double *str1,double *str2,int dim)
{
	double dis=0;
	int i;
	for(i=0;i<dim;i++)
		dis=dis+(double)(str1[i]-str2[i])*(str1[i]-str2[i]);
	return dis;
}

int	SAA::Compare(double *p1,double *p2,int dim)
{
	int i;
	for(i=0;i<dim;i++)
		if(p1[i]!=p2[i])
			return 1;
	return 0;
}		



//******************************************
//              模拟退火算法   *************
//**************             ***************
//******************************************

void SAA::SA( )
{ 
	int i;
	int HALT,Count=1,TopNum=1,GetNum=0,AcceptNum=0;
	int InnerLoop,OuterLoop,TempNum=0;
	double gap,gap1,Temp0,TempRate=0;

	HALT=FALSE;
//	cout<<endl<<"temp: "<<Temp;
	CurrentStatus=new DataType[DataNum];
	NewStatus=new DataType[DataNum];
	CurrentCluster=new ClusterType[K];
	NewCluster=new ClusterType[K];

	Temp0=Temp=(double)CO*AimFunction(DataMember,ClusterMember);
//	cout<<endl<<"temp: "<<Temp;
//	Generate(DataMember,ClusterMember);

	for(i=0;i<DataNum;i++)
	{
		NewStatus[i].data=new double[Dimension];	
		CurrentStatus[i].data=new double[Dimension];	

		NewStatus[i].uncle=new double[K];
		CurrentStatus[i].uncle=new double[K];
	}
	
	for(i=0;i<K;i++)
	{
		CurrentCluster[i].center=new double[Dimension];
		NewCluster[i].center=new double[Dimension];
	}

	cout<<endl<<endl<<" ***********开 始 退 火！***************"
		<<endl<<" 退火起始温度： "<<Temp<<endl;


	OuterLoop=1;
	while(HALT==FALSE)
	{
		cout<<endl<<"              已连续降温次数: "<<OuterLoop<<"  T温度："<<Temp;
		CopyStatus(CurrentStatus,CurrentCluster,DataMember,ClusterMember);
		CopyStatus(NewStatus,NewCluster,DataMember,ClusterMember);
//		cout<<endl<<"  Newfun: "<<AimFunction(NewStatus,NewCluster)<<" Old Aim: "<<AimFunction(CurrentStatus,CurrentCluster);
		InnerLoop=1;

		for(i=0;i<MarkovLengh;i++)
		{
		    if(TopNum%100000==0)
			{
				cout<<endl<<endl<<"  "<<TopNum<<"次   系统正在计算，请耐心等待！"<<endl;
//				cout<<endl<<"gap:"<<gap<<" Temp: "<<Temp<<"  Rate: "<<Temp/Temp0;
			}
			TopNum+=1;
			Generate(NewStatus,NewCluster);
			GetNum=GetNum+1;
//			cout<<endl<<"  NewStatus: "<<AimFunction(NewStatus,NewCluster)<<" OldStatus: "<<AimFunction(CurrentStatus,CurrentCluster);

			gap=AimFunction(NewStatus,NewCluster)-AimFunction(CurrentStatus,CurrentCluster);
			if(gap<0)
			{
				CopyStatus(CurrentStatus,CurrentCluster,NewStatus,NewCluster);
				InnerLoop=+1;
				AcceptNum=AcceptNum+1;
//				cout<<endl<<"****************第"<<TopNum<<"次运算";
			}
			else if(exp(-gap/Temp)>(double)rand()/(double)RAND_MAX)
			{
				CopyStatus(CurrentStatus,CurrentCluster,NewStatus,NewCluster);
				InnerLoop+=1;
				AcceptNum=AcceptNum+1;
				TempRate=TempRate+exp(-gap/Temp);
				TempNum+=1;
//				cout<<endl<<"****************第"<<TopNum<<"次运算*";
			}
			else 
			{
				CopyStatus(NewStatus,NewCluster,CurrentStatus,CurrentCluster);
				InnerLoop+=1;
//				cout<<endl<<"****************第"<<TopNum<<"次运算";
			}

			gap1=AimFunction(NewStatus,NewCluster)-AimFunction(DataMember,ClusterMember);
			if(gap1<0)
			{
			CopyStatus(DataMember,ClusterMember,NewStatus,NewCluster);
			CopyStatus(CurrentStatus,CurrentCluster,NewStatus,NewCluster);
			OuterLoop=1;
			InnerLoop=1;
			cout<<endl<<" ******************* 成功改进"<<Count<<"次聚类*************************************";
//			cout<<endl<<"  ChangeTemp:"<<Temp
//				<<"  Rate:"<<Temp/Temp0<<"  接受率"<<(double)AcceptNum/(double)GetNum;
			Count=Count+1;
			cout<<endl<<" 改进后准则函数: "<<AimFunction(DataMember,ClusterMember)<<endl;
			}
			if(GetNum%1000==0)
			{
//				cout<<endl<<"  接受率"<<(double)AcceptNum/(double)GetNum;
				 if((double)AcceptNum/(double)GetNum>AcceptRate)
					 Temp=Temp*ForceDecline;
				 GetNum=0;
				 AcceptNum=0;
			}

//		cout<<endl<<" Innter Loop:********************:  "<<InnerLoop;
//		cout<<endl<<" gap:****************:  "<<gap1<<endl;
			if(InnerLoop>MaxInnerLoop)
			break;

		}// end of InnerLoop/for

		OuterLoop+=1;
//		cout<<endl<<endl<<" before: "<<Temp;
		Temp=Temp*DeclineRate;
//		cout<<endl<<" After: "<<Temp;

		if(OuterLoop>MaxOuterLoop)
		HALT=1;
//		if(GetNum==1000)
//			break;
	}//end of while
	cout<<endl<<"  TempRate:"<<TempRate/TempNum;


	cout<<endl<<"  退火次数: "<<OuterLoop-1
		<<endl<<"  终止温度："<<Temp<<endl;
	

}

//将样本分成一类，将其准则函数作为初始温度TEMP。

double SAA::MaxFunc()
{
	int i,j;
	double *mean,*dis;
	Temp=0;

	mean=new double[Dimension];
	dis=new double[Dimension];
	for(i=0;i<Dimension;i++)
	{
		mean[i]=0;
		dis[i]=0;
	}

	for(i=0;i<DataNum;i++)
		for(j=0;j<Dimension;j++)
			mean[j]=mean[j]+DataMember[i].data[j];

	for(i=0;i<Dimension;i++)
		mean[i]=mean[i]/DataNum;

	for(i=0;i<DataNum;i++)
		for(j=0;j<Dimension;j++)
			dis[j]=dis[j]+(DataMember[i].data[j]-mean[j])*(DataMember[i].data[j]-mean[j]);
	for(i=0;i<Dimension;i++)
		Temp=Temp+dis[i];	
	return Temp;
}

//按照某种规则，生成新的状态空间

void SAA::Generate(DataType *p1,ClusterType *c1)
{
	int ra,secfa,fa;

	ra=(int)(((double)rand()/(double)RAND_MAX)*DataNum);
	fa=p1[ra].father;
//	cout<<"father: "<<c1[fa].sonnum;

	while(ra==DataNum||c1[fa].sonnum==1)    
	{
		ra=(int)( (double)rand()/(double)RAND_MAX*(double)DataNum);
		fa=p1[ra].father;
	}

//	secfa=SecondFather(p1,ra,K);
	secfa=(int)(((double)rand()/(double)RAND_MAX)*K);
	while(secfa==K||secfa==fa)
		secfa=(int)(((double)rand()/(double)RAND_MAX)*K);

//	cout<<endl<<"sera: "<<secfa<<endl;
//	cout<<endl<<" "<<p1[ra].data[0]<<" in "<<p1[ra].father+1<<""<<" to secfa:"<<secfa+1<<endl;

	p1[ra].father=secfa;
	c1[fa].sonnum-=1;
	c1[secfa].sonnum+=1;
	NewCenterReduce(c1,fa,p1[ra].data,Dimension);
//	cout<<endl<<"new father's center: "<<c1[secfa].center[0];
	NewCenterPlus(c1,secfa,p1[ra].data,Dimension);
//	cout<<endl<<"new father's center: "<<c1[secfa].center[0];

}



void SAA::CopyStatus(DataType *p1,ClusterType *c1,DataType *p2,ClusterType *c2)
{
	int i,j;
	for(i=0;i<DataNum;i++)
	{
		for(j=0;j<Dimension;j++)
			p1[i].data[j]=p2[i].data[j];
		p1[i].father=p2[i].father;
		for(j=0;j<K;j++)
			p1[i].uncle[j]=p2[i].uncle[j];
	}

	for(i=0;i<K;i++)
	{
		for(j=0;j<Dimension;j++)
			c1[i].center[j]=c2[i].center[j];
		c1[i].sonnum=c2[i].sonnum;
	}

}

int SAA::SecondFather(DataType *p,int t,int k)
{
	int i,fa,secfa;
	double min;
	secfa=fa=p[t].father;
	min=p[t].uncle[fa];
	for(i=0;i<K;i++)
	{
//		cout<<"  "<<t<<"->:"<<i+1<<": "<<p[t].uncle[i];
		if(secfa==fa&&p[t].uncle[i]>min)
			secfa=i;
		else
			if(p[t].uncle[i]>min&&p[t].uncle[i]<p[t].uncle[secfa])
				secfa=i;	
	}
	
//	cout<<endl<<"  fa:"<<fa+1<<" secfa:"<<secfa+1<<endl<<endl;
//	cout<<endl<<endl;

	return secfa;
}

double SAA::FRand(double a,double b)
{
	return a+(double)(((double)rand()/(double)RAND_MAX)*(b-a));

}

void SAA::DisPlay()
{
	int i,N,j,t;
	ofstream  result("聚类过程结果显示.txt",ios::ate);
	for(i=0;i<K;i++)
	{
		N=0;
		cout<<endl<<endl<<"******************** 第 "<<i+1<<" 类样本:*******************"<<endl;
		result<<endl<<endl<<"******************** 第 "<<i+1<<" 类样本:*******************"<<endl;
		for(j=0;j<DataNum;j++)
			if(DataMember[j].father==i)
			{
				cout<<" [";
				for(t=0;t<Dimension;t++)
				cout<<" "<<setw(5)<<DataMember[j].data[t];
				cout<<" ]  ";				
				if((N+1)%Row==0)
					cout<<endl;

				result<<" [";
				for(t=0;t<Dimension;t++)
				result<<" "<<setw(5)<<DataMember[j].data[t];
				result<<" ]  ";				
				if((N+1)%Row==0)
					result<<endl;

				N=N+1;
			}
	}//end for

	cout<<endl<<endl<<"  聚类结果，总体误差准则函数："<<AimFunction(DataMember,ClusterMember)<<endl;
	result<<endl<<"  聚类结果，总体误差准则函数："<<AimFunction(DataMember,ClusterMember)<<endl;
  
	result.close();

}//end of Display

double SAA::AimFunction(DataType *q,ClusterType *c)
{
	int i,j;
	double *p;
	p=new double[K];
	for(i=0;i<K;i++)
		p[i]=0;
	for(i=0;i<K;i++)
	{
		for(j=0;j<DataNum;j++)
			if(q[j].father==i)
			{
				p[i]=p[i]+SquareDistance(c[i].center,q[j].data,Dimension);
			}
	}

			
	AimFunc=0;
	for(i=0;i<K;i++)
		AimFunc=AimFunc+p[i];
	return AimFunc;

}



//************************************
//            主函数入口          ****   
//************************************


void main()
{
	//用户输入数据
	srand((unsigned)time(NULL));
	GETDATA getdata;
	getdata.Initial();
	ofstream file("聚类过程结果显示.txt",ios::trunc);

	//k-均值聚类方法聚类
	SAA saa;    //****此行不可与上行互换。
	
	saa.KMeans();
	cout<<endl<<"***********************K-均值聚类结果：**********************";
	file<<endl<<"***********************K-均值聚类结果：**********************"<<endl;

	file.close();
	saa.DisPlay();

	//模拟退火，进一步聚类
	saa.SA();
	saa.DisPlay();

//	saa.KMeans1();
//	saa.DisPlay();
	
//	cout<<endl<<"Temp="<<saa.MaxFunc();

	/*cout<<endl<<" 数据输入完毕  开始聚类 按B键   退出  按Q键：";		 
		cin>>begin;
		if(begin=='b'||begin=='B')
			saa.SA();
		else
			if(begin=='q'||begin=='Q')
				cout<<endl<<" 程序已退出！"<<endl;
			else 
				cout<<endl<<" 输入错误，程序停止！"<<endl;
	*/

	cout<<endl<<"  程序运行结束！"<<endl;
			
}