mat.cpp

这是个人脸识别程序

// $mat\mat.cpp 1.5 milbo$
// Initiallly derived from gslwrap-0.2\matrix_double.cc
// Warning: this is raw research code -- expect it to be quite messy.
// milbo jun 05, petaluma

#include <iomanip>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#pragma warning(disable:4786)		// MSC compiler: disable warnings about extremely long variable names in STL
#include "cvec.hpp"
#include "mcommon.hpp"
#include "mfile.hpp"
#include "mat.hpp"
#include "matview.hpp"
#include "mchol.hpp"
#include "memstate.hpp"

#define USE_EQ_DOUBLES	1	// 1 to for match within tolerance when testing for div by zero etc.
#define TRACE			0	// for debug printfs

#define PRINT_FORMAT 	"% 8.2f"	// default print format

namespace GslMat
{
// Matrix printing format, global setting for all matrices

char Mat::sPrintFormat[20] = PRINT_FORMAT;	// for function print
int  Mat::nPrintWidth = 8;					// for operator<<, derived from sgPrintFormat
int  Mat::nPrintPrecision = 2;				// ditto

bool Mat::fAutoRedimOnAssign = false;

//-----------------------------------------------------------------------------
Mat::Mat (const double Vals[], size_t nrows1, size_t ncols1)
{
m = gsl_matrix_alloc(nrows1, ncols1);
int iVal = 0;
for (int iRow = 0; iRow < nrows1; iRow++)
	for (int iCol = 0; iCol < ncols1; iCol++)
		gsl_matrix_set(m, iRow, iCol, Vals[iVal++]);
}

//-----------------------------------------------------------------------------
Mat::Mat (const double Vals[], size_t nelems, const char sIsRowVec[]) // declare matrix as vector and init
{
int iVal = 0;
if (sIsRowVec)								// create row vector if ANY string
	{										// column matrix: one column, many rows
	m = gsl_matrix_alloc(1, nelems);
	for (int iCol = 0; iCol < nelems; iCol++)
		gsl_matrix_set(m, 0, iCol, Vals[iVal++]);
	}
else										// row matrix: one row, many columns
	{
	m = gsl_matrix_alloc(nelems, 1);
	for (int iRow = 0; iRow < nelems; iRow++)
		gsl_matrix_set(m, iRow, 0, Vals[iVal++]);
	}
}

//-----------------------------------------------------------------------------
void Mat::getDim (int *pnrows, int *pncols) const
{
if (m)
	{
	*pnrows = nrows();
	*pncols = ncols();
	}
else
	{
	*pnrows = 0;
	*pncols = 0;
	}
}

//-----------------------------------------------------------------------------
// If dimensions are same then matrix not touched.
// Else redimensions the matrix.  The new contents are unspecified.
// This function is quicker than dimKeep() or dimClear()

void Mat::dim (size_t nrows1, size_t ncols1)
{
if (nrows1 == 0 || ncols1 == 0)	// zero size?
	{
	if (m)
		{
		DASSERT(m->owner);		// only the owner can change the size (a view can't)
		gsl_matrix_free(m);
		m = NULL;
		}
	}
else if (!m)
	m = gsl_matrix_alloc(nrows1, ncols1);
else if (nrows() != nrows1 || ncols() != ncols1)
	{
	DASSERT(m->owner);		// only the owner can change the size (a view can't)
	gsl_matrix_free(m);
	m = gsl_matrix_alloc(nrows1, ncols1);
	}
}

void Mat::dim(size_t nelems, const char *sIsRowVec) 	// redimension, new contents unspecified
{
if (sIsRowVec)								// create row vector if ANY string
	dim(1, nelems);
else
	dim(nelems, 1);
}

void Mat::dim (const Mat &other)		// make this the same size as other
{
dim(other.nrows(), other.ncols());
}

//-----------------------------------------------------------------------------
// If dimensions are same then matrix not touched.
// Else clears the matrix when it actually redimensions it.

void Mat::dimClear (size_t nrows1, size_t ncols1)
{
if (nrows1 == 0 || ncols1 == 0)	// zero size?
	{
	if (m)
		{
		DASSERT(m->owner);	// only the owner can change the size (a view can't)
		gsl_matrix_free(m);
		m = NULL;
		}
	}
else if (!m)
	m = gsl_matrix_calloc(nrows1, ncols1);
else if (nrows() != nrows1 || ncols() != ncols1)
	{
	DASSERT(m->owner);
	gsl_matrix_free(m);
	m = gsl_matrix_calloc(nrows1, ncols1);
	}
}

//-----------------------------------------------------------------------------
// This preserves as much of the old data as possible.  If new matrix is
// bigger than or same size as the old matrix then all data will be preserved.
// Unused entries in the new matrix are cleared i.e. set to 0.

void Mat::dimKeep (size_t nrows1, size_t ncols1)
{
if (nrows1 == 0 || ncols1 == 0)	// zero size?
	{
	if (m)
		{
		DASSERT(m->owner);	// only the owner can change the size (a view can't)
		gsl_matrix_free(m);
		m = NULL;
		}
	}
else if (!m)
	m = gsl_matrix_calloc(nrows1, ncols1);
else if (nrows() != nrows1 || ncols() != ncols1)
	{
	DASSERT(m->owner);
	gsl_matrix *p = gsl_matrix_calloc( nrows1, ncols1);

	// copy as much of the data as will fit in the new matrix

	size_t nMinrows = __min(nrows1, nrows());
	size_t nMincols = __min(ncols1, ncols());
	for (size_t i = 0; i < nMinrows; i++)
		for (size_t j = 0; j < nMincols; j++)
			gsl_matrix_set(p, i, j, gsl_matrix_get(m, i, j));
	gsl_matrix_free(m);
	m = p;					// copy top structure only, not values, for speed
	}
}

//-----------------------------------------------------------------------------
// Redimension this to the same size as other and copy other into this
// If no redimensioning is needed this is about as fast as operator=

void Mat::assign (const Mat &other)
{
dim(other.nrows(), other.ncols());
copy(other);
}

//-----------------------------------------------------------------------------
// Helper function

const Mat Mat::reshapeAux (const Mat &A, size_t nNewRows, size_t nNewCols,
			bool fRowMajor,	size_t nStartRow, size_t nStartCol, size_t nrows1, size_t ncols1)
{
if (nrows1 == 0)
	nrows1 = A.nrows() - nStartRow;
if (ncols1 == 0)
	ncols1 = A.ncols() - nStartCol;

#if GSL_RANGE_CHECK
DASSERT(nStartRow >= 0 && nStartCol >= 0);

if (nStartRow + nrows1 > A.nrows() || nStartCol + ncols1 > A.ncols())
	SysErr("Mat::reshape rows %d+%d > %d or columns %d+%d > %d",
		nStartRow, nrows1, A.nrows(), nStartCol, ncols1, A.ncols());

if (nrows1 * ncols1 < nNewRows * nNewCols)
	SysErr("Mat::reshape old dimensions %dx%d < new dimensions %dx%d",
		nrows1, ncols1, nNewRows, nNewCols);
#endif

Mat m(nNewRows, nNewCols);
int i0 = 0, j0 = 0;
if (fRowMajor)
	{
	for (size_t j = nStartCol; j < nStartCol + ncols1; j++)
		for (size_t i = nStartRow; i < nStartRow + nrows1; i++)
			{
			m(i0, j0) = A(i,j);
			if (++i0 == nNewRows)
				{
				i0 = 0;
				j0++;
				if (j0 == nNewCols)
					return m;
				}
			}
	}
else
	{
	for (size_t i = nStartRow; i < nStartRow + nrows1; i++)
		for (size_t j = nStartCol; j < nStartCol + ncols1; j++)
			{
			m(i0, j0) = A(i,j);
			if (++j0 == nNewCols)
				{
				j0 = 0;
				i0++;
				if (i0 == nNewRows)
					return m;
				}
			}
	}
DASSERT(false);	// should never get here
return m;
}

//-----------------------------------------------------------------------------
// Return a submatrix of this, redimensioned to nNewRows and nNewCols
// If you just want a submatrix then use submatrix().
//
// nStartRow nStartCol nrows1 ncols1 specify the submatrix in the original matrix
//
// nNewRows and nNewCols specify the layout in the new matrix
//
// If nrows1=0 it will be set to max possible, ditto for ncols1
//
// Thus if just want to relayout the rows and cols of a matrix then set all
// of nStartRow nStartRow nStartCol nrows1 ncols1 to 0 and nNewRows and
// nNewCols to the new layout.
//
// If you want to access parts of a matrix without data copying, perhaps a
// view will suit you better -- see MatView in matview.hpp.  This is usually
// more efficient.

Mat	Mat::reshape (size_t nNewRows, size_t nNewCols, bool fRowMajor,
				size_t nStartRow, size_t nStartCol, size_t nrows1, size_t ncols1)
{
return reshapeAux(*this, nNewRows, nNewCols, fRowMajor, nStartRow, nStartCol, nrows1, ncols1);
}

//-----------------------------------------------------------------------------
// Make me a submatrix of myself, redimensioned to nNewRows and nNewCols
// See header comments in reshape().

Mat Mat::reshapeMe (size_t nNewRows, size_t nNewCols, bool fRowMajor,
				size_t nStartRow, size_t nStartCol, size_t nrows1, size_t ncols1)
{
Mat result(reshapeAux(*this, nNewRows, nNewCols, fRowMajor, nStartRow, nStartCol, nrows1, ncols1));
this->assign(result);
return *this;
}

//-----------------------------------------------------------------------------
void Mat::printAux (const char sMsg[], const char sFormat[], 	// helper for print()
					const char *sFile, FILE *pFile, int nMax, bool fPrintRowIndex) const
{
bool fCloseFile = false;

ASSERT(!(pFile && !sFile));	// must give a file name for error reporting if pFile is given

if (sFile && pFile == NULL)
	{
	pFile = fopen(sFile, "w");
	if (!pFile)
		Err("Can't open %s", sFile);
	fCloseFile = true;
	}
if (!sFile)	// make sure we have something to print in user messages
	sFile = "file";

if (nMax == 0)
	nMax = CONF_nMaxMatDim;
int nrows1 = __min(nMax, nrows());
int ncols1 = __min(nMax, ncols());

if (sMsg)
	if (pFile)
		fprintf(pFile, sMsg, nrows(), ncols());
	else
		lprintf(sMsg, nrows(), ncols());

for (size_t iRow = 0; iRow < nrows1; iRow++)
	{
	if (fPrintRowIndex)
		if (pFile)
			fprintf(pFile, "%3d: ", iRow);
		else
			lprintf("%3d: ", iRow);
   	for (size_t iCol = 0; iCol < ncols1; iCol++)
		if (pFile)
			if (sFormat)
				fprintf(pFile, sFormat, getElem(iRow, iCol));
			else
				fprintf(pFile, sPrintFormat, getElem(iRow, iCol));
		else
			if (sFormat)
				lprintf(sFormat, getElem(iRow, iCol));
			else
				lprintf(sPrintFormat, getElem(iRow, iCol));
	if (iRow < nrows1)
		{
		if (pFile)
			fprintf(pFile, "\n");
		else
			lprintf("\n");
		}
	}
if (pFile)
	fflush(pFile);	// flush it out so it can be seen even if application crashes later
else
	fflush(stdout);

if (fCloseFile)
	fclose(pFile);
}

void Mat::print (const char sMsg[], const char sFormat[],
					const char *sFile, FILE *pFile, int nMax, bool fPrintRowIndex) const
{
printAux(sMsg, sFormat, sFile, pFile, nMax, fPrintRowIndex);
}

void Mat::print (const char sMsg[], const char sFormat[],
					const char *sFile, FILE *pFile, int nMax, bool fPrintRowIndex)
{
printAux(sMsg, sFormat, sFile, pFile, nMax, fPrintRowIndex);
}

//-----------------------------------------------------------------------------
#if OSTREAM_SUPPORT
ostream& operator<< (ostream& os, const Mat &m)
{
using namespace std;
os.setf(ios::fixed);
os << setprecision(Mat::nPrintPrecision);
for (size_t iRow = 0; iRow < m.nrows(); iRow++)
	{
   	for (size_t iCol = 0; iCol < m.ncols(); iCol++)
		os << setw(Mat::nPrintWidth) << m.getElem(iRow, iCol) << " ";
	if (iRow < m.nrows()-1)
		os << endl;
	}
return os;
}
#endif

//-----------------------------------------------------------------------------
// You should init nrows() and ncols() before calling this routine

Mat& Mat::operator= (const double Vals[])	// assign an array of doubles
{
int iVal = 0;

for (int iRow = 0; iRow < nrows(); iRow++)
		for (int iCol = 0; iCol < ncols(); iCol++)
			gsl_matrix_set(m, iRow, iCol, Vals[iVal++]);

return *this;
}

//-----------------------------------------------------------------------------
bool Mat::operator== (const Mat &other) const
{
if (nrows() != other.nrows() || ncols() != other.ncols())
  return false;

for (size_t iRow = 0; iRow < nrows(); iRow++)
 	for (size_t iCol = 0; iCol < ncols(); iCol++)
		if (this->getElem(iRow, iCol) != other.getElem(iRow, iCol))
			return false;

return true;
}

//-----------------------------------------------------------------------------
Mat Mat::operator+ (const Mat &other) const
{
Mat result(*this);
gsl_matrix_add(result.m, other.m);
return result;
}

//-----------------------------------------------------------------------------
Mat Mat::operator+ (const double &d) const
{
Mat result(*this);
gsl_matrix_add_constant(result.m, d);
return result;
}

//-----------------------------------------------------------------------------
Mat operator+ (const double &d, const Mat &other)
{
Mat result(other);
gsl_matrix_add_constant(result.m, d);
return result;
}

//-----------------------------------------------------------------------------
Mat &Mat::operator+= (const double &d)
{
gsl_matrix_add_constant(m, d);
return *this;
}

//-----------------------------------------------------------------------------
Mat &Mat::operator+= (const Mat &other)
{
gsl_matrix_add(m, other.m);