mat.cpp

这是个人脸识别程序

return i;
}

//-----------------------------------------------------------------------------
double Mat::trace () const
{
double sum = 0;
#if GSL_RANGE_CHECK
void CheckIsSquare(const Mat &A, const char *sMsg);	// forward declaration
CheckIsSquare(*this, "Mat trace");
#endif
for (int iRow = 0; iRow < nrows(); iRow++)
	sum += getElem(iRow, iRow);
return sum;
}

//-----------------------------------------------------------------------------
Mat Mat::t () const 			// transpose of matrix
{
Mat result(ncols(), nrows());

for (int i = 0; i < nrows(); i++)
	for (int j = 0; j < ncols(); j++)
		gsl_matrix_set(result.m, j, i, gsl_matrix_get(m, i, j));

return result;
}

//-----------------------------------------------------------------------------
Mat Mat::transposeMe ()
{
assign(t());
return *this;
}

//-----------------------------------------------------------------------------
// Returns sum of log of absolute value of diag elems. If an elem is 0 returns -DBL_MAX
// Equivalent to product of diag elems, but without overflow
//
// TODO not yet tested in test.cpp

double Mat::sumLnDiag () const
{
#if GSL_RANGE_CHECK
CheckIsSquare(*this, "Mat sumLnDiag");
#endif
double sum = 0;
for (int iRow = 0; iRow < nrows(); iRow++)
	{
	double Elem = fabs(getElem(iRow, iRow));
	if (fEqual(Elem, 0))
		return -DBL_MAX;
	else
		sum += log(Elem);
	}
return sum;
}

double Mat::prodDiag () const	// product of diag elements
{
#if GSL_RANGE_CHECK
CheckIsSquare(*this, "Mat prodDiag");
#endif
double prod = 1;
for (int iRow = 0; iRow < nrows(); iRow++)
	prod *= getElem(iRow, iRow);
return prod;
}

//-----------------------------------------------------------------------------
double Mat::det (int sign) const
{
Mat a = *this;
a = a.LU_decomp();
return gsl_linalg_LU_det(a.m, sign);
}

//-----------------------------------------------------------------------------
double Mat::lndet () const
{
Mat a = *this;
a = a.LU_decomp();
return gsl_linalg_LU_lndet(a.m);
}

//-----------------------------------------------------------------------------
Mat Mat::inverse () const
{
gsl_permutation *pPerm = gsl_permutation_alloc(nrows());
Mat A = *this;
A = A.LU_decomp(pPerm);
Mat AInverse(nrows(), ncols());
int iRet = gsl_linalg_LU_invert(A.m, pPerm, AInverse.m);
DASSERT(iRet == 0);
gsl_permutation_free(pPerm);
return AInverse;
}

//-----------------------------------------------------------------------------
Mat Mat::invertMe ()
{
assign(inverse());
return *this;
}

//-----------------------------------------------------------------------------
Mat Mat::LU_decomp (gsl_permutation *pPerm, int *pSign) const
{
bool fRetPerm = (pPerm != NULL);
if (!pPerm)
	pPerm = gsl_permutation_alloc(nrows());
int Sign;
Mat result = *this;
int iRet = gsl_linalg_LU_decomp(result.m, pPerm, &Sign);

// check if matrix is singular i.e. has a zero pivot
for (int i = result.nrows()-1; i >= 0; i--)		// smallest last so look there 1st
	ASSERT(!fEqual(result(i,i), FLT_EPSILON));	// use FLT_EPSILON as a convenient small nbr

DASSERT(iRet == 0);

if (!fRetPerm)
	gsl_permutation_free(pPerm);
if (pSign)
	*pSign = Sign;
return result;
}

//-----------------------------------------------------------------------------
// Set matrix print format
// Expects a printf style format string like %g or %8.2f
// The format applies to all matrices, not just this one

void Mat::setPrintFormat (const char *sPrintFormat)
{
if (sPrintFormat == NULL)		// if null parameter then use default format
	sPrintFormat = PRINT_FORMAT;

DASSERT(strlen(sPrintFormat) < sizeof(Mat::sPrintFormat)-1);
sprintf(Mat::sPrintFormat, "%s ", sPrintFormat);

// get nPrintWidth and nPrintPrecision from sPrintFormat, also checks sPrintFormat
// TODO this may longer work with new PRINT_FORMAT, added space after %

DASSERT(sPrintFormat[0] == '%');

char c0 = 'x', c1 = 'x';
if (sPrintFormat[1] >= '0' && sPrintFormat[1] <= '9')
	sscanf(sPrintFormat+1, "%d%c%d%c", &nPrintWidth, &c0, &nPrintPrecision, &c1);
else
	sscanf(sPrintFormat+1, "%c", &c0);

#if 0	//TODO fix this some time
// sanity check format helps prevent crash later when used in a printf
DASSERT(c0 == 'f' || c0 == 'g' || c0 == 'e' || c1 == 'f' || c1 == 'g' || c1 == 'e');
#endif
}

//-----------------------------------------------------------------------------
// If pFile is NULL then it is opened.  If pFile is already open then we append
// to it (and sFile isn't used except for user messages).
//
// File format is:
//
// # optional comments each prefixed by #
// "optional mat tag string, can be read back or ignored by Mat::read, all chars except tilde allowed"
// # more optional comments
// {
// nrows ncols			// comments NOT allowed between { and }
// elements...
// }
//
// { and } are header and footer chars (for sanity checking when reading).

void Mat::write (const char *sFile, FILE *pFile, bool fVerbose, const char sComment[], const char sFormat[], bool fLimitOutputWidth) const
{
bool fCloseFile = false;

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

if (pFile == NULL)
	{
	pFile = fopen(sFile, "w");
	if (!pFile)
		Err("Can't open %s", sFile);
	fCloseFile = true;
	}
size_t nrows1 = nGetRows(), ncols1 = nGetCols();

if (!sFile)		// make sure we have something to print in user messages
	sFile = "file";

if (fVerbose)
	lprintf("Write %s %dx%d\n", sFile, nrows(), ncols());

if (sComment)
	fprintf(pFile, "# %s\n", sComment);

char sStart[] = "{ ", sEnd[] = "}\n";
if (2 != ::fwrite(sStart, 1, 2, pFile))
	Err("Can't write matrix to %s", sFile);

// we don't use gsl_matrix_fprintf because debugging is easier if matrix is laid out on a per row basis

if (0 == fprintf(pFile, "%d %d", nrows1, ncols1))
	Err("Can't write matrix to %s", sFile);

if (nrows1 == 1 && ncols1 < 20)		// a little tweak so small vectors can be seen more easily
	fprintf(pFile, " ");
else
	fprintf(pFile, "\n");

int ncolsMax = ncols1;
if (fLimitOutputWidth)
	{
	if (ncolsMax > 20)
		ncolsMax = int(sqrt(ncols1));
	if (ncolsMax > 20)
		ncolsMax = 20;
	}
for (int iRow = 0; iRow < nrows1; iRow++)
	for (int iCol = 0; iCol < ncols1; iCol++)
		{
		if (sFormat)
			{
			if (fprintf(pFile, sFormat, getElem(iRow, iCol)) < 1)
				Err("Can't write matrix to %s", sFile);
			}
		else if (fprintf(pFile, "%13g", getElem(iRow, iCol)) < 13)
			Err("Can't write matrix to %s", sFile);
		if (iCol % ncolsMax == (ncolsMax-1) || iCol == ncols1-1)	// print \n every 20th number and at end of row
			{
			if (fprintf(pFile, "\n") != 1)
				Err("Can't write matrix to %s", sFile);
			}
		else if (fprintf(pFile, " ") != 1)
			Err("Can't write matrix to %s", sFile);
		}

if (2 != ::fwrite(sEnd, 1, 2, pFile))
	Err("Can't write matrix to %s", sFile);

if (fCloseFile)
	fclose(pFile);
}

//-----------------------------------------------------------------------------
// Like Mat::write but writes the matrix as a C language array
// If sArrayName is NULL then we don't write the matrix name or terminating ";".
// This is used in WriteMatVecACArray to write multiple matrices as big one array.

void Mat::writeAsCArray (const char sArrayName[], const char *sFile, FILE *pFile,
						 const char *sType, bool fVector, bool fVerbose) const
{
bool fCloseFile = false;

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

if (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";

size_t nrows1 = nGetRows(), ncols1 = nGetCols();

if (fVerbose)
	lprintf("Write %s %s[%d];\n", sType, sArrayName, (nrows1>1? nrows1: ncols1));

char *sFormat;
double Max;
if (0 == strcmp(sType, "double"))
	{ sFormat = "%g"; Max = DBL_MAX; }
else if (0 == strcmp(sType, "float"))
	{ sFormat = "%#gf"; Max = FLT_MAX; } // will print as 1.23f, the # forces decimal point (needed for f suffix)
else if (0 == strcmp(sType, "int"))
	{ sFormat = "%d"; Max = INT_MAX; }
else if (0 == strcmp(sType, "short"))
	{ sFormat = "%d"; Max = SHRT_MAX; }
else if ( 0 == strcmp(sType, "char"))
	{ sFormat = "%d"; Max = SCHAR_MAX; }
else
	SysErr("WriteMatVecAsCArray: sType must be one of double, float, int, short, char");

if (sArrayName)
	{
	Fprintf(sFile, pFile, "#ifdef DEFS_ONLY\n");
	if (fVector)
		Fprintf(sFile, pFile, "extern %s %s[%d];\n", sType, sArrayName, (nrows1>1? nrows1: ncols1));
	else
		Fprintf(sFile, pFile, "extern %s %s[%d][%d];\n", sType, sArrayName, nrows1, ncols1);
	Fprintf(sFile, pFile, "#else\n");
	if (fVector)
		Fprintf(sFile, pFile, "%s %s[%d] = {\n", sType, sArrayName, (nrows1>1? nrows1: ncols1));
	else
		Fprintf(sFile, pFile, "%s %s[%d][%d] = {\n", sType, sArrayName, nrows1, ncols1);
	}
for (int iRow = 0; iRow < nrows1; iRow++)
	{
	Fprintf(sFile, pFile, "\t");
	if (!fVector)
		Fprintf(sFile, pFile, "{");
	for (int iCol = 0; iCol < ncols1; iCol++)
		{
		double Elem = getElem(iRow, iCol);
		if (Elem < -Max+1 || Elem > Max-1)
			{
			Warn("writeAsCArray %s overflow %g range %g %g sFormat %s", sType, Elem, -Max+1, Max-1, sFormat);
			Elem = (Elem < 0? -Max+1:Max-1);
			}
		if (sType[0] == 'i' || sType[0] == 's' || sType[0] == 'c')	// int or short or char type?
			fprintf(pFile, sFormat, iRound(Elem));					// fprintf not Fprintf for speed
		else
			fprintf(pFile, sFormat, Elem);
		if (iCol < ncols1-1)
			fprintf(pFile, ",");
		if ((iCol+1) % 12 == 0 && iCol < ncols1-1)				// limit output width (up to 12 numbers per line)
			fprintf(pFile, "\n\t ");
		}
	if (!fVector)
		Fprintf(sFile, pFile, "},");
	Fprintf(sFile, pFile, "\n");
	}
if (sArrayName)
	{
	Fprintf(sFile, pFile, "};\n");
	Fprintf(sFile, pFile, "#endif // DEFS_ONLY\n");
	}
if (fCloseFile)
	fclose(pFile);
}

//-----------------------------------------------------------------------------
// If pFile is NULL then it is opened.  If pFile is already open then we read
// from it (and sFile isn't used except for user messages).
//
// If we try to read a matrix and can't find the matrix header '{' then we return
// with an error message, if fExitOnErr is false.
//
// We copy the matrix tag string into sTag, if sTag is not null.
// The caller must ensure that there is enough space in sTag.
// A tag string is an optional string (with no ~) enclosed in quotes before the matrix first delimiter '{'.
// It can be used for storing filenames and other info for the matrix.
// The leading and closing quotes are removed before copying into sTag.
// If there is no tag string, sTag[0] is set to 0.
//
// File format: see Mat::write

char *Mat::sread (const char *sFile, FILE *pFile, char *sTag, bool fVerbose, bool fExitOnErr)
{
static char sErr[256];
char sTag1[256]; sTag1[0] = 0;	// mat prefix we read from file, possibly copied into sTag
if (sTag)
	sTag[0] = 0;

bool fCloseFile = false;

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

if (pFile == NULL)
	{
	DASSERT(sFile);
	pFile = fopen(sFile, "r");
	if (!pFile)
		if (fExitOnErr)
			Err("Can't open %s", sFile);
		else
			return "Can't open file";
	fCloseFile = true;
	}
size_t nrows, ncols;

if (!sFile)	// make sure we have something to print for user messages
	sFile = "file";


if (fVerbose)
	lprintf("Read %s ", sFile);

// search for first delimiter '{', skipping optional comments prefixed by # or "
char c;
do {
	if (1 != ::fread(&c, 1, 1, pFile))
		{
		sprintf(sErr, "Can't read matrix header from %s (no '{') [last valid tag string \"%s\"]", sFile, sTag1);
		if (fExitOnErr)
			Err(sErr);
		else
			{
			if (fCloseFile)
				fclose(pFile);
			return sErr;
			}
		}
	if (c == '#')						// comment?
		{
		while (c != '\n')				// if so, then skip it
			if (1 != ::fread(&c, 1, 1, pFile))
				{
				sprintf(sErr, "Can't read matrix header from %s (no '{') [last valid tag string \"%s\"]", sFile, sTag1);
				if (fExitOnErr)
					Err(sErr);
				else
					{
					if (fCloseFile)
						fclose(pFile);
					return sErr;
					}
				}
		}
	else if (c == '"')			// embedded string?
		{
		if (!fgets(sTag1, 256-1, pFile))
			{
			sprintf(sErr, "Can't read matrix header from %s (no '{') [last valid tag string \"%s\"]", sFile, sTag1);
			if (fExitOnErr)
				Err(sErr);
			else
				{
				if (fCloseFile)
					fclose(pFile);
				return sErr;
				}
			}
		if (sTag)
			{
			char *pDest = strchr(sTag1, '"');	// there may or may not be a terminating ", be lenient
			int iLen = strlen(sTag1)-1;
			if (pDest)
				{
				iLen--;
				*pDest = 0;
				}