gmatrix.cpp

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

/*
	Copyright (C) 2006, Mike Gashler

	This library is free software; you can redistribute it and/or
	modify it under the terms of the GNU Lesser General Public
	License as published by the Free Software Foundation; either
	version 2.1 of the License, or (at your option) any later version.

	see http://www.gnu.org/copyleft/lesser.html
*/

#include "GMatrix.h"
#include "GMacros.h"
#include "GVector.h"
#include "GBits.h"
#include "GArff.h"
#include <math.h>

// These tell which library to use for computing eigen vectors. You should uncomment only one of them. Octave is known to give good results, but it's GPL'ed. I got ARPACK working once, but I don't know if it's still in working shape.
//#define OCTAVE
//#define ARPACK


GMatrix::GMatrix(int nRows, int nColumns)
{
	m_nRows = nRows;
	m_nColumns = nColumns;
	int nSize = nRows * nColumns;
	if(nSize > 0)
	{
		m_pData = new double[nSize];
		int n;
		for(n = 0; n < nSize; n++)
			m_pData[n] = 0;
	}
	else
		m_pData = NULL;
}

GMatrix::~GMatrix()
{
	delete(m_pData);
}

void GMatrix::SetToIdentity()
{
	int nRow, nColumn;
	for(nRow = 0; nRow < m_nRows; nRow++)
	{
		for(nColumn = 0; nColumn < m_nColumns; nColumn++)
			Set(nRow, nColumn, 0);
	}
	int nMin = MIN(m_nColumns, m_nRows);
	for(nColumn = 0; nColumn < nMin; nColumn++)
		Set(nColumn, nColumn, 1);
}

void GMatrix::Resize(int nRows, int nColumns)
{
	if(m_nRows == nRows && m_nColumns == nColumns)
		return;
	delete(m_pData);
	m_pData = new double[nRows * nColumns];
	m_nRows = nRows;
	m_nColumns = nColumns;
}

void GMatrix::Copy(const GMatrix* pMatrix)
{
	Resize(pMatrix->m_nRows, pMatrix->m_nColumns);
	memcpy(m_pData, pMatrix->m_pData, m_nRows * m_nColumns * sizeof(double));
}

void GMatrix::Multiply(double dScalar)
{
	int nElements = m_nColumns * m_nRows;
	int i;
	for(i = 0; i < nElements; i++)
		m_pData[i] *= dScalar;
}

void GMatrix::Multiply(const double* pVectorIn, double* pVectorOut)
{
	int i;
	for(i = 0; i < m_nRows; i++)
		pVectorOut[i] = GVector::ComputeDotProduct(GetRow(i), pVectorIn, m_nColumns);
}

void GMatrix::Multiply(const GMatrix* pA, const GMatrix* pB)
{
	if(pA->m_nColumns != pB->m_nRows)
	{
		GAssert(false, "Incompatible sizes--can't multiply");
		return;
	}
	Resize(pA->m_nRows, pB->m_nColumns);
	double dSum;
	int nPos;
	int nRow;
	int nColumn;
	for(nRow = 0; nRow < pA->m_nRows; nRow++)
	{
		for(nColumn = 0; nColumn < pB->m_nColumns; nColumn++)
		{
			dSum = 0;
			for(nPos = 0; nPos < pA->m_nColumns; nPos++)
				dSum += pA->Get(nRow, nPos) * pB->Get(nPos, nColumn);
			Set(nRow, nColumn, dSum);
		}
	}
}

void GMatrix::GetColumn(int nCol, double* pVector)
{
	int i;
	for(i = 0; i < m_nRows; i++)
	{
		pVector[i] = m_pData[nCol];
		nCol += m_nColumns;
	}
}

double GMatrix::ComputeColumnMean(int nCol) const
{
	return ComputeColumnSum(nCol) / m_nRows;
}

double GMatrix::ComputeColumnSum(int nCol) const
{
	double dSum = 0;
	int i;
	for(i = 0; i < m_nRows; i++)
		dSum += Get(i, nCol);
	return dSum;
}

/*
// todo: this algorithm isn't very numerically stable because the
// values will tend to get too large or too small and lose data
double GMatrix::GetDeterminant()
{
	if(m_nRows != m_nColumns)
	{
		GAssert(false, "expected a square matrix");
		return 0;
	}
	double dDet = 0;
	double d;
	int col, i;
	for(col = 0; col < m_nColumns; col++)
	{
		d = 1;
		for(i = 0; i < m_nRows; i++)
			d *= Get(i, (col + i) % m_nColumns);
		dDet += d;
		d = 1;
		for(i = 0; i < m_nRows; i++)
			d *= Get(i, (col - i + m_nColumns) % m_nColumns);
		dDet -= d;
	}
	return dDet;
}
*/
void GMatrix::Print()
{
	int nRow;
	int nColumn;
	for(nRow = 0; nRow < m_nRows; nRow++)
	{
		printf("%f", Get(nRow, 0));
		for(nColumn = 1; nColumn < m_nColumns; nColumn++)
			printf("\t%f", Get(nRow, nColumn));
		printf("\n");
	}
}

void GMatrix::PrintCorners(int n)
{
	int nRow, nCol;
	for(nRow = 0; nRow < n && nRow < m_nRows; nRow++)
	{
		printf("%f", Get(nRow, 0));
		for(nCol = 1; nCol < n && nCol < m_nColumns; nCol++)
			printf("\t%f", Get(nRow, nCol));
		if(nCol < m_nColumns - n)
			printf("\t...");
		for(nCol = MAX(nCol, m_nColumns - n); nCol < m_nColumns; nCol++)
			printf("\t%f", Get(nRow, nCol));
		printf("\n");
	}
	if(nRow < m_nRows - n)
		printf("...\n");
	for(nRow = MAX(nRow, m_nRows - n); nRow < m_nRows; nRow++)
	{
		printf("%f", Get(nRow, 0));
		for(nCol = 1; nCol < n && nCol < m_nColumns; nCol++)
			printf("\t%f", Get(nRow, nCol));
		if(nCol < m_nColumns - n)
			printf("\t...");
		for(nCol = MAX(nCol, m_nColumns - n); nCol < m_nColumns; nCol++)
			printf("\t%f", Get(nRow, nCol));
		printf("\n");
	}
}

void GMatrix::Transpose()
{
	double* pNewData = new double[m_nRows * m_nColumns];
	int nRow, nCol;
	for(nRow = 0; nRow < m_nRows; nRow++)
	{
		for(nCol = 0; nCol < m_nColumns; nCol++)
			pNewData[nCol * m_nRows + nRow] = m_pData[nRow * m_nColumns + nCol];
	}
	delete(m_pData);
	m_pData = pNewData;
	int nTmp = m_nRows;
	m_nRows = m_nColumns;
	m_nColumns = nTmp;
}

double GMatrix::ComputeTrace()
{
	int nMin = m_nColumns;
	if(m_nRows < nMin)
		nMin = m_nRows;
	double dSum = 0;
	int n;
	for(n = 0; n < nMin; n++)
		dSum += Get(n, n);
	return dSum;
}
/*
void GMatrix::Invert()
{
	GAssert(m_nRows == m_nColumns, "only square matrices supported");
	GAssert(m_nRows > 1, "needs to be at least 2x2");
	for(int i=1; i < actualsize; i++)
		data[i] /= data[0]; // normalize row 0
	for(int i=1; i < actualsize; i++)
	{ 
		for(int j=i; j < actualsize; j++)
		{ // do a column of L
			D sum = 0.0;
			for (int k = 0; k < i; k++)  
				sum += data[j*maxsize+k] * data[k*maxsize+i];
			data[j*maxsize+i] -= sum;
		}
		if(i == actualsize-1)
			continue;
		for(int j=i+1; j < actualsize; j++)
		{  // do a row of U
			D sum = 0.0;
			for (int k = 0; k < i; k++)
				sum += data[i*maxsize+k]*data[k*maxsize+j];
			data[i*maxsize+j] = (data[i*maxsize+j]-sum) / data[i*maxsize+i];
		}
	}
	for( int i = 0; i < actualsize; i++ )  // invert L
	{
		for( int j = i; j < actualsize; j++ )
		{
			D x = 1.0;
			if ( i != j )
			{
				x = 0.0;
				for ( int k = i; k < j; k++ ) 
					x -= data[j*maxsize+k]*data[k*maxsize+i];
			}
			data[j*maxsize+i] = x / data[j*maxsize+j];
		}
	}
	for( int i = 0; i < actualsize; i++ )   // invert U
	{
		for ( int j = i; j < actualsize; j++ )
		{
			if ( i == j ) continue;
			D sum = 0.0;
			for ( int k = i; k < j; k++ )
				sum += data[k*maxsize+j]*( (i==k) ? 1.0 : data[i*maxsize+k] );
			data[i*maxsize+j] = -sum;
		}
	}
	for( int i = 0; i < actualsize; i++ )   // final inversion
	{
		for ( int j = 0; j < actualsize; j++ )
		{
			D sum = 0.0;
			for ( int k = ((i>j)?i:j); k < actualsize; k++ )  
				sum += ((j==k)?1.0:data[j*maxsize+k])*data[k*maxsize+i];
			data[j*maxsize+i] = sum;
		}
	}
}
*/

void GMatrix::Solve(double* pVector)
{
	// Subtract out the non-diagonals
	GAssert(m_nRows == m_nColumns, "Expected a square matrix");
	int nRow, nCol, i;
	double dVal;
	for(nCol = 0; nCol < m_nColumns; nCol++)
	{
		for(nRow = 0; nRow < m_nRows; nRow++)
		{
			if(nRow == nCol)
				continue;
			dVal = m_pData[nRow * m_nColumns + nCol] / m_pData[nCol * m_nColumns + nCol];
			for(i = 0; i < m_nColumns; i++)
				m_pData[nRow * m_nColumns + i] -= dVal * m_pData[nCol * m_nColumns + i];
			pVector[nRow] -= dVal * pVector[nCol];
		}
	}

	// Normalize to make the diagonals one
	for(nRow = 0; nRow < m_nRows; nRow++)