matrix.cpp

模式识别的作业代码

	for(i=0;i<KP;i++) //单位化
	{
		tmp = -A.GetElement(i,IP)*CQ;
		B.SetElement(i,IP,tmp);
	}

	for(i=0;i<KP;i++)
	{
		tmp = A.GetElement(IP,i)*CQ; //单位化
		B.SetElement(IP,i,tmp);
	}
	B.SetElement(IP,IP,CQ);

	A = B;
}

////////////////////////////
// 对当前矩阵求逆:
// 矩阵求逆 B = inv( /*this */ )  (高斯消元法)
// 输出参数 B --this的逆
CMatrix CMatrix::inv() const
{
	int nDim;
	int i,j;
	CMatrix invA(*this);
	
	if((nDim=this->GetNumRow()) != this->GetNumColumns())
	{
		TRACE("不是方阵");
		return invA;
	}

	int KP = nDim+1;
	CMatrix G(KP); //创建一个nDim+1维的方阵

	for(i=0;i<nDim;i++) //
	for(j=0;j<nDim;j++)
	{
		G.SetElement(i,j,this->GetElement(i,j));
	}
	G.SetElement(KP,KP,1.0e5);

	for(i=0;i<nDim;i++)
		WGAU(i,G);

	for(i=0;i<nDim;i++)
	for(j=0;j<nDim;j++)
		invA.SetElement(i,j,G.GetElement(i,j));

	return invA;
}

int CMatrix::GetRowVector(int nRow,double* pVector) const
{
	if(pVector == NULL) return -1;
	int COL = this->m_nNumColumns ;
	for(int i=0;i< COL;i++)
	{
		ASSERT(&pVector[i]);
		pVector[i] = this->GetElement(nRow,i);
	}
	return nRow;
}

int CMatrix::GetColVector(int nCol,double* pVector) const
{
	if(pVector == NULL) return -1;
	int ROW = this->m_nNumRows ;
	for(int i=0;i< ROW;i++)
	{
		ASSERT(&pVector[i]);
		pVector[i] = this->GetElement(i,nCol);
	}
	return nCol;
}

int CMatrix::GetSize() const
{
	return (this->m_nNumColumns*this->m_nNumRows);
}

CMatrix CMatrix::mean(int dim) const //求矩阵的均值
{
	CMatrix clsMean;
	int nSize = this->GetSize();
	double dMean = 0;
	int i,j;

	if(m_nNumColumns == 1 || m_nNumRows == 1) //为一维数据
	{
		clsMean.Init(1,1); //只有一个值;

		dMean = 0;
		for(i=0;i<nSize;i++)
			dMean += m_pData[i];
		dMean /= nSize;
		clsMean.SetElement(0,0,dMean);
	}
	else //为矩阵
	{
		if(dim == 2) //对每行向量求均值，返回一个均值列向量 n*1
		{
			clsMean.Init(m_nNumRows,1);
			for(i=0;i<m_nNumRows;i++)
			{
				dMean = 0;
				for(j=0;j<m_nNumColumns;j++)
				{
					dMean += this->GetElement(i,j);
				}
				dMean /= m_nNumColumns;
				clsMean.SetElement(i,0,dMean);
			}
		}
		else if(dim == 3) //对整个矩阵求均值 ,返回 1*1矩阵
		{
			clsMean.Init(1,1);

			dMean = 0;
			for(i=0;i<m_nNumRows;i++)
			{
				for(j=0;j<m_nNumColumns;j++)
				{
					dMean += this->GetElement(i,j);
				}
			}
			dMean /= nSize;
			clsMean.SetElement(0,0,dMean);
		}
		else ////对每列向量求均值，返回一个均值行向量 1*n
		{
			clsMean.Init(1,m_nNumColumns);
			for(i=0;i<m_nNumColumns;i++)
			{
				dMean = 0;
				for(j=0;j<m_nNumRows;j++)
				{
					dMean += this->GetElement(j,i);
				}
				dMean /= m_nNumRows;
				clsMean.SetElement(0,i,dMean);
			}
		}
	}
	return clsMean;
}

CMatrix CMatrix::std(int dim) const //求矩阵的方差
{
	CMatrix clsStd;

	CMatrix clsMean;
	clsMean = this->mean(dim); //先求均值

	int nSize = this->GetSize();
	double dStd = 0 , dMean;
	int i,j;

	if(m_nNumColumns == 1 || m_nNumRows == 1) //为一维数据
	{
		clsStd.Init(1,1); //只有一个值;
		dMean = clsMean.GetElement(0,0);

		dStd = 0;
		for(i=0;i<nSize;i++)
			dStd += (m_pData[i] - dMean)*(m_pData[i] - dMean);
		dStd = sqrt(dStd/(nSize-1));

		clsStd.SetElement(0,0,dStd);
	}
	else //为矩阵
	{
		if(dim == 2) //对每行向量求均值，返回一个均值列向量 n*1
		{
			clsStd.Init(m_nNumRows,1);
			for(i=0;i<m_nNumRows;i++)
			{
				dMean = clsMean.GetElement(i,0);

				dStd = 0;
				for(j=0;j<m_nNumColumns;j++)
				{
					dStd += (GetElement(i,j) - dMean)*(GetElement(i,j) - dMean);
				}
				dStd = sqrt(dStd/(m_nNumColumns-1));
				clsStd.SetElement(i,0,dStd);
			}
		}
		else if(dim == 3) //对整个矩阵求均值 ,返回 1*1矩阵
		{
			clsStd.Init(1,1);
			dMean = clsMean.GetElement(0,0);

			dStd = 0;
			for(i=0;i<m_nNumRows;i++)
			{
				for(j=0;j<m_nNumColumns;j++)
				{
					dStd += (GetElement(i,j) - dMean)*(GetElement(i,j) - dMean);
				}
			}
			dStd =sqrt(dStd/(nSize-1));
			clsStd.SetElement(0,0,dStd);
		}
		else ////对每列向量求均值，返回一个均值行向量 1*n
		{
			clsStd.Init(1,m_nNumColumns);
			for(i=0;i<m_nNumColumns;i++)
			{
				dStd = 0;
				dMean = clsMean.GetElement(0,i);
				for(j=0;j<m_nNumRows;j++)
				{
					dStd += (GetElement(j,i) - dMean)*(GetElement(j,i) - dMean);
				}
				dStd =sqrt(dStd/(m_nNumRows-1));
				clsStd.SetElement(0,i,dStd);
			}
		}
	}
	return clsStd;
}

void CMatrix::multipVec(double *dMultipVec, double *dFinalVec)
{
	/*
	
		这个程序的作用是this矩阵乘上一个向量dMultipVec（长为矩阵列数），
		获得结果向量dFinalVec（长为矩阵行数）的意思。

		gived by ZWC 
		2003 年 10 月 6 日
	
	*/
	CMatrix  dMatrixTem(*this);
	double	 dTem;
	int		 i,j;

	for(i=0;i<m_nNumRows;i++)
	{
		dTem=0;
		for(j=0;j<m_nNumColumns;j++)
		{
			dTem+=m_pData[m_nNumColumns*i+j]*dMultipVec[j];
		}
		dFinalVec[i]=dTem;
	}
		

}




//////////////////////////////////////////////////////////////////////
// 求实对称矩阵特征值与特征向量的雅可比过关法
//
// 参数：
// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
//    数组dblEigenValue中第j个特征值对应的特征向量
// 3. double eps - 计算精度，默认值为0.000001
//
// 返回值：BOOL型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::JacobiEigenv2(double dblEigenValue[], CMatrix& mtxEigenVector, double eps /*= 0.000001*/)
{ 
	int i,j,p,q,u,w,t,s;
    double ff,fm,cn,sn,omega,x,y,d;
    
	if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
		return FALSE;

	for (i=0; i<=m_nNumColumns-1; i++)
    { 
		mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;
        for (j=0; j<=m_nNumColumns-1; j++)
			if (i!=j) 
				mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;
    }
    
	ff=0.0;
    for (i=1; i<=m_nNumColumns-1; i++)
	{
		for (j=0; j<=i-1; j++)
		{ 
			d=m_pData[i*m_nNumColumns+j]; 
			ff=ff+d*d; 
		}
	}

    ff=sqrt(2.0*ff);

Loop_0:
    
	ff=ff/(1.0*m_nNumColumns);

Loop_1:

    for (i=1; i<=m_nNumColumns-1; i++)
	{
		for (j=0; j<=i-1; j++)
        { 
			d=fabs(m_pData[i*m_nNumColumns+j]);
            if (d>ff)
            { 
				p=i; 
				q=j;
                goto Loop_2;
            }
        }
	}
        
	if (ff<eps) 
	{
		for (i=0; i<m_nNumColumns; ++i)
			dblEigenValue[i] = GetElement(i,i);
		return TRUE;
	}
    
	goto Loop_0;

Loop_2: 
		
	u=p*m_nNumColumns+q; 
	w=p*m_nNumColumns+p; 
	t=q*m_nNumColumns+p; 
	s=q*m_nNumColumns+q;
    x=-m_pData[u]; 
	y=(m_pData[s]-m_pData[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=m_pData[w];
    m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData[u]*omega;
    m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega;
    m_pData[u]=0.0; m_pData[t]=0.0;
    
	for (j=0; j<=m_nNumColumns-1; j++)
	{
		if ((j!=p)&&(j!=q))
		{ 
			u=p*m_nNumColumns+j; 
			w=q*m_nNumColumns+j;
			fm=m_pData[u];
			m_pData[u]=fm*cn+m_pData[w]*sn;
			m_pData[w]=-fm*sn+m_pData[w]*cn;
		}
	}

    for (i=0; i<=m_nNumColumns-1; i++)
    {
		if ((i!=p)&&(i!=q))
        { 
			u=i*m_nNumColumns+p; 
			w=i*m_nNumColumns+q;
			fm=m_pData[u];
			m_pData[u]=fm*cn+m_pData[w]*sn;
			m_pData[w]=-fm*sn+m_pData[w]*cn;
        }
	}
    
	for (i=0; i<=m_nNumColumns-1; i++)
    { 
		u=i*m_nNumColumns+p; 
		w=i*m_nNumColumns+q;
        fm=mtxEigenVector.m_pData[u];
        mtxEigenVector.m_pData[u]=fm*cn+mtxEigenVector.m_pData[w]*sn;
        mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn;
	}

	goto Loop_1;
}



BOOL CMatrix::SymFullRankOrNot()
{
//JacobiEigenv2(double dblEigenValue[], CMatrix& mtxEigenVector, double eps /*= 0.000001*/)

	BOOL    bTem;

	double*  dblEigenValue;

	double dTem0,dTem1;

	dblEigenValue=new double[m_nNumRows];

	CMatrix mtxEigenVector(m_nNumRows,m_nNumRows);

	//行列不等
	if(m_nNumColumns!=m_nNumRows)
		return FALSE;

	bTem=JacobiEigenv2(dblEigenValue,mtxEigenVector);

	//求特征值失败
	if(! bTem)
		return FALSE;

	dTem1=fabs(dblEigenValue[m_nNumRows-1] );

	//具有为零的特征值

	if(dTem1<1.0e-10 )
		return FALSE;

	dTem0=fabs(dblEigenValue[0] );

	//最大最小特征值相差太大
	if((dTem0/dTem1)>1.0e+6)
		return FALSE;

	delete []dblEigenValue;

	return TRUE;
}