matrix.cpp

变分独立因子分析C++代码,H.Attias说比独立分量分析要好,但是这个程序分析效果不好,可能是程序问题,也可能是对理论理解不透

// Matrix.cpp: implementation of the CMatrix class.
// Shijun Sun    Tutor: Peng Chenglin
// Chongqing University, Bioengineering College   2007/10/18
//////////////////////////////////////////////////////////////////////

#include "Matrix.h"
#include <memory.h>
#include <math.h>
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CMatrix::CMatrix()
{
	m_nr=m_nc=0;
	m_matrix=0;
}

CMatrix::~CMatrix()
{
	clear();
}

void CMatrix::clear()
{
	if(m_matrix!=0)
	{
		for(int i=0;i<m_nr;i++) delete[] m_matrix[i];
		delete[] m_matrix;
		m_nr=m_nc=0;
		m_matrix=0;
	}
}

void CMatrix::create(const int nr, const int nc)
{
	m_nr=nr;
	m_nc=nc;
	m_matrix=new double*[nr];
	for(int i=0;i<nr;i++) m_matrix[i]=new double[nc];
}

CMatrix::CMatrix(const int nr, const int nc)
{
	create(nr,nc);
}

CMatrix::CMatrix(double **matrix, const int nr, const int nc)
{	
	create(nr,nc);
	for(int i=0;i<nr;i++)
	{
		memcpy(m_matrix[i], matrix[i], nc*sizeof(double));
	}
}

CMatrix::CMatrix(const CMatrix &matrix)
{
	create( matrix.m_nr,matrix.m_nc );
	for(int i=0;i<m_nr;i++)
	{
		memcpy(m_matrix[i],matrix[i],m_nc*sizeof(double));
	}
}

CMatrix & CMatrix::operator =(const CMatrix &matrix)
{	
    if(&matrix==this) return *this;  
	
	clear();
	create(matrix.m_nr,matrix.m_nc);
	for(int i=0;i<m_nr;i++) 
	{
		memcpy(m_matrix[i], matrix.m_matrix[i],m_nc*sizeof(double) );
	}
	return *this;
}

CMatrix CMatrix::operator -() const
{
	CMatrix m=*this;
	for(int i=0;i<m_nr;i++) 
		for(int j=0;j<m_nc;j++) 
			m[i][j]=-m[i][j];
	return m;
}

CMatrix CMatrix::operator +() const
{
	return CMatrix(*this);
}

CMatrix CMatrix::operator ~() const
{
	CMatrix m(m_nc,m_nr);
	for(int i=0;i<m_nc;i++)
	{
		for(int j=0;j<m_nr;j++)
		{
			m[i][j]=m_matrix[j][i];
		}
	}
	return m;
}

CMatrix CMatrix::operator +(const CMatrix &matrix) 
{
	CMatrix m=*this;
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++) 
		{
			m[i][j]+=matrix[i][j];
		}
	}	
	return m;
}

CMatrix CMatrix::operator -(const CMatrix &matrix)
{
	CMatrix m=*this;
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++) 
		{
			m[i][j]-=matrix[i][j];
		}
	}
	return m;
}

CMatrix CMatrix::operator *(const CMatrix &matrix)
{
	CMatrix m(m_nr,matrix.m_nc);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<matrix.m_nc;j++)
		{
			m[i][j]=0;
			for(int k=0;k<m_nc;k++)
			{
				m[i][j]+=m_matrix[i][k]*matrix[k][j];
			}
		}
	}
	return m;
}

CMatrix CMatrix::operator *(const double x)
{
	CMatrix m=*this;
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++) m[i][j]*=x;
	}
	return m;
}

CMatrix CMatrix::operator /(const double x)
{	
	CMatrix m=*this;
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++) m[i][j]/=x;
	}
	return m; //return *this*(1/x);
}

////////////////////////////////////////////////////////////////
// friend function and operator

double scalarvalue(const CMatrix &matrix)
{
	return matrix[0][0];
}

CMatrix operator *(const double x, const CMatrix & matrix)
{
	CMatrix m=matrix;
	for(int i=0;i<m.getnr();i++)
		for(int j=0;j<m.getnc();j++) 
			m[i][j]=x*matrix[i][j];
	return m;
}


double diagdeterminant(const CMatrix & matrix)
{
	double det=1;
	for(int i=0;i<matrix.m_nr;i++)
	{
		det*=matrix[i][i];
	}
	return det;
}

CMatrix diaginverse(const CMatrix & matrix)
{
	CMatrix m=matrix;
	for(int i=0;i<m.m_nr;i++)
		m[i][i]=1/matrix[i][i];
	return m;
}

CMatrix symmlowtridecompose(const CMatrix & matrix)
{
	int n=matrix.m_nr;
	CMatrix m(n,n);
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<i;j++)
		{
			m[i][j]=matrix[i][j];
			for(int k=0;k<j;k++)
			{
				m[i][j]-=m[i][k]*m[j][k];
			}
			m[i][j]/=m[j][j];
		}

		m[i][i]=matrix[i][i];
		for(int k=0;k<i;k++)
		{
			m[i][i]-=m[i][k]*m[i][k];
		}
		m[i][i]=sqrt(m[i][i]);

		for(k=i+1;k<n;k++)
		{
			m[i][k]=0;
		}
	}
	return m;
}

double symmdeterminant(const CMatrix & matrix)
{
	double det=lowtrideterminat(symmlowtridecompose(matrix));
	return det*det;
}

double lowtrideterminat(const CMatrix & matrix)
{
	return diagdeterminant(matrix);
}

CMatrix lowtriinverse(const CMatrix & matrix)
{
	int n=matrix.m_nr;
	CMatrix m(n,n);
	for(int i=0;i<n;i++)
	{	
		m[i][i]=1/matrix[i][i];
		for(int j=0;j<i;j++)
		{	
			m[i][j]=0;
			for(int k=j;k<i;k++) m[i][j]-=matrix[i][k]*m[k][j];
			m[i][j]/=matrix[i][i];
		}
	}
	for(i=0;i<n;i++)
		for(int j=0;j<i;j++)
			m[j][i]=0;
	return m;
}

CMatrix symminverse(const CMatrix & matrix)
{
	CMatrix low=symmlowtridecompose(matrix);
	return ~lowtriinverse(low)*lowtriinverse(low);
}

CMatrix linearsimultaneousequations(const CMatrix & A, const CMatrix & B)
{
	CMatrix a=A;
	CMatrix b=B;
	int n=a.m_nr;
	CMatrix x(n,1);
	for(int i=0;i<n-1;i++)// ok
	{
		int pos=i;
		double temp=a[i][i];
		for(int j=i+1;j<n;j++)
		{
			if(temp<a[j][i])
			{
				temp=a[j][i];
				pos=j;
			}
		}
		if(i!=pos)
		{
			for(int k=i;k<n;k++)
			{
				swap(a[i][k],a[pos][k]);
			}
			swap(b[i][0],b[pos][0]);
		}

		for(j=i+1;j<n;j++)
		{
			double temp=a[j][i]/a[i][i];
			for(int k=i+1;k<n;k++)
			{
				a[j][k]-=temp*a[i][k];
			}
			b[j][0]-=temp*b[i][0];
		}
	}


	for(i=n-1;i>=0;i--)
	{
		x[i][0]=b[i][0];
		for(int j=i+1;j<n;j++)
		{
			x[i][0]-=a[i][j]*x[j][0];
		}
		x[i][0]/=a[i][i];
	}
	return x;
}

void swap(double &x, double &y)
{
	double temp=x;
	x=y;
	y=temp;
}