29 #include "bias_config.h"
31 #define _USE_MATH_DEFINES
34 #include <Base/Debug/Error.hh>
36 #include <Base/Math/tnt/tnt.h>
37 #include <Base/Math/tnt/cmat.h>
38 #include <Base/Math/tnt/fmat.h>
42 #ifdef BUILD_ALLOW_NON_INTEGRAL_MATRIX_ELEMENTS
43 #include "Base/Math/Polynom.hh"
48 #ifndef DOUBLE_EPSILON
49 # define DOUBLE_EPSILON 1E-12
52 template <
class T>
class Vector;
86 template <
class T=
double>
102 Matrix(
int rows,
int cols,
const T value)
103 : TNT::
Matrix<T>(rows, cols, value) {};
108 : TNT::
Matrix<T>(rows, cols, T(0)) {
112 for (
int r=0; r<rows; r++) {
114 (*this)[r][r] = T(1);
118 default: BIASERR(
"undefined MatrixInitType");
128 explicit Matrix(
const int rows,
const int cols,
const std::string & s);
138 for (
int i=0; i<v.
size(); i++)
139 (*
this)[i][0] = v[i];
186 {
return this->v_; };
189 {
return this->v_; };
193 {
return (
const T **)this->row_; }
196 {
return this->row_; }
200 {
return this->rowm1_; };
202 inline unsigned int GetRows()
const {
return this->m_; }
204 inline unsigned int GetCols()
const {
return this->n_; }
215 BIASASSERT( GetNumElements()>0 );
216 return GetData() + GetNumElements() - 1;
221 BIASASSERT( GetNumElements()>0 );
222 return GetData() + GetNumElements() - 1;
238 void GetSubMatrix(
const size_t startRow,
const size_t startCol,
239 const size_t numRows,
const size_t numCols,
272 void SetSubMatrix(
const size_t startRowInThis,
const size_t startColInThis,
274 const size_t startRowInSub,
const size_t startColInSub,
275 const size_t numRows,
const size_t numCols);
281 void SetSubMatrix(
const size_t startRowInThis,
const size_t startColInThis,
288 void SetSubMatrix(
const size_t startRowInThis,
const size_t startColInThis,
293 void SetRow(
const int row,
const Vector<T>& data);
297 void SetCol(
const int row,
const Vector<T>& data);
302 void Set(
const int row,
const int col,
const Matrix<T> &data);
308 void Set(
const int row,
const int col,
const Vector<T> &data);
314 void SetTranspose(
const int row,
const int col,
const Vector<T> &data);
322 {
return this->num_cols() * this->num_rows(); };
337 void GetMaxMin(T& max, T& min)
const;
340 void GetAbsMaxMin(T& max, T& min)
const;
353 inline T NormL1()
const;
357 inline double NormL2()
const;
361 inline double NormFrobenius()
const;
364 inline T Trace()
const;
370 BIASASSERT((*this).num_cols()==(*this).num_rows());
371 const unsigned int dim = (*this).num_rows();
372 const unsigned int subdim = dim-1;
375 return ((*
this)[0][0]);
377 else if (subdim==1)
return ((*
this)[0][0]*(*
this)[1][1]-
378 (*
this)[1][0]*(*
this)[0][1]);
381 for (
register unsigned int sub=0; sub<dim; sub++) {
385 for (
register unsigned int i=0; i<subdim; i++) {
388 for (
register unsigned int j=0; j<subdim; j++) {
389 SubMatrix[j][i] = (*this)[(sub+j+1)%dim][i+1];
393 d += SubMatrix.
DetSquare() * (*this)[sub][0];
416 { Add(scalar, *
this); };
432 inline void Add(
const T& scalar,
Matrix<T> &dest)
const;
438 { Sub(scalar, *
this); };
443 inline void Sub(
const T& scalar,
Matrix<T> &dest)
const;
449 { Multiply(scalar, *
this); };
454 inline void Multiply(
const T& scalar,
Matrix<T> &dest)
const;
468 inline void MultLeft(
const Matrix<T>& arg);
487 { Divide(scalar, *
this); };
492 inline void Divide(
const T& scalar,
Matrix<T> &dest)
const;
497 inline void DivideElementwise(
const Matrix<T> &arg);
506 inline void GetSystemMatrix(
Matrix<T> &dest)
const;
517 inline void SetZero();
523 BIASASSERT(GetCols()==GetRows());
524 const int num = GetCols();
525 for (
int r=0; r<num; r++){
526 for (
int c=r; c<num; c++){
527 (*this)[c][r] = (*this)[r][c] =
528 T(((*
this)[c][r] + (*
this)[r][c]) / 2.0);
561 void Fill(
const T& scalar);
577 void NormalizeColsToOne(
const int norm_row_index);
587 void NormalizeRowsToOne(
const int norm_col_index);
591 void NormalizeRows();
595 void NormalizeCols();
604 void ScaleRow(
int NoRow,T scale);
607 void ScaleCol(
int NoCol,T scale);
619 bool Load(
const std::string & filename);
626 bool Save(
const std::string & filename)
const;
633 int WriteMatlab(std::ostream& ostr,
const std::string& name)
const;
640 bool BinRead(
const char * filename)
const;
646 bool BinWrite(
const char * filename)
const;
653 bool IsZero(
double eps=0.0 )
const;
659 bool IsIdentity(
double eps=0.0 )
const;
669 void SwapRows(
const int i,
const int r);
692 { mat.
AddIP(scalar);
return mat; }
699 { mat.
AddIP(arg);
return mat; }
706 { mat.
SubIP(arg);
return mat; }
713 { mat.
AddIP(scalar);
return mat; }
732 BIASASSERT( argl.
size() == argr.
size() );
733 register const T* argP = argl.
GetData();
734 register const T* matP = argr.
GetData();
746 {
return !(argl==argr); }
755 if (this->num_rows() != dest.
num_rows() ||
756 this->num_cols() != dest.
num_cols())
757 dest.
newsize(this->num_rows(), this->num_cols());
758 register T* destP = dest.
GetData();
759 register const T* matP = GetData();
760 for(;matP <= GetDataLast(); matP++, destP++) {
761 *destP = *matP + scalar;
768 if (this->num_rows() != dest.
num_rows() ||
769 this->num_cols() != dest.
num_cols())
770 dest.
newsize(this->num_rows(), this->num_cols());
771 register T* destP = dest.
GetData();
772 register const T* matP = GetData();
773 for(;matP <= GetDataLast(); matP++, destP++) {
774 *destP = *matP - scalar;
781 if (this->num_rows() != dest.
num_rows() || this->num_cols() != dest.
num_cols())
782 dest.
newsize(this->num_rows(), this->num_cols());
783 register T* destP = dest.
GetData();
784 register const T* matP = GetData();
785 for(;matP <= GetDataLast(); matP++, destP++) {
786 *destP = *matP * scalar;
793 if (this->num_rows() != dest.
num_rows() || this->num_cols() != dest.
num_cols())
794 dest.
newsize(this->num_rows(), this->num_cols());
796 BIASERR(
"Division by Zero");
800 register T* destP = dest.
GetData();
801 register const T* matP = GetData();
802 for(;matP <= GetDataLast(); matP++, destP++) {
803 *destP = *matP / scalar;
811 if (this->num_rows() != arg.
num_rows() || this->num_cols() != arg.
num_cols()){
812 BIASERR(
"Unequal matrix dimensions");
815 register const T* argP = arg.
GetData();
816 register T* matP = GetData();
817 for(;matP <= GetDataLast(); matP++, argP++) {
818 *matP = *matP / *argP;
825 Matrix<T> Out(this->num_cols(),this->num_rows());
826 for (
int i=0; i<this->num_rows(); i++)
827 for (
int j=0; j<this->num_cols(); j++)
828 Out[j][i]= (*
this)[i][j];
835 unsigned int num_rows = this->num_rows();
836 unsigned int num_cols = this->num_cols();
837 BIASASSERT (num_rows==num_cols);
839 Matrix<T> SubMatrix(num_rows-1,num_cols-1);
840 for (
unsigned int i=0; i<num_rows; i++) {
841 for (
unsigned int j=0; j<num_cols; j++) {
842 for (
unsigned int k=0; k<num_rows-1; k++) {
843 unsigned int kk = k<i?k:k+1;
844 for (
unsigned int r=0; r<num_rows-1; r++) {
845 unsigned int rr = r<j?r:r+1;
846 SubMatrix[k][r] = (*this)[kk][rr];
849 Out[j][i]= (((i+j)%2==0)?1:-1) * SubMatrix.
DetSquare();
858 #ifdef BUILD_ALLOW_NON_INTEGRAL_MATRIX_ELEMENTS
859 for(
register T* dataP = GetData(); dataP <= GetDataLast(); dataP++)
862 memset(GetData(), 0, GetNumElements()*
sizeof(T));
870 for(
register const T* dataP = GetDataLast(); dataP >= GetData(); dataP--) {
871 result += *dataP > 0 ? *dataP : *dataP * -1;
881 for (
register const T* dataP = GetDataLast(); dataP >= GetData(); dataP--){
882 result += (double)((*dataP) * (*dataP));
886 #ifdef BUILD_ALLOW_NON_INTEGRAL_MATRIX_ELEMENTS
889 BIASERR(
"NormL2 is only pollsible for integral types");
905 BIASASSERT(this->num_cols() == this->num_rows())
907 for (
int x=0; x<this->num_rows(); x++)
908 result = result + (*
this)[x][x];
916 result.
num_rows() != this->num_rows()){
920 if (this->num_cols() != arg.
num_rows() ||
922 result.
num_rows() != this->num_rows()){
923 BIASERR(
"invalid matrix sizes "<<this->num_rows()<<
"x"<< this->num_cols()
931 register int x, y, i;
934 for (i=0; i<this->num_cols(); i++)
935 result[y][x]+=(*
this)[y][i]*arg[i][x];
950 arg.
Mult(*
this,result);
957 if (result.
size() != this->num_rows()) {
958 result.
newsize(this->num_rows());
961 if (this->num_cols() != arg.
size() ||
962 result.
size() != this->num_rows()) {
963 BIASERR(
"invalid matrix/vector sizes "<<this->num_rows()<<
"x"<< this->num_cols()
964 <<
" * "<<arg.
size()<<
"x1"
965 <<
" = "<<result.
size()<<
"x1");
972 for (y=0; y<result.
size(); y++)
973 for (i=0; i<this->num_cols(); i++)
974 result[y]+=(*
this)[y][i]*arg[i];
981 if (result.
size() != this->num_cols()) {
982 result.
newsize(this->num_cols());
985 if (this->num_rows() != arg.
size() ||
986 result.
size() != this->num_cols()) {
987 BIASERR(
"invalid matrix/vector sizes "
988 <<arg.
size()<<
"x1"<<
" * "
989 <<this->num_rows()<<
"x"<< this->num_cols()
990 <<
" = "<<result.
size()<<
"x1");
994 if (result.
size() != this->num_cols()) {
995 result.
newsize(this->num_cols());
999 for (y=0; y<this->num_cols(); y++)
1000 for (i=0; i<this->num_rows(); i++)
1001 result[y]+=(*
this)[i][y]*arg[i];
1009 result.
num_rows() != this->num_rows()){
1012 BIASWARN(
"Result-matrix has been resized!");
1016 if (this->num_cols() != arg.
num_cols()){
1017 BIASERR(
"invalid matrix sizes "<<this->num_rows()<<
"x"<< this->num_cols()
1025 register int x, y, i;
1026 for (x=0; x<result.
num_cols(); x++)
1027 for (y=0; y<result.
num_rows(); y++)
1028 for (i=0; i<this->num_cols(); i++)
1029 result[y][x]+=(*
this)[y][i]*arg[x][i];
1036 if (this->num_cols() != arg.
num_cols() ||
1037 this->num_rows() != arg.
num_rows()) {
1038 BIASERR(
"invalid matrix sizes "<<this->num_rows()<<
"x"<< this->num_cols()
1044 T* thisData=GetData();
1045 const T* argData=arg.
GetData();
1046 register unsigned int size=this->num_rows()*this->num_cols();
1047 for (
unsigned int i=0;i<size;i++) {
1048 *thisData = *thisData + T(*argData);
1049 thisData++; argData++;
1057 if (this->num_cols() != arg.
num_cols() ||
1058 this->num_rows() != arg.
num_rows()) {
1059 BIASERR(
"invalid matrix sizes "<<this->num_rows()<<
"x"<< this->num_cols()
1065 T* thisData=GetData();
1066 const T* argData=arg.
GetData();
1067 register unsigned int size=this->num_rows()*this->num_cols();
1068 for (
unsigned int i=0;i<size;i++) {
1069 *thisData = *thisData - T(*argData);
1070 thisData++; argData++;
1077 const unsigned int colsJ = this->num_cols(), rowsJ = this->num_rows();
1083 for (
unsigned int col = 0; col<colsJ; col++) {
1084 for (
unsigned int k = 0; k<rowsJ; k++) {
1085 dest[col][col]+=this->row_[k][col]*this->row_[k][col];
1089 for (
unsigned int r=1; r<colsJ; r++) {
1090 for (
unsigned int c=0; c<r; c++) {
1091 for (
unsigned int k = 0; k<rowsJ; k++) {
1092 const T* pJ_k_ = this->row_[k];
1093 dest[r][c] += pJ_k_[r] * pJ_k_[c];
1095 dest[c][r]=dest[r][c];
1116 #include <Base/Math/Operators.hh>
1118 #endif // _Matrix_hh_
Matrix(int rows, int cols, const T value)
assignment with a constant values for all elements (=set)
void SubIP(const T &scalar)
in place subtraction function
Matrix< T > & operator=(const Matrix< T > &A)
MatrixInitType
can be passed to matrix constructors to init the matrix with the most often used values ...
Subscript num_cols() const
class for column vectors with arbitrary size
Matrix< T > Transpose() const
transpose function, storing data destination matrix
void DivideElementwise(const Matrix< T > &arg)
elementwise division function in place
Matrix< T > & newsize(Subscript M, Subscript N)
int GetNumElements() const
Returns the total number of elements.
bool operator!=(const Matrix< T > &argl, const Matrix< T > &argr)
Matrix< T > & operator=(const TNT::Matrix< T > &mat)
assignment operators calling corresponding operator from base class if appropriate ...
Matrix(const Matrix< T > &A)
double NormL2() const
Return the L2 norm: a^2 + b^2 + c^2 + ...
void SetZero()
equivalent to matrix call
Matrix(const Vector< T > &v)
constructs Nx1 Matrix from N-Vector
class BIASMathBase_EXPORT Matrix
virtual void clear()
stl conform interface destroys Matrix JW
Matrix< T > & operator+=(Matrix< T > &mat, const Matrix< T > &arg)
const T * GetDataLast() const
Get a pointer to the last data element Do not use this on unitilized matrices with <= on pointers beca...
Vector< T > & newsize(Subscript N)
void SetZero()
Sets all values to zero.
T ** GetDataArray1() const
returns 1 based array to data access
unsigned int GetRows() const
Matrix< T > & operator+=(Matrix< T > &mat, const T &scalar)
operators
T * GetData()
get the pointer to the data array of the matrix (for faster direct memeory access) ...
is a 'fixed size' quadratic matrix of dim.
Matrix(int rows, int cols)
void DivideIP(const T &scalar)
in place division function
Matrix< T > & operator*=(Matrix< T > &mat, const T &scalar)
Matrix(const TNT::Matrix< T > &A)
class Vector3 contains a Vector of fixed dim.
bool operator==(const Matrix< T > &argl, const Matrix< T > &argr)
unsigned int GetCols() const
void Mult(const Matrix< T > &arg, Matrix< T > &result) const
matrix multiplication, result is not allocated
matrix class with arbitrary size, indexing is row major.
const T ** GetDataArray() const
returns zero based arry for data access
Matrix(int rows, int cols, const MatrixInitType &i)
init with standard form
const T * GetData() const
Matrix< T > Adjoint() const
computes the adjoint matrix
void Multiply(const T &scalar, Matrix< T > &dest) const
multiplication function, storing data destination matrix
void MultiplyIP(const T &scalar)
in place multiplication function
void Sub(const T &scalar, Matrix< T > &dest) const
substraction function, storing data destination matrix
void MultLeft(const Matrix< T > &arg)
in Place matrix multiplication this is equal to M = arg*M, but faster
void AddIP(const T &scalar)
in place addition function
void Add(const T &scalar, Matrix< T > &dest) const
addition function, storing data destination matrix
class BIASMathBase_EXPORT Vector
Subscript num_rows() const
void Divide(const T &scalar, Matrix< T > &dest) const
division function, storing data destination matrix
void SubIP(const Matrix< T > &arg)
Subtracts arg from this this -= arg.
void MultiplyWithTransposeOf(const Matrix< T > &arg, Matrix< T > &result) const
matrix matrix multiplication for multiplication with the transpose of the given matrix, result=this*arg^T.
T NormL1() const
Return the L1 norm: |a| + |b| + |c| + ...
void MakeSymmetric()
componentwise: this = 0.5(this + this^T) yields symmetric matrix only allowed for square shaped matri...
double NormFrobenius() const
Return Frobenius norm = sqrt(trace(A^t * A)).
Matrix< T > & operator-=(Matrix< T > &mat, const Matrix< T > &arg)
Matrix< T > & operator-=(Matrix< T > &mat, const T &scalar)
void GetSystemMatrix(Matrix< T > &dest) const
compute square system matrix dest = A^T * A
void Fill(const T &scalar)
fills complete Vector with scalar value