mat.cpp

这是个人脸识别程序

return *this;
}

//-----------------------------------------------------------------------------
Mat Mat::operator- (const Mat &other) const
{
Mat result(*this);
gsl_matrix_sub(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(-1 * 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_sub(m, other.m);
return *this;
}

//-----------------------------------------------------------------------------
// Note: this isn't consistent with operator+ because we don't use
// matrix/gsl_matrix_mul here.  This is because for matrix multiply we have to
// redimension the output and therefore we use linalg/gsl_linalg_matmult.

Mat Mat::operator* (const Mat &other) const
{
Mat result(nrows(), other.ncols());
gsl_linalg_matmult(m, other.m, result.m);
return result;
}

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

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

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

//-----------------------------------------------------------------------------
Mat &Mat::operator*= (const Mat &other)
{
*this = (*this) * other;
return *this;
}

//-----------------------------------------------------------------------------
Mat Mat::operator/(const double &d) const
{
Mat result(*this);

#if USE_EQ_DOUBLES
DASSERT(!fEqual(d, 0));
#else
DASSERT(d != 0);
#endif

gsl_matrix_scale(result.m, 1.0 / d);
return result;
}

//-----------------------------------------------------------------------------
Mat &Mat::operator/= (const double &d)
{
#if USE_EQ_DOUBLES
DASSERT(!fEqual(d, 0));
#else
DASSERT(d != 0.0);
#endif
gsl_matrix_scale(m, 1.0 / d);
return *this;
}

//-----------------------------------------------------------------------------
Mat Mat::mulElems (const Mat &other) const
{
Mat result(*this);
gsl_matrix_mul_elements(result.m, other.m);
return result;
}

//-----------------------------------------------------------------------------
Mat Mat::mulElems (const Mat &other)
{
Mat result(*this);
gsl_matrix_mul_elements(result.m, other.m);
return result;
}

//-----------------------------------------------------------------------------
Mat Mat::divElems (const Mat &other) const
{
Mat result(*this);
gsl_matrix_div_elements(result.m, other.m);
return result;
}

//-----------------------------------------------------------------------------
Mat Mat::divElems (const Mat &other)
{
Mat result(*this);
gsl_matrix_div_elements(result.m, other.m);
return result;
}

//-----------------------------------------------------------------------------
// Apply the given function to each element in matrix.
// See the definitions of tpMatFunc at top of mat.hpp to see meaning of parameters.
//
// Examples:
//
// A.map(sqrt) returns a mat where each elem is the root of the coresponding element in A.
//
// double times2(double Elem) { return 2 * x; }
// A.map(times2) returns a mat where each elem is twice the coresponding element in A.

Mat Mat::map (tpMatFunc pFunc) const
{
Mat result(*this);
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		gsl_matrix_set(result.m, iRow, iCol, pFunc(getElem(iRow, iCol)));

return result;
}

//-----------------------------------------------------------------------------
// Utility function for squareElems

static inline double __cdecl square1 (double x)
{
return x * x;
}

//-----------------------------------------------------------------------------
Mat  Mat::squareElems () const
{
return map(square1);
}

//-----------------------------------------------------------------------------
inline Mat  Mat::absElems () const
{
return map(fabs);
}

//-----------------------------------------------------------------------------
// Utility function for roundElems

inline static double __cdecl round1 (double x)
{
// this is symmetrical for ".5" in that e.g. 1.5 becomes 1 and -1.5 becomes -1
if (x < 0)
	return (int)(x - 0.5);
else
	return (int)(x + 0.5);
}

//-----------------------------------------------------------------------------
Mat Mat::roundElems () const
{
return map(round1);
}

//-----------------------------------------------------------------------------
bool Mat::fIeeeNormal () const	// true if fIeeeNormal true for all elements of mat
{
extern inline bool fIeeeNormal(double x);	// mathutil.cpp

for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		if (!fIeeeNormal(getElem(iRow, iCol)))
			return false;
return true;
}

//-----------------------------------------------------------------------------
bool Mat::fAnElemIsZero (int i) const	// true if an element is zero, or equal to i
{
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		if (getElem(iRow, iCol) == i)
			return true;
return false;
}

//-----------------------------------------------------------------------------
bool Mat::fARowIsZero (int i) const		// true if a row is all zeroes, or all equal to i
{
for (int iRow = 0; iRow < nrows(); iRow++)
	{
	bool fAllZero = true;
	for (int iCol = 0; iCol < ncols(); iCol++)
		if (getElem(iRow, iCol) != i)
			{
			fAllZero = false;
			break;
			}
	if (fAllZero)
		return true;
	}
return false;
}

//-----------------------------------------------------------------------------
bool Mat::fAllZeroes (double Val) const		// true all values zero, or equal to Val
{
bool fAllZero = true;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		if (getElem(iRow, iCol) != Val)
			{
			fAllZero = false;
			break;
			}

return fAllZero;
}

//-----------------------------------------------------------------------------
// Clear (i.e. set to 0) all elems that are absolute smaller than MinVal
// If MinVal==0 then nothing will be changed

void Mat::clearSmallElems (double MinVal) const
{
DASSERT(MinVal >= 0);

if (MinVal == 0.0)		// for efficiency, not strictly needed
	return;

for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		if (fabs(gsl_matrix_get(m, iRow, iCol)) < MinVal)
			gsl_matrix_set(m, iRow, iCol, 0);
}

//-----------------------------------------------------------------------------
Mat Mat::identityMe (double DiagVal)
{
gsl_matrix_set_zero(m);
return diagMe(DiagVal);
}

//-----------------------------------------------------------------------------
Mat Mat::normalizeMe ()
{
double Len = this->len();

// We want to be able to normalize zero length vectors without worrying about
// about div by zero or checks exterior to this function. So we do a check here.

if (!fEqual(Len, 0))
	gsl_matrix_scale(this->m, 1.0 / Len);

return *this;
}

Mat Mat::normalize () const
{
Mat result(*this);
double Len = this->len();
if (!fEqual(Len, 0))
	gsl_matrix_scale(result.m, 1.0 / Len);
return result;
}

//-----------------------------------------------------------------------------
Mat Mat::diagMe (double d)
{
int n;
if (nrows() < ncols())
	n = nrows();
else
	n = ncols();
for (size_t i = 0; i < n; i++)
	gsl_matrix_set(m, i, i, d);
return *this;
}

//-----------------------------------------------------------------------------
double Mat::sum () const
{
double sum = 0;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		sum += getElem(iRow, iCol);
return sum;
}

//-----------------------------------------------------------------------------
// Dot product
// This is an efficient implementation of:	mulElems(other)).sum()

double Mat::dot (const Mat &other) const
{
const int nrows1 = this->nrows();
const int ncols1 = this->ncols();
ASSERT(other.nrows() == nrows1 && other.ncols() == ncols1);

double result = 0;
const int thisTda   = this->m->tda;
const int otherTda  = other.m->tda;
const double *thisData = this->m->data;
const double *otherData = other.m->data;

// for speed we check for special cases

if (nrows1 == 1)						// column vector?
	for (int j = 0; j < ncols1; j++)
		result += thisData[j] * otherData[j];

else if (ncols1 == 1)					// row vector?
	for (int i = 0; i < nrows1; i++)
		result += thisData[i * thisTda] * otherData[i * otherTda];

else 									// not a vector, must be a matrix
	for (int i = 0; i < nrows1; i++)
		{
		const int iThis  = i * thisTda;
		const int iOther = i * otherTda;
		for (int j = 0; j < ncols1; j++)
			result += thisData[iThis + j] * otherData[iOther + j];
		}

return result;
}

//-----------------------------------------------------------------------------
double Mat::absSum () const
{
double sum = 0;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		sum += fabs(getElem(iRow, iCol));
return sum;
}

//-----------------------------------------------------------------------------
double Mat::maxAbsElem () const
{
double max = -DBL_MIN;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		{
		double Elem = fabs(getElem(iRow, iCol));
		if (Elem  > max)
			max = Elem;
		}
return max;
}

//-----------------------------------------------------------------------------
int Mat::iMaxElem () const
{
double max = -DBL_MIN;
int i = 0;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		{
		double Elem = getElem(iRow, iCol);
		if (Elem  > max)
			{
			max = Elem;
			i = iCol + iRow * ncols();
			}
		}
return i;
}

int Mat::iMaxAbsElem () const
{
double max = -DBL_MIN;
int i = 0;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		{
		double Elem = fabs(getElem(iRow, iCol));
		if (Elem  > max)
			{
			max = Elem;
			i = iCol + iRow * ncols();
			}
		}
return i;
}

int Mat::iMinElem () const
{
double min = DBL_MAX;
int i = 0;
for (int iRow = 0; iRow < nrows(); iRow++)
	for (int iCol = 0; iCol < ncols(); iCol++)
		{
		double Elem = getElem(iRow, iCol);
		if (Elem  < min)
			{
			min = Elem;
			i = iCol + iRow * ncols();
			}
		}