mat.hpp

这是个人脸识别程序

// $mat\mat.hpp 1.5 milbo$ matrix class, elements are doubles.  Built for Microsoft VC 6.0.
//
// Warning: this is raw research code -- expect it to be quite messy.
//
// This file needs a few external functions: SysErr() Err() fEqual() lprintf(), ...
//
// Facilities Provided
// -------------------
//
// This is a C++ wrapper class for the Gnu Scientific library.  It provides
// two classes: Mat and MatView.  It allows matrix expressions like A + B
// where A and B are matrices: the + operator is overloaded.  It also has
// facilities for the inverse and other linear algebra functions. This library
// is good for non-sparse matrices with n in dozens or hundreds.  It is not good
// for the specialized matrix functions used in graphics.
//
// The MatView class allows different views of a matrix. For example, you
// can view a matrix column as a vector.  Or create a view that is a
// sub matrix of a larger matrix.  Matrix views point to the same underlying
// data as the matrix they are viewing -- they are an alias to the same data.
// A speed advantage of views is that they allow you to manipulate a matrix
// without copying the actual  elements of the matrix.
//
// A variable of the Mat class has a single element m which is a GSL
// matrix of doubles.  There are also class static variables (common to the
// entire class).  Complex matrices aren't supported -- thus when extracting
// eigenvalues, matrices must be symmetric.
//
// The approach is to make programmer errors as visible as possible, although
// this can often only be done at run time by issuing an error msg (for example
// for matrix conformance mismatches).  Convenience extensions like negative
// indexing aren't allowed.  This seems appropriate in this compiled
// non-interactive environment.
//
// Assignment
// ----------
//
// You can't assign different sized matrices using "=".
// Thus you will get a runtime error message if the size of mat1 size is not the
// same as the size as mat1 in the assignment
//		mat1 = mat2;
// There is an exception though: if mat1's size is 0 then it will automatically
// be resized by the above statement.  This allows code like the following:
//
//		Mat mat1; 		// default constructor
//		...
//		mat1 = mat2;	// this assignment could be in a different function
//
// You can change the default behaviour and make _all_ matrix
// redimensioning-on-assignment automatic by calling the member function
//		setAutoRedimensionOnAssignment(true);
//  Using this is not recommended because, in practice, the runtime check of matrix
//	conformance helps find bugs.
//
// If you want to assign different sized matrices, use assign().
//
// The () operator
// ---------------
//
// Thus matrix indexing is (iRow,iCol) or (y,x). This is the mathematical/linear algebra
// convention, but is opposite to what you may expect: (x,y).  Matrices are stored
// using row-major ordering i.e. elements in each row are stored contigously.
//
// Indexed matrices like A(1,2) or A(3) always refer to a single element of A,
// except when defining a matrix.
// If you want to refer to a whole column or row conveniently, use views e.g. row()
// Vectors have the same type as matrices, so you can mix and match
// vectors and matrices in expressions,  as long as the dimensions conform.
//
// dyadic ()	X(1,2) refers to an element of X if X is a matrix -- row 1 col 2
// 				X(1,2) is illegal if X is a vector (will get runtime SysErr msg)
// 				X(0,2) is okay if X is a row vector
//
// monadic ()	X(1)   refers to an element of X if X is a vector
// 				X(1)   is illegal if X is an two-dimesional matrix (runtime SysErr)
//
// Accessing a one-dimensional matrix using 2 subscripts (one at 0) is slightly
// faster than accessing it using one subscript.  Using just one subscript can make
// the programmer's intention clearer though.
//
//
// Declaring Vectors
// -----------------
//
// Vec v(3);  		creates a column vector e.g.  1		// identical to Mat v(3,1)
// 				    							  2
// 											      3
//
// Vec v(3,"");		creates a row vector e.g. [1 2 3]	// identical to Mat v(1,3)
//
// In the above declaration, "" can be any string.  A string type
// was chosen so (3,"") could be disambiguated from (3,i) where i is an int.
// The #define ROWVEC should be used for the string.
//
// Views
// -----
//
// TODO fill in
// Caveat: mat.srow(A.t()) gives view of temporary object, not what you may think
//
// Error handling, and speed
// -------------------------
//
// If you this use module and the GSL library with GSL_RANGE_CHECK enabled,
// then all matrix indexes are checked at run time.  To ease bug location this class
// sometimes precheck matrix indexes (redundantly) before calling GSL
// functions.  These range checks help a lot.  but of course all this
// checking slows things down.  So what you should do is build two binaries for
// this class and the GSL object files: one with range checking on, one
// with it off.  You also want to HAVE_INLINE defined to 1 in the GSL library (although
//
// The code runs 2..4 times faster with compiler optimizations enabled and
// GSL_RANGE_CHECK disabled (tested with MS compiler on test.cpp and the test.cpp parameter
// NBIG varying from 100 to 500).  Compiling takes noticeably longer with optimizations enabled.
//
// This class needs some exernal functions e.g. SysErr().  Look at the make file
// for test.cpp to see where these functions are defined.
// It is a good idea to point the GSL error handler to SysErr so error handling
// is uniform across the application. See mat/test.cpp
//
// Zero sized matrices
// -------------------
//
// Matrices with no rows or no cols have m==NULL.  This is convenient for new()
// and related funcs on arrays of matrices.
// For efficiency, not all member functions can be used on zero sized matrices.
// Accessing array elements using (x,y) is illegal on zero size arrays -- you will get a
// null pointer crash.  Changing this would introduce a not-worth-it inefficiency in operator().
// Some that functions that are legal on zero sized matrices:
//		dim() dimClear() getDim() nelems() nGetRows() nGetCols() free()
//
// Notes
// -----
//
// Matrix expressions like A + B are somewhat inefficient, because we have to create a
// hidden temporary matrix for the result.  But these expressions are notationally convenient.
//
// TO DO
// -----
//
// add matView(anotherMat) which gives a view of a mat using anotherMat's dimensions
// use func pointer instead of fprintf/lprintf in printAux
// remove operators from class so dyadic doubles ops can be in either order
//
// Naming
// ------
//
// TODO fill this in, f prefix etc.
// A space is inserted before the parameter list "(" when defining a function.
// There isn't a space when declaring but not defining the function.
//
//
// History
// -------
//
// The following people wrote the original gslwrap class. This class is a modification
// of that -- specifically from gslwrap-0.2\matrix_double.h
// 		Ramin Nakisa
// 		Brian Gough, rodney sparapani, Marcel Bosc, Nigel Wetters
//
// milbo jun 05, petaluma

#if !defined(mat_hpp)
#define mat_hpp

#include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include "gsl/config.h"		// must be included for correct build with inline correctly defined etc.
#include "gsl/gsl_nan.h"
#include "gsl/gsl_errno.h"
#include "gsl/gsl_math.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_eigen.h"
#include "mcommon.hpp"

#define OSTREAM_SUPPORT	1	// turn off if get ostream compiler err msgs when using this file

#define CONF_nMaxMatDim		10000	// mostly used for parameter sanity checking but also used sometimes for a max temp alloc
#define CONF_nMaxMatVec		30000	// maximum number of matrices in a matrix vec
#define MAX_TEMP_ARRAY_LEN	1000	// maximum length of temporary arrays used internally in some functions
#define ROWVEC				""		// used to define a row vec (as opposed to a column vec): use as 2nd index when defining the vec
#define NO_CLEAR			0		// for fClear parameter

#define Vec 	Mat					// vectors are matrices with 1 row or 1 col
#define VecView	MatView				// ditto

#define VX 0						// X,Y,Z coords for indexing
#define VY 1
#define VZ 2

// Utility for printing matrices during debugging. Prints matrix name and the matrix.
#define PrintMat(e)	 { if ((e).nrows() == 1) (e).print(#e " %dx%d="); else (e).print(#e " %dx%d=\n"); }
#define PrintMatT(e) { if ((e).ncols() == 1) (e).t().print(#e " %dx%d="); else (e).t().print(#e " %dx%d=\n"); }

static const char sMatOperatorSingleArgIllegal[] = \
	"Mat::operator() single arg illegal on mat with ncols>1 or nrows>1";

// Function argument for map function, see map function header for details. Arg is current_elem
typedef double (__cdecl *tpMatFunc)(double);

//-----------------------------------------------------------------------------
namespace GslMat
{
#if OSTREAM_SUPPORT	// added because can get VC err C2874: using-declaration causes a multiple declaration of 'ostream'
using std::ostream;
#endif

class Mat;
class MatView;

// externs for functions defined in mat.cpp
extern double OneElemToDouble(const Mat &A);
extern bool fEqualMat(const Mat &A, const Mat &B, double MaxDiff=1e-10);
extern bool fSameDim(const Mat &A, const Mat &B);
extern void CheckReasonableDim(size_t nrows, size_t ncols, int nMin, const char *sMsg, const char *sFile=NULL);
extern void CheckSameDim(const Mat &A, const Mat &B, const char *sMsg);
extern void CheckDim(const Mat &A, size_t nrows, size_t ncols, const char *sMsg);
extern void Check2Cols(const Mat &A, const char *sMsg);
extern void CheckIsVector(const Mat &A, const char *sMsg);
extern void CheckIsSquare(const Mat &A, const char *sMsg);
extern void CheckSameNbrRows(const Mat &A, const Mat &B, const char *sMsg);
extern void CheckSameNbrCols(const Mat &A, const Mat &B, const char *sMsg);
extern Vec	GetElemsInRange(const Mat &A, double GreaterThan, double LessThan, bool fInRange);
extern Vec	GetElemsInRange(const Mat &A, const Mat &GreaterThan, const Mat &LessThan, bool fInRange);
extern Vec 	SolveWithLU(const Mat &A, Vec &b, double *pSmallestPivot, Mat &LU, gsl_permutation **ppPerm);
extern Vec 	SolveWithLU(const Mat &A, Vec &b, double *pSmallestPivot=NULL);
extern Vec	SolveWithLUAgain(Mat &LU, Vec &b, gsl_permutation *pPerm);
extern Vec	RefineLU(Mat &A, Mat &LU, gsl_permutation *pPerm, Vec &x, Vec &b);
extern Vec	SolveWithSVD(Mat &A, const Vec &b, bool fVerbose=false);
extern void SingularValDecomp(Mat &A, Mat &VTranspose, Vec &S);
extern Vec  SolveWithCholesky(Mat &A, const Vec &b, bool fVerbose, bool *pfPosDef);
extern bool fPosDef(const Mat &A);						// true if A is positive definite
extern bool fPosDef1(const Mat &A, int nAllowedZeroEigs=0); 	// true if A is positive definite (alternative slower routine)
extern Mat  GetEigsForSymMat(const Mat &A, Vec &EigVals, bool fSort=true, bool fTestForSymmetry=true);
extern Mat  NormalizeCols(const Mat &A);
extern Mat	HilbertMat(size_t n);
extern Mat	IdentMat(size_t n);		// returns an identity mat of the given size
extern Vec	AllOnesMat(size_t n, const char *sIsRowVec=NULL);
extern double xAx(const Vec &x, const Mat &A);
extern bool fIsMatrixSymmetric(const Mat &A, double Min=0.0);
extern double SumElems(const Mat &m);
extern double AbsSumElems(const Mat &m);

//-----------------------------------------------------------------------------
class Mat
{
public:

// For convenient reference, this is what a gsl_matrix is (defined in gsl_matrix_double.h)
//	typedef struct
//	{
//	  size_t size1;			// nrows
//	  size_t size2;			// ncols
//	  size_t tda;
//	  double *data;
//	  gsl_block *block;
//	  int owner;			// non-zero if own the block pointed to data
//	} gsl_matrix;

gsl_matrix *m;	// the underlying GSL matrix

// the copy function doesn't actually need templates but will need so later when I extend this class

template<typename T> void copy(const T &other)		// used for copy constructor and assignment
	{
	if (static_cast<const void *>(this) == static_cast<const void *>(&other))
		return;
	if (other.m == NULL)	// zero sized other?
		free();
	else
		{
		if (m == NULL)		// zero sized this?
			m = gsl_matrix_alloc(other.nrows(), other.ncols());
		else if (nrows() != other.nrows() || ncols() != other.ncols())
			{
			if (fAutoRedimOnAssign)
				{
				gsl_matrix_free(m);
				m = gsl_matrix_alloc(other.nrows(), other.ncols());
				}
			else
				SysErr("mats not conformable assigning %dx%d to %dx%d",
					other.nrows(), other.ncols(), nrows(), ncols());
			}
		gsl_matrix_memcpy(m, other.m);
		}
	}
Mat() :m(NULL) { }									// default constructor

Mat(size_t nrows1, size_t ncols1, bool fClear=true)
{
if (nrows1 == 0 || ncols1 == 0)
	m = NULL;
else
	m = (fClear? gsl_matrix_calloc(nrows1, ncols1): gsl_matrix_alloc(nrows1, ncols1));
}