path: root/htk_io/src/KaldiLib/Matrix.h

#ifndef TNet_Matrix_h
#define TNet_Matrix_h

#include <stddef.h>
#include <stdlib.h>
#include <stdexcept>
#include <iostream>

#ifdef HAVE_ATLAS
extern "C"{
  #include <cblas.h>
  #include <clapack.h>
}
#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<typename _ElemT> class Vector;
  template<typename _ElemT> class SubVector;
  template<typename _ElemT> class Matrix;
  template<typename _ElemT> class SubMatrix;

  // we need to declare the friend << operator here
  template<typename _ElemT>
    std::ostream & operator << (std::ostream & rOut, const Matrix<_ElemT> & rM);

  // we need to declare the friend << operator here
  template<typename _ElemT>
    std::istream & operator >> (std::istream & rIn, Matrix<_ElemT> & rM);

  // we need to declare this friend function here
  template<typename _ElemT>
   _ElemT TraceOfProduct(const Matrix<_ElemT> &A, const Matrix<_ElemT> &B); // tr(A B)

  // we need to declare this friend function here
  template<typename _ElemT> 
	 _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<typename _ElemT>
    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<typename _ElemU>
      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<typename _ElemU> 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<Rows(); row++) {
          for(size_t col=0; col<Cols(); col++) {
            if(isnan((*this)(row,col)) || isinf((*this)(row,col))) {
              std::ostringstream os;
              os << "Invalid value: " << (*this)(row,col)
                 << " in matrix row: " << row
                 << " col: " << col
                 << " file: " << file;
              Error(os.str());
            }
          }
        }
      }

      /**
       *  **********************************************************************
       *  **********************************************************************
       *  @defgroup RESHAPE Matrix reshaping rutines
       *  **********************************************************************
       *  **********************************************************************
       * @{
       */

      /**
       *  @brief Removes one row from the matrix. The memory is not reallocated.
       */
      ThisType &
      RemoveRow(size_t i);      

      /** @} */

      /**
       *  **********************************************************************
       *  **********************************************************************
       *  @defgroup ACCESS Access functions and operators
       *  **********************************************************************
       *  **********************************************************************
       * @{
       */

      /**
       *  @brief Gives access to a specified matrix row without range check
       *  @return Subvector object representing the row
       */
      inline