classificationmqdf.cpp

该程序包实现了几个常用的模式识别分类器算法

/* CMQdf类实现改进二次判别方程算法（MQDF）。
   使用方法如下：创建一个对象，调用train方法训练分类器，调用test方法进行分类。 */
#include "stdlib.h"
#include <iostream>
#include <string.h>
#include <fstream>
#include <math.h>
using namespace std;
#include "global.h"

CMQdf::CMQdf()
{
	paramU = NULL;
	paramSigma = NULL;
	resultTotal = NULL;
	resultRight = NULL;
	paramInvSigma = NULL;
	paramLogSigma = NULL;
}

CMQdf::~CMQdf()
{
	if (this->paramU != NULL)
		delete[] this->paramU;
	if (this->paramSigma != NULL)
		delete[] this->paramSigma;
	if (this->paramInvSigma != NULL)
		delete[] this->paramInvSigma;
	if (this->paramLogSigma != NULL)
		delete[] this->paramLogSigma;
	if (this->resultTotal != NULL)
		delete[] this->resultTotal;
	if (this->resultRight != NULL)
		delete[] this->resultRight;
}

/* 训练分类器，输入训练样本的文件名，返回true表示训练成功，返回false表示训练失败。 */
bool CMQdf::train(char* fileName)
{
	bool ret = true;
	int i, j, k, l, index, cnum;
	ret = this->sdata.readFile(fileName);//读入训练数据
	if (!ret)
		return false;//如果数据格式不正确，退出程序。
	this->sdata.normalize();
	int numClass = this->sdata.numClass;
	int numFeature = this->sdata.numFeature;
	int numSample = this->sdata.numSample;
	DOUBLE *sumX = new DOUBLE[numClass*numFeature];
	DOUBLE *sumXX = new DOUBLE[numClass*numFeature*numFeature];
	int* classSample = new int[numClass];
	int total;
	DOUBLE temp1, temp2, maxv;
	for (i=0; i<numClass; i++)
	{
		classSample[i]=0;
	}
	total = numClass * numFeature;
	for (i=0; i<total; i++)
	{
		sumX[i]=0;
	}
	total = numClass * numFeature * numFeature;
	for (i=0; i<total; i++)
	{
		sumXX[i]=0;
	}

	for (i=0; i<numSample; i++)
	{
		cnum = this->sdata.ydata[i];
		classSample[cnum] ++;
		for (j=0; j<numFeature; j++)
		{
			sumX[cnum*numFeature+j] += this->sdata.xdata[i*numFeature+j];
			for (k=j; k<numFeature; k++)
			{
				sumXX[(cnum*numFeature+j)*numFeature+k] +=this->sdata.xdata[i*numFeature+j] * this->sdata.xdata[i*numFeature+k];
			}
		}
	}
	this->paramU = new DOUBLE[numClass*numFeature];
	this->paramSigma = new DOUBLE[numClass*numFeature*numFeature];
	for (i=0; i<numClass; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			this->paramU[i*numFeature+j] = sumX[i*numFeature+j] / classSample[i];
		}
	}
	for (i=0; i<numClass; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			for (k=j; k<numFeature; k++)
			{
				index = (i*numFeature+j)*numFeature+k;
				this->paramSigma[index] = sumXX[index]/classSample[i]-(this->paramU[i*numFeature+j]*this->paramU[i*numFeature+k]);
			}
		}
	}
	for (i=0; i<numClass; i++)
	{
		for (j=1; j<numFeature; j++)
		{
			for (k=0; k<j; k++)
			{
				this->paramSigma[(i*numFeature+j)*numFeature+k] = this->paramSigma[(i*numFeature+k)*numFeature+j];
			}
		}
	}
	/*for (i=0; i<numClass; i++)//限制协方差矩阵的对角元素的下限
	{
		for (j=0; j<numFeature; j++)
		{
			if (this->paramSigma[(i*numFeature+j)*numFeature+j]<0.001)
			{
				this->paramSigma[(i*numFeature+j)*numFeature+j] = 0.001;
			}
		}
	}*/
	CMatrix2 theMat;
	float** Cova;
	Cova = theMat.allocMat(numFeature);		// n: dimensionality of vector
	float** P;
	P = theMat.allocMat(numFeature);
	for (i=0; i<numClass; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			for (k=0; k<numFeature; k++)
			{
				Cova[j][k] = (float)this->paramSigma[(i*numFeature+j)*numFeature+k];
			}
		}
		theMat.diagonalize(Cova, numFeature, P);
		//显示计算的特征值
		for (j=0; j<numFeature; j++)
		{
			cout<<endl;
			for (k=0; k<numFeature; k++)
			{
				temp2 = 0;
				for (l=0; l<numFeature; l++)
				{
					temp2 += ((DOUBLE)Cova[j][l])*((DOUBLE)P[l][k]);
				}
				cout<<temp2/P[j][k]<<"  ";
			}
		}
		cout<<endl;
		DOUBLE* featureValue = new DOUBLE[numFeature];
		for (j=0; j<numFeature; j++)
		{
			featureValue[j] = 0;
		}
		for (k=0; k<numFeature; k++)//计算特征值
		{
			temp1 = 0;
			for (j=0; j<numFeature; j++)
			{
				temp2 = 0;
				for (l=0; l<numFeature; l++)
				{
					temp2 += ((DOUBLE)Cova[j][l])*((DOUBLE)P[l][k]);
				}
				if (P[j][k]>0)
				{
					temp1 += (DOUBLE)P[j][k];
					featureValue[k] += temp2;
				}
				else
				{
					temp1 -= (DOUBLE)P[j][k];
					featureValue[k] -= temp2;
				}
			}
			featureValue[k] = featureValue[k]/temp1;
			cout<<featureValue[k]<<endl;
		}
		maxv = 0;
		for (j=0; j<numFeature; j++)
		{
			if (featureValue[j] > maxv)
				maxv = featureValue[j];
		}
		int numLittle = 0;
		DOUBLE sunV = 0;
		for (j=0; j<numFeature; j++)
		{
			if (featureValue[j] < maxv/100)
			{
				numLittle ++;
				sunV += featureValue[j];
			}
		}
		for (j=0; j<numFeature; j++)
		{
			if (featureValue[j] < maxv/100)
			{
				featureValue[j] = sunV / ((DOUBLE)numLittle);
			}
		}
		for (j=0; j<numFeature; j++)
		{
			for (k=0; k<numFeature; k++)
			{
				Cova[j][k] = P[j][k]*((float)featureValue[k]);
			}
		}
		for (j=0; j<numFeature; j++)
		{
			for (k=0; k<numFeature; k++)
			{
				this->paramSigma[(i*numFeature+j)*numFeature+k] = 0;
				for (l=0; l<numFeature; l++)
				{
					this->paramSigma[(i*numFeature+j)*numFeature+k] += ((DOUBLE)Cova[j][l]*P[k][l]);
				}
			}
		}
	}

	//模型参数预处理
	this->paramInvSigma = new DOUBLE[numClass*numFeature*numFeature];
	this->paramLogSigma = new DOUBLE[numClass];
	CMatrix sigma;
	CMatrix invSigma;
	int size = numFeature * numFeature;
	int totalSample = this->sdata.numSample;
	sigma.nline = numFeature;
	sigma.ncol = numFeature;
	sigma.pdata = new DOUBLE[size];
	for (i=0; i<numClass; i++)
	{
		for (j=0; j<size; j++)
		{
			sigma.pdata[j] = this->paramSigma[i*size+j];
		}
		sigma.inv(invSigma);
		for (j=0; j<size; j++)
		{
			this->paramInvSigma[i*size+j] = invSigma.pdata[j];
		}
		this->paramLogSigma[i] = log(sigma.det())-2*log(((DOUBLE)classSample[i])/((DOUBLE)totalSample));//这里log表示自然对数
	}
	delete[] sumX;
	delete[] sumXX;
	delete[] classSample;
	return ret;
}

/* 使用分类器进行分类，输入测试样本的文件名称，返回true表示程序执行正常，false表示程序执行错误。 */
bool CMQdf::test(char* fileName)
{
	bool ret = true;
	ifstream ifs(fileName);
	int bufSize = 1024;
	char* buf=new char[bufSize];
	char* temp=new char[bufSize];
	int numClass = this->sdata.numClass;
	int numFeature = this->sdata.numFeature;
	DOUBLE* sample=new DOUBLE[numFeature];
	DOUBLE* value=new DOUBLE[numClass];
	int n=0;
	int i, j, k;
	int cnum;
	DOUBLE max;
	int index;

	if (NULL==this->resultTotal)
		this->resultTotal = new DOUBLE[numClass];
	if (NULL==this->resultRight)
		this->resultRight = new DOUBLE[numClass];
	for (i=0; i<numClass; i++)
	{
		this->resultTotal[i] = 0;
		this->resultRight[i] = 0;
	}
	while (ifs.good()) {
		ifs.getline(buf, bufSize);
		n=0;
		while (*(buf+n)!=',' && n<bufSize)
		{
			*(temp+n)=*(buf+n);
			n++;
		}
		*(temp+n)=0;
		cnum=this->sdata.searchClassName(temp);
		if (cnum>-1)
		{
			this->resultTotal[cnum] += 1;
			for (j=0; j<numFeature; j++)
			{
				n++;
				int n2=n;
				while (*(buf+n)!=',' && n<bufSize)
				{
					*(temp+n-n2)=*(buf+n);
					n++;
				}
				*(temp+n-n2)=0;
				*(sample+j)=atof(temp);
				if (this->sdata.maxValue[j] > this->sdata.minValue[j])
				{
					*(sample+j) = (*(sample+j) - this->sdata.minValue[j])/(this->sdata.maxValue[j] - this->sdata.minValue[j]);
				}
				else
					*(sample+j) = 0.5;
			}
			for (i=0; i<numClass; i++)
			{
				*(value+i) = -this->paramLogSigma[i];
				for (j=0; j<numFeature; j++)
				{
					for (k=0; k<numFeature; k++)
					{
						*(value+i) -= (*(sample+j) - this->paramU[i*numFeature+j])*(this->paramInvSigma[(i*numFeature+j)*numFeature+k])*(*(sample+k) - this->paramU[i*numFeature+k]);
					}
				}
			}
			max = *value;
			index = 0;
			for (i=1; i<numClass; i++)
			{
				if (*(value+i)>max)
				{
					max = *(value+i);
					index = i;
				}
			}
			if (index == cnum)
				this->resultRight[cnum] += 1;
		}
	}
	delete[] buf;
	delete[] temp;
	delete[] sample;
	delete[] value;
	return ret;
}