mat.hpp

这是个人脸识别程序

Mat(size_t nelems, const char *sIsRowVec=NULL, bool fClear=true)		// declare mat as a vector TODO make fClear default false some time
	{												// sIsRowVec is a string instead of a bool to disambiguate type
	if (nelems == 0)
		m = NULL;
	else
		{
		if (sIsRowVec)								// create row vector if any string (use ROWVEC defined above)
			m = (fClear? gsl_matrix_calloc(1, nelems): gsl_matrix_alloc(1, nelems)); // row matrix: one row, many columns
		else
			m = (fClear? gsl_matrix_calloc(nelems, 1): gsl_matrix_alloc(nelems, 1)); // column matrix: one column, many rows
		}
	}
Mat(const double Vals[], size_t nrows1, size_t ncols1);					// declare and init from an array

Mat(const double Vals[], size_t nelems, const char *sIsRowVec=NULL);	// declare as vector and init

Mat(const Mat &other)		// copy constructor
{
if (other.m)
	{
	m = gsl_matrix_alloc(other.nrows(), other.ncols());
	copy(other);
	}
else
	m = NULL;
}
~Mat() { free(); }

size_t nelems() const 	{ if (m == NULL) return 0; else return m->size1 * m->size2; }

size_t nGetRows() const { if (m == NULL) return 0; else return m->size1; }
size_t nGetCols() const { if (m == NULL) return 0; else return m->size2; }

size_t nrows() const 	{ return m->size1; }
size_t ncols() const 	{ return m->size2; }

size_t size1() const 	{ return m->size1; }
size_t size2() const 	{ return m->size2; }

size_t M() 	const 		{ return m->size1; }
size_t N() 	const 		{ return m->size2; }

void   getDim(int *pnrows, int *pncols) const;

// release underlying GSL mat (gsl_matrix_free checks if owner before doing the free)

void   free() { if (m) gsl_matrix_free(m); m = NULL; }

// (i,j) operator -- for accessing mat elements

const double &operator()(size_t iRow, size_t iCol) const { return *gsl_matrix_ptr(m, iRow, iCol); }

double 		 &operator()(size_t iRow, size_t iCol) 	     { return *gsl_matrix_ptr(m, iRow, iCol); }

// (i) operator -- for accessing mat elements when matrix is treated like a vector

const double &operator()(size_t i) const
	{
	if (nrows() == 1)
		return *gsl_matrix_ptr(m, 0, i);
#if GSL_RANGE_CHECK
	if (ncols() != 1) 	// must have either 1 col or 1 row to be a vec
		SysErr(sMatOperatorSingleArgIllegal);
#endif
	return *gsl_matrix_ptr(m, i, 0);
	}

double &operator()(size_t i)
	{
	if (nrows() == 1)
		return *gsl_matrix_ptr(m, 0, i);
#if GSL_RANGE_CHECK
	if (ncols() != 1) 	// must have either 1 col or 1 row to be a vec
		SysErr(sMatOperatorSingleArgIllegal);
#endif
	return *gsl_matrix_ptr(m, i, 0);
	}

double	getElem(size_t iRow, size_t iCol) const 	 	 { return  gsl_matrix_get(m, iRow, iCol); }
void	setElem(size_t row, size_t col, const double &v) { gsl_matrix_set(m, row, col, v); }

void	dim(size_t nrows1, size_t ncols1);				// redimension, new contents unspecified
void	dim(size_t nelems, const char *sIsRowVec=NULL);	// redimension as vector, new contents unspecified
void	dim(const Mat &other);					// redimension to same size as other, new contents unspecified
void	dimKeep(size_t nrows1, size_t ncols1);	// redimension, keep old data
void	dimClear(size_t nrows1, size_t ncols1);	// redimension the mat and zero it
void	assign(const Mat &other);				// redimension to same size as other and copy other into this

// return a redimensioned submatrix of this matrix
Mat		reshape(size_t nNewRows, size_t nNewCols, bool fRowMajor=false,
				size_t nStartRow=0, size_t nStartCol=0, size_t nrows1=0, size_t ncols1=0);

Mat  	reshapeMe(size_t nNewRows, size_t nNewCols, bool fRowMajor=false,	// $ not yet in test.cpp
				size_t nStartRow=0, size_t nStartCol=0, size_t nrows1=0, size_t ncols1=0);

// return submatrix of this matrix
Mat		submat(size_t nStartRow, size_t nStartCol, size_t nrows1=0, size_t ncols1=0)
				{ return reshape(nrows1, ncols1, false, nStartRow, nStartCol); }

Mat 	t() const;							// transpose of this matrix

// assign copies all elements of other (use views to refer to another matrix without copying)
Mat&     	operator=(const Mat &other)	{ copy(other); return *this; }
Mat&       	operator=(const double Vals[]);	// assign an array of doubles

#if 0	// deprecated
// copyStruct is essentially a quick assign which doesn't copy elements (similar to a view)
// careful: after a copyStruct the data is shared and destruction of this or m will release the data
void 		copyStruct(const Mat &other)	{ this->m = other.m; }
#endif

bool 	    operator==(const Mat &other) const;	// does exact compare, use fEqualMat if want tolerance
bool 	    operator!=(const Mat &other) const { return !(*this == other); }

Mat			operator+(const Mat &other) const;
Mat 	    operator+(const double &f) const;
friend Mat  operator+(const double &f, const Mat &other);
Mat 	    &operator+=(const double &f);
Mat 	    &operator+=(const Mat &other);

Mat 	    operator-(const Mat &other) const;
Mat 	    operator-(const double &f) const;
friend Mat 	operator-(const double &f, const Mat &other);
Mat 	    &operator-=(const double &f);
Mat 	    &operator-=(const Mat &other);

Mat 	    operator*(const Mat &other) const;
Mat 	    operator*(const double &f) const;
friend Mat 	operator*(const double &f, const Mat &other);
Mat 	    &operator*=(const double &f);
Mat 	    &operator*=(const Mat &other);

Mat			operator/(const double &) const;
Mat 		&operator/=(const double &);

Mat 		mulElems(const Mat &other) const;
Mat 		mulElems(const Mat &other);

Mat 		divElems(const Mat &other) const;
Mat 		divElems(const Mat &other);

inline Mat  absElems() const;					// returns abs value of elems

Mat  		squareElems() const;				// returns square of elems

Mat 		roundElems() const;					// round to nearest integer

Mat			map(tpMatFunc pFunc) const;			// applies pFunc to each elem

double sum() const;								// scalar sum of all elements
double absSum() const;							// scalar sum of absolute value of elements
double trace() const;							// sum of diag elements
double dot(const Mat &other) const;				// dot product
double firstElem() const			{ return getElem(0,0); };
double lastElem() const				{ return getElem(nrows()-1, ncols()-1); };
double distSquared(Mat &other) 	  	{ return (*this - other).squareElems().sum(); }
double distSquared(const Mat &other) const { return (*this - other).squareElems().sum(); }
double maxElem() const 				{ return gsl_matrix_max(m); }
double maxAbsElem() const;						// max of absolute value of elements
double minElem() const 				{ return gsl_matrix_min(m); }
int    iMaxElem() const;						// index of max element
int    iMaxAbsElem() const;						// index of max of absolute value of elements
int    iMinElem() const;						// index of min element
double l1Norm() const 				{ return this->map(fabs).sum(); }
double l2Norm() const 				{ return sqrt(dot(*this)); }
double supNorm() const 				{ return this->map(fabs).maxElem(); }
double len() const 					{ return sqrt(dot(*this)); }				// same as l2Norm()
double lenSquared() const 			{ return dot(*this); }
Mat    fill(const double &d) 		{ gsl_matrix_set_all(m, d); return *this; }	// set all elems to a scalar value
void   fill(double Val, int iStartRow, int iStartCol, int nEndRow, int nEndCol);// set specified elems to a scalar value
Mat	   normalizeMe();							// uses l2Norm
Mat	   normalize() const;						// uses l2Norm
Mat    identityMe(double DiagVal=1); 			// turn mat into an diag mat (DiagVal on diag, others 0)
Mat	   transposeMe();				  			// transposes the mat, redims if necessary
Mat	   diagMe(double f);			  			// set diag elems to f, doesn't clear other elems
Mat	   zeroMe() 					{ gsl_matrix_set_zero(m); return *this;}	// set all elems to 0
bool   fSquare() const 				{ return nrows() == ncols(); }				// true if mat is square
bool   fVector() const 				{ return nrows() == 1 || ncols() == 1;} 	// true if mat is vec
bool   fIeeeNormal() const;						// true if fIeeeNormal() true for all elements
bool   fAnElemIsZero(int i = 0) const;			// true if an element is equal to i
bool   fARowIsZero(int i = 0) const;			// true if a row is all equal to i
bool   fAllZeroes(double Val = 0) const;		// true if matrix is all equal to Val
void   clearSmallElems(double MinVal) const;	// clear all elems that are absolute smaller than MinVal
double sumLnDiag() const;						// sum of log of abs value of diag elems, if elem is 0 return -DBL_MAX
double prodDiag() const;						// product of diag elements
double det(int sign=1) const;					// determinant
double lndet() const;							// log of determinant
Mat	   inverse() const;							// return inverse, bad results if mat is singular
Mat	   invertMe();								// invert this, bad results if mat is singular
Vec	   mat33TimesVec2(const Vec &v2) const;		// returns Vec(2)
Vec	   colAverages() const;						// vector of column averages
Vec	   colSums() const;							// vector of column sums
double colSum(int iCol) const;
double rowSum(int iRow) const;

// print() prints the mat, all parameters are optional
//
// If sMsg is specified it gets printed before the matrix
// nrows() and ncols() are passed as printf params to sMsg which you
// can use if you want (use %d in sMsg).
//
// sFormat is a printf style format string for printing matrice elements e.g. "%g "
// Default sFormat is the class static variable sPrintFormat,
// change the default with setPrintFormat().
//
// If sFile and pFile both null      prints using lprintf
// If pFile specified                writes to it (and uses sFile only for user msgs)
// If sFile specified but pFile null then opens sFile and writes to it and then closes it
//
// nMax is max rows and/or cols to print, 0 for all.
// If you specify nMax, will print min of nMax, nrows, ncols.
//
// If fPrintRowIndex is true then prefixes each row with the row index

void	print(const char sMsg[]=NULL, const char sFormat[]=NULL,
			  const char *sFile=NULL, FILE *pFile=NULL, int nMax=0, bool fPrintRowIndex=false);

void	print(const char sMsg[]=NULL, const char sFormat[]=NULL,
			  const char *sFile=NULL, FILE *pFile=NULL, int nMax=0, bool fPrintRowIndex=false) const;

#if OSTREAM_SUPPORT
friend ostream& operator<<(ostream& os, const Mat& m);
#endif

// Following two settings apply to all matrices, not just this one
void 	setPrintFormat(const char *sPrintFormat=NULL);
void 	setAutoRedimensionOnAssignment(bool f) const { fAutoRedimOnAssign = f; }

// For the file reading and writing routines, if pFile parameter is not null then
// the file is assumed to be already open and pFile is used for appending (and
// sFile is just for msgs to the user).
//

void write(const char *sFile, FILE *pFile=NULL, bool fVerbose=false,	// write to a file
			    const char sComment[]=NULL, const char sFormat[] = NULL, bool fLimitOutputWidth=false) const;

void writeAsCArray(const char sArrayName[], const char *sFile, FILE *pFile=NULL,
				const char *sType="double", bool fVector=false, bool fVerbose=false) const;

// sread() resizes the mat to match the mat in the file (prefix "s" means sread returns a possible err msg string)

char *sread(const char *sFile, FILE *pFile=NULL, char *sTag=NULL, bool fVerbose=false, bool fExitOnErr=true);	// read from a file

MatView		  view			(size_t iStartRow, size_t iStartCol, size_t nrows1, size_t ncols1, int Tda=0);
const MatView view			(size_t iStartRow, size_t iStartCol, size_t nrows1, size_t ncols1, int Tda=0) const;

VecView       row			(size_t iRow);				// select a single row of mat
const VecView row			(size_t iRow) const;
VecView       rowAsCol		(size_t iRow);				// view a single column as a row
const VecView rowAsCol		(size_t iRow) const;

VecView       col			(size_t iCol);				// select a single column of mat
const VecView col			(size_t iCol) const;
#if 0	// not possible: i*m->tda + j, no multiplier on j (need i*m->tda + j*m->tda2)
VecView       colAsRow		(size_t iCol);
const VecView colAsRow		(size_t iCol) const;
#endif

VecView       viewAsRow		();							// view entire mat as a row vector
const VecView viewAsRow		() const;
VecView       viewAsCol		();							// view entire mat as a col vector
const VecView viewAsCol		() const;

VecView       viewDiagAsCol	();							// view diag of mat as a col vector
const VecView viewDiagAsCol	() const;

MatView       viewAsSquare	();							// view matrix as a square matrix
const MatView viewAsSquare	() const;

private:

// Default matrix printing format, global setting for all matrices
// Change this with setPrintFormat()
// Use sFormat param of print() to change print format for a single case

static char sPrintFormat[20];	// for function print, inited in m.cpp to "% 8.2f"
static int  nPrintWidth;		// for operator<<, automatically derived from sgPrintFormat
static int  nPrintPrecision;	// ditto

// If you set fAutoRedimOnAssign, you can have statements like A = B where A is of
// different dimensions to B.  A will be automatically redimensioned to match B.
// This is a global setting for all matrices.
// There is an exception: matrices of size 0 are always redimensioned as needed
// if they are on the left hand side of an assignment (regardless of the setting
// of the flag).  This allows, for example: Mat a; A = anotherMat;

static bool	fAutoRedimOnAssign;

Mat	LU_decomp(gsl_permutation *pPerm=NULL, int *pSign=NULL) const;

void printAux(const char sMsg[], const char sFormat[], 			// helper for print()
			  const char *sFile, FILE *pFile, int nMax, bool fPrintRowIndex) const;

const 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);

};	// end class Mat
}   // end namespace GslMat

#endif // mat_hpp