pca.cpp

主元分析PCA的C代码

#include "PCA.h"

//求实对称矩阵的特征值及特征向量的雅格比法
//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量
//返回值小于0表示超过迭代jt次仍未达到精度要求
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值
//n-矩阵的阶数
//u-长度为n*n的数组，返回特征向量(按列存储)
//eps-控制精度要求
//jt-整型变量，控制最大迭代次数
double eejcb(double a[],int n,double v[],double eps,int jt)
{ 
	int i,j,p,q,u,w,t,s,l;
    double fm,cn,sn,omega,x,y,d;
    l=1;
    for (i=0; i<=n-1; i++)
	{ 
		v[i*n+i]=1.0;
		for (j=0; j<=n-1; j++)
		{
			if (i!=j) 
			{
				v[i*n+j]=0.0;
			}
		}
	}
    while (1==1)
	{ 
		fm=0.0;
        for (i=0; i<=n-1; i++)
		{
			for (j=0; j<=n-1; j++)
			{ 
				d=fabs(a[i*n+j]);
				if ((i!=j)&&(d>fm))
				{ 
					fm=d; 
					p=i; 
					q=j;
				}
			}
		}
        if (fm<eps)  
		{
			return l;
		}
        if (l>jt)  
		{
			return l;
		}
        l=l+1;
        u=p*n+q; 
		w=p*n+p; 
		t=q*n+p; 
		s=q*n+q;
        x=-a[u];
		y=(a[s]-a[w])/2.0;
        omega=x/sqrt(x*x+y*y);
        if (y<0.0)
		{
			omega=-omega;
		}
        sn=1.0+sqrt(1.0-omega*omega);
        sn=omega/sqrt(2.0*sn);
        cn=sqrt(1.0-sn*sn);
        fm=a[w];
        a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;
        a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;
        a[u]=0.0;
		a[t]=0.0;
        for (j=0; j<=n-1; j++)
		{
			if ((j!=p)&&(j!=q))
			{ 
				u=p*n+j;
				w=q*n+j;
				fm=a[u];
				a[u]=fm*cn+a[w]*sn;
				a[w]=-fm*sn+a[w]*cn;
			}
		}
        for (i=0; i<=n-1; i++)
		{
			if ((i!=p)&&(i!=q))
            { 
				u=i*n+p; 
				w=i*n+q;
				fm=a[u];
				a[u]=fm*cn+a[w]*sn;
				a[w]=-fm*sn+a[w]*cn;
            }
		}
        for (i=0; i<=n-1; i++)
		{ 
		   u=i*n+p; 
		   w=i*n+q;
           fm=v[u];
           v[u]=fm*cn+v[w]*sn;
           v[w]=-fm*sn+v[w]*cn;
		}
	}
	return l;
}

void AllocPCAData(PCA_Data &data)
{
	//FreePCAData(data);
	if(data.number_sample>0)
	{
		data.x=(double **)malloc(data.number_sample*sizeof(double *));
		data.y=(double *)malloc(data.number_sample*sizeof(double));
	}
	for(int i=0; i<data.number_sample; i++)
		data.x[i]=(double *)malloc(data.number_input*sizeof(double));
}

void FreePCAData(PCA_Data &data)
{
	if(data.x)
	{
		for(int i=0; i<data.number_sample; i++)
			if(data.x[i])free(data.x[i]);
		free(data.x);
	}
	if(data.y)free(data.y);
}

int PCA(PCA_Param param, PCA_Data o_data)
{
	int i, j, k;
	//协方差矩阵
	double *cov_ma=(double *)malloc(o_data.number_input*o_data.number_input*sizeof(double));
	//特征向量
	double *eigen_vector=(double *)malloc(o_data.number_input*o_data.number_input*sizeof(double));
	//样本集均值
	double *mean_value=(double *)malloc(o_data.number_input*sizeof(double));
	FILE *fp;
	char str[MAX_NAME_LEN];

	//求各特征均值
	for(i=0; i<o_data.number_input; i++)
	{
		for(j=0, mean_value[i]=0.; j<o_data.number_sample; j++)
			mean_value[i]+=o_data.x[j][i];
		mean_value[i]/=o_data.number_sample;
	}
	//求协方差矩阵
	for(i=0; i<o_data.number_input; i++)
	{
		for(j=i; j<o_data.number_input; j++)
		{
			for(k=0, cov_ma[i*o_data.number_input+j]=0.; k<o_data.number_sample; k++)
				cov_ma[i*o_data.number_input+j]+=(o_data.x[k][i]-mean_value[i])*(o_data.x[k][j]-mean_value[j]);
			cov_ma[i*o_data.number_input+j]/=(o_data.number_sample-1);
		}
	}
	//协方差矩阵左下部
	for(i=0; i<o_data.number_input; i++)
	{
		for(j=0; j<i; j++)
			cov_ma[i*o_data.number_input+j]=cov_ma[j*o_data.number_input+i];
	}
	//求协方差矩阵特征值与特征向量
	printf("迭代精度：%lf\n", eejcb(cov_ma, o_data.number_input, eigen_vector, param.eps, param.max_iter));
	sprintf(str, "result/%s/evalue.txt", param.data_set_name);
	fp=fopen(str, "w");
	for(i=0; i<o_data.number_input; i++)
		fprintf(fp, "%lf\n", cov_ma[i*o_data.number_input+i]);
	fclose(fp);
	//挑出前n个特征值及相应的特征向量
	unsigned short *index=(unsigned short *)malloc(param.number_pc*sizeof(unsigned short));
	double max_v;
	for(i=0; i<param.number_pc; i++)
	{
		//一趟冒泡，找出当前最大特征值
		for(j=0, max_v=0.; j<o_data.number_input; j++)
		{
			if(cov_ma[j*o_data.number_input+j]>max_v)
			{
				index[i]=j;								//记录该最大特征值的序号，即特征的序号
				max_v=cov_ma[j*o_data.number_input+j];
			}
		}
		//将本次找大的最大特征值置零，即下一趟不再参与查找
		cov_ma[index[i]*o_data.number_input+index[i]]=0.;
	}

	//保存结果
	sprintf(str, "result/%s/model.txt", param.data_set_name);
	fp=fopen(str, "w");
	//输入数、主元数
	fprintf(fp, "%d %d\n", param.number_input, param.number_pc);
	//各输入的均值
	for(i=0; i<o_data.number_input; i++)
		fprintf(fp, "%lf ", mean_value[i]);
	fprintf(fp, "\n");
	//选出的特征向量
	for(j=0; j<param.number_pc; j++)
	{
		for(i=0; i<o_data.number_input; i++)
			fprintf(fp, "%lf ", eigen_vector[i*o_data.number_input+index[j]]);
		fprintf(fp, "\n");
	}
	fclose(fp);

	return 1;
}

void DataTransform(PCA_Param param, PCA_Data o_data)
{
	FILE *fp;
	int number_input, number_pc;
	char str[MAX_NAME_LEN];
	sprintf(str, "result/%s/model.txt", param.data_set_name);
	fp=fopen(str, "r");
	if(!fp)
	{
		printf("模型文件打开错误！\n");
		return;
	}
	//输入数、主元数
	fscanf(fp, "%d", &number_input);
	fscanf(fp, "%d", &number_pc);
	if(number_input!=o_data.number_input)
	{
		printf("模型文件中的特征数与给定数据集的特征数不符！");
		fclose(fp);
		return;
	}

	int i, j, k;
	//样本集均值
	double *mean_value=(double *)malloc(o_data.number_input*sizeof(double));
	for(i=0; i<o_data.number_input; i++)
		fscanf(fp, "%lf ", &mean_value[i]);
	//特征向量
	double **eigen_vector=(double **)malloc(number_pc*sizeof(double *));
	for(j=0; j<param.number_pc; j++)
	{
		eigen_vector[j]=(double *)malloc(o_data.number_input*sizeof(double));
		for(i=0; i<o_data.number_input; i++)
			fscanf(fp, "%lf ", &eigen_vector[j][i]);
	}
	fclose(fp);

	PCA_Data res_data;
	res_data.number_input=number_pc;
	res_data.number_sample=o_data.number_sample;
	res_data.x=NULL;
	res_data.y=NULL;
	AllocPCAData(res_data);
	//计算生成新的数据集
	double *max_v=(double *)malloc(res_data.number_input*sizeof(double));
	double *min_v=(double *)malloc(res_data.number_input*sizeof(double));
	for(i=0; i<res_data.number_input; i++)
	{
		max_v[i]=-1000000.;
		min_v[i]=1000000.;
	}
	for(i=0; i<res_data.number_sample; i++)
	{
		for(j=0; j<res_data.number_input; j++)
		{
			res_data.x[i][j]=0.;
			for(k=0; k<o_data.number_input; k++)
				res_data.x[i][j]+=(o_data.x[i][k]-mean_value[k])*eigen_vector[j][k];
			max_v[j]=(max_v[j]>res_data.x[i][j])?max_v[j]:res_data.x[i][j];
			min_v[j]=(min_v[j]<res_data.x[i][j])?min_v[j]:res_data.x[i][j];
		}
	}
	memcpy(res_data.y, o_data.y, o_data.number_sample*sizeof(double));
	free(mean_value);
	for(j=0; j<param.number_pc; j++)
		free(eigen_vector[j]);
	free(eigen_vector);

	sprintf(str, "result/%s/%s.data", param.data_set_name, param.data_set_name);
	fp=fopen(str, "w");
	fprintf(fp, "%d\n", res_data.number_sample);
	//线性归一化
	for(i=0; i<res_data.number_sample; i++)
	{
		for(j=0; j<res_data.number_input; j++)
		{
			res_data.x[i][j]=(res_data.x[i][j]-min_v[j])/(max_v[j]-min_v[j]);
			fprintf(fp, "%lf ", res_data.x[i][j]);
		}
		fprintf(fp, "%lf\n", res_data.y[i]);
	}
	fclose(fp);
	FreePCAData(res_data);
}