gmatrix.cpp

一个由Mike Gashler完成的机器学习方面的includes neural net, naive bayesian classifier, decision tree, KNN, a genet

	{
		dVal = m_pData[m_nColumns * nRow + nRow];
		for(nCol = 0; nCol < m_nColumns; nCol++)
			m_pData[m_nColumns * nRow + nCol] /= dVal;
		pVector[nRow] /= dVal;
	}
}

int GMatrix::CountNonZeroElements()
{
	int nCount = 0;
	int nRow, nCol;
	for(nRow = 0; nRow < m_nRows; nRow++)
	{
		for(nCol = 0; nCol < m_nColumns; nCol++)
		{
			if(Get(nRow, nCol) != 0)
				nCount++;
		}
	}
	return nCount;
}

bool GMatrix::Cholesky(const GMatrix* pIn)
{
	Resize(pIn->m_nRows, pIn->m_nColumns);
	int i, j, k;
	double d;
	for(j = 0; j < m_nRows; j++)
	{
		for(i = 0; i < j; i++)
		{
			d = 0;
			for(k = 0; k < i; k++)
				d += (Get(i, k) * Get(j, k));
			Set(j, i, (1.0 / Get(i, i)) * (pIn->Get(i, j) - d));
		}
		d = 0;
		for(k = 0; k < i; k++)
			d += (Get(i, k) * Get(j, k));
		d = pIn->Get(j, i) - d;
		if(d < 0)
		{
			if(d > -1e-12)
				d = 0; // it's probably just round error
			else
				return false; // pIn is not positive definite
		}
		Set(j, i, sqrt(d));
		for(i++; i < m_nColumns; i++)
			Set(j, i, 0);
	}
	return true;
}

void GMatrix::ComputeEigenVectors(int nCount, const GMatrix* pInputMatrix)
{
	// Check values and resize this matrix
	int nDims = pInputMatrix->GetColumnCount();
	GAssert(nDims == pInputMatrix->GetRowCount(), "not a covariance matrix");
	GAssert(nCount <= nDims, "Can't have more eigenvectors than columns");
	Resize(nCount, nDims);
	GArffData data(nDims * 2);

	// Compute the deviation matrix. (Deviation is just the square root of variance. Covariance
	// is the multi-dimensional equivalent of variance, and eigenvectors have meaning in respect to a
	// covariance matrix, so we simply take the square root of the input matrix.)
	GMatrix deviation;
	deviation.Cholesky(pInputMatrix);
	deviation.Transpose();
	GMatrix negdev;
	negdev.Copy(&deviation);
	negdev.Multiply(-1);

	// Generate representative data. (If your covariance matrix was computed from
	// original data then obviously it's a waste of time to generating new data
	// to represent your original data. In that case you should have just called
	// GArffData::ExtractPrincipleComponent directly on the original data and never
	// bothered to compute the covariance matrix.)
	int i;
	for(i = 0; i < nDims; i++)
	{
		data.AddVector(negdev.GetRow(i));
		data.AddVector(deviation.GetRow(i));
	}

	// Extract the principle components
	for(i = 0; i < nCount; i++)
		data.ComputePrincipleComponent(nDims, GetRow(i), 10, true);

	data.DropAllVectors();
}

void GMatrix::ComputeEigenValues(double* pOutValues, GMatrix* pEigenVectors)
{
	GTEMPBUF(double, pVec, pEigenVectors->GetColumnCount());
	double* pEigVec;
	int row, col;
	for(row = 0; row < pEigenVectors->GetRowCount(); row++)
	{
		GMatrix tmp;
		tmp.Copy(this);
		pEigVec = pEigenVectors->GetRow(row);
		tmp.Multiply(pEigVec, pVec);
		pOutValues[row] = 0;
		for(col = 0; col < pEigenVectors->GetColumnCount(); col++)
			pOutValues[row] += pVec[col] / pEigVec[col];
		pOutValues[row] /= pEigenVectors->GetRowCount();
	}
}

/*
#ifdef ARPACK

#	include <arpack++/armat.h>
#	include <arpack++/arcomp.h>
#	include <arpack++/arlsmat.h>
#	include <arpack++/arlnsmat.h>
#	include <arpack++/arlssym.h>
#	include <arpack++/arlgsym.h>
#	include <arpack++/arlgnsym.h>
#	include <arpack++/arlscomp.h>
#	include <arpack++/arlgcomp.h>

//#	include <arpack++/arbsnsym.h>
//#	include <arpack++/ardsnsym.h>
#	include <arpack++/arlsnsym.h>
//#		include <arpack++/arusnsym.h>

void ArpackError::Set(ArpackError::ErrorCode eCode, char const* szMessage)
{
	fprintf(stderr, "Got an error: [%d] %s\n", eCode, szMessage);
}

void MemoryOverflow()
{
	fprintf(stderr, "memory overflow\n");
}

typedef double ARFLOAT;

void GMatrix::ComputeEigenVectors(GVector* pOutEigenValues, GMatrix* pOutEigenVectors)
{
	// Allocate space for a sparse matrix
	printf("allocating space for the matrix...\n");
	int nRows = m_nRows;
	GAssert(m_nRows == m_nColumns, "expected a square matrix");
	int nNonZeroElements = CountNonZeroElements();
	int* pRowIndexes = new int[nNonZeroElements];
	int* pColIndexes = new int[nRows + 1];
	ARFLOAT* pMatrix = new ARFLOAT[nNonZeroElements];

	// Make the matrix
	printf("making the matrix...\n");
	int col, row;
	int j = 0; // matrix position
	pColIndexes[0] = 0;
	for(col = 0; col < m_nColumns; col++)
	{
		for(row = 0; row < m_nRows; row++)
		{
			if(Get(row, col) != 0)
			{
				pRowIndexes[j] = row;
				pMatrix[j++] = Get(row, col);
			}
		}
		pColIndexes[col + 1] = j;
	}

	// Allocate space for the output
	int nMaxEigenVals = 100;
	printf("allocating space for the output...\n");
	ARFLOAT* pEigValReal = new ARFLOAT[nRows];
	ARFLOAT* pEigValImag = new ARFLOAT[nRows];
	ARFLOAT* pEigVec = new ARFLOAT[nMaxEigenVals * nRows];

	// Create an ARPACK++ matrix
	printf("creating an arpack++ matrix...\n");
	ARluNonSymMatrix<ARFLOAT, ARFLOAT> matrix(nRows, nNonZeroElements, pMatrix, pRowIndexes, pColIndexes);

	// Construct the eigenvalue solver
	printf("constructing eigen solver...\n");
	ARluNonSymStdEig<ARFLOAT> solver(nMaxEigenVals, matrix, "SR", 0, 0.00001, 500000, NULL, true);

	// Find the eigenvalues and vectors
	printf("solving for eigen values and vectors...\n");
	int nEigenValues = solver.EigenValVectors(pEigVec, pEigValReal, pEigValImag);
	printf("Number of eigen values computed: %d\n", nEigenValues);

	// Copy the data into the output space
	pOutEigenVectors->Resize(m_nRows, m_nColumns);
	pOutEigenValues->Resize(m_nRows);
	for(j = 0; j < nEigenValues; j++)
	{
		row = m_nRows - 1 - j;
		for(col = 0; col < m_nColumns; col++)
			pOutEigenVectors->Set(row, col, pEigVec[m_nRows * j + col]);
		pOutEigenValues->Set(row, pEigValReal[j]);
	}

	// Clean up
	delete(pEigVec);
	delete(pEigValImag);
	delete(pEigValReal);
}
#else // ARPACK
#	ifdef OCTAVE

#		define EIGHT_BYTE_INT long long int
#		define GCC_ATTR_NORETURN
#		define SIZEOF_SHORT 2
#		define SIZEOF_INT 4
//#		define SIZEOF_LONG 4
#		include <octave/octave/Matrix.h>
#		include <octave/octave/EIG.h>

void GMatrix::ComputeEigenVectors(GVector* pOutEigenValues, GMatrix* pOutEigenVectors)
{
	Matrix mTmp(m_nRows, m_nColumns);
	int row, col;
	for(row = 0; row < m_nRows; row++)
	{
		for(col = 0; col < m_nColumns; col++)
			mTmp.elem(row, col) = Get(row, col);
	}
	EIG eig(mTmp);
	ComplexColumnVector eigenvalues = eig.eigenvalues();
	ComplexMatrix eigenvectors = eig.eigenvectors();
	pOutEigenVectors->Resize(m_nRows, m_nColumns);
	pOutEigenValues->Resize(m_nRows);
	for(row = 0; row < m_nRows; row++)
	{
		for(col = 0; col < m_nColumns; col++)
			pOutEigenVectors->Set(row, col, eigenvectors.elem(col, row).real());
		pOutEigenValues->Set(row, (double)eigenvalues.elem(row).real());
	}
}

#	else // OCTAVE

void GMatrix::ComputeEigenVectors(GVector* pOutEigenValues, GMatrix* pOutEigenVectors)
{
	const char* szMessage = "No method is available for computing eigen vectors. (You need to uncomment one of the #defines at the top of this file and link to the corresponding library.)";
	GAssert(false, szMessage);
	throw(szMessage);
}

#	endif // !OCTAVE
#endif // !ARPACK
*/
double GMatrix::ComputeSumSquaredDifference(GMatrix* other)
{
	GAssert(m_nRows == other->m_nRows && m_nColumns == other->m_nColumns, "mismatching sizes");
	int i, j;
	double dError = 0;
	double d;
	for(j = 0; j < m_nRows; j++)
	{
		for(i = 0; i < m_nRows; i++)
		{
			d = Get(j, i) - other->Get(j, i);
			dError += (d * d);
		}	
	}
	return dError;
}

#ifndef NO_TEST_CODE
/*static*/ void GMatrix::Test()
{
	// Test matrix multiplication
	GMatrix m1(5, 5);
	m1.Set(0, 0, 3);	m1.Set(0, 1, 2);	m1.Set(0, 2, 7);	m1.Set(0, 3, 2);	m1.Set(0, 4, 1);
	m1.Set(1, 0, 4);	m1.Set(1, 1, 6);	m1.Set(1, 2, 2);	m1.Set(1, 3, 6);	m1.Set(1, 4, 6);
	m1.Set(2, 0, 0);	m1.Set(2, 1, 2);	m1.Set(2, 2, 3);	m1.Set(2, 3, 2);	m1.Set(2, 4, 9);
	m1.Set(3, 0, 0);	m1.Set(3, 1, 3);	m1.Set(3, 2, 3);	m1.Set(3, 3, 3);	m1.Set(3, 4, 9);
	m1.Set(4, 0, 0);	m1.Set(4, 1, 6);	m1.Set(4, 2, 3);	m1.Set(4, 3, 2);	m1.Set(4, 4, 5);
	GMatrix m2(3, 3);
	m2.Copy(&m1);
	m2.Set(2, 2, 8.1);
	m2.Set(2, 3, 8.8);
	GMatrix m3(8, 8);
	m3.Multiply(&m1, &m2);
	if(m3.Get(4, 3) < 78.39999999 || m3.Get(4, 3) > 78.400000001)
		throw "matrix multiplication return the wrong answer";

	// Test Cholesky decomposition
	m1.Resize(3, 3);
	m1.Set(0, 0, 3);	m1.Set(0, 1, 0);	m1.Set(0, 2, 0);
	m1.Set(1, 0, 1);	m1.Set(1, 1, 4);	m1.Set(1, 2, 0);
	m1.Set(2, 0, 2);	m1.Set(2, 1, 2);	m1.Set(2, 2, 7);
	m2.Copy(&m1);
	m2.Transpose();
	m3.Multiply(&m1, &m2);
	GMatrix m4;
	if(!m4.Cholesky(&m3))
		throw "huh? not positive definate?";
	if(m1.ComputeSumSquaredDifference(&m4) >= .0001)
		throw "Cholesky decomposition didn't work right";

	// Test eigenvector computation
	m1.Resize(2, 2);
	m1.Set(0, 0, 1);	m1.Set(0, 1, 1);
	m1.Set(1, 0, 1);	m1.Set(1, 1, 4);
	m2.ComputeEigenVectors(2, &m1);
	if(ABS(m2.Get(0, 0) * m2.Get(0, 0) + m2.Get(0, 1) * m2.Get(0, 1) - 1) > .0001)
		throw "answer not normalized";
	if(ABS(m2.Get(0, 0) * m2.Get(0, 1) - .27735) > .0001)
		throw "wrong answer";
	if(ABS(m2.Get(1, 0) * m2.Get(1, 0) + m2.Get(1, 1) * m2.Get(1, 1) - 1) > .0001)
		throw "answer not normalized";
	if(ABS(m2.Get(1, 0) * m2.Get(1, 1) + .27735) > .0001)
		throw "wrong answer";
}
#endif // !NO_TEST_CODE