sofm.cpp

自组织的算法 学神经网络的必会算法 很重要

#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#include "time.h"

//求两个2维向量之间的模
double dist(double *x1,double*x2)
{
	return(pow(pow(x1[0]-x2[0],2)+pow(x1[1]-x2[1],2),0.5));
}

//将映射层神经元代表的样本聚类结果显示为
//第几类（代表的此类样本个数）
//第一次迭代时创建新的名为name的结果文件，之后继续写
//M映射层神经元行数
//N映射层神经元列数
//classnum样本总类数
//x存放神经元分类结果的矩阵指针
//c是否新建保存结果的文件
void display(int M,int N,int classnum,int ***x,char name[20],char c)
{
	int i,j,k,cache;
	FILE *ep;

	if(c=='c')
	{
		if((ep=fopen(name,"wt+"))==NULL)
		{
			printf("\nSorry can not write to the file %s!\n",name);
			exit(0);
		}
	}
	else
	{
		if((ep=fopen(name,"at+"))==NULL)
		{
			printf("\nSorry can not write to the file %s!\n",name);
			exit(0);
		}
	}

	printf("\n\n----------------------------------------------");
	fprintf(ep,"\n\n----------------------------------------------");
	printf("\n     ");
	fprintf(ep,"\n     ");
	for(j=0;j<N;j++)
	{
		printf("-%d-   ",j+1);
		fprintf(ep,"-%d-   ",j+1);
	}
	for(i=0;i<M;i++)
	{
		printf("\n-%d-  ",i+1);
		fprintf(ep,"\n-%d-  ",i+1);
		
		for(j=0;j<N;j++)
		{
			//求出该神经元代表哪类样本最多
			cache=0;
			for(k=0;k<classnum;k++)
				if(x[i][j][k]>x[i][j][cache]) cache=k;
			//输出竞争优胜的神经元代表的结果
			if(x[i][j][cache]!=0)
			{
				printf("%d(%d)  ",cache+1,x[i][j][cache]);
				fprintf(ep,"%d(%d)  ",cache+1,x[i][j][cache]);
			}
			//不是竞争优胜的神经元不显示
			else
			{
				printf("----  ");
				fprintf(ep,"----  ");
			}
		}
	}
	fclose(ep);
}

void save_w(int M,int N,double ***w,char name[20])
{
	int i,j;
	FILE *ep;

	if((ep=fopen(name,"wt+"))==NULL)
	{
		printf("\nSorry can not write to the file %s!\n",name);
		exit(0);
	}
	
	for(i=0;i<M;i++)
	{
		for(j=0;j<N;j++)
			fprintf(ep,"\nw[%d][%d][0]=%f  w[%d][%d][1]=%f",i,j,w[i][j][0],i,j,w[i][j][1]);
		fprintf(ep,"\n");
	}
	fclose(ep);
}
//自组织映射网络函数
//samplenum样本个数
//classnum样本总类数
//M映射层神经元行数
//N映射层神经元列数
//iteration最大迭代次数
//x存放样本数据的数组指针
void sofm(int samplenum,int classnum,int M,int N,int iteration,double **x)
{
	int i,j,k,l,p,q,temp[2],range0,range,***xcla,cla,a,b;
	double ***w,**d,cache,cache1,cache2,dlast,dnew,c1,c2;

	//w是映射层神经元系数矩阵，每次都赋值为0到1之间的随机数
	//xcla存放映射层神经元代表每类样本个数
	w=(double***)malloc(M*sizeof(double**));
	xcla=(int***)malloc(M*sizeof(int**));
	{
		for(i=0;i<M;i++)
		{
			w[i]=(double**)malloc(N*sizeof(double*));
			xcla[i]=(int**)malloc(N*sizeof(int*));
			for(j=0;j<N;j++)
			{
				w[i][j]=(double*)malloc(2*sizeof(double));
				xcla[i][j]=(int*)malloc(classnum*sizeof(int));//////////
			}
		}
	}
	//d存放映射层每个神经元到1个样本的距离
	d=(double**)malloc(M*sizeof(double*));
	for(i=0;i<M;i++)
		d[i]=(double*)malloc(N*sizeof(double));
	
	//权值赋初值
	srand( (unsigned)time( NULL ) );
	for(j=0;j<M;j++)
		for(i=0;i<N;i++)
		{
			w[j][i][0]=(double)rand()/32767;
			w[j][i][1]=(double)rand()/32767;
			for(k=0;k<classnum;k++) xcla[j][i][k]=0;
		}
	
	//确定调整范围最大值range0
	frexp(N*M,&range0);
	range0--;
	//迭代开始
	for(l=0;l<iteration;l++)
	{
		for(j=0;j<M;j++)
			for(i=0;i<N;i++)
				for(k=0;k<classnum;k++) xcla[i][j][k]=0;
		//本次迭代映射层所有神经元到所有样本的距离总合
		dnew=0;
		//本次迭代调整范围range
		frexp(l/100+1,&range);
		if(range<=range0) range=range0-range;
		else range=0;
		for(k=0;k<samplenum;k++)
		{
			//映射层所有神经元到一个样本的距离
			cla=(k)/5;
			for(i=0;i<M;i++)
				for(j=0;j<N;j++)
					d[i][j]=dist(w[i][j],x[k]);
				
			//找对该样本竞争获胜的神经元编号
			cache=d[0][0];
			temp[0]=0;temp[1]=0;
			for(i=0;i<M;i++)
				for(j=0;j<N;j++)
					if(d[i][j]<cache)
					{
						cache=d[i][j];
						temp[0]=i;temp[1]=j;
					}
			//记录竞争获胜的神经元对应该类样本的个数
			xcla[temp[0]][temp[1]][cla]++;
			//权值调整
			//范围是以获胜神经元为中心，range为半径的正方形
			//若调整范围超过M*N则到边界为止
			for(p=-range;p<range+1;p++)
			{
				//c1,c2是针对调整范围内每个神经元调整量大小的系数
				//以中心点向外逐渐递减
				if(range!=0) c1=1-(double)(fabs(p)/range);
				else c1=1;
				for(q=-range;q<range+1;q++)
				{
					if(range!=0) c2=1-(double)(fabs(q)/range);
					else c2=1;
					a=temp[0]+p;
					b=temp[1]+q;
					if(a<0) continue;
					else if(a>=M) continue;
					if(b<0) continue;
					else if(b>=N) continue;
					cache1=c1*c2*0.1/pow(l+1,0.2)*(x[k][0]-w[a][b][0]);
					cache2=c1*c2*0.1/pow(l+1,0.2)*(x[k][1]-w[a][b][1]);
					dnew+=fabs(cache1)+fabs(cache2);
					w[a][b][0]+=cache1;
					w[a][b][1]+=cache2;
				}
			}
		}
				
		printf("\n%d dnew=%f",l+1,dnew);
		//判断本次和上次迭代结果总距离变化量是否足够小
		//是，结束迭代，退出，否，继续迭代
		if(l==0)
		{
			display(M,N,classnum,xcla,"result_course.dat",'c');
			dlast=dnew;
		}
		else
		{
			display(M,N,classnum,xcla,"result_course.dat",'1');
	//		if(fabs(dnew-dlast)<10e-7)
			if(fabs(dnew)<10e-7)
			{
				printf("\nThe iteration has convergenced!\n");
				display(M,N,classnum,xcla,"result.dat",'c');
				save_w(M,N,w,"w.dat");
				exit(0);
			}
			else dlast=dnew;
		}
	}
	printf("\nThe iteration abort for not convergenced!\n");
	display(M,N,classnum,xcla,"result.dat",'c');
	save_w(M,N,w,"w.dat");
}

void main()
{
	int i;
	double **sample;

	sample=(double**)malloc(15*sizeof(double*));
	for(i=0;i<15;i++)
		sample[i]=(double*)malloc(2*sizeof(double));
	
	sample[0][0]=5;sample[0][1]=0;
	sample[1][0]=6;sample[1][1]=0;
	sample[2][0]=5;sample[2][1]=1;
	sample[3][0]=6;sample[3][1]=1;
	sample[4][0]=5.5;sample[4][1]=0.5;
	sample[5][0]=0;sample[5][1]=9;
	sample[6][0]=0;sample[6][1]=10;
	sample[7][0]=1;sample[7][1]=9;
	sample[8][0]=1;sample[8][1]=10;
	sample[9][0]=0.5;sample[9][1]=9.5;
	sample[10][0]=5;sample[10][1]=9;
	sample[11][0]=5;sample[11][1]=10;
	sample[12][0]=6;sample[12][1]=9;
	sample[13][0]=6;sample[13][1]=10;
	sample[14][0]=5.5;sample[14][1]=9.5;

	sofm(15,3,5,5,50000,sample);
}