From a74183ddb4ab8383bfe214b3745eb8a0a99ee47a Mon Sep 17 00:00:00 2001 From: Determinant Date: Thu, 25 Jun 2015 12:56:45 +0800 Subject: let HTK I/O implementation be a single package --- htk_io/src/KaldiLib/Matrix.h | 677 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 677 insertions(+) create mode 100644 htk_io/src/KaldiLib/Matrix.h (limited to 'htk_io/src/KaldiLib/Matrix.h') diff --git a/htk_io/src/KaldiLib/Matrix.h b/htk_io/src/KaldiLib/Matrix.h new file mode 100644 index 0000000..d33cb0c --- /dev/null +++ b/htk_io/src/KaldiLib/Matrix.h @@ -0,0 +1,677 @@ +#ifndef TNet_Matrix_h +#define TNet_Matrix_h + +#include +#include +#include +#include + +#ifdef HAVE_ATLAS +extern "C"{ + #include + #include +} +#endif + +#include "Common.h" +#include "MathAux.h" +#include "Types.h" +#include "Error.h" + +//#define TRACE_MATRIX_OPERATIONS +#define CHECKSIZE + +namespace TNet +{ + + + // class matrix_error : public std::logic_error {}; + // class matrix_sizes_error : public matrix_error {}; + + // declare the class so the header knows about it + template class Vector; + template class SubVector; + template class Matrix; + template class SubMatrix; + + // we need to declare the friend << operator here + template + std::ostream & operator << (std::ostream & rOut, const Matrix<_ElemT> & rM); + + // we need to declare the friend << operator here + template + std::istream & operator >> (std::istream & rIn, Matrix<_ElemT> & rM); + + // we need to declare this friend function here + template + _ElemT TraceOfProduct(const Matrix<_ElemT> &A, const Matrix<_ElemT> &B); // tr(A B) + + // we need to declare this friend function here + template + _ElemT TraceOfProductT(const Matrix<_ElemT> &A, const Matrix<_ElemT> &B); // tr(A B^T)==tr(A^T B) + + + /** ************************************************************************** + ** ************************************************************************** + * @brief Provides a matrix class + * + * This class provides a way to work with matrices in TNet. + * It encapsulates basic operations and memory optimizations. + * + */ + template + class Matrix + { + public: + /// defines a transpose type + + struct HtkHeader + { + INT_32 mNSamples; + INT_32 mSamplePeriod; + INT_16 mSampleSize; + UINT_16 mSampleKind; + }; + + + /** + * @brief Extension of the HTK header + */ + struct HtkHeaderExt + { + INT_32 mHeaderSize; + INT_32 mVersion; + INT_32 mSampSize; + }; + + + + + /// defines a type of this + typedef Matrix<_ElemT> ThisType; + + // Constructors + + /// Empty constructor + Matrix<_ElemT> (): + mpData(NULL), mMCols(0), mMRows(0), mStride(0) +#ifdef STK_MEMALIGN_MANUAL + , mpFreeData(NULL) +#endif + {} + + /// Copy constructor + Matrix<_ElemT> (const Matrix<_ElemT> & rM, MatrixTrasposeType trans=NO_TRANS): + mpData(NULL) + { if(trans==NO_TRANS){ Init(rM.mMRows, rM.mMCols); Copy(rM); } else { Init(rM.mMCols,rM.mMRows); Copy(rM,TRANS); } } + + /// Copy constructor from another type. + template + explicit Matrix<_ElemT> (const Matrix<_ElemU> & rM, MatrixTrasposeType trans=NO_TRANS): + mpData(NULL) + { if(trans==NO_TRANS){ Init(rM.Rows(), rM.Cols()); Copy(rM); } else { Init(rM.Cols(),rM.Rows()); Copy(rM,TRANS); } } + + /// Basic constructor + Matrix(const size_t r, const size_t c, bool clear=true) + { mpData=NULL; Init(r, c, clear); } + + + Matrix<_ElemT> &operator = (const Matrix <_ElemT> &other) { Init(other.Rows(), other.Cols()); Copy(other); return *this; } // Needed for inclusion in std::vector + + /// Destructor + ~Matrix() + { Destroy(); } + + + /// Initializes matrix (if not done by constructor) + ThisType & + Init(const size_t r, + const size_t c, bool clear=true); + + /** + * @brief Dealocates the matrix from memory and resets the dimensions to (0, 0) + */ + void + Destroy(); + + + ThisType & + Zero(); + + ThisType & + Unit(); // set to unit. + + /** + * @brief Copies the contents of a matrix + * @param rM Source data matrix + * @return Returns reference to this + */ + template ThisType & + Copy(const Matrix<_ElemU> & rM, MatrixTrasposeType Trans=NO_TRANS); + + + + /** + * @brief Copies the elements of a vector row-by-row into a matrix + * @param rV Source vector + * @param nRows Number of rows of returned matrix + * @param nCols Number of columns of returned matrix + * + * Note that rV.Dim() must equal nRows*nCols + */ + ThisType & + CopyVectorSplicedRows(const Vector<_ElemT> &rV, const size_t nRows, const size_t nCols); + + /** + * @brief Returns @c true if matrix is initialized + */ + bool + IsInitialized() const + { return mpData != NULL; } + + /// Returns number of rows in the matrix + inline size_t + Rows() const + { + return mMRows; + } + + /// Returns number of columns in the matrix + inline size_t + Cols() const + { + return mMCols; + } + + /// Returns number of columns in the matrix memory + inline size_t + Stride() const + { + return mStride; + } + + + /** + * @brief Gives access to a specified matrix row without range check + * @return Pointer to the const array + */ + inline const _ElemT* __attribute__((aligned(16))) + pData () const + { + return mpData; + } + + + /** + * @brief Gives access to a specified matrix row without range check + * @return Pointer to the non-const data array + */ + inline _ElemT* __attribute__((aligned(16))) + pData () + { + return mpData; + } + + + /** + * @brief pData_workaround is a workaround that allows SubMatrix to get a + * @return pointer to non-const data even though the Matrix is const... + */ + protected: + inline _ElemT* __attribute__((aligned(16))) + pData_workaround () const + { + return mpData; + } + public: + + + /// Returns size of matrix in memory + size_t + MSize() const + { + return mMRows * mStride * sizeof(_ElemT); + } + + /// Checks the content of the matrix for nan and inf values + void + CheckData(const std::string file = "") const + { + for(size_t row=0; row + operator [] (size_t i) const + { + assert(i < mMRows); + return SubVector<_ElemT>(mpData + (i * mStride), Cols()); + } + + inline SubVector<_ElemT> + operator [] (size_t i) + { + assert(i < mMRows); + return SubVector<_ElemT>(mpData + (i * mStride), Cols()); + } + + /** + * @brief Gives access to a specified matrix row without range check + * @return pointer to the first field of the row + */ + inline _ElemT* + pRowData(size_t i) + { + assert(i < mMRows); + return mpData + i * mStride; + } + + /** + * @brief Gives access to a specified matrix row without range check + * @return pointer to the first field of the row (const version) + */ + inline const _ElemT* + pRowData(size_t i) const + { + assert(i < mMRows); + return mpData + i * mStride; + } + + /** + * @brief Gives access to matrix elements (row, col) + * @return reference to the desired field + */ + inline _ElemT& + operator () (size_t r, size_t c) + { +#ifdef PARANOID + assert(r < mMRows && c < mMCols); +#endif + return *(mpData + r * mStride + c); + } + + /** + * @brief Gives access to matrix elements (row, col) + * @return pointer to the desired field (const version) + */ + inline const _ElemT + operator () (size_t r, size_t c) const + { +#ifdef PARANOID + assert(r < mMRows && c < mMCols); +#endif + return *(mpData + r * mStride + c); + } + + /** + * @brief Returns a matrix sub-range + * @param ro Row offset + * @param r Rows in range + * @param co Column offset + * @param c Coluns in range + * See @c SubMatrix class for details + */ + SubMatrix<_ElemT> + Range(const size_t ro, const size_t r, + const size_t co, const size_t c) + { return SubMatrix<_ElemT>(*this, ro, r, co, c); } + + const SubMatrix<_ElemT> + Range(const size_t ro, const size_t r, + const size_t co, const size_t c) const + { return SubMatrix<_ElemT>(*this, ro, r, co, c); } + /** @} */ + + + /** + * ********************************************************************** + * ********************************************************************** + * @defgroup MATH ROUTINES + * ********************************************************************** + * ********************************************************************** + * @{ + **/ + + /** + * @brief Returns sum of all elements + */ + _ElemT& + Sum() const; + + ThisType & + DotMul(const ThisType& a); + + ThisType & + Scale(_ElemT alpha); + + ThisType & + ScaleCols(const Vector<_ElemT> &scale); // Equivalent to (*this) = (*this) * diag(scale). + + ThisType & + ScaleRows(const Vector<_ElemT> &scale); // Equivalent to (*this) = diag(scale) * (*this); + + /// Sum another matrix rMatrix with this matrix + ThisType& + Add(const Matrix<_ElemT>& rMatrix); + + + /// Sum scaled matrix rMatrix with this matrix + ThisType& + AddScaled(_ElemT alpha, const Matrix<_ElemT>& rMatrix); + + /// Apply log to all items of the matrix + ThisType& + ApplyLog(); + + /** + * @brief Computes the determinant of this matrix + * @return Returns the determinant of a matrix + * @ingroup MATH + * + */ + _ElemT LogAbsDeterminant(_ElemT *DetSign=NULL); + + + /** + * @brief Performs matrix inplace inversion + */ + ThisType & + Invert(_ElemT *LogDet=NULL, _ElemT *DetSign=NULL, bool inverse_needed=true); + + /** + * @brief Performs matrix inplace inversion in double precision, even if this object is not double precision. + */ + ThisType & + InvertDouble(_ElemT *LogDet=NULL, _ElemT *DetSign=NULL, bool inverse_needed=true){ + double LogDet_tmp, DetSign_tmp; + Matrix dmat(*this); dmat.Invert(&LogDet_tmp, &DetSign_tmp, inverse_needed); if(inverse_needed) (*this).Copy(dmat); + if(LogDet) *LogDet = LogDet_tmp; if(DetSign) *DetSign = DetSign_tmp; + return *this; + } + + + /** + * @brief Inplace matrix transposition. Applicable only to square matrices + */ + ThisType & + Transpose() + { + assert(Rows()==Cols()); + size_t M=Rows(); + for(size_t i=0;i(const Matrix<_ElemT> &A, const Matrix<_ElemT> &B); // tr(A B) + friend _ElemT TNet::TraceOfProductT<_ElemT>(const Matrix<_ElemT> &A, const Matrix<_ElemT> &B); // tr(A B^T)==tr(A^T B) + friend class SubMatrix<_ElemT>; // so it can get around const restrictions on the pointer to mpData. + + /** ********************************************************************** + * ********************************************************************** + * @defgroup BLAS_ROUTINES BLAS ROUTINES + * @ingroup MATH + * ********************************************************************** + * ********************************************************************** + **/ + + ThisType & + BlasGer(const _ElemT alpha, const Vector<_ElemT>& rA, const Vector<_ElemT>& rB); + + ThisType & + Axpy(const _ElemT alpha, const Matrix<_ElemT> &rM, MatrixTrasposeType transA=NO_TRANS); + + ThisType & + BlasGemm(const _ElemT alpha, + const ThisType& rA, MatrixTrasposeType transA, + const ThisType& rB, MatrixTrasposeType transB, + const _ElemT beta = 0.0); + + + /** @} */ + + + /** ********************************************************************** + * ********************************************************************** + * @defgroup IO Input/Output ROUTINES + * ********************************************************************** + * ********************************************************************** + * @{ + **/ + + friend std::ostream & + operator << <> (std::ostream & out, const ThisType & m); + + void PrintOut(char *file); + void ReadIn(char *file); + + + bool + LoadHTK(const char* pFileName); + + /** @} */ + + + protected: +// inline void swap4b(void *a); +// inline void swap2b(void *a); + + + protected: + /// data memory area + _ElemT* mpData; + + /// these atributes store the real matrix size as it is stored in memory + /// including memalignment + size_t mMCols; ///< Number of columns + size_t mMRows; ///< Number of rows + size_t mStride; ///< true number of columns for the internal matrix. + ///< This number may differ from M_cols as memory + ///< alignment might be used + +#ifdef STK_MEMALIGN_MANUAL + /// data to be freed (in case of manual memalignment use, see Common.h) + _ElemT* mpFreeData; +#endif + }; // class Matrix + + template<> Matrix & Matrix::Invert(float *LogDet, float *DetSign, bool inverse_needed); // state that we will implement separately for float and double. + template<> Matrix & Matrix::Invert(double *LogDet, double *DetSign, bool inverse_needed); + + + + /** ************************************************************************** + ** ************************************************************************** + * @brief Sub-matrix representation + * + * This class provides a way to work with matrix cutouts in STK. + * + * + */ + template + class SubMatrix : public Matrix<_ElemT> + { + typedef SubMatrix<_ElemT> ThisType; + + public: + /// Constructor + SubMatrix(const Matrix<_ElemT>& rT, // Input matrix cannot be const because SubMatrix can change its contents. + const size_t ro, + const size_t r, + const size_t co, + const size_t c); + + + /// The destructor + ~SubMatrix<_ElemT>() + { +#ifndef STK_MEMALIGN_MANUAL + Matrix<_ElemT>::mpData = NULL; +#else + Matrix<_ElemT>::mpFreeData = NULL; +#endif + } + + /// Assign operator + ThisType& operator=(const ThisType& rSrc) + { + //std::cout << "[PERFORMing operator= SubMatrix&^2]" << std::flush; + this->mpData = rSrc.mpData; + this->mMCols = rSrc.mMCols; + this->mMRows = rSrc.mMRows; + this->mStride = rSrc.mStride; + this->mpFreeData = rSrc.mpFreeData; + return *this; + } + + + + /// Initializes matrix (if not done by constructor) + ThisType & + Init(const size_t r, + const size_t c, bool clear=true) + { Error("Submatrix cannot do Init"); return *this; } + + /** + * @brief Dealocates the matrix from memory and resets the dimensions to (0, 0) + */ + void + Destroy() + { Error("Submatrix cannot do Destroy"); } + + + + }; + + + + //Create useful shortcuts + typedef Matrix BfMatrix; + typedef SubMatrix BfSubMatrix; + + /** + * Function for summing matrices of different types + */ + template + void Add(Matrix<_ElemT>& rDst, const Matrix<_ElemU>& rSrc) { + assert(rDst.Cols() == rSrc.Cols()); + assert(rDst.Rows() == rSrc.Rows()); + + for(size_t i=0; i + void AddScaled(Matrix<_ElemT>& rDst, const Matrix<_ElemU>& rSrc, _ElemT scale) { + assert(rDst.Cols() == rSrc.Cols()); + assert(rDst.Rows() == rSrc.Rows()); + + Vector<_ElemT> tmp(rDst[0]); + + for(size_t i=0; i