CLAM-Development  1.4.0
MatrixTmplDec.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 #ifndef _MatrixTmplDec_
24 #define _MatrixTmplDec_
25 
26 #include <iosfwd>
27 #include "Array.hxx"
28 #include "CLAM_Math.hxx"
29 
30 namespace CLAM
31 {
32 
33  template <class T>
34  class MatrixTmpl
35  {
36  public:
37  MatrixTmpl();
38  MatrixTmpl(unsigned int dim1, unsigned int dim2);
39  MatrixTmpl(const MatrixTmpl<T>& originalMatrix);
40  ~MatrixTmpl();
41  int GetNumRows() const {return mNumRows;}
42  int GetNumColumns() const {return mNumColumns;}
43  int GetNumElements() const {return mNumRows*mNumColumns;}
44  const Array <T>& GetBuffer() const {return MatrixBuffer();}
46 
47  T Sum() const
48  {
49  T sum=0;
50  for (unsigned int i=0; i<(mNumRows * mNumColumns); i++)
51  sum += MatrixBuffer()[i];
52  return sum;
53  }
54 
55  T Max() const
56  {
57  T max = MatrixBuffer()[0];
58  for (unsigned int i=1; i<(mNumRows*mNumColumns); i++)
59  if (MatrixBuffer()[i] > max)
60  max = MatrixBuffer()[i];
61  return max;
62  }
63 
64  T Min() const
65  {
66  T min = MatrixBuffer()[0];
67  for (unsigned int i=1; i<(mNumRows*mNumColumns); i++)
68  if (MatrixBuffer()[i] < min)
69  min = MatrixBuffer()[i];
70  return min;
71  }
72 
73  float Mean() const {return (float)(Sum())/(mNumRows*mNumColumns);}
74 
75  // Print the matrix
76  void Print() const;
77 
78  // Determinant
79  float GetDet() const
80  {
82  "Determinant only valid for square matrices");
83  float sum = 0;
84 
85  if (mNumColumns == 1)
86  return (*this)(0,0);
87 
88  for (unsigned int i=0; i< mNumColumns;i++) // For the first row
89  {
90  sum += MatrixBuffer()[i] * pow((TData)-1,(TData)i) * (GetSubmatrix(0,i).GetDet());
91  }
92  // i+1 : the matrix index are from [0..N-1]
93  return sum;
94  }
95 
96  // Get Transposed Matrix
98  {
100  for (unsigned int i=0; i<mNumColumns; i++)
101  for(unsigned int j=0; j<mNumRows; j++)
102  ret(i,j) = (*this)(j,i);
103  return ret;
104  }
105 
106  // Transpose the matrix
107  void Trans()
108  {
110  for (unsigned int i=0; i<mNumColumns; i++)
111  for(unsigned int j=0; j<mNumRows; j++)
112  ret(i,j) = (*this)(j,i);
113  (*this) = ret;
114  }
115  // Reset the matrix (set all the elements to 0)
116  void Reset()
117  {
119  for (unsigned int i=0; i<mNumRows; i++)
120  for(unsigned int j=0; j<mNumColumns; j++)
121  ret(i,j) = 0;
122  (*this) = ret;
123  }
124 
125 
126  // Invert the matrix
127  void Invert()
128  {
130  "Inverse can only be calculated for square matrix");
133  for (unsigned int i=0; i<mNumRows; i++)
134  for (unsigned int j=0; j<mNumRows; j++)
135  {
136  aux = GetSubmatrix(i,j);
137  ret(i,j) = pow((TData)-1,(TData)i+j) * aux.GetDet();
138  }
139  ret = (1 / GetDet()) * (ret.GetTrans());
140  (*this) = ret;
141  }
142 
143  // Get the inverse of a matrix
145  {
147  "Inverse can only be calculated for square matrix");
149  for (unsigned i=0; i<mNumRows; i++)
150  for (unsigned j=0; j<mNumColumns; j++)
151  ret(i,j) = pow((TData)-1,(TData)i+j) * ( GetSubmatrix(i,j).GetDet() );
152  ret = (1/GetDet()) * (ret.GetTrans());
153  return ret;
154  }
155 
156  // Delete a Row and return a new Matrix Oject with reduced dims
157  MatrixTmpl<T> GetDelRow(unsigned int row) const
158  {
160 
161  for (unsigned i=0;i<row;i++) // copy unshifted part
162  {
163  for(unsigned j=0;j<mNumColumns;j++)
164  ret(i,j) = MatrixBuffer()[i*mNumColumns+j];
165  }
166  for (unsigned i=row;i<mNumRows-1;i++) // shift rest one row up
167  {
168  for(unsigned j=0;j<mNumColumns;j++)
169  ret(i,j) = MatrixBuffer()[(i+1)*mNumColumns+j];
170  }
171  return ret;
172  }
173 
174  // Delete Column
175  MatrixTmpl<T> GetDelColumn(unsigned int column) const
176  {
178 
179  for (unsigned i=0;i<mNumRows;i++) // shift matrix one column left
180  {
181  for(unsigned j=0;j<column;j++)
182  ret(i,j) = MatrixBuffer()[i*mNumColumns+j];
183  }
184 
185  for (unsigned i=0;i<mNumRows;i++) // shift matrix one column left
186  {
187  for(unsigned j=column;j<mNumColumns-1;j++)
188  ret(i,j) = MatrixBuffer()[i*mNumColumns+j+1];
189  }
190  return ret;
191  }
192 
193  // Get a Submatrix(i,j) deleting the i row and the j column
194  MatrixTmpl<T> GetSubmatrix(unsigned int i, unsigned int j) const
195  {
196  MatrixTmpl<T> aux( GetDelRow(i) );
197  MatrixTmpl<T> ret( aux.GetDelColumn(j) );
198  return ret;
199  }
200 
201  // Convert matrix to his Submatrix(i,j).
202  void Submatrix(unsigned int i, unsigned int j)
203  {
204  MatrixTmpl<T> aux( GetDelRow(i) );
205  (*this) = aux.GetDelColumn(j);
206  }
207 
208  // Solving linear equations [A][x]=[b]
209  // friend MatrixTmpl<T> Solve(const MatrixTmpl<T>& A, const MatrixTmpl<T>& b);
210 
211  // LU Decomposition
212  // Decompose the matrix into a product of lower and upper triangular matrices.
213  // Decompose();
214  // Compute forward/backward substitution on decomposed LU matrix product
215  // Substitution(LU);
216 
217  // Eigenvalues and Eigen
218  // friend MatrixTmpl<T> Eigenvalues(const MatrixTmpl<T>& m);
219  // friend MatrixTmpl<T> Eigenvectors(const MatrixTmpl<T>& m);
220 
221  // Gram Schmidt Orthonormalization
222  // friend MatrixTmpl<T> Orth(const MatrixTmpl<T>& m);
223 
224  /**** Get/Set an element *****/
225  void SetAt (unsigned int iPosition, unsigned int jPosition, T element)
226  {
227  CLAM_ASSERT(iPosition < mNumRows,
228  "Index beyond matrix dimension");
229  CLAM_ASSERT(iPosition >= 0,
230  "Index beyond matrix dimension");
231  CLAM_ASSERT(jPosition < mNumColumns,
232  "Index beyond matrix dimension");
233  CLAM_ASSERT(jPosition >= 0,
234  "Index beyond matrix dimension");
235 
236  MatrixBuffer()[mNumColumns*iPosition+jPosition] = element;
237  }
238 
239  T GetAt (unsigned int iPosition, unsigned int jPosition) const
240  {
241  CLAM_ASSERT(iPosition < mNumRows,
242  "Index beyond matrix dimension");
243  CLAM_ASSERT(iPosition >= 0,
244  "Index beyond matrix dimension");
245  CLAM_ASSERT(jPosition < mNumColumns,
246  "Index beyond matrix dimension");
247  CLAM_ASSERT(jPosition >= 0,
248  "Index beyond matrix dimension");
249 
250  return MatrixBuffer()[mNumColumns*iPosition+jPosition];
251  }
252 
253  // Get one column
254  friend MatrixTmpl<T> GetColumn(unsigned int column, MatrixTmpl<T>& m )
255  {
256  CLAM_ASSERT(column<m.mNumColumns,
257  "Column beyond matrix dimensions");
258  CLAM_ASSERT(column >= 0,
259  "Column beyond matrix dimensions");
260 
261  MatrixTmpl<T> ret(m.mNumRows, 1); // Column vector
262  for (unsigned int i=0; i<m.mNumRows; i++)
263  ret(i,0) = m(i,column);
264  return ret;
265  }
266 
267  // Get one row
268  friend MatrixTmpl<T> GetRow(unsigned int row, MatrixTmpl<T>& m)
269  {
270  CLAM_ASSERT(row < m.mNumRows,
271  "Row beyond matrix dimensions");
272  CLAM_ASSERT(row >= 0,
273  "Row beyond matrix dimensions");
274 
275  MatrixTmpl<T> ret(1, m.mNumColumns); // Row vector
276  for (unsigned int i=0; i<m.mNumColumns; i++)
277  ret(0,i) = m(row,i);
278  return ret;
279  }
280 
281  /* Apply an unary function to all the elements of the matrix */
282  /* sin, cos, sqrt, cbrt, exp, log, log10, asin, acos, tan, atan */
283  void Apply( T (*f)(T) )
284  {
285  for (unsigned int i=0; i<mNumRows; i++)
286  for(unsigned int j=0; j<mNumColumns; j++)
287  (*this)(i,j) = f( (*this)(i,j) );
288  }
289 
290  void Apply( T (*f)(T,int),int parameter )
291  {
292  for (unsigned int i=0; i<mNumRows; i++)
293  for(unsigned int j=0; j<mNumColumns; j++)
294  (*this)(i,j) = f( (*this)(i,j), parameter );
295  }
296 
297  friend MatrixTmpl<T> GetApply( const MatrixTmpl<T> &m, double f(double) )
298  {
300  for (unsigned int i=0; i<m.mNumRows; i++)
301  for(unsigned int j=0; j<m.mNumColumns; j++)
302  ret(i,j) = T(f( m(i,j) ));
303  return ret;
304  }
305 
306  friend MatrixTmpl<T> AbsMatrix(const MatrixTmpl<T> &m) // absolute value
307  {
309  for (unsigned int i=0; i<m.mNumRows; i++)
310  for(unsigned int j=0; j<m.mNumColumns; j++)
311  ret(i,j) = abs( m(i,j) );
312  return ret;
313  }
314 
315  /**** Operators ****/
316  T& operator () (unsigned int iPosition,unsigned int jPosition) const {
317  CLAM_ASSERT(iPosition < mNumRows,
318  "Index beyond matrix dimension");
319  CLAM_ASSERT(iPosition >= 0,
320  "Index beyond matrix dimension");
321  CLAM_ASSERT(jPosition < mNumColumns,
322  "Index beyond matrix dimension");
323  CLAM_ASSERT(jPosition >= 0,
324  "Index beyond matrix dimension");
325 
326  return ( MatrixBuffer()[mNumColumns*iPosition+jPosition] );
327  }
328 
329  const MatrixTmpl<T>& operator = (const MatrixTmpl<T>& originalMatrix)
330  {
331  // MatrixBuffer() = originalMatrix.MatrixBuffer();
332  *mpMatrixBuffer = *(originalMatrix.mpMatrixBuffer);
333  //*mpMatrixBuffer = originalMatrix.MatrixBuffer();
334  mNumRows = originalMatrix.mNumRows;
335  mNumColumns = originalMatrix.mNumColumns;
336  return *this;
337  }
338 
339  const MatrixTmpl<T>& operator = (const T element)
340  {
341  MatrixBuffer().SetSize(1);
342  MatrixBuffer()[0] = element;
343  mNumRows = mNumColumns = 1;
344  return *this;
345  }
346 
347  const MatrixTmpl<T>& operator += (const MatrixTmpl<T>& newMatrix)
348  {
349  // Adding element by element
350  CLAM_ASSERT(mNumRows == newMatrix.mNumRows,"Adding matrix with different dimensions");
351  CLAM_ASSERT(mNumColumns == newMatrix.mNumColumns,"Adding matrix with different dimensions");
352  for (unsigned int i=0; i< (mNumRows*mNumColumns) ; i++)
353  MatrixBuffer()[i] += newMatrix.MatrixBuffer()[i];
354  return *this;
355  }
356 
357  const MatrixTmpl<T>& operator -= (const MatrixTmpl<T>& newMatrix)
358  {
359  // Substract element by element
360  CLAM_ASSERT(mNumRows == newMatrix.mNumRows,"Substracting matrix with different dimensions");
361  CLAM_ASSERT(mNumColumns == newMatrix.mNumColumns,"Substracting matrix with different dimensions");
362  for (unsigned int i=0; i< (mNumRows*mNumColumns) ; i++)
363  MatrixBuffer()[i] -= newMatrix.MatrixBuffer()[i];
364  return *this;
365  }
366 
368  {
369  CLAM_ASSERT(m1.mNumRows == m2.mNumRows,"Adding matrix with different dimensions");
370  CLAM_ASSERT(m1.mNumColumns == m2.mNumColumns,"Adding matrix with different dimensions");
371  MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); // Construction of an Vector object
372  for (unsigned int i=0; i<ret.mNumRows; i++)
373  for (unsigned int j=0; j<ret.mNumColumns; j++)
374  ret(i,j) = m1(i,j) + m2(i,j);
375  return ret;
376  }
377 
378  // add an element to all the matrix elements
379  friend MatrixTmpl<T> operator + (const MatrixTmpl<T>& m1, const T & element)
380  {
381  MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns);
382  for (unsigned int i=0; i<ret.mNumRows; i++)
383  for (unsigned int j=0; j<ret.mNumColumns; j++)
384  ret(i,j) = m1(i,j) + element;
385  return ret;
386  }
387 
389  {
390  CLAM_ASSERT(m1.mNumRows == m2.mNumRows,"Substracting matrix with different dimensions");
391  CLAM_ASSERT(m1.mNumColumns == m2.mNumColumns,"Substracting matrix with different dimensions");
392  MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns); // Construction of an Vector object
393  for (unsigned int i=0; i<ret.mNumRows; i++)
394  for (unsigned int j=0; j<ret.mNumColumns; j++)
395  ret(i,j) = m1(i,j) - m2(i,j);
396  return ret;
397  }
398 
399  // substract an element to all the matrix elements
400  friend MatrixTmpl<T> operator - (MatrixTmpl<T>& m1, const T element)
401  {
402  MatrixTmpl<T> ret(m1.mNumRows, m1.mNumColumns);
403  for (unsigned int i=0; i<ret.mNumRows; i++)
404  for (unsigned int j=0; j<ret.mNumColumns; j++)
405  ret(i,j) = m1(i,j) - element;
406  return ret;
407  }
408 
409  friend MatrixTmpl<T> operator * (T scalar,const MatrixTmpl<T>& m)
410  {
412  for (unsigned int i=0; i<(mult.mNumRows*mult.mNumColumns); i++)
413  mult.MatrixBuffer()[i] = scalar * m.MatrixBuffer()[i];
414  return mult;
415  }
416 
418  {
419  CLAM_ASSERT( m1.mNumColumns == m2.mNumRows,"Multiplications with incompatible arrays");
420  MatrixTmpl<T> ret(m1.mNumRows, m2.mNumColumns);
421  for (unsigned int i=0; i<ret.mNumRows; i++)
422  for (unsigned int j=0; j<ret.mNumColumns; j++) {
423  ret(i,j) = 0;
424  for( unsigned int k=0; k<m1.mNumColumns; k++)
425  ret(i,j) += m1(i,k)*m2(k,j);
426  }
427  return ret;
428  }
429 
430  friend MatrixTmpl<T> operator / (const MatrixTmpl<T>& m, T scalar)
431  {
433  for (unsigned int i=0; i<(ret.mNumRows*ret.mNumColumns); i++)
434  ret.MatrixBuffer()[i] = m.MatrixBuffer()[i] / scalar;
435  return ret;
436  }
437 
438  friend bool operator == (const MatrixTmpl<T>& m1, const MatrixTmpl<T>& m2)
439  {
440  if ( (m1.mNumRows == m2.mNumRows) && (m1.mNumColumns == m2.mNumColumns) && (m1.MatrixBuffer() ==
441  m2.MatrixBuffer()) )
442  return true;
443  else
444  return false;
445  }
446 
447 
448  Array <T>& MatrixBuffer() const { return *mpMatrixBuffer;}
449 
450  protected:
451 
452  // Dimensions
453  unsigned int mNumRows;
454  unsigned int mNumColumns;
455 
456  // Buffer
458 
459  };
460 
461  template <class T>
462  std::istream& operator >> (std::istream & stream, MatrixTmpl<T> & a);
463 
464  template <class T>
465  std::ostream& operator << (std::ostream & stream, const MatrixTmpl<T> & a);
466 
467 } // namespace CLAM
468 
469 #endif // _MatrixTmplDec_
470