CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 00023 #ifndef _MatrixTmplDec_ 00024 #define _MatrixTmplDec_ 00025 00026 #include <iosfwd> 00027 #include "Array.hxx" 00028 #include "CLAM_Math.hxx" 00029 00030 namespace CLAM 00031 { 00032 00033 template <class T> 00034 class MatrixTmpl 00035 { 00036 public: 00037 MatrixTmpl(); 00038 MatrixTmpl(unsigned int dim1, unsigned int dim2); 00039 MatrixTmpl(const MatrixTmpl<T>& originalMatrix); 00040 ~MatrixTmpl(); 00041 int GetNumRows() const {return mNumRows;} 00042 int GetNumColumns() const {return mNumColumns;} 00043 int GetNumElements() const {return mNumRows*mNumColumns;} 00044 const Array <T>& GetBuffer() const {return MatrixBuffer();} 00045 Array <T>& GetBuffer() {return MatrixBuffer();} 00046 00047 T Sum() const 00048 { 00049 T sum=0; 00050 for (unsigned int i=0; i<(mNumRows * mNumColumns); i++) 00051 sum += MatrixBuffer()[i]; 00052 return sum; 00053 } 00054 00055 T Max() const 00056 { 00057 T max = MatrixBuffer()[0]; 00058 for (unsigned int i=1; i<(mNumRows*mNumColumns); i++) 00059 if (MatrixBuffer()[i] > max) 00060 max = MatrixBuffer()[i]; 00061 return max; 00062 } 00063 00064 T Min() const 00065 { 00066 T min = MatrixBuffer()[0]; 00067 for (unsigned int i=1; i<(mNumRows*mNumColumns); i++) 00068 if (MatrixBuffer()[i] < min) 00069 min = MatrixBuffer()[i]; 00070 return min; 00071 } 00072 00073 float Mean() const {return (float)(Sum())/(mNumRows*mNumColumns);} 00074 00075 // Print the matrix 00076 void Print() const; 00077 00078 // Determinant 00079 float GetDet() const 00080 { 00081 CLAM_ASSERT(mNumRows == mNumColumns, 00082 "Determinant only valid for square matrices"); 00083 float sum = 0; 00084 00085 if (mNumColumns == 1) 00086 return (*this)(0,0); 00087 00088 for (unsigned int i=0; i< mNumColumns;i++) // For the first row 00089 { 00090 sum += MatrixBuffer()[i] * pow((TData)-1,(TData)i) * (GetSubmatrix(0,i).GetDet()); 00091 } 00092 // i+1 : the matrix index are from [0..N-1] 00093 return sum; 00094 } 00095 00096 // Get Transposed Matrix 00097 MatrixTmpl<T> GetTrans() 00098 { 00099 MatrixTmpl<T> ret(mNumColumns, mNumRows); 00100 for (unsigned int i=0; i<mNumColumns; i++) 00101 for(unsigned int j=0; j<mNumRows; j++) 00102 ret(i,j) = (*this)(j,i); 00103 return ret; 00104 } 00105 00106 // Transpose the matrix 00107 void Trans() 00108 { 00109 MatrixTmpl<T> ret(mNumColumns, mNumRows); 00110 for (unsigned int i=0; i<mNumColumns; i++) 00111 for(unsigned int j=0; j<mNumRows; j++) 00112 ret(i,j) = (*this)(j,i); 00113 (*this) = ret; 00114 } 00115 // Reset the matrix (set all the elements to 0) 00116 void Reset() 00117 { 00118 MatrixTmpl<T> ret(mNumRows, mNumColumns); 00119 for (unsigned int i=0; i<mNumRows; i++) 00120 for(unsigned int j=0; j<mNumColumns; j++) 00121 ret(i,j) = 0; 00122 (*this) = ret; 00123 } 00124 00125 00126 // Invert the matrix 00127 void Invert() 00128 { 00129 CLAM_ASSERT(mNumRows == mNumColumns, 00130 "Inverse can only be calculated for square matrix"); 00131 MatrixTmpl<T> ret(mNumRows, mNumColumns); 00132 MatrixTmpl<T> aux(mNumRows-1, mNumColumns-1); 00133 for (unsigned int i=0; i<mNumRows; i++) 00134 for (unsigned int j=0; j<mNumRows; j++) 00135 { 00136 aux = GetSubmatrix(i,j); 00137 ret(i,j) = pow((TData)-1,(TData)i+j) * aux.GetDet(); 00138 } 00139 ret = (1 / GetDet()) * (ret.GetTrans()); 00140 (*this) = ret; 00141 } 00142 00143 // Get the inverse of a matrix 00144 MatrixTmpl<T> GetInverse() const 00145 { 00146 CLAM_ASSERT(mNumRows == mNumColumns, 00147 "Inverse can only be calculated for square matrix"); 00148 MatrixTmpl<T> ret(mNumRows, mNumColumns); 00149 for (unsigned i=0; i<mNumRows; i++) 00150 for (unsigned j=0; j<mNumColumns; j++) 00151 ret(i,j) = pow((TData)-1,(TData)i+j) * ( GetSubmatrix(i,j).GetDet() ); 00152 ret = (1/GetDet()) * (ret.GetTrans()); 00153 return ret; 00154 } 00155 00156 // Delete a Row and return a new Matrix Oject with reduced dims 00157 MatrixTmpl<T> GetDelRow(unsigned int row) const 00158 { 00159 MatrixTmpl<T> ret (mNumRows-1,mNumColumns); 00160 00161 for (unsigned i=0;i<row;i++) // copy unshifted part 00162 { 00163 for(unsigned j=0;j<mNumColumns;j++) 00164 ret(i,j) = MatrixBuffer()[i*mNumColumns+j]; 00165 } 00166 for (unsigned i=row;i<mNumRows-1;i++) // shift rest one row up 00167 { 00168 for(unsigned j=0;j<mNumColumns;j++) 00169 ret(i,j) = MatrixBuffer()[(i+1)*mNumColumns+j]; 00170 } 00171 return ret; 00172 } 00173 00174 // Delete Column 00175 MatrixTmpl<T> GetDelColumn(unsigned int column) const 00176 { 00177 MatrixTmpl<T> ret (mNumRows,mNumColumns-1); 00178 00179 for (unsigned i=0;i<mNumRows;i++) // shift matrix one column left 00180 { 00181 for(unsigned j=0;j<column;j++) 00182 ret(i,j) = MatrixBuffer()[i*mNumColumns+j]; 00183 } 00184 00185 for (unsigned i=0;i<mNumRows;i++) // shift matrix one column left 00186 { 00187 for(unsigned j=column;j<mNumColumns-1;j++) 00188 ret(i,j) = MatrixBuffer()[i*mNumColumns+j+1]; 00189 } 00190 return ret; 00191 } 00192 00193 // Get a Submatrix(i,j) deleting the i row and the j column 00194 MatrixTmpl<T> GetSubmatrix(unsigned int i, unsigned int j) const 00195 { 00196 MatrixTmpl<T> aux( GetDelRow(i) ); 00197 MatrixTmpl<T> ret( aux.GetDelColumn(j) ); 00198 return ret; 00199 } 00200 00201 // Convert matrix to his Submatrix(i,j). 00202 void Submatrix(unsigned int i, unsigned int j) 00203 { 00204 MatrixTmpl<T> aux( GetDelRow(i) ); 00205 (*this) = aux.GetDelColumn(j); 00206 } 00207 00208 // Solving linear equations [A][x]=[b] 00209 // friend MatrixTmpl<T> Solve(const MatrixTmpl<T>& A, const MatrixTmpl<T>& b); 00210 00211 // LU Decomposition 00212 // Decompose the matrix into a product of lower and upper triangular matrices. 00213 // Decompose(); 00214 // Compute forward/backward substitution on decomposed LU matrix product 00215 // Substitution(LU); 00216 00217 // Eigenvalues and Eigen 00218 // friend MatrixTmpl<T> Eigenvalues(const MatrixTmpl<T>& m); 00219 // friend MatrixTmpl<T> Eigenvectors(const MatrixTmpl<T>& m); 00220 00221 // Gram Schmidt Orthonormalization 00222 // friend MatrixTmpl<T> Orth(const MatrixTmpl<T>& m); 00223 00224 /**** Get/Set an element *****/ 00225 void SetAt (unsigned int iPosition, unsigned int jPosition, T element) 00226 { 00227 CLAM_ASSERT(iPosition < mNumRows, 00228 "Index beyond matrix dimension"); 00229 CLAM_ASSERT(iPosition >= 0, 00230 "Index beyond matrix dimension"); 00231 CLAM_ASSERT(jPosition < mNumColumns, 00232 "Index beyond matrix dimension"); 00233 CLAM_ASSERT(jPosition >= 0, 00234 "Index beyond matrix dimension"); 00235 00236 MatrixBuffer()[mNumColumns*iPosition+jPosition] = element; 00237 } 00238 00239 T GetAt (unsigned int iPosition, unsigned int jPosition) const 00240 { 00241 CLAM_ASSERT(iPosition < mNumRows, 00242 "Index beyond matrix dimension"); 00243 CLAM_ASSERT(iPosition >= 0, 00244 "Index beyond matrix dimension"); 00245 CLAM_ASSERT(jPosition < mNumColumns, 00246 "Index beyond matrix dimension"); 00247 CLAM_ASSERT(jPosition >= 0, 00248 "Index beyond matrix dimension"); 00249 00250 return MatrixBuffer()[mNumColumns*iPosition+jPosition]; 00251 } 00252 00253 // Get one column 00254 friend MatrixTmpl<T> GetColumn(unsigned int column, MatrixTmpl<T>& m ) 00255 { 00256 CLAM_ASSERT(column<m.mNumColumns, 00257 "Column beyond matrix dimensions"); 00258 CLAM_ASSERT(column >= 0, 00259 "Column beyond matrix dimensions"); 00260 00261 MatrixTmpl<T> ret(m.mNumRows, 1); // Column vector 00262 for (unsigned int i=0; i<m.mNumRows; i++) 00263 ret(i,0) = m(i,column); 00264 return ret; 00265 } 00266 00267 // Get one row 00268 friend MatrixTmpl<T> GetRow(unsigned int row, MatrixTmpl<T>& m) 00269 { 00270 CLAM_ASSERT(row < m.mNumRows, 00271 "Row beyond matrix dimensions"); 00272 CLAM_ASSERT(row >= 0, 00273 "Row beyond matrix dimensions"); 00274 00275 MatrixTmpl<T> ret(1, m.mNumColumns); // Row vector 00276 for (unsigned int i=0; i<m.mNumColumns; i++) 00277 ret(0,i) = m(row,i); 00278 return ret; 00279 } 00280 00281 /* Apply an unary function to all the elements of the matrix */ 00282 /* sin, cos, sqrt, cbrt, exp, log, log10, asin, acos, tan, atan */ 00283 void Apply( T (*f)(T) ) 00284 { 00285 for (unsigned int i=0; i<mNumRows; i++) 00286 for(unsigned int j=0; j<mNumColumns; j++) 00287 (*this)(i,j) = f( (*this)(i,j) ); 00288 } 00289 00290 void Apply( T (*f)(T,int),int parameter ) 00291 { 00292 for (unsigned int i=0; i<mNumRows; i++) 00293 for(unsigned int j=0; j<mNumColumns; j++) 00294 (*this)(i,j) = f( (*this)(i,j), parameter ); 00295 } 00296 00297 friend MatrixTmpl<T> GetApply( const MatrixTmpl<T> &m, double f(double) ) 00298 { 00299 MatrixTmpl<T> ret(m.mNumRows, m.mNumColumns); 00300 for (unsigned int i=0; i<m.mNumRows; i++) 00301 for(unsigned int j=0; j<m.mNumColumns; j++) 00302 ret(i,j) = T(f( m(i,j) )); 00303 return ret; 00304 } 00305 00306 friend MatrixTmpl<T> AbsMatrix(const MatrixTmpl<T> &m) // absolute value 00307 { 00308 MatrixTmpl<T> ret(m.mNumRows, m.mNumColumns); 00309 for (unsigned int i=0; i<m.mNumRows; i++) 00310 for(unsigned int j=0; j<m.mNumColumns; j++) 00311 ret(i,j) = abs( m(i,j) ); 00312 return ret; 00313 } 00314 00315 /**** Operators ****/ 00316 T& operator () (unsigned int iPosition,unsigned int jPosition) const { 00317 CLAM_ASSERT(iPosition < mNumRows, 00318 "Index beyond matrix dimension"); 00319 CLAM_ASSERT(iPosition >= 0, 00320 "Index beyond matrix dimension"); 00321 CLAM_ASSERT(jPosition < mNumColumns, 00322 "Index beyond matrix dimension"); 00323 CLAM_ASSERT(jPosition >= 0, 00324 "Index beyond matrix dimension"); 00325 00326 return ( MatrixBuffer()[mNumColumns*iPosition+jPosition] ); 00327 } 00328 00329 const MatrixTmpl<T>& operator = (const MatrixTmpl<T>& originalMatrix) 00330 { 00331 // MatrixBuffer() = originalMatrix.MatrixBuffer(); 00332 *mpMatrixBuffer = *(originalMatrix.mpMatrixBuffer); 00333 //*mpMatrixBuffer = originalMatrix.MatrixBuffer(); 00334 mNumRows = originalMatrix.mNumRows; 00335 mNumColumns = originalMatrix.mNumColumns; 00336 return *this; 00337 } 00338 00339 const MatrixTmpl<T>& operator = (const T element) 00340 { 00341 MatrixBuffer().SetSize(1); 00342 MatrixBuffer()[0] = element; 00343 mNumRows = mNumColumns = 1; 00344 return *this; 00345 } 00346 00347 const MatrixTmpl<T>& operator += (const MatrixTmpl<T>& newMatrix) 00348 { 00349 // Adding element by element 00350 CLAM_ASSERT(mNumRows == newMatrix.mNumRows,"Adding matrix with different dimensions"); 00351 CLAM_ASSERT(mNumColumns == newMatrix.mNumColumns,"Adding matrix with different dimensions"); 00352 for (unsigned int i=0; i< (mNumRows*mNumColumns) ; i++) 00353 MatrixBuffer()[i] += newMatrix.MatrixBuffer()[i]; 00354 return *this; 00355 } 00356 00357 const MatrixTmpl<T>& operator -= (const MatrixTmpl<T>& newMatrix) 00358 { 00359 // Substract element by element 00360 CLAM_ASSERT(mNumRows == newMatrix.mNumRows,"Substracting matrix with different dimensions"); 00361 CLAM_ASSERT(mNumColumns == newMatrix.mNumColumns,"Substracting matrix with different dimensions"); 00362 for (unsigned int i=0; i< (mNumRows*mNumColumns) ; i++) 00363 MatrixBuffer()[i] -= newMatrix.MatrixBuffer()[i]; 00364 return *this; 00365 } 00366 00367 friend MatrixTmpl<T> operator + (MatrixTmpl<T>& m1, MatrixTmpl<T>& m2) 00368 { 00369 CLAM_ASSERT(m1.mNumRows == m2.mNumRows,"Adding matrix with different dimensions"); 00370 CLAM_ASSERT(m1.mNumColumns == m2.mNumColumns,"Adding matrix with different dimensions"); 00371 MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); // Construction of an Vector object 00372 for (unsigned int i=0; i<ret.mNumRows; i++) 00373 for (unsigned int j=0; j<ret.mNumColumns; j++) 00374 ret(i,j) = m1(i,j) + m2(i,j); 00375 return ret; 00376 } 00377 00378 // add an element to all the matrix elements 00379 friend MatrixTmpl<T> operator + (const MatrixTmpl<T>& m1, const T & element) 00380 { 00381 MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); 00382 for (unsigned int i=0; i<ret.mNumRows; i++) 00383 for (unsigned int j=0; j<ret.mNumColumns; j++) 00384 ret(i,j) = m1(i,j) + element; 00385 return ret; 00386 } 00387 00388 friend MatrixTmpl<T> operator - (MatrixTmpl<T>& m1, MatrixTmpl<T>& m2) 00389 { 00390 CLAM_ASSERT(m1.mNumRows == m2.mNumRows,"Substracting matrix with different dimensions"); 00391 CLAM_ASSERT(m1.mNumColumns == m2.mNumColumns,"Substracting matrix with different dimensions"); 00392 MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); // Construction of an Vector object 00393 for (unsigned int i=0; i<ret.mNumRows; i++) 00394 for (unsigned int j=0; j<ret.mNumColumns; j++) 00395 ret(i,j) = m1(i,j) - m2(i,j); 00396 return ret; 00397 } 00398 00399 // substract an element to all the matrix elements 00400 friend MatrixTmpl<T> operator - (MatrixTmpl<T>& m1, const T element) 00401 { 00402 MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); 00403 for (unsigned int i=0; i<ret.mNumRows; i++) 00404 for (unsigned int j=0; j<ret.mNumColumns; j++) 00405 ret(i,j) = m1(i,j) - element; 00406 return ret; 00407 } 00408 00409 friend MatrixTmpl<T> operator * (T scalar,const MatrixTmpl<T>& m) 00410 { 00411 MatrixTmpl<T> mult(m.mNumRows, m.mNumColumns); 00412 for (unsigned int i=0; i<(mult.mNumRows*mult.mNumColumns); i++) 00413 mult.MatrixBuffer()[i] = scalar * m.MatrixBuffer()[i]; 00414 return mult; 00415 } 00416 00417 friend MatrixTmpl<T> operator * (const MatrixTmpl<T>& m1, const MatrixTmpl<T>& m2) 00418 { 00419 CLAM_ASSERT( m1.mNumColumns == m2.mNumRows,"Multiplications with incompatible arrays"); 00420 MatrixTmpl<T> ret(m1.mNumRows, m2.mNumColumns); 00421 for (unsigned int i=0; i<ret.mNumRows; i++) 00422 for (unsigned int j=0; j<ret.mNumColumns; j++) { 00423 ret(i,j) = 0; 00424 for( unsigned int k=0; k<m1.mNumColumns; k++) 00425 ret(i,j) += m1(i,k)*m2(k,j); 00426 } 00427 return ret; 00428 } 00429 00430 friend MatrixTmpl<T> operator / (const MatrixTmpl<T>& m, T scalar) 00431 { 00432 MatrixTmpl<T> ret(m.mNumRows, m.mNumColumns); 00433 for (unsigned int i=0; i<(ret.mNumRows*ret.mNumColumns); i++) 00434 ret.MatrixBuffer()[i] = m.MatrixBuffer()[i] / scalar; 00435 return ret; 00436 } 00437 00438 friend bool operator == (const MatrixTmpl<T>& m1, const MatrixTmpl<T>& m2) 00439 { 00440 if ( (m1.mNumRows == m2.mNumRows) && (m1.mNumColumns == m2.mNumColumns) && (m1.MatrixBuffer() == 00441 m2.MatrixBuffer()) ) 00442 return true; 00443 else 00444 return false; 00445 } 00446 00447 00448 Array <T>& MatrixBuffer() const { return *mpMatrixBuffer;} 00449 00450 protected: 00451 00452 // Dimensions 00453 unsigned int mNumRows; 00454 unsigned int mNumColumns; 00455 00456 // Buffer 00457 Array <T>* mpMatrixBuffer; 00458 00459 }; 00460 00461 template <class T> 00462 std::istream& operator >> (std::istream & stream, MatrixTmpl<T> & a); 00463 00464 template <class T> 00465 std::ostream& operator << (std::ostream & stream, const MatrixTmpl<T> & a); 00466 00467 } // namespace CLAM 00468 00469 #endif // _MatrixTmplDec_ 00470